Repository: quantumlib/OpenFermion
Pull request #91 (2017-10-30 02:40:58)

Fixes for Slater determinant preparation
The main things I changed were:
- The canonical form for antisymmetric matrices should have nonnegative entries in the upper right block; I forgot to do this initially.
- Particle-hole transformations used in the preparation of fermionic Gaussian states are now on the last site, rather than the first. The particle-hole transformation on the last site is easier to implement on an actual device.


Patch:
diff --git a/src/openfermion/utils/_slater_determinants.py b/src/openfermion/utils/_slater_determinants.py
index 5341c96..5ebecd6 100644
--- a/src/openfermion/utils/_slater_determinants.py
+++ b/src/openfermion/utils/_slater_determinants.py
@@ -22,7 +22,7 @@ from openfermion.config import EQ_TOLERANCE
 
 def fermionic_gaussian_decomposition(unitary_rows):
     """Decompose a matrix into a sequence of Givens rotations and
-    particle-hole transformations on the first fermionic mode.
+    particle-hole transformations on the last fermionic mode.
 
     The input is an n x (2 * n) matrix W with orthonormal rows.
     Furthermore, W has the block form::
@@ -36,17 +36,19 @@ def fermionic_gaussian_decomposition(unitary_rows):
 
     W can be decomposed as::
 
-        V * W * U^\dagger = [ 0  |  L ]
+        V * W * U^\dagger = [ 0  |  D ]
 
-    where V and U are unitary matrices and L is an antidiagonal unitary matrix.
+    where V and U are unitary matrices and D is a diagonal unitary matrix.
     Furthermore, we can decompose U as a sequence of Givens rotations
-    and particle-hole transformations on the first fermionic mode.
+    and particle-hole transformations on the last fermionic mode.
+    This particle-hole transformation maps a^\dagger_n to a^n and vice
+    versa, while leaving the other ladder operators invariant.
 
     The decomposition of U is returned as a list of tuples of objects
     describing rotations and particle-hole transformations. The list looks
-    something like [('p-h', ), (G_1, ), ('p-h', G_2), ... ].
-    The objects within a tuple are either the string 'p-h', which indicates
-    a particle-hole transformation on the first fermionic mode, or a tuple
+    something like [('pht', ), (G_1, ), ('pht', G_2), ... ].
+    The objects within a tuple are either the string 'pht', which indicates
+    a particle-hole transformation on the last fermionic mode, or a tuple
     of the form (i, j, theta, phi), which indicates a Givens roation
     of rows i and j by angles theta and phi.
 
@@ -57,7 +59,7 @@ def fermionic_gaussian_decomposition(unitary_rows):
     Returns:
         left_unitary(ndarray): An n x n matrix representing V.
         decomposition(list[tuple]): The decomposition of U.
-        antidiagonal(ndarray): A list of the nonzero entries of L.
+        diagonal(ndarray): A list of the nonzero entries of D.
     """
     current_matrix = numpy.copy(unitary_rows)
     n, p = current_matrix.shape
@@ -85,9 +87,9 @@ def fermionic_gaussian_decomposition(unitary_rows):
 
     # Compute left_unitary using Givens rotations
     left_unitary = numpy.eye(n, dtype=complex)
-    for k in reversed(range(1, n)):
+    for k in range(n - 1):
         # Zero out entries in column k
-        for l in range(k):
+        for l in range(n - 1 - k):
             # Zero out entry in row l
             givens_rotation = givens_matrix_elements(current_matrix[l, k],
                                                      current_matrix[l + 1, k])
@@ -104,46 +106,46 @@ def fermionic_gaussian_decomposition(unitary_rows):
         parallel_ops = list()
 
         # Perform a particle-hole transformation if necessary
