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 semidefinite), 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 very large problems, you might
want to check out the dedicated software package
Gloptipoly
(the solution time will be the same, but the setup time might be reduced).
For the underlying theory of moment relaxations, the reader is referred to
[Lasserre].
Solving polynomial problems by relaxations
The following code calculates a lower bound on a concave quadratic
optimization problem. As you can see, the only difference compared to
solving the problem using a standard solver, such as
fmincon or
penbmi, is that we call solvemoment instead of
solvesdp.
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 (or several) 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 that was needed to extract the
solution.
[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).
assign([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|
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
Extracting sum of squares decompositions.
Moment based approaches and sum of squares
computations are directly linked through duality. Given the solution
to the moment problem, a sum of squares decomposition is easily derived
from dual variables. Consider a simple unconstrained problem.
sdpvar x1 x2
p = x1^4+x2^4-x1*x2+1
|
A sum of squares decomposition of p can be obtained by using
solvesos
[sol,v,Q] = solvesos(set(sos(p)));
sdisplay(clean(p-v{1}'*Q{1}*v{1},1e-6))
0
|
The decomposition can alternatively be computed from the moment
results. If we minimize the polynomial using moments,
solvemoment will automatically
extract t, v and Q in the decomposition p(x)-t=vTQv
[sol,x,moments,sosdec] = solvemoment([],p);
t = sosdec.t;
v = sosdec.v0;
Q = sosdec.Q0;
sdisplay(clean(p-t-v'*Q*v,1e-6))
0
|
In the more general constrained case, the polynomial multipliers for
the constraints will also be returned.
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);
assign([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.rceftol',1e-8);
[sol,xe,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.rceftol',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));
assign(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|
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
The
semidefinite programs rapidly grows when the number of variables and the
polynomial degree increase, so be careful when you model your problem.
The two
most useful practical tips when working semidefinite relaxations seem to
be to de-symmetrize your objective function, and compactify your
feasible region. These two tricks typically increase the likelihood that
you will be able to extract global solutions. By adding a perturbation
to the polynomial, symmetry is lost, which generically means that there
will not be infinitely many optima, and the extraction algorithm is more
likely to work. Most theory in moment relaxations assumes that the
feasible set is compact, and this is also affecting practical
performance. By adding redundant compactifying constraints, you
typically increase the likelihood of success. As an example, a simple
redundant constraint which often work well in practice is an upper bound on
the objective function based on a known feasible sub-optimal solution. |