Repository: quantumlib/OpenFermion
Pull request #104 (2017-11-07 16:22:48)

Procedure for ground state of quadratic Hamiltonians
- Added `jw_get_quadratic_hamiltonian_ground state`, which computes the ground energy and state of a quadratic Hamiltonian. It can be used instead of `get_ground_state` when possible; it should be faster and more memory efficient. It also serves as a test for the Slater determinant preparation circuits.
- Some minor optimizations and rewriting of Givens rotations code and tests
- Fixed a bug in `get_ground_state` which caused it to return the complex conjugate of the ground state instead of the true ground state.

kevinsung: I think `jw_get_quadratic_hamiltonian_ground_state` still fits better in `_slater_determinants.py` for now. However, I took your second suggestion.

Patch:
diff --git a/src/openfermion/ops/_quadratic_hamiltonian.py b/src/openfermion/ops/_quadratic_hamiltonian.py
index fdd42be..3d8cf0b 100644
--- a/src/openfermion/ops/_quadratic_hamiltonian.py
+++ b/src/openfermion/ops/_quadratic_hamiltonian.py
@@ -40,6 +40,10 @@ class QuadraticHamiltonian(PolynomialTensor):
 
     We separate the chemical potential \mu from M so that we can use it
     to adjust the expectation value of the total number of particles.
+
+    Attributes:
+        constant(float): A constant term in the operator.
+        chemical_potential(float): The chemical potential \mu.
     """
 
     def __init__(self, constant, hermitian_part,
@@ -67,6 +71,7 @@ class QuadraticHamiltonian(PolynomialTensor):
                     hermitian_part -
                     chemical_potential * numpy.eye(n_qubits))
 
+        # Initialize the PolynomialTensor
         if antisymmetric_part is None:
             super(QuadraticHamiltonian, self).__init__(
                     {(): constant,
@@ -78,6 +83,8 @@ class QuadraticHamiltonian(PolynomialTensor):
                      (1, 1): .5 * antisymmetric_part,
                      (0, 0): -.5 * antisymmetric_part.conj()})
 
+        # Add remaining attributes
+        self.constant = self.n_body_tensors[()]
         self.chemical_potential = chemical_potential
 
     def combined_hermitian_part(self):
@@ -113,8 +120,8 @@ class QuadraticHamiltonian(PolynomialTensor):
 
         where the s_i are normalized Majorana fermion operators:
 
-            s_j = i / sqrt(2) (a^\dagger_j - a_j)
-            s_{j + n_qubits} = 1 / sqrt(2) (a^\dagger_j + a_j)
+            s_j = 1 / sqrt(2) (a^\dagger_j + a_j)
+            s_{j + n_qubits} = i / sqrt(2) (a^\dagger_j - a_j)
 
         and A is a (2 * n_qubits) x (2 * n_qubits) real antisymmetric matrix.
         This function returns the matrix A and the constant.
diff --git a/src/openfermion/utils/__init__.py b/src/openfermion/utils/__init__.py
index d75656f..bb46d80 100644
--- a/src/openfermion/utils/__init__.py
+++ b/src/openfermion/utils/__init__.py
@@ -22,7 +22,8 @@ from ._operator_utils import (commutator, count_qubits,
 
 from ._slater_determinants import (fermionic_gaussian_decomposition,
                                    givens_decomposition,
-                                   ground_state_preparation_circuit)
+                                   ground_state_preparation_circuit,
+                                   jw_get_quadratic_hamiltonian_ground_state)
 
 from ._sparse_tools import (expectation,
                             expectation_computational_basis_state,
@@ -34,6 +35,7 @@ from ._sparse_tools import (expectation,
                             jw_hartree_fock_state,
                             jw_get_ground_states_by_particle_number,
                             jw_number_restrict_operator,
+                            jw_slater_determinant,
                             qubit_operator_sparse,
                             sparse_eigenspectrum)
 
diff --git a/src/openfermion/utils/_slater_determinants.py b/src/openfermion/utils/_slater_determinants.py
index d0d7a83..0485665 100644
--- a/src/openfermion/utils/_slater_determinants.py
+++ b/src/openfermion/utils/_slater_determinants.py
@@ -16,9 +16,13 @@ from __future__ import absolute_import
 
 import numpy
 from scipy.linalg import schur
+from scipy.sparse import csc_matrix, eye
 
 from openfermion.config import EQ_TOLERANCE
 from openfermion.ops import QuadraticHamiltonian
+from openfermion.utils._sparse_tools import (jw_hartree_fock_state,
+                                             kronecker_operators,
+                                             pauli_matrix_map)
 
 
 def ground_state_preparation_circuit(quadratic_hamiltonian):
@@ -38,6 +42,11 @@ def ground_state_preparation_circuit(quadratic_hamiltonian):
             transformation on the last fermionic mode, or a tuple of
             the form (i, j, theta, phi), indicating a Givens rotation
             of modes i and j by angles theta and phi.
+        n_electrons(int):
+            The number of electrons to start with. This describes the
+            initial state that the circuit should be applied to: it should
+            be a Slater determinant (in the computational basis) with the
+            first n_electrons orbitals filled.
     """
     if not isinstance(quadratic_hamiltonian, QuadraticHamiltonian):
         raise ValueError('Input must be an instance of QuadraticHamiltonian.')
@@ -59,7 +68,8 @@ def ground_state_preparation_circuit(quadratic_hamiltonian):
         # Get the circuit description
         left_unitary, decomposition, diagonal = givens_decomposition(
                 slater_determinant_matrix)
-        circuit_description = decomposition
+        circuit_description = list(reversed(decomposition))
+        n_electrons = num_negative_energies
     else:
         # The Hamiltonian does not conserve particle number, so we
         # need to use the most general procedure.
@@ -75,8 +85,62 @@ def ground_state_preparation_circuit(quadratic_hamiltonian):
         left_unitary, decomposition, diagonal = (
                 fermionic_gaussian_decomposition(
                     slater_determinant_matrix))
-        circuit_description = decomposition
-    return circuit_description
+        circuit_description = list(reversed(decomposition))
+        n_electrons = 0
+    return circuit_description, n_electrons
+
+
+def jw_get_quadratic_hamiltonian_ground_state(quadratic_hamiltonian):
+    """Compute the lowest eigenvalue and eigenstate of a quadratic Hamiltonian.
+
+    Returns:
+        ground_energy(float): The lowest eigenvalue.
+        state(sparse): The lowest eigenstate in scipy.sparse csc format.
+    """
+    if not isinstance(quadratic_hamiltonian, QuadraticHamiltonian):
+        raise ValueError('Input must be an instance of QuadraticHamiltonian.')
+
+    n_qubits = quadratic_hamiltonian.n_qubits
+
+    # Compute the ground energy
+    if quadratic_hamiltonian.conserves_particle_number():
+        hermitian_matrix = quadratic_hamiltonian.combined_hermitian_part()
+        energies, diagonalizing_unitary = numpy.linalg.eigh(
+                hermitian_matrix)
+        num_negative_energies = numpy.count_nonzero(
+                energies < -EQ_TOLERANCE)
+        negative_energies = energies[:num_negative_energies]
+        ground_energy = (numpy.sum(negative_energies) +
+                         quadratic_hamiltonian.constant)
+    else:
+        majorana_matrix, majorana_constant = (
+                quadratic_hamiltonian.majorana_form())
+        canonical, orthogonal = antisymmetric_canonical_form(majorana_matrix)
+        negative_diagonal_entries = canonical[range(n_qubits, 2 * n_qubits),
+                                              range(n_qubits)]
+        ground_energy = (.5 * numpy.sum(negative_diagonal_entries) +
+                         majorana_constant)
+
+    # Obtain the circuit that prepares the ground state
+    circuit_description, n_electrons = ground_state_preparation_circuit(
+            quadratic_hamiltonian)
+
+    # Initialize the starting state
+    state = jw_hartree_fock_state(n_electrons, n_qubits)
+
+    # Apply the circuit
+    particle_hole_transformation = (
+            jw_sparse_particle_hole_transformation_last_mode(n_qubits))
+    for parallel_ops in circuit_description:
+        for op in parallel_ops:
+            if op == 'pht':
+                state = particle_hole_transformation.dot(state)
+            else:
+                i, j, theta, phi = op
+                state = jw_sparse_givens_rotation(
+                            i, j, theta, phi, n_qubits).dot(state)
+
+    return ground_energy, state
 
 
 def fermionic_gaussian_decomposition(unitary_rows):
