sdpvar variable Matrix Inverse workaround in YALMIP?

[Mark L. Stone]

I would like to optimize Mahalanobis Distance squared, where the symmetric n by n matrix A being inverted is the optimization variable, subject to semidefinite constraint A>=0 and linear constraints on the entries in A.

Let x be an n by 1 vector of numerical; data for some specified value of n.  Then I would like to minimize x'*inv(A)*x,but :Function 'inv' is not defined for values of class 'sdpvar' ".I know that if I really wanted to minimize this, I cold get rid of the inverse by using the Schur complement and introducing a Linear Matrix Inequality converting it into an SDP..  And that would work for minimization.  However, I also want to maximize it (which would be a non-convex objective) and feed it into a global solver (e.g. BMIBNB) subject to the convex senidefinite and linear constraints. That Schur complement trick won't work for maximization when A >= 0 (although ti would if A <= 0, which it is not in my case).   Is there a way to do this in YALMIP?  Is there a nonlinear evaluation capability which can be added to do this, or some other workaround/approach?  Specifying 'full' (with my intention to add other constraints explicitly to force symmetric semidefinite) still results in the same error message,

Thanks.

[Johan Löfberg]

If you are willing to create a nonlinear nonconvex program, and thus accept that you most likely won't be able to solve the problem for anyhting but trivial instances, you can simply introduce a new matrix B, minimize x'Bx and add the bilinear equality A*B = eye(n).

[Johan Löfberg]

BTW, a more compact representation might be to minimize x^Tz and add the bilinear constraint x == Az. A bit too early in the morning to see if it is correct. I believe it is if you know A is definite, i.e., the inverse actually exist.

[Johan Löfberg]

Anyhow, internally in YALMIP, the third-order expression x'*B*x would be reduced to a series of bilinear/quadratic terms, so you probably benefit from going bilinear manually, since you might be able to do this more efficiently as you know the form, i.e., x^Ty, y==Bx, AB=I in case you work with the form where you explicitly ensure the inverse exist and create a model for it

[Mark L. Stone]

Thanks Johan, I should have thought of that workaround. Actually, for the maximization problem I believe there would be a minus in either the objective or the bilinear equality.  Your approach goes the opposite of the usual direction of moving non-convexities from constraints to the objective, as I had nice convex constraints in the "natural" formulation, but I guess there's no free lunch, so you have to pay the piper somewhere.

To be explicit, the matrices I'm interested in are 6 by 6.  In one version of the problem, I have only 9 free variables,with the remaining 12 unique elements constrained to known values, but in other variants, I have 9 fully free variables (other than semidefinite constraint) and only linear constraints (leaving some degrees of freedom) on the other 12 unique elements. I'm prepared that the run time could be quite long (hours to maybve even days)?

In case the above problem is not challenging enough, I believe your approach could also be applied in principle (leaving aside how well solvers could actually deal with it) to a related problem, in which, subject to the same constraints as before, the objective function is a linear function of the product of A and inv(A +C), where C is a known symmetric psd matrix, which could be formulated as a the objective function being a linear function of the quadratic A*B subject to (A + C)*B=eye(n) plus my other constraints on A  Again 6 dimensional.  Do you have any estimate as to the prospects of solving the Mahalanobis Distacnce Squared problem with 9 free variables, or the quadratic objective problem I just showed?

Thanks again for your prompt and very helpful response!!!

[Johan Löfberg]

Yes, if you want to maximize, you negate the objective.

I would not agree on the statement that it is standard to move nonconvexities from constraints to the objective. In deterministic global optimization, it is most common to untangle the model to simple canonical operators (univariate functions and bilinear terms), and then connect all auxiliary variables through equality constraints. The reason is that we have simple methods to derive cutting planes for these univariate functions and bilinear terms

9 variables does not sound out of reach. The proof is in the pudding. Try to hack up some simple model, and I can help you to see if it can be attacked using YALMIP and its interfaced solvers, and see if the model can be improved upon,

[Mark L. Stone]

I had only seen your first response at the time of my first response.  Based on what you have written, I shall endeavor to solve the Mahalanobis Distance Squared maximization via the formulation x'*y s.t. y==B*x,A*B==eye(n).  I am interpreting your comments to mean that doing it this way rather than without the intermediate variable y may result in an easier to solve problem for the solver.  Anyhow, it shouldn;t be hard to try both (other than the run time) and see which works better.  Is there a way of examining the model that YALMIP will pass to the solver to make a reasonable assessment as to which will be easier to solve?

For my quadratic objective problem (when bilinear constraint is introduced), it looks like I could reformulate the problem the same way as you suggested for the Mahalanobis Distance Squared problem.  Specifically, I have some linear function applied to C*B*F, where F is a known positive semidefinite matrix.  So I could make that C*Y, s.t. Y==B*F, (C+D)*B==eye(n).  I presume that this should benefit for a similar reason as the Mahalanobis Distance Squared problem.

Thanks again.

[Johan Löfberg]

Hmm, looks like YALMIP is more clever than I am

x = sdpvar(3,1);

