diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md
index 66dc415a..09b945ba 100644
--- a/.github/CHANGELOG.md
+++ b/.github/CHANGELOG.md
@@ -15,6 +15,8 @@
 
 ### Bug fixes
 
+* Added tolerance parameter to `deecompositions.takagi` to prevent from incorrect Autonne-Takagi decomposition, which happened with some matrices. [(#393)](https://github.com/XanaduAI/thewalrus/pull/393)
+
 ### Documentation
 
 ### Contributors
diff --git a/thewalrus/decompositions.py b/thewalrus/decompositions.py
index 78f2575c..3ca4dd29 100644
--- a/thewalrus/decompositions.py
+++ b/thewalrus/decompositions.py
@@ -152,7 +152,7 @@ def blochmessiah(S):
     return O, D, Q
 
 
-def takagi(A, svd_order=True):
+def takagi(A, svd_order=True, tol=1e-16):
     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 +162,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.
+        tol (float): tolerance parameter to kill small values in the matrix A. This was installed because takagi decomposition did not work at some matrix A.
 
     Returns:
         tuple[array, array]: (r, U), where r are the singular values,
@@ -174,7 +175,9 @@ def takagi(A, svd_order=True):
     # 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)
+    # Kill small values in the matrix A. This was installed because takagi decomposition did not work at some matrix A with small values.
+    rtol = tol * np.amax(np.abs(A))
+    A = np.where(np.abs(A.real) <= rtol, 0, A.real) + np.where(np.abs(A.imag) <= rtol, 0, A.imag)
 
     if np.allclose(A, 0):
         return np.zeros(n), np.eye(n)
diff --git a/thewalrus/internal_modes/prepare_cov.py b/thewalrus/internal_modes/prepare_cov.py
index fc7e12a2..b4b368ec 100644
--- a/thewalrus/internal_modes/prepare_cov.py
+++ b/thewalrus/internal_modes/prepare_cov.py
@@ -230,7 +230,7 @@ def state_prep(eps, W, thresh=1e-4, hbar=2):
 
 def LO_overlaps(chis, LO_shape):
     r"""
-    Computes the overlap integral between the orthonormal moes and the local oscillator shape
+    Computes the overlap integral between the orthonormal modes and the local oscillator shape
 
     Args:
         chis (list[array]): list of the temporal functions for the new orthonormal basis