[tomdodd4598]

Context:

1. The current computation of the W matrix in takagi does not always produce a unitary matrix (up to floating point noise). The cause seems to be the existence of eigenvalues of the matrix to be square-rooted very close to the branch cut of the complex square root through an angle of ±π. We avoid this source of inaccuracy by temporarily adding an additional phase to the matrix which best clears the branch cut of any eigenvalues, computing the square root, and then removing half of the same phase.

2. (Not important as discussed below) The current implementation uses (vh @ np.conjugate(u)).T in the matrix square root, when the correct expression to use seems to be (u.T @ vh).conj() = vh @ u.conj() as shown in box 1 on page 5 of the paper referenced in the docstring.

An example of the current implementation failing can be found here.

Description of the Change:
Improve the current general-case Takagi decomposition algorithm.

Benefits:
Makes takagi correctly return unitary matrix.

[nquesada]

What does this fix?

[tomdodd4598]

Currently there are certain matrices for which the takagi function returns a W matrix which is not accurately unitary. For example, this test prints False.

Should we expect that (Uᵀ V)* is equal to (V U*)ᵀ?

[nquesada]

Yes. They should be the same. The issue is likely a subtle numerical problem

[tomdodd4598]

If so, then perhaps, but the numerical instability does seem rather bad. For that particular example, the absolute value of WW† is
[[1.00212044 0.01949431 0.00495123 0.0017475 ]
[0.01949431 0.99789623 0.0030558 0.00418326]
[0.00495123 0.0030558 1.00077574 0.00122061]
[0.0017475 0.00418326 0.00122061 0.99964825]].

With the new implementation, it is
[[1.00000000e+00 6.15535818e-15 7.44900262e-16 1.43805873e-15]
[6.15535818e-15 1.00000000e+00 1.74707629e-15 6.43570492e-16]
[7.44900262e-16 1.74707629e-15 1.00000000e+00 1.39291723e-15]
[1.43805873e-15 6.43570492e-16 1.39291723e-15 1.00000000e+00]].

[heltluke]

Could you add your test here @tomdodd4598?

[tomdodd4598]

I've moved this PR back to a draft while I work out if/when this causes the function to return a less accurate decomposition for certain inputs. I'll also add a suitable test once I better understand what's happening.

[nquesada]

Currently there are certain matrices for which the takagi function returns a W matrix which is not accurately unitary. For example, this test prints False.

Should we expect that (Uᵀ V)* is equal to (V U*)ᵀ?

This might help to understand:

from https://cdnsciencepub.com/doi/epdf/10.1139/cjp-2024-0070

[heltluke]

So what's going wrong with this example is that V^T U is not diagonal? Is there a threshold below which the singular values of \Lambda should be considered "machine zero"? complex256 will allow for some very small numbers...

[heltluke]

I think it's important to remember that numpy.linalg.svd returns the hermitian conjugate of V (as in Eq. (34) of the paper).

Suggested change

	u, d, v = np.linalg.svd(A)
	u, d, vh = np.linalg.svd(A)

Thus, I'm not convinced that either of these expressions matches the expression in Box 1. There the desired expression is (U^T V)* = (V^T U)* -> vh @ u.conj(). What was previously written on line 221 has an extra transpose (.T), and the new expression on line 221 has not taken into account the relationship between V and vh.

[nquesada]

I think it's important to remember that numpy.linalg.svd returns the hermitian conjugate of V (as in Eq. (34) of the paper).

Thus, I'm not convinced that either of these expressions matches the expression in Box 1. There the desired expression is (U^T V)* = (V^T U)* -> vh @ u.conj(). What was previously written on line 221 has an extra transpose (.T), and the new expression on line 221 has not taken into account the relationship between V and vh.

what you call the last matrix in the SVD is immaterial, it is just that, the last unitary:

[heltluke]

Hmm... if all of these products are equivalent, but sqrtm doesn't seem to treat them as such, then indeed something funny seems to be going on.

[tomdodd4598]

@nquesada As well as the question Luke has above (perhaps there is a reason I'm not seeing for why it's not an issue), I think I found a way around the problem with the matrix square root. The issue is that it can sometimes be non-unitary with a relatively large error (well beyond machine epsilon).

A solution that appears to work in testing is to compute the SVD of this square root, yielding two other unitary matrices and another diagonal matrix which holds all of the non-unitary error. We can then multiply the two new unitaries together to get the "corrected" unitary square root.

[nquesada]

@nquesada As well as the question Luke has above (perhaps there is a reason I'm not seeing for why it's not an issue), I think I found a way around the problem with the matrix square root. The issue is that it can sometimes be non-unitary with a relatively large error (well beyond machine epsilon).

