({{, (RefLink[SemidefiniteOptimization, paclet:ref/SemidefiniteOptimization])[f,cons,vars]
finds values of variables vars that minimize the linear objective f subject to semidefinite constraints cons.}, {, (RefLink[SemidefiniteOptimization, paclet:ref/SemidefiniteOptimization])[c,{(Subscript[a, 0]),(Subscript[a, 1]),…,(Subscript[a, k])}]
finds a vector x that minimizes the quantity c.x subject to the linear matrix inequality constraint (Subscript[a, k]) (Subscript[x, k])+(Subscript[a, 1]) (Subscript[x, 1])+(Subscript[a, 0])+…⪰0.}, {, (RefLink[SemidefiniteOptimization, paclet:ref/SemidefiniteOptimization])[…,"prop"]
specifies what solution property "prop" should be returned. }})

OpenersNotesSection(Spacer1)Details and Options

(RefLink[SemidefiniteOptimization, paclet:ref/SemidefiniteOptimization]) is also known as semidefinite programming (SDP).

Semidefinite optimization is a convex optimization problem that can be solved globally and efficiently.

Semidefinite optimization finds x that solves the primal problem:

The matrices (Subscript[a, i]) must be symmetric n ×n matrices.

The constraints cons can be specified by:

With (RefLink[SemidefiniteOptimization, paclet:ref/SemidefiniteOptimization])[f,cons,vars], parameter equations of the form par==val, where par is not in vars and val is numerical or an array with numerical values, may be included in the constraints to define parameters used in f or cons. »

The primal minimization problem has a related maximization problem that is the Lagrangian dual problem. The dual maximum value is always less than or equal to the primal minimum value, so it provides a lower bound. The dual maximizer provides information about the primal problem, including sensitivity of the minimum value to changes in the constraints. »

The semidefinite optimization has a dual: »

The possible solution properties "prop" include:

The following options may be given:

The option (RefLink[Method, paclet:ref/Method])->method may be used to specify the method to use. Available methods include:

Computations are limited to (RefLink[MachinePrecision, paclet:ref/MachinePrecision]).

OpenersPrimaryExamplesSection(Spacer1)Examples  (73)

OpenersExampleSection0(Spacer1)Basic Examples  (2)

