Karush-Kuhn-Tucker (KKT) Conditions

Adam Rumpf, 5/7/2018

Introduction

This is a graphical representation of the KKT conditions applied to a 2D constrained optimization problem of the form:

min  f(x,y)
s.t. g(x,y) ≤ 0

In the main Manipulate environment defined below, the contour lines of the objective function f and the constraint function g are shown in gray and blue, respectively, with the feasible set being shaded blue. g is adjustable with the “low”, “mid”, and “high” controls, which move it up and down the y-axis. The solution to this optimization problem is the point ((SuperStar[x]),(SuperStar[y])) within the feasible region that achieves the lowest value of f((SuperStar[x]),(SuperStar[y])), and is shown as a red dot in the upper left figure. Naturally the optimal solution may change as g changes.

The user may click and drag the top left plot to select a candidate solution. The other parts of the figure will display whether or not the chosen point satisfies the various KKT conditions, displaying “✓” if the condition is satisfied and “X” if not. Try moving the point around to watch how this effects the evaluations of the various conditions. Also try coming up with solutions that satisfy some conditions but not others, which should help to explain why all of them are necessary for optimality.

To explain the different parts of the display:

The top left shows the chosen point and the feasible set, and corresponds to the primal feasibility condition, which simply states that the point must lie inside the feasible set. The chosen point is also accompanied by a black and a blue vector which show the gradients of f and g, respectively.

The top right shows an expanded view of the ∇f and ∇g vectors from the top left, and corresponds to the stationarity condition, which is that ∇f=-μ ∇g for some constant μ. This is satisfied either if ∇f=0 (in which case μ=0) or if ∇f and ∇g are collinear (in which case μ is the scaling factor of ∇g relative to ∇f).

The bottom left shows a combination of the dual feasibility and complementary slackness constraints. Both concern the constant μ defined in the stationarity condition, which only exists if stationarity is satisfied. If stationarity is not satisfied then μ is undefined and nothing is displayed. Otherwise, the display shows the point (μ,g) on the μ versus g axis. Complementary slackness states that μ g=0, which is equivalent to requiring that at least one out of μ and g be zero. Dual feasibility is μ≥0 while primal feasibility is g≤0, and combining all of this tells us that the optimal solution must have (μ,g) somewhere on either the positive μ-axis or the negative g-axis.

The bottom right does not correspond to any KKT condition, and simply prints the objective values of the global optimum and of the current solution.

Code

Initialization

(* objective *)
f[x_,y_]:=(((x-0.5))^2)+0.5(y^4)-0.5(((y+0.5))^2)

(* unconstrained min *)
fmin={x,y}/.NMinimize[f[x,y],{x,y}][[2]];

(* boundary function with adjustable offset *)
g[x_,y_,c_]:=y-c+(x^2)

(* constraint offset to put min on boundary *)
co=c/.NSolve[g[fmin[[1]],fmin[[2]],c]==0,c][[1,1]];

(* global min for different boundary offsets *)
sol[c_]:={x,y}/.NMinimize[{f[x,y],g[x,y,c]==0},{x,y}][[2]]

(* global min for specific boundary offsets *)
sol1=sol[0];

sol2=sol[co];

sol3=sol[2];

(* contour levels necessary to line up with global mins *)
cont1=f[sol1[[1]],sol1[[2]]];

cont3=f[sol3[[1]],sol3[[2]]];

contdiff=cont1-cont3;

(* calculate smallest angle between two vectors *)
angle[u_,v_]:=N[If[Norm[u]≠0&&Norm[v]≠0,ArcCos[Abs[((u.v)/(Norm[u]Norm[v]))]],0]]

Demonstration