@@ -149,12 +213,13 @@ def fermionic_gaussian_decomposition(unitary_rows):
     for k in range(n - 1):
         # Zero out entries in column 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])
-            # Apply Givens rotation
-            givens_rotate(current_matrix, givens_rotation, l, l + 1)
-            givens_rotate(left_unitary, givens_rotation, l, l + 1)
+            # Zero out entry in row l if needed
+            if abs(current_matrix[l, k]) > EQ_TOLERANCE:
+                givens_rotation = givens_matrix_elements(
+                        current_matrix[l, k], current_matrix[l + 1, k])
+                # Apply Givens rotation
+                givens_rotate(current_matrix, givens_rotation, l, l + 1)
+                givens_rotate(left_unitary, givens_rotation, l, l + 1)
 
     # Initialize list to store decomposition of current_matrix
     decomposition = list()
@@ -191,7 +256,7 @@ def fermionic_gaussian_decomposition(unitary_rows):
                                                          right_element)
 
                 # Add the parameters to the list
-                theta = numpy.arccos(numpy.real(givens_rotation[0, 0]))
+                theta = numpy.arcsin(numpy.real(givens_rotation[1, 0]))
                 phi = numpy.angle(givens_rotation[1, 1])
                 parallel_ops.append((j, j + 1, theta, phi))
 
@@ -257,12 +322,13 @@ def givens_decomposition(unitary_rows):
     for k in reversed(range(n - m + 1, n)):
         # Zero out entries in column k
         for l in range(m - n + k):
-            # Zero out entry in row l
-            givens_rotation = givens_matrix_elements(current_matrix[l, k],
-                                                     current_matrix[l + 1, k])
-            # Apply Givens rotation
-            givens_rotate(current_matrix, givens_rotation, l, l + 1)
-            givens_rotate(left_unitary, givens_rotation, l, l + 1)
+            # Zero out entry in row l if needed
+            if abs(current_matrix[l, k]) > EQ_TOLERANCE:
+                givens_rotation = givens_matrix_elements(
+                        current_matrix[l, k], current_matrix[l + 1, k])
+                # Apply Givens rotation
+                givens_rotate(current_matrix, givens_rotation, l, l + 1)
+                givens_rotate(left_unitary, givens_rotation, l, l + 1)
 
     # Compute the decomposition of current_matrix into Givens rotations
     givens_rotations = list()
@@ -313,14 +379,11 @@ def givens_decomposition(unitary_rows):
                 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), :]
+                    givens_rotation = givens_matrix_elements(
+                            left_element, right_element, which='right')
 
                     # Add the parameters to the list
-                    theta = numpy.arccos(numpy.real(givens_rotation[0, 0]))
+                    theta = numpy.arcsin(numpy.real(givens_rotation[1, 0]))
                     phi = numpy.angle(givens_rotation[1, 1])
                     parallel_rotations.append((j - 1, j, theta, phi))
 
@@ -440,38 +503,45 @@ def antisymmetric_canonical_form(antisymmetric_matrix):
     return canonical, orthogonal.T
 
 
-def givens_matrix_elements(a, b):
+def givens_matrix_elements(a, b, which='left'):
     """Compute the matrix elements of the Givens rotation that zeroes out one
     of two row entries.
 
-    Returns a matrix G such that
+    If `which='left'` then returns a matrix G such that
 
         G * [a  b]^T= [0  r]^T
 
+    otherwise, returns a matrix G such that
+
+        G * [a  b]^T= [r  0]^T
+
     where r is a complex number.
 
     Args:
-        a: A complex number representing the upper row entry
-        b: A complex number representing the lower row entry
+        a(complex or float): A complex number representing the upper row entry
+        b(complex or float): A complex number representing the lower row entry
+        which(string): Either 'left' or 'right', indicating whether to
+            zero out the left element (first argument) or right element
+            (second argument). Default is `left`.
     Returns:
-        G: A 2 x 2 numpy array representing the matrix G. The numbers in the
-            first column of G are real.
+        G(ndarray): A 2 x 2 numpy array representing the matrix G.
+            The numbers in the first column of G are real.
     """
     # Handle case that a is zero
     if abs(a) < EQ_TOLERANCE:
-        c = 1.
-        s = 0.
+        cosine = 1.
+        sine = 0.
         phase = 1.
     # Handle case that b is zero and a is nonzero
     elif abs(b) < EQ_TOLERANCE:
-        c = 0.
-        s = 1.
+        cosine = 0.
+        sine = 1.
         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
+        cosine = abs(b) / denominator
+        sine = abs(a) / denominator
         sign_b = b / abs(b)
         sign_a = a / abs(a)
         phase = sign_a * sign_b.conjugate()
@@ -480,8 +550,24 @@ def givens_matrix_elements(a, b):
             phase = numpy.real(phase)
 
     # Construct matrix and return
-    givens_rotation = numpy.array([[c, -phase * s],
-                                  [s, phase * c]])
+    if which == 'left':
+        # We want to zero out a
+        if numpy.isreal(a) and numpy.isreal(b):
+            givens_rotation = numpy.array([[cosine, -phase * sine],
+                                           [phase * sine, cosine]])
+        else:
+            givens_rotation = numpy.array([[cosine, -phase * sine],
+                                           [sine, phase * cosine]])
+    elif which == 'right':
+        # We want to zero out b
+        if numpy.isreal(a) and numpy.isreal(b):
+            givens_rotation = numpy.array([[sine, phase * cosine],
+                                           [-phase * cosine, sine]])
+        else:
+            givens_rotation = numpy.array([[sine, phase * cosine],
+                                           [cosine, -phase * sine]])
+    else:
+        raise ValueError('"which" must be equal to "left" or "right".')
     return givens_rotation
 
 
@@ -495,7 +581,7 @@ def givens_rotate(operator, givens_rotation, i, j, which='row'):
                        givens_rotation[0, 1] * row_j)
         operator[j] = (givens_rotation[1, 0] * row_i +
                        givens_rotation[1, 1] * row_j)
-    else:
+    elif which == 'col':
         # Rotate columns i and j
         col_i = operator[:, i].copy()
         col_j = operator[:, j].copy()
@@ -503,6 +589,8 @@ def givens_rotate(operator, givens_rotation, i, j, which='row'):
                           givens_rotation[0, 1].conj() * col_j)
         operator[:, j] = (givens_rotation[1, 0] * col_i +
                           givens_rotation[1, 1].conj() * col_j)
+    else:
+        raise ValueError('"which" must be equal to "row" or "col".')
 
 
 def double_givens_rotate(operator, givens_rotation, i, j, which='row'):
@@ -523,7 +611,7 @@ def double_givens_rotate(operator, givens_rotation, i, j, which='row'):
         givens_rotate(operator[:n], givens_rotation, i, j, which='row')
         # Rotate rows n + i and n + j
         givens_rotate(operator[n:], givens_rotation.conj(), i, j, which='row')
