nonlinear convex optimization

[Dylan]

Hi All,

I have a problem which is convex, but nonlinear, and Yalmip returns "no suitable solver" (I am running yalmip('version')=3).  Code looks something like:

p=sdpvar(4,4,'hermitian','complex');

j=0;

len=length(counts);

for k=1:len

    j=j-counts(k)*log(trace(o(:,:,k)*p));

end

ops = sdpsettings('sedumi.eps',1e-16,'sedumi.numtol',1e-16,'sedumi.maxiter',100000,'sedumi.bigeps',1e-16,'verbose',2);

solvesdp([p>0,trace(p)==1],j)

The variable o is a list of 4X4 hermitian matrices, and the variable counts is a list of integers>0.  The key to the problem is in the log function here.  It was my understanding that fmincon should just take over?  I tried forcing yalmip to use it explicitly but still got the same error.

Thanks for any help you can offer!

-Dylan

[Johan Löfberg]

You have to use the logdet operator, and this problem probably requires a solver which has native logdet support (sdpt3)

[Dylan]

Hi Johan,

Thanks for the reply, but I don't want to minimize the log determinant, but rather the log trace, which I do not think are equivalent.  

-Dylan

[Dylan]

Also, logdet can "only be applied to hermitian SDVAR objects" (I'm quoting Matlab here), and the product of o and p is not necessarily hermitian.

Really at a loss as to how to solve this problem.  I know for a fact this function is well behaved and convex.  For example I could parametrize p in a positive hermitian form and then just use fminsearch, but this would be crazy slow and as I am interested in solving this problem many times this is not the greatest option.

[Johan Löfberg]

Aha, didn't see the nonsymmetric product with o.

Then the problem is not structurally convex, and there are no direct solvers available.

[Johan Löfberg]

I take that back, I was confused by your claim that trace destroyed the setup. The problem is convex and standard logdet problem (assuming o*p is real). log(tr()) is directly solvable by log(det()) since the logarithm acts on the scalar trace result, i.e., log(tr(X)) = log(det(tr(X))

X = sdpvar(3);
O1 = randn(3);
O2 = randn(3);
solvesdp([X >= 0, trace(X)==1],-logdet(trace(O1*X))-logdet(trace(O2*X)))

[Johan Löfberg]

BTW, fmincon would not take over since it doesn't support SDP constraints

[Dylan]

Ah, thank you very much!  

-Dylan

[Dylan]

Hmm, now I'm getting some funny results.  If I solve the problem (Yalmip automatically selects sedumi for the task):

p=sdpvar(4,4,'hermitian','complex');

j=0;

len=length(counts);

for k=1:len

    j=j+counts(k)*logdet(trace(o(:,:,k)*p));

end

solvesdp([p>0,trace(p)==1],-j)

and output j with double(j), and then I simply calculate what j is from the solution (using double(p) in place of p and log in place of logdet) I get two different answers!  They are significantly different (for one specific example I got about a factor of 4 difference).  Also, if I explicitly force sdpt3 to solve the problem, I get NaN for all the elements of double(p).

-Dylan

[Johan Löfberg]

Things like this are best discussed with a fully reproducible model (create code I can run which behaves in the same way)

Note, when using sedumi to solve this, you are not solving the logdet problem but a problem involving geomean operators. This can lead to very different solutions when you have several terms.http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Commands.Logdet. Non-logdet capable solvers are only called for historical reasons, it should really warn about this...

[Dylan]

To make the reproducible code simpler, I scaled the problem down a bit, here is a function: 

function rx = MaximumLikelihoodR(o,counts)

p=sdpvar(2,2,'hermitian','complex');

j=0;

len=length(counts);

for k=1:len

    j=j+counts(k)*logdet(trace(o(:,:,k)*p));

end

ops = sdpsettings('solver','sdpt3','verbose',2);

solvesdp([p>0,trace(p)==1],-j,ops)

rx=double(p);

which you can run with the following command:

MaximumLikelihoodR(cat(3,[1 0; 0 0],[0 0; 0 1],1/2*[1 1;1 1],1/2*[1 i; -i 1]),[100 1 60 40])

which yields NaN.  

The answer should be very close to [1 0; 0 0].

Any thoughts?

Thanks again so much.

-Dylan

[Johan Löfberg]

Does that code really run? Fails here, since trace(o(:,:,4)*p) is complex

K>> sdisplay(trace(o(:,:,4)*p))

0.5000*p(1,1)+0.5000*p(2,2)+(0-3.0000e+03i)*p(2,1)

[Dylan]

Hmmm, not sure how you get that, tr(o(:,:,4)*p) = 0.5*p(1,1)+0.5*p(2,2)+0.5*i*(p(2,1)-p(1,2)).

Since p is explicitly hermitian, that last bit is 0.5*i*2*i*Imag(P(2,1)) which is real.

-Dylan

[Johan Löfberg]

I had something messed up in my workspace.

Solves nicely with both sdpt3 (and sedumi but then you are solving another problem). The optimal solution is


   0.8711             0.2700 + 0.1985i
   0.2700 - 0.1985i   0.1289         

Is the problem really solved when you run SDPT3, i.e, is the iteration log visible. If not, turn on the option debug in sdpsettings

However, you did catch a bug in the display of logdet terms.

>> sdpvar x
>> assign(x,2)
>> 2*logdet(x)+logdet(1+x)
Logdet-term 1x2 (current value: 3.5835)
>> 2*log(det(x))+log(det(1+x))
Linear scalar (real, derived, 2 variables, current value : 2.4849)

Edit logdet/display and change X.gain(1) to X.gain(i)

[Dylan]

Hmmm... I tried turning debug on and received the following error:

??? Cell contents indices must be greater

than 0

Error in ==> callsdpt34 at 42

            options.sdpt3.parbarrier{end-i+1,1}

            = -K.maxdetgain(i);

Error in ==> solvesdp at 291

    eval(['output = ' solver.call

    '(interfacedata);']);

Error in ==> MaximumLikelihoodR at 20

solvesdp([p>0,trace(p)==1],-j,ops)

-Dylan

[Johan Löfberg]

There is no such code on line 42 in the latest YALMIP version. Are you really running the latest YALMIP version?

[Dylan]

Sorry about that!  I have now installed the latest version of Yalmip, and the program in fact runs, but it yields a solution which is *visibly* false. 
It is still giving me solutions which I do not expect, however.  The solution you give in a previous quote is the same one I get, but its not the one I expect from the problem.  I would have expected something much closer to [1 0 ;0 0].

Where is logdet/display located?  I can't find it anywhere.

Thanks
-Dylan

[Dylan]

Nevermind, I just found my bug, which is not related to Yalmip. 

Thank you so much for such a great tool, and for all your help!
-Dylan