fix merge conflict

.github/CHANGELOG.md

@@ -2,17 +2,30 @@
 
 ### New features
 
+* Adds the Takagi decomposition [(#363)](https://github.com/XanaduAI/thewalrus/pull/338)
+
 ### Breaking changes
 
 ### Improvements
 
+
+* Tighten power-trace bound of odd loop Hafnian. [(#362)](https://github.com/XanaduAI/thewalrus/pull/362)
+
+* Simplifies the internal working of Bloch-Messiah decomposition [(#363)](https://github.com/XanaduAI/thewalrus/pull/338). 
+
+* Simplifies the internal working of Williamson decomposition [(#366)](https://github.com/XanaduAI/thewalrus/pull/338). 
+
+* Improves the handling of an edge case in Takagi [(#373)](https://github.com/XanaduAI/thewalrus/pull/373).
+
 ### Bug fixes
 
 ### Documentation
 
 ### Contributors
 
-This release contains contributions from (in alphabetical order):
+This release contains contributions from (in alphabetical order): 
+
+Gregory Morse, Nicolas Quesada
 
 ---
 

.readthedocs.yml

@@ -9,19 +9,18 @@ version: 2
 sphinx:
   configuration: docs/conf.py
 
-# Build documentation with MkDocs
-#mkdocs:
-#  configuration: mkdocs.yml
-
 # Optionally build your docs in additional formats such as PDF and ePub
 formats:
   - htmlzip
 
-# Optionally set the version of Python and requirements required to build your docs
+# Set the version of Python and container image
+build:
+  os: ubuntu-22.04
+  tools:
+    python: "3.10"
+
+# Set the requirements required to build your docs
 python:
-  version: 3.8
   install:
     - requirements: docs/requirements.txt
 
-build:
-  image: latest

docs/requirements.txt

@@ -1,13 +1,13 @@
-docutils==0.16
+docutils
 sphinxcontrib-bibtex
-sphinx==3.5.3
+sphinx
 ipykernel
-nbsphinx==0.7
+nbsphinx
 numba>=0.49.1
 numpy>=1.9
 sympy>=1.5.1
 scipy>=1.8.0
-Jinja2==2.11.3
+Jinja2
 version_information
 sphinx-copybutton
 dask[delayed]

thewalrus/_hafnian.py

@@ -264,7 +264,7 @@ def f_loop_odd(AX, AX_S, XD_S, D_S, n, oddloop, oddVX_S):  # pragma: no cover
     count = 0
     comb = np.zeros((2, n + 1), dtype=np.complex128)
     comb[0, 0] = 1
-    powtrace = charpoly.powertrace(AX, n + 1)
+    powtrace = charpoly.powertrace(AX, n // 2 + 1)
     for i in range(1, n + 1):
         if i == 1:
             factor = oddloop

thewalrus/decompositions.py

@@ -70,36 +70,31 @@ def williamson(V, rtol=1e-05, atol=1e-08):
     omega = sympmat(n)
     vals = np.linalg.eigvalsh(V)
 
-    for val in vals:
-        if val <= 0:
-            raise ValueError("Input matrix is not positive definite")
+    if not np.all(vals > 0):
+        raise ValueError("Input matrix is not positive definite")
 
     Mm12 = sqrtm(np.linalg.inv(V)).real
     r1 = Mm12 @ omega @ Mm12
     s1, K = schur(r1)
-    X = np.array([[0, 1], [1, 0]])
-    I = np.identity(2)
-    seq = []
-
-    # In what follows I construct a permutation matrix p  so that the Schur matrix has
+    # In what follows a permutation matrix perm1 is constructed so that the Schur matrix has
     # only positive elements above the diagonal
-    # Also the Schur matrix uses the x_1,p_1, ..., x_n,p_n  ordering thus I permute using perm
+    # Also the Schur matrix uses the x_1,p_1, ..., x_n,p_n  ordering thus a permutation perm2 is used
     # to go to the ordering x_1, ..., x_n, p_1, ... , p_n
-
+    perm1 = np.arange(2 * n)
     for i in range(n):
-        if s1[2 * i, 2 * i + 1] > 0:
-            seq.append(I)
-        else:
-            seq.append(X)
-    perm = np.array([2 * i for i in range(n)] + [2 * i + 1 for i in range(n)])
-    p = block_diag(*seq)
-    Kt = K @ p
-    Ktt = Kt[:, perm]
-    s1t = p @ s1 @ p
-    dd = [1 / s1t[2 * i, 2 * i + 1] for i in range(n)]
-    Db = np.diag(dd + dd)
-    S = Mm12 @ Ktt @ sqrtm(Db)
-    return Db, np.linalg.inv(S).T
+        if s1[2 * i, 2 * i + 1] <= 0:
+            (perm1[2 * i], perm1[2 * i + 1]) = (perm1[2 * i + 1], perm1[2 * i])
+
+    perm2 = np.array([perm1[2 * i] for i in range(n)] + [perm1[2 * i + 1] for i in range(n)])
+
+    Ktt = K[:, perm2]
+    s1t = s1[:, perm1][perm1]
+
+    dd = np.array([1 / s1t[2 * i, 2 * i + 1] for i in range(n)])
+    dd = np.concatenate([dd, dd])
+    ddsqrt = np.sqrt(dd)
+    S = Mm12 @ Ktt * ddsqrt
+    return np.diag(dd), np.linalg.inv(S).T
 
 
 def symplectic_eigenvals(cov):
@@ -112,8 +107,8 @@ def symplectic_eigenvals(cov):
         (array): symplectic eigenvalues
     """
     M = int(len(cov) / 2)
-    D, _ = williamson(cov)
-    return np.diag(D)[:M]
+    Omega = sympmat(M)
+    return np.real_if_close(-1j * np.linalg.eigvals(Omega @ cov))[::2]
 
 
 def blochmessiah(S):
@@ -216,7 +211,10 @@ def takagi(A, svd_order=True):
         sorted_ls, permutation = zip(*list_vals)
         return np.array(sorted_ls), Uc[:, np.array(permutation)]
 
-    phi = np.angle(A[0, 0])
+    # Find the element with the largest absolute value
+    pos = np.unravel_index(np.argmax(np.abs(A)), (n, n))
+    # Use it to find whether the input is a global phase times a real matrix
+    phi = np.angle(A[pos])
     Amr = np.real_if_close(np.exp(-1j * phi) * A)
     if np.isrealobj(Amr):
         vals, U = takagi(Amr, svd_order=svd_order)

thewalrus/symplectic.py

@@ -438,39 +438,6 @@ def is_symplectic(S, rtol=1e-05, atol=1e-08):
     return np.allclose(S.T @ Omega @ S, Omega, rtol=rtol, atol=atol)
 
 
-def autonne(A, rtol=1e-05, atol=1e-08, svd_order=True):
-    r"""Autonne-Takagi decomposition of a complex symmetric (not Hermitian!) matrix.
-
-    Args:
-        A (array): square, symmetric matrix
-        rtol (float): the relative tolerance parameter between ``A`` and ``A.T``
-        atol (float): the absolute tolerance parameter between ``A`` and ``A.T``
-        svd_order (boolean): whether to return result by ordering the singular values of ``A`` in descending (``True``) or asceding (``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")
-    if not np.allclose(A, A.T, rtol=rtol, atol=atol):
-        raise ValueError("The input matrix is not symmetric")
-    Areal = A.real
-    Aimag = A.imag
-
-    B = np.empty((2 * n, 2 * n))
-    B[:n, :n] = Areal
-    B[n : 2 * n, :n] = Aimag
-    B[:n, n : 2 * n] = Aimag
-    B[n : 2 * n, n : 2 * n] = -Areal
-    vals, vects = np.linalg.eigh(B)
-    U = vects[:n, n : 2 * n] + 1j * vects[n : 2 * n, n : 2 * n]
-    if svd_order:
-        return (vals[n : 2 * n])[::-1], U[:, ::-1]
-    return vals[n : 2 * n], U
-
-
 def xxpp_to_xpxp(S):
     """Permutes the entries of the input from xxpp ordering to xpxp ordering.
 

thewalrus/tests/test_decompositions.py

@@ -396,8 +396,6 @@ def test_real_degenerate():
     assert np.allclose(U @ np.diag(rl) @ U.T, mat)
 
 
-
-
 @pytest.mark.parametrize("n", [5, 10, 50])
 @pytest.mark.parametrize("datatype", [np.complex128, np.float64])
 @pytest.mark.parametrize("svd_order", [True, False])
@@ -415,3 +413,24 @@ def test_autonne_takagi(n, datatype, svd_order):
         assert np.all(np.diff(r) <= 0)
     else:
         assert np.all(np.diff(r) >= 0)
+
+@pytest.mark.parametrize("size", [10, 20, 100])
+def test_flat_phase(size):
+    """Test that the correct decomposition is obtained even if the first entry is 0"""
+    A = np.random.rand(size, size) + 1j * np.random.rand(size, size)
+    A += A.T
+    A[0, 0] = 0
+    l, u = takagi(A)
+    assert np.allclose(A, u * l @ u.T)
+
+
+def test_real_input_edge():
+    """Adapted from https://math.stackexchange.com/questions/4418925/why-does-this-algorithm-for-the-takagi-factorization-fail-here"""
+    rng = np.random.default_rng(0)  # Important for reproducibility
+    A = (rng.random((100, 100)) - 0.5) * 114
+    A = A * A.T  # make A symmetric
+    l, u = takagi(A)
+    # Now, reconstruct A, see
+    Ar = u * l @ u.T
+    assert np.allclose(A, Ar)
+