Rounding issue in Takagi when complex

[nquesada]

When using decomposition.takagi there is a problem when the matrix is complex and has singular value close to zero. The funcion takagi has an optional argument called rounding which uses to decide how many decimal digits it should keep when deciding if a two singular values are the same, the default value is rounding=13. Note in the following example how this fails:

from strawberryfields.utils import random_interferometer
from strawberryfields.decompositions import takagi

real = False # Whether to use real or complex symmetric matrices
rounding = 13 # 13 is the default value in Takagi
np.random.seed(137)
n = 10 # Half the number of nonzero singular values
num_zeros = 3 # Number of zero singular values
abs_vals = 10 * np.random.rand(n)
diag_vals = np.concatenate([abs_vals, -abs_vals, np.zeros([num_zeros])])
U = random_interferometer(2 * n + num_zeros, real=real)
A = U @ np.diag(diag_vals) @ U.T
r, V = takagi(A, rounding = rounding)
test1 = np.allclose(V @ V.T.conj(), np.identity(2*n + num_zeros))
test2 = np.allclose(V @ np.diag(r) @ V.T, A)
print(test1, test2)

which gives False False. Note that if rounding=12 of real=True everything works as expected.

[josh146]

@nquesada, what do you think the best solution to this is? If we change rounding=12, are we just delaying this bug (e.g., is it possible another edge case will fail later?)

[nquesada]

It is certainly tempting to just change rounding, but I'd prefer a more satisfactory solution. Will try to find sometime this week to think about this issue.

[nquesada]

I did some research into this issue. There are at least two other ways to implement Takagi-Autonne.
The first outlined in wikipedia https://en.wikipedia.org/wiki/Symmetric_matrix#Complex_symmetric_matrices relies on jointly diagonalizing two normal matrices that commute. As it turns out the QuTiP people have faced this problem before and found more or less the same issues that we found already, i.e., that you run into trouble with degenerancies (cf. scipy/scipy#8445).

The second approach follows what is done here https://pypi.org/project/takagi-fact/ but unless you are working with mpmath it will fail when dealing with matrices that are degenerate near zero, like the one constructed above. I suspect that, surprisingly, our current implementation is the best we will be able to get.

[josh146]

Nice! So I guess the solution is stick with the existing Takagi decomposition?

We could do something similar to what happens in SciPy/Matlab etc. when you ask for the matrix square root of a matrix that is near-singular; a warning is raised that says 'Matrix is ill-conditioned, results might be inaccurate'

[mariaschuld]

@nquesada can we close this issue or is this on the todo list?

[nquesada]

Not sure. It is still a bug, but as far as I can tell there is no fix to it.