GEVP treated with bisection / penlab / bmibnb

[mar sag]

Dear Johan, 

I'm trying to implement the following GEVP by using both the bisection algorithm and alternatively penlab, as you suggest in the decay rate example. The GEVP I want to solve is the following 

given a full-rank matrix M, determine P>0, Q>0, t >1, such that 

t is minimized, and 

t*Q>=M'*P*M>=Q, t>1.

I tried with the bisection and alternatively I used penlab as well. The solvers state that a solution is found, but when I check the numerical values of P and Q, they are not positive definite (they are not positive semi-definite either). 

Am I coding anything wrong?

Enclosed you can find the two minimal working examples. I would be glad if you could tell me what´s going on, as I rechecked the code several times, and it seems consistent to me. In any case for a quick look I'm copying the main functions here too. (the subfunction only hardcodes the input)

Thanks in advance, 

Best Regards,

Marco

%%%%%% EXAMPLE 1 - use of bisection 

function GEVP_testing

% solving GEVP problem s.t.

% P>0, Q>0, (t)*Q>=M'*P*M>=Q, t>1

% load the matrix

M = load_M;

% check the rank

if (rank(M)-size(M,1))<0

    disp('matrix is not full-ranked');

    return

end

% cond number

c0 = cond(M);

% define the dimensions of the matrices P and Q

dim = size(M);

p = dim(1);

q = dim(2);

% define in Yalmip the matrices P and Q, and the variable t; they must be PSD matrices

t = sdpvar(1,1);

P = sdpvar(p,p,'diagonal');

Q = sdpvar(q,q,'diagonal');

% define matrices to ensure P and Q are PD

eps_p = 1e-15*eye(p);

eps_q = 1e-15*eye(q);

% define the LMI problem

Objective = t;

% the LMI is defined here - eps_p and eps_q ensure the PD of P and Q. since

% t has to be larger than 1, but don't explicitly bound it, we include a

% +1, and t_true = t+1;

F = [P-eps_p>=0, Q-eps_q>=0, Q<=M'*P*M, (t+1)*Q>=M'*P*M];

% set the solver and use the bisection command

h_opt = bisection(F, Objective,sdpsettings('solver','sedumi'))

% check the positive definiteness

% check the positive definiteness

positivedefinite = prod((eig(value(P))>0)) &&  prod((eig(value(Q))>0));

if positivedefinite==0

    disp('not positive definite')

end

% end of test

%%%%% EXAMPLE 2 - use of penlab

 

function GEVP_testing_alternative

% solving GEVP problem s.t.

% P>0, Q>0, (t)*Q>=M'*P*M>=Q, t>1

% load the matrix

M = load_M;

% check the rank

if (rank(M)-size(M,1))<0

    disp('matrix is not full-ranked');

    return

end

% cond number

c0 = cond(M);

% define the dimensions of the matrices P and Q

dim = size(M);

p = dim(1);

q = dim(2);

% define in Yalmip the matrices P and Q, and the variable t; they must be PSD matrices

t = sdpvar(1,1);

P = sdpvar(p,p,'diagonal');

Q = sdpvar(q,q,'diagonal');

% define matrices to ensure P and Q are PD

eps_p = 1e-15*eye(p);

eps_q = 1e-15*eye(q);

% version without eps_p and eps_q - not working

% F = [P>=0, Q>=0, Q<=M'*P*M, t*Q>=M'*P*M, c0 >= t >= 2];

F = [P>=eps_p, Q>=eps_q, Q<=M'*P*M, (t)*Q>=M'*P*M, c0 >= t >= 2];

Objective = t;

% version with bmibnb - not working

% sol = optimize(F,t,sdpsettings('solver','bmibnb'));

sol = optimize(F,t,sdpsettings('solver','penlab'));

% check the positive definiteness

positivedefinite = prod((eig(value(P))>0)) &&  prod((eig(value(Q))>0));

if positivedefinite==0

    disp('not positive definite')

end

end

[Mark L. Stone]

I think this is a case of solver tolerances.  On your examples, the min eigenvalue (equivalently, diagonal value) of P was -2e-14, and positive for Q using bisection with sedumi. Min eigenvalues for P and Q of -1e-7 and -2e-8 respectively using penlab. Your eps values of 1e-15 are insufficient to guarantee actual positive definiteness given solver tolerances.

[Johan Löfberg]

bisection maximizes the scalar which is specified as objective. Nope, not consistent with other commands...Note that you cannot simply switch the sign of objective, as the code is very limited and assumes that the objective is a simple scalar with no coefficients or anything. Hence, you have to move the sign convention to the constraints

Additionally, you must mean P - eps_p*eye(p) >= 0 etc, and as Mark says using a margin of 1e-15 makes no difference as it is lost in the general numerical tolerances of the solver

[Johan Löfberg]

Disregard my comment about the missing identity matrix. Saw that you had that already

[mar sag]

Thank you Johan, and thanks to Mark as well. 

Is there any way to provide initial guesses when I use bisection or penlab for my problem? Since the unitary matrices Q and P correspond to the solution t = c0, this would mean that I could always provide a valid initial guess, and every improvement from this condition would be welcome. 

Moreover, about the tolerances you were mentioning, can I change it by using sdpsettings?

Thanks!

[Johan Löfberg]

With bisecton, initial values are not relevant (as the problems are solved using a standard sdp solver, and yalmip automatically derives lower/upper bounds on the scalar which is bisected). For penlab, use standard assign/usex0 logic

You'll have to look at solver specific settings, but you can typically not improve precision to any significant degree. The solvers do as good as they can. 

[Mark L. Stone]

In order to improve precision significantly, presuming you can get it to work (under YALMIP) and have the computational resources (time), perhaps VSDP as the LMI solver called by bisecrtion with tight tol?

[mar sag]

thanks Mark. Yalmip does not support VSDP anymore

http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Solvers.VSDP

In any case I need to solve the problem on a normal machine in a CPU time, which is in the order of some seconds. This is what happens at the moment with the  penlab / bisection schemes, so I think I´ll try to refine those schemes on the basis of your suggestions. Once it is working I'll update the post in case it can be useful for other users as well,

Thanks a lot!

Marco

[Johan Löfberg]

With reasonable numerics in your data, and reasonable size on the problem, this should be a trivial problem. Your main flaw at the moment though is that the problem isn't feasible.

First note that the model is homogenuous, i.e., if P and Q are feasible, then so are x*P and x*Q. Hence, you might just as well say Q>=I

Now try to solve the following relaxation of your problem

optimize([P>=0, Q>=eye(q), Q<=M'*P*M])

It's infeasible when you use sedumi. Sedumi has problems with this model. With Mosek, it works (of course, as M'*M is psd, so just pick P=s*I for s large enough)

However, now you want to add t*Q >= M'*P*M. We can directly compute an upper bound on a feasible t to be 

K>> max(eig(value(M'*P*M)))/min(eig(value(Q)))

ans =

     1.323372942778177e+07

I.e., very ill-conditioned. bisection gives up if the upper bound isn't found when testing t = 1e6 (to be precise, it tests 2^19 which fails, and then gives up as 2^20 is larger than 1e6)

If you edit bisection.m and change 1e6 to 1e8, and use mosek, you will see that bisection finds the optimal value around 9.85e5

F = [P>=1e-7*eye(p), Q>=eye(q), Q<=M'*P*M, (-t+1)*Q>=M'*P*M];

h_opt = bisection(F, Objective,sdpsettings('solver','mosek'))