Moment based relaxation of polynomials programs


YALMIP comes with a built-in module for polynomial programming using moment relaxations. This can be used for finding lower bounds on constrained polynomial programs (inequalities and equalities, element-wise and in the semidefinite cone), and to extract the corresponding optimizers. The implementation is entirely based on high-level YALMIP code, and can be somewhat inefficient for large problems (the inefficiency would then show in the setup of the problem, not actually solving the semidefinite resulting program). For larger problems, you might want to check out the dedicated software package Gloptipoly. For the underlying theory of moment relaxations, the reader is referred to [Lasserre].

The following code calculates a lower bound on a concave quadratic optmization problem.

sdpvar x1 x2 x3
obj = -2*x1+x2-x3;
F = set(x1*(4*x1-4*x2+4*x3-20)+x2*(2*x2-2*x3+9)+x3*(2*x3-13)+24>0);
F = F + set(4-(x1+x2+x3)>0);
F = F + set(6-(3*x2+x3)>0);
F = F + set(x1>0);
F = F + set(2-x1>0);
F = F + set(x2>0);
F = F + set(x3>0);
F = F + set(3-x3>0);
solvemoment(F,obj);
double(obj)

 ans =
     -6.0000

Notice that YALMIP does not recover variables by default, a fact showing up in the difference between lifted variables and actual nonlinear variables (lifted variables are the variables used in the semidefinite relaxation to model nonlinear variables) The lifted variables can be obtained by using the command relaxdouble . The quadratic constraint above is satisfied in the lifted variables, but not in the true variables, as the following code illustrates.

relaxdouble(x1*(4*x1-4*x2+4*x3-20)+x2*(2*x2-2*x3+9)+x3*(2*x3-13)+24)

 ans =

   23.2648
double(x1*(4*x1-4*x2+4*x3-20)+x2*(2*x2-2*x3+9)+x3*(2*x3-13)+24)

 ans =

  -2.0000

An tighter relaxation can be obtained by using a higher order relaxation (the lowest possible is used if it is not specified).

solvemoment(F,obj,[],2);
double(obj)

 ans =
     -5.6593

The obtained bound can be used iteratively to improve the bound by adding dynamically generated cuts.

solvemoment(F+set(obj>double(obj)),obj,[],2);
double(obj)

 ans =
     -5.3870
solvemoment(F+set(obj>double(obj)),obj,[],2);
double(obj)

 ans =

     -5.1270

The known true minima, -4, is found in the fourth order relaxation.

solvemoment(F,obj,[],4);
double(obj)

 ans =
     -4.0000

The true global minima is however not recovered with the lifted variables, as we can see if we check the current solution (still violates the nonlinear constraint).

checkset(F)

+++++++++++++++++++++++++++++++++++++++++++++++++++
| ID| Constraint   |         Type| Primal residual|
+++++++++++++++++++++++++++++++++++++++++++++++++++
| #1| Numeric value| Element-wise|        -0.88573|
| #2| Numeric value| Element-wise|           1.834|
| #3| Numeric value| Element-wise|           5.668|
| #4| Numeric value| Element-wise|           1.834|
| #5| Numeric value| Element-wise|         0.16599|
| #6| Numeric value| Element-wise|     2.0873e-006|
| #7| Numeric value| Element-wise|         0.33198|
| #8| Numeric value| Element-wise|           2.668|
+++++++++++++++++++++++++++++++++++++++++++++++++++

Extracting solutions

To extract a globally optimal solution, we need two output arguments. The first output is a diagnostic structure (standard solution structure from the semidefinite solver), the second output is the (hopefully) extracted globally optimal solutions and the third output is a data structure containing all data needed to extract variables.

[sol,x] = solvemoment(F,obj,[],4);
x{1}

ans =

    0.5000
    0.0000
    3.0000
x{2}

ans =

    2.0000
   -0.0000
    0.0000

We can easily check that these are globally optimal solutions since they reach the lower bound -4 and satisfy the constraints (up to numerical precision).

setsdpvar([x1;x2;x3],x{1});
double(obj)
ans =
   -4.0000
checkset(F)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|   ID|      Constraint|                       Type|   Primal residual|
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|   #1|   Numeric value|   Element-wise (quadratic)|      -1.1034e-006|
|   #2|   Numeric value|               Element-wise|               0.5|
|   #3|   Numeric value|               Element-wise|                 3|
|   #4|   Numeric value|               Element-wise|               0.5|
|   #5|   Numeric value|               Element-wise|               1.5|
|   #6|   Numeric value|               Element-wise|        5.956e-007|
|   #7|   Numeric value|               Element-wise|                 3|
|   #8|   Numeric value|               Element-wise|       4.6084e-007|
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Polynomial semidefinite constraints

Nonlinear semidefinite constraints can be added using exactly the same notation and syntax. The following example is taken from [D. Henrion, J. B. Lasserre].

sdpvar x1 x2
obj = -x1^2-x2^2;
F = set([1-4*x1*x2 x1;x1 4-x1^2-x2^2]);
[sol,x] = solvemoment(F,obj,[],2);
setsdpvar([x1;x2],x{1});
double(obj)
ans =
   -4.00003
checkset(F)
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|   ID|      Constraint|              Type|   Primal residual|
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|   #1|   Numeric value|   LMI (quadratic)|       -0.00034633|
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Advanced features

A number of advanced features are available. We just briefly introduce these here by a quick example where we refine the extracted solution using a couple of Newton steps on an algebraic systems defining the global solutions given the optimal moment matrices, and change the numerical tolerance in the extraction algorithm. Finally, we calculate some different global solutions using the optimal moment matrices. Please check the moment relaxation example in yalmipdemo for details.

sdpvar x1 x2
obj = -x1^2-x2^2;
F = set([1-4*x1*x2 x1;x1 4-x1^2-x2^2]);
ops = sdpsettings('moment.refine',5','moment.creftol',1e-8);
[sol,data] = solvemoment(F,obj,ops);
xopt1 = extractsolution(data,sdpsettings('moment.refine',0));
xopt2 = extractsolution(data,sdpsettings('moment.refine',1));
xopt3 = extractsolution(data,sdpsettings('moment.refine',10));
xopt4 = extractsolution(data,sdpsettings('moment.creftol',1e-3,'moment.refine',5));

The moment relaxation solver can also be called using a more standard YALMIP notation, by simply defining the solver as 'moment'.

sdpvar x1 x2
obj = -x1^2-x2^2;
F = set([1-4*x1*x2 x1;x1 4-x1^2-x2^2]);
sol = solvesdp(F,obj,sdpsettings('solver','moment','moment.order',2));
setsdpvar(sol.momentdata.x,sol.xoptimal{1});
double(obj)
ans =
   -4.00003
checkset(F)
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|   ID|      Constraint|              Type|   Primal residual|
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|   #1|   Numeric value|   LMI (quadratic)|       -0.00034633|
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++