Quadratic objective with SDP constraints

[ericphanson]

Is it possible to use Alpine for such a problem? It seems like JuMP 18.5 does not accept such problems when I try to formulate one for use with Alpine:

ERROR: JuMP does not support quadratic objectives for conic problems

[kaarthiksundar]

Currently Alpine.jl does not support SDP constraints

[ccoffrin]

If your problem is convex give SCS, Mosek or Pajarito a try.

[ericphanson]

Ah ok, thanks for letting me know. Unfortunately SCS and Mosek don't support quadratic objectives, and Pajarito yields the same problem with "ERROR: JuMP does not support quadratic objectives for conic problems" since it seems to be a limitation of MathProgBase (though maybe I can try to fool it with a dummy variable). Thanks again!

[ccoffrin]

Oh, I didn't realize. You can always add lifted variables to turn the objective into a linear one.

[harshangrjn]

I don't think JuMP supports quadratic and conic constraints too. So lifting may not help.

[harshangrjn]

@ericphanson But a work around is that you can write a (simple) quadratic constraint in the SDP form. For example: x^2 <= y*z can be written as an SDP constraint as follows: [x y; z x] \succeq 0

[ccoffrin]

In JuMP v0.18 you can use norm in the conic api. An example in PowerModels, https://github.com/lanl-ansi/PowerModels.jl/blob/master/src/core/objective.jl#L178

[ericphanson]

Thanks for the suggestions :). I'm not sure whether or not these things can work in my case. I don't think the problem is convex; for example, I tried COSMO.jl on JuMP 19, and it transforms my problem into the COSMO.jl form shown here: https://oxfordcontrol.github.io/COSMO.jl/stable/. After that transformation, the matrix P I get out is not PSD, which I think reflects a non-convexity in the original problem (and COSMO.jl gives a very poor solution, likely because my problem thus doesn't fit their assumptions). It it quadratic in the sense that two variables are multiplied, but not in the form e.g. x^T Q x.

Instead, my objective function is a sum of terms like

sum( tr( A[i] * X ) * Phi[i,j] * c[j] for i=1:n for j = 1:n )

where the A[i] are constant (PSD) matrices, X is a PSD variable (with some linear constraints), c is a constant vector, and Phi is a permutation matrix (or in the relaxed version, a doubly stochastic one). So it's basically an inner product of a vector whose components are traces against my PSD variable X with a permuted vector c, optimizing over both the permutation and the PSD variable.

So JuMP calls it a quadratic expression, since you have components of X multiplied by components of Phi. But I'm not sure if it can be cast as the norm of some vector.

One thing I can do is fix Phi and relatively quickly solve the resulting SDP (with e.g. Mosek, or CSDP which seems even faster for my problem), and then try to do some heuristic optimization over Phi, e.g. simulated annealing as an outer optimization. This gets me something (but not as good as certain other heuristic things I've tried), and it seems like that at least the relaxed version where Phi is doubly stochastic should be not so hard to solve in a better way.

[harshangrjn]

@ericphanson Are the elements of your Phi matrix binary variables since you mention that it is a permutation matrix? If they are indeed binaries, based on your objective function description, you should be able to get rid of all your nonlinearities by replacing every Phi_ij * X_ij with a set of auxiliary variables, say \hat{phi}_ij, and the following constraints:

\hat{phi}_ij >= Phi_ij * X_ij_L
\hat{phi}_ij >= X_ij + Phi_ij * X_ij_U - X_ij_U
\hat{phi}_ij <= Phi_ij * X_ij_U
\hat{phi}_ij <= X_ij + Phi_ij * X_ij_L - X_ij_L

where X_ij \in [X_ij_L, X_ij_U]

Note that the above linearization is exact, i.e., you should be able to use any MISDP solver, like Pajarito.jl, and obtain global optimum values for your problem. Of course, the convergence of solvers will now depend on the bounds of X_ij variables. If the bounds are tight, you should be able to solve it fairly fast.

If your Phi_ij-s are in the continuous set [0,1] instead of {0,1}, the above linearization, solved using Mosek, guarantees you a valid lower bound but not necessarily global optimal values.

Hope this helps!

P.S: If you are interested in a more formal description, we do this above linearization for a similar nonlinear MISDP problem in this paper: Eqns (22)-(24) in https://arxiv.org/pdf/1903.08354.pdf

[ericphanson]

Wow, that's great! Thank you, I will try this (and of course cite the paper if indeed this makes it into our work). Yes, they are binary ideally-- I figured I might have to relax them to be in [0,1] since I thought that might be easier than an integer program, but maybe not then! Based on the earlier comments I already tried replacing the products with auxiliary variables, but I didn't think of anything like this clever argument.

[harshangrjn]

👍 I will be now curious to see how well Pajarito.jl would do on your problem. :)

[ericphanson]

Pajarito solved my test problem exactly, in only 28 seconds! I am impressed. Well, with CSDP as the SDP solver it ran for a few minutes and then segfaulted, but with Mosek it worked perfectly (edit: reported here: jump-dev/CSDP.jl#39 (comment)). I'm now trying on a bigger version of the problem, but I might just come back in the morning and see how it did, since it will likely take awhile.

[ericphanson]

@harshangrjn just to say the paper is out now, and we used this linearization trick: https://arxiv.org/abs/2001.03598 (acknowledging you). Thanks very much for the help! The MISDP ends up being hard to solve for non-small problems but it’s already interesting to me for small problems; there’s an SDP formulation (with exponentially many variables) solvable for small problems as well, but the misdp formulation can prescribe a number of non-zero variables which leads to more informative solutions.

[harshangrjn]

@ericphanson Thats great! Thanks for keeping me posted. In general MISDPs are quite hard to solve, but I have been off-late working on developing relatively fast methods to obtain tight lower/upper bounds. We can take this offline.