Convex Optimization

Rob Knapp, Wolfram Research

Join the Conversation #WolframTechConf

ABSTRACT

This presentation will provide a general overview of convex optimization with a focus on modeling techniques for optimization problems that take advantage of new solvers to be included in the Wolfram Language.

Why the fuss?

I’ll try to explain this via a rough timeline (taken from Wikipedia and other sources)

1800’s Fourier publishes a method for solving systems of  linear inequalities

Late 1930’s.

Soviet economist Leonid Kantorovich, formulated the general linear programming problem and a method for solving it

American economist  T. C. Koopmans formulated classic economic problems as linear programming problems.

Kantorovich and Koopmans shared the 1975 Nobel prize in economics.

Late 1940’s   George Dantzig independently developed a general linear programming formulation

He developed the simplex algorithm as an efficient way to solve them.

When Dantzig presented the simplex algorithm to John von Neumann he conjectured the theory of duality.

Dantzig published a formal proof of the duality theory in 1948.

1979  Leonid Khachiyan showed linear optimization to be solvable in polynomial time

1984  Narendra Karmarkar introduced a new interior-point (polynomial time) method for solving linear-programming problems

Late 1980’s   development of interior point methods for nonlinear convex optimization

1994 Nesterov & Nemirovski show that convex optimization problems can be solved in polynomial time

Late 90’s - now

Development of convex problem classes and solvers like second order cone optimization, semidefinite optimization

Advances in modeling techniques developed by many

Packages CVX and Yalmip provide somewhat comprehensible interfaces to solver libraries

Courses and textbooks (e.g. Boyd and Vanderberghe) become more mainstream.

The big upshot is that convex optimization opens up a whole class of nonlinear problems that can now be considered solved with the help of efficient solvers and appropriate modeling transformations.

What’s new?

Special Convex Optimization Solvers

LinearOptimization — minimize {c.x|a.x+b \[VectorGreaterEqual]0}