Minimize  2(Subscript[x, 1])+3(Subscript[x, 2]) subject to the linear matrix inequality constraint  (Subscript[a, 0])+(Subscript[a, 1])(Subscript[x, 1])+(Subscript[a, 2]) (Subscript[x, 2])(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0:

(Subscript[a, 0])=(({{0, 1}, {1, 0}}));(Subscript[a, 1])=(({{1, 0}, {0, 0}}));(Subscript[a, 2])=(({{0, 0}, {0, 1}}));
res=SemidefiniteOptimization[2(Subscript[x, 1])+3(Subscript[x, 2]),(Subscript[a, 0])+(Subscript[a, 1])(Subscript[x, 1])+(Subscript[a, 2]) (Subscript[x, 2])(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,{(Subscript[x, 1]),(Subscript[x, 2])}]

The optimal point is where 2(Subscript[x, 1])+3(Subscript[x, 2]) is smallest within the region defined by the constraints:

Show[Plot3D[2(Subscript[x, 1])+3(Subscript[x, 2]),{(Subscript[x, 1]),(Subscript[x, 2])}∈ImplicitRegion[(Subscript[x, 1])(Subscript[x, 2])-1≥0&&(Subscript[x, 1])≥0&&(Subscript[x, 2])≥0,{(Subscript[x, 1]), (Subscript[x, 2])}], ],Graphics3D[{Blue, PointSize[0.04],Point[{(Subscript[x, 1]),(Subscript[x, 2]),2(Subscript[x, 1])+3(Subscript[x, 2])}/.res]}]]

Minimize  2(Subscript[x, 1])+3(Subscript[x, 2]) subject to the linear matrix inequality constraint  (Subscript[a, 0])+(Subscript[a, 1])(Subscript[x, 1])+(Subscript[a, 2]) (Subscript[x, 2])(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0:

(Subscript[a, 0])=(({{0, 1}, {1, 0}}));(Subscript[a, 1])=(({{1, 0}, {0, 0}}));(Subscript[a, 2])=(({{0, 0}, {0, 1}}));
res=SemidefiniteOptimization[2(Subscript[x, 1])+3(Subscript[x, 2]),(Subscript[a, 0])+(Subscript[a, 1])(Subscript[x, 1])+(Subscript[a, 2]) (Subscript[x, 2])(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,{(Subscript[x, 1]),(Subscript[x, 2])}]

Use the equivalent formulation with the objective vector and constraint matrices:

xmin = SemidefiniteOptimization[{2,3},{(Subscript[a, 0]),(Subscript[a, 1]),(Subscript[a, 2])}]

OpenersExampleSection1(Spacer1)Scope  (23)

OpenersExampleSubsection0(Spacer1)Basic Uses  (13)

Minimize x+2y subject to constraints (x^2)+2(y^2)≤3,x+y==2,x≥1:

res=SemidefiniteOptimization[x+2y,{(x^2)+2(y^2)≤ 3,x+y==2,x≥ 1},{x,y}]

Minimize x+y when the matrix (({{x, 1}, {1, y}})) is positive semidefinite:

constraint= VectorGreaterEqual[{{{x,1},{1,y}},0}, {"SemidefiniteCone",2}]

Find the solution:

SemidefiniteOptimization[x+y,constraint,{x,y}]

Minimize c.{x,y} subject to the linear matrix inequality constraint (Subscript[a, 1]) x+(Subscript[a, 2]) y+(Subscript[a, 0])(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0:

c = {1,1};(Subscript[a, 0])=(({{0, -1}, {-1, 0}}));(Subscript[a, 1])=(({{1, 1}, {1, 1}}));(Subscript[a, 2])=(({{1, 0}, {0, -1}}));

SemidefiniteOptimization[c.{x,y},(Subscript[a, 0])+(Subscript[a, 1])x+(Subscript[a, 2])y(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,{x,y}]

The left-hand side of the constraint can be given in evaluated form:

(Subscript[a, 0])+(Subscript[a, 1]) x+(Subscript[a, 2]) y//MatrixForm

SemidefiniteOptimization[x+y,(({{x+y, -1+x}, {-1+x, x-y}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,{x,y}]

Express the problem with the objective vector and constraint matrices:

SemidefiniteOptimization[c,{(Subscript[a, 0]),(Subscript[a, 1]),(Subscript[a, 2])}]

Use a vector variable v={x,y}:

c = {1,2};(Subscript[a, 0])=(({{0, 1}, {1, 0}}));(Subscript[a, 1])=(({{1, 0}, {0, 0}}));(Subscript[a, 2])=(({{0, 0}, {0, 1}}));

SemidefiniteOptimization[c.v,{(Subscript[a, 1]),(Subscript[a, 2])}.v(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])-(Subscript[a, 0]),v]

Use a vector variable v={x,y} and (RefLink[Inactive, paclet:ref/Inactive])[(RefLink[Plus, paclet:ref/Plus])] to avoid unintended threading:

c = {1,2};(Subscript[a, 0])=(({{0, 1}, {1, 0}}));(Subscript[a, 1])=(({{1, 0}, {0, 0}}));(Subscript[a, 2])=(({{0, 0}, {0, 1}}));

SemidefiniteOptimization[c.v,(Subscript[a, 0])+{(Subscript[a, 1]),(Subscript[a, 2])}.v(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,v]

Use a vector variable v={x,y} and parameter equations to avoid unintended threading:

parEqs = {(Subscript[a, 0])==(({{0, 1}, {1, 0}})),(Subscript[a, 1])==(({{1, 0}, {0, 0}})),(Subscript[a, 2])==(({{0, 0}, {0, 1}}))};

SemidefiniteOptimization[{1,2}.v,{parEqs,(Subscript[a, 0])+{(Subscript[a, 1]),(Subscript[a, 2])}.v(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0},v]

Use a vector variable x and (RefLink[Indexed, paclet:ref/Indexed])[x,i] to specify individual components:

SemidefiniteOptimization[{1,2}.x,(({{(IndexedDefault[x, 1]), 1}, {1, (IndexedDefault[x, 2])}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,x]

Use (RefLink[Vectors, paclet:ref/Vectors])[n] to specify the dimension of a vector variable when it is ambiguous:

SemidefiniteOptimization[(IndexedDefault[x, 1])+2(IndexedDefault[x, 2]),(({{(IndexedDefault[x, 1]), 1}, {1, (IndexedDefault[x, 2])}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,x]

SemidefiniteOptimization[(IndexedDefault[x, 1])+(IndexedDefault[2x, 2]),(({{(IndexedDefault[x, 1]), 1}, {1, (IndexedDefault[x, 2])}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,x∈Vectors[2]]

Several linear inequality constraints can be expressed with (RefLink[VectorGreaterEqual, paclet:ref/VectorGreaterEqual]):

SemidefiniteOptimization[x+y,VectorGreaterEqual[{{x+2y,x+1},{3,-1}}],{x,y}]

Use (Key0[esc]) v>= (Key0[esc]) or (RefLink[\[VectorGreaterEqual], paclet:ref/character/VectorGreaterEqual]) to enter the vector inequality sign \[VectorGreaterEqual]:

SemidefiniteOptimization[x+y,{x+2y,x}\[VectorGreaterEqual]{3,-1},{x,y}]

An equivalent form using scalar inequalities:

SemidefiniteOptimization[x+y,{x+2y≥3,x≥-1},{x,y}]

Use a vector variable v={x,y} and vector inequality:

SemidefiniteOptimization[{1,1}.v, {{1,2},{1,0}}.v\[VectorGreaterEqual]{3,-1},v]

Specify non-negative constraints using (RefLink[NonNegativeReals, paclet:ref/NonNegativeReals]) ((NonNegativeReals)):

SemidefiniteOptimization[{1,1}.v, {1,2}.v\[VectorGreaterEqual]3,v∈Vectors[2, NonNegativeReals]]

An equivalent form using vector inequalities:

SemidefiniteOptimization[{1,1}.v, {{{1,2},{1,0}}.v\[VectorGreaterEqual]{3,0},v\[VectorGreaterEqual]0},v]

Second-order cone constraints of the form (Norm[a.x+b])≤α.x+β can be used:

SemidefiniteOptimization[{1,1}.x,Norm[x]≤1,x]

"NormCone" constraints of the form a.x+b(Subscript[\[VectorGreaterEqual], (NormConeString)])0 can be used:

constraint = VectorGreaterEqual[{{x,y,3},0},{"NormCone",3}]

SemidefiniteOptimization[x+y,constraint,{x,y}]

OpenersExampleSubsection1(Spacer1)Primal Model Properties  (3)

Minimize the function x+y subject to the constraint (({{x, 1}, {1, y}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0:

SemidefiniteOptimization[x+y,(({{x, 1}, {1, y}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,{x,y}]

Get the primal minimizer as a vector:

{xmin,ymin}=SemidefiniteOptimization[x+y,(({{x, 1}, {1, y}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,{x,y},"PrimalMinimizer"]

Get the minimal value:

minval =SemidefiniteOptimization[x+y,(({{x, 1}, {1, y}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,{x,y},"PrimalMinimumValue"]

Show[Plot3D[x+y,{x,y}∈ImplicitRegion[x y≥1,{{x,0,5},{y,0,5}}],MeshFunctions->{#3&}], Graphics3D[{Red,PointSize[0.05],Point[{xmin,ymin,minval}]}]]

Extract the objective vector:

c=SemidefiniteOptimization[x+y,(({{x, 1}, {1, y}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,{x,y},"ObjectiveVector"]

Extract the constraint matrices:

mats=SemidefiniteOptimization[x+y,(({{x, 1}, {1, y}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,{x,y},"ConstraintMatrices"];

Map[MatrixForm, mats]

Use the extracted objective vector and constraint matrices for direct input:

SemidefiniteOptimization[c,mats]

The slack for an inequality g(x)(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[n])])0 at the minimizer (SuperStar[x]) is given by g((SuperStar[x])):

s=SemidefiniteOptimization[(Subscript[x, 1])+2(Subscript[x, 2]),(({{(Subscript[x, 1]), 1}, {1, (Subscript[x, 2])}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,{(Subscript[x, 1]),(Subscript[x, 2])},"Slack"]//Normal

Extract the minimizer and constraint matrices:

{(SuperStar[x]), a}=SemidefiniteOptimization[(Subscript[x, 1])+2(Subscript[x, 2]),
(({{(Subscript[x, 1]), 1}, {1, (Subscript[x, 2])}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,{(Subscript[x, 1]),(Subscript[x, 2])},{"PrimalMinimizer", "ConstraintMatrices"}]//Normal

Verify that the slack matrix satisfies (Subscript[a, 0])+(Subscript[a, 1]) (SuperStar[(Subscript[x, 1])])+(Subscript[a, 2])(SuperStar[(Subscript[x, 2])])-s==0:

Prepend[(SuperStar[x]),1].a-s

OpenersExampleSubsection2(Spacer1)Dual Model Properties  (3)

Minimize c.x subject to (Subscript[a, 0])+(Subscript[a, 1]) (Subscript[x, 1])+…+(Subscript[a, k]) (Subscript[x, k])⪰0:

c = {2,3};(Subscript[a, 0])=(({{0, 1}, {1, 0}}));(Subscript[a, 1])=(({{1, 0}, {0, 0}}));(Subscript[a, 2])=(({{0, 0}, {0, 1}}));
xv=Table[Indexed[x,i],{i,2}];

psol=SemidefiniteOptimization[c.xv,(Subscript[a, 0])+{(Subscript[a, 1]),(Subscript[a, 2])}.xv(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,xv]

The dual problem is to maximize -Tr[(Subscript[a, 0]).y] subject to Tr[(Subscript[a, i]).y]==(Subscript[c, i]),i=1,..., k,y\[VectorGreaterEqual]0:

constraints = {Table[Tr[(Subscript[a, i]).y]==c[[i]],{i,2}],y(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0};

dsol=SemidefiniteOptimization[Tr[(Subscript[a, 0]).y],constraints,y∈Matrices[2]]

The primal minimum value and the dual maximum value coincide because of strong duality:

{c.xv /. psol, -Tr[(Subscript[a, 0]).y]/.dsol}

That is the same as having a duality gap of zero. In general, c.(SuperStar[x])≤Tr[(Subscript[a, 0]).y] at optimal points:

SemidefiniteOptimization[c.xv,(Subscript[a, 0])+{(Subscript[a, 1]),(Subscript[a, 2])}.xv(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,xv,"DualityGap"]

Construct the dual problem using constraint matrices extracted from the primal problem:

obj=2x+3y;
constraints=(({{x, 1}, {1, y}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0;

SemidefiniteOptimization[obj,constraints,{x,y}]

Extract the objective vector and constraint matrices:

{c,a}=SemidefiniteOptimization[obj,constraints,{x,y},{"ObjectiveVector","ConstraintMatrices"}];

The dual problem is to maximize -Tr[(Subscript[a, 0]).y] subject to Tr[(Subscript[a, i]).y]==(Subscript[c, i]),y\[VectorGreaterEqual]0:

SemidefiniteOptimization[Tr[a[[1]].y],{Table[Tr[a[[i+1]].y]==c[[i]],{i,2}],y(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0},y∈Matrices[2]]

Get the dual maximum value and dual maximizer directly using solution properties:

obj=2x+3y;
constraints=(({{x, 1}, {1, y}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0;

The "DualMaximumValue" is:

SemidefiniteOptimization[obj,constraints,{x,y},"DualMaximumValue"]

The "DualMaximizer" can be obtained with:

SemidefiniteOptimization[obj,constraints,{x,y},"DualMaximizer"]//Normal

OpenersExampleSubsection3(Spacer1)Sensitivity Properties  (4)

Use "ConstraintSensitivity" to find the change in optimal value due to constraint perturbations:

c={2,3};(Subscript[a, 0])=(({{0, 3}, {3, 0}}));(Subscript[a, 1])=(({{1, 0}, {0, 0}}));(Subscript[a, 2])=(({{0, 0}, {0, 1}}));

(SuperStar[p])=SemidefiniteOptimization[c,{(Subscript[a, 0]),(Subscript[a, 1]),(Subscript[a, 2])},"PrimalMinimumValue"]

The sensitivity is a matrix:

sens=SemidefiniteOptimization[c,{(Subscript[a, 0]),(Subscript[a, 1]),(Subscript[a, 2])},"ConstraintSensitivity"]//Normal

Consider new constraints (Subscript[a, 0])+(Subscript[a, 1]) (Subscript[x, 1])+…+(Subscript[a, k]) (Subscript[x, k])+δ⪰0 where δ is the perturbation:

δ=0.01*{{1,2},{2,3}};

The approximate new optimal value is:

(SuperStar[p])+Tr[δ.sens]

Compare to directly solving the perturbed problem:

SemidefiniteOptimization[c,{(Subscript[a, 0])+δ,(Subscript[a, 1]),(Subscript[a, 2])},"PrimalMinimumValue"]

The optimal value changes according to the signs of the sensitivity matrix elements:

c={2,3};(Subscript[a, 0])=(({{0, 3}, {3, 0}}));(Subscript[a, 1])=(({{1, 0}, {0, 0}}));(Subscript[a, 2])=(({{0, 0}, {0, 1}}));

{(SuperStar[p]),sens}=SemidefiniteOptimization[c,{(Subscript[a, 0]),(Subscript[a, 1]),(Subscript[a, 2])},{"PrimalMinimumValue","ConstraintSensitivity"}];

sens//Normal

At negative sensitivity element position, a positive perturbation will decrease the optimal value:

(SuperStar[p])>SemidefiniteOptimization[c,{(Subscript[a, 0])+{{0.01,0},{0,0}},(Subscript[a, 1]),(Subscript[a, 2])},"PrimalMinimumValue"]

At positive sensitivity element position, a positive perturbation will increase the optimal value:

(SuperStar[p])<SemidefiniteOptimization[c,{(Subscript[a, 0])+{{0,0.01},{0.01,0}},(Subscript[a, 1]),(Subscript[a, 2])},"PrimalMinimumValue"]

Express the perturbed constraints symbolically using Sylvester's criterion for semidefiniteness:

c={2,3};(Subscript[a, 0])=(({{0, 3}, {3, 0}}));(Subscript[a, 1])=(({{1, 0}, {0, 0}}));(Subscript[a, 2])=(({{0, 0}, {0, 1}}));

parameters=(({{(Subscript[p, 1,1]), (Subscript[p, 1,2])}, {(Subscript[p, 1,2]), (Subscript[p, 2,2])}}));

matrix = (Subscript[a, 0])+parameters+(Subscript[x, 1]) (Subscript[a, 1])+(Subscript[x, 2]) (Subscript[a, 2]);
constraints = {matrix[[1,1]]≥0,matrix[[2,2]]≥ 0, Det[matrix]≥0}

With this form, the minimum value can be found exactly as a function of the parameters around 0:

Assuming[-1<(Subscript[p, 1,2])<1,exact=Refine[MinValue[2 (Subscript[x, 1]) + 3 (Subscript[x, 2]),
 {(Subscript[x, 1])+(Subscript[p, 1,1]) ≥0,(Subscript[x, 2])+(Subscript[p, 2,2]) ≥ 0, -1+(Subscript[x, 1]) (Subscript[x, 2])+(Subscript[x, 2]) (Subscript[p, 1,1])-2 (Subscript[p, 1,2])-(Subscript[p, 1,2]^2)+(Subscript[x, 1]) (Subscript[p, 2,2])+(Subscript[p, 1,1]) (Subscript[p, 2,2])≥ 0},{(Subscript[x, 1]),(Subscript[x, 2])}]]]

Make a symmetric matrix with the derivatives of the minimum with respect to the parameters:

MatrixForm[N[smat=(({{D[exact,(Subscript[p, 1,1])], (1/2)D[exact,(Subscript[p, 1,2])]}, {(1/2)D[exact,(Subscript[p, 1,2])], D[exact,(Subscript[p, 2,2])]}}))]]

This is the sensitivity to perturbation given by the "ConstraintSensitivity" property:

SemidefiniteOptimization[{2,3},{(Subscript[a, 0]),(Subscript[a, 1]),(Subscript[a, 2])},"ConstraintSensitivity"]//MatrixForm

The constraint sensitivity can also be obtained as the negative of the dual maximizer:

c={2,3};(Subscript[a, 0])=(({{0, 3}, {3, 0}}));(Subscript[a, 1])=(({{1, 0}, {0, 0}}));(Subscript[a, 2])=(({{0, 0}, {0, 1}}));

-SemidefiniteOptimization[c,{(Subscript[a, 0]),(Subscript[a, 1]),(Subscript[a, 2])},"DualMaximizer"]//MatrixForm

SemidefiniteOptimization[c,{(Subscript[a, 0]),(Subscript[a, 1]),(Subscript[a, 2])},"ConstraintSensitivity"]//MatrixForm

OpenersExampleSection2(Spacer1)Options  (8)

OpenersExampleSubsection4(Spacer1)Method  (5)

The default method "CSDP" is an interior point method:

SemidefiniteOptimization[(Subscript[x, 1])+(Subscript[x, 2]),(({{(Subscript[x, 1]), 1}, {1, (Subscript[x, 2])}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[n])])0,{(Subscript[x, 1]),(Subscript[x, 2])}, Method->Automatic]

SemidefiniteOptimization[(Subscript[x, 1])+(Subscript[x, 2]),(({{(Subscript[x, 1]), 1}, {1, (Subscript[x, 2])}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[n])])0,{(Subscript[x, 1]),(Subscript[x, 2])}, Method->"CSDP"]

"DSDP" is an alternative interior point method:

SemidefiniteOptimization[(Subscript[x, 1])+(Subscript[x, 2]),(({{(Subscript[x, 1]), 1}, {1, (Subscript[x, 2])}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[n])])0,{(Subscript[x, 1]),(Subscript[x, 2])}, Method->"DSDP"]

"SCS" uses a splitting conic solver method:

SemidefiniteOptimization[(Subscript[x, 1])+(Subscript[x, 2]),(({{(Subscript[x, 1]), 1}, {1, (Subscript[x, 2])}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[n])])0,{(Subscript[x, 1]),(Subscript[x, 2])}, Method->"SCS"]

Different methods have different default tolerances, which affects the accuracy and precision:

{obj,cons,var}={(Subscript[x, 1])+(Subscript[x, 2]),(({{(Subscript[x, 1]), 1}, {1, (Subscript[x, 2])}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,{(Subscript[x, 1]),(Subscript[x, 2])}};

Compute exact and approximate solutions:

{minval,argmin}=Minimize[obj,cons,var];
exact={minval,var/.argmin}

solCSDP=SemidefiniteOptimization[obj,cons,var, {"PrimalMinimumValue", "PrimalMinimizer"}, Method->"CSDP"];

solDSDP=SemidefiniteOptimization[obj,cons,var, {"PrimalMinimumValue", "PrimalMinimizer"}, Method->"DSDP"];

solSCS=SemidefiniteOptimization[obj,cons,var, {"PrimalMinimumValue", "PrimalMinimizer"}, Method->"SCS"];

"CSDP" and "DSDP" have default tolerances of (10^(-8)):

solCSDP-exact

solDSDP-exact

"SCS" has a default tolerance of (10^(-3)):

solSCS-exact

When the default method "CSDP" produces a message, try "DSDP" first:

a = {{{0,2},{2,0}},{{1,2},{2,3}},{{5,0},{0,1}}};
c = {2,3};

{xmin, min}=SemidefiniteOptimization[c,a, {"PrimalMinimizer", "PrimalMinimumValue"},Method->"CSDP"]

In this case, "DSDP" succeeds in finding a good solution:

{xmin, min}=SemidefiniteOptimization[c,a, {"PrimalMinimizer", "PrimalMinimumValue"},Method->"DSDP"]

"SCS" with a default tolerance of (10^(-3)) is an alternative method to try:

{xmin, min}=SemidefiniteOptimization[c,a, {"PrimalMinimizer", "PrimalMinimumValue"},Method->"SCS"]

The quality of the result with "SCS" can often be improved with a smaller (RefLink[Tolerance, paclet:ref/Tolerance]):

{xmin, min}=SemidefiniteOptimization[c,a, {"PrimalMinimizer", "PrimalMinimumValue"},Method->{"SCS", Tolerance->10^-8}]

OpenersExampleSubsection5(Spacer1)Tolerance  (2)

A smaller (RefLink[Tolerance, paclet:ref/Tolerance]) setting gives a more precise result:

c={1,1};(Subscript[a, 0])=(({{0, 1}, {1, 0}}));(Subscript[a, 1])=(({{1, 0}, {0, 0}}));(Subscript[a, 2])=(({{0, 0}, {0, 1}}));

Compute the exact minimum value with (RefLink[Minimize, paclet:ref/Minimize]):

exact =Minimize[{(Subscript[x, 1])+(Subscript[x, 2]),(Subscript[a, 0])+(Subscript[a, 1])(Subscript[x, 1])+(Subscript[a, 2])(Subscript[x, 2])(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0}, {(Subscript[x, 1]),(Subscript[x, 2])}][[1]]

Compute the error in the minimum value with different (RefLink[Tolerance, paclet:ref/Tolerance]) settings:

err=Table[{tol,Abs[SemidefiniteOptimization[c,{(Subscript[a, 0]),(Subscript[a, 1]),(Subscript[a, 2])},"PrimalMinimumValue",Tolerance->tol]-exact]},{tol,10.^-Range[7]}]

Visualize the change in minimum value error with respect to tolerance:

ListLogLogPlot[err, PlotRange->All,Joined->True,PlotMarkers->Automatic]

A smaller (RefLink[Tolerance, paclet:ref/Tolerance]) setting gives a more precise result, but may take longer to compute:

n=2000;c = ConstantArray[1,n];
(Subscript[a, 0])=SparseArray[{{1,n}->1, {n,1}->1},{n,n}];Table[(Subscript[a, i])=SparseArray[{{i,i}->1},{n,n}], {i,1,n}];

A smaller tolerance takes longer:

Timing[ptight=SemidefiniteOptimization[c,Table[(Subscript[a, i]), {i,0,n}],"PrimalMinimumValue",Tolerance->(10^(-8))];]

Timing[ploose=SemidefiniteOptimization[c,Table[(Subscript[a, i]), {i,0,n}],"PrimalMinimumValue",Tolerance->(10^(-1))];]

The tighter tolerance gives a more precise answer:

{ptight,ploose}

OpenersExampleSubsection6(Spacer1)PerformanceGoal  (1)

The default value of the option (RefLink[PerformanceGoal, paclet:ref/PerformanceGoal]) is (RefLink[$PerformanceGoal, paclet:ref/$PerformanceGoal]):

$PerformanceGoal

Use (RefLink[PerformanceGoal, paclet:ref/PerformanceGoal])->"Quality" to get a more accurate result:

c={1,1};(Subscript[a, 0])=(({{0, 1}, {1, 0}}));(Subscript[a, 1])=(({{1, 0}, {0, 0}}));(Subscript[a, 2])=(({{0, 0}, {0, 1}}));
SemidefiniteOptimization[c,{(Subscript[a, 0]),(Subscript[a, 1]),(Subscript[a, 2])}, PerformanceGoal->"Quality"]

SemidefiniteOptimization[c,{(Subscript[a, 0]),(Subscript[a, 1]),(Subscript[a, 2])}, PerformanceGoal->"Speed"]

Use (RefLink[PerformanceGoal, paclet:ref/PerformanceGoal])->"Speed" to get a result faster, but at the cost of quality:

n=2000;c = ConstantArray[1,n];
(Subscript[a, 0])=SparseArray[{{1,n}->1, {n,1}->1},{n,n}];Table[(Subscript[a, i])=SparseArray[{{i,i}->1},{n,n}], {i,1,n}];

Compare the timings:

Timing[resSpeed=SemidefiniteOptimization[c,Table[(Subscript[a, i]), {i,0,n}],"PrimalMinimumValue",PerformanceGoal->"Speed"];]

Timing[resQuality=SemidefiniteOptimization[c,Table[(Subscript[a, i]), {i,0,n}],"PrimalMinimumValue",PerformanceGoal->"Quality"];]

The "Speed" goal gives a less accurate result:

{resSpeed, resQuality}

OpenersExampleSection3(Spacer1)Applications  (29)

OpenersExampleSubsection7(Spacer1)Basic Modeling Transformations  (13)

Maximize x+y subject to (({{-x, 1}, {1, -y}}))⪰0. Solve a maximization problem by negating the objective function:

SemidefiniteOptimization[-(x+y),(({{-x, 1}, {1, -y}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,{x,y}]

Negate the primal minimum value to get the corresponding maximal value:

-SemidefiniteOptimization[-(x+y),(({{-x, 1}, {1, -y}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,{x,y},"PrimalMinimumValue"]

Find {x,y} that minimizes the largest eigenvalue of a symmetric matrix that depends linearly on the decision variables {x,y}, A(x,y) =(Subscript[a, 0])+(Subscript[a, 1])x+(Subscript[a, 2])y. The problem can be formulated as linear matrix inequality, since λℐ⪰A(x,y) is equivalent to λ≥(Subscript[λ, 1]),…,λ≥(Subscript[λ, n]) where (Subscript[λ, i]) is the i(^th) eigenvalue of A(x,y). Define the linear matrix function A(x,y):

A[x_,y_]:=(({{2, -3}, {-3, 4}}))+(({{1, 2}, {2, 3}})) x+(({{4, 5}, {5, 6}})) y

A[x,y]//MatrixForm

A real symmetric matrix A can be diagonalized with an orthogonal matrix Q so (Q^T).A.Q=diag({(Subscript[λ, 1]),…,(Subscript[λ, n])}). Hence (Q^T).(λ ℐ-A).Q=diag({λ-(Subscript[λ, 1]),…,λ-(Subscript[λ, n])})⪰0 iff λ≥(Subscript[λ, 1]),…,λ≥(Subscript[λ, n]). Since any u=Q .(Q^T)u, taking v=(Q^T).u, (u^T).(λ ℐ-A). u=(v^T).(Q^T).(λ ℐ-A). Q .v≥0, hence λ ℐ-A⪰0 iff λ≥(Subscript[λ, 1]),…,λ≥(Subscript[λ, n]). Numerically simulate to show that these formulations are equivalent:

eigs=Eigenvalues[A[1.,2.]]

data=Table[{λ,λ (({{1, 0}, {0, 1}})) (Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])A[1.,2.]},{λ,RandomReal[{0,35},100]}];
{true,false}={Cases[data,{λ_,True}:>λ],Cases[data,{λ_,False}:>λ]};

NumberLinePlot[{eigs,true,false},]

The resulting problem:

res=SemidefiniteOptimization[λ,λ (({{1, 0}, {0, 1}})) (Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])A[x,y],{λ,x,y}]

Run a Monte Carlo simulation to check the plausibility of the result:

Min[Table[Max[Eigenvalues[A@@RandomReal[{-10000,10000},2]]],10000]]

Find {x,y} that maximizes the smallest eigenvalue of a symmetric matrix A(x,y) that depends linearly on the decision variables {x,y}. Define the linear matrix function A(x,y):

A[x_,y_]:=(({{x+4 y+2, 2 x+5 y-3}, {2 x+5 y-3, 3 x+6 y+4}}));

The problem can be formulated as linear matrix inequality, since λℐ⪯A(x,y) is equivalent to λ⪯(Subscript[λ, 1]),…,λ⪯(Subscript[λ, n]) where (Subscript[λ, i]) is the i(^th) eigenvalue of A(x,y). To maximize λ, minimize -λ:

res=SemidefiniteOptimization[-λ,
λ (({{1, 0}, {0, 1}})) (Underscript[\[VectorLessEqual], (SemidefiniteConeList[2])])A[x,y],{λ,x,y}]

Run a Monte Carlo simulation to check the plausibility of the result:

Max[Table[Min[Eigenvalues[A@@RandomReal[{-10000,10000},2]]],10000]]

Find {x,y} that minimizes the difference between the largest and the smallest eigenvalues of a symmetric matrix A(x,y) that depends linearly on the decision variables {x,y}. Define the linear matrix function A(x,y):

A[x_,y_]:=(({{x+4 y+2, 2 x+5 y-3}, {2 x+5 y-3, 3 x+6 y+4}}));

The problem can be formulated as linear matrix inequality, since λℐ⪯A(x,y)⪯μℐ is equivalent to λ⪯(Subscript[λ, 1])⪯μ,…,λ⪯(Subscript[λ, n])⪯μ where (Subscript[λ, i]) is the i(^th) eigenvalue of A(x,y). Solve the resulting problem:

res=SemidefiniteOptimization[μ-λ,{λ (({{1, 0}, {0, 1}})) (Underscript[\[VectorLessEqual], (SemidefiniteConeList[2])])A[x,y](Underscript[\[VectorLessEqual], (SemidefiniteConeList[2])])μ (({{1, 0}, {0, 1}}))},{λ,μ,x,y}]

In this case, the minimum and maximum eigenvalues coincide and the difference is 0:

{Eigenvalues[A[x,y]/.res],Chop[μ-λ/.res, 10^-7]}

Minimize the largest by absolute value eigenvalue of a linear in {x,y} symmetric matrix A(x,y):

A[x_,y_]:=(({{x+4 y+2, 2 x+5 y-3}, {2 x+5 y-3, 3 x+6 y+4}}));

The largest eigenvalue satisfies λ I-A(x,y)⪰0. The largest by absolute value negative eigenvalue of A(x,y) is the largest eigenvalue of -A(x,y) and satisfies λ I+A(x,y)⪰0:

res=SemidefiniteOptimization[λ,{-λ (({{1, 0}, {0, 1}}))(Underscript[\[VectorLessEqual], (SemidefiniteConeList[2])])A[x,y](Underscript[\[VectorLessEqual], (SemidefiniteConeList[2])])λ (({{1, 0}, {0, 1}}))},{λ,x,y}]

Eigenvalues[A[x,y]/.res]

Find {x,y} that minimizes the largest singular value σ of a linear in {x,y} matrix A(x,y):

A[x_,y_]:=(({{1+x+2 y, 2+3 y, x+1}, {3+5 y, 4-x+8 y, 2y}}));

The largest singular value σ of A(x,y)  is the square root of the largest eigenvalue of ((A )^T)A and from a preceding example it satisfies (σ^2) I-((A(x,y))^T)A(x,y)⪰0  or equivalently (({{σ I, A(x,y)}, {((A(x,y))^T), σ I}}))⪰0:

res=SemidefiniteOptimization[σ,{{σ (({{1, 0}, {0, 1}})),A[x,y]},{A[x,y]\[Transpose],σ (({{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}))}} (Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[5])])0,{x,y,σ}]

Plot the result:

f[x,y] := Module[{v = SingularValueList[A[x,y]]},Max[v]];
Show[Plot3D[f[x,y],{x,-10,10},{y,-2,2}],Graphics3D[{Red,PointSize[.05],Point[res[[All,2]]]}]]

Minimize ||a.x-b||. Using an auxiliary variable t, transform the problem to minimize t such that ||a.x-b||≤t. This is the same as (a.x-b)ᵀ.(a.x-b)≤(t^2)⟺t-(a.x-b)ᵀ.(((t I))^(-1))(a.x-b)≥0:

a={{1,2},{3,4},{5,6}};
b={7,8,10};

A Schur complement condition says that if C\[VectorGreater]0, a block matrix (({{A, B}, {(B^T), C}}))\[VectorGreaterEqual]0 iff A-B.(C^(-1)).(B^T)\[VectorGreaterEqual]0. Thus ||a.x-b||≤t ⧦ (({{t, (((a.x-b))^T)}, {a.x-b, t I}}))\[VectorGreaterEqual]0 for t>0 and for t=0, since then a.x-b must be 0:

SemidefiniteOptimization[t,{ℐ==IdentityMatrix[3],a1==a,b1==b,(({{t, a1.x-b1}, {a1.x-b1, t ℐ}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeString)])0},{t,x}]

Use the constraint ||a.x-b||≤t directly and it will automatically convert into semidefinite form:

SemidefiniteOptimization[t,{ℐ==IdentityMatrix[3],a1==a,b1==b,Norm[a1.x-b1]≤t},{t,x}]

Minimize (((α.x+β))^2)/(γ.x+δ) subject to a.x+b≥0, assuming γ.x+δ>0 when a.x+b≥0. Using the auxiliary variable t, the objective is to minimize t such that (((α.x+β))^2)/(γ.x+δ)≤t:

{α,β,γ,δ}={{1,2},1,{-3,4},2};
{a,b}={{{1,-2},{-3,4}},{-5,-1}};

Check that a.x+b≥0 implies γ.x+δ>0:

Reduce[Implies[a.{(IndexedDefault[x, 1]), (IndexedDefault[x, 2])}+b>=0, γ.{(IndexedDefault[x, 1]), (IndexedDefault[x, 2])}+δ>0],Reals]

Using the Schur complement condition as in the preceding, (((α.x+β))^2)/(γ.x+δ)≤t iff (({{t, α.x+β}, {α.x+β, γ.x+δ }}))⪰0. Use (RefLink[Inactive, paclet:ref/Inactive])[(RefLink[Plus, paclet:ref/Plus])] for constructing the constraints to avoid threading:

SemidefiniteOptimization[t,{(({{t, α.x+β}, {α.x+β, γ.x+δ}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0, a.x+b\[VectorGreaterEqual]0},{t,x∈Vectors[2]}]

For quadratic sets  (Subscript[𝒬, i])={x∈(ℝ^n)|(x^T).(Subscript[a, i]) .x+2 x.(Subscript[b, i]^T)+(Subscript[c, i])≤0}, which include ellipsoids, quadratic cones and paraboloids, determine whether (Subscript[𝒬, 1])⊆(Subscript[𝒬, 2]), where (Subscript[a, i]) are symmetric matrices, (Subscript[b, i]) are vectors and (Subscript[c, i]) scalars:

{a1,b1,c1}={{{0.7,0.68},{0.68,0.94}},{0.35,.42},0.12};{a2,b2,c2}={{{0.13,0.2},{0.2,0.62}},{.2,.4},.01};

With[{x={x1,x2}},RegionPlot[{x.a1.x+2 b1.x+c1≤0,x.a2.x+2 b2.x+c2≤0},{x1,-2,2},{x2,-2,2},PlotLegends->{(Subscript[𝒬, 1]),(Subscript[𝒬, 2])},ImageSize->100]]

Assuming that the sets (Subscript[𝒬, i]) are full dimensional, the S-procedure says that  (Subscript[𝒬, 1])⊆(Subscript[𝒬, 2]) iff there exists some non-negative number λ such that λ (({{(Subscript[a, 1]), (Subscript[b, 1])}, {(Subscript[b, 1]^T), (Subscript[c, 1])}})) ⪰ (({{(Subscript[a, 2]), (Subscript[b, 2])}, {(Subscript[b, 2]^T), (Subscript[c, 2])}})). Visually see that there exists a non-negative λ:

m1=ArrayFlatten[{{a1,{b1}\[Transpose]},{{b1},{{c1}}}}];
m2=ArrayFlatten[{{a2,{b2}\[Transpose]},{{b2},{{c2}}}}];
Plot[Boole[λm1(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[3])])m2],{λ,0,5},ImageSize->Small]

Since λ≥0, (Subscript[𝒬, 1])⊆(Subscript[𝒬, 2]):

SemidefiniteOptimization[0, λ(({{a1, b1}, {b1, c1}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[3])])(({{a2, b2}, {b2, c2}})),λ]

Minimize |x+y| subject to (({{x, 1}, {1, y}}))⪰0. Convert the objective|x+y| into a linear function t with the additional constraint |x+y|≤t, which is equivalent to  -t≤x+y≤t:

SemidefiniteOptimization[t,{(({{x, 1}, {1, y}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0, -t≤x+y≤t},{x,y,t}]

Minimize ||{x+y,x-y}(Subscript[||, 1]) subject to (({{x, 1}, {1, y}}))⪰0. Convert the objective into a linear function using ||{x+y,x-y}(Subscript[||, 1])≤(Subscript[s, 1])+(Subscript[s, 2]) and the additional constraints -(Subscript[s, 1])≤x+y≤(Subscript[s, 1]),-(Subscript[s, 2])≤x-y≤(Subscript[s, 2]):

SemidefiniteOptimization[s1+s2,{(({{x, 1}, {1, y}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0, -s1≤x+y≤s1,-s2≤x-y≤s2},{x,y,s1,s2}]

Minimize ||{x+y,x-y}(Subscript[||, ∞]) subject to (({{x, 1}, {1, y}}))⪰0. Convert the objective into a linear function s with the additional constraint ||{x+y,x-y}(Subscript[||, ∞])≤ s, which is equivalent to -s≤x+y≤s,-s≤x-y≤s:

SemidefiniteOptimization[s,{(({{x, 1}, {1, y}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0, -s≤x+y≤s,-s≤x-y≤s},{x,y,s}]

Minimize f(c.x) subject to (Subscript[a, 0])+(∑_(i=1))(Subscript[a, i])(Subscript[x, i])≥0, where f is a nondecreasing function, by instead minimizing c.x. The primal minimizer (SuperStar[x]) will remain the same for both problems. Consider minimizing f(z)=(Sqrt[z]), z=x+y, z≥0, subject to  (({{x, 1}, {1, y}}))⪰0:

f=Function[{s},Sqrt[s]];

{minval, (Subscript[z, min])}=SemidefiniteOptimization[x+y,{(({{x, 1}, {1, y}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,x+y≥0},{x,y},{"PrimalMinimumValue","PrimalMinimizerRules"}]

The true minimum value can be obtained by applying f to the minimum value of  z=x+y:

f[minval]

OpenersExampleSubsection8(Spacer1)Data-Fitting Problems  (3)

Find the coefficients c of a fifth-order polynomial by minimizing (Norm[a.c-b]) that fits a discrete dataset:

data=;
dataplot=ListPlot[data, PlotRange->All]

Select the polynomial bases and construct the input matrix using (RefLink[DesignMatrix, paclet:ref/DesignMatrix]) and output vector:

bases=Map[(x^(#))&,Range[0,5]];
input=DesignMatrix[data,bases,x];
output=data[[All,2]];

Using an auxiliary variable t, the objective is transformed to minimize t such that ||a.c-b||≤t, which is equivalent to (({{t, (((a.c-b))^T)}, {a.c-b, t I}}))⪰0  as shown under Basic Modeling Transformations:

res=SemidefiniteOptimization[t,{ℐ==IdentityMatrix[Length[input]],a==input, b==output,(({{t, a.c-b}, {a.c-b, t ℐ }}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeString)])0},{t,c}]

Compare fit with data:

Show[dataplot, Plot[Evaluate[bases.(c/.res)],{x,0,1}]]

Find an approximating function to discrete data that varies on a logarithmic scale by minimizing (Norm2[a.x-b, ∞])=(max_i)|(Subscript[a, i]).x-b| using Chebyshev bases:

data = ;

Select Chebyshev basis functions and compute their values at the random data points:

basis=Table[ChebyshevT[k,s],{k,0,5}];
a=DesignMatrix[data,basis,s];
b = data[[All,2]];

Since the data is on a logarithmic scale, direct data-fitting is not ideal. Instead, transform the problem to minimize (max_i)|log((Subscript[a, i]^T).x)-log ((Subscript[b, i]))|=log ((max_(i ))((Subscript[a, i])ᵀ.x/(Subscript[b, i]), (Subscript[b, i])/(Subscript[a, i])ᵀ.x)). Using auxiliary variable t, minimize t such that 1/t≤(Subscript[a, i]).x/(Subscript[b, i])≤t. This constraint is equivalent to (({{t-(Subscript[a, i]).x/(Subscript[b, i]), 0, 0}, {0, (Subscript[a, i]).x/(Subscript[b, i]), 1}, {0, 1, t}}))⪰0:

constraints=MapThread[(({{t-((#1.x)/(#2)), 0, 0}, {0, ((#1.x)/(#2)), 1}, {0, 1, t}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeString)])0&,{a,b}];

Find the coefficients of the approximating function:

res=SemidefiniteOptimization[t,constraints,{t,x}]

The resulting fit is:

chebyshevFit=(basis.x)/. res

Visualize the fit:

Show[ListLogPlot[data],LogPlot[chebyshevFit,{s,-1,1}]]

The data-fitting can also be obtained directly using the function (RefLink[Fit, paclet:ref/Fit]). However, without the log transformation, there are significant oscillations in the approximating function:

fit1=basis.Fit[{a,b},NormFunction->Infinity];

Show[ListLogPlot[data],LogPlot[{chebyshevFit,fit1},{s,-1,1}]]

Represent a given bivariate polynomial p(x,y) in terms of sum-of-squares polynomial (∑_(i=1))(((Subscript[g, i])(x,y))^2):

poly=2+4 x+6 (x^2)+4 (x^3)+(x^4)-2 y+(y^2);

The objective is to find Q such that p(x,y)=v(x,y).Q.v(x,y), where v(x,y) is a vector of monomials:

v={1,x,y,(x^2),(y^2),x y};

Construct the symmetric matrix Q:

n=Length[v];
mat = Array[Indexed[q, {##}] &, {n, n}];
Q = UpperTriangularize[mat] + Transpose[UpperTriangularize[mat, 1]];

Find the polynomial coefficients of p(x,y) and v(x,y).Q.v(x,y) and make sure they are equal:

coeffs =Map[#==0&,DeleteCases[Flatten[CoefficientList[v.Q.v-poly,{x,y}]],0]]

Find the elements of Q:

res=SemidefiniteOptimization[0,{coeffs,q (Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[n])]) 0},q]

The quadratic term v.Q.v = (v.L).((L^T)v)=g.g, where L is a lower-triangular matrix obtained from the Cholesky decomposition of Q:

L=CholeskyDecomposition[(Q/.res)+(10^(-8))IdentityMatrix[n]];
g=L.v;
Chop[Expand[Sum[((g[[i]])^2),{i,Length[g]}]],(10^(-7))]

Compare the sum-of-squares polynomial to the given polynomial:

Chop[Expand[Total[g^2]] - poly,(10^(-7))]

OpenersExampleSubsection9(Spacer1)Geometry Problems  (6)

Find the smallest disk centered at c=(x,y) of radius r that encloses a set of points:

SeedRandom[123];
pts = RandomReal[{-1,1},{100,2}];

For each point (Subscript[p, i])={(Subscript[x, i]),(Subscript[y, i])}, the constraint (Norm[c-(Subscript[p, i])])≤r, i=1,2,…,n must be satisfied. This constraint is equivalent to (({{r, (((c-(Subscript[p, i])))^Τ)}, {c-(Subscript[p, i]), r I}})) ⪰0, i=1,2,…,n. Use (RefLink[Inactive, paclet:ref/Inactive]) when forming the constraints:

ℐ=IdentityMatrix[2,SparseArray];
constraints=Table[(({{r, -c+pi}, {-c+pi, r ℐ}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeString)])0,{pi,pts}];

Find the enclosing disk by minimizing the radius r:

res=SemidefiniteOptimization[r,constraints,{r,c∈Vectors[2]}]

Visualize the enclosing region:

Graphics[{Circle[c,r]/.res, {Red,PointSize[0.02],Point[pts]}}]

The minimal area bounding disk can also be found using (RefLink[BoundingRegion, paclet:ref/BoundingRegion]):

BoundingRegion[pts,"MinDisk"]

Find the smallest ellipse parametrized as {x:(Norm[a.x+b])≤1} that encompasses a set of points by minimizing the area:

pts = ;

For each point (Subscript[p, i])={(Subscript[x, i]),(Subscript[y, i])}, the constraint (Norm[a.(Subscript[p, i])+b ])≤1, i=1,2,…,n must be satisfied:

constraints=Table[Norm[a.pi+b]≤1,{pi,pts}];

The area is proportional to det((a^(-1))). Applying the monotone function (RefLink[Log, paclet:ref/Log]), the function to minimize is -log (det(a))=-log ((Subscript[Π, i=1])(Subscript[λ, i]))=-(∑_(i=1))log((Subscript[λ, i])). This in turn is equivalent to minimizing -Tr(a):

res=SemidefiniteOptimization[-Tr[a],{constraints, a(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0},{a,b∈Vectors[2]}]

Convert the parameterized ellipse into the explicit form (x-(Subscript[x, 0])).(Σ^(-1)).(x-(Subscript[x, 0]))≤1:

{Σ,x0}={Inverse[Transpose[a].a /.res], -LinearSolve[a,b]/.res};
minAreaEllipse=Ellipsoid[x0,Σ]

Graphics[{{Opacity[0.5],Red,Ellipsoid[x0,Σ]},{PointSize[0.025],Point[pts]}}]

A bounding ellipse, not necessarily minimal area, can be found using (RefLink[BoundingRegion, paclet:ref/BoundingRegion]):

fastEllipse=BoundingRegion[pts,"FastEllipse"]

The optimal ellipse has a smaller area:

{Area[minAreaEllipse],Area[fastEllipse]}

Graphics[{{Lighter[Blue,0.6],fastEllipse},{Opacity[0.3],Red,minAreaEllipse}}, ImageSize->Tiny]

Find the smallest ellipsoid parametrized as {x:(Norm[a.x+b])≤1} that encompasses a set of points in 3D by minimizing the volume:

pts = ;

For each point (Subscript[p, i])={(Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])}, the constraint (Norm[a.(Subscript[p, i])+b ])≤1, i=1,2,…,n must be satisfied:

constraints=Table[Norm[a.pi+b]≤1,{pi,pts}];

Minimizing the volume is equivalent to minimizing -log (det(a)), which is equivalent to minimizing -Tr(a):

res=SemidefiniteOptimization[-Tr[a],{constraints, a(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[3])])0},{a,b∈Vectors[3]}]

Convert the parameterized ellipse into the explicit form (x-(Subscript[x, 0])).(Σ^(-1)).(x-(Subscript[x, 0]))≤1:

{Σ,x0}={Inverse[Transpose[a].a /.res], -LinearSolve[a,b]/.res};
minVolEllipsoid=Ellipsoid[x0,Σ]

Graphics3D[{{Opacity[0.3],Red,Ellipsoid[x0,Σ]},{PointSize[0.025],Point[pts]}}]

A bounding ellipsoid, not necessarily minimum volume, can also be found using (RefLink[BoundingRegion, paclet:ref/BoundingRegion]):

fastEllipsoid=BoundingRegion[pts,"FastEllipsoid"]

{Volume[minVolEllipsoid],Volume[fastEllipsoid]}

Find the maximum area ellipse parametrized as {x=C.y+d :(Norm[y])≤1} that can be fitted into a convex polygon:

cmesh=ConvexHullMesh[];

Each segment of the convex polygon can be represented as intersections of half-planes (Subscript[a, i]).x≤(Subscript[b, i]),i=1,2,…,n:

segs=MeshPrimitives[cmesh,1][[All,1]];
a=segs/.{p_,q_}:> {-1,1}*Reverse[p-q];
b=MapThread[#1.#2&,{a,segs[[All,1]]}];

Applying the parametrization to the half-planes gives sup{((Subscript[a, i]).C).y:(Norm[y])≤1}+(Subscript[a, i]).d≤(Subscript[b, i]). The term sup{((Subscript[a, i]).C).y:(Norm[y])≤1}=sup{(Norm[(Subscript[a, i]).C])(Norm[y])Cos(θ):(Norm[y])≤1}=(Norm[(Subscript[a, i]).C]). Thus, the constraints are  (Norm[(Subscript[a, i]).C])+(Subscript[a, i]).d≤ (Subscript[b, i]), which is equivalent to (({{(Subscript[b, i])-(Subscript[a, i]).d, ((Subscript[a, i]).C)ᵀ}, {(Subscript[a, i]).C,  ((Subscript[b, i])-(Subscript[a, i]).d) I}})) ⪰0, i=1,2,…,n:

ℐ=IdentityMatrix[2];
constraints=MapThread[(({{-#1.d+#2, c.#1}, {c.#1, (#2-#1.d) ℐ}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeString)])0&,{a,b}];

Minimizing the area is equivalent to minimizing -log (det(C)), which is equivalent to minimizing -Tr(C):

res=SemidefiniteOptimization[-Tr[c],{constraints,c(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0},{c,d∈Vectors[2]}]

Convert the parameterized ellipse into the explicit form as (x-d).((((C^Τ).C))^(-1)).(x-d)≤1:

inscribedEllipse=Ellipsoid[d,Transpose[c].c] /.res

Show[cmesh,Graphics[{LightRed,inscribedEllipse}]]

Find the center p and radius r of a disk given by (x^Τ).x-2x.p+γ≤0 that encloses three ellipses of the form x.(Subscript[a, i]).x+2(Subscript[b, i]).x+(Subscript[c, i])≤0:

ellipse[1]={{{0.13,0.02},{0.02,0.62}},{.2,.4},.1}; 
ellipse[2]={{{.18,-.16},{-.16,.35}},{.1,.2},0.1}; 
ellipse[3] = {{{0.7,0.68},{0.68,0.94}},{0.4,.3},0.12};

Using S-procedure, it can be shown that the disk contains the ellipses iff  (({{γ, -p}, {-(p^Τ), I}})) \[VectorLessEqual](Subscript[τ, k])(({{c, b}, {(b^Τ), a}})), (Subscript[τ, k])≥0:

ℐ=IdentityMatrix[2];
ellipseConstraint=Table[{a,b,c}=ellipse[i];(({{γ, -p}, {-p, ℐ}}))(Underscript[\[VectorLessEqual], (SemidefiniteConeList[3])])Indexed[τ,i](({{c, b}, {b, a}})),{i,3}];

The goal is to minimize the radius r given by r=√(p^T).p-γ. Using auxiliary variable t=(p^Τ).p-γ, the objective is to minimize t such that (r^2)≤t, which can be written as (({{t+γ, p}, {(p^Τ), I}}))\[VectorGreaterEqual]0:

auxConstraint=(({{t+γ, p}, {p, ℐ}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[3])])0;

Find the center and radius of the disk:

res=SemidefiniteOptimization[t,{auxConstraint,ellipseConstraint,τ\[VectorGreaterEqual]0},{t,γ,p∈Vectors[2],τ∈Vectors[3]}]

The disk is given by:

enclosingDisk=Disk[p,Sqrt[t]]/.res

Convert the quadratic form of the ellipse to the explicit form (x-(Subscript[x, 0]))(Σ^(-1))(x-(Subscript[x, 0]))≤1:

explicitEllipses=Table[{a,b,c}=ellipse[i];
With[{ai=Inverse[a]},Ellipsoid[-ai.b,ai(b.ai.b-c)]],{i,3}];

Visualize the result:

Graphics[{{Lighter[Red,0.8],enclosingDisk},{Opacity[0.5],Blue,explicitEllipses}},Frame->True]

Test whether an ellipsoid is a subset of another ellipsoid of the form v.(Subscript[a, i]).v+2(Subscript[b, i]).v+(Subscript[c, i])≤0:

{a1,b1,c1}={{{0.13,0.2},{0.2,0.62}},{.2,.4},.01}; 
{a2,b2,c2}={{{0.7,0.68},{0.68,0.94}},{0.35,.42},0.12};

Using S-procedure, it can be shown that ellipse 2 is a subset of ellipse 1 iff  (({{(Subscript[a, 1]), (Subscript[b, 1])}, {((Subscript[b, 1])^Τ), (Subscript[c, 1])}})) \[VectorLessEqual]τ (({{(Subscript[a, 2]), (Subscript[b, 2])}, {((Subscript[b, 2])^Τ), (Subscript[c, 2])}})), τ≥0:

subsetConstraint=(({{a1, b1}, {b1, c1}}))(Underscript[\[VectorLessEqual], (SemidefiniteConeList[3])])τ(({{a2, b2}, {b2, c2}}));

Check if the condition is satisfied:

SemidefiniteOptimization[0,{subsetConstraint,τ≥0},τ]

Convert the ellipsoids into explicit form and confirm that ellipse 2 is within ellipse 1:

ellipse1=With[{ai=Inverse[a1]},Ellipsoid[-ai.b1,ai(b1.ai.b1-c1)]];
ellipse2=With[{ai=Inverse[a2]},Ellipsoid[-ai.b2,ai(b2.ai.b2-c2)]];

Graphics[{{Red,Opacity[0.5],ellipse1},{Blue,Opacity[0.5],ellipse2}}]

Move ellipsoid 2 such that it overlaps with ellipsoid 1:

{a3,b3,c3}={{{0.7,0.68},{0.68,0.94}},{0.45,.42},0.12};
ellipse3=With[{ai=Inverse[a3]},Ellipsoid[-ai.b3,ai(b3.ai.b3-c3)]];

A test now shows that the problem is infeasible, indicating that ellipsoid 2 is not a subset of ellipsoid 1:

SemidefiniteOptimization[0,{(({{a1, b1}, {b1, c1}}))(Underscript[\[VectorLessEqual], (SemidefiniteConeList[3])])τ(({{a3, b3}, {b3, c3}})),τ≥0},τ];

Graphics[{{Red,Opacity[0.5],ellipse1},{Blue,Opacity[0.5],ellipse3}}]

OpenersExampleSubsection10(Spacer1)Classification Problems  (2)

Find an ellipse that separates two groups of points P={(Subscript[p, 1]),(Subscript[p, 2]),..(Subscript[p, n])} and Q={(Subscript[q, 1]),(Subscript[q, 2]),..(Subscript[q, m])}:

{set1,set2}=;
dataplot = ListPlot[{set1,set2},PlotMarkers->{Automatic,PointSize[0.025]}]

For separation, set 1 must satisfy ((Subscript[p, i])^Τ).A.(Subscript[p, i])+b.(Subscript[p, i])+c≥0 and set 2 must satisfy ((Subscript[q, j])^Τ).A.(Subscript[q, j])+b.(Subscript[q, j])+c≤0:

ellipsefun = Function[{x,y},c+(x^2) (IndexedDefault[a, 1, 1])+2 x y (IndexedDefault[a, 1, 2])+(y^2) (IndexedDefault[a, 2, 2])+x (IndexedDefault[b, 1])+y (IndexedDefault[b, 2])];
constraints = {Apply[ellipsefun,Transpose[set1]]\[VectorLessEqual]0,
Apply[ellipsefun,Transpose[set2]]\[VectorGreaterEqual]0,a(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0};

Find the coefficients of the separating ellipsoid:

res=SemidefiniteOptimization[0,constraints, {a,b∈Vectors[2,Reals],c}]

Visualize the result:

Show[Region[With[{v={x,y}},ImplicitRegion@@{v.a.v+b.v+c≤0/.res,v}]],dataplot]

Find an ellipse that is as close as possible to a circle that separates two groups of points P={(Subscript[p, 1]),(Subscript[p, 2]),..(Subscript[p, n])} and Q={(Subscript[q, 1]),(Subscript[q, 2]),..(Subscript[q, m])}:

{set1,set2}=;
dataplot = ListPlot[{set1,set2},PlotMarkers->{Automatic,PointSize[0.025]}]

For separation, set 1 must satisfy ((Subscript[p, i])^Τ).A.(Subscript[p, i])+b.(Subscript[p, i])+c≥0 and set 2 must satisfy ((Subscript[q, j])^Τ).A.(Subscript[q, j])+b.(Subscript[q, j])+c≤0:

ellipsefun = Function[{x,y},c+(x^2) (IndexedDefault[a, 1, 1])+2 x y (IndexedDefault[a, 1, 2])+(y^2) (IndexedDefault[a, 2, 2])+x (IndexedDefault[b, 1])+y (IndexedDefault[b, 2])];
constraints = {Apply[ellipsefun,Transpose[set1]]\[VectorLessEqual]0,
Apply[ellipsefun,Transpose[set2]]\[VectorGreaterEqual]0};

For the ellipsoid to be as close as possible to a circle, the constraint I\[VectorLessEqual]a\[VectorLessEqual]γ I, γ≥1:

ℐ=IdentityMatrix[2];
circleConstraints={ℐ(Underscript[\[VectorLessEqual], (SemidefiniteConeList[2])])a (Underscript[\[VectorLessEqual], (SemidefiniteConeList[2])])γ ℐ,γ≥1};

Find the coefficients of the separating ellipsoid by minimizing γ:

res=SemidefiniteOptimization[γ,{constraints,circleConstraints,a(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0}, {a,b∈Vectors[2,Reals],c,γ}]

Visualize the result:

Show[Region[With[{v={x,y}},ImplicitRegion@@{v.a.v+b.v+c≤0/.res,v}]],dataplot]

OpenersExampleSubsection11(Spacer1)Graph Problems  (2)

The Lovász number, computable using semidefinite optimization, is used as a bound for hard compute graph invariants:

{ShannonC, LovaszNum, cliqueNum, cromaticNum}=GraphData[{"Cycle", 5},{ "ShannonCapacity","LovaszNumber", "CliqueNumber", "ChromaticNumber"}]

The Lovász number is an upper bound for the Shannon capacity of a graph:

ShannonC≤ LovaszNum

According to the Lovász sandwich theorem:

cliqueNum≤ LovaszNum≤ cromaticNum

The Lovász number ϑ(G) for a graph G=(V,E) is given by ϑ(G) =(max_Y)(∑_(i,j∈V))(Subscript[Y, i j]), where Y\[VectorGreaterEqual]0, Tr(Y)=1 and (Subscript[Y, i j])\[VectorGreaterEqual]0 for i<->j ∈ E. It can be written in a dual semidefinite form: (max_Y)-Tr(J.Y), (Subscript[J, i j])=-1, subject to Tr(I.Y)=1 and Tr((Subscript[S, i j]).Y)=0, (Subscript[S, i j])=(Subscript[S, j i])=1,i<->j ∈ E , and 0 elsewhere:

LovaszTheta[g_?GraphQ]:=Module[{n=VertexCount[g],m=AdjacencyMatrix[g],amats,adjvs},
adjvs=UpperTriangularize[m]["NonzeroPositions"];amats=Join[{ConstantArray[-1,{n,n}],IdentityMatrix[n]},Apply[Function[SparseArray[{{#1,#2}->1,{#2,#1}->1},{n,n}]],adjvs,1]];
SemidefiniteOptimization[Flatten[{1, ConstantArray[0,Length[adjvs]]}],amats,"DualMaximumValue"]]

Compare to the exact Lovász number values from (RefLink[GraphData, paclet:ref/GraphData]):

{G=CycleGraph[5], N[GraphData[{"Cycle", 5},"LovaszNumber"]],LovaszTheta[G] }

{G=GraphData[{"Complete",4}], GraphData[{"Complete",4},"LovaszNumber"], LovaszTheta[G]}

{G=GraphData["PetersenGraph"],GraphData["PetersenGraph","LovaszNumber"], LovaszTheta[G] }

Find the approximate result for ϑ(G) when the exact result is not available:

{G=GraphData[{"Rook",{2,4}}],GraphData[{"Rook",{2,4}},"LovaszNumber"], LovaszTheta[G]}

Find ϑ(G) for a random graph:

{G=RandomGraph[{300,500}], LovaszTheta[G]}

The max-cut problem determines a subset S of the vertices V of a graph, for which the sum of the weights (Subscript[w, i j]) of the edges that cross from S to its complement V \ S is maximized. Let (Subscript[x, k])=1 for k∈S and (Subscript[x, k])=-1 for k∈V \ S. Maximize (Underoverscript[(1/4)∑, (i,j=1), n])(Subscript[w, i j ])(1-(Subscript[x, i])(Subscript[x, j]))=Tr ((1/4)  L · (x·xᵀ)), where (Subscript[x, k])∈{-1, 1}, k=1,…, n and L is the Laplacian matrix of the graph:

laplacianMatrix[G_?GraphQ]:=(DiagonalMatrix[Total[#,{1}]]-#)&[WeightedAdjacencyMatrix[G]]

For smaller cases, the max-cut problem can be solved exactly, but this is impractical for larger graphs since in general the problem has NP-complete complexity:

maxCutExact[G_?GraphQ]:=Module[{L=laplacianMatrix[G],n,x},
n=Length[First[L]];
x=Array["var",n];
MaxValue[{1/4x.L.x,x^2==ConstantArray[1,n]},x,Integers]
]

The problem minimizes (1/4) Tr[L·Y], where Y=x·(Transpose[x]) is a symmetric rank-1 positive semidefinite matrix, with (Subscript[Y, i i])=1 for each i, equivalent to Tr[(Subscript[E, i i])·Y]=1, where (Subscript[E, ii]) is the matrix with 1 at the i(^th) diagonal position and 0 everywhere else. To make the solution practical, solve a relaxed problem where the rank-1 condition is eliminated. For such Y, a cut is constructed by randomized rounding: decompose Y=U (Transpose[U]), let ρ be a uniformly distributed random vector of the unit norm and let (Subscript[x, k])=1⟺ρ·(Subscript[U, k])>0. For demonstration, a function is defined that shows the relaxed value, the rounded value and the graph with the vertices in S shown as red:

maxCutRelaxed[G_?GraphQ, trials_:10]:=Module[{L=laplacianMatrix[G],Y,U,n,sdpmax,x},
n=Length[First[L]];{Y,sdpmax}=SemidefiniteOptimization[ConstantArray[1,n], Prepend[Table[SparseArray[{i,i}->1.,{n,n},0.], {i,n}], -L/4],{"DualMaximizer","DualMaximumValue"}];
U=CholeskyDecomposition[Y];
x=First[MaximalBy[Table[Sign[RandomPoint[Sphere[n]].U],trials],1/4#.L.#&]];
{sdpmax,1/4x.L.x, HighlightGraph[G,Pick[Range[n],x,1]]}
]

Find an approximate max cut using the previously shown procedure, and compare with the exact result:

triang=Graph[{1\[UndirectedEdge]2,2\[UndirectedEdge]3,3\[UndirectedEdge]1},EdgeWeight->{2,3,4}, VertexSize->Small,VertexLabels->Placed["Name",Center],EdgeLabels->"EdgeWeight"];
{maxCutExact[triang],maxCutRelaxed[triang]}

Find the max cut for a grid graph:

grid=GridGraph[{3,4},PlotLabel->"(Subscript[G, 3, 4])",VertexLabels->"Name",VertexSize->Medium];
{maxCutExact[grid],maxCutRelaxed[grid]}

Find the max cut for a random graph:

g1=RandomGraph[{15,24}];
{maxCutExact[g1], maxCutRelaxed[g1]}

Compare timings for the relaxed and exact algorithms for a Peterson graph:

pet=PetersenGraph[8,2,VertexSize->Medium];
AbsoluteTiming[maxCutExact[pet]]

AbsoluteTiming[maxCutRelaxed[pet]]

OpenersExampleSubsection12(Spacer1)Control & Dynamic Systems Problems  (3)

Show that a linear dynamical system x'(t)=A.x(t) will be asymptotically stable for any initial condition. The system is said to be stable iff there exists a positive definite matrix P such that V'(t)<0, where V(x(t))=x(t)ᵀ.P.x(t) is called the Lyapunov function:

A={{-1,2,0},{-3,-4,1},{0,0,-2}};

Differentiating the Lyapunov function V(t) gives V'(t)=x(t).(Aᵀ.P+P.A).x(t)<0. Therefore, the stability conditions are Aᵀ.P+P.A\[VectorLess]0, P\[VectorGreater]0:

stabilityConditions={A\[Transpose].P+P.A(Underscript[\[VectorLess], (SemidefiniteConeList[3])])0,P(Underscript[\[VectorGreater], (SemidefiniteConeList[3])])0};

Find a matrix P:

res=SemidefiniteOptimization[0,stabilityConditions,P]

The eigenvalues of Aᵀ.P+P.A are all negative, making the matrix negative definite, which proves stability:

Eigenvalues[A\[Transpose].P+P.A/.res]

Since the analytic solution to the system is x(t)=(e^(A.t))x(0), numerically verify that any the system will go to zero for any initial condition. Take t=20 for this simulation:

ics=RandomReal[{-100,100},{3,100000}];
Norm[MatrixExp[A*20.].ics,"Frobenius"]

Find a controller u(t)=Q.x(t), such that the closed-loop system x'(t)=A.x(t)+B.u(t) is stable:

A={{2,-0.1},{-0.5,2}};
B={{1,1,0},{0,1,1}};

Using Lyapunov stability theorem, the objective is to find a matrices Q,P such that the stability constraints (A+B.Q).P + P.(A+B.Q)ᵀ\[VectorLess]0,P\[VectorGreater]0 are satisfied. Letting Z=Q.P, the first constraint becomes a proper semidefinite constraint A.P+B.Z+P.Aᵀ+(B.Z)ᵀ\[VectorLess]0:

stabilityConditions={A.P+B.Z+P.A\[Transpose]+(B.Z)\[Transpose](Underscript[\[VectorLess], (SemidefiniteConeList[2])])0,P(Underscript[\[VectorGreater], (SemidefiniteConeList[2])])0};

Find the matrices P,Z:

res=SemidefiniteOptimization[0,stabilityConditions,{P,Z∈Matrices[{3,2}]}]

The control matrix can be computed as Q=Z.(P^(-1)):

Q=Z.Inverse[P] /. res;

The closed-loop system is stable if the real parts of the eigenvalues of A+B.Q are negative:

Eigenvalues[A+B.Q]

Perform a numerical run to see that the system is stable:

sol=NDSolveValue[{x'[t]==(A+B.Q).x[t],x[0]=={1,1}},x,{t,0,10}];
Plot[sol[t],{t,0,10}, PlotRange->All]

Find a Lyapunov function of the dynamical system x'(t)=-y-(3/2)((x(t))^2)-(1/2)((x(t))^3),y'(t)=3x(t)-y(t):

f={-y-(3/2)x^2-(1/2)x^3,3x-y};

The objective is to find P such that V(x,y)=v(x,y).P.v(x,y), where v(x,y) is a vector of monomials:

v={x,y,x^2,y^2,x y};

Construct the matrix P:

n=Length[v];
mat = Array[Indexed[p, {##}] &, {n, n}];
P = UpperTriangularize[mat] + Transpose[UpperTriangularize[mat, 1]];

For stability, V(t)≥0,V'(t)≤0:

V=v.P.v;
Vdot=Expand[D[V,{{x,y}}].f];

Match the coefficients such that they are all positive for V(t) and negative for V'(t):

coeffs1 =Map[#>=0&,DeleteCases[Flatten[CoefficientList[V-0.(x^2+y^2),{x,y}]],0]];
coeffs2 = Map[#<=0&,DeleteCases[Flatten[CoefficientList[Vdot,{x,y}]],0]];

Find the matrix P:

res=SemidefiniteOptimization[0,{coeffs1,coeffs2, p (Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[6])])0},p ];

The Lyapunov function is given by:

lyap=Expand@Chop[Simplify[V/.res],10^-9]

Visualize the Lyapunov function. The minimal location of the function matches with the location of the attractor:

Show[StreamPlot[f,{x,-5,5},{y,-5,5},StreamStyle->"Segment"], ContourPlot[lyap,{x,-5,5},{y,-5,5},],AspectRatio->1/GoldenRatio, ImageSize->350]

OpenersExampleSection4(Spacer1)Properties & Relations  (8)

(RefLink[SemidefiniteOptimization, paclet:ref/SemidefiniteOptimization]) gives the global minimum of the objective function:

res=SemidefiniteOptimization[2(Subscript[x, 1])+3(Subscript[x, 2]),(({{(Subscript[x, 1]), 1}, {1, (Subscript[x, 2])}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,{(Subscript[x, 1]),(Subscript[x, 2])}]

Plot the objective function with the minimum value over the feasible region:

Show[Plot3D[2(Subscript[x, 1])+3(Subscript[x, 2]),{(Subscript[x, 1]),(Subscript[x, 2])}∈], ]

(RefLink[Minimize, paclet:ref/Minimize]) gives global exact results for semidefinite problems:

Minimize[{2(Subscript[x, 1])+3(Subscript[x, 2]),(({{(Subscript[x, 1]), 1}, {1, (Subscript[x, 2])}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0},{(Subscript[x, 1]),(Subscript[x, 2])}]

%//N

Compare to (RefLink[SemidefiniteOptimization, paclet:ref/SemidefiniteOptimization]):

SemidefiniteOptimization[2(Subscript[x, 1])+3(Subscript[x, 2]),(({{(Subscript[x, 1]), 1}, {1, (Subscript[x, 2])}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,{(Subscript[x, 1]),(Subscript[x, 2])}, {"PrimalMinimumValue","PrimalMinimizerRules"}]

(RefLink[NMinimize, paclet:ref/NMinimize]) can be used to get inexact results using global methods:

NMinimize[{2(Subscript[x, 1])+3(Subscript[x, 2]),(({{(Subscript[x, 1]), 1}, {1, (Subscript[x, 2])}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0},{(Subscript[x, 1]),(Subscript[x, 2])}]

(RefLink[FindMinimum, paclet:ref/FindMinimum]) can be used to obtain inexact results using local methods:

FindMinimum[{2x+3y,(({{x, 1}, {1, y}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0},{x,y}]

(RefLink[ConicOptimization, paclet:ref/ConicOptimization]) is more general than (RefLink[SemidefiniteOptimization, paclet:ref/SemidefiniteOptimization]):

ConicOptimization[2(Subscript[x, 1])+3(Subscript[x, 2]),(({{(Subscript[x, 1]), 1}, {1, (Subscript[x, 2])}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,{(Subscript[x, 1]),(Subscript[x, 2])}]

SemidefiniteOptimization[2(Subscript[x, 1])+3(Subscript[x, 2]),(({{(Subscript[x, 1]), 1}, {1, (Subscript[x, 2])}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,{(Subscript[x, 1]),(Subscript[x, 2])}]

(RefLink[SecondOrderConeOptimization, paclet:ref/SecondOrderConeOptimization]) is a special case of (RefLink[SemidefiniteOptimization, paclet:ref/SemidefiniteOptimization]):

SecondOrderConeOptimization[10x+11 y,{x+y≥1,{x,y}∈Disk[]},{x,y}]

SemidefiniteOptimization[10x+11 y,{x+y≥1,{x,y}∈Disk[]},{x,y}]

(RefLink[QuadraticOptimization, paclet:ref/QuadraticOptimization]) is a special case of (RefLink[SemidefiniteOptimization, paclet:ref/SemidefiniteOptimization]):

obj=2(x^2)+20(y^2)+6x y+5x;
cons={-x+y≥2,y≥0};

QuadraticOptimization[obj,cons,{x,y}]

Use auxiliary variable t and minimize t with additional constraint 2(x^2)+20(y^2)+6x y+5x≤t:

SemidefiniteOptimization[t,{obj≤t,cons},{x,y,t}]

(RefLink[LinearOptimization, paclet:ref/LinearOptimization]) is a special case of (RefLink[SemidefiniteOptimization, paclet:ref/SemidefiniteOptimization]):

LinearOptimization[10.x+11.y,{x,y}∈RegularPolygon[4],{x,y}]

SemidefiniteOptimization[10.x+11.y,{x,y}∈RegularPolygon[4],{x,y}]

OpenersExampleSection5(Spacer1)Possible Issues  (3)

The constraints at the optimal point are satisfied up to some tolerance:

res=SemidefiniteOptimization[2(Subscript[x, 1])+3(Subscript[x, 2]),(({{(Subscript[x, 1]), 1}, {1, (Subscript[x, 2])}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,{(Subscript[x, 1]),(Subscript[x, 2])}]

(({{(Subscript[x, 1]), 1}, {1, (Subscript[x, 2])}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0/.res

With default options, the constraint violation tolerance is (10^(-8)):

Eigenvalues[(({{(Subscript[x, 1]), 1}, {1, (Subscript[x, 2])}}))/.res]

The minimum value of an empty set or infeasible problem is defined to be ∞:

SemidefiniteOptimization[x,{x≤0,x≥1},x,"PrimalMinimumValue"]

The minimizer is (RefLink[Indeterminate, paclet:ref/Indeterminate]):

SemidefiniteOptimization[x,{x≤0,x≥1},x]

The minimum value for an unbounded set or unbounded problem is -∞:

SemidefiniteOptimization[-(x+y),x+y≥1,{x,y},"PrimalMinimumValue"]

The minimizer is (RefLink[Indeterminate, paclet:ref/Indeterminate]):

SemidefiniteOptimization[-(x+y),x+y≥1,{x,y}]

({{DocumentationFooterIconsRelatedFunction.png, ({{See Also}, {(RefLink[QuadraticOptimization, paclet:ref/QuadraticOptimization])  ▪  (RefLink[SecondOrderConeOptimization, paclet:ref/SecondOrderConeOptimization])  ▪  (RefLink[ConicOptimization, paclet:ref/ConicOptimization])  ▪  (RefLink[FindMinimum, paclet:ref/FindMinimum])  ▪  (RefLink[NMinimize, paclet:ref/NMinimize])  ▪  (RefLink[Minimize, paclet:ref/Minimize])}})}})

({{DocumentationFooterIconsRelatedGuide.png, ({{Related Guides}, {▪(RefLinkPlain[Optimization, paclet:guide/Optimization])}, {▪(RefLinkPlain[Convex Optimization, paclet:guide/ConvexOptimization])}})}})