A couple of years ago I was complaining about the lack of direct support for SDPs with complex numbers, and documenting my efforts to get it working on SeDuMi via YALMIP. That worked, but in the meantime I’ve stopped using MATLAB (and thus YALMIP), and started solving SDPs in Julia via JuMP.
No problem, I did once, I can do it again. Funnily enough, the work was pretty much the complement of what I did for YALMIP: YALMIP already had a complex interface for SeDuMi, it was just disabled. JuMP didn’t have one, I had to write it from scratch. On the other hand, I had to write the dualization support for YALMIP, but not JuMP, because there we get it for free once the interface is working. The codebases couldn’t be more different: not only Julia is a much better programming language than MATLAB, but also YALMIP is a one-man project which has been growing by accretion for decades1. In contrast, JuMP is a much newer project, and has been written by many hands. I has a sane, modular design by necessity.
That was done, but it was pointless, as SeDuMi is just not a competitive solver. To get a real benefit I’d need to add complex support to the best solvers available. My favourite one is Hypatia, which was in a similar situation: the only thing lacking was the JuMP interface. I contacted the devs, and they quickly wrote the interface themselves. Amazing! But that was also the end of the line: to the best of my knowledge no other solver could handle complex numbers natively.
I thought maybe I can hack the solver myself, as the algorithm is really simple. I looked a bit around, and found the perfect target: COSMO, an ADMM solver written in Julia. It was harder than I thought because COSMO fiddled directly with LAPACK, and that requires dealing with a FORTRAN interface from the 70s. Nevertheless I did it, and my PR was accepted. Luckily the interface came for free, as COSMO is written to work with JuMP directly. But it was also a bit pointless, because it turns out COSMO is not really good either, and it’s mostly abandoned.
I decided to go for the best ADMM solver available: SCS. The problem is, it’s written in C2. But I had a strong motivation to do it: I was working on a paper that needed to solve some huge complex SDPs that only SCS could handle. Getting a speedup there was really worth it. And it turned out to be much easier than I thought: the codebase is rather clean and well-organized. The only difficulty was dealing with LAPACK, but I had already learned how to do it from my work with COSMO. It just took ages for my PR to be accepted, so long that my paper was already done by then, and I had no intention of running the calculations again.
With the solver done, I went to work on the JuMP interface, which was just released today. As part of my evil plan to get people to switch to Julia, I didn’t write the YALMIP interface. I mean, writing MATLAB code for something that I’ll never use personally? No way. In any case, the new version of SCS showed a 4x speedup on an artificial problem designed to display the advantage, and a 75% speedup on a real complex SDP I had lying around.
Is that the end of my crusade? I really wanted to add support to a second-order solver, as these are faster and more precise (at the cost of not being able to handle the largest problems). The best one available is MOSEK. But it’s a closed source solver owned by a private, for-profit corporation, so nope. I actually tried to contribute a bug report to them once, and as a result I got personally insulted by MOSEK’s CEO.
An interesting possibility is Clarabel. It is written in Julia, so it should be really easy to hack. And in principle it should be much faster than Hypatia, as it uses the specialized algorithm for symmetric cones, whereas Hypatia uses the generic one that can also handle non-symmetric cones. But my benchmarks showed that it is not competitive with Hypatia, so until that changes there is no point, and I did enough pointless work.
Thank you for your efforts to improve the status of SDP solvers for quantum information problems! I have a couple of quick questions.
1) Do I understand correctly that the current state of the art (in terms of speed) for large quantum information SDPs is the native complex solver you have written for SCS? And for the moment, the only interface available to use this solver (without manually deriving the canonical form) is the one you have written for JuMP.jl, right? Is anybody working on a CVXPY interface perhaps?
2) Regarding MOSEK: if I understand correctly, it does not directly support SDPs with complex numbers, but only through the standard conversion into real numbers (doubling the dimension), right? So something like what you have done for SCS to natively support complex numbers would need to be implemented by the company, I guess.
1) Indeed, to the best of my knowledge there is no better solver for large-scale complex SDPs (for small-scale ones MOSEK is better). And the only interface for it is the one I wrote. I don’t have any contact with the CVXPY devs, sorry, it’s better to ask yourself. It’s not looking promising, though, because they support a solver that has a native complex PSD cone, QICS, but don’t seem to actually use it, but do the inefficient real mapping instead.
2) Yes, exactly.
Hi, a quick update.
I opened a feature request on cvxpy’s Github repo, and one of the developers replied that he likes the idea of supporting the hermitian PSD cone https://github.com/cvxpy/cvxpy/issues/3045 (I mentioned you, this blog post and your comment, I hope it’s fine).
Let’s see what happens. Thanks again for your work on this.
Thanks for doing that, I hope they add it. Maybe with increasing interest we’ll get proper complex support in Clarabel and MOSEK as well.