Pushing the Math Frontier

Breakthroughs in Symbolic Optimization

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 Exact Optimization

General Analytic Method

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

Methods Using Periodicity

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

Methods Using Function Property Information

Use information about properties of built-in mathematical functions, including domain, singularities, monotonicity, limits at domain endpoints, local and global minima and more.

Multivariate Exact 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}]

Analytic Methods

Find all the local extrema, using the Lagrange multipliers or the Karush–Kuhn–Tucker (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}]

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]

Convex optimization (new)

Multiparametric linear programming (new)

Multiparametric quadratic programming (new)

Exact Convex Optimization

KKT Conditions for Convex Optimization Problems

Theorem (sufficiency of KKT conditions)

Let f,(Subscript[g, 1]),…,(Subscript[g, m]) be continuously differentiable convex functions over (ℝ^n), let (Subscript[h, 1]),…,(Subscript[h, p]) be affine functions, and let

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

Suppose 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]) may not be polynomial, so you cannot solve the KKT conditions directly using symbolic methods. However, you 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.

Exact Solutions of Square Systems of Equations

Root[{{(Subscript[f, 1]),…,(Subscript[f, n])},{(Subscript[c, 1]),…,(Subscript[c, n])}},j] represents the 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 the 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]

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

N[%,100]

A similar procedure is now used used by FindInstance to find solutions of transcendental systems:

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}]

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]]}]]

Multiparametric Linear Programming

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}}}

Multiparametric Quadratic Programming

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 Version 12.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 five 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;

Version 12 .2 does not finish in five minutes:

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

Version 12 .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!