Hessian?

[Sébastien Loisel]

Hi,

I'm working on some large-scale optimization problems. Long story short: does YALMIP calculate the Hessian of the Lagrangian in a smart way or does it resort to numerical differentiation or something?

I've got a custom interior point method for my problem so I have a good idea of what should happen. I've also solved my problem using CVX and SeDuMi and SDPT3; I'm just testing a wider set of solvers using YALMIP. YALMIP accepts that my model is convex with 'allownonconvex',0 in the options. Nevertheless, with 66 unknowns (tiny for this class of problems), it takes IPOPT 896 iterations but what baffles me even more is that it says that it evaluated the Lagrangian of the Hessian 0 times. I read the IPOPT paper and it's based on a backtracking line search in the Newton search direction so I'm not sure what is happening. Only thing I can think of is: does it fall back on some sort of gradient descent or BFGS thing, which would explain something maybe?

So does YALMIP provide the Hessian of the Lagrangian?

Thanks,

S

[Johan Löfberg]

YALMIP only supplies gradients

ipopt is very sensitive to conditioning of objective if it doesn't have any hessian. Other solvers such as fmincon often performs better in this case. Alternatively, clever coordinate changes or moving stuff from objective to constraints can change things significantly

A = randn(150);

S = diag(1:150)

Q = A'*S'*S*A

x = sdpvar(150,1);
% Ill-conditioned hessian

optimize(-1 <= x <= 1, x'*Q*x+sum(S*A*x),sdpsettings('solver','ipopt'))

e = sdpvar(150,1);
% Nice hessian

optimize([e == S*A*x,-1 <= x <= 1], e'*e+sum(S*A*x),sdpsettings('solver','ipopt'))

[Johan Löfberg]

and if you have a problem that you can solve using sdpt3/sedumi, using ipopt sounds weird. It means the problem is an LP/QP/SOCP (meaing cplex,gurobi,mosek would be most interesting) or an sdp (meaning mosek would be most interesting, but you also have solvers like csdp,sdpa,sdpnal etc)

[Sébastien Loisel]

> YALMIP only supplies gradients

Well that answers the original question. However, switching the solver to 'fmincon' results in bad juju that I don't understand. I'll reply to your other things after I show you the bad juju code:

yalmip('clear');

% Set some options for YALMIP and solver

options = sdpsettings('verbose',1,'solver','fmincon','allownonconvex',0);

% Define variables

Dx = randn(5,3); Dy = randn(5,3);

cx = randn(5,1); cy = randn(5,1);

f = randn(3,1);

p = 1.1;

u = sdpvar(3,1);

y = [Dx*u+cx, Dy*u+cy];

z = sdpvar(5,1);

C = cone([z';y']);

O = dot(ones(5,1),z.^p)/p+dot(f,u);

% Solve the problem

sol = optimize(C,O,options);

This results in "Error using barrier Too many input arguments." Any clue what this means?

Responses to your other things.

1. This is a research problem which is "badly conditioned by design". It's a weighted average of p-th powers of Euclidian norms, plus a linear term ``forcing''. I can't find any way of phrasing it so the Hessian is good.

2. CVX is able to reduce p-th powers to quadratic cone programming but this is some nontrivial reduction described in a paper I haven't read. My custom barrier method is primal only because I can prove things about the primal method. I'd like to use a "canned" primal-dual barrier method that someone else has already tested a lot. My understanding is we don't have a good analysis of long-step primal-dual methods when the cone is not self-dual so this is just for numerical comparison.

3. I have not tried the commercial offerings but it would indeed make sense to enable the Gurobi and Mosek solvers in CVX to see what happens. SDPT3 chokes on my problem but SeDuMi ends up doing okay after this nontrivial reduction of pth powers to some sort of quadratic cone program.

[Johan Löfberg]

which versions are you using? ipopt runs without issues here and easily solves the problem

If you want to us an socp representaton and socp solver, you use cpower(z,p) instead of z.^p

[Johan Löfberg]

...and fmincon works too