search

cvxgrp / scs

Splitting Conic Solver
MIT License
544 stars 134 forks source link

Support for complex PSD cone #290

Open araujoms opened 2 weeks ago

araujoms commented 2 weeks ago

Are you interested in adding support for the complex PSD cone? It shows up all the time in quantum information problems, and is much more efficient than mapping the problem to the real PSD cone via

A \mapsto \begin{pmatrix}\Re(A) & \Im(A)\\ -\Im(A) & \Re(A)\end{pmatrix}

I just did it for another ADMM solver, COSMO: oxfordcontrol/COSMO.jl#194

If you're interested I volunteer to do it for SCS as well. I would need help with the interface, though.

bodono commented 2 weeks ago

Sure, that would be very interesting! I'm happy to help with the interface. If you get something working we can run some tests in python to see the advantage the approach has over converting into a regular PSD cone, I guess it should be something like 4x faster?

araujoms commented 2 weeks ago

I'm glad you're interested! Naïvely I would expect a 8x speedup, as the complexity in ADMM is dominated by diagonalization, which is roughly $d^3$ complexity, and we're dividing dimension by two. This is roughly what I observed when benchmarking COSMO.

My proposal is that I take care of core functions like set_up_sd_cone_work_space and proj_semi_definite_cone, and you connect that to the rest of SCS so that they can actually be used. In Julia I used some metaprogramming to avoid code duplication, but for C that's a dubious proposition, so I'd just accept the inelegance.

A design choice that needs to made now is how to vectorize the complex matrix. I'd use the same vectorization as Hypatia, adapted to the lower triangular for consistency with your real psd cone. That means, x[0,0], real(x[1,0]), imag(x[1,0]), x[1,1], real(x[2,0]), imag(x[2,0]), etc. Is that fine?

Also, is there any particular reason why you're using syev instead of syevr? The latter is significantly faster.

bodono commented 1 week ago

Great, thanks for the explanation! That vectorization seems fine, though do you need to scale the diagonals by sqrt(2)? See details here for why it's necessary for the usual sdp cone: https://www.cvxgrp.org/scs/api/cones.html#semidefinite-cones

I switched from syevr to syev a while back, I think for robustness issues, though perhaps syevr is better: https://github.com/cvxgrp/scs/commit/e07d5abb255c798d174dabc268eba122a718cb20#diff-3ff078d764f5668dad6750cc5905768ecdfb80e76f67d3e8fa2f6e66fd1c2bc9R562

araujoms commented 1 week ago

Great, I'll use this vectorization then, together with syevr. I know one needs to scale the off-diagonals by sqrt(2), I just didn't feel it was worth mentioning.