Advances in Equation Solving and Symbolic Optimization

Adam Strzeboński, Wolfram Research Inc.

(RowDefault[, (Spacer1), , (Spacer1), ])  #WolframTechConf

ABSTRACT

Recent and upcoming releases of Mathematica include significant functionality extensions in functions for finding exact solutions of systems of equations and inequalities and for solving exact optimization problems. The extensions include methods for solving systems of transcendental function equations, exact convex optimization, and solving optimization problems that depend on symbolic parameters. In my talk, I will discuss the newly added algorithms and show some examples.

New methods in equation solving and symbolic optimization

Methods based on interval arithmetic

Mathematica 12.3 implements real and complex interval arithmetic for more than 100 mathematical functions (based on the Arb library). Availability of interval arithmetic allowed to develop several new equation solving and exact optimization methods.

• Representation of roots of systems of analytic function equations.

• Finding solution instances of systems of transcendental equations and inequalities.

• Solving systems of transcendental equations in bounded regions.

• Exact convex optimization.

• Counting roots of univariate meromorphic functions in bounded regions (upcoming version).

New special case parametric optimization methods

• Linear optimization (MPLP, multiparametric linear programming).

• Quadratic optimization (MPQP, multiparametric quadratic programming).

Center-radius interval arithmetic

({{, }, {real interval, complex interval}})

f[(Subscript[Δ, 1]),…,(Subscript[Δ, n])] yields a center-radius interval Δ which contains f[(Subscript[a, 1]),…,(Subscript[a, n])] for any (Subscript[a, i])∈(Subscript[Δ, i]).

In the upcoming version of Mathematica center-radius intervals will be available directly as CenteredInterval objects (the typesetting presented here is a prototype).

{a,b}={CenteredInterval[1,10^-10],CenteredInterval[2+3I,10^-10+10^-10I]}

BesselJ[Log[a]+1,a/b^2]

Representation of roots of systems of analytic function equations

Root[{{(Subscript[f, 1]),…,(Subscript[f, n])},{(Subscript[c, 1]),…,(Subscript[c, n])}},j] represents j-th coordinate of the unique root of the system of equations {(Subscript[f, 1])[(Subscript[x, 1]),…,(Subscript[x, n])]==0,…,(Subscript[f, n])[(Subscript[x, 1]),…,(Subscript[x, n])]==0} in the product of intervals/rectangles represented by {(Subscript[c, 1]),…,(Subscript[c, n])}. Interval arithmetic is used to prove existence and uniqueness of the root.