B = sdpvar(3,3);

solvesdp([-1<=[x;B(:)]<=1,sum(x)==1],x'*B*x,sdpsettings('solver','bmibnb'))

y = sdpvar(3,1);

solvesdp([y == B*x, -1<=[x;B(:)]<=1,sum(x)==1],x'*y,sdpsettings('solver','bmibnb'))

The first model will have 27 variables of which 12 are nonlinear, and YALMIP will branch on 6 variables, while the second model will have 24 variables, out of which 12 are nonlinear and YALMIP will branch on all of these

[Johan Löfberg]

No, there is no way to see the intermediate model at the moment (unless one puts a break in the code etc)

[Mark L. Stone]

Perhaps there is a misunderstanding.  My x is an input vector (of numbers), not an sdpvar.

Also, do you have suggestions for the best free lower, upper, and LP (if applicable?) solvers to use on these problens?

[Johan Löfberg]

Ah, right. A bit easier then.

Depends if you are in academia or not, since cplex, gurobi, mosek, and xpress all are free for academia.

Otherwise, I am not particularly biased to anyone of the free solvers. scip has worked well recently when I have tried it (compiled variants available in the opti toolbox)

[Johan Löfberg]

Not sure I understand what you are saying.

If you define A as

Q=sdpvar(3,3,'full');A = [data Q;Q' data]

YALMIP not introduce complicating nonlinear expression if you then introduce B and A*B==I (you might get some trivial linear equalities compared to a manually crafted model, but that will not be anything important in the big picture)

[Mark L. Stone]

Sorry, my posts were getting a bit clobbered.  I was suggesting defining a 3 by 3 full sdpvar 'full' instead of a 6 by 6 sdpvar, but still making B  be a 6 by 6 sdpvar, and then using 3 by 3 block multiplication to implement the blinear equality.  So I would have explicitly eliminated 12 variables, which in the other formulation would be equality constrained ti numerical values.  If there's no advantage to the solver of doing it this way (maybe no difference), then I won't, since it just reduces model transparency and won;t work for other cases in which the 12 other elements are not known individually, but only by linear constraints.  I realize I have the sane nonlinearity both ways, so it's really just a matter of writing things out differently.

I am not in academia (well technically, I am an inactive student, but not with the active status I would need for free academic versions.

Thanks again for your help and suggestions.  I should probably let you get dome work done now.

[Mark L. Stone]

It's getting a bit late for me, so I'm not getting things quite clear the first tine.  I was suggesting exactly what you have interpreted, namely Q=sdpvar(3,3,'full');A = [data Q;Q' data], and of course A>=0, B==sdpvar(6,6).  Both Data blocks are symmetric psd, so A>=0 implies B>=0 in the block formulation.  Would it be a good idea to still specify B>=0 (I'm guessing yes)?  I don't really understand your answer as to whether this would be better than making A be a 6 by 6 sdpvar, A>=0, A(1:3,1:3)==upper block diagonal data, A(4:6,4:6)==lower block diagonal data.


[Johan Löfberg]

You mean B=sdpvar(6,6) I guess

If A is psd, and A*B==I then natuarally B>=0, if that is what you mean

It is a waste of variables to define A=sdpvar(6), and then immediately constraint large portions of the matrix to constant data using equalitiy constraints. Most often much better to define the model without these variables.

Remember, x=sdpvar(1,2) with the constraint x(2)==3 is an optimization problem in two variables with one equality constraint, whereas x = [sdpvar(1) 3] is an optimization problem with only 1 variable. I.e., there is a big difference between assignments and equality constraints

[Yulin]

Hello! 

I want to try a similar problem, which means there is a inverse of some sdp variable involving in the problem. 

Here is the very simple code:

clear all

x0=[ 2.7183;

   13.5914;

    1.4969;

   12.3700;

    1.2214;

    1.2214]

X=sdpvar(6,6,'full')

P=sdpvar(6,6)

P1=sdpvar(6,6) %P1=inv(P)

LL=sdpvar(6,6,'full') %LL=inv(P)*X

alpha=sdpvar(1)

G=[alpha,x0';x0,X'*LL]

ops = sdpsettings('verbose',1,'solver','bmibnb','debug',1);

solvesdp([P>0,G>0,P*P1==eye(6),P1*X==LL],[],ops)

However, it can not get a solution, and I don't  understand the debug part. Is it something wrong with my code or is it unsolvable itself? 

[Johan Löfberg]

I don't get any crash with debug messages. I do however get a diagnostic code that the search-space is unbounded. All variables involved nonlinearly must be explicitly bounded
http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Solvers.BMIBNBTheory

However, I would not waste to much energy on this. Spend your time on trying to get a better (convex) model instead.

BTW, I think you have something wrong in the model. It looks like you want the constraint G>0 to be interpreted as a semidefinite constraint. That will not be the case though, as X'*LL is a non-symmetric matrix. Hence, it is interpreted as a condition that all elements in G are non-negative.

And as I constantly repeat, strict inequalities has no relevance in YALMIP and associated solvers
http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Blog.Prepare-your-code