-    else:
+    elif which == 'col':
         if p % 2 != 0:
             raise ValueError('To apply a double Givens rotation on columns, '
                              'the number of columns must be even.')
@@ -533,6 +621,45 @@ def double_givens_rotate(operator, givens_rotation, i, j, which='row'):
         # Rotate cols n + i and n + j
         givens_rotate(operator[:, n:], givens_rotation.conj(), i, j,
                       which='col')
+    else:
+        raise ValueError('"which" must be equal to "row" or "col".')
+
+
+def jw_sparse_givens_rotation(i, j, theta, phi, n_qubits):
+    """Return the matrix (acting on a full wavefunction) that performs a
+    Givens rotation of modes i and j in the Jordan-Wigner encoding."""
+    if j != i + 1:
+        raise ValueError('Only adjacent modes can be rotated.')
+    if j > n_qubits - 1:
+        raise ValueError('Too few qubits requested.')
+
+    cosine = numpy.cos(theta)
+    sine = numpy.sin(theta)
+    phase = numpy.exp(1.j * phi)
+
+    # Create the two-qubit rotation matrix
+    rotation_matrix = csc_matrix(
+            ([1., phase * cosine, -phase * sine, sine, cosine, phase],
+             ((0, 1, 1, 2, 2, 3), (0, 1, 2, 1, 2, 3))),
+            shape=(4, 4))
+
+    # Initialize identity operators
+    left_eye = eye(2 ** i, format='csc')
+    right_eye = eye(2 ** (n_qubits - 1 - j), format='csc')
+
+    # Construct the matrix and return
+    givens_matrix = kronecker_operators([left_eye, rotation_matrix, right_eye])
+
+    return givens_matrix
+
+
+def jw_sparse_particle_hole_transformation_last_mode(n_qubits):
+    """Return the matrix (acting on a full wavefunction) that performs a
+    particle-hole transformation on the last mode in the Jordan-Wigner
+    encoding.
+    """
+    left_eye = eye(2 ** (n_qubits - 1), format='csc')
+    return kronecker_operators([left_eye, pauli_matrix_map['X']])
 
 
 def swap_rows(M, i, j):
diff --git a/src/openfermion/utils/_sparse_tools.py b/src/openfermion/utils/_sparse_tools.py
index 84943bd..a1d8c2f 100644
--- a/src/openfermion/utils/_sparse_tools.py
+++ b/src/openfermion/utils/_sparse_tools.py
@@ -26,7 +26,8 @@ import warnings
 
 from openfermion.config import *
 from openfermion.ops import (FermionOperator, hermitian_conjugated,
-                             normal_ordered, number_operator, QubitOperator)
+                             normal_ordered, number_operator,
+                             QuadraticHamiltonian, QubitOperator)
 from openfermion.utils import commutator, fourier_transform, Grid
 from openfermion.hamiltonians._jellium import (momentum_vector,
                                                position_vector,
@@ -385,25 +386,35 @@ def is_hermitian(sparse_operator):
     return True
 
 
-def get_ground_state(sparse_operator):
+def get_ground_state(operator):
     """Compute lowest eigenvalue and eigenstate.
 
     Returns:
         eigenvalue: The lowest eigenvalue, a float.
         eigenstate: The lowest eigenstate in scipy.sparse csc format.
     """
-    if not is_hermitian(sparse_operator):
-        raise ValueError('sparse_operator must be Hermitian.')
-
-    values, vectors = scipy.sparse.linalg.eigsh(
-        sparse_operator, 2, which='SA', maxiter=1e7)
+    if isinstance(operator, QuadraticHamiltonian):
+        from openfermion.utils._slater_determinants import (
+                jw_get_quadratic_hamiltonian_ground_state)
+        eigenvalue, eigenstate = (
+                jw_get_quadratic_hamiltonian_ground_state(operator))
+    elif isinstance(operator, scipy.sparse.spmatrix):
+        if not is_hermitian(operator):
+            raise ValueError('operator must be Hermitian.')
+
+        values, vectors = scipy.sparse.linalg.eigsh(
+            operator, 2, which='SA', maxiter=1e7)
+
+        order = numpy.argsort(values)
+        values = values[order]
+        vectors = vectors[:, order]
+        eigenvalue = values[0]
+        eigenstate = scipy.sparse.csc_matrix(vectors[:, 0]).T
+    else:
+        raise ValueError('operator must be a sparse matrix or '
+                         'QuadraticHamiltonian.')
 
-    order = numpy.argsort(values)
-    values = values[order]
-    vectors = vectors[:, order]
-    eigenvalue = values[0]
-    eigenstate = scipy.sparse.csc_matrix(vectors[:, 0])
-    return eigenvalue, eigenstate.getH()
+    return eigenvalue, eigenstate
 
 
 def sparse_eigenspectrum(sparse_operator):


Test patch:
diff --git a/src/openfermion/utils/_slater_determinants_test.py b/src/openfermion/utils/_slater_determinants_test.py
index da94ad1..d7f9444 100644
--- a/src/openfermion/utils/_slater_determinants_test.py
+++ b/src/openfermion/utils/_slater_determinants_test.py
@@ -19,58 +19,199 @@ from scipy.linalg import qr
 
 from openfermion.config import EQ_TOLERANCE
 from openfermion.ops import QuadraticHamiltonian
+from openfermion.transforms import get_fermion_operator
 from openfermion.utils import (fermionic_gaussian_decomposition,
+                               get_ground_state,
                                givens_decomposition,
-                               ground_state_preparation_circuit)
+                               ground_state_preparation_circuit,
+                               jordan_wigner_sparse,
+                               jw_get_quadratic_hamiltonian_ground_state)
 from openfermion.utils._slater_determinants import (
         antisymmetric_canonical_form,
         diagonalizing_fermionic_unitary,
         double_givens_rotate,
         givens_rotate,
+        givens_matrix_elements,
         swap_rows)
 
 
 class GroundStatePreparationCircuitTest(unittest.TestCase):
 
+    def test_bad_input(self):
+        """Test bad input."""
+        with self.assertRaises(ValueError):
+            description, n_electrons = ground_state_preparation_circuit('a')
+
+
+class GetQuadraticHamiltonianGroundStateTest(unittest.TestCase):
+
     def setUp(self):
