strictly sos certificate with sos.model = 2;

[Yalmiper]

In this tutorial, http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Blog.Strictly-feasible-SOS-solutions

I only need set sos.model = 2 in the last, feasibility problem, and if it returns Q with strictly positive eigenvalue, that means it is a strict sos certificate, is that correct? Because in my example, it seems setting sos.model in primary form in the original optimization gives better result.

So when setting sos.model =2, is the equality never violated?

Thank you!

[Johan Löfberg]

When you set model=2, it explicitly parametrizes Q such that p = v'*Q*v, hence if Q is psd, p is psd. In model=1, the condition that p = v'*Q*v would be encoded through equalities, which only will hold up to solver precision/tolerance. As you've noticed though, model=1 is typically better, but then you have to relate the errors in the equalities with the positive definiteness of Q to draw any conclusions, as described in the paper

[Yalmiper]

So if I understand correctly I can first set sos.model = 1, solve optimization problem (A) to get an optimal value (say value(b)); then solve a new feasibility problem, with the original constraints plus b<1.01*value(b) (I am minimizing b), and this time set sos.model =2 and solve feasibility problem (B). As long as the (B) is feasible then I get a strict sos certificate, or to put another way I validate the optimal solution of the problem (A) solved with sos.model = 1, within 1% of course.

I mean, in my problem, if I set sos.model = 2 in problem (A), i.e., directly minimize b with image representation, it gave a crappy result. But somehow the feasibility problem with sos.model, with 1% tolerance is indeed feasible.

[Johan Löfberg]

Yes. 

This sounds reasonable, when you push some variable to optimality, you will have a non-strict Q, and things possibly going towards large/small values, causing issues

[Yalmiper]

Hi Prof Lofberg,

I think you are talking mainly about polynomials without parameters?
My problem actually has one parameter in the polynomial and that is the one I am minimizing (i.e. "b" in previous post).
I realize that even when I set model=2, and the optimal Q is strict positive, the residue is still very bad. MOSEK report optimal and YALMIP suspect it is infeasible. I checked it manually, by replacing "b" in the left hand side as "value(b)" and some of the coefficients of the difference can be in the order of 1e-2. I guess this is a problem that is very insensitive to "b". Is there anyway to deal problems like this?

Thank you!

[Johan Löfberg]

You will have to supply an example

[Yalmiper]

Dear Professor Lofberg,

In fact, I now try set b to a fix value, now it is a polynomial without parameter, and the residue is bad.

Here is the code,

x = sdpvar(2,1);

lhs = 1357500.122+2.5588e+07*x(2)^2-3.0876e+07*x(2)^4-8.7309e+06*x(1)^2+3.0120e+07*x(1)^4-4.4903e+07*x(1)^6-1.8107e+08*x(2)^6+3.4204e+10*x(2)^8+3.0441e+09*x(2)^10-8.1294e+09*x(2)^12-9.5647e+09*x(2)^14+6.1053e+09*x(2)^16-9.0845e+08*x(2)^18+4.0129e+07*x(2)^20+3.0961e+07*x(1)^8-9.3812e+06*x(1)^10+1.0472e+06*x(1)^12-2.0073e+08*x(1)^2*x(2)^2-4.1267e+09*x(1)^2*x(2)^4+1.1283e+09*x(1)^4*x(2)^2+5.9385e+08*x(1)^2*x(2)^6+2.7951e+10*x(1)^4*x(2)^4-4.1687e+09*x(1)^6*x(2)^2+1.3196e+11*x(1)^2*x(2)^8-1.4124e+11*x(1)^4*x(2)^6-4.1472e+10*x(1)^6*x(2)^4-2.8262e+10*x(1)^2*x(2)^10-4.9590e+11*x(1)^2*x(2)^12+3.0234e+10*x(1)^2*x(2)^14+1.3769e+12*x(1)^2*x(2)^16-6.8587e+11*x(1)^2*x(2)^18+1.2039e+11*x(1)^2*x(2)^20-3.5009e+11*x(1)^4*x(2)^8+1.1319e+11*x(1)^4*x(2)^10+3.4424e+12*x(1)^4*x(2)^12-1.7411e+13*x(1)^4*x(2)^14+5.5342e+12*x(1)^4*x(2)^16-2.4789e+12*x(1)^4*x(2)^18+2.0473e+11*x(1)^6*x(2)^6+6.5129e+11*x(1)^6*x(2)^8-5.0836e+11*x(1)^6*x(2)^10-5.1200e+12*x(1)^6*x(2)^12+2.5276e+13*x(1)^6*x(2)^14+3.0111e+13*x(1)^6*x(2)^16+6.9852e+09*x(1)^8*x(2)^2+1.7716e+10*x(1)^8*x(2)^4-8.0929e+10*x(1)^8*x(2)^6-1.3639e+11*x(1)^8*x(2)^8+8.6015e+11*x(1)^8*x(2)^10+2.1871e+12*x(1)^8*x(2)^12-9.9913e+12*x(1)^8*x(2)^14-5.1169e+09*x(1)^10*x(2)^2+1.2670e+08*x(1)^10*x(2)^4-1.1501e+10*x(1)^10*x(2)^8-6.3172e+11*x(1)^10*x(2)^10+1.5642e+10*x(1)^10*x(2)^12+1.3491e+09*x(1)^12*x(2)^2+1.2929e+08*x(1)^12*x(2)^8+1.6656e+11*x(1)^12*x(2)^10;
C = sos(lhs);
ops = sdpsettings('solver','mosek','sos.model',2);
[sol,v,Q,res] = solvesos(C,[],ops);
rhs = v{1}'*Q{1}*v{1};
df = clean(lhs-rhs,1e-6);
sdisplay(df)

[Johan Löfberg]

You are doomed from the start here

> norm(getbase(lhs),inf)

ans =

   3.0111e+13

That can of numerics is impossible for any solver, and in the compilation of the SOS problem there are some numerical steps (finding a basis for Q for instance) where this kind of horrible numerics will destroy things.

[Yalmiper]

Hi Johan,

I remember you pointed out similiar thing once, although I am not sure they are the same issue.

However, you are calculating the max cofficient here, so I normalize everything with the minimum coefficient, here is 1047200, this time is works fine. Does this sounds a reasonable approach to you? Basically I need to try to scale my variable so that the coefficients should roughly in the same magnitude, is that correct?

Thank you!

[Johan Löfberg]

YALMIP scales internally anyways, but you essentially have the same bad numerics, as you now have very small elements instead. No matter what you do once you have these coefficients, you will have a small/large or huge spread in coefficients in the model, which is a recipe for numerical problems

[Johan Löfberg]

Yes, you have to go back to the polynomial and use a better representation, and not just work with the given coefficients

[Yalmiper]

Thank you! I also realized normalization does not change anything. Working on the variable scaling now.