Elucidating the Geometry of Constrained Nonlinear Optimization With Symbolic Linear Algebra

Frank J. Kampas, Physicist at Large Consulting
frank @ physicistatlarge.com

Abstract

The mathematics behind various nonlinear optimization techniques can be obscured by the details of the numerical methods used to reach the solution.  Using symbolic matrix functions and symbolic solution of the optimization conditions it is possible to make the fundamentals of various techniques more transparent. Four different techniques are described.  Two of them, the Karush- Kuhn-Tucker equations and a Jacobian Minors approach, use the fact that the objective function gradient lies in the space of the constraint functions gradients.  The other two, using a tangent space basis set or a tangent space projection operator, use the fact that the objective function gradient has no component in the constraint tangent space.

Introduction

The condition for a point to be a local minimum or maximum of an objective function subject to constraints is straightforward:  the gradient of the objective function at that point must be orthogonal to any vector in the tangent space of the constraints at that point.  The tangent space can be thought of as a flattened-out version of the surface satisfying the constraints at that point.  Mathematically, any vector in the tangent space is orthogonal to the constraint gradients. Since both the objective function and constraint gradients are orthogonal to the tangent space, it is not surprising that the objective function gradient must be a linear combination of the constraint gradients at a minimum or maximum of the objective function.  The coefficients are the negatives of the Lagrange multipliers, assuming that the Lagrangian is defined as the objective function plus the multipliers times the constraints (Luenberger, p. 300).  So, at a point satisfying the constraints, one has two subspaces:  the tangent space and constraint gradient space.  For that point to be a local minimum or maximum, the gradient of the objective function must be completely in the constraint gradient space, having no component in the tangent space.

Optimization Conditions

In this section, the mathematical conditions for the maximum or minimum are constructed for the four different approaches.  First it is useful to defined some auxiliary functions.

Auxiliary Functions

gradObj calculates the gradient of the objective function with respect to the variables

gradObj[obj_,vars_List]:= D[obj,{vars}]

consconvrule converts constraints to the form expr == 0 or expr <= 0.  This includes constraints of the form a <= expr <= b and a >= expr >= b

consconvrule= {
GreaterEqual[x_,y_]->LessEqual[y-x,0],
Equal[x_,y_]->Equal[x-y,0],
LessEqual[x_,y_]->LessEqual[x-y,0],
LessEqual[lb_,x_,ub_]-> LessEqual[(x-lb)(x-ub),0], 
GreaterEqual[ub_,x_,lb_]-> LessEqual[(x-lb)(x-ub),0]
};

lhsCons uses consconvrule to put constraints the form expr == 0 and expr <= 0 and returns the left-hand side, giving the expr terms, in two separate lists for the equality and inequality constraints.