-        self.n_qubits = 5
-        self.constant = 1.7
-        self.chemical_potential = 2.
+        self.n_qubits_range = range(2, 10)
+
+    def test_particle_conserving(self):
+        """Test the case that the Hamiltonian conserves particle number."""
+        constant = 1.7
+        chemical_potential = 2.4
+        for n_qubits in self.n_qubits_range:
+            # Obtain a random Hermitian matrix
+            rand_mat_A = numpy.random.randn(n_qubits, n_qubits)
+            rand_mat_B = numpy.random.randn(n_qubits, n_qubits)
+            rand_mat = rand_mat_A + 1.j * rand_mat_B
+            hermitian_mat = rand_mat + rand_mat.T.conj()
+
+            # Initialize a particle-number-conserving Hamiltonian
+            quadratic_hamiltonian = QuadraticHamiltonian(
+                    constant, hermitian_mat,
+                    chemical_potential=chemical_potential)
+
+            # Compute the true ground state
+            fermion_operator = get_fermion_operator(quadratic_hamiltonian)
+            sparse_operator = jordan_wigner_sparse(fermion_operator)
+            ground_energy, ground_state = get_ground_state(sparse_operator)
+
+            # Compute the ground state using the circuit
+            circuit_energy, circuit_state = (
+                    jw_get_quadratic_hamiltonian_ground_state(
+                        quadratic_hamiltonian))
+
+            # Check that the energies match
+            self.assertAlmostEqual(ground_energy, circuit_energy)
+
+            # Check that the state obtained using the circuit is a ground state
+            difference = (sparse_operator * circuit_state -
+                          ground_energy * circuit_state)
+            discrepancy = 0.
+            if difference.nnz:
+                discrepancy = max(abs(difference.data))
+
+            self.assertTrue(discrepancy < EQ_TOLERANCE)
+
+    def test_particle_nonconserving(self):
+        """Test the case that the Hamiltonian does not conserve particle
+        number."""
+        constant = 1.7
+        chemical_potential = 2.4
+        for n_qubits in self.n_qubits_range:
+            # Obtain random Hermitian and antisymmetric matrices
+            rand_mat_A = numpy.random.randn(n_qubits, n_qubits)
+            rand_mat_B = numpy.random.randn(n_qubits, n_qubits)
+            rand_mat = rand_mat_A + 1.j * rand_mat_B
+            hermitian_mat = rand_mat + rand_mat.T.conj()
+            rand_mat_A = numpy.random.randn(n_qubits, n_qubits)
+            rand_mat_B = numpy.random.randn(n_qubits, n_qubits)
+            rand_mat = rand_mat_A + 1.j * rand_mat_B
+            antisymmetric_mat = rand_mat - rand_mat.T
+
+            # Initialize a non-particle-number-conserving Hamiltonian
+            quadratic_hamiltonian = QuadraticHamiltonian(
+                    constant, hermitian_mat, antisymmetric_mat,
+                    chemical_potential)
+
+            # Compute the true ground state
+            fermion_operator = get_fermion_operator(quadratic_hamiltonian)
+            sparse_operator = jordan_wigner_sparse(fermion_operator)
+            ground_energy, ground_state = get_ground_state(sparse_operator)
+
+            # Compute the ground state using the circuit
+            circuit_energy, circuit_state = (
+                    jw_get_quadratic_hamiltonian_ground_state(
+                        quadratic_hamiltonian))
+
+            # Check that the energies match
+            self.assertAlmostEqual(ground_energy, circuit_energy)
+
+            # Check that the state obtained using the circuit is a ground state
+            difference = (sparse_operator * circuit_state -
+                          ground_energy * circuit_state)
+            discrepancy = 0.
+            if difference.nnz:
+                discrepancy = max(abs(difference.data))
+
+            self.assertTrue(discrepancy < EQ_TOLERANCE)
+
+    def test_bad_input(self):
+        """Test bad input."""
+        with self.assertRaises(ValueError):
+            energy, state = jw_get_quadratic_hamiltonian_ground_state('a')
 
-        # Obtain random Hermitian and antisymmetric matrices
-        rand_mat_A = numpy.random.randn(self.n_qubits, self.n_qubits)
-        rand_mat_B = numpy.random.randn(self.n_qubits, self.n_qubits)
-        rand_mat = rand_mat_A + 1.j * rand_mat_B
-        self.hermitian_mat = rand_mat + rand_mat.T.conj()
-        rand_mat_A = numpy.random.randn(self.n_qubits, self.n_qubits)
-        rand_mat_B = numpy.random.randn(self.n_qubits, self.n_qubits)
-        rand_mat = rand_mat_A + 1.j * rand_mat_B
-        self.antisymmetric_mat = rand_mat - rand_mat.T
 
-        self.combined_hermitian = (
-                self.hermitian_mat -
-                self.chemical_potential * numpy.eye(self.n_qubits))
+class GivensDecompositionTest(unittest.TestCase):
 
-        # Initialize a particle-number-conserving Hamiltonian
-        self.quad_ham_pc = QuadraticHamiltonian(
-                self.constant, self.hermitian_mat)
+    def setUp(self):
+        self.test_dimensions = [(3, 4), (3, 5), (3, 6), (3, 7), (3, 8), (3, 9),
+                                (4, 7), (4, 8), (4, 9)]
+
+    def test_main_procedure(self):
+        for m, n in self.test_dimensions:
+            # Obtain a random matrix of orthonormal rows
+            x = numpy.random.randn(n, n)
+            y = numpy.random.randn(n, n)
+            A = x + 1.j*y
+            Q, R = qr(A)
+            Q = Q[:m, :]
+
+            # Get Givens decomposition of Q
+            V, givens_rotations, diagonal = givens_decomposition(Q)
 
-        # Initialize a non-particle-number-conserving Hamiltonian
-        self.quad_ham_npc = QuadraticHamiltonian(
-                self.constant, self.hermitian_mat, self.antisymmetric_mat,
-                self.chemical_potential)
+            # Compute U
+            U = numpy.eye(n, dtype=complex)
+            for parallel_set in givens_rotations:
+                combined_givens = numpy.eye(n, dtype=complex)
+                for i, j, theta, phi in reversed(parallel_set):
+                    c = numpy.cos(theta)
+                    s = numpy.sin(theta)
+                    phase = numpy.exp(1.j * phi)
+                    G = numpy.array([[c, -phase * s],
+                                    [s, phase * c]], dtype=complex)
+                    givens_rotate(combined_givens, G, i, j)
+                U = combined_givens.dot(U)
+
+            # Compute V * Q * U^\dagger
+            W = V.dot(Q.dot(U.T.conj()))
+
+            # Construct the diagonal matrix
+            D = numpy.zeros((m, n), dtype=complex)
+            D[numpy.diag_indices(m)] = diagonal
+
+            # Assert that W and D are the same
+            for i in range(m):
+                for j in range(n):
+                    self.assertAlmostEqual(D[i, j], W[i, j])
+
+    def test_real_numbers(self):
+        for m, n in self.test_dimensions:
+            # Obtain a random matrix of orthonormal rows
+            A = numpy.random.randn(n, n)
+            Q, R = qr(A)
+            Q = Q[:m, :]
+
+            # Get Givens decomposition of Q
+            V, givens_rotations, diagonal = givens_decomposition(Q)
+
+            # Compute U
+            U = numpy.eye(n, dtype=complex)
+            for parallel_set in givens_rotations:
+                combined_givens = numpy.eye(n, dtype=complex)
+                for i, j, theta, phi in reversed(parallel_set):
+                    c = numpy.cos(theta)
+                    s = numpy.sin(theta)
+                    phase = numpy.exp(1.j * phi)
+                    G = numpy.array([[c, -phase * s],
+                                    [s, phase * c]], dtype=complex)
+                    givens_rotate(combined_givens, G, i, j)
+                U = combined_givens.dot(U)
 
-    def test_no_error(self):
-        """Test that the procedure runs without error."""
-        # Test a particle-number-conserving Hamiltonian
-        circuit_description = (
-                ground_state_preparation_circuit(self.quad_ham_pc))
-        # Test a non-particle-number-conserving Hamiltonian
-        circuit_description = (
-                ground_state_preparation_circuit(self.quad_ham_npc))
+            # Compute V * Q * U^\dagger
+            W = V.dot(Q.dot(U.T.conj()))
 
+            # Construct the diagonal matrix
+            D = numpy.zeros((m, n), dtype=complex)
+            D[numpy.diag_indices(m)] = diagonal
 
-class GivensDecompositionTest(unittest.TestCase):
+            # Assert that W and D are the same
+            for i in range(m):
+                for j in range(n):
+                    self.assertAlmostEqual(D[i, j], W[i, j])
 
     def test_bad_dimensions(self):
         m, n = (3, 2)
@@ -122,7 +263,7 @@ class GivensDecompositionTest(unittest.TestCase):
             for j in range(n):
                 self.assertAlmostEqual(VQ[i, j], D[i, j])
 
-    def test_3_by_3(self):
+    def test_square(self):
         m, n = (3, 3)
         # Obtain a random matrix of orthonormal rows
         x = numpy.random.randn(n, n)
