Exact Global Constrained Optimization

Adam Strzebonski

Wolfram Research

Abstract

Mathematica functions Minimize and Maximize allow users to compute exact global extrema of functions on sets constrained by systems of equations and inequalities. When the objective function and the constraints are real algebraic functions, methods based on cylindrical algebraic decomposition allow users to always compute global extrema (or extremal values, if they are not attained). If parameters are present, the developmental version of Mathematica can compute the extrema as functions of parameters. In addition to cylindrical algebraic decomposition, in certain special cases, like linear objective and constraints, bounded constraints, or univariate objective and constraints, faster methods can be used. A method based on Lagrange multipliers, used for (sufficiently regular) bounded constraints, can also succeed for some transcendental problems, thanks to the transcendental equation and inequality solving functionality. Mathematica can also solve some optimization problems over the integers. Bounded linear problems are solved by integer linear programming methods. Several other integer optimization problems can be solved by a combination of real optimization methods and integer solution finding. In my talk I will show 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]) consisting of equations and inequalities.
• Domain specifications for the variables (Subscript[x, 1]),  ... ,  (Subscript[x, n]) (either ℝ or ℤ).

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 satifying 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.

Optimization by Cylindrical Algebraic Decomposition

The most general method of exact global optimization is based on the cylindrical algebraic decomposition (CAD) algorithm. It applies when the objective function and the constraints are real algebraic functions. The method allows to always compute global extrema (or extremal values, if the extrema are not attained). If parameters are present, the extrema can be computed as piecewise-algebraic functions of the parameters. A downside of the method is its high, doubly exponential in the number of variables, complexity.

This finds the point on the cubic curve (x^3) + (y^2) - x == (1/4) which is closest to the origin.

Minimize[{(x^2) + (y^2), (x^3) + (y^2) - x == (1/4)}, {x, 
        y}]

The method 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}] // N[#, 20] &

Here the input differs only in the 21-st decimal digit, but the answer is quite different, especially 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}] // N[#, 20] &

Strict Polynomial Inequality Constraints

If there are no parameters, all constraints are strict inequalities and we only need the extremum value, a significantly simpler version of the CAD algorithm can be used. (We can safely make inequality constraints strict if we know  that \(TraditionalForm\`C ⊆ \(int(C)\)\&_\), where C is the solution set of the constraints.) In this case we can disregard lower dimensional parts of the solution set of the constraints. This means that we can use a simpler projection operator, need to construct only full-dimensional cells, and the cell construction phase of CAD does not require algebraic number computations.

This finds the smallest ball centered at the origin which contains the set bounded by the surface (x^4) + 2 (y^4) + 3 (z^4) - x y z == 1. A full Minimize call with the same input did not finish in 10 minutes.

Experimental`Supremum[{(x^2) + (y^2) + (z^2), 
          (x^4) + 2 (y^4) + 3 (z^4) - x y z < 1}, {x, y, z}] // 
      Timing

\(N[\@%[[2]], 20]\)

Linear Optimization

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]

Univariate Optimization

For univariate problems equation and inequality solving methods are used to find the constraint solution set and discontinuity points and zeros of the derivative of the objective function within the set.

Minimize[(x^2) + (2^x), x]

N[%, 20]

Plot[(x^2) + (2^x), {x, -1, 1}]

Here Mathematica recognizes that the objective functions and the constraints are periodic.

Minimize[{Tan[2 x - (π/2)]^2, -(1/2) ≤ 
          Sin[x] ≤ (1/2)}, x]

Optimization by Finding Stationary and Singular Points

Another approach to finding global extrema is to find all the local extrema, using the Lagrange multipliers or the KKT conditions, and pick the smallest (or the greatest). However, for a fully automatic method, there are additional complications. In addition to solving the necessary conditions for local extrema, it needs to check smoothness of the objective function and smoothness and nondegeneracy of the constraints. It also needs to check for extrema at the boundary of the set defined by the constraints and at infinity, if the set is unbounded. This in general requires exact solving of systems of equations and inequalities over the reals, which may lead to CAD computations that are harder than in the optimization-by-CAD algorithm. Currently Mathematica uses Lagrange multipliers only for equational constraints within a bounded box. The method also requires that the number of stationary points and the number of singular points of the constraints be finite. An advantage of this method over the CAD-based algorithm is that it can solve some transcendental problems, as long as they lead to systems of equations that Mathematica can solve.

Here is an example where the minimum is attained at a singular point of the constraints.

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

<< Graphics`ImplicitPlot`;
ImplicitPlot[(y^3) == (x^2), {x, -2, 2}]

The maximum of the same objective function is attained on the boundary of the set defined by the constraints.

Maximize[{y, 
        (y^3) == (x^2) && -2 ≤ x ≤ 
            2 && -2 ≤ y ≤ 2}, {x, y}]

There are no stationary points in this example.

Reduce[
      (y^3) == (x^2) && -2 x λ == 0 && 
        1 + 3 (y^2) λ == 0, {x, y, λ}]

Here is a set of 3-dimensional examples with the same constraints. Depending on the objective function, the maximum is attained at a stationary point of the objective function on the solution set of the constraints, at a stationary point of the restriction of the objective function to the boundary of  the solution set of the constraints, and at the boundary of the boundary of the solution set of the constraints.

Maximize[x + y + z, 
      (x^2) + (y^2) + (z^2) == 9 && -2 ≤ x ≤ 
          2 && -2 ≤ y ≤ 2 && -2 ≤ 
          z ≤ 2, {x, y, z}]

Maximize[x + y + 2 z, 
      (x^2) + (y^2) + (z^2) == 9 && -2 ≤ x ≤ 
          2 && -2 ≤ y ≤ 2 && -2 ≤ 
          z ≤ 2, {x, y, z}]

Maximize[x + 2 y + 2 z, 
      (x^2) + (y^2) + (z^2) == 9 && -2 ≤ x ≤ 
          2 && -2 ≤ y ≤ 2 && -2 ≤ 
          z ≤ 2, {x, y, z}]

<< Graphics`ContourPlot3D`;
ContourPlot3D[
      (x^2) + (y^2) + (z^2) - 9, {x, -2, 2}, {y, -2, 2}, {z, -2, 2}, 
      PlotPoints -> 5]

The Lagrange multiplier method works for some optimization problems involving transcendental functions.

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

N[%,  20]

Optimization over the Integers

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];
A = Table[Random[Integer, {-1000, 1000}], {3}, {7}];
a = Table[Random[Integer, {-1000, 1000}], {3}];
B = 
        Table[Random[
            Integer, {-1000, 
              1000}], {3}, {7}];
