Merge branch 'master' into master

.github/CHANGELOG.md

@@ -17,6 +17,10 @@
 
 ### Bug fixes
 
+* Add the calculation method of `takagi` when the matrix is diagonal. [(#394)](https://github.com/XanaduAI/thewalrus/pull/394)
+
+* Add the lines for avoiding the comparison of np.ndarray and list. [(#395)](https://github.com/XanaduAI/thewalrus/pull/395)
+
 ### Documentation
 
 ### Contributors

thewalrus/decompositions.py

@@ -152,7 +152,8 @@ def blochmessiah(S):
     return O, D, Q
 
 
-def takagi(A, svd_order=True):
+def takagi(A, svd_order=True, rtol=1e-16):
+    # pylint: disable=too-many-return-statements
     r"""Autonne-Takagi decomposition of a complex symmetric (not Hermitian!) matrix.
     Note that the input matrix is internally symmetrized by taking its upper triangular part.
     If the input matrix is indeed symmetric this leaves it unchanged.
@@ -162,6 +163,7 @@ def takagi(A, svd_order=True):
     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.
+        rtol (float): the relative tolerance parameter used in ``np.allclose`` when judging if the matrix is diagonal or not. Default to 1e-16.
 
     Returns:
         tuple[array, array]: (r, U), where r are the singular values,
@@ -202,6 +204,19 @@ 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)), rtol=rtol):
+        d = np.diag(A)
+        l = np.abs(d)
+        idx = np.argsort(l)
+        d = d[idx]
+        l = l[idx]
+        U = np.diag(np.exp(1j * 0.5 * np.angle(d)))
+        U = U[::-1, :]
+        if svd_order:
+            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/symplectic.py

@@ -171,6 +171,9 @@ def reduced_state(mu, cov, modes):
     """
     N = len(mu) // 2
 
+    if type(modes) == np.ndarray:
+        modes = modes.tolist()
+
     if modes == list(range(N)):
         # reduced state is full state
         return mu, cov

thewalrus/tests/test_decompositions.py

@@ -324,6 +324,35 @@ def test_takagi_error():
         takagi(A)
 
 
+@pytest.mark.parametrize("svd_order", [True, False])
+def test_takagi_diagonal_matrix(svd_order):
+    """Test the takagi decomposition works well for a specific matrix that was not decomposed accurately in a previous implementation.
+    See more info in PR #393 (https://github.com/XanaduAI/thewalrus/pull/393)"""
+    A = np.array(
+        [
+            [
+                -8.4509484628125742e-01 + 1.0349426984742664e-16j,
+                6.3637197288239186e-17 - 7.4398922703555097e-33j,
+                2.6734481396039929e-32 + 1.7155650257063576e-35j,
+            ],
+            [
+                6.3637197288239186e-17 - 7.4398922703555097e-33j,
+                -2.0594021562561332e-01 + 2.2863956908382538e-17j,
+                -5.8325863096557049e-17 + 1.6949718400585382e-18j,
+            ],
+            [
+                2.6734481396039929e-32 + 1.7155650257063576e-35j,
+                -5.8325863096557049e-17 + 1.6949718400585382e-18j,
+                4.4171453199503476e-02 + 1.0022350742842835e-02j,
+            ],
+        ]
+    )
+    d, U = takagi(A, svd_order=svd_order)
+    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
     correct decomposition when input a real but highly degenerate matrix. This test uses the

thewalrus/tests/test_symplectic.py

@@ -367,6 +367,15 @@ def test_tms(self, hbar, tol):
         assert np.allclose(res[0], expected[0], atol=tol, rtol=0)
         assert np.allclose(res[1], expected[1], atol=tol, rtol=0)
 
+    def test_ndarray(self, hbar, tol):
+        """Test numpy.ndarray in the third argument of `reduced_state` is converted to list correctly"""
+        mu, cov = symplectic.vacuum_state(4, hbar=hbar)
+        res = symplectic.reduced_state(mu, cov, np.array([0, 1, 2, 3]))
+        expected = np.zeros([8]), np.identity(8) * hbar / 2
+
+        assert np.allclose(res[0], expected[0], atol=tol, rtol=0)
+        assert np.allclose(res[1], expected[1], atol=tol, rtol=0)
+
 
 class TestLossChannel:
     """Tests for the loss channel"""