A crude hack is: trap the exception and check if the code returned (info) satisfies info>0. If so, add normally distributed random numbers with standard deviation 1e-20 to the matrix and try to compute the eigenvalues again. That hack works in this case.

It could be that it's not just that the matrix is almost the identity. The degeneracy in the eigenvalues, or the diagonal, is necessary to cause the error... in this case at least. Changing one of the diagonal entries so that it ends in 9 rather than 8 is enough to see success.

In this case, multiplying the matrix by $\sqrt{-1}$ also avoids the error.

A crude hack is: trap the exception and check if the code returned (info) satisfies info>0. If so, add normally distributed random numbers with standard deviation 1e-20 to the matrix and try to compute the eigenvalues again. That hack works in this case.

The changes in #6896 were intended to avoid the randomized algorithm in TwoQubitWeylDecomposition, which after a lot of battle-hardening seems finally 🤞 immune to such pathologies. I guess rather than trapping the error and ad-hoc randomizing, I'd rather trap the error and fall back to the TwoQubitWeylDecomposition alogorithm. Although having two paths like this worries me, and if the general eigensolver is unreliable maybe its better to simply use the TwoQubitWeylDecomposition algorithm everywhere unless the performance hit is substantial.

I think the use of randomization is quite different. And pretty small in terms of code and complexity. But, I understand that introducing another piece of randomness might result in a series of tweaks as failing cases trickle in under general use.

If there are understandable characteristics of offending matrices, we could make some kind of test suite in order to be more confident (if not perfectly) that a solution will work. I suppose a small subset of these would go in the big test suite.

EDIT: There already is something like this

The terra testsuite exercises a bunch of the pathological cases. It would be interesting to see if the testsuite passes on your setup. Something like python -m unittest test.python.quantum_info.test_synthesis -cb should run the relevant pieces

Apropos environment-specific features: The example above succeeds when using MKL (on my machine). I'm not 100% sure because I used Julia (I didnt look up yet how to use mkl with scipy.linalg). But, it looks like python and julia are calling the same lapack routines. Furthermore, with openblas, Julia fails with the same error as numpy.

It's definitely a function of the environment beyond just the specific blas implementation. On my linux system I do not encounter any failures with openblas, but on my m1 mac also using openblas it does reliably fail.

Just to follow up here, we've moved most of the two qubit synthesis code to rust and the eigensolver in faer seems seems to be immune from this issue. I have generally observed that running on arm64 mac numpy/openblas has a lot of these numeric issues, besides this specific error we've seen a lot of weird numeric stability issues and failures on arm64 (I suspect it might be an apple clang issue, but lack the time and knowledge to really dig into it) and moving the numerics to rust has been a pretty blunt tool to avoid them. That being said @jlapeyre wrote the weyl_coordinates function in rust already, but we're not leveraging it from python yet, although the only remaining use of it is for tests. I'll push up a small PR to remove the scipy implementation and just expose the rust function to python and have that close this issue.

It looks like we could move the stuff in local_invariance.py to Rust. Then test this from Rust rather than Python. But exposing weyl_coordinates from Python is reasonable in the meantime.