@@ -149,305 +290,65 @@ class GivensDecompositionTest(unittest.TestCase):
             for j in range(n):
                 self.assertAlmostEqual(D[i, j], W[i, j])
 
-    def test_3_by_4(self):
-        m, n = (3, 4)
-        # Obtain a random matrix of orthonormal rows
-        x = numpy.random.randn(n, n)
-        y = numpy.random.randn(n, n)
-        A = x + 1.j*y
-        Q, R = qr(A)
-        Q = Q[:m, :]
-
-        # Get Givens decomposition of Q
-        V, givens_rotations, diagonal = givens_decomposition(Q)
-
-        # Compute U
-        U = numpy.eye(n, dtype=complex)
-        for parallel_set in givens_rotations:
-            combined_givens = numpy.eye(n, dtype=complex)
-            for i, j, theta, phi in parallel_set:
-                c = numpy.cos(theta)
-                s = numpy.sin(theta)
-                phase = numpy.exp(1.j * phi)
-                G = numpy.array([[c, -phase * s],
-                                [s, phase * c]], dtype=complex)
-                givens_rotate(combined_givens, G, i, j)
-            U = combined_givens.dot(U)
-
-        # Compute V * Q * U^\dagger
-        W = V.dot(Q.dot(U.T.conj()))
-
-        # Construct the diagonal matrix
-        D = numpy.zeros((m, n), dtype=complex)
-        D[numpy.diag_indices(m)] = diagonal
-
-        # Assert that W and D are the same
-        for i in range(m):
-            for j in range(n):
-                self.assertAlmostEqual(D[i, j], W[i, j])
-
-    def test_3_by_5(self):
-        m, n = (3, 5)
-        # Obtain a random matrix of orthonormal rows
-        x = numpy.random.randn(n, n)
-        y = numpy.random.randn(n, n)
-        A = x + 1.j*y
-        Q, R = qr(A)
-        Q = Q[:m, :]
-
-        # Get Givens decomposition of Q
-        V, givens_rotations, diagonal = givens_decomposition(Q)
-
-        # Compute U
-        U = numpy.eye(n, dtype=complex)
-        for parallel_set in givens_rotations:
-            combined_givens = numpy.eye(n, dtype=complex)
-            for i, j, theta, phi in parallel_set:
-                c = numpy.cos(theta)
-                s = numpy.sin(theta)
-                phase = numpy.exp(1.j * phi)
-                G = numpy.array([[c, -phase * s],
-                                [s, phase * c]], dtype=complex)
-                givens_rotate(combined_givens, G, i, j)
-            U = combined_givens.dot(U)
-
-        # Compute V * Q * U^\dagger
-        W = V.dot(Q.dot(U.T.conj()))
-
-        # Construct the diagonal matrix
-        D = numpy.zeros((m, n), dtype=complex)
-        D[numpy.diag_indices(m)] = diagonal
-
-        # Assert that W and D are the same
-        for i in range(m):
-            for j in range(n):
-                self.assertAlmostEqual(D[i, j], W[i, j])
-
-    def test_3_by_6(self):
-        m, n = (3, 6)
-        # Obtain a random matrix of orthonormal rows
-        x = numpy.random.randn(n, n)
-        y = numpy.random.randn(n, n)
-        A = x + 1.j*y
-        Q, R = qr(A)
-        Q = Q[:m, :]
-
-        # Get Givens decomposition of Q
-        V, givens_rotations, diagonal = givens_decomposition(Q)
-
-        # Compute U
-        U = numpy.eye(n, dtype=complex)
-        for parallel_set in givens_rotations:
-            combined_givens = numpy.eye(n, dtype=complex)
-            for i, j, theta, phi in parallel_set:
-                c = numpy.cos(theta)
-                s = numpy.sin(theta)
-                phase = numpy.exp(1.j * phi)
-                G = numpy.array([[c, -phase * s],
-                                [s, phase * c]], dtype=complex)
-                givens_rotate(combined_givens, G, i, j)
-            U = combined_givens.dot(U)
-
-        # Compute V * Q * U^\dagger
-        W = V.dot(Q.dot(U.T.conj()))
-
-        # Construct the diagonal matrix
-        D = numpy.zeros((m, n), dtype=complex)
-        D[numpy.diag_indices(m)] = diagonal
-
-        # Assert that W and D are the same
-        for i in range(m):
-            for j in range(n):
-                self.assertAlmostEqual(D[i, j], W[i, j])
-
-    def test_3_by_7(self):
-        m, n = (3, 7)
-        # Obtain a random matrix of orthonormal rows
-        x = numpy.random.randn(n, n)
-        y = numpy.random.randn(n, n)
-        A = x + 1.j*y
-        Q, R = qr(A)
-        Q = Q[:m, :]
-
-        # Get Givens decomposition of Q
-        V, givens_rotations, diagonal = givens_decomposition(Q)
-
-        # Compute U
-        U = numpy.eye(n, dtype=complex)
-        for parallel_set in givens_rotations:
-            combined_givens = numpy.eye(n, dtype=complex)
-            for i, j, theta, phi in parallel_set:
-                c = numpy.cos(theta)
-                s = numpy.sin(theta)
-                phase = numpy.exp(1.j * phi)
-                G = numpy.array([[c, -phase * s],
-                                [s, phase * c]], dtype=complex)
-                givens_rotate(combined_givens, G, i, j)
-            U = combined_givens.dot(U)
-
-        # Compute V * Q * U^\dagger
-        W = V.dot(Q.dot(U.T.conj()))
-
-        # Construct the diagonal matrix
-        D = numpy.zeros((m, n), dtype=complex)
-        D[numpy.diag_indices(m)] = diagonal
-
-        # Assert that W and D are the same
-        for i in range(m):
-            for j in range(n):
-                self.assertAlmostEqual(D[i, j], W[i, j])
-
-    def test_3_by_8(self):
-        m, n = (3, 8)
-        # Obtain a random matrix of orthonormal rows
-        x = numpy.random.randn(n, n)
-        y = numpy.random.randn(n, n)
-        A = x + 1.j*y
-        Q, R = qr(A)
-        Q = Q[:m, :]
-
-        # Get Givens decomposition of Q
-        V, givens_rotations, diagonal = givens_decomposition(Q)
-
-        # Compute U
-        U = numpy.eye(n, dtype=complex)
-        for parallel_set in givens_rotations:
-            combined_givens = numpy.eye(n, dtype=complex)
-            for i, j, theta, phi in parallel_set:
-                c = numpy.cos(theta)
-                s = numpy.sin(theta)
-                phase = numpy.exp(1.j * phi)
-                G = numpy.array([[c, -phase * s],
-                                [s, phase * c]], dtype=complex)
-                givens_rotate(combined_givens, G, i, j)
-            U = combined_givens.dot(U)
-
-        # Compute V * Q * U^\dagger
-        W = V.dot(Q.dot(U.T.conj()))
-
-        # Construct the diagonal matrix
-        D = numpy.zeros((m, n), dtype=complex)
-        D[numpy.diag_indices(m)] = diagonal
-
-        # Assert that W and D are the same
-        for i in range(m):
-            for j in range(n):
-                self.assertAlmostEqual(D[i, j], W[i, j])
-
-    def test_3_by_9(self):
-        m, n = (3, 9)
-        # Obtain a random matrix of orthonormal rows
-        x = numpy.random.randn(n, n)
-        y = numpy.random.randn(n, n)
-        A = x + 1.j*y
-        Q, R = qr(A)
-        Q = Q[:m, :]
-
-        # Get Givens decomposition of Q
-        V, givens_rotations, diagonal = givens_decomposition(Q)
-
-        # Compute U
-        U = numpy.eye(n, dtype=complex)
-        for parallel_set in givens_rotations:
-            combined_givens = numpy.eye(n, dtype=complex)
-            for i, j, theta, phi in parallel_set:
-                c = numpy.cos(theta)
-                s = numpy.sin(theta)
-                phase = numpy.exp(1.j * phi)
-                G = numpy.array([[c, -phase * s],
-                                [s, phase * c]], dtype=complex)
-                givens_rotate(combined_givens, G, i, j)
-            U = combined_givens.dot(U)
-
-        # Compute V * Q * U^\dagger
-        W = V.dot(Q.dot(U.T.conj()))
-
-        # Construct the diagonal matrix
-        D = numpy.zeros((m, n), dtype=complex)
-        D[numpy.diag_indices(m)] = diagonal
-
-        # Assert that W and D are the same
-        for i in range(m):
-            for j in range(n):
-                self.assertAlmostEqual(D[i, j], W[i, j])
-
-    def test_4_by_5(self):
-        m, n = (4, 5)
-        # Obtain a random matrix of orthonormal rows
-        x = numpy.random.randn(n, n)
-        y = numpy.random.randn(n, n)
-        A = x + 1.j*y
-        Q, R = qr(A)
-        Q = Q[:m, :]
-
-        # Get Givens decomposition of Q
-        V, givens_rotations, diagonal = givens_decomposition(Q)
-
-        # Compute U
-        U = numpy.eye(n, dtype=complex)
-        for parallel_set in givens_rotations:
-            combined_givens = numpy.eye(n, dtype=complex)
-            for i, j, theta, phi in parallel_set:
-                c = numpy.cos(theta)
-                s = numpy.sin(theta)
-                phase = numpy.exp(1.j * phi)
-                G = numpy.array([[c, -phase * s],
-                                [s, phase * c]], dtype=complex)
-                givens_rotate(combined_givens, G, i, j)
-            U = combined_givens.dot(U)
-
-        # Compute V * Q * U^\dagger
-        W = V.dot(Q.dot(U.T.conj()))
-
-        # Construct the diagonal matrix
-        D = numpy.zeros((m, n), dtype=complex)
-        D[numpy.diag_indices(m)] = diagonal
-
-        # Assert that W and D are the same
-        for i in range(m):
-            for j in range(n):
-                self.assertAlmostEqual(D[i, j], W[i, j])
-
-    def test_4_by_9(self):
-        m, n = (4, 9)
-        # Obtain a random matrix of orthonormal rows
-        x = numpy.random.randn(n, n)
-        y = numpy.random.randn(n, n)
-        A = x + 1.j*y
-        Q, R = qr(A)
-        Q = Q[:m, :]
-
-        # Get Givens decomposition of Q
-        V, givens_rotations, diagonal = givens_decomposition(Q)
-
-        # Compute U
-        U = numpy.eye(n, dtype=complex)
-        for parallel_set in givens_rotations:
-            combined_givens = numpy.eye(n, dtype=complex)
-            for i, j, theta, phi in parallel_set:
-                c = numpy.cos(theta)
-                s = numpy.sin(theta)
-                phase = numpy.exp(1.j * phi)
-                G = numpy.array([[c, -phase * s],
-                                [s, phase * c]], dtype=complex)
-                givens_rotate(combined_givens, G, i, j)
-            U = combined_givens.dot(U)
-
-        # Compute V * Q * U^\dagger
-        W = V.dot(Q.dot(U.T.conj()))
-
-        # Construct the diagonal matrix
-        D = numpy.zeros((m, n), dtype=complex)
-        D[numpy.diag_indices(m)] = diagonal
-
-        # Assert that W and D are the same
-        for i in range(m):
-            for j in range(n):
-                self.assertAlmostEqual(D[i, j], W[i, j])
-
 
 class FermionicGaussianDecompositionTest(unittest.TestCase):
 
