An issue with the penlab solver

[Eike F]

Greetings,

I have the following problem and would appreciate your advice:

I wish to find the greatest ball around 0 on which a polynomial P is positive. Of course in 1 dimension this could be done by finding its roots but I'm interested in the general case.

Therefore I wrote this code:

x=sdpvar(1);

P=2*x^2-1./3.*x^4;

%polynomial describing the ball
radius_squared=sdpvar(1);
ball=(x'*x-radius_squared);

%multiplier for Positivstellensatz
L=polynomial(x,6);

constraints=[sos(L),sos(P+L*ball)];

ops=sdpsettings('sos.model',2);

%The tiny quadratic part in the objectice is a "convexifier".
[diag,mons,grams,resids]=solvesos(constraints,-radius_squared+0.01*radius_squared^2,ops,[radius_squared;coefficients(L,x)])

The only bmi solver I have is PENLAB, which returns:

PenLab didn't converge: unconstrained minimization failed

or falsely that radius_squared is about 15^2.

If, on the other hand I choose the fixed value radius_squared = 2^2, CSDP correctly states feasibility.
Lastly, if I choose radius_squared=3^2, CSDP correctly states infeasibility.


How can I get PENLAB to solve my problem?

Best,

Eike

[Johan Löfberg]

I would not use penlab (for the exact reasons you just experienced) for this simple quasi-convex problem. A simple bisection solves the problem easily

P = optimizer(constraints,[],sdpsettings('solver','sedumi','sos.model',2),radius_squared,coefficients(L,x))

lower = 0;
upper = 100;
while upper-lower > 1e-3
    test = (lower + upper)/2;
    [~,fail] = P{test};
    if fail
        upper = test;
    else
        lower = test
    end
end

[Eike F]

A pity. Is this due to theoretical problems which apply to all bmi solvers or is it just an issue with penlab?

Thanks for the code!

[Mark L. Stone]

From limited experience, PENBMI, which is the non-free version, seems to converge more reliably than PENLAB (the open source version from the same authors). I had several nonlinear SDPs for which PENLAB had a propensity to diverge (even playing with algorithm parameter settings), but PENBMI reliably converged to a valid local optimum, but none for which the reverse occurred.

[Eike F]

Thanks for the hint! Could I put you to the trouble of trying my code with penbmi?

[Johan Löfberg]

Both. BMIs are intractable so all solvers will have problems on some problems. Of course, they might behave better or worse.

However, with your simple structure, I would use bisection as above with a linear solver. You then leverage on the solid numerical performance of linear sdp solvers, and you know you get the globally optimal solution.