Symbolic Optimization

Adam Strzeboński, Wolfram Research

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

ABSTRACT

The Mathematica functions Minimize, Maximize, MinValue, MaxValue, ArgMin and ArgMax compute exact global extrema of univariate or multivariate functions, constrained by systems of equations and inequalities. The variables can range over the reals or over the integers. If symbolic parameters are present, the extrema are given as functions of parameters. In the recent releases of Mathematica, the exact global optimization capabilities have been significantly extended, and more optimization methods are currently under development. In my talk, I will show some examples and discuss the algorithms used.

The global optimization problem

Input

• A function f((Subscript[x, 1]),..., (Subscript[x, n])).
• A set of constraints on (Subscript[x, 1]),..., (Subscript[x, n]).
• Domain specifications for the variables (Subscript[x, 1]),..., (Subscript[x, n]) (either ℝ or ℤ).

Constraints are constructed from equations, inequalities, and region element statements, using Boolean operators and quantifiers .

Goals

• Find the global optimum (infimum or supremum) of values of f((Subscript[x, 1]),..., (Subscript[x, n])) on the solution set of the constraints.
• If possible, find a point satisfying the constraints at which the optimum is attained.
• If f((Subscript[x, 1]),..., (Subscript[x, n])) or the constraints contain parameters, find formulas for the optimum value and the optimum point in terms of the parameters.

Univariate Optimization

General analytic method

Find all zeros of the derivative, singular points, and limits at domain endpoints.

For some classes of functions, like rational functions, exp-log functions, or tame elementary functions, there are finitely many zeros and singularities over unbounded reals and we can compute them.

f=((E^x)/x)-(x^2) (Log[x]-(2/3));

Minimize[f,x]

Plot[f,{x,0,5},Epilog->{Red,PointSize[Large],Point[{x,f}/.%[[2]]]}]

For general analytic functions we can solve optimization problems over bounded intervals.

f=((BesselJ[2,x])/(Gamma[x]))+(((x+1))^(Sin[x]));

Maximize[{f,0≤x≤25},x]

Plot[f,{x,0,25},Epilog->{Red,PointSize[Large],Point[{x,f}/.%[[2]]]}]

Univariate Optimization

Methods using periodicity

For periodic functions it is sufficient to solve the optimization problem within a single period.

Here we have a composition of a periodic function with a function whose range includes a full period.

f=Log[Mod[(x^3)-1,2]+1]-Sin[π ((x^3)-2)];

Minimize[f,x]

Plot[f,{x,-3,3},Epilog->{Red,PointSize[Large],Point[{x,f}/.%[[2]]]}]

Find the maximum value of a combination of periodic functions with incommensurable periods. The maximum value is not attained.

f=Sin[Sin[(Sqrt[2]) x]]+2 Cos[3 x]-Erf[2+Cos[2 (Sqrt[2]) x]];

MaxValue[f,x]

Plot[f,{x,-50,50},Epilog->{Red,Line[{{-50,%},{50,%}}]}]

Find the minimum value of a combination of periodic functions with incommensurable periods over the integers.

f=(((Sin[π x]+Cos[Cos[E π x]]))/(Cos[(Sqrt[3]) π x]+2));

MinValue[f,x,Integers]

ListPlot[Table[f,{x,0,10000}],Epilog->{Red,Line[{{0,%},{10000,%}}]}]

Univariate Optimization

Methods using function property information

Mathematica contains information about properties of built-in mathematical functions, including domain, singularities, monotonicity, limits at domain endpoints, local and global minima, and more. This information is used to solve optimization problems.

This combines information of global minima of DawsonF and monotonicity of SinhIntegral.

f=SinhIntegral[DawsonF[x]];

Minimize[f,x]

Plot[f,{x,-5,5},Epilog->{Red,PointSize[Large],Point[{x,f}/.%[[2]]]}]

This maxes use of information about global maxima of BesselJ.

f=BesselJ[(7/4),(x^4)+x+1];

Maximize[f,x]

Plot[f,{x,-5,5},Epilog->{Red,PointSize[Large],Point[{x,f}/.%[[2]]]}]

Multivariate Optimization

Methods based on Cylindrical Algebraic Decomposition

Advantages

This is the most general complete method. The objective function and the constraints may involve arbitrary polynomials, rational functions, or real algebraic functions.

• No convexity requirements
• No degree bounds
• Constraints may involve quantifiers
• Works with parameters
• No numerical stability problems

Disadvantages

• Doubly-exponential complexity in the number of variables.

Examples

MaxValue[{x^3-a x y^2,x^2+b y^2≤1},{x,y}]

MinValue[{(x+1)^2+(y+7)^2,(Subscript[∀, z])(Subscript[∃, w])3 (z^2) w+(x+y) (z^4)-1==((x^3)-x y+(y^2)-1) (w^2)},{x,y}]

Numerical stability

CAD always finds a global extremum, even in cases that are numerically unstable. The left-hand side of the constraint here is ((((x^2)+(y^2)-(10^10)))^2) ((x^2)+(y^2)).

Maximize[{x+y,100000000000000000000 (x^2)-20000000000 (x^4)+(x^6)+100000000000000000000 (y^2)-40000000000 (x^2) (y^2)+3 (x^4) (y^2)-20000000000 (y^4)+3 (x^2) (y^4)+(y^6)≤1},{x,y}]//RootReduce