-        if k % 2 == 0 and abs(current_matrix[k // 2, 0]) > EQ_TOLERANCE:
-            parallel_ops.append('p-h')
-            swap_columns(current_matrix, 0, n)
+        if k % 2 == 0 and abs(current_matrix[k // 2, n - 1]) > EQ_TOLERANCE:
+            parallel_ops.append('pht')
+            swap_columns(current_matrix, n - 1, 2 * n - 1)
 
         # Get the (row, column) indices of elements to zero out in parallel.
         if k < n:
             end_row = k
-            end_column = k
+            end_column = n - 1 - k
         else:
             end_row = n - 1
-            end_column = 2 * (n - 1) - k
-        column_indices = range(end_column, 0, -2)
+            end_column = k - (n - 1)
+        column_indices = range(end_column, n - 1, 2)
         row_indices = range(end_row, end_row - len(column_indices), -1)
         indices_to_zero_out = zip(row_indices, column_indices)
 
         for i, j in indices_to_zero_out:
-            # Compute the Givens rotation to zero out the (i, j) element
-            left_element = current_matrix[i, j - 1].conj()
-            right_element = current_matrix[i, j].conj()
-            givens_rotation = givens_matrix_elements(left_element,
-                                                     right_element)
-            # Need to switch the rows to zero out right_element
-            # rather than left_element
-            givens_rotation = givens_rotation[(1, 0), :]
-
-            # Add the parameters to the list
-            theta = numpy.arccos(numpy.real(givens_rotation[0, 0]))
-            phi = numpy.angle(givens_rotation[1, 1])
-            parallel_ops.append((j - 1, j, theta, phi))
-
-            # Update the matrix
-            double_givens_rotate(current_matrix, givens_rotation,
-                                 j - 1, j, which='col')
+            # Compute the Givens rotation to zero out the (i, j) element,
+            # if needed
+            left_element = current_matrix[i, j].conj()
+            if abs(left_element) > EQ_TOLERANCE:
+                # We actually need to perform a Givens rotation
+                right_element = current_matrix[i, j + 1].conj()
+                givens_rotation = givens_matrix_elements(left_element,
+                                                         right_element)
+
+                # Add the parameters to the list
+                theta = numpy.arccos(numpy.real(givens_rotation[0, 0]))
+                phi = numpy.angle(givens_rotation[1, 1])
+                parallel_ops.append((j, j + 1, theta, phi))
+
+                # Update the matrix
+                double_givens_rotate(current_matrix, givens_rotation,
+                                     j, j + 1, which='col')
 
         # Append the current list of parallel rotations to the list
         decomposition.append(tuple(parallel_ops))
 
-    # Get the antidiagonal entries
-    antidiagonal = current_matrix[range(n), range(2 * n - 1, n - 1, -1)]
-    return left_unitary, decomposition, antidiagonal
+    # Get the diagonal entries
+    diagonal = current_matrix[range(n), range(n, 2 * n)]
+    return left_unitary, decomposition, diagonal
 
 
 def diagonalizing_fermionic_unitary(antisymmetric_matrix):
@@ -202,7 +204,7 @@ def antisymmetric_canonical_form(antisymmetric_matrix):
         [  0     D ]
         [ -D     0 ]
 
-    where D is a diagonal matrix.
+    where D is a diagonal matrix with nonnegative entries.
 
     Args:
         antisymmetric_matrix(ndarray): An antisymmetric matrix with even
@@ -239,6 +241,13 @@ def antisymmetric_canonical_form(antisymmetric_matrix):
             swap_columns(canonical, num_blocks - 1, num_blocks + i)
             swap_columns(orthogonal, num_blocks - 1, num_blocks + i)
 
+    # Now we permute so that the upper right block is non-negative
+    for i in range(num_blocks):
+        if canonical[i, num_blocks + i] < -EQ_TOLERANCE:
+            swap_rows(canonical, i, num_blocks + i)
+            swap_columns(canonical, i, num_blocks + i)
+            swap_columns(orthogonal, i, num_blocks + i)
+
     return canonical, orthogonal.T
 
 
@@ -340,23 +349,26 @@ def givens_decomposition(unitary_rows):
 
             parallel_rotations = list()
             for i, j in indices_to_zero_out:
-                # Compute the Givens rotation to zero out the (i, j) element
-                left_element = current_matrix[i, j - 1].conj()
+                # Compute the Givens rotation to zero out the (i, j) element,
+                # if needed
                 right_element = current_matrix[i, j].conj()
-                givens_rotation = givens_matrix_elements(left_element,
-                                                         right_element)
-                # Need to switch the rows to zero out right_element
-                # rather than left_element
-                givens_rotation = givens_rotation[(1, 0), :]
-
-                # Add the parameters to the list
-                theta = numpy.arccos(numpy.real(givens_rotation[0, 0]))
-                phi = numpy.angle(givens_rotation[1, 1])
-                parallel_rotations.append((j - 1, j, theta, phi))
-
-                # Update the matrix
-                givens_rotate(current_matrix, givens_rotation,
-                              j - 1, j, which='col')
+                if abs(right_element) > EQ_TOLERANCE:
+                    # We actually need to perform a Givens rotation
+                    left_element = current_matrix[i, j - 1].conj()
+                    givens_rotation = givens_matrix_elements(left_element,
+                                                             right_element)
+                    # Need to switch the rows to zero out right_element
+                    # rather than left_element
+                    givens_rotation = givens_rotation[(1, 0), :]
+
+                    # Add the parameters to the list
+                    theta = numpy.arccos(numpy.real(givens_rotation[0, 0]))
+                    phi = numpy.angle(givens_rotation[1, 1])
+                    parallel_rotations.append((j - 1, j, theta, phi))
+
+                    # Update the matrix
+                    givens_rotate(current_matrix, givens_rotation,
+                                  j - 1, j, which='col')
 
             # Append the current list of parallel rotations to the list
             givens_rotations.append(tuple(parallel_rotations))
@@ -393,12 +405,10 @@ def givens_matrix_elements(a, b):
     elif abs(b) < EQ_TOLERANCE:
         c = 0.
         s = 1.
-        sign_a = a / abs(a)
-        phase = sign_a
+        phase = 1.
     # Handle case that a and b are both nonzero
     else:
         denominator = numpy.sqrt(abs(a) ** 2 + abs(b) ** 2)
-
         c = abs(b) / denominator
         s = abs(a) / denominator
         sign_b = b / abs(b)


Test patch:
diff --git a/src/openfermion/utils/_slater_determinants_test.py b/src/openfermion/utils/_slater_determinants_test.py
index 73ced28..d0c838c 100644
--- a/src/openfermion/utils/_slater_determinants_test.py
+++ b/src/openfermion/utils/_slater_determinants_test.py
@@ -17,10 +17,14 @@ import numpy
 import unittest
 from scipy.linalg import qr
 
+from openfermion.config import EQ_TOLERANCE
 from openfermion.utils import (fermionic_gaussian_decomposition,
                                givens_decomposition)
 from openfermion.utils._slater_determinants import (
-        diagonalizing_fermionic_unitary, double_givens_rotate, givens_rotate,
+        antisymmetric_canonical_form,
+        diagonalizing_fermionic_unitary,
+        double_givens_rotate,
+        givens_rotate,
         swap_rows)
 
 
@@ -427,16 +431,23 @@ class FermionicGaussianDecompositionTest(unittest.TestCase):
         lower_unitary = ferm_unitary[n:]
 
         # Get fermionic Gaussian decomposition of lower_unitary
-        left_unitary, decomposition, antidiagonal = (
+        left_unitary, decomposition, diagonal = (
                 fermionic_gaussian_decomposition(lower_unitary))
 
+        # Check that left_unitary zeroes out the correct entries of
+        # lower_unitary
+        product = left_unitary.dot(lower_unitary)
+        for i in range(n - 1):
+            for j in range(n - 1 - i):
+                self.assertAlmostEqual(product[i, j], 0.)
+
         # Compute right_unitary
         right_unitary = numpy.eye(2 * n, dtype=complex)
         for parallel_set in decomposition:
             combined_op = numpy.eye(2 * n, dtype=complex)
             for op in parallel_set:
-                if op == 'p-h':
-                    swap_rows(combined_op, 0, n)
+                if op == 'pht':
+                    swap_rows(combined_op, n - 1, 2 * n - 1)
                 else:
                     i, j, theta, phi = op
                     c = numpy.cos(theta)
@@ -451,13 +462,13 @@ class FermionicGaussianDecompositionTest(unittest.TestCase):
         # Compute left_unitary * lower_unitary * right_unitary^\dagger
         product = left_unitary.dot(lower_unitary.dot(right_unitary.T.conj()))
 
-        # Construct the antidiagonal matrix
-        anti_diag = numpy.zeros((n, 2 * n), dtype=complex)
-        anti_diag[range(n), range(2 * n - 1, n - 1, -1)] = antidiagonal
+        # Construct the diagonal matrix
+        diag = numpy.zeros((n, 2 * n), dtype=complex)
+        diag[range(n), range(n, 2 * n)] = diagonal
 
         # Assert that W and D are the same
         for i in numpy.ndindex((n, 2 * n)):
-            self.assertAlmostEqual(anti_diag[i], product[i])
+            self.assertAlmostEqual(diag[i], product[i])
 
     def test_n_equals_4(self):
         n = 4
@@ -470,16 +481,23 @@ class FermionicGaussianDecompositionTest(unittest.TestCase):
         lower_unitary = ferm_unitary[n:]
 
         # Get fermionic Gaussian decomposition of lower_unitary
-        left_unitary, decomposition, antidiagonal = (
+        left_unitary, decomposition, diagonal = (
                 fermionic_gaussian_decomposition(lower_unitary))
 
+        # Check that left_unitary zeroes out the correct entries of
+        # lower_unitary
+        product = left_unitary.dot(lower_unitary)
+        for i in range(n - 1):
+            for j in range(n - 1 - i):
+                self.assertAlmostEqual(product[i, j], 0.)
+
         # Compute right_unitary
         right_unitary = numpy.eye(2 * n, dtype=complex)
         for parallel_set in decomposition:
             combined_op = numpy.eye(2 * n, dtype=complex)
             for op in parallel_set:
-                if op == 'p-h':
-                    swap_rows(combined_op, 0, n)
+                if op == 'pht':
+                    swap_rows(combined_op, n - 1, 2 * n - 1)
                 else:
                     i, j, theta, phi = op
                     c = numpy.cos(theta)
@@ -494,13 +512,13 @@ class FermionicGaussianDecompositionTest(unittest.TestCase):
         # Compute left_unitary * lower_unitary * right_unitary^\dagger
         product = left_unitary.dot(lower_unitary.dot(right_unitary.T.conj()))
 
-        # Construct the antidiagonal matrix
-        anti_diag = numpy.zeros((n, 2 * n), dtype=complex)
-        anti_diag[range(n), range(2 * n - 1, n - 1, -1)] = antidiagonal
+        # Construct the diagonal matrix
+        diag = numpy.zeros((n, 2 * n), dtype=complex)
+        diag[range(n), range(n, 2 * n)] = diagonal
 
         # Assert that W and D are the same
         for i in numpy.ndindex((n, 2 * n)):
-            self.assertAlmostEqual(anti_diag[i], product[i])
+            self.assertAlmostEqual(diag[i], product[i])
 
     def test_n_equals_5(self):
         n = 5
@@ -513,16 +531,23 @@ class FermionicGaussianDecompositionTest(unittest.TestCase):
         lower_unitary = ferm_unitary[n:]
 
         # Get fermionic Gaussian decomposition of lower_unitary
-        left_unitary, decomposition, antidiagonal = (
+        left_unitary, decomposition, diagonal = (
                 fermionic_gaussian_decomposition(lower_unitary))
 
+        # Check that left_unitary zeroes out the correct entries of
+        # lower_unitary
+        product = left_unitary.dot(lower_unitary)
+        for i in range(n - 1):
+            for j in range(n - 1 - i):
+                self.assertAlmostEqual(product[i, j], 0.)
+
         # Compute right_unitary
         right_unitary = numpy.eye(2 * n, dtype=complex)
         for parallel_set in decomposition:
             combined_op = numpy.eye(2 * n, dtype=complex)
             for op in parallel_set:
-                if op == 'p-h':
-                    swap_rows(combined_op, 0, n)
+                if op == 'pht':
+                    swap_rows(combined_op, n - 1, 2 * n - 1)
                 else:
                     i, j, theta, phi = op
                     c = numpy.cos(theta)
@@ -537,13 +562,13 @@ class FermionicGaussianDecompositionTest(unittest.TestCase):
         # Compute left_unitary * lower_unitary * right_unitary^\dagger
         product = left_unitary.dot(lower_unitary.dot(right_unitary.T.conj()))
 
-        # Construct the antidiagonal matrix
-        anti_diag = numpy.zeros((n, 2 * n), dtype=complex)
-        anti_diag[range(n), range(2 * n - 1, n - 1, -1)] = antidiagonal
+        # Construct the diagonal matrix
+        diag = numpy.zeros((n, 2 * n), dtype=complex)
+        diag[range(n), range(n, 2 * n)] = diagonal
 
         # Assert that W and D are the same
         for i in numpy.ndindex((n, 2 * n)):
-            self.assertAlmostEqual(anti_diag[i], product[i])
+            self.assertAlmostEqual(diag[i], product[i])
 
     def test_n_equals_6(self):
         n = 6
@@ -556,16 +581,23 @@ class FermionicGaussianDecompositionTest(unittest.TestCase):
         lower_unitary = ferm_unitary[n:]
 
         # Get fermionic Gaussian decomposition of lower_unitary
-        left_unitary, decomposition, antidiagonal = (
+        left_unitary, decomposition, diagonal = (
                 fermionic_gaussian_decomposition(lower_unitary))
 
+        # Check that left_unitary zeroes out the correct entries of
+        # lower_unitary
+        product = left_unitary.dot(lower_unitary)
+        for i in range(n - 1):
+            for j in range(n - 1 - i):
+                self.assertAlmostEqual(product[i, j], 0.)
+
         # Compute right_unitary
         right_unitary = numpy.eye(2 * n, dtype=complex)
         for parallel_set in decomposition:
             combined_op = numpy.eye(2 * n, dtype=complex)
             for op in parallel_set:
-                if op == 'p-h':
-                    swap_rows(combined_op, 0, n)
+                if op == 'pht':
+                    swap_rows(combined_op, n - 1, 2 * n - 1)
                 else:
                     i, j, theta, phi = op
                     c = numpy.cos(theta)
@@ -580,13 +612,13 @@ class FermionicGaussianDecompositionTest(unittest.TestCase):
         # Compute left_unitary * lower_unitary * right_unitary^\dagger
         product = left_unitary.dot(lower_unitary.dot(right_unitary.T.conj()))
 
-        # Construct the antidiagonal matrix
-        anti_diag = numpy.zeros((n, 2 * n), dtype=complex)
-        anti_diag[range(n), range(2 * n - 1, n - 1, -1)] = antidiagonal
+        # Construct the diagonal matrix
+        diag = numpy.zeros((n, 2 * n), dtype=complex)
+        diag[range(n), range(n, 2 * n)] = diagonal
 
         # Assert that W and D are the same
         for i in numpy.ndindex((n, 2 * n)):
-            self.assertAlmostEqual(anti_diag[i], product[i])
+            self.assertAlmostEqual(diag[i], product[i])
 
     def test_n_equals_7(self):
         n = 7
@@ -599,16 +631,23 @@ class FermionicGaussianDecompositionTest(unittest.TestCase):
         lower_unitary = ferm_unitary[n:]
 
         # Get fermionic Gaussian decomposition of lower_unitary
-        left_unitary, decomposition, antidiagonal = (
+        left_unitary, decomposition, diagonal = (
                 fermionic_gaussian_decomposition(lower_unitary))
 
+        # Check that left_unitary zeroes out the correct entries of
+        # lower_unitary
+        product = left_unitary.dot(lower_unitary)
+        for i in range(n - 1):
+            for j in range(n - 1 - i):
+                self.assertAlmostEqual(product[i, j], 0.)
+
         # Compute right_unitary
         right_unitary = numpy.eye(2 * n, dtype=complex)
         for parallel_set in decomposition:
             combined_op = numpy.eye(2 * n, dtype=complex)
             for op in parallel_set:
-                if op == 'p-h':
-                    swap_rows(combined_op, 0, n)
+                if op == 'pht':
+                    swap_rows(combined_op, n - 1, 2 * n - 1)
                 else:
                     i, j, theta, phi = op
                     c = numpy.cos(theta)
@@ -623,13 +662,13 @@ class FermionicGaussianDecompositionTest(unittest.TestCase):
         # Compute left_unitary * lower_unitary * right_unitary^\dagger
         product = left_unitary.dot(lower_unitary.dot(right_unitary.T.conj()))
 
-        # Construct the antidiagonal matrix
-        anti_diag = numpy.zeros((n, 2 * n), dtype=complex)
-        anti_diag[range(n), range(2 * n - 1, n - 1, -1)] = antidiagonal
+        # Construct the diagonal matrix
+        diag = numpy.zeros((n, 2 * n), dtype=complex)
+        diag[range(n), range(n, 2 * n)] = diagonal
 
         # Assert that W and D are the same
         for i in numpy.ndindex((n, 2 * n)):
-            self.assertAlmostEqual(anti_diag[i], product[i])
+            self.assertAlmostEqual(diag[i], product[i])
 
     def test_n_equals_8(self):
         n = 8
@@ -642,16 +681,23 @@ class FermionicGaussianDecompositionTest(unittest.TestCase):
         lower_unitary = ferm_unitary[n:]
 
         # Get fermionic Gaussian decomposition of lower_unitary
-        left_unitary, decomposition, antidiagonal = (
+        left_unitary, decomposition, diagonal = (
                 fermionic_gaussian_decomposition(lower_unitary))
 
+        # Check that left_unitary zeroes out the correct entries of
+        # lower_unitary
+        product = left_unitary.dot(lower_unitary)
+        for i in range(n - 1):
+            for j in range(n - 1 - i):
+                self.assertAlmostEqual(product[i, j], 0.)
+
         # Compute right_unitary
         right_unitary = numpy.eye(2 * n, dtype=complex)
         for parallel_set in decomposition:
             combined_op = numpy.eye(2 * n, dtype=complex)
             for op in parallel_set:
-                if op == 'p-h':
-                    swap_rows(combined_op, 0, n)
+                if op == 'pht':
+                    swap_rows(combined_op, n - 1, 2 * n - 1)
                 else:
                     i, j, theta, phi = op
                     c = numpy.cos(theta)
@@ -666,13 +712,13 @@ class FermionicGaussianDecompositionTest(unittest.TestCase):
         # Compute left_unitary * lower_unitary * right_unitary^\dagger
         product = left_unitary.dot(lower_unitary.dot(right_unitary.T.conj()))
 
-        # Construct the antidiagonal matrix
-        anti_diag = numpy.zeros((n, 2 * n), dtype=complex)
-        anti_diag[range(n), range(2 * n - 1, n - 1, -1)] = antidiagonal
+        # Construct the diagonal matrix
+        diag = numpy.zeros((n, 2 * n), dtype=complex)
+        diag[range(n), range(n, 2 * n)] = diagonal
 
         # Assert that W and D are the same
         for i in numpy.ndindex((n, 2 * n)):
-            self.assertAlmostEqual(anti_diag[i], product[i])
+            self.assertAlmostEqual(diag[i], product[i])
 
     def test_n_equals_9(self):
         n = 9
@@ -685,16 +731,23 @@ class FermionicGaussianDecompositionTest(unittest.TestCase):
         lower_unitary = ferm_unitary[n:]
 
         # Get fermionic Gaussian decomposition of lower_unitary
-        left_unitary, decomposition, antidiagonal = (
+        left_unitary, decomposition, diagonal = (
                 fermionic_gaussian_decomposition(lower_unitary))
 
+        # Check that left_unitary zeroes out the correct entries of
+        # lower_unitary
+        product = left_unitary.dot(lower_unitary)
+        for i in range(n - 1):
+            for j in range(n - 1 - i):
+                self.assertAlmostEqual(product[i, j], 0.)
+
         # Compute right_unitary
         right_unitary = numpy.eye(2 * n, dtype=complex)
         for parallel_set in decomposition:
             combined_op = numpy.eye(2 * n, dtype=complex)
             for op in parallel_set:
-                if op == 'p-h':
-                    swap_rows(combined_op, 0, n)
+                if op == 'pht':
+                    swap_rows(combined_op, n - 1, 2 * n - 1)
                 else:
                     i, j, theta, phi = op
                     c = numpy.cos(theta)
@@ -709,13 +762,13 @@ class FermionicGaussianDecompositionTest(unittest.TestCase):
         # Compute left_unitary * lower_unitary * right_unitary^\dagger
         product = left_unitary.dot(lower_unitary.dot(right_unitary.T.conj()))
 
-        # Construct the antidiagonal matrix
-        anti_diag = numpy.zeros((n, 2 * n), dtype=complex)
-        anti_diag[range(n), range(2 * n - 1, n - 1, -1)] = antidiagonal
+        # Construct the diagonal matrix
+        diag = numpy.zeros((n, 2 * n), dtype=complex)
+        diag[range(n), range(n, 2 * n)] = diagonal
 
         # Assert that W and D are the same
         for i in numpy.ndindex((n, 2 * n)):
-            self.assertAlmostEqual(anti_diag[i], product[i])
+            self.assertAlmostEqual(diag[i], product[i])
 
 
 class DiagonalizingFermionicUnitaryTest(unittest.TestCase):
@@ -756,3 +809,34 @@ class DiagonalizingFermionicUnitaryTest(unittest.TestCase):
         for i in numpy.ndindex((n, n)):
             self.assertAlmostEqual(identity[i], constraint_matrix_1[i])
             self.assertAlmostEqual(0., constraint_matrix_2[i])
+
+
+class AntisymmetricCanonicalFormTest(unittest.TestCase):
+
+    def test_equality(self):
+        """Test that the decomposition is valid."""
+        n = 7
+        rand_mat = numpy.random.randn(2 * n, 2 * n)
+        antisymmetric_matrix = rand_mat - rand_mat.T
+        canonical, orthogonal = antisymmetric_canonical_form(
+                antisymmetric_matrix)
+        result_matrix = orthogonal.dot(antisymmetric_matrix.dot(orthogonal.T))
+        for i in numpy.ndindex(result_matrix.shape):
+            self.assertAlmostEqual(result_matrix[i], canonical[i])
+
+    def test_canonical(self):
+        """Test that the returned canonical matrix has the right form."""
+        n = 7
+        # Obtain a random antisymmetric matrix
+        rand_mat = numpy.random.randn(2 * n, 2 * n)
+        antisymmetric_matrix = rand_mat - rand_mat.T
+        canonical, orthogonal = antisymmetric_canonical_form(
+                antisymmetric_matrix)
+        for i in range(2 * n):
+            for j in range(2 * n):
+                if i < n and j == n + i:
+                    self.assertTrue(canonical[i, j] > -EQ_TOLERANCE)
+                elif i >= n and j == i - n:
+                    self.assertTrue(canonical[i, j] < -EQ_TOLERANCE)
+                else:
+                    self.assertAlmostEqual(canonical[i, j], 0.)