lhsCons[cons_List]:= Cases[cons /. consconvrule,#[_,0]][[All,1]]& /@ {Equal,LessEqual}

gradeqCons uses lhsCons to find the gradient of the equality constraint expressions with respect to the variables.

gradEqCons[cons_List, vars_List]:= D[lhsCons[cons][[1]],{vars}]

The clearDen function clears the denominators of vectors to avoid problems where the denominator terms are zero

clearDen[l_List]:= l*Times @@ Denominator[l]

addSlackVars converts inequality constraints of the form expr <= 0 to equality constraints of the form expr + (ζ^2)==0, for approaches that do not handle inequality constraints directly.

addSlackVars[obj_, cons_List, vars_List]:=
Module[{eqcons, ineqcons, ζs, newineqcons, newcons,newvars},
{eqcons,ineqcons} = lhsCons[cons];
Ζ = (Subscript[ζ, #]) & /@ Range[Length[ineqcons]];
newineqcons = ineqcons + (Ζ^2);
newcons = Thread[Join[eqcons,newineqcons]==0];
newvars = Join[vars,Ζ];
{obj, newcons, newvars}
]

Karush Kuhn Tucker (KKT) Equations

In this approach the optimality condition is that the gradient of the objective function be a linear combination of the gradients of the equality and inequality constraints, with the coefficients being - (Subscript[λ, i]) for the equality constraints and - (Subscript[μ, j]) for the inequality constraints.  (The negatives signs are used to be consistent with the usual definition of the Lagrangian as being the objective function plus the multipliers times the constraints.)  The inactive inequality constraints are eliminated by equations in which the product of the multiplier and the constraint value are set to zero.  This forces the multiplier to be zero if the inequality constraint is nonzero (inactive).  The condition that the μ multipliers be non-negative was not included, in order to obtain both minima and maxima.  See Luenberger, p. 315.

In 2005 I described solving the KKT equations using Reduce (Kampas).  Reduce is much more powerful now than it was then, so more difficult problems can be handled.

KKTConds[obj_, cons_List, vars_List]:= 
Module[
{eqcons,ineqcons,λs,μs,gradeqs, gradineqs, mults},
{gradeqs, gradineqs} = D[{eqcons, ineqcons} = lhsCons[cons],{vars}];
Λ = (Subscript[λ, #])& /@ Range[Length[eqcons]];
Μ = (Subscript[μ, #])& /@ Range[Length[ineqcons]];
{And @@ Join[
Thread[gradObj[obj,vars] == -Λ.gradeqs - Μ.gradineqs],
Thread[Μ*ineqcons==0]], 
obj, cons, Join[vars,Λ,Μ]}
]

Jacobian Minors Approach

As described in the documentation for the Minors function, Jacobian Minors are zero where a function is singular.  The combination of the objective function and constraints is singular at a maximum or minimum, since they cannot be inverted.  Therefore the Jacobian of the system will be 0 (Edwards, p. 133).  Another way to view this is noting that the determinant of a matrix with linearly dependent rows or columns is zero, so that the Jacobian Minors condition finds points where the objective and constraint gradients are linearly dependent.

jacMinors constructs the minors of the combined objective and constraint gradients.

jacMinors[obj_, cons_List, vars_List]:= Minors[Join[{gradObj[obj,vars]},gradEqCons[cons,vars]],1+Length[cons]]

jacMinorsConds sets the result of jacMinors equal to 0.  The function has the condition that it only applies to equality constraints.  It returns the optimality conditions, the objective function, constraints, and variables, for use by a function defined later which computes the solution to the optimality conditions and the constraints and computes the value of the objective function.

jacMinorsConds[obj_,cons:{HoldPattern[Equal[__]]...}, vars_List]:= {And @@ Thread[Flatten @ jacMinors[obj,cons,vars] ==0],obj,cons,vars}

Tangent Space Approach

tanSpaceBasis constructs a basis for the constraint tangent space, using the NullSpace function, as the tangent space is, by definition, the null space of the constraint gradient matrix.  The NullSpace function uses row reduction resulting in vectors which have a 1 as one component and terms in the denominator of the other components, so the clearDen function is used to avoid problems that would otherwise occur where denominators equal 0.

tanSpaceBasis[cons_List, vars_List]:= clearDen /@ NullSpace[gradEqCons[cons,vars]]

The tanSpaceConds function sets the inner producst of the tangent space basis vectors with the objective function gradient equal to zero and returns those conditions along with the objective function, constraints and variables.

tanSpaceConds[obj_,cons:{HoldPattern[Equal[__]]...}, vars_List]:=  
{And @@Thread[tanSpaceBasis[cons,vars].gradObj[obj,vars] == 0], obj, cons, vars}

Projection Matrix Approach

In this approach, a projection matrix onto the tangent space is constructed using the pseudoinverse of the constraint gradients (Luenberger, p. 340).  The clearDen function is also used here to simplify the expressions.

pseudoInvConsGrad[cons_List, vars_List]:=  Simplify @ Assuming[vars ∈Reals,Refine @  PseudoInverse[gradEqCons[cons,vars]]]

projMat[cons_List, vars_List]:= Simplify[clearDen /@(IdentityMatrix[Length[vars]] - pseudoInvConsGrad[cons,vars].gradEqCons[cons,vars])]

The optimality condition is that the projection of the objective function gradient onto the tangent space is zero.

projMatConds[obj_, cons:{HoldPattern[Equal[__]]...}, vars_List]:= 
{And @@ Thread[projMat[cons,vars].gradObj[obj,vars]==0],obj,cons,vars}

Slack Variables For Inequality Constraints

Only the KKT equations handle inequality constraints directly.  The other three approaches can be extended to problems with inequality constraints by the use of slack variables.

jacMinorsConds[obj_,cons_List, vars_List]:=jacMinorsConds[Sequence @@ addSlackVars[obj,cons,vars]]

tanSpaceConds[obj_,cons_List,vars_List]:= tanSpaceConds[Sequence @@ addSlackVars[obj,cons,vars]]

projMatConds[obj_,cons_List, vars_List]:=projMatConds[Sequence @@ addSlackVars[obj,cons,vars]]

Examples of the Application of the Optimality Conditions

Some functions need to be defined to combine the optimality conditions with the constraints and calculate the value of the objective function.

Auxiliary Functions

getResfromConds takes the results of the various functions and uses Reduce to solve the equations and then uses Reduce again to find the value of the objective function at the various solutions.  Solutions with negative values of the slack variables are eliminated to avoid duplication.

getResfromConds[{conds_, obj_, cons_,vars_}]:=
Module[{r1, Ζ= Cases[vars,Subscript[ζ,_]]},
r1 = Reduce[conds && And @@ cons && And @@ Thread[Ζ ≥ 0], vars, Reals];
Reduce[r1 && objval == obj,Join[{objval},vars],Reals,Backsubstitution -> True, Cubics -> True, Quartics -> True]
]

getResbyMeth returns the time required for the calculation, the optimality conditions and the results from getResfromConds

getResbyMeth[method_, obj_, cons_List, vars_List]:=
Module[{mres = method[obj,cons,vars],gc, l},
{t,gc} = AbsoluteTiming @ getResfromConds[mres];
TableForm @ Flatten[{method,ToString[t]<> " seconds",mres[[1]]/. And -> List,gc/. Or -> List}]
]

testall solves a problem with all the methods and displays the results for each.

testall[obj_,cons_List, vars_List]:= TableForm[getResbyMeth[#,obj,cons,vars]& /@ {KKTConds,jacMinorsConds,tanSpaceConds,projMatConds },TableSpacing -> {3,1}]

Tests

In this section results from the built-in functions Minimize and Maximize are presented for various test problems and compared to the four methods.

Simple Problem:  Optimize x + y subject to (x^2)+(y^2)= 1

The results from the built-in functions Minimize and Maximize:

#[{x+y,(x^2)+(y^2)==1},{x,y}]& /@ {Minimize,Maximize}//ColumnForm

All four methods give identical results. The KKT equations use a Lagrange multiplier (Subscript[λ, 1]); the other approaches use no auxiliary variables.

testall[x+y,{(x^2)+(y^2)==1},{x,y}]

Constraint Gradients Linearly Dependent at Solution
 Optimize y subject to (x^6) = z and (y^3)= z

#[{y,(x^6)-z==0, (y^3)-z==0},{x,y,z}]& /@ {Minimize,Maximize}//ColumnForm

In this problem, the gradients of the constraints are linearly dependent at the solution, causing the KKT equations to fail.

gradEqCons[{(x^6)==z, (y^3)==z},{x,y,z}]

% /. {x->0,y->0,z->0}

testall[y,{(x^6)==z,(y^3)==z},{x,y,z}]

Constraint Gradients Linearly Dependent Everywhere
 Optimize x+y subject to (x^2)+(y^2)= 1 and 2(x^2)+2(y^2)= 2

#[{x+y,(x^2)+(y^2)==1, 2 (x^2)+ 2 (y^2)==2},{x,y}]& /@ {Minimize,Maximize}//ColumnForm

In this problem, the constraint gradients are linearly dependent everywhere, causing the Jacobian Minors approach to fail.  The KKT approach works, with the Lagrange multipliers being linearly dependent.

testall[x+y,{(x^2)+(y^2)==1, 2 (x^2)+ 2 (y^2)==2},{x,y}]

Simple Inequality Constraint 
 Optimize x + y subject to (x^2)+(y^2)<=  1

Here the first problem has been modified by changing the equality constraint to an inequality constraint.

#[{x+y,(x^2)+(y^2)≤ 1},{x,y}]& /@ {Minimize,Maximize}//ColumnForm

Results are identical for all four methods although the optimality conditions are different.

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

Simple Inequality Constraint Problem Where Slack Variable is NonZero at Minimum
 Optimize (x - (((1/2)))^2) + (y^2) subject to (x^2)+(y^2)<=  1

In this problem, the objective function has been modified so that the minimum occurs where the inequality constraint is inactive.

#[{(((x-(1/2)))^2)+(y^2),(x^2)+(y^2)≤ 1},{x,y}]& /@ {Minimize,Maximize}//ColumnForm

The multiplier (Subscript[μ, 1  ])is 0 in the solution from the KKT conditions for the minimum and the slack variable is zero in the other two approaches for the minimum.

testall[(((x-(1/2)))^2)+(y^2),{(x^2)+(y^2)≤ 1},{x,y}]

The point at {1,0} with objective function value of (1/4)is a local maximum.

Plot3D[(((x-(1/2)))^2)+(y^2),{x,-1.5,1.5},{y,-1.5,1.5},RegionFunction -> Function[{x,y},(x^2)+(y^2)≤ 1]]

Problem where the solution lies along a line
 Optimize x + y subject to x + y <= 1 and (x^2) + (y^2) <= 1

#[{x+y,(x^2)+(y^2) ≤ 1, x+y ≤ 1},{x,y}]& /@ {Minimize,Maximize}//ColumnForm

In this problem the maximum actually lies along a line from {0,1} to {1,0}.  
The tangent space conditions give an incorrect result -(Sqrt[2])<objval<1&&x==(objval/2)&&y==(objval/2)&&(Subscript[ζ, 1])==((Sqrt[2-(objval^2)])/(Sqrt[2]))&&(Subscript[ζ, 2])==(Sqrt[1-objval])
along with the correct results.  This may be due to the fact that the objective function and one of the constraints have the same gradient.

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

Problem With Boundary Constraints
 Optimize (x^2)-(y^2) subject to Cos[x-y] ≥ (1/2) and x and y between -5 and 5

There are a lot of local minima and maxima.

Plot3D[(x^2)-(y^2),{x,-5,5},{y,-5,5},RegionFunction -> Function[{x,y}, Cos[x-y] ≥ (1/2)]]

The built-in symbolic functions do not do this problem

#[{(x^2)-(y^2),Cos[x-y]≥ (1/2), -5≤x≤5,-5≤y<=5},{x,y}]& /@ {Minimize,Maximize}//ColumnForm

but it can be done numerically.

#[{(x^2)-(y^2),Cos[x-y]≥ (1/2), -5≤x≤5,-5≤y<=5},{x,y}, Method -> "DifferentialEvolution"]& /@ {NMinimize,NMaximize}//ColumnForm

Only the result of the KKT equations is shown here for brevity.  The Jacobian Minors and Projection Matrix methods give similar results but take longer.
The result obtained agree with the results of NMinimize and NMaximize.

getResbyMeth[KKTConds,(x^2)-(y^2),{Cos[x-y]≥ (1/2),-5≤x≤5,-5≤y<=5},{x,y}]

TableForm @ N[%]

Conclusions

The geometry and mathematics of constrained nonlinear optimization can be illustrated by using a symbolic linear algebra to illustrate various approaches to the problem, based both in the constrant tangent space and in the constraint gradient space.

References

Edwards, Harold M. “Advanced Calculus:  A Differential Forms Approach”, Birkhauser, 1969

Kampas, Frank J. “Mathematica Journal, vol. 9, Issue 4”, 2005, http://www.mathematica-journal.com/issue/v9i4/contents/Tricks9-4/Tricks9-4_2.html

Luenberger, David G., “Linear and Nonlinear Programming”, 2nd edition, Addison-Wesley, 1989