{f1,f2}={Sin[#1+Cos[#2]]-#1 #2+1&,Exp[#1 #2]-2#1Erf[#2]-3&};

c=SetPrecision[{x,y}/.FindRoot[{f1[x,y],f2[x,y]},{x,1},{y,1}],10]

Evaluation of Root verifies existence and uniqueness of the root.

Root[{{f1,f2},c},1]

N[%,100]

Finding solution instances of transcendental systems

FindInstance uses a new method which finds an approximate solution using FindRoot and then proves that a box around the approximate solution contains a unique root using interval arithmetic.

A real solution:

FindInstance[Gamma[x+y+z]==Sin[x]+Beta[x,y]+3&&SinIntegral[x y]==x-Cos[z]+2&&BesselJ[2,x-y]≥3 z,{x,y,z},Reals]

A complex solution:

FindInstance[Sinh[x y]-Exp[I y]==0&&2 x Log[y]-y Log[x]==1,{x,y}]

Solving systems of transcendental equations in bounded regions

Interval arithmetic provides criteria for proving that a box contains a unique root and for proving that a box contains no roots. Hence if a system contains a finite number of isolated simple roots in a box, then the roots can be found using a subdivision method.

A real system:

Solve[Gamma[x+y+1]-Sin[x y]==1&&Erf[x^2-y]-E^y-x+4==0&&0<x<3&&0<y<3,{x,y}, Reals]

A complex system:

Reduce[FresnelS[x-y]-AiryAi[x]+y==(3+I)/2&&CosIntegral[x] y-Sin[x y]==-5(1-I)/4&&1≤Re[x]≤3&&0≤Im[x]≤1&&1≤Re[y]≤3&&0≤Im[y]≤1,{x,y}]//AbsoluteTiming

Exact convex optimization

Convex optimization problems

Find the global minimum of a convex function f over the set

({{({{C={x∈}})(ℝ^n): (Subscript[g, 1])(x)≤0∧…∧(Subscript[g, m])(x)≤0∧(Subscript[h, 1])(x)==0∧…∧(Subscript[h, p])(x)==0}}})

where  (Subscript[g, 1]),…,(Subscript[g, m]) are convex functions and (Subscript[h, 1]),…,(Subscript[h, p]) are affine functions.

KKT conditions for convex optimization problems

Theorem (Sufficiency of KKT conditions)

Suppose that f,(Subscript[g, 1]),…,(Subscript[g, m]) are continuously differentiable over (ℝ^n), 
and there exist a∈C, (Subscript[λ, 1]),…,(Subscript[λ, m])≥0 , and (Subscript[μ, 1]),…,(Subscript[μ, p])∈ℝ such that

({{({{({{ (Subscript[∇, x])f(a)+(∑_(i=1)^m)(Subscript[λ, i])(Subscript[∇, x])(Subscript[g, i])(a)+(∑_(j=1)^p)(Subscript[μ, j])(Subscript[∇, x])(Subscript[h, j])(a)==0 }})
(Subscript[∀, 1≤i≤m])(Subscript[λ, i])(Subscript[g, i])(a)==0  }})}})

Then f(a)  is the global minimum of f on C.

Solving the KKT conditions

Functions f,(Subscript[g, 1]),…,(Subscript[g, m]) need not to be polynomial, so we cannot solve the KKT conditions directly using symbolic methods. However, we can solve convex optimization problems numerically. The new method finds an exact minimum point by using an approximate minimum point as a hint and then proving that the obtained exact point satisfies the KKT conditions. The algorithm proceeds as follows.

Find an approximate minimum point (OverTilde[a]).

From the approximate solution guess the set I⊆{1,…,m} such that (Subscript[∀, i∈I])(Subscript[g, i])(a)==0.

Find an exact real solution (a,λ,μ) with a near (OverTilde[a]) of the square system of equations

({{({{({{ (Subscript[∇, x])f(x)+(∑_(i∈I))(Subscript[λ, i])(Subscript[∇, x])(Subscript[g, i])(x)+(∑_(j=1)^p)(Subscript[μ, j])(Subscript[∇, x])(Subscript[h, j])(x)==0 }})
(Subscript[∀, i∈I])(Subscript[g, i])(x)==0 
(Subscript[∀, 1≤j≤p])(Subscript[h, j])(x)==0 }})}})

Verify that (Subscript[g, i])(a)≤0 for i∉I and (Subscript[λ, i])≥0 for i∈I.

Example

Find an exact solution to a convex optimization problem:

{min,pt}=Minimize[{Log[(E^(x+y))+(E^(2 x-y))],{x,y,x+y+1}(Underscript[\[VectorGreaterEqual], (ExponentialConeString)])0},{x,y}]

Show[Plot3D[Log[(E^(x+y))+(E^(2 x-y))],{x,-2,(5/4)},{y,-(1/4),1},MeshFunctions->{#3&},RegionFunction->Function[{x,y},{x,y,x+y+1}(Underscript[\[VectorGreaterEqual], (ExponentialConeString)])0]],Graphics3D[{Red,PointSize[0.03],Point[N[{x,y,min}/. pt]]}]]

Counting roots of univariate meromorphic functions in bounded regions

The algorithm

The following method finds the roots of a univariate meromorphic function in a bounded box:
1) Find approximate roots using a numeric method.
2) Prove existence and uniqueness of each root.
3) Prove that all roots have been found by counting the number of roots in the box.

Counting the number of roots can be done either by computing the winding number of the function along the boundary of the box, or by computing the integral of the logarithmic derivative of the function along the boundary of the box (and by computing the number of poles, if the function is not analytic).

The winding number can be computed using interval arithmetic, which in earlier versions of Mathematica was available only for elementary functions. For non-elementary functions the only available method was numerical computation of the integral of the logarithmic derivative.

Examples

Mathematica 12.3 uses numeric integration to count the roots here, so it issues a message warning that the completeness of the result depends on the correctness of the numeric integration result:

Reduce[Gamma[x]==Sin[Cos[x]]&&-5/2<=Re[x]<=5/2&&-1<=Im[x]<=1,x]

The development version result uses interval arithmetic for root counting, so the result has been proven to be correct:

Reduce[Gamma[x]==Sin[Cos[x]]&&-5/2<=Re[x]<=5/2&&-1<=Im[x]<=1,x]

Here the numeric integration used by Mathematica 12.3 times out:

f=Hypergeometric2F1[x,2x,1+x,3+x]-E^Sin[4x]+1+I;

Reduce[f==0&&0≤Re[x]≤1&&1/2≤Im[x]≤1,x]//AbsoluteTiming

Root counting based on interval arithmetic (development version) is much faster:

CountRoots[f,{x,I/2,1+I}]//AbsoluteTiming

Reduce[f==0&&0≤Re[x]≤1&&1/2≤Im[x]≤1,x]//AbsoluteTiming

Parametric linear optimization

The objective function has the form c.x and all constraints are linear in x and t, where x is a vector of optimization variables, and t is a vector of parameters.

Clear[a,b];m=MinValue[{-8x-7y,-6 x-4 y≤8-5 a-7 b&&7 x+5 y≤a+4 b&&-x+y≤6+4a-5 b&&-4 x+7 y≤-1-2 a-7 b&&5 y≤6-9 b},{x,y}]

Plot3D[m,{a,-3,3},{b,0,3}]

Piecewise solution regions:

regs=ImplicitRegion@@{#,{a,b}}&/@Internal`FromPiecewise[m][[1]];
RegionPlot@@{regs,PlotRange->{{-3,3},{0,3}}}

Parametric quadratic optimization

The objective function has the form (((Q.x+H.t+c))^T).x and all constraints are linear in x and t, where (Subscript[Q,  ]) is a positive-definite matrix, x is a vector of optimization variables, and t is a vector of parameters.

m=MinValue[{(x-1)^2+(2y-1)^2,x+2y≤a+b&&2x-y≤a-b+1&&x-2y≤2a-b+1},{x,y}]

Plot3D[m,{a,-5,5},{b,-5,5}]

Piecewise solution regions:

regs=ImplicitRegion@@{#,{a,b}}&/@Internal`FromPiecewise[m][[1]];
RegionPlot@@{regs,PlotRange->{{-5,5},{-5,5}}}

Timing comparison with V12.2

Arbitrary parametric polynomial optimization problems can be solved using cylindrical decomposition methods. However, the new specialized methods  for linear and quadratic problems are much faster.

Distance from the point {-1, 1} to the polygon given by 5 inequalities that depend on parameters a and b.

constr=x+2y≤a&&2x-y≤a-b+1&&x-2y≤2a-b+1&&-x-y≤2a+b&&y-x≤2a+b;

V12 .2 does not finish in 5 minutes:

TimeConstrained[(m=MinValue[{(x+1)^2+(y-1)^2,constr&&-1/2≤a≤1/2&&0≤b≤2},{x,y}])//AbsoluteTiming,300]

V12 .3:

(m=MinValue[{(x+1)^2+(y-1)^2,constr&&-1/2≤a≤1/2&&0≤b≤2},{x,y}])//AbsoluteTiming

Thank you for listening!