Animating Optimization
With the Gradient Projection Method

Frank J. Kampas
Frank @ PhysicistatLarge.com

Abstract

Implementation of the gradient projection optimization method as a differential equation, rather than than a series of discrete steps, enables calculation of optimization paths which are sufficiently smooth for the creation of animations.  This is illustrated by animating the packing of different size circles inside a circumscribing circle.

Introduction

The gradient projection mimization technique (see Luenberger), can be understood in terms of the tangent space of the active constraints at a point.  Since any vector in the tangent space is orthogonal to the gradients of the constraints, a displacement in the direction of a vector in the tangent space does not change the value of the constraints, to first order.  The gradient projection method starts at a feasible point and moves in the direction which is the projection of the negative of the objective function gradient onto the tangent space.  The resulting point has a lower value for the objective function but is no longer feasible.  A second step is then taken, moving perpendicular to the tangent space of the orginal point, to get back to feasibility.  This process is repeated until a local minimum is found. 

In the approach taken here, the gradient projection method is implemented using a differential equation approach, in which the time derivative of a point is equal to the projected negative objective function gradient plus the correction term for moving to feasibility, starting from a non-feasible point.  This is slower than the usual step-wise approach, but results in a smooth path from a non-feasible, non-optimal point to a feasible local minimum, since NDSolve controls the step size.  This path is suitable for animation.  The method is illustrated with an animation of the packing of different size circles inside a circumscribing circle.

Mathematical Approach

Following Luenberger (section 11.4), we consider the problem of minimizing f(x), subject to the equality constraints h(x) = 0, where x is the point in question, f(x) is the objective function and h(x) is the vector of constraints.  The projection matrix (Subscript[P, k]) onto the constraint tangent space at point (Subscript[x, k]) is given by (Luenberger p. 336)

(Subscript[P, k]) = I - Transpose[ ∇h((Subscript[x, k])) ] . Inverse[ ∇h((Subscript[x, k])) . Transpose[ ∇h((Subscript[x, k])) ] ] . ∇h((Subscript[x, k]))

where  ∇h((Subscript[x, k])) is the gradient (Jacobian) matrix of the constraints and I is an identity matrix.  In Mathematica notation, the gradient is given by

∇h((Subscript[x, k])) = ( D[h(x), {x}] ) /. x -> (Subscript[x, k])

The projection matrix onto the constraint tangent space can also be written as

(Subscript[P, k]) = I - PseudoInverse[ ∇h((Subscript[x, k])) ]  .  ∇h((Subscript[x, k]))

since the projection matrix onto the constraint gradient space at point (Subscript[x, k]) is  PseudoInverse[ ∇h((Subscript[x, k])) ]  .  ∇h((Subscript[x, k])) (see Strang, p. 404).  As described by Strang, PseudoInverse[A ]. A is the projection matrix onto the row space of A.  For ∇h((Subscript[x, k])), each row is a gradient of a constraint.  The projection matrix onto the tangent space is the identity matrix minus the projection matrix onto the constraint gradient space since the two spaces are complementary.

To see the equivalence of the PseudoInverse formula to Luenberger’s formula, consider the following example which shows that for a symbolic 2 x 2 matrix, the PseudoInverse[m] == Transpose[m].Inverse[m.Transpose[m]]:

m = {{a,b},{c,d}};

Reduce[Assuming[{a,b,c,d}∈Reals,Refine[PseudoInverse[m]]]== Transpose[m].Inverse[m.Transpose[m]]]

The reader is encouraged to read Strang’s description of projections onto a subspace (p. 209).  The formula used by Luenberger assumes that the constraint gradients are linearly independent, whereas the PseudoInverse approach is more general.

As described previously, the gradient projection technique involves starting at the point (Subscript[x, k]), moving in the direction of the negative objective function gradient projected on the constraint tangent subspace using the projection method to a point (Subscript[y, 1]), and then moving back to feasibility in a series of steps, in a direction orthogonal to the constraint tangent space at point x (not (Subscript[y, j])) to reduce the computations required.  Luenberger’s equation (p. 336) for moving back to feasibility is

(Subscript[y, j+1])= (Subscript[y, j]) - Transpose[ ∇h((Subscript[x, k])) ] . Inverse[ ∇h((Subscript[x, k])) . Transpose[ ∇h((Subscript[x, k])) ] ] . h((Subscript[y, j]))

where   (Subscript[y, j]) and  (Subscript[y, j+1]) are successive points on the path back to feasibility.

As before, this is rewritten as

(Subscript[y, j+1])= (Subscript[y, j]) - PseudoInverse[ ∇h((Subscript[x, k])) ]   . h((Subscript[y, j]))