Manipulate[Module[{flag=True,reddot,boundary,snap,fgrad,ggrad,fgrado,ggrado,μ=0,gval,tol=0.05,tempsnap,temp,plot1,plot2,plot3,plot4,pfeas=False,static=False},
reddot=Switch[c,0,sol1,co,sol2,2,If[flag==False,sol3,fmin]];
(* snap to boundary if barely outside it *)
If[g[pt[[1]],pt[[2]],c]≤0,
snap=pt;
pfeas=True,
tempsnap=RegionNearest[ImplicitRegion[g[x,y,c]==0,{x,y}],pt];
If[EuclideanDistance[tempsnap,pt]≤0.15,
snap=tempsnap;
pfeas=True,
snap=pt;
];
];
(* snap to exact solution if near it *)
If[EuclideanDistance[snap,reddot]≤0.15,
snap=reddot;
];
fgrado=fgrad=Grad[f[x,y],{x,y}]/.{x->snap[[1]],y->snap[[2]]};
pt=snap;
(* normalize ∇f, unless it is zero *)
If[Norm[fgrad]>0.1,
fgrad/=3Norm[fgrad],
fgrad={0,0};
static=True;
];
(* normalize ∇g, unless it is zero *)
ggrado=ggrad=Grad[g[x,y,c],{x,y}]/.{x->snap[[1]],y->snap[[2]]};
If[Norm[ggrad]>0.1,
ggrad/=2Norm[ggrad],
ggrad={0,0};
];
(* draw feasible curve/region *)
If[flag==False,
boundary=ContourPlot[g[x,y,c]==0,{x,-4,2},{y,-4,2},ContourStyle->{Thick,Blue}],
boundary={RegionPlot[g[x,y,c]≤0,{x,-4,2},{y,-4,2},BoundaryStyle->{Thick,Blue}],
ContourPlot[g[x,y,c],{x,-4,2},{y,-4,2},ContourShading->None,ContourStyle->Blue,Contours->Table[k,{k,-0.2,-10,-0.2}]]};
];
(* calculate μ and gval *)
gval=g[snap[[1]],snap[[2]],c];
If[Norm[fgrad]==0,
μ=0,
If[angle[fgrad,ggrad]≤tol,
static=True;
temp=1;
If[Abs[ggrad[[2]]]>Abs[ggrad[[1]]],
temp=2;
];
μ=-((fgrado[[temp]])/(ggrado[[temp]])),
μ=Infinity;
];
];
(* display combined plot *)
plot1=Show[ContourPlot[f[x,y],{x,-1,2},{y,-1,2},Contours->Table[cont3+contdiff k,{k,0,20}],ContourShading->None,ContourStyle->LightGray],boundary,Graphics[{LightGray,PointSize[Large],Point[fmin]}],Graphics[{Red,PointSize[Large],Point[reddot]}],Graphics[{Blue,Arrow[{snap,snap+ggrad}],Black,Arrow[{snap,snap+fgrad}],PointSize[Large],Point[snap]}],Graphics[{If[pfeas==True,Text[Style["✓",Darker@Green,14],{1.75,1.75}],Text[Style["Χ",Red],{1.75,1.75}]]}],PlotLabel->"primal feasibility"];
plot2=Show[Graphics[{Thick,Blue,Arrow[{{0,0},If[Max[Abs[ggrado]]≤3,ggrado,3(ggrado/(Max[Abs[ggrado]]))]}],Black,Arrow[{{0,0},If[Max[Abs[fgrado]]≤3,fgrado,3(fgrado/(Max[Abs[fgrado]]))]}]}],Graphics[If[static==True,{Text[Style["✓",Darker@Green,14],{2.5,2.5}],If[Norm[fgrad]==0,Text[Row[{"∇",Style["f",Italic]," = 0"}],{2.5,2.0}],Text[Row[{"∇",Style["f",Italic]," = "<>ToString[PaddedForm[-μ,{2,1}]]<>"∇",Style["g",Italic]}],{3.0,2.0},{1,0}]]},{Text[Style["Χ",Red],{2.5,2.5}]}]],Axes->True,PlotRange->{{-3,3},{-3,3}},AspectRatio->1,PlotLabel->"stationarity",AxesStyle->Gray];
plot3=Show[ListPlot[{},PlotRange->{{-3,3},{-3,3}},AspectRatio->1,AxesLabel->{Row[{Style["g",Italic],"(",Style["x",Italic],",y)"}],"μ"},PlotLabel->"dual feasibility and\ncomplementary slackness",AxesStyle->If[μ==Infinity,LightGray,Gray],TicksStyle->If[μ==Infinity,LightGray,Gray]],Graphics[{Thickness[0.03],Arrowheads[0.1],If[μ==Infinity,LightGray,Gray],Arrow[{{0,0},{-3,0}}],Arrow[{{0,0},{0,3}}],PointSize[0.03],Point[{0,0}],Black,PointSize[0.075],Point[{gval,Min[μ,100]}]}],Graphics[If[μ==Infinity,{Text[Style["Χ",Red],{2.5,2.5}],Text["μ DNE",{3.0,2.0},{1,0}]},If[μ<0,{Text[Style["Χ",Red],{2.5,2.5}],Text["μ < 0",{3.0,2.0},{1,0}]},
If[gval>0.01,
{Text[Style["X",Red],{2.5,2.5}],Text["g > 0",{3.0,2.0},{1,0}]},
If[μ≠0&&Abs[gval]>0.01,
{Text[Style["X",Red],{2.5,2.5}],Text["μg ≠ 0",{3.0,2.0},{1,0}]},
{Text[Style["✓",Darker@Green,14],{2.5,2.5}]}
]]]]]];
plot4=Show[Graphics[{Text["optimal solution: ",{0,0.75},{1,0}],Text[ToString[PaddedForm[f[reddot[[1]],reddot[[2]]],{3,2}]],{0.25,0.75},{1,0}],Text["current solution: ",{0,0.6},{1,0}],Text[ToString[PaddedForm[f[snap[[1]],snap[[2]]],{3,2}]],{0.25,0.6},{1,0}]}],PlotRange->{{-1,0.45},{-1,1}},PlotLabel->"objective values"];
Grid[{{plot1,plot2},{plot3,plot4}},Spacings->{3,3}]
],{{c,0,"boundary height"},{0->"low",co->"mid",2->"high"},RadioButton},{{pt,{0,-0.5}},Locator},SaveDefinitions->True]