+    def setUp(self):
+        self.test_dimensions = [3, 4, 5, 6, 7, 8, 9]
+
+    def test_main_procedure(self):
+        for n in self.test_dimensions:
+            # Obtain a random antisymmetric matrix
+            rand_mat = numpy.random.randn(2 * n, 2 * n)
+            antisymmetric_matrix = rand_mat - rand_mat.T
+
+            # Get the diagonalizing fermionic unitary
+            ferm_unitary = diagonalizing_fermionic_unitary(
+                    antisymmetric_matrix)
+            lower_unitary = ferm_unitary[n:]
+
+            # Get fermionic Gaussian decomposition of lower_unitary
+            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 reversed(parallel_set):
+                    if op == 'pht':
+                        swap_rows(combined_op, n - 1, 2 * n - 1)
+                    else:
+                        i, j, theta, phi = op
+                        c = numpy.cos(theta)
+                        s = numpy.sin(theta)
+                        phase = numpy.exp(1.j * phi)
+                        givens_rotation = numpy.array(
+                                [[c, -phase * s],
+                                 [s, phase * c]], dtype=complex)
+                        double_givens_rotate(
+                                combined_op, givens_rotation, i, j)
+                right_unitary = combined_op.dot(right_unitary)
+
+            # Compute left_unitary * lower_unitary * right_unitary^\dagger
+            product = left_unitary.dot(lower_unitary.dot(
+                                       right_unitary.T.conj()))
+
+            # 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(diag[i], product[i])
+
     def test_bad_dimensions(self):
         n, p = (3, 7)
         rand_mat = numpy.random.randn(n, p)
@@ -462,358 +363,49 @@ class FermionicGaussianDecompositionTest(unittest.TestCase):
             left_unitary, decomposition, antidiagonal = (
                     fermionic_gaussian_decomposition(ones_mat))
 