The differential equation approach starts at an infeasible point and applies the sum of the two corrections at that point to obtain the time derivative of x.  Therefore the correction to improve the objective function does not affect the constraint violations and the correction to reduce the constraint violations does not affect the objective function value, both to first order.  

The correction to improve the objective function is

-  (Subscript[P, k]) . ∇f(x) = - ( I - PseudoInverse[ ∇h(x) ]  .  ∇h(x) ) .  ∇f(x) =

-∇f(x) +  PseudoInverse[ ∇h(x) ]  .  ∇h(x)  .  ∇f(x)

where ∇f(x) is the gradient of the objective function.

The correction to reduce the infeasibility is

- PseudoInverse[ ∇h(x) ]   . h(x)

Therefore the total correction, which becomes ((dx)/dt)is

PseudoInverse[ ∇h(x) ] . ( ∇h(x)  .  ∇f(x)  -  h(x) )  - ∇f(x)

Code

The function cg is a compiled function which takes the gradient of the objective function, the gradient of the constraints and the constraints themselves as arguments and returns the total correction.

cg= Compile[
{
{gradobj,_Real,1},(* gradient of the objective function *)
{gradcons,_Real,2},(* gradient of the constraints *)
{cons,_Real,1} (* constraints *)
},
PseudoInverse[gradcons].(gradcons.gradobj - cons)-gradobj, RuntimeOptions -> "Speed"];

The function cgn will be used to require that the argument for cg be numeric.  This is necessary to insure numeric computation of the Pseudoinverse, which can be calculated symbolically.  For problems with a small number of variables, calculating the PseudoInverse once symbolically is faster than calculating it numerically at every point, but that is not the case for most problems.

cgn[f1_,f2_,f3_][x:{_?NumberQ..}]:= cg[f1[x],f2[x],f3[x]]

The function projgraddeq takes as arguments the objective function, the equality constraints (in the form of the the left-hand side of expr == 0), the variables, the initial point and the “time” ending point of the differential equation which is generated from the inputs.  The correction term is multiplied by 1 + (t^2) to improve convergence.  The discrete points found by NDSolve are returned, using Reap and Sow, rather than the interpolating function generated from them.

