Quadratic Programs

Code

Test Matrices etc.

Clear[RandomSym, RandomVec, RandomMat]
RandomSym[n_, Distribution_:UniformDistribution[{1,2}]]:=Module[
{Q,Dia},
Dia=SparseArray[Band[{1,1}]:>RandomReal[Distribution],{n,n}];
Q=QRDecomposition[RandomReal[{-1,1},{n,n}]][[1]];
Q.Dia.Transpose[Q]
]
RandomVec[n_,Distribution_:UniformDistribution[{-1,1}]]:= RandomReal[Distribution,n]
RandomMat[{m_,n_},Distribution_:UniformDistribution[{-1,1}]]:= RandomReal[Distribution,{m,n}]

n=1200;
{
Histogram[Eigenvalues[RandomSym[n,NormalDistribution[2,0.1]]]],
Histogram[Eigenvalues[RandomSym[n]]]
}

KKTQuadratic Solver: Increment version

The following takes in matrices and vectors for the Quadratic program and solves the linear system of KKT equations.  It assumes that A has full rank and that G is symmetric.  It gives the only KKT point x+p if G is positive definate.  It is the increment version 16.5 with

KKTSolve[x_,{G_,c_},{A_,b_}]:= Module[
{KKTMat,KKTRHS,pλStar,n=Length[c]},
If[Length[b]==0, 
{-LinearSolve[G,c+G.x],{}},
KKTMat=ArrayFlatten[{{G, Transpose[A]},{A,0}}];
KKTRHS=Join[c+G.x, A.x-b];
pλStar=LinearSolve[KKTMat, KKTRHS];
{-Take[pλStar,n],Drop[pλStar,n]}
]
]

A 2D test.

G= {
{1.0,0},
{0,3.0}
};
c= {1.0,2.0};
A= {{1.0,1.0}};
b = {1.2};
x={1,-1};
{p,λStar}=KKTSolve[x,{G,c},{A,b}]
Show[
ContourPlot[
(1/2){x1,x2}.G.{x1,x2} + c.{x1,x2},
{x1, -2,2},
{x2,-2,2},
ColorFunction->Hue,
Contours->120
],
ContourPlot[A.{x1,x2}-b==0,
{x1, -2,2},
{x2,-2,2},
ContourStyle->Thickness[0.01]
],
Graphics[{Red, Arrow[{x,x+p}],PointSize[0.02],Blue,Point[x],Point[x+p]}]
]

G= {
{1.0,0},
{0,3.0}
};
c= {1.0,2.0};
A= {{}};
b = {};
x={1,-1};
{p,λStar}=KKTSolve[x,{G,c},{A,b}]
Show[
ContourPlot[
(1/2){x1,x2}.G.{x1,x2} + c.{x1,x2},
{x1, -2,2},
{x2,-2,2},
ColorFunction->Hue,
Contours->120
],
Graphics[{Black, Arrow[{x,x+p}],PointSize[0.02],Blue,Point[x],Point[x+p]}]
]

Sequential Quadratic Programs.

Notation etc.

Min (1/2)x.G.x+c.x satisfying (Subscript[a, i]).x=(Subscript[b, i]) for i=1(Subscript[…m, e])  and (Subscript[a, i+(Subscript[m, e])]).x≥(Subscript[b, i+(Subscript[m, e])]) for i=1(Subscript[…m, i]). 
where G is symmetric.

The Lagrangian is 
	ℒ(x, λ) = 1/2 x.G.x + x.c - (∑_(i=1)^((Subscript[m, e])+(Subscript[m, i])))(Subscript[λ, i])(x.(Subscript[a, i])-(Subscript[b, i]))
The active set at a feasible point x is 
	𝒜(x)={i: 1≤i≤(Subscript[m, i])+(Subscript[m, e]) and x.(Subscript[a, i])-(Subscript[b, i])=0}