-    def test_n_equals_3(self):
-        n = 3
-        # Obtain a random antisymmetric matrix
-        rand_mat = numpy.random.randn(2 * n, 2 * n)
-        antisymmetric_matrix = rand_mat - rand_mat.T
-
-        # Get the diagonalizing fermionic unitary
-        ferm_unitary = diagonalizing_fermionic_unitary(antisymmetric_matrix)
-        lower_unitary = ferm_unitary[n:]
-
-        # Get fermionic Gaussian decomposition of lower_unitary
-        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 == 'pht':
-                    swap_rows(combined_op, n - 1, 2 * n - 1)
-                else:
-                    i, j, theta, phi = op
-                    c = numpy.cos(theta)
-                    s = numpy.sin(theta)
-                    phase = numpy.exp(1.j * phi)
-                    givens_rotation = numpy.array(
-                            [[c, -phase * s],
-                             [s, phase * c]], dtype=complex)
-                    double_givens_rotate(combined_op, givens_rotation, i, j)
-            right_unitary = combined_op.dot(right_unitary)
-
-        # Compute left_unitary * lower_unitary * right_unitary^\dagger
-        product = left_unitary.dot(lower_unitary.dot(right_unitary.T.conj()))
-
-        # 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(diag[i], product[i])
-
-    def test_n_equals_4(self):
-        n = 4
-        # Obtain a random antisymmetric matrix
-        rand_mat = numpy.random.randn(2 * n, 2 * n)
-        antisymmetric_matrix = rand_mat - rand_mat.T
-
-        # Get the diagonalizing fermionic unitary
-        ferm_unitary = diagonalizing_fermionic_unitary(antisymmetric_matrix)
-        lower_unitary = ferm_unitary[n:]
-
-        # Get fermionic Gaussian decomposition of lower_unitary
-        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 == 'pht':
-                    swap_rows(combined_op, n - 1, 2 * n - 1)
-                else:
-                    i, j, theta, phi = op
-                    c = numpy.cos(theta)
-                    s = numpy.sin(theta)
-                    phase = numpy.exp(1.j * phi)
-                    givens_rotation = numpy.array(
-                            [[c, -phase * s],
-                             [s, phase * c]], dtype=complex)
-                    double_givens_rotate(combined_op, givens_rotation, i, j)
-            right_unitary = combined_op.dot(right_unitary)
-
-        # Compute left_unitary * lower_unitary * right_unitary^\dagger
-        product = left_unitary.dot(lower_unitary.dot(right_unitary.T.conj()))
-
-        # 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(diag[i], product[i])
-
-    def test_n_equals_5(self):
-        n = 5
-        # Obtain a random antisymmetric matrix
-        rand_mat = numpy.random.randn(2 * n, 2 * n)
-        antisymmetric_matrix = rand_mat - rand_mat.T
-
-        # Get the diagonalizing fermionic unitary
-        ferm_unitary = diagonalizing_fermionic_unitary(antisymmetric_matrix)
-        lower_unitary = ferm_unitary[n:]
-
-        # Get fermionic Gaussian decomposition of lower_unitary
-        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 == 'pht':
-                    swap_rows(combined_op, n - 1, 2 * n - 1)
-                else:
-                    i, j, theta, phi = op
-                    c = numpy.cos(theta)
-                    s = numpy.sin(theta)
-                    phase = numpy.exp(1.j * phi)
-                    givens_rotation = numpy.array(
-                            [[c, -phase * s],
-                             [s, phase * c]], dtype=complex)
-                    double_givens_rotate(combined_op, givens_rotation, i, j)
-            right_unitary = combined_op.dot(right_unitary)
-
-        # Compute left_unitary * lower_unitary * right_unitary^\dagger
-        product = left_unitary.dot(lower_unitary.dot(right_unitary.T.conj()))
-
-        # 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(diag[i], product[i])
-
-    def test_n_equals_6(self):
-        n = 6
-        # Obtain a random antisymmetric matrix
-        rand_mat = numpy.random.randn(2 * n, 2 * n)
-        antisymmetric_matrix = rand_mat - rand_mat.T
-
-        # Get the diagonalizing fermionic unitary
-        ferm_unitary = diagonalizing_fermionic_unitary(antisymmetric_matrix)
-        lower_unitary = ferm_unitary[n:]
-
-        # Get fermionic Gaussian decomposition of lower_unitary
-        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 == 'pht':
-                    swap_rows(combined_op, n - 1, 2 * n - 1)
-                else:
-                    i, j, theta, phi = op
-                    c = numpy.cos(theta)
-                    s = numpy.sin(theta)
-                    phase = numpy.exp(1.j * phi)
-                    givens_rotation = numpy.array(
-                            [[c, -phase * s],
-                             [s, phase * c]], dtype=complex)
-                    double_givens_rotate(combined_op, givens_rotation, i, j)
-            right_unitary = combined_op.dot(right_unitary)
-
-        # Compute left_unitary * lower_unitary * right_unitary^\dagger
-        product = left_unitary.dot(lower_unitary.dot(right_unitary.T.conj()))
-
-        # 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(diag[i], product[i])
-
-    def test_n_equals_7(self):
-        n = 7
-        # Obtain a random antisymmetric matrix
-        rand_mat = numpy.random.randn(2 * n, 2 * n)
-        antisymmetric_matrix = rand_mat - rand_mat.T
-
-        # Get the diagonalizing fermionic unitary
-        ferm_unitary = diagonalizing_fermionic_unitary(antisymmetric_matrix)
-        lower_unitary = ferm_unitary[n:]
-
-        # Get fermionic Gaussian decomposition of lower_unitary
-        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 == 'pht':
-                    swap_rows(combined_op, n - 1, 2 * n - 1)
-                else:
-                    i, j, theta, phi = op
-                    c = numpy.cos(theta)
-                    s = numpy.sin(theta)
-                    phase = numpy.exp(1.j * phi)
-                    givens_rotation = numpy.array(
-                            [[c, -phase * s],
-                             [s, phase * c]], dtype=complex)
-                    double_givens_rotate(combined_op, givens_rotation, i, j)
-            right_unitary = combined_op.dot(right_unitary)
-
-        # Compute left_unitary * lower_unitary * right_unitary^\dagger
-        product = left_unitary.dot(lower_unitary.dot(right_unitary.T.conj()))
-
-        # 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(diag[i], product[i])
-
-    def test_n_equals_8(self):
-        n = 8
-        # Obtain a random antisymmetric matrix
-        rand_mat = numpy.random.randn(2 * n, 2 * n)
-        antisymmetric_matrix = rand_mat - rand_mat.T
-
-        # Get the diagonalizing fermionic unitary
-        ferm_unitary = diagonalizing_fermionic_unitary(antisymmetric_matrix)
-        lower_unitary = ferm_unitary[n:]
-
-        # Get fermionic Gaussian decomposition of lower_unitary
-        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 == 'pht':
-                    swap_rows(combined_op, n - 1, 2 * n - 1)
-                else:
-                    i, j, theta, phi = op
-                    c = numpy.cos(theta)
-                    s = numpy.sin(theta)
-                    phase = numpy.exp(1.j * phi)
-                    givens_rotation = numpy.array(
-                            [[c, -phase * s],
-                             [s, phase * c]], dtype=complex)
-                    double_givens_rotate(combined_op, givens_rotation, i, j)
-            right_unitary = combined_op.dot(right_unitary)
-
-        # Compute left_unitary * lower_unitary * right_unitary^\dagger
-        product = left_unitary.dot(lower_unitary.dot(right_unitary.T.conj()))
-
-        # 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(diag[i], product[i])
-
-    def test_n_equals_9(self):
-        n = 9
-        # Obtain a random antisymmetric matrix
-        rand_mat = numpy.random.randn(2 * n, 2 * n)
-        antisymmetric_matrix = rand_mat - rand_mat.T
 
-        # Get the diagonalizing fermionic unitary
-        ferm_unitary = diagonalizing_fermionic_unitary(antisymmetric_matrix)
-        lower_unitary = ferm_unitary[n:]
-
-        # Get fermionic Gaussian decomposition of lower_unitary
-        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 == 'pht':
-                    swap_rows(combined_op, n - 1, 2 * n - 1)
-                else:
-                    i, j, theta, phi = op
-                    c = numpy.cos(theta)
-                    s = numpy.sin(theta)
-                    phase = numpy.exp(1.j * phi)
-                    givens_rotation = numpy.array(
-                            [[c, -phase * s],
-                             [s, phase * c]], dtype=complex)
-                    double_givens_rotate(combined_op, givens_rotation, i, j)
-            right_unitary = combined_op.dot(right_unitary)
-
-        # Compute left_unitary * lower_unitary * right_unitary^\dagger
-        product = left_unitary.dot(lower_unitary.dot(right_unitary.T.conj()))
+class DiagonalizingFermionicUnitaryTest(unittest.TestCase):
 
-        # Construct the diagonal matrix
-        diag = numpy.zeros((n, 2 * n), dtype=complex)
-        diag[range(n), range(n, 2 * n)] = diagonal
+    def setUp(self):
+        self.n_qubits = 5
+        self.constant = 1.7
+        self.chemical_potential = 2.
 
-        # Assert that W and D are the same
-        for i in numpy.ndindex((n, 2 * n)):
-            self.assertAlmostEqual(diag[i], product[i])
+        # Obtain random Hermitian and antisymmetric matrices
+        rand_mat_A = numpy.random.randn(self.n_qubits, self.n_qubits)
+        rand_mat_B = numpy.random.randn(self.n_qubits, self.n_qubits)
+        rand_mat = rand_mat_A + 1.j * rand_mat_B
+        self.hermitian_mat = rand_mat + rand_mat.T.conj()
+        rand_mat_A = numpy.random.randn(self.n_qubits, self.n_qubits)
+        rand_mat_B = numpy.random.randn(self.n_qubits, self.n_qubits)
+        rand_mat = rand_mat_A + 1.j * rand_mat_B
+        self.antisymmetric_mat = rand_mat - rand_mat.T
 