projgraddeq[obj_, cons_List, vars_List,inits_List,timelimit_]:=
Module[
{rl,(* A rule to convert the variables to the vector xx *)
gradobj,(* The gradient of the objective function *)
gradcons,(* The gradient of the constraints *)
gradobjx,(* gradobj as a function of xx *)
gradconsx, (* gradcons as a function of xx *)
consx,(* the constraints as a function of xx *)
cgradobj,(* compiled gradobjx *)
cgradcons,(* compiled gradconsx *)
ccons,(* compiled cconsx *)
var,(* the variable for NDSolve *)
res (* the results from the Reap on NDSolve *)
},
Quiet[rl = Thread[vars -> Table[xx[[i]],{i,Length[vars]}]]];
{gradobj,gradcons} = D[{obj,cons},{vars}];
{gradobjx,gradconsx,consx} = {gradobj,gradcons,cons} /. rl;
cgradobj= Compile[{{xx,_Real,1}},Evaluate@gradobjx, RuntimeOptions -> "Speed"];
cgradcons = Compile[{{xx,_Real,1}},Evaluate @ gradconsx, RuntimeOptions -> "Speed"];
ccons = Compile[{{xx,_Real,1}},Evaluate @ consx, RuntimeOptions -> "Speed"]; res = Reap @ NDSolve[{var'[t]== (1+(t^2))*cgn[cgradobj,cgradcons,ccons][var[t]],var[0]==inits},var[t],{t,0,timelimit},MaxSteps -> (10^5), EvaluationMonitor :> Sow @ var[t]] [[1]];
res[[2,1]]
]

A Simple Example

In this example, minimize x + y subject to (x^2)+(y^2)==1 is solved, starting at the point {1.0, 1.5}.

AbsoluteTiming[sln=projgraddeq[x+y,{(x^2)+(y^2)-1},{x,y},{1,1.5},5];]

The path taken shows the effect of projected objective function gradient as well as the movement towards feasibility.  Initially both are in about the same direction.  As feasibility is approaches, the path moves closer and closer to the circle satisfying the constraints as it moves to the minimum.

p1=Show[ContourPlot[(x^2)+(y^2)==1,{x,-2,2},{y,-2,2}],ListPlot[sln, AspectRatio -> 1]]

The last point is close to the constrained minimum.

sln[[-1]]

N @ Minimize[{x+y,(x^2)+(y^2)==1},{x,y}]

This problem is simple enough that the calculation can be done symbolically.

f = x+y;h = {(x^2)+(y^2)-1};vars = {x,y};

{gradf,gradh} = D[{f,h},{vars}]

pi = PseudoInverse[gradh]/. Conjugate[x_] -> x

The correction term for the improving the objective function

fcorr[{x_,y_}] = -gradf + pi.gradh.gradf//Simplify

The correction term for reducing the constraint violation

hcorr[{x_,y_}] = -pi.h

As expected, the two terms are orthogonal

fcorr[{x,y}].hcorr[{x,y}]//Simplify

The correction terms are shown graphically with a red line for the objective function correction and a green line for the constraint violation correction.

fcorra[{x_,y_}]= {Red,Line[{{x,y},{x,y}+fcorr[{x,y}]}]};

hcorra[{x_,y_}]= {Green,Line[{{x,y},{x,y}+hcorr[{x,y}]}]};

fhcorra[{x_,y_}] = Sequence[fcorra[{x,y}],hcorra[{x,y}]];

npts = Length[sln];

corrs = Table[fhcorra[sln[[Floor[r*npts]]]],{r,.05,.6,.05}];

Show[p1,Graphics @ corrs]

Circle Packing

The function circlesv takes the number of circles and the function describing the radii of the circles as arguments.  It returns the objective function (the radius of the circumscribing circle), the constraints, the variables and an initial point, generated randomly.  The inequality constraints in the problem are converted to equality constraints using surplus variables t[i] and t[i,j].  Converting all the constraints to equality constraints works better for the differential equation approach, as compared to an active set approach.  Note that the initial value for the radius of the circumscribing circle is set to 2.0.  The initial values of the surplus variables are set to 0.01, as some of them will end up close to 0.0.

circlessv[n_Integer,func_]:=
Module[
{
r, (* the radii of the circles *)
circrad, (* the radius of the circumscribing circle *)
cons1,(* the separation of the centers of the circles is ≥ to sum of their radii *)
cons2,(* the radius of the circumscribing circles is ≥ the points on the circles most distant from the origin *)
cons,(* all the constraints *)
vars,(* all the variables, including the surplus variables *)
inits (* the randomly generated initial values *)
},
r = func /@ Range[n];
cons1=Flatten [Table[  (Sqrt[(((x[i]-x[j]))^2)+(((y[i]-y[j]))^2)])- (r[[i]]+r[[j]])-((t[i,j])^2),{i,n-1},{j,i+1,n}]];
cons2 = Table[circrad -r[[i]]-   (Sqrt[ ((x[i])^2)+((y[i])^2)])-((t[i])^2),{i,n}];
cons = Join[cons1,cons2];
vars=Join[Flatten[Table[{x[i],y[i]},{i,n}]],Flatten[Table[t[i,j],{i,n-1},{j,i+1,n}]],Table[t[i],{i,n}],{circrad}];
SeedRandom[0];
inits = Join[RandomReal[{-1,1},2*n],Table[0.01,{n+ ((n(n-1))/2)}],{2}];
{circrad,cons,vars,inits}
]

c10sv is the result of circlessv for 10 circles with radii 1, (Sqrt[(1/2)]), (Sqrt[(1/3)]) etc.

c10sv = circlessv[10,((#)^((-1)/2))&];

AbsoluteTiming[sln10 = projgraddeq[Sequence @@ c10sv,15];]

Several thousand circle configurations are generated.

nconfigs = Length[sln10]

plt generates a graphic of the circle configuration from the j’th configuaration.

plt10[j_]:=Show @ Graphics[Join[Table[{Hue[(i/10),1,1,.75],Disk[{x[i],y[i]},(i^((-1)/2))]},{i,10}],{Circle[{0,0},c10sv[[3,-1]]]}]/. Thread[c10sv[[3]] -> sln10[[j]]],PlotRange -> {{-2.3,2.3},{-2.3,2.4}}]

The initial configuration

plt10[1]

The final configuration

plt10[nconfigs]

An animation of the optimization. (Click on the + to reveal the animation menu)

Manipulate[plt10[j],{j,1,nconfigs,1}]

The radius of the circumscribing circle as a function of step number.

ListPlot[sln10[[All,-1]]]

It is also possible to generate an animated GIF

Export["circle_packing10.gif",Table[plt10[j],{j,1,nconfigs,10}]]

Conclusion

Implenting the gradient projection optimization method as a differential equation enables generating solutions with a smooth path, which can be used in animations.

References

“Linear and Nonlinear Programming” (2nd edition) by David G. Luenberger, Addison-Wesley Publishing Company, 1989
“Introduction to Linear Algebra” (4th edition) by Gilbert Strang, Wellesley - Cambridge Press, 2009