Update decompositions.py and readds the Takagi decomposition

thewalrus/decompositions.py

@@ -125,9 +125,9 @@ def blochmessiah(S):
         S (array[float]): 2N x 2N real symplectic matrix
 
     Returns:
-        tupple(array[float],  : orthogonal symplectic matrix uff
-               array[float],  : diagional matrix dff
-               array[float])  : orthogonal symplectic matrix vff
+        tuple(array[float],  : orthogonal symplectic matrix uff
+              array[float],  : diagonal matrix dff
+              array[float])  : orthogonal symplectic matrix vff
     """
 
     N, _ = S.shape
@@ -148,16 +148,11 @@ def blochmessiah(S):
     alpha = Unitary[0 : N // 2, 0 : N // 2]
     beta = Sig[0 : N // 2, N // 2 : N]
     # Bloch-Messiah in this Basis
-    u2, d2, v2 = np.linalg.svd(beta)
+    d2, takagibeta = takagi(beta)
     sval = np.arcsinh(d2)
-    takagibeta = u2 @ sqrtm(np.conjugate(u2).T @ (v2.T))
-    uf = np.block([[takagibeta, 0 * takagibeta], [0 * takagibeta, np.conjugate(takagibeta)]])
-    vf = np.block(
-        [
-            [np.conjugate(takagibeta).T @ alpha, 0 * takagibeta],
-            [0 * takagibeta, np.conjugate(np.conjugate(takagibeta).T @ alpha)],
-        ]
-    )
+    uf = block_diag(takagibeta, takagibeta.conj())
+    blc = np.conjugate(takagibeta).T @ alpha
+    vf = block_diag(blc, blc.conj())
     df = np.block(
         [
             [np.diag(np.cosh(sval)), np.diag(np.sinh(sval))],
@@ -172,3 +167,67 @@ def blochmessiah(S):
     vff = np.real_if_close(vff)
     uff = np.real_if_close(uff)
     return uff, dff, vff
+
+
+def takagi(A, svd_order=True):
+    r"""Autonne-Takagi decomposition of a complex symmetric (not Hermitian!) matrix.
+    Note that the input matrix is internally symmetrized. If the input matrix is indeed symmetric this leaves it unchanged.
+    See `Carl Caves note. <http://info.phys.unm.edu/~caves/courses/qinfo-s17/lectures/polarsingularAutonne.pdf>`_
+
+    Args:
+        A (array): square, symmetric matrix
+        svd_order (boolean): whether to return result by ordering the singular values of ``A`` in descending (``True``) or ascending (``False``) order.
+
+    Returns:
+        tuple[array, array]: (r, U), where r are the singular values,
+        and U is the Autonne-Takagi unitary, such that :math:`A = U \diag(r) U^T`.
+    """
+
+    n, m = A.shape
+    if n != m:
+        raise ValueError("The input matrix is not square")
+    # Here we force symmetrize the matrix
+    A = 0.5 * (A + A.T)
+
+    A = np.real_if_close(A)
+
+    if np.allclose(A, 0):
+        return np.zeros(n), np.eye(n)
+
+    if np.isrealobj(A):
+        # If the matrix A is real one can be more clever and use its eigendecomposition
+        ls, U = np.linalg.eigh(A)
+        U = U / np.exp(1j * np.angle(U)[0])
+        vals = np.abs(ls)  # These are the Takagi eigenvalues
+        phases = -np.ones(vals.shape[0], dtype=np.complex128)
+        for j, l in enumerate(ls):
+            if np.allclose(l, 0) or l > 0:
+                phases[j] = 1
+        phases = np.sqrt(phases)
+        Uc = U @ np.diag(phases)  # One needs to readjust the phases
+        signs = np.sign(Uc.real)[0]
+        for k, s in enumerate(signs):
+            if np.allclose(s, 0):
+                signs[k] = 1
+        Uc = np.real_if_close(Uc / signs)
+        list_vals = [(vals[i], i) for i in range(len(vals))]
+        # And also rearrange the unitary and values so that they are decreasingly ordered
+        list_vals.sort(reverse=svd_order)
+        sorted_ls, permutation = zip(*list_vals)
+        return np.array(sorted_ls), Uc[:, np.array(permutation)]
+
+    phi = np.angle(A[0, 0])
+    Amr = np.real_if_close(np.exp(-1j * phi) * A)
+    if np.isrealobj(Amr):
+        vals, U = takagi(Amr, svd_order=svd_order)
+        return vals, U * np.exp(1j * phi / 2)
+
+    u, d, v = np.linalg.svd(A)
+    U = u @ sqrtm((v @ np.conjugate(u)).T)
+    # The line above could be simplifed to the line below if the product v @ np.conjugate(u) is diagonal
+    # Which it should be according to Caves http://info.phys.unm.edu/~caves/courses/qinfo-s17/lectures/polarsingularAutonne.pdf
+    # U = u * np.sqrt(0j + np.diag(v @ np.conjugate(u)))
+    # This however breaks test_degenerate
+    if svd_order is False:
+        return d[::-1], U[:, ::-1]
+    return d, U

thewalrus/tests/test_decompositions.py

@@ -1,4 +1,4 @@
-# Copyright 2019 Xanadu Quantum Technologies Inc.
+# Copyright 2023 Xanadu Quantum Technologies Inc.
 
 # Licensed under the Apache License, Version 2.0 (the "License");
 # you may not use this file except in compliance with the License.
@@ -19,7 +19,7 @@
 
 from thewalrus.random import random_interferometer as haar_measure
 from thewalrus.random import random_symplectic
-from thewalrus.decompositions import williamson, blochmessiah
+from thewalrus.decompositions import williamson, blochmessiah, takagi
 from thewalrus.symplectic import sympmat as omega
 from thewalrus.quantum.gaussian_checks import is_symplectic
 
@@ -253,3 +253,144 @@ def test_transform(self, passive, create_transform, tol):
         assert np.allclose(O1.T @ O @ O1, O, atol=tol, rtol=0)
         assert np.allclose(O2.T @ O @ O2, O, atol=tol, rtol=0)
         assert np.allclose(S @ O @ S.T, O, atol=tol, rtol=0)
+
+
+@pytest.mark.parametrize("n", [5, 10, 50])
+@pytest.mark.parametrize("datatype", [np.complex128, np.float64])
+@pytest.mark.parametrize("svd_order", [True, False])
+def test_takagi(n, datatype, svd_order):
+    """Checks the correctness of the Takagi decomposition function"""
+    if datatype is np.complex128:
+        A = np.random.rand(n, n) + 1j * np.random.rand(n, n)
+    if datatype is np.float64:
+        A = np.random.rand(n, n)
+    A += A.T
+    r, U = takagi(A, svd_order=svd_order)
+    assert np.allclose(A, U @ np.diag(r) @ U.T)
+    assert np.all(r >= 0)
+    if svd_order is True:
+        assert np.all(np.diff(r) <= 0)
+    else:
+        assert np.all(np.diff(r) >= 0)
+
+
+@pytest.mark.parametrize("n", [5, 10, 50])
+@pytest.mark.parametrize("datatype", [np.complex128, np.float64])
+@pytest.mark.parametrize("svd_order", [True, False])
+@pytest.mark.parametrize("half_rank", [0, 1])
+@pytest.mark.parametrize("phase", [0, 1])
+def test_degenerate(n, datatype, svd_order, half_rank, phase):
+    """Tests Takagi produces the correct result for very degenerate cases"""
+    nhalf = n // 2
+    diags = [half_rank * np.random.rand()] * nhalf + [np.random.rand()] * (n - nhalf)
+    if datatype is np.complex128:
+        U = haar_measure(n)
+    if datatype is np.float64:
+        U = np.exp(1j * phase) * haar_measure(n, real=True)
+    A = U @ np.diag(diags) @ U.T
+    r, U = takagi(A, svd_order=svd_order)
+    assert np.allclose(A, U @ np.diag(r) @ U.T)
+    assert np.allclose(U @ U.T.conj(), np.eye(n))
+    assert np.all(r >= 0)
+    if svd_order is True:
+        assert np.all(np.diff(r) <= 0)
+    else:
+        assert np.all(np.diff(r) >= 0)
+
+
+def test_zeros():
+    """Verify that the Takagi decomposition returns a zero vector and identity matrix when
+    input a matrix of zeros"""
+    dim = 4
+    a = np.zeros((dim, dim))
+    rl, U = takagi(a)
+    assert np.allclose(rl, np.zeros(dim))
+    assert np.allclose(U, np.eye(dim))
+
+
+def test_takagi_error():
+    """Tests the value errors of Takagi"""
+    n = 10
+    m = 11
+    A = np.random.rand(n, m)
+    with pytest.raises(ValueError, match="The input matrix is not square"):
+        takagi(A)
+
+
+def test_real_degenerate():
+    """Verify that the Takagi decomposition returns a matrix that is unitary and results in a
+    correct decomposition when input a real but highly degenerate matrix. This test uses the
+    adjacency matrix of a balanced tree graph."""
+
+    vals = [
+        1,
+        2,
+        31,
+        34,
+        35,
+        62,
+        67,
+        68,
+        94,
+        100,
+        101,
+        125,
+        133,
+        134,
+        157,
+        166,
+        167,
+        188,
+        199,
+        200,
+        220,
+        232,
+        233,
+        251,
+        265,
+        266,
+        283,
+        298,
+        299,
+        314,
+        331,
+        332,
+        346,
+        364,
+        365,
+        377,
+        397,
+        398,
+        409,
+        430,
+        431,
+        440,
+        463,
+        464,
+        472,
+        503,
+        535,
+        566,
+        598,
+        629,
+        661,
+        692,
+        724,
+        755,
+        787,
+        818,
+        850,
+        881,
+        913,
+        944,
+    ]
+    mat = np.zeros([31 * 31])
+    mat[vals] = 1
+    mat = mat.reshape(31, 31)
+    # The lines above are equivalent to:
+    # import networkx as nx
+    # g = nx.balanced_tree(2, 4)
+    # a = nx.to_numpy_array(g)
+    rl, U = takagi(mat)
+    assert np.allclose(U @ U.conj().T, np.eye(len(mat)))
+    assert np.allclose(U @ np.diag(rl) @ U.T, mat)