diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md
index 48c70c10..6eef08c9 100644
--- a/.github/CHANGELOG.md
+++ b/.github/CHANGELOG.md
@@ -5,6 +5,7 @@
 * Adds the Takagi decomposition [(#363)](https://github.com/XanaduAI/thewalrus/pull/363)
 
 * Adds the Montrealer and Loop Montrealer functions [(#363)](https://github.com/XanaduAI/thewalrus/pull/374).
+
 ### Breaking changes
 
 ### Improvements
@@ -18,6 +19,8 @@
 
 * Improves the handling of an edge case in Takagi [(#373)](https://github.com/XanaduAI/thewalrus/pull/373).
 
+* Adds extra tests for the Takagi decomposition [(#377)](https://github.com/XanaduAI/thewalrus/pull/377)
+
 ### Bug fixes
 
 ### Documentation
diff --git a/thewalrus/decompositions.py b/thewalrus/decompositions.py
index ad0e767b..f272b052 100644
--- a/thewalrus/decompositions.py
+++ b/thewalrus/decompositions.py
@@ -166,7 +166,8 @@ def blochmessiah(S):
 
 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.
+    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.
     See `Carl Caves note. <http://info.phys.unm.edu/~caves/courses/qinfo-s17/lectures/polarsingularAutonne.pdf>`_
 
     Args:
@@ -181,8 +182,8 @@ def takagi(A, svd_order=True):
     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)
+    # Here we build a Symmetric matrix from the top right triangular part
+    A = np.triu(A) + np.triu(A, k=1).T
 
     A = np.real_if_close(A)
 
@@ -192,24 +193,16 @@ def takagi(A, svd_order=True):
     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)]
+        signs = (-1) ** (1 + np.heaviside(ls, 1))
+        phases = np.sqrt(np.complex128(signs))
+        Uc = U * phases  # One needs to readjust the phases
+        # Find the permutation to sort in decreasing order
+        perm = np.argsort(vals)
+        # if svd_order reverse it
+        if svd_order:
+            perm = perm[::-1]
+        return vals[perm], Uc[:, perm]
 
     # Find the element with the largest absolute value
     pos = np.unravel_index(np.argmax(np.abs(A)), (n, n))
@@ -222,10 +215,6 @@ def takagi(A, svd_order=True):
 
     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
diff --git a/thewalrus/tests/test_decompositions.py b/thewalrus/tests/test_decompositions.py
index ab7deeff..d0f15d07 100644
--- a/thewalrus/tests/test_decompositions.py
+++ b/thewalrus/tests/test_decompositions.py
@@ -274,23 +274,30 @@ def test_takagi(n, datatype, svd_order):
         assert np.all(np.diff(r) >= 0)
 
 
+# pylint: disable=too-many-arguments
 @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):
+@pytest.mark.parametrize("null_space", [0, 5, 10])
+@pytest.mark.parametrize("offset", [0, 0.5])
+def test_degenerate(n, datatype, svd_order, half_rank, phase, null_space, offset):
     """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)
+    diags = (
+        [half_rank * np.random.rand()] * nhalf
+        + [np.random.rand() - offset] * (n - nhalf)
+        + [0] * null_space
+    )
     if datatype is np.complex128:
-        U = haar_measure(n)
+        U = haar_measure(n + null_space)
     if datatype is np.float64:
-        U = np.exp(1j * phase) * haar_measure(n, real=True)
+        U = np.exp(1j * phase) * haar_measure(n + null_space, 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.allclose(U @ U.T.conj(), np.eye(n + null_space))
     assert np.all(r >= 0)
     if svd_order is True:
         assert np.all(np.diff(r) <= 0)