b = Table[Random[Integer, {-1000, 1000}], {3}];
c = Table[Random[Integer, {-1000, 1000}], {7}];
X = x /@ Range[7];
eqns = And @@ Thread[A . X == a];
ineqs = And @@ Thread[B . X ≤ b];
bds = 
        And @@ Thread[X ≥ 0] && 
          Total[X] ≤ (10^100);
Minimize[{c . X, eqns && ineqs && bds && X ∈ Integers}, 
      X]

Optimization over the Reals + Integer Solution Finding

Suppose we need to minimize a function f∈ℤ[vars] over the integer solution set of constraints constr. Let m be the minimum of f over the real solution set of constr. If there exists an integer point satisfying f == Ceiling[m]∧constr, then Ceiling[m] is the minumum of f over the integer solution set of constr. Otherwise we try to find an integer solution of f == Ceiling[m] + 1∧constr etc. This heuristic works if we can solve the real optimization problem and all the integer solution finding problems, and we can find an integer solution within a predefined number of attempts. (By default Mathematica tries 10 candidate optium values. This can be changed using the IntegerOptimumCandidates system option.)

Maximize[{x + y, 
        (x^3) + (y^3) ≤ 1000 && (x | y) ∈ 
            Integers}, {x, y}]

Show[{Plot[Root[(x^3) + #(1^3) - 1000 &, 1], {x, -20, 20}, 
          PlotRange -> All, DisplayFunction -> Identity], 
        Graphics[{PointSize[0.03], Red, Point[{6, 9}]}]}, 
      DisplayFunction -> $DisplayFunction]

Customized CAD Algorithm for Optimization

The CAD problem setup

Suppose we need to minimize an algebraic function f over the solutions set of algebraic constraints constr. f and constr depend on variables vars and parameters params. Let y be a new variable. The minimization of f over constr is equivalent to the minimization of the coordinate function y over the semialgebraic set given by y == f∧constr.

If f happens to be a monotonic function of one variable x we do not need to introduce a new variable, we can minimize x instead.

The projection phase of CAD

We project the variables vars first, then the new variable y, and then the parameters params.

If algebraic functions are present, they are replaced with new variables and equations and inequalities satisfied by the new variables are added. The variables replacing algebraic functions are projected first. They also require special handling in the lifting phase of the algorithm.

Projection operator improvements by Hong, McCallum, and Brown can be used here, including the use of equational constraints. Note, that if we needed to introduce a new variable, there is at least one equational constraint, namely y == f.

The lifting phase of CAD

We lift the parameters first, constructing the algebraic function equation and inequality description of the cells. If there are constraints that depend only on parameters and we can determine that constr is identically false over a parameter cell, we don't need to lift this cell further.

When lifting the minimization variable y, we start with the smallest values of y, and proceed (lifting the remaining variables in the depth-first manner) until we find the first cell for which the constraints are satisfied. If this cell corresponds to a root of a projection polynomial in y, the root is the minimum value of f, and the coordinates corresponding to vars of any point in the cell give a point at which the minimum is attained. If parameters are present, we get a parametric description of a point in the cell in terms of roots of polynomials bounding the cell. If there are no parameters we can simply give the sample point used by the CAD algorithm. If the first cell satisfying the constraints corresponds to an interval (r, s), where r is a root of a projection polynomial in y, then r  is the infimum of values of f, and the infimum value is not attained. Finally, if the first cell satisfying the constraints corresponds to an interval (-∞, s), f is unbounded from below.