A solution that appears to work in testing is to compute the SVD of this square root, yielding two other unitary matrices and another diagonal matrix which holds all of the non-unitary error. We can then multiply the two new unitaries together to get the "corrected" unitary square root.

This works in that it forces a unitary to be returned. However, it will not guarantee that the decomposition is correct, since if the singular values of the second SVD are not ones, then something went truly amiss. On the other hand the solution to this problem escapes me: it seems to be a subtle issue from numerical precision.

@tomdodd4598 do add your name here: https://github.com/XanaduAI/thewalrus/blob/master/.github/ACKNOWLEDGMENTS.md

[tomdodd4598]

This works in that it forces a unitary to be returned. However, it will not guarantee that the decomposition is correct, since if the singular values of the second SVD are not ones, then something went truly amiss.

@nquesada You are correct - unfortunately I initially missed this, as the specific situation I was working with only threw an error for non-unitary U, but the more direct test of takagi indeed continued to throw assertion errors.

However, I've just made another commit which, in testing, seems to fix the error completely, even for pathological inputs. The idea is the following: once we have a unitary U, the associated D matrix (equal to U†AU*) may be non-diagonal, leading to an inaccuracy. We try to fix this by iteratively putting D through more Takagi decompositions and rotating U in response.

Non-pathological matrices will not go through any of these iterations. In our tests, it seems like close-to-real matrices with imaginary component noise of order ~1e-10 go through just one.

[nquesada]

@tomdodd4598 this sound very promosing. As @heltluke advised, it would be great to have a strong battery of tests, specially for rank deficient and also for degenerate symmetric matrices, and also for matrices with a small imaginary parts on the background.

[nquesada]

and of course: thanks for catching this insidious bug and working on squashing it!

[tomdodd4598]

@arsalan-motamedi Has found a much nicer solution which doesn't require a loop to refine the accuracy, but instead uses some random phases to move any bad matrix elements away from a complex argument close to ±π. I have tried thousands of repeats of tests across all of the Takagi tests already in thewalrus (rank deficient, degenerate, etc.), as well as real random matrices with small imaginary noise (just committed here), and not seen a failure:
=== 327000 passed, 51000 deselected in 1067.55s (0:17:47) ===

[nquesada]

Would there be a way to make that phase a function of the input matrix and not a random number?

[tomdodd4598]

@nquesada Done, and all tests passed again. Instead of a random angle, the best angle is computed by looking at the complex arguments of the eigenvalues of the matrix to be square-rooted, and finding the angle of the midpoint of the largest arc not containing any of them. The matrix entries are multiplied by the negative of this phase so that the problematic region near ±π is clear of points. After the safe square root, the result is then rotated back by half of the angle.

[nquesada]

Was reading the docs of numpy.linalg.eig . They explicitly mention that if you want diagonalize normal matrices (the ones that commute with their dagger, which includes hermitian as well as anti-hermitian and unitary you are much better off diagonalizing them using the scipy.linalg.schur decomposition:

Not sure how much more time you want to spend on this PR, but this might be something worth considering, especially if you have good reasons to believe that what is supposed to be unitary is indeed unitary.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 100.00%. Comparing base (12633ca) to head (1285ade).
⚠️ Report is 4 commits behind head on master.

Additional details and impacted files

@@            Coverage Diff            @@
##            master      #403   +/-   ##
=========================================
  Coverage   100.00%   100.00%           
=========================================
  Files           28        28           
  Lines         1922      1936   +14     
=========================================
+ Hits          1922      1936   +14     

Files with missing lines	Coverage Δ
thewalrus/decompositions.py	100.00% <100.00%> (ø)

---

Continue to review full report in Codecov by Sentry.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update deb45fd...1285ade. Read the comment docs.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[ziofil]

I think it would be best to avoid sqrtm altogether: you could take the sqrt of the eigenvalues instead no?

[nquesada]

I think it would be best to avoid sqrtm altogether: you could take the sqrt of the eigenvalues instead no?

See my comment above schur vs eig for unitaries.

[arsalan-motamedi]

I think it would be best to avoid sqrtm altogether: you could take the sqrt of the eigenvalues instead no?

See my comment above schur vs eig for unitaries.

I agree with using a better matrix square root method. But the problem here was with the usual square root. You can see that conventions are so that at around -1, the square root funciton is discontinuous (and in general around any point on -x axis):

Our solution here was to rotate the eigen-values away from there (by finding some optimal angle), and then rotate back after.

[ziofil]

sure sure, i agree. my suggestion was about efficiency