takagi fix

.github/CHANGELOG.md

@@ -15,6 +15,8 @@
 
 ### Bug fixes
 
+* Add the calculation method of `takagi` when the matrix is diagonal. 
+
 ### Documentation
 
 ### Contributors

thewalrus/decompositions.py

@@ -202,6 +202,15 @@ def takagi(A, svd_order=True):
         vals, U = takagi(Amr, svd_order=svd_order)
         return vals, U * np.exp(1j * phi / 2)
 
+    # If the matrix is diagonal, Takagi decomposition is easy
+    if np.allclose(A, np.diag(np.diag(A))):
+        d = np.diag(A)
+        U = np.diag(np.exp(1j * 0.5 * np.angle(d)))
+        l = np.abs(d)
+        if svd_order is False:
+            return l[::-1], U[:, ::-1]
+        return l, U
+
     u, d, v = np.linalg.svd(A)
     U = u @ sqrtm((v @ np.conjugate(u)).T)
     if svd_order is False:

thewalrus/tests/test_decompositions.py

@@ -323,6 +323,14 @@ def test_takagi_error():
     with pytest.raises(ValueError, match="The input matrix is not square"):
         takagi(A)
 
+def test_takagi_sepcific_matrix():
+    """Test the takagi decomposition works well for a specific matrix that was not deecomposed accuratelyin a previous version.
+    See more info in PR #393 (https://github.com/XanaduAI/thewalrus/pull/393)"""
+    A = np.load('test_matrix_for_takagi.npy')
+    d, U = takagi(A)
+    assert np.allclose(A, U @ np.diag(d) @ U.T)
+    assert np.allclose(U @ np.conjugate(U).T, np.eye(len(U)))
+    assert np.all(d >= 0)
 
 def test_real_degenerate():
     """Verify that the Takagi decomposition returns a matrix that is unitary and results in a