Since I’ve first used MATLAB I have dreamt of finding a replacement for it. Not only it is expensive, proprietary software, but also a terrible programming language. Don’t get me wrong, I’m sure it was amazing when it was invented, but this was in the 70s. We know better now. I’ve had to deal with so many fascinating bugs due to its poor design decisions!
Most recently, I had code that was failing because 'asdf' and "asdf" are very similar, but not exactly the same. The former is a character vector, and the latter is a string. Almost always you can use them interchangeably, but as it turns out, not always. Another insane design decision is that you don’t need to declare variables to work on them. I declared a matrix called constraints, worked on it a bit, and then made an assignment with a typo contraints(:,1) = v. Instead of throwing an error like any sane programming language, MATLAB just silently created a new variable contraints. Perhaps more seriously, MATLAB does not support namespaces. If you are using two packages that both define a function called square, you have to be careful about the order in which they appear in the MATLAB path to get the correct one. If you need both versions? You’re just out of luck.
Perhaps I should stop ranting at this point, but I just can’t. Another thing that drives me mad is that loop indices are always global, so you must be very careful about reusing index names. This interacts greatly with another “feature” of MATLAB, that i is both the imaginary unity and a valid variable name. If you have for i=1:3, i, end followed by a = 2+3*i you’re not getting a complex number, you’re getting 111. The parser is downright stone age, it can’t handle simple operators like +=, or double indexing like a(2)(4). To vectorize a matrix there’s no function, just the operator :, so if you want to vectorize a(2) you have to either call reshape(a(2),[],1), or define x = a(2) and then do x(:). Which of course leads to everyone and their dog defining a function vec() for convenience, which then all conflict with each other because of the lack of namespaces.
I wouldn’t be surprised to find out that the person that designed the function length() tortured little animals as a child. If you call it on a vector, it works as expected. But if you call it on an $m \times n$ matrix, what should it do? I think the most sensible option is to give the number of elements, $mn$, but it’s also defensible to give $m$ or $n$. MATLAB takes of course the fourth option, $\max(m,n)$. I could also mention the lack of support for types, Kafkaesque support for optional function arguments, mixing row and column vectors… It would keep me ranting forever. But enough about MATLAB. What are the alternatives?
The first one I looked at is Octave. It is open source, great, but its fundamental goal is to be compatible with MATLAB, so it cannot fix its uncountable design flaws. Furthermore, it isn’t 100% compatible with MATLAB2, so almost always when I have to use MATLAB because of a library, this library doesn’t work with Octave. If I give up on compatibility, then I can use the Octave extensions that make the programming language more tolerable. But it’s still a terrible programming language, and is even slower than MATLAB, so there isn’t much point.
Then came Python. No hope of compatibility here, but I accept that, no pain no gain. The language is a joy to program with, but I absolutely need some optimization libraries (that I’m not going to write myself). There are two available, CVXPY and PICOS. Back when I first looked at them, about a decade ago, neither supported complex numbers, so Python was immediately discarded. In the meanwhile they have both added support, so a few years ago I gave it a shot. It turns out both are unbearably slow. CVXPY gets an extra negative point for demanding its own version of partial trace and partial transposition3, but that’s beside the point, I can’t use them for any serious problem anyway. I did end up publishing a paper using Python code, but it was only because the optimization problem I was solving was so simple that performance wasn’t an issue.
After that I gave up for several years, resigned to my fate of programming with MATLAB until it drove me to suicide. But then Sébastien Designolle came to visit Vienna, and told me of a programming language called Julia that was even nicer to program than Python, almost as fast as C++, and had an optimization library supporting every solver under the Sun, JuMP. I couldn’t believe my ears. Had the promised land been there all along? After all, I knew Julia, it just had never occurred to me that I could do optimization with it.
I immediately asked Sébastien if it supported complex numbers, and if it needed funny business to accept partial transpositions. Yes, and no, respectively. Amazing! To my relief JuMP had just added support for complex numbers, so I hadn’t suffered all these years for nothing. I started testing the support for complex numbers, and it turned out to be rather buggy. However, the developers Benoît Legat and Oscar Dowson fixed the bugs as fast as I could report them, so now it’s rock solid. Dowson in particular seemed to never sleep, but as it turned out he just lives in New Zealand.
Since then I have been learning Julia and writing serious code with it, and I can confirm, the language is all that Sébastien promised and more. Another big advantage is the extensive package ecosystem, where apparently half of academia has been busy solving the problems I need. The packages can be easily installed from within Julia itself and have proper support for versions and dependencies4. Also worth mentioning is the powerful type system, that makes it easy to write functions that work differently for different types, and switch at runtime between floats and complex floats and double floats and quadruple floats and arbitrary precision floats. This makes it easy to do optimization with arbitrary precision, as JuMP in fact allows for the solvers that support it (as far as I know they are Hypatia, COSMO, and Clarabel). As you might know this is a nightmare in MATLAB.5
Now Julia is not perfect. It has some design flaws. Some are because it wants to be familiar to MATLAB users, such as having 1-indexed arrays and col-major orientation. Some are incomprehensible (why is Real not a subtype of Complex? Why is M[:,1] a copy instead of a view?). It’s not an ideal language, it’s merely the best that exists. Maybe in a couple of decades someone will release a 0-indexed version called Giulia and we’ll finally have flying cars and world peace.
It’s a bit ironic to write this blog post after I have released a paper that is based on a major MATLAB library that I wrote together with Andy, Moment6. In my defence, Andy wrote almost all of the code, the vast majority is in C++, and we started it before Sébastien’s visit. And it demonstrated beyond any doubt that MATLAB is completely unsuitable for any serious programming7. I promise that when I have time (haha) I’ll rewrite it in Julia.
But the time for irony is over. My new projects are all in Julia, and I’ll start releasing them very soon. In the meanwhile, I wrote a tutorial to help refugees from MATLAB to settle in the promised land.
Hi Matteus,
Enjoy your musings, as always.
Like you, I have had a similar trajectory of finding the right combination of tooling to invoke SDPs for problems in quantum information. I wouldn’t necessarily call myself a MATLAB refugee, but I used it extensively before Python solvers supported complex-valued matrices (which I’m glad they now do). These days, I’ve been building a lot of tooling around the solvers in Python in my toqito Python package.
Out of curiosity, have you tried solving your SDPs in Python using a solver like MOSEK? As you point out, the CVXOPT solver can be painfully slow but I’ve noticed that MOSEK is much better. I wonder how that compares to the SDP solvers in Julia–I’d be curious to know as that would be an interesting datapoint for me to seriously consider switching over (or at least to start using Julia more in my workflow).
I’ve gone back and forth on this for Julia. Indeed, when I started with `toqito`, the complex-valued thing was still an issue, and Julia was still quite cumbersome and challenging to work with (as you point out). Plus, many of the bugs that you requested the developers fix weren’t fixed at that time either.
Perhaps I should give the Julia thing another go. It seems we are interested in a fairly wide overlap of problems, so if you happen to have any Julia code that showcases some neat SDPs for quantum information, I’d love to play with it!
Cheers!
And I appreciate your comments, as always.
My complaint about Python’s performance was not about the solvers, but about Python itself: I did try MOSEK, for instance, and as expected it doesn’t make any difference, as MOSEK offloads the work to C anyway. But the pre-processing part is so painfully slow that it overshadows the solution time. Often it doesn’t matter, but I need to do a lot of non-commutative polynomial optimization, which needs rather nontrivial pre-processing. In fact it was my frustration with the Python package for it, ncpol2sdpa, that drove me to write Moment.
I did a couple of benchmarks with Julia, nothing serious, and the result was as expected: little difference in the solving times (except the difference from using primal versus dual forms), and much faster pre-processing.
As for Julia code, the really cool stuff is yet to be published, but I wrote an example for bounding random robustness at the end of the tutorial I linked above, and another for doing quantum state discrimination here.
To be fair to python, ncpol2sdpa could be much faster in many of the commonly used cases but it would require some overhaul of the package design. E.g., we encounter sets of operators {{A1, A2, …, Ar}, {B1,B2,…,Br},…} which do not necessarily commute within the subsets but two operators from different subsets will always commute. Ncpol2sdpa will create monomials, ignoring the commutation relations and then simplify at the global level — forgetting about the commutativity structure. However, for many applications, simplifications need only be made on the submonomials which are built from the noncommuting subsets and at the end you just multiply your submonomials to build your total monomial. This would likely be a huge speedup for ncpol2sdpa.
Thanks for letting me know about ncpol2sdpa, I had assumed that it used the specialized algorithm for Bell scenarios (which Moment does, btw).
But Python is still slow, nevertheless. I tried implementing my algorithm for the local bound in several languages, and Python was the slowest by far, orders of magnitude slower than C++. MATLAB was a bit better than Python but not much, and Julia had performance similar to C++.
It’s good to hear that Julia gains popularity it deserves!
I’ve been using Julia for quite a while now. It was not a love from first sight, Julia is harder that Python and has some uncommon design choices, after all. But eventually I realized how powerful it is. And it’s not only about speed. The killer feature is the multiple dispatch and its type system. It instantly makes Julia superior to other popular languages. I can write algorithms using abstract types, someone else can independently introduce specific types and it could work together flawlessly without much interaction between developers. E.g. think of combining the interval arithmetic and dual numbers. How much pain would it cost to combine the two in other languages, especially if developed independently?