📄 sos.htm
字号:
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<html>
<head>
<meta http-equiv="Content-Language" content="en-us">
<title>YALMIP Example : Sum-of-squares decompositions</title>
<meta http-equiv="Content-Type" content="text/html; charset=windows-1251">
<meta content="Microsoft FrontPage 6.0" name="GENERATOR">
<meta name="ProgId" content="FrontPage.Editor.Document">
<link href="yalmip.css" type="text/css" rel="stylesheet">
<base target="_self">
</head>
<body leftMargin="0" topMargin="0">
<div align="left">
<table border="0" cellpadding="4" cellspacing="3" style="border-collapse: collapse" bordercolor="#000000" width="100%" align="left" height="100%">
<tr>
<td width="100%" align="left" height="100%" valign="top">
<h2>Sum of squares decompositions</h2>
<hr noShade SIZE="1">
<p>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
<strong>p(x)</strong> and tries to find a real-valued polynomial vector function
<strong>h(x)</strong> such
that <strong>p(x)=h<sup>T</sup>(x)h(x)</strong> (or equivalently <strong>p(x)=v<sup>T</sup>(x)Qv(x)</strong>,
where <b>Q</b> is positive semidefinite and <strong>v(x)</strong> is a vector of monomials), hence proving non-negativity of the polynomial <strong>p(x)</strong>.
Read more about sum of squares decompositions in
<a href="readmore.htm#PARRILO">[Parrilo]</a>.</p>
<h3>Options</h3>
<p>In the examples below, 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:</p>
<table border="1" cellspacing="1" style="border-collapse: collapse" width="100%" bordercolor="#000000" bgcolor="#FFFFE0">
<tr>
<td width="301"><code>sdpsettings('sos.model',[0|1|2])</code></td>
<td>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).<p>
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.</p> </td>
</tr>
<tr>
<td width="301"><code>sdpsettings('sos.newton',[0|1])</code></td>
<td>To reduce the size of the resulting SDP, a Newton polytope
reduction algorithm is applied by default. For efficiency, you must
have <a href="solvers.htm#cdd">CDDMEX</a> or <a href="solvers.htm#glpk">
GLPKMEX</a> installed.</td>
</tr>
<tr>
<td width="301"><code>sdpsettings('sos.congruence',[0|1])</code></td>
<td>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.</td>
</tr>
<tr>
<td width="301"><code>sdpsettings('sos.inconsistent',[0|1])</code></td>
<td>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.</td>
</tr>
<tr>
<td width="301"><code>sdpsettings('sos.extlp',[0|1])</code></td>
<td>When a kernel representation model is used, the SDP problem is
derived using the <a href="dual.htm#dualize">dualization</a>
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 <b>set(t>-100)</b>).</td>
</tr>
<tr>
<td width="301"><code>sdpsettings('sos.postprocess',[0|1])</code></td>
<td>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.</td>
</tr>
</table>
<h3>Simple sum of squares problems</h3>
<p>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 <code>yalmipdemo.m</code>). </p>
<p>The most improtant commands are <a href="reference.htm#sos">
sos</a> (to define a SOS constraint) and
<a href="reference.htm#solvesos">
solvesos</a> (to solve the problem).</p>
<table cellPadding="10" width="100%">
<tr>
<td class="xmpcode">
<pre>x = sdpvar(1,1);y = sdpvar(1,1);
p = (1+x)^4 + (1-y)^2;
F = set(sos(p));
solvesos(F);</pre>
</td>
</tr>
</table>
<p>The sum of squares decomposition is extracted with the command
<a href="reference.htm#sosd">
sosd</a>.</p>
<table cellPadding="10" width="100%">
<tr>
<td class="xmpcode">
<pre>h = sosd(F);
sdisplay(h)
<font color="#000000">
ans = </font></pre>
<pre><font color="#000000"> '-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'</font></pre>
</td>
</tr>
</table>
<p>To see if the decomposition was successful, we simply calculate <strong>
p(x)-h<sup>T</sup>(x)h(x) </strong>which should be 0. However, due to numerical errors,
the difference will not be zero. A useful command then is
<a href="reference.htm#clean">
clean</a>. Using
<a href="reference.htm#clean">
clean</a>, we remove all monomials with coefficients smaller than, e.g., 1e-6.</p>
<table cellPadding="10" width="100%">
<tr>
<td class="xmpcode">
<pre>clean(p-h'*h,1e-6)
<font color="#000000">
ans =</font></pre>
<pre><font color="#000000"> 0</font></pre>
</td>
</tr>
</table>
<p>The decomposition <strong>
p(x)-v<sup>T</sup>(x)Qv(x) </strong>can also be obtained easily.</p>
<table cellPadding="10" width="100%">
<tr>
<td class="xmpcode">
<pre>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)</pre>
<pre><font color="#000000"> ans =
0</font></pre>
</td>
</tr>
</table>
<p><font color="#FF0000">Note</font> : Even though <strong>
h<sup>T</sup>(x)h(x) </strong>and <strong>
v<sup>T</sup>(x)Qv(x)</strong> 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
<b>Q</b>
not is positive definite (often the case due to numerical issues in the
SDP solver). The decomposition is in theory defined as <strong>h(x)=chol(Q)v(x)</strong>
. If <strong>
Q</strong> is indefinite, YALMIP uses an approximation of the Cholesky
factor based on a singular value decomposition.</p>
<p>The quality of the decomposition can alternatively be evaluated using
<a href="reference.htm#checkset">checkset</a>. The value reported is the
largest coefficient in the polynomial <strong>
p(x)-h<sup>T</sup>(x)h(x)</strong></p>
<table cellPadding="10" width="100%" id="table1">
<tr>
<td class="xmpcode">
<pre>checkset(F)
<font color="#000000">++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
| ID| Constraint| Type| Primal residual|
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
| #1| Numeric value| SOS constraint (polynomial)| 7.3674e-012|
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++</font></pre>
<pre>e = checkset(F(is(F,'sos')))
<font color="#000000">e =</font></pre>
<pre><font color="#000000"> 7.3674e-012</font></pre>
</td>
</tr>
</table>
<p>It is very rare (except in contrived academic examples) that one finds a
decomposition with a positive definite <strong>
Q </strong>and zero mismatch<strong>
</strong>between <strong>
v<sup>T</sup>(x)Qv(x) </strong>and the polynomial <strong>
p(x)</strong> .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).</p>
<h3>Parameterized problems</h3>
<p>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).</p>
<table cellPadding="10" width="100%">
<tr>
<td class="xmpcode">
<pre>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)
<font color="#000000">
ans =</font></pre>
<pre><font color="#000000"> 0.75</font></pre>
</td>
</tr>
</table>
<p>YALMIP automatically declares all variables appearing in the objective
function and in non-SOS constraints as parametric variables. Hence, the
following code is equivalent.</p>
<table cellPadding="10" width="100%">
<tr>
<td class="xmpcode">
<pre>solvesos(F,-lower);
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -