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)=h^T(x)h(x) (or equivalently p(x)=v^T(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)-h^T(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)-v^T(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 h^T(x)h(x) and v^T(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)-h^T(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 v^T(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)²-xy+(1-y)² and (1-xy)²+xy+t(1+y)² 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) = v^T(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: