Takagi tol

[RyosukeNORO]

Context: There were issues with Takagi decomposition, specifically sqrtm, with some matrix.

Description of the Change: Implemented the tolerance parameter to kill small values in the elements of the matrix to be Takagi-decomposed.   Benefits: Takagi decomposition works well.

Possible Drawbacks: Takagi decomposition does not work well in other matrices because I have not fully understood why the decomposition have not done well before the change.

Related GitHub Issues:

[nquesada]

Could you provide the matrix for which takagi fails @RyosukeNORO ?

[RyosukeNORO]

Could you provide the matrix for which takagi fails @RyosukeNORO ?

Here you are. matrix_inaccurate_Takagi.zip

The zip file includes .npy file of a matrix for which takagi fails. You can simply load the .npy file A = np.load('matrix_inaccurate_Takagi.npy') and run takagi(A). When I printed the matrix and copied it to create a new matrix, takagi with the new matrix works well.

[nquesada]

Thanks. That is a very interesting edge case you found. I suggest you use tol to check if your matrix is diagonal:

np.allclose(A, np.diag(np.diag(A)), rtol, atol) is True

If this the case, then the Takagi is trivial:

d = np.diag(A)
U = np.diag(np.exp(1j*0.5*np.angle(d)))
l = np.abs(d)

[nquesada]

The issue seems to be related to sqrtm having problems with diagonal matrices :| .

You might want to check if an issue exists for this on the scipy repo and if not make one.

[elib20]

Thanks Nico! I think it would be fine if we can just hardcode the diagonal edge case rather than suppress these numerical zeros. There are already multiple different cases that are treated differently in the function. What do you think?

[nquesada]

Agreed! I don't like the idea of suppressing things. One thing that requires some thinking is whether we need both an rtol and atol. I think you don't need rtol if you are doing np.allclose(something, 0) : https://numpy.org/doc/stable/reference/generated/numpy.allclose.html `

I do suggest to call the function input rtol following the usage from np.allclose. Also, please add this matrix as a test into into test_decompositions. That way we know we are covering this edge case and the tests bot won't yell at you.

[nquesada]

Also, please make sure to order the diagonal elements based on their absolute value (in increasing or decreasing order, following svd_order)

[nquesada]

And also, extra also, why don't you make a PR directly into main?