+        # Initialize a non-particle-number-conserving Hamiltonian
+        self.quad_ham_npc = QuadraticHamiltonian(
+                self.constant, self.hermitian_mat, self.antisymmetric_mat,
+                self.chemical_potential)
 
-class DiagonalizingFermionicUnitaryTest(unittest.TestCase):
+    def test_diagonalizes_quadratic_hamiltonian(self):
+        """Test that the unitary returned indeed diagonalizes a
+        quadratic Hamiltonian."""
+        hermitian_part = self.quad_ham_npc.combined_hermitian_part()
+        antisymmetric_part = self.quad_ham_npc.antisymmetric_part()
+        block_matrix = numpy.zeros((2 * self.n_qubits, 2 * self.n_qubits),
+                                   dtype=complex)
+        block_matrix[:self.n_qubits, :self.n_qubits] = antisymmetric_part
+        block_matrix[:self.n_qubits, self.n_qubits:] = hermitian_part
+        block_matrix[self.n_qubits:, :self.n_qubits] = -hermitian_part.conj()
+        block_matrix[self.n_qubits:, self.n_qubits:] = (
+                -antisymmetric_part.conj())
+
+        majorana_matrix, majorana_constant = self.quad_ham_npc.majorana_form()
+        canonical, orthogonal = antisymmetric_canonical_form(majorana_matrix)
+        ferm_unitary = diagonalizing_fermionic_unitary(majorana_matrix)
+        diagonalized = ferm_unitary.conj().dot(
+                block_matrix.dot(ferm_unitary.T.conj()))
+        for i in numpy.ndindex((2 * self.n_qubits, 2 * self.n_qubits)):
+            self.assertAlmostEqual(diagonalized[i], canonical[i])
 
     def test_bad_dimensions(self):
         n, p = (3, 4)
@@ -879,6 +471,69 @@ class AntisymmetricCanonicalFormTest(unittest.TestCase):
                 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)
+                    self.assertTrue(canonical[i, j] < EQ_TOLERANCE)
                 else:
                     self.assertAlmostEqual(canonical[i, j], 0.)
+
+
+class GivensMatrixElementsTest(unittest.TestCase):
+
+    def setUp(self):
+        self.num_test_repetitions = 5
+
+    def test_already_zero(self):
+        """Test when some entries are already zero."""
+        # Test when left entry is zero
+        v = numpy.array([0., numpy.random.randn()])
+        G = givens_matrix_elements(v[0], v[1])
+        self.assertAlmostEqual(G.dot(v)[0], 0.)
+        # Test when right entry is zero
+        v = numpy.array([numpy.random.randn(), 0.])
+        G = givens_matrix_elements(v[0], v[1])
+        self.assertAlmostEqual(G.dot(v)[0], 0.)
+
+    def test_real(self):
+        """Test that the procedure works for real numbers."""
+        for _ in range(self.num_test_repetitions):
+            v = numpy.random.randn(2)
+            G_left = givens_matrix_elements(v[0], v[1], which='left')
+            G_right = givens_matrix_elements(v[0], v[1], which='right')
+            self.assertAlmostEqual(G_left.dot(v)[0], 0.)
+            self.assertAlmostEqual(G_right.dot(v)[1], 0.)
+
+    def test_bad_input(self):
+        """Test bad input."""
+        with self.assertRaises(ValueError):
+            v = numpy.random.randn(2)
+            G = givens_matrix_elements(v[0], v[1], which='a')
+
+
+class GivensRotateTest(unittest.TestCase):
+
+    def test_bad_input(self):
+        """Test bad input."""
+        with self.assertRaises(ValueError):
+            v = numpy.random.randn(2)
+            G = givens_matrix_elements(v[0], v[1])
+            givens_rotate(v, G, 0, 1, which='a')
+
+
+class DoubleGivensRotateTest(unittest.TestCase):
+
+    def test_odd_dimension(self):
+        """Test that it raises an error for odd-dimensional input."""
+        A = numpy.random.randn(3, 3)
+        v = numpy.random.randn(2)
+        G = givens_matrix_elements(v[0], v[1])
+        with self.assertRaises(ValueError):
+            double_givens_rotate(A, G, 0, 1, which='row')
+        with self.assertRaises(ValueError):
+            double_givens_rotate(A, G, 0, 1, which='col')
+
+    def test_bad_input(self):
+        """Test bad input."""
+        A = numpy.random.randn(3, 3)
+        v = numpy.random.randn(2)
+        G = givens_matrix_elements(v[0], v[1])
+        with self.assertRaises(ValueError):
+            double_givens_rotate(A, G, 0, 1, which='a')
diff --git a/src/openfermion/utils/_sparse_tools_test.py b/src/openfermion/utils/_sparse_tools_test.py
index 02073f9..da4d3b9 100644
--- a/src/openfermion/utils/_sparse_tools_test.py
+++ b/src/openfermion/utils/_sparse_tools_test.py
@@ -21,7 +21,8 @@ from scipy.sparse import csc_matrix
 
 from openfermion.hamiltonians import jellium_model, wigner_seitz_length_scale
 from openfermion.ops import FermionOperator, normal_ordered, number_operator
-from openfermion.transforms import get_sparse_operator, jordan_wigner
+from openfermion.transforms import (get_fermion_operator, get_sparse_operator,
+                                    jordan_wigner)
 from openfermion.utils import fourier_transform, Grid
 from openfermion.utils._jellium_hf_state import (
     lowest_single_particle_energy_states)
@@ -268,12 +269,43 @@ class GroundStateTest(unittest.TestCase):
 
         self.assertAlmostEqual(ground[0], -2)
         self.assertAlmostEqual(
-            numpy.absolute(expected_state.T.dot(ground[1].A))[0, 0], 1.0)
+            numpy.absolute(
+                expected_state.T.conj().dot(ground[1].A))[0, 0], 1.0)
 
     def test_get_ground_state_nonhermitian(self):
         with self.assertRaises(ValueError):
             get_ground_state(get_sparse_operator(1j * QubitOperator('X1')))
 
+    def test_get_ground_state_quadratic_hamiltonian(self):
+        """Test when input is QuadraticHamiltonian."""
+        n_qubits = 3
+        # Obtain a random Hermitian matrix
+        rand_mat_A = numpy.random.randn(n_qubits, n_qubits)
+        rand_mat_B = numpy.random.randn(n_qubits, n_qubits)
+        rand_mat = rand_mat_A + 1.j * rand_mat_B
+        hermitian_mat = rand_mat + rand_mat.T.conj()
+
+        # Initialize a particle-number-conserving Hamiltonian
+        quadratic_hamiltonian = QuadraticHamiltonian(
+                0., hermitian_mat)
+
+        # Compute the ground energy by passing in a QuadraticHamiltonian
+        energy_1, state_1 = get_ground_state(quadratic_hamiltonian)
+
+        # Compute the ground energy by passing in a sparse operator
+        fermion_operator = get_fermion_operator(quadratic_hamiltonian)
+        sparse_operator = jordan_wigner_sparse(fermion_operator)
+        energy_2, state_2 = get_ground_state(sparse_operator)
+
+        # Check that the energies match
+        self.assertAlmostEqual(energy_1, energy_2)
+
+    def test_bad_input(self):
+        """Test bad input."""
+        H = numpy.eye(3)
+        with self.assertRaises(ValueError):
+            get_ground_state(H)
+
 
 class ExpectationTest(unittest.TestCase):
     def test_expectation_correct(self):