The KKT conditions for a minimizer x are that there is a Lagrange multiplier vector λ satisfying
	({{G.x+c-A.λ, =, 0,  }, {(Subscript[a, i]).x, =, (Subscript[b, i]), for i∈𝒜(x)}, {(Subscript[a, i]).x, ≥, (Subscript[b, i]),        for i∈ℐ\𝒜(x)}, {(Subscript[λ, i]), ≥, 0,          for i∈ℐ⋂𝒜(x)}})
The key to progress for QPs is this system.

Equality constrained problems

For positive definite QPs with equality constraints you simply solve the linear system and you get a single KKT point which is the global minima. The code for this is in two equivalent flavors: one giving the optimal (SuperStar[x]) (eq 16.4) and the other giving the shift p (eq 16.5) required to get to the optimal (SuperStar[x]) from some initial condition x i.e. (SuperStar[x])=x+p.    Both are easy to code as long as we just use the linear algebra in matlab or Mathematica.   The function KKTSolve[x_, {G_, c_}, {A_, b_}] implements algorithm 16.5.

n=2;
m=1;
G= RandomSym[n];
c= RandomVec[n];
A=RandomMat[{m,n}];
b = RandomVec[m];
x=RandomVec[n];
{p,λStar}=KKTSolve[x,{G,c},{A,b}];
M=1.5;
Show[
ContourPlot[
(1/2){x1,x2}.G.{x1,x2} + c.{x1,x2},
{x1, -M,M},
{x2,-M,M},
ColorFunction->Hue,
Contours->120
],
ContourPlot[A.{x1,x2}-b==0,
{x1, -M,M},
{x2,-M,M},
ContourStyle->Thickness[0.01]
],
Graphics[{ Arrow[{x,x+p}],PointSize[0.02],Blue,Point[x],Point[x+p]}]
]

Inequality constrained problems

For positive definite QPs with equality and inequality constraints you can use the active set algorithm 16.3 on p472.  The idea is simple add some (includes none and all) of the inequality constraints to the equality constraints to form a working set of constraints, solve the equailty constrained problem with the working set, step as far as you can (and remain feasible) in the direction specified by the solution of the working set problem and adjust the working set. You adjust the working set by adding the constraint (or possibly constraints) that blocked you from going further and/or deleting constraints that do not have the correct sign of multiplier.

We are looking at primal algorithms that always use feasible points. So we will assume we have a feasible point to start with.

DeactivateConstraints[w_,λ_]:= Module[
{j},
(* Empty constraint special case *) 
If[w=={},Return[w]];
(* Deletes Active Constraint with the most negative multiplier *) 
j=Ordering[λ,1][[1]];
If[λ[[j]]>0, Return[w]];
Drop[w,{j}]
]