Here the input differs only in the 21-st decimal digit, but the answer is quite different, both in the minimum value and in the location of the minimum point.

Maximize[{x+y,100000000000000000001 (x^2)-20000000000 (x^4)+(x^6)+100000000000000000000 (y^2)-40000000000 (x^2) (y^2)+3 (x^4) (y^2)-20000000000 (y^4)+3 (x^2) (y^4)+(y^6)≤1},{x,y}]//RootReduce

Numeric optimization fails, even with very high working precision.

AbsoluteTiming[NMaximize[{x+y,100000000000000000000 (x^2)-20000000000 (x^4)+(x^6)+100000000000000000000 (y^2)-40000000000 (x^2) (y^2)+3 (x^4) (y^2)-20000000000 (y^4)+3 (x^2) (y^4)+(y^6)≤1},{x,y},WorkingPrecision->1000]]//N

Reduce`CADPlot[((((x^2)+(y^2)-10))^2) ((x^2)+(y^2))≤1,{x,-5,5},{y,-5,5}]

Reduce`CADPlot[((((x^2)+(y^2)-10))^2) ((x^2)+(y^2))+(x^2)≤1,{x,-5,5},{y,-5,5}]

Multivariate Optimization

Analytic methods

• Find all the local extrema, using the Lagrange multipliers or the KKT conditions, and pick the smallest (or the greatest).

• Needs to check smoothness of the objective function and smoothness and nondegeneracy of the constraints.

• Needs to check for extrema at the boundary of the set defined by the constraints and at infinity, if the set is unbounded.

Advantages

• Can solve some transcendental problems, as long as they lead to systems of equations that Mathematica can solve.

Disadvantages

• Requires ability to solve (exactly) systems of equations and inequalities over the reals. A general algorithm (CAD) exists only for algebraic problems, and in this case solving the generated systems may be harder than running the optimization-by-CAD algorithm.

• The implemented methods require bounded solution sets of constraints with simple to compute boundaries and finite sets of stationary and singular points.

Examples

Minimize[{y+Sin[2x],(y^3)==Cos[5 x]&&-5≤x≤5&&-5≤y≤5},{x,y}]

Minimize[{y+(E^x),(x^2)+(y^2)≤1},{x,y}]

Multivariate Optimization

Special case methods

Special case methods solve optimization problems with much larger numbers of variables and constraints than the general methods could handle.

Linear programming

When the objective function and all constraints are linear, global extrema can be computed exactly using the simplex algorithm.

SeedRandom[1];n=10;
A=Table[Random[Integer,{-1000,1000}],{n/2},{n}];
B=Table[Random[Integer,{-1000,1000}],{n/2},{n}];
a=Table[Random[Integer,{-1000,1000}],{n/2}];
b=Table[Random[Integer,{-1000,1000}],{n/2}];
c=Table[Random[Integer,{-1000,1000}],{n}];
X=x/@Range[n];
Minimize[{c.X, And@@Thread[A.X==a]&& And@@Thread[b≤ B.X≤ b+100]},X]

Linear-fractional programming

The objective function is a ratio of linear functions and all constraints are linear.

SeedRandom[1];n=10;
A=Table[Random[Integer,{-1000,1000}],{n/2},{n}];
B=Table[Random[Integer,{-1000,1000}],{n/2},{n}];
a=Table[Random[Integer,{-1000,1000}],{n/2}];
b=Table[Random[Integer,{-1000,1000}],{n/2}];
c=Table[Random[Integer,{-1000,1000}],{n}];
d=Table[Random[Integer,{-1000,1000}],{n}];
X=x/@Range[n];
Maximize[{c.X/d.X, And@@Thread[A.X==a]&& And@@Thread[b≤ B.X≤ b+100]},X]

Integer linear programming

To solve an integer linear programming problem Mathematica first solves the equational constraints, reducing the problem to one containing inequality constraints only. Then it uses lattice reduction techniques to put the inequality system in a simpler form. Finally, it solves the simplified optimization problem using a branch-and-bound method.

SeedRandom[1];n=7;
A=Table[Random[Integer,{-1000,1000}],{3},{n}];
a=Table[Random[Integer,{-1000,1000}],{3}];
B=Table[Random[Integer,{-1000,1000}],{3},{n}];
b=Table[Random[Integer,{-1000,1000}],{3}];
c=Table[Random[Integer,{-1000,1000}],{n}];
X=x/@Range[n];
eqns=And@@Thread[A.X==a];
ineqs=And@@Thread[B.X≤b];
bds=And@@Thread[X≥0]&&Total[X]≤(10^100);
Maximize[{c.X,eqns&&ineqs&&bds&&X∈Integers},X]

Other methods

• Monotonic decomposition
• Piecewise subdivision
• ...

f=SinhIntegral[ScorerHi[LogisticSigmoid[Gudermannian[Erf[ArcTan[ArcSinh[(x^3)+(y^3)]]]]]]];

Maximize[{f,x^4+y^4≤1},{x,y}]//RootReduce

f=Max[Sin[x],Cos[x]];

Minimize[f,x]

Plot[f,{x,-5Pi,5Pi},Epilog->{Red,PointSize[Large],Point[{x,f}/.%[[2]]]}]

Thank you for listening!