diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md
index 66dc415a..392e5d28 100644
--- a/.github/CHANGELOG.md
+++ b/.github/CHANGELOG.md
@@ -15,6 +15,8 @@
 
 ### Bug fixes
 
+* Add the calculation method of `takagi` when the matrix is diagonal. [(#394)](https://github.com/XanaduAI/thewalrus/pull/394)
+
 ### Documentation
 
 ### Contributors
diff --git a/thewalrus/decompositions.py b/thewalrus/decompositions.py
index 46585429..a0e65733 100644
--- a/thewalrus/decompositions.py
+++ b/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:
diff --git a/thewalrus/tests/test_decompositions.py b/thewalrus/tests/test_decompositions.py
index 122c0821..18e3adf6 100644
--- a/thewalrus/tests/test_decompositions.py
+++ b/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