ActivateConstraints[w_,{x_,p_},{AI_,bI_}]:= Module[
{wComp,mI=Length[bI],αs,jBlock,α=1,Tol=(10^(-14))},
(* Select all inactive inequality constraints *)wComp=Complement[Range[mI],w];
(* Select only downhill constraints *)
wComp=Select[ wComp, AI[[#]].p<0&];
(* Deal with no constraints to test *)
If[wComp=={},Return[{x+ α p,w}]];
(* compute αs that would violate constraints *)
αs=((bI[[wComp]]-AI[[wComp]].x)/(AI[[wComp]].p));
jBlock=Ordering[αs,1][[1]];
(* Constraint blocks only if the smallest α is <1 *)
α = Min[1,αs[[jBlock]]];
If[α<1,
(* return modified step with additional constraint if blocking*)
{x+α p,Append[w,wComp[[jBlock]]] },
(* return full step with no additional constraint if not blocking*)
{x+α p,w}
]
]
ActiveSetQP[{x_,w_},{G_,c_},{{AE_,bE_},{AI_,bI_}}]:= Module[
{A,b,p,λ,Tol=(10^(-6))},
{A,b}={Join[AE,AI[[w]]],Join[bE,bI[[w]]]};
{p,λ}=KKTSolve[x,{G,c},{A,b}];
(*  *)
If[ Norm[p]<Tol,
{x,DeactivateConstraints[w,λ]},
ActivateConstraints[w,{x,p},{AI,bI}]
]
]

Now I need to check which of the λs have the correct sign.  I also need to check which of the constraints the ps are consistent with

n=2;
mE=0;
mI=2;
G= RandomSym[n];
c= RandomVec[n];
AE=RandomMat[{mE,n}];
AI=RandomMat[{mI,n}];
bE = RandomVec[mE];
bI = RandomVec[mI];

Clear[x1,x2]
M=1.5;
x0={x1,x2}/.FindInstance[
{Min[AI.{x1,x2}-bI]≥0.2,-M<x1<M, -M<x2<M},{x1,x2}][[1]]

w0={};
Data=FixedPointList[
ActiveSetQP[#,{G,c},{{AE,bE},{AI,bI}}]&,
 {x0,w0}]
Show[
ContourPlot[
(1/2){x1,x2}.G.{x1,x2} + c.{x1,x2},
{x1, -M,M},
{x2,-M,M},
ColorFunction->ColorData["TemperatureMap"],
Contours->120,
RegionFunction->Function[{x1,x2,f},Min[AI.{x1,x2}-bI]≥0]
],
Graphics[{ 
{Red,Thickness[0.01],Arrow[Data[[All,1]]]},
{Cyan,PointSize[0.02],Map[Point,Data[[All,1]]]},
{Red, Point[Data[[-1,1]]]}
}],
PlotRange->All
]

Just to demonstrate that this can be very slow.

n=123;
{AI,bI}= Transpose[Table[{{Cos[θ],Sin[θ]},-1},{θ,0.0,(1-(1/n))2π ,2π/n}]];
G = {{3.4,1.2},{1.2,1.4}};
c={3.1,4.2};
M=1.5;
{x0,w0}={{0.99,0},{}};
Data=FixedPointList[
ActiveSetQP[#,{G,c},{{AE,bE},{AI,bI}}]&,
 {x0,w0}]
Show[
ContourPlot[
(1/2){x1,x2}.G.{x1,x2} + c.{x1,x2},
{x1, -M,M},
{x2,-M,M},
ColorFunction->ColorData["TemperatureMap"],
Contours->120,
RegionFunction->Function[{x1,x2,f},Min[AI.{x1,x2}-bI]≥0]
],
Epilog->{
{Red,Thickness[0.005], Arrow[Data[[All,1]]]},
{Green,PointSize[0.01], Map[Point,Data[[All,1]]]}
}
]

Interior Point Methods for Quadratic Programs

We are going to implement the algorithm 16.4 on p484.  This algorithm needs a feasible initial point along with the slack varaiables and the multipliers.

KKTSolve1658[{x_,y_,λ_},{σ_,μ_},{{G_,c_},{AI_,bI_}}]:=Module[
{KKTMat, KKTRHS,m=Length[bI],n=Length[c],e,xyλ},
e=ConstantArray[1,{m}];
KKTMat= ArrayFlatten[(({{G, 0, -Transpose[AI]}, {AI, -IdentityMatrix[m], 0}, {0, DiagonalMatrix[λ], DiagonalMatrix[y]}}))];
KKTRHS=Join[-(G.x-Transpose[AI].λ+c),-(AI.x-y-bI),-λ y + σ μ e];
xyλ=LinearSolve[KKTMat,KKTRHS];
{Take[xyλ,n],Take[xyλ,{n+1,n+m}],Take[xyλ,{n+m+1,n+2m}]}
]

αBlockingComputation[{y_,λ_},{Δy_,Δλ_}]:=Module[
{αs},
αs=Join[-(y/(Δy)),-(λ/(Δλ)),{1}];
αs=Select[αs,0<#≤1&];
Min[αs]
]

InteriorPCStepQP[{x_,y_,λ_}, {τ_,{G_,c_},{AI_,bI_}}]:= Module[
{μ,KKTMat,αAff,μAff,σ,m=Length[bI],Δx,Δy,Δλ,αHat,xNew,yNew,λNew},
{σ,μ}={0,((y.λ)/m)};
{Δx,Δy,Δλ}=KKTSolve1658[{x,y,λ},{σ,μ},{{G,c},{AI,bI}}];
αAff= αBlockingComputation[{y,λ},{Δy,Δλ}];
μAff=(((y+αAff Δy).(λ+αAff Δλ))/m);
σ=(((((μAff)/μ)))^3);
{Δx,Δy,Δλ}=KKTSolve1658[{x,y,λ},{σ,μ},{{G,c},{AI,bI}}];
αHat=αBlockingComputation[{τ y,τ λ},{Δy,Δλ}]; (* see 16.66 *)
{x,y,λ}+αHat {Δx,Δy, Δλ}
]

n=2;
mI=2;
G= RandomSym[n];
c= RandomVec[n];
AI=RandomMat[{mI,n}];
bI = RandomVec[mI];

Clear[x1,x2]
M=1.5;
x0={x1,x2}/.FindInstance[
{Min[AI.{x1,x2}-bI]≥0.2,-M<x1<M, -M<x2<M},{x1,x2}][[1]]
y0=AI.x0-bI;
λ0=ConstantArray[1,{mI}];

τ=0.99;
{MaxIter,Tol}={100,(10^(-5))};
Data=FixedPointList[
InteriorPCStepQP[#,{τ,{G,c},{AI,bI}}]&, {x0,y0,λ0},
MaxIter,
SameTest->Function[{xyλ1,xyλ2},Norm[xyλ1[[1]]-xyλ2[[1]]]<Tol]
];
ContourPlot[
(1/2){x1,x2}.G.{x1,x2} + c.{x1,x2},
{x1, -M,M},
{x2,-M,M},
ColorFunction->ColorData["TemperatureMap"],
Contours->120,
RegionFunction->Function[{x1,x2,f},Min[AI.{x1,x2}-bI]≥0],
Epilog->{ 
{Red,Thickness[0.01],Arrow[Data[[All,1]]]},
{Cyan,PointSize[0.02],Map[Point,Data[[All,1]]]},
{Red, Point[Data[[-1,1]]]}
},
PlotLabel->Length[Data]
]

For our many edged problem.

n=23;
{AI,bI}= Transpose[Table[{{Cos[θ],Sin[θ]},-1},{θ,0.0,(1-(1/n))2π ,2π/n}]];
G = {{3.4,1.2},{1.2,1.4}};
c={3.1,4.2};
M=1.5;
x0={-0.2,0.4};
y0=AI.x0-bI;
λ0=ConstantArray[1,{Length[bI]}];
τ=0.5;
{MaxIter,Tol}={100,(10^(-5))};
Data=FixedPointList[
InteriorPCStepQP[#,{τ,{G,c},{AI,bI}}]&, {x0,y0,λ0},
MaxIter,
SameTest->Function[{xyλ1,xyλ2},Norm[xyλ1[[1]]-xyλ2[[1]]]<Tol]
];
ContourPlot[
(1/2){x1,x2}.G.{x1,x2} + c.{x1,x2},
{x1, -M,M},
{x2,-M,M},
ColorFunction->ColorData["TemperatureMap"],
Contours->120,
RegionFunction->Function[{x1,x2,f},Min[AI.{x1,x2}-bI]≥0],
Epilog->{ 
{Red,Thickness[0.01],Arrow[Data[[All,1]]]},
{Cyan,PointSize[0.02],Map[Point,Data[[All,1]]]},
{Red, Point[Data[[-1,1]]]}
},
PlotLabel->Length[Data]
]