diff --git a/.github/ACKNOWLEDGMENTS.md b/.github/ACKNOWLEDGMENTS.md
index c9ba3d427..7281c8060 100644
--- a/.github/ACKNOWLEDGMENTS.md
+++ b/.github/ACKNOWLEDGMENTS.md
@@ -62,4 +62,6 @@
 
 * [Fabian Laudenbach](https://github.com/fab1an-q) (Xanadu) - :revolving_hearts: Entangler of hearts
 
-* [Martin Houde](https://github.com/MHoude2) (Polytechnique Montréal) - :upside_down: Minister of amplification
\ No newline at end of file
+* [Martin Houde](https://github.com/MHoude2) (Polytechnique Montréal) - 🙃 Minister of amplification
+
+* Will McCutcheon (Heriot-Watt University) - 🧅 Gaussian Onion Merchant
diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md
index 5dfc61bc5..a8875f9bf 100644
--- a/.github/CHANGELOG.md
+++ b/.github/CHANGELOG.md
@@ -2,6 +2,8 @@
 
 ### New features
 
+* New function to produce Bloch-Messiah decomposition of symplectic matrices.
+
 ### Breaking changes
 
 ### Improvements
@@ -22,7 +24,7 @@
 
 This release contains contributions from (in alphabetical order):
 
-Mikhail Andrenkov, Antonín Hoskovec
+Mikhail Andrenkov, Antonín Hoskovec, Martin Houde
 
 ---
 
diff --git a/thewalrus/decompositions.py b/thewalrus/decompositions.py
index f1f5c462f..e179ec43f 100644
--- a/thewalrus/decompositions.py
+++ b/thewalrus/decompositions.py
@@ -28,6 +28,7 @@
 .. autosummary::
     williamson
     symplectic_eigenvals
+    blochmessiah
 
 Code details
 ------------
@@ -36,6 +37,7 @@
 
 from scipy.linalg import block_diag, sqrtm, schur
 from thewalrus.symplectic import sympmat
+from thewalrus.quantum.gaussian_checks import is_symplectic
 
 
 def williamson(V, rtol=1e-05, atol=1e-08):
@@ -112,3 +114,64 @@ def symplectic_eigenvals(cov):
     M = int(len(cov) / 2)
     D, _ = williamson(cov)
     return np.diag(D)[:M]
+
+
+def blochmessiah(S):
+    """Returns the Bloch-Messiah decomposition of a symplectic matrix S = uff @ dff @ vff
+       where uff and vff are orthogonal symplectic matrices and dff is a diagonal matrix
+       of the form diag(d1,d2,...,dn,d1^-1, d2^-1,...,dn^-1),
+
+    Args:
+        S (array[float]): 2N x 2N real symplectic matrix
+
+    Returns:
+        tupple(array[float],  : orthogonal symplectic matrix uff
+               array[float],  : diagional matrix dff
+               array[float])  : orthogonal symplectic matrix vff
+    """
+
+    N, m = S.shape
+
+    if N != m:
+        return False
+    if N % 2 != 0:
+        return False
+    if not is_symplectic(S):
+        return False
+    # Changing Basis
+    R = (1 / np.sqrt(2)) * np.block(
+        [[np.eye(N // 2), 1j * np.eye(N // 2)], [np.eye(N // 2), -1j * np.eye(N // 2)]]
+    )
+    Sc = R @ S @ np.conjugate(R).T
+    # Polar Decomposition
+    u1, d1, v1 = np.linalg.svd(Sc)
+    Sig = u1 @ np.diag(d1) @ np.conjugate(u1).T
+    Unitary = u1 @ v1
+    # Blocks of Unitary and Hermitian symplectics
+    alpha = Unitary[0 : N // 2, 0 : N // 2]
+    beta = Sig[0 : N // 2, N // 2 : N]
+    # Bloch-Messiah in this Basis
+    u2, d2, v2 = np.linalg.svd(beta)
+    sval = np.arcsinh(d2)
+    takagibeta = u2 @ sqrtm(np.conjugate(u2).T @ (v2.T))
+    uf = np.block([[takagibeta, 0 * takagibeta], [0 * takagibeta, np.conjugate(takagibeta)]])
+    vf = np.block(
+        [
+            [np.conjugate(takagibeta).T @ alpha, 0 * takagibeta],
+            [0 * takagibeta, np.conjugate(np.conjugate(takagibeta).T @ alpha)],
+        ]
+    )
+    df = np.block(
+        [
+            [np.diag(np.cosh(sval)), np.diag(np.sinh(sval))],
+            [np.diag(np.sinh(sval)), np.diag(np.cosh(sval))],
+        ]
+    )
+    # Rotating Back to Original Basis
+    uff = np.conjugate(R).T @ uf @ R
+    vff = np.conjugate(R).T @ vf @ R
+    dff = np.conjugate(R).T @ df @ R
+    dff = np.real_if_close(dff)
+    vff = np.real_if_close(vff)
+    uff = np.real_if_close(uff)
+    return uff, dff, vff
diff --git a/thewalrus/tests/test_decompositions.py b/thewalrus/tests/test_decompositions.py
index 77b47b02b..d3882da71 100644
--- a/thewalrus/tests/test_decompositions.py
+++ b/thewalrus/tests/test_decompositions.py
@@ -18,8 +18,10 @@
 from scipy.linalg import block_diag
 
 from thewalrus.random import random_interferometer as haar_measure
-from thewalrus.decompositions import williamson
+from thewalrus.random import random_symplectic
+from thewalrus.decompositions import williamson, blochmessiah
 from thewalrus.symplectic import sympmat as omega
+from thewalrus.quantum.gaussian_checks import is_symplectic
 
 
 class TestWilliamsonDecomposition:
@@ -147,3 +149,33 @@ def test_mixed_state(self, create_cov, hbar, tol):
         assert np.allclose(sorted(nbar[:n]), sorted(nbar_in), atol=tol, rtol=0)
         # check S is symplectic
         assert np.allclose(S @ O @ S.T, O, atol=tol, rtol=0)
+
+
+class TestBlochMessiahDecomposition:
+    """Tests for the Williamson decomposition"""
+
+    @pytest.mark.parametrize("N", range(50, 500, 50))
+    def test_blochmessiah_rand(self, N):
+        """Tests blochmessiah function for different matrix sizes."""
+        S = random_symplectic(N)
+        u, d, v = blochmessiah(S)
+        assert np.allclose(u @ d @ v, S)
+        np.allclose(u.T @ u, np.eye(len(u)))
+        np.allclose(v.T @ v, np.eye(len(v)))
+        is_symplectic(u)
+        is_symplectic(v)
+
+    def test_blochmessiah_odd(self):
+        """Tests that odd matrices return False in blochmessiah."""
+        S = np.random.rand(5, 5)
+        assert not blochmessiah(S)
+
+    def test_blochmessiah_rect(self):
+        """Tests that rectangular matrices return False in blochmessiah"""
+        S = np.random.rand(4, 5)
+        assert not blochmessiah(S)
+
+    def test_blochmessiah_false(self):
+        """Tests that non-symplectic mattrices return False in blochmessiah"""
+        S = np.random.rand(4, 4)
+        assert not blochmessiah(S)