LinearFractionalOptimization — minimize {(α.x+β)/(γ.x+δ|a.x+b\[VectorGreaterEqual]0}

QuadraticOptimization — minimize {(1/2)x.q.x+c.x|a.x+b\[VectorGreaterEqual]0}

SecondOrderConeOptimization — minimize {c.x|‖(Subscript[a, j]).x+(Subscript[b, j])‖≤(Subscript[α, j])+(Subscript[β, j]),j=1… k}

SemidefiniteOptimization — minimize {c.x|(Subscript[a, 0])+(Subscript[x, 1])(Subscript[a, 1])+…+(Subscript[x, k])(Subscript[a, k])(Subscript[\[VectorGreaterEqual], (Subscript[𝒮, +]^n)])0}

ConicOptimization — minimize {c.x|(Subscript[a, j]).x+(Subscript[b, j])(Subscript[\[VectorGreaterEqual], (Subscript[κ, j])])0, j=1… k}

Nina and Paritosh have presentations that will show how to effectively use  these solvers immediately following this presentation.

... What’s new?

General Optimization Solvers

FindMinimum — numerical local constrained optimization

NMinimize — numerical global constrained optimization

These are not new ....
... what is now is that when a problem is symbolically detected as convex and can be transformed to a form useable by the solvers, they are used to solve the problem efficiently.

Fitting

FitRegularization — regularization for Fit and FindFit

... What’s new?

Modeling tools

Vector inequalities

Vector (more generally array) valued variables with dimensional inference

Parameter “constraint” equations

Transformations to efficient canonical forms

Vector inequalities

VectorLess  {(Subscript[x, 1]),…,(Subscript[x, n])}\[VectorLess]{(Subscript[y, 1]),…,(Subscript[y, n])} if  (Subscript[x, i])<(Subscript[y, i]) for  i=1 … n

VectorLessEqual  {(Subscript[x, 1]),…,(Subscript[x, n])}\[VectorLessEqual]{(Subscript[y, 1]),…,(Subscript[y, n])} if  (Subscript[x, i])≤(Subscript[y, i]) for  i=1 … n

VectorGreaterEqual  {(Subscript[x, 1]),…,(Subscript[x, n])}\[VectorGreaterEqual]{(Subscript[y, 1]),…,(Subscript[y, n])} if  (Subscript[x, i])≥(Subscript[y, i]) for  i=1 … n

VectorGreater  {(Subscript[x, 1]),…,(Subscript[x, n])}\[VectorGreater]{(Subscript[y, 1]),…,(Subscript[y, n])} if  (Subscript[x, i])>(Subscript[y, i]) for  i=1 … n

These are only partial orderings

Manipulate[Quiet[RegionPlot[{x,y}\[VectorGreaterEqual]p,{x,-2,2},{y,-2,2}]],{{p,{0,0}},Locator}]

{1,2}\[VectorLess]{3,4}

If any element are not Less, then it is False

{1,2}\[VectorLess]{1,3}

Either side can be scalar (or of any dimensions compatible with the other side):

{1,2}\[VectorLess]a

{1,2}\[VectorLess]a /. a->3

{1,2}\[VectorLess]a /. a->{{1,2,3},{4,5}}

Generalized Vector inequalities

The generalized inequality x(Underscript[\[VectorGreaterEqual], κ])y gives True if x-y∈κ where κ is a proper convex cone. 
Some important examples of proper cones are

“NonNegativeCone” ((NonNegativeConeList[n])) x∈(ℝ^n) such that (Subscript[x, i])≥0

“NormCone” ((NormConeList[n])) x∈(ℝ^n) such that Norm[{(Subscript[x, 1]),…,(Subscript[x, n-1])}]≤(Subscript[x, n])

“SemidefiniteCone” ((SemidefiniteConeList[n]))  symmetric positive semidefinite matrices x∈(ℝ^(n×n))

Norm cone inequalities

“NormCone” ((NormConeList[n])) x∈(ℝ^n) such that Norm[{(Subscript[x, 1]),…,(Subscript[x, n-1])}]≤(Subscript[x, n])

Quiet[RegionPlot[{x,y,1}(Underscript[\[VectorGreaterEqual], (NormConeString)])0,{x,-1,1},{y,-1,1}]]

RegionPlot3D[{x,y,z}(Underscript[\[VectorGreaterEqual], (NormConeString)])0,{x,-1,1},{y,-1,1},{z,0,1}]

Semidefinite cone inequalities

“SemidefiniteCone” ((SemidefiniteConeList[n]))  symmetric positive semidefinite matrices x∈(ℝ^(n×n))

{{E,1},{1,π}}(Underscript[\[VectorGreaterEqual], (SemidefiniteConeString)])IdentityMatrix[2]

N[Eigenvalues[{{E,1},{1,π}}-IdentityMatrix[2]]]

Quiet[RegionPlot[{{x,1},{1,y}}(Underscript[\[VectorGreaterEqual], (SemidefiniteConeString)])0,{x,0,2},{y,0,2}]]

Manipulate[Quiet[RegionPlot[{{x,1},{1,y}}(Underscript[\[VectorLessEqual], (SemidefiniteConeString)])DiagonalMatrix[d],{x,-1,1},{y,-1,1}]],{{d,{0,0}},Slider2D}]

Parameter equation “constraints”

Here is a simple example

parameters = {a=={{1,2},{3,4},{5,6}},b=={7,8,9},λ==3};

NMinimize[{((Norm[a.x-b])^2)+   λ Norm[x, 1],parameters},x]

Now suppose we’d try to do this by assignment first:

SetAttributes[withAssignment, HoldAll];
withAssignment[e_] := Block[{a={{1,2},{3,4},{5,6}},b={7,8,9},λ=3},e];

withAssignment[NMinimize[((Norm[a.x-b])^2)+ λ Norm[x, 1],x]]

withAssignment[NMinimize[((Norm[Inactive[Plus][a.x,-b]])^2)+ λ Norm[x, 1],x]]

I’ll note that this can be done now with

withAssignment[Fit[{a,b},FitRegularization->{"L1",λ/2}]]

Normal[%]

Model Transformation

The goal here is to transform an optimization model (a set of objective function, constraints, and variables) into a form that can be used by the solvers that I showed you above.

Lets consider how to do this on the simple example above

Introduce a new vector variable t such that t==(Abs[x]).  This is done by adding constraints

tconstraints = {t\[VectorGreaterEqual]0,-t\[VectorLessEqual]x\[VectorLessEqual]t}; (* Note Vector inequalities *)

Since ((Norm[.])^2) is quadratic, we can use

QuadraticOptimization[((Norm[a.x-b])^2)+λ Total[t], {tconstraints,parameters}, {x,t}]

One aspect of model transformation that is important in this context is maintaining information about convexity

Convexity

A set C is affine if for x,y in C then θ x+(1-θ) y is also in C

A set C is convex if the above holds with θ≥0.

A function f is convex if f(θ x + (1-θ)y) ≤ θ f(x) + (1-θ)f(y) for θ≥0 and all x,y in the domain of f.

The epigraph of a function is the set {x,t} such that x is in the domain of f and f(x)≤ t

f[x_]:=x^2;
constrain[{x_,_}] := {x,f[x]};
plot = Plot[f[x],{x,-1,1},Filling->Top];
DynamicModule[{px={-1,1},py={1,1}},Show[plot,Graphics[{Dynamic[Line[{px,py}]],Locator[Dynamic[px,Function[px = constrain[#]]]],Locator[Dynamic[py,Function[py = constrain[#]]]]}]]]

Convex Optimization

A convex optimization problem is to

 ({{minimize, f(x)}, {subject to, ({{g(x)≤0 }, {h(x)==0}})}}) 
 
 where f and g={(Subscript[g, 1]),…}  are convex functions and  h is affine (i.e.h(x)== (Subscript[a, 0]).x+(Subscript[b, 0]))

Using the “epigraph transformation”, the problem may be expressed as

 ({{minimize, t}, {subject to, ({{f(x)≤t }, {g(x)\[VectorLessEqual]0}, {h(x)==0}})}})

So any convex optimization problem may be expressed as optimization of a linear function over a convex set.

This problem is convex:

FindMinimum[{x.{3,4},x\[VectorGreaterEqual]0,x∈Disk[]}, x]

This one is not:

FindMinimum[{x.{3,4},x\[VectorGreaterEqual]0,x∈Circle[]}, x]

plot = RegionPlot[Disk[{0,0},1,{0,π/2}]];
Manipulate[
xd = Flatten[FindArgMin[{x.c,x\[VectorGreaterEqual]0,x∈Disk[]}, x]]; 
xc=Flatten[FindArgMin[{x.c,x\[VectorGreaterEqual]0,x∈Circle[]}, x]];
Show[plot,Graphics[{PointSize[0.1],Point[xd],PointSize[.05],Blue,Point[xc],Text["c="<>ToString[c],{.8,.8}]}]],
{c,{-1,-1},{1,1}},ControlType->Slider2D,ControlPlacement->Left]

Conic Optimization

A convex optimization problem is to  ({{minimize, c.x}, {subject to, ({{h(x)==0}, {(Subscript[g, j])(x)≤0}})}}), j=1,…,k

Yet another way of expressing it is

({{minimize, c.x}, {subject to, ({{(Subscript[a, 0]).x+(Subscript[b, 0])==0}, {(Subscript[a, j]).x+(Subscript[b, j])(Underscript[\[VectorGreaterEqual], (Subscript[κ, j])])0}})}}) 

where (Subscript[κ, j]) is a convex proper cone,, j=1,…,k

This is called a conic linear optimization problem.

Note [http://www.seas.ucla.edu/~vandenbe/236C/lectures/conic.pdf] for a convex function (Subscript[f, i]) there is a proper convex cone (Subscript[κ, i]) such that (Subscript[f, i])(x) ≤ 0 is equivalent to (Subscript[a, i]).x+(Subscript[b, i])(Underscript[\[VectorGreaterEqual], (Subscript[κ, i])])0

For our solvers, the cones are restricted to a smaller set of cones:

Don’t worry; these are sufficient to handle a rich set of functions.

Model Transformation revisited

The general goal for the convex model transformation is to reduce a given optimization problem into a conic linear optimization problem.

An alternative to transforming the fitting problem (minimize ((Norm[a.x-b])^2)+λ Norm[x,1]) to quadratic optimization form is to get the linear conic form.  Since both Norm[a.x-b] and Norm[x,1] are convex, the problem can be expressed as

({{minimize, s + λ t}, {subject to, ({{(u^2)≤s}, {Norm[a.x+b]≤u}, {Norm[x,1]≤t}})}})

and we can solve it with

ConicOptimization[s+λ t,{u^2≤s,Norm[a.x-b]≤u,Norm[x,1]≤t,parameters},{x,s,t,u}]

Internally the ConicOptimization function has done more transformations to get it to linear conic form.  You can see this using

conicConstraints = ConicOptimization[s+3t,{u^2≤s,Norm[a.x+b]≤u,Norm[x,1]≤t,parameters},{x,s,t,u},"ConicConstraints"]

conicConstraints[[All,1,1,1,1]]=Map[MatrixForm,conicConstraints[[All,1,1,1,1]]];
conicConstraints

Another Example

Find the minimum volume ellipsoid that encloses a set of n points in dimension d:

SeedRandom[1234];
n = 100; d = 3;
points = RandomReal[1,{n,d}];

An ellipsoid in dimension d can be described by the set of points such that Norm[a.x+b]<=1 for a positive semidefinite symmetric d×d matrix a and vector b.

constraints = {VectorGreaterEqual[{a,0},{"SemidefiniteCone",d}],Map[(Norm[a.#+b]≤1)&, points]};

The thing to minimize is the volume of the ellipsoid.  For a symmetric positive semidefinite matrix, the volume of the ellipsoid is proportional to Log[Det[Inverse[a]]] = -Log[Det[a]].  To minimize the volume set the objective function to

objective = -Log[Det[a]];

and voilà:

{minval, rules} = FindMinimum[{objective, constraints},{a,Element[b,Vectors[d]]}]

The representation of the ellipsoid above is not the same as the graphics primitive, so we can do:

ainv = Inverse[a /. rules];
center = -ainv.(b /. rules);
Σ = ainv.ainv;
RegionMeasure[ellipsoid = Ellipsoid[center,Σ]]

approx = BoundingRegion[points,"FastEllipsoid"];
RegionMeasure[approx]

Show[RegionPlot3D[{ellipsoid,approx},PlotStyle->Opacity[.5],Mesh->None],Graphics3D[{PointSize[.025],Point[points]}]]

... Another Example Transformations

Constraint transformation use Schur complements. If the block matrix (({{A, B}, {C, D}})) is positive semidefinite, then so are A-B (D^(-1))C and D-C (A^(-1))B.

Thus (({{1, (Norm[a.(Subscript[p, i])+b])}, {(Norm[a.(Subscript[p, i])+b]), 1}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeString)])0 => 1-((Norm[a.(Subscript[p, i])+b])^2)≥0 => (Norm[a.(Subscript[p, i])+b])≤1

So all of the constraints can be expressed as semidefinite cone constraints.

The objective function uses a few more tricks.  Consider the problem of maximizing Log[Det[a]]
(equivalent to minimizing -Log[Det[a]])

Monotonicity: Log[x] is monotone increasing so if we maximize x we maximize Log[x]

Plot[Log[x],{x,.01,10}]

Determinant:  For a symmetric matrix, the determinant is the product of the eigenvalues, (Subscript[λ, 1])×…×(Subscript[λ, d])

Geometric mean,  The geometric mean of the eigenvalues is ((((Subscript[λ, 1])×…×(Subscript[λ, d])))^(1/d)). Since all the eigenvalues are positive, we can use the monotonicity of the power function, (x^(1/d)) for d > 0, and maximize the geometric mean.

Plot[(x^(1/d)),{x,0,10}]

Arithmetic mean.  The arithmetic mean is greater than or equal to the geometric mean: ((((Subscript[λ, 1])×…×(Subscript[λ, d])))^(1/d))≤(((Subscript[λ, 1])+…+(Subscript[λ, d]))/d)=((Tr[a])/d) and the product is increasing in each (Subscript[λ, i]) so if we maximize Tr[a] we also maximize the geometric mean

So to maximize Tr[a], we minimize -Tr[a] -- this is now just a linear objective

Are we done yet?

No ... 
Upcoming things we plan to do

Integration of some commercial libraries that can be used with the commercial license

Mixed integer-real convex optimization

Multiobjective convex optimization

Convex optimization with uncertainty

Summary

Convex optimization is a powerful tool that integrates very well in the Wolfram Language environment.

Extended class of problems that can be considered “solved”

Convex optimization solver functions with access solution details

Addition of modeling tools to the Wolfram language

Generalized vector inequalities

Dimensionality inference

Parameter specifications

Automated transformation engine to get to specific solvers