Another day working with MATLAB, another nasty surprise. I was solving a SDP that was already working well, did a minor change, and suddenly it started taking minutes to reach a solution, when usually it took seconds.

After a long investigation, I realized that the problem was that my input was not exactly Hermitian anymore. I had switched from writing manually the components of a matrix to defining the matrix as $A = \ket{\psi}\bra{\psi}$ for a complex vector $\ket{\psi}$. Now how could that possibly be a problem? It is a very simple theorem to prove that the outer product of a vector with itself is always Hermitian. Well, not with MATLAB. I can’t figure out why, but probably MATLAB uses some smart algorithm for the outer product that ruins this property due to floating-point arithmetic. What I could determine is that if you compute the outer product the naïve way, then even with MATLAB the result will be exactly Hermitian. Also, with Python the output is exactly Hermitian.

To solve this problem, one can simply redefine $A$ as $(A+A^\dagger)/2$, which is very fast and gives an exactly Hermitian output. I don’t like doing that, though, as $(A+A^\dagger)/2$ is always exactly Hermitian, even when $A$ is not even close to Hermitian. If I fucked up the definition of $A$ this will hide my mistake, and make debugging harder. Instead, what I did was write my own outer product function, computing it as $\ket{\psi}\otimes\bra{\phi}$. It is slower than whatever black magic MATLAB is doing, sure, but it is fast enough, and leaves my mistakes easily visible.

It dawned on me that this subtlety was probably the source of many bugs and slowdowns in my codes over the years. I decided to go hunting, and find out the fascinating properties of MATLAB algebra that doesn’t preserve Hermiticity where it should. It turns out that if $A$ and $B$ are exactly Hermitian, then both $A+B$ and $A \otimes B$ are exactly Hermitian, as they should. The problem is really when you do matrix multiplication.

Which shouldn’t produce further problems, right? After all, $AB$ is in general not Hermitian, so we don’t have anything to worry about. Except, of course, that the Hilbert-Schmidt inner product $\operatorname{tr}(A^\dagger B)$ is real for Hermitian $A,B$, and this is not true in MATLAB. Argh! $\operatorname{tr}(A^\dagger B)$ appears very often in the objective of a SDP, and it really needs to be a real number, as you can’t minimize a complex function. It turns out that this is not MATLAB’s fault, it is a fundamental problem of floating-point arithmetic.

An well-known but often-forgotten fact is the floating-point addition is not associative. Let $a=10^{20}$ and $b=10^{-20}$. Then with floating point numbers we see that $b+(a-a) = b$ but $(b+a)-a = 0$. Which is the issue with the Hilbert-Schmidt inner product: we get stuff like $z+\bar{z}+w+\bar{w}+\ldots$, but not in order. $z+\bar{z}$ is of course real, and $w+\bar{w}$ as well, but $z+(\bar{z}+w)+\bar{w}$? Nope, not in floating-point arithmetic.

Here I don’t think a proper solution is possible. One can of course write $\Re[\operatorname{tr}(A^\dagger B)]$, but that will hide your mistakes when $A$ and $B$ are not in fact Hermitian. A better idea is to vectorize $A$ and $B$ to avoid matrix multiplication, writing $\Re[\langle A|B\rangle]$. This only helps with speed, though, it doesn’t touch the real problem.