Sum of squares decompositions
YALMIP has a built-in module for sum of squares calculations. In its most basic formulation, a
sum of squares (SOS) problem takes a polynomial
p(x) and tries to find a real-valued polynomial vector function
h(x) such
that p(x)=hT(x)h(x) (or equivalently p(x)=vT(x)Qv(x),
where Q is positive semidefinite and v(x) is a vector of monomials), hence proving non-negativity of the polynomial p(x).
Read more about standard sum of squares decompositions in
[Parrilo].
In addition to standard SOS decompositions, YALMIP also supports
linearly and nonlinearly
parameterized problems, decomposition of matrix
valued polynomials and computation of low rank
decompositions.
The examples below only introduce the basic features related to
sum-of-squares in YALMIP. For more features and details, run the example
sosex.m
Doing it your self the hard way
Before we introduce the efficiently implemented SOS module in YALMIP,
let us briefly mention how one could implement a SOS solver in
high-level YALMIP code. Of course, you would never use this approach,
but it might be useful to see the basic idea. Define a polynomial.
x = sdpvar(1,1);y = sdpvar(1,1);
p = (1+x)^4 + (1-y)^2;
|
Introduce a decomposition
v = monolist([x y],degree(p)/2);
Q = sdpvar(length(v));
p_sos = v'*Q*v;
|
The polynomials have to match, hence all coefficient in the polynomial
describing the difference of the two polynomials have to be zero.
F = set(coefficients(p-p_sos,[x y]) == 0) + set(Q > 0);
solvesdp(F)
|
The problem above is typically more efficiently solved when interpreted
in primal form, hence we dualize it.
F = set(coefficients(p-p_sos,[x y]) == 0) + set(Q > 0);
[F,obj] = dualize(F,[]);
solvesdp(F,-obj)
|
Adding parameterizations does not change the code significantly. Here is the
code to find a lower bound on the polynomial
sdpvar t
F = set(coefficients((p-t)-p_sos,[x y]) == 0) + set(Q > 0);
solvesdp(F,-t)
|
Matrix valued decompositions are a bit trickier, but still
straightforward.
sdpvar x y
P = [1+x^2 -x+y+x^2;-x+y+x^2 2*x^2-2*x*y+y^2];
m = size(P,1);
v = monolist([x y],degree(P)/2);
Q = sdpvar(length(v)*m);
R = kron(eye(m),v)'*Q*kron(eye(m),v)-P;
s = coefficients(R(find(triu(R))),[x y]);
solvesdp(set(Q > 0) + set(s==0));
sdisplay(clean(kron(eye(m),v)'*double(Q)*kron(eye(m),v),1e-6))
|
Once again, this is the basic idea behind the SOS module in YALMIP.
However, the implementation is far more efficient and uses various
approaches to reduce complexity, hence the approaches above should never be
used in practice.
Simple sum of squares problems
The following lines of code presents some typical manipulation when working
with SOS-calculations (a more detailed description is available if you run
the sum of squares example in yalmipdemo.m ).
The most important commands are
sos (to define a SOS constraint) and
solvesos (to solve the problem).
x = sdpvar(1,1);y = sdpvar(1,1);
p = (1+x)^4 + (1-y)^2;
F = set(sos(p));
solvesos(F);
|
The sum of squares decomposition is extracted with the command
sosd.
h = sosd(F);
sdisplay(h)
ans =
'-1.203-1.9465x+0.22975y-0.97325x^2'
'0.7435-0.45951x-0.97325y-0.22975x^2'
'0.0010977+0.00036589x+0.0010977y-0.0018294x^2'
|
To see if the decomposition was successful, we simply calculate
p(x)-hT(x)h(x) which should be 0. However, due to numerical errors,
the difference will not be zero. A useful command then is
clean. Using
clean, we remove all monomials with coefficients smaller than, e.g., 1e-6.
clean(p-h'*h,1e-6)
ans =
0
|
The decomposition
p(x)-vT(x)Qv(x) can also be obtained easily.
x = sdpvar(1,1);y = sdpvar(1,1);
p = (1+x)^4 + (1-y)^2;
F = set(sos(p));
[sol,v,Q] = solvesos(F);
clean(p-v{1}'*Q{1}*v{1},1e-6)
ans =
0
|
Note : Even though
hT(x)h(x) and
vT(x)Qv(x) should be the same in theory, they are
typically not. The reason is partly numerical issues with floating point
numbers, but more importantly due to the way YALMIP treats the case when
Q
not is positive definite (often the case due to numerical issues in the
SDP solver). The decomposition is in theory defined as h(x)=chol(Q)v(x)
. If
Q is indefinite, YALMIP uses an approximation of the Cholesky
factor based on a singular value decomposition.
The quality of the decomposition can alternatively be evaluated using
checkset. The value reported is the
largest coefficient in the polynomial
p(x)-hT(x)h(x)
checkset(F)
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
| ID| Constraint| Type| Primal residual|
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
| #1| Numeric value| SOS constraint (polynomial)| 7.3674e-012|
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
e = checkset(F(is(F,'sos')))
e =
7.3674e-012
|
It is very rare (except in contrived academic examples) that one finds a
decomposition with a positive definite
Q and zero mismatch
between
vT(x)Qv(x) and the polynomial
p(x) .The reason is simply that all SDP solver uses real
arithmetics. Hence, in principle, the decomposition has no theoretical
value as a certificate for non-negativity unless additional post-analysis
is performed (relating the size of the residual with the eigenvalues of
the Gramian).
Parameterized problems
The involved polynomials can be parameterized, and we can optimize the parameters
under the constraint that the polynomial is a sum of squares. As an example,
we can calculate a lower bound on a polynomial. The second argument can be used for an objective function to be
minimized. Here, we maximize the lower bound, i.e. minimize the negative
lower bound. The third argument is an options structure while the fourth argument is a vector
containing all parametric variables in the polynomial (in this example
we only have one variable, but several parametric variables can be defined
by simply concatenating them).
sdpvar x y lower
p = (1+x*y)^2-x*y+(1-y)^2;
F = set(sos(p-lower));
solvesos(F,-lower,[],lower);
double(lower)
ans =
0.75
|
YALMIP automatically declares all variables appearing in the objective
function and in non-SOS constraints as parametric variables. Hence, the
following code is equivalent.
solvesos(F,-lower);
double(lower)
ans =
0.75
|
Multiple SOS constraints can also be used. Consider the following problem of
finding the smallest possible t such that the polynomials t(1+xy)2-xy+(1-y)2
and (1-xy)2+xy+t(1+y)2 are both
sum of squares.
sdpvar x y t
p1 = t*(1+x*y)^2-x*y+(1-y)^2;
p2 = (1-x*y)^2+x*y+t*(1+y)^2;
F = set(sos(p1))+set(sos(p2));
solvesos(F,t);
double(t)
ans =
0.2500
sdisplay(sosd(F(1)))
ans =
'-1.102+0.95709y+0.14489xy'
'-0.18876-0.28978y+0.47855xy'
sdisplay(sosd(F(2)))
ans =
'-1.024-0.18051y+0.76622xy'
'-0.43143-0.26178y-0.63824xy'
'0.12382-0.38586y+0.074568xy'
|
If you have
parametric variables, bounding the feasible region typically
improves numerical behavior. Having lower bounds will
additionally decrease the size of the optimization problem
(variables bounded from below can be treated as translated cone
variables in dualization, which
is used by
solvesos).
One of the
most common mistake people make when using the sum of squares module
is to forget to declare some parametric variables. This will
typically lead to a (of-course erroneous) huge sum of squares
problem which easily freezes MATLAB and/or give strange error
messages (trivially infeasible, nonlinear parameterization, etc).
Make sure to initially run the module in verbose mode to see how
YALMIP characterizes the problem (most importantly to check the
number of parametric variables and independent variables). If you
use nonlinear operators (min, max, abs,...) on
parametric variables in your problem, it is recommended to always
declare the parametric variables.
When you use
a kernel representation (sos.model=1 and typically the case also
for sos.model=0 ), YALMIP will derive and solve a
problem which is related to the dual of your original problem. This
means that warnings about infeasibility after solving the SDP actually means
unbounded objective, and vice versa. If you use sos.model=2 ,
a primal problem is solved, and YALMIP error messages are directly related to
your problem.
The quality
of the SOS approximation is typically improved substantially if the
tolerance and precision options of the semidefinite solver is
decreased. As an example, having sedumi.eps
less than 10-10 when solving sum of squares problems is
typically recommended for anything but trivial problems. There is a
higher likelihood that the semidefinite solver will complain about
numerical problems in the end-phase, but the resulting solutions are
typically much better. This seem to be even more important in
parameterized problems.
To evaluate
the quality and success of the sum of squares decomposition, do not
forget to study the discrepancy between the decomposition and the
original polynomial. No problems in the
semidefinite solver is no guarantee for a successful decomposition
(due to numerical reasons, tolerances in the solvers etc.)
Polynomial parameterizations
A special feature of the sum of squares package in YALMIP is the
possibility to work with nonlinear SOS parameterizations, i.e. SOS problems
resulting in PMIs (polynomial matrix inequalities) instead of LMIs. The following piece of code solves a
nonlinear control synthesis problem using sum of squares. Note
that this requires the solver PENBMI.
clear all
yalmip('clear');
% States...
sdpvar x1 x2
x = [x1;x2];
% Non-quadratic Lyapunov z'Pz
z = [x1;x2;x1^2];
P = sdpvar(3,3);
V = z'*P*z;
% Non-linear state feedback
v = [x1;x2;x1^2];
K = sdpvar(1,3);
u = K*v;
% System x' = f(x)+Bu
f = [1.5*x1^2-0.5*x1^3-x2; 3*x1-x2];
B = [0;1];
% Closed loop system, u = Kv
fc = f+B*K*v;
% Stability and performance constraint dVdt < -x'x-u'u
% NOTE : This polynomial is bilinear in P and K
F = set(sos(-jacobian(V,x)*fc-x'*x-u'*u));
% P is positive definite, bound P and K for numerical reasons
F = F + set(P>0)+set(25>P(:)>-25)+set(25>K>-25);
% Minimize trace(P)
% Parametric variables P and K automatically detected
% by YALMIP since they are both constrained
solvesos(F,trace(P));
|
Matrix valued problems
In the same sense that the moment implementation in YALMIP supports
optimization over nonlinear semidefiniteness constraint
using moments, YALMIP also supports the dual SOS approach. Instead of computing a
decomposition
p(x) =
vT(x)Qv(x), the matrix SOS decomposition is P(x) = (kron(I,v(x))TQ(kron(I,v(x)).
sdpvar x1 x2
P = [1+x1^2 -x1+x2+x1^2;-x1+x2+x1^2 2*x1^2-2*x1*x2+x2^2];
[sol,v,Q] = solvesos(set(sos(P)));
sdisplay(v{1})
ans =
'1' '0'
'x2' '0'
'x1' '0'
'0' '1'
'0' 'x2'
'0' 'x1'
clean(v{1}'*Q{1}*v{1}-P,1e-8)
ans =
0 0
0 0
|
Of course, parameterized problems etc can also be solved with matrix
valued SOS constraints.
At the moment, YALMIP extends some of the reduction techniques from
the scalar case to exploit symmetry and structure of the polynomial
matrix, but there is room for obvious improvements. If you think you
might need this, make a feature request.
Keep in mind, that a simple scalarization can be more efficient in
many cases.
w = sdpvar(length(P),1);
[sol,v,Q] = solvesos(set(sos(w'*P*w)));
clean(v{1}'*Q{1}*v{1}-w'*P*w,1e-8)
ans =
0
|
This approach will however only prove positive semidefiniteness, it
will not provide a decomposition of the polynomial matrix.
Low-rank sum-of-squares (experimental)
By using the capabilities of the solver LMIRANK,
we can pose sum-of-squares problems where we search for decompositions
with few terms (low-rank Gramian). Consider the following problem where
a trace heuristic leads to an effective rank of 5, perhaps 6.
x = sdpvar(1,1);
y = sdpvar(1,1);
f = 200*(x^3 - 4*x)^2+200 * (y^3 - 4*y)^2+(y - x)*(y + x)*x*(x + 2)*(x*(x - 2)+2*(y^2 - 4));
g = 1 + x^2 + y^2;
p = f * g;
F = set(sos(p));
[sol,v,Q] = solvesos(F,[],sdpsettings('sos.traceobj',1));
eig(Q{1})
ans =
1.0e+003 *
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0124
0.3977
3.3972
3.4000
6.7972
|
We solve the problem using LMIRANK
instead, and aim for a rank less than or equal to 4. The desired rank is
specified easily in the sos construct.
x = sdpvar(1,1);
y = sdpvar(1,1);
f = 200*(x^3 - 4*x)^2+200 * (y^3 - 4*y)^2+(y - x)*(y + x)*x*(x + 2)*(x*(x - 2)+2*(y^2 - 4));
g = 1 + x^2 + y^2;
p = f * g;
F = set(sos(p,4));
[sol,v,Q] = solvesos(F,[],sdpsettings('lmirank.eps',0));
eig(Q{1})
ans =
1.0e+003 *
-0.0000
-0.0000
-0.0000
-0.0000
-0.0000
-0.0000
-0.0000
-0.0000
0.0000
0.0000
0.4634
4.2674
4.6584
7.1705
|
The resulting rank is indeed effectively 4. Note though that LMIRANK
works on the dual problem side, and can return slightly infeasible
solutions (in terms of positive definiteness.) Keep in mind that
sum-of-squares decompositions almost always be slightly wrong, in
one way or the other. If a dual (also called image) approach is used
(corresponding to sos.model=2), positive
definiteness may be violated, and if a primal approach (also called
kernel) approach is used (corresponding to
sos.model=1), there is typically a difference between the
polynomial and its decomposition. This simply due to the way SDP solvers
and floating point arithmetic work. See more in the example sosex.m
Remember that LMIRANK is a local
solver, hence there are no guarantees that it will find a low rank
solution even though one is known to exist. Moreover, note that LMIRANK
does not support objective functions. Options
In the examples above, we are mainly using the default settings when
solving the SOS problem, but there are a couple of options that can be
changed to alter the computations. The most important are:
sdpsettings('sos.model',[0|1|2]) |
The SDP formulation of a SOS problem is not unique but can be
done in several ways. YALMIP supports two version, here called image
and kernel representation. If you set the value to 1, a kernel
representation will be used, while 2 will result in an image
representation. If the option is set to the default value 0, YALMIP
will automatically select the representation (kernel by default).
The kernel representation will most often give a smaller and
numerically more robust semidefinite program, but cannot be used for
nonlinearly parameterized SOS programs (i.e. problems resulting in
BMIs etc) or problems with integrality
constraints on parametric variables. |
sdpsettings('sos.newton',[0|1]) |
To reduce the size of the resulting SDP, a Newton polytope
reduction algorithm is applied by default. For efficiency, you must
have CDDMEX or
GLPKMEX installed. |
sdpsettings('sos.congruence',[0|1]) |
A useful feature in YALMIP is the use of symmetry of the
polynomial to block-diagonalize the SDP. This can make a huge difference
for some SOS problems and is applied by default. |
sdpsettings('sos.inconsistent',[0|1]) |
This options can be used to further reduce the size of the SDP. It
is turned off by default since it typically not gives any major
reduction in problem size once the Newton polytope reduction has been
applied. In some situations, it might however be useful to use this
strategy instead of the linear programming based Newton reduction (it
cannot suffer from numerical problems and does not require any
efficient LP solver), and for some problems, it can reduce models that
are minimal in the Newton polytope sense, leading to a more
numerically robust solution of the resulting SDP. |
sdpsettings('sos.extlp',[0|1]) |
When a kernel representation model is used, the SDP problem is
derived using the dualization
function. For some problems, the strategy in the dualization may
affect the numerical conditioning of the problem. If you encounter
numerical problems, it can sometimes be resolved by setting this
option to 0. This will typically yield a slightly larger problem, but
can improve the numerical performance. Note that this option only is of use
if you have parametric variables with explicit non-zero lower bounds (constraints like set(t>-100)). |
sdpsettings('sos.postprocess',[0|1]) |
In practice, the SDP computations will never give a completely
correct decomposition (due to floating point numbers and in many
cases numerical problems in the solver). YALMIP can try to recover
from these problems by applying a heuristic post-processing
algorithm. This can in many cases improve the results. |
|