Repository: quantumlib/OpenFermion
Pull request #107 (2017-11-13 00:16:27)

Preparation of arbitrary Gaussian states
- Generalized Slater determinant preparation to arbitrary fermionic Gaussian states. Previously you could only prepare ground states, but now you can specify which orbitals to be filled.
- Added some methods to QuadraticHamiltonian.
- Minor optimization in conversion from FermionOperator to InteractionOperator or QuadraticHamiltonian (avoids having to call `compress()` before conversion sometimes).
- Moved some code from `utils/_slater_determinants.py` to `ops/_quadratic_hamiltonian.py`.


Patch:
diff --git a/examples/givens_rotations.ipynb b/examples/givens_rotations.ipynb
index 24eea98..0a14ff3 100644
--- a/examples/givens_rotations.ipynb
+++ b/examples/givens_rotations.ipynb
@@ -102,7 +102,7 @@
    ],
    "source": [
     "# Get Givens decomposition of U\n",
-    "V, givens_rotations, diagonal = givens_decomposition(Q)\n",
+    "givens_rotations, V, diagonal = givens_decomposition(Q)\n",
     "print(V)\n",
     "print()\n",
     "print(diagonal)"
diff --git a/src/openfermion/ops/_quadratic_hamiltonian.py b/src/openfermion/ops/_quadratic_hamiltonian.py
index 24935fc..a3c97f9 100644
--- a/src/openfermion/ops/_quadratic_hamiltonian.py
+++ b/src/openfermion/ops/_quadratic_hamiltonian.py
@@ -15,6 +15,7 @@ in the fermionic ladder operators."""
 from __future__ import absolute_import
 
 import numpy
+from scipy.linalg import schur
 
 from openfermion.config import EQ_TOLERANCE
 from openfermion.ops import FermionOperator, PolynomialTensor
@@ -114,6 +115,45 @@ class QuadraticHamiltonian(PolynomialTensor):
         discrepancy = numpy.max(numpy.abs(self.antisymmetric_part))
         return discrepancy < EQ_TOLERANCE
 
+    def add_chemical_potential(self, chemical_potential):
+        """Add a chemical potential."""
+        self.n_body_tensors[1, 0] -= (chemical_potential *
+                                      numpy.eye(self.n_qubits))
+        self.chemical_potential += chemical_potential
+
+    def orbital_energies(self):
+        """Return the orbital energies.
+
+        Any quadratic Hamiltonian is unitarily equivalent to a Hamiltonian
+        of the form::
+
+            \sum_{j} \epsilon_j a^\dagger_j a_j + constant.
+
+        We call the \epsilon_j the orbital energies.
+        The eigenvalues of the Hamiltonian are sums of subsets of the
+        orbital energies (up to the additive constant).
+
+        Returns:
+            orbital_energies(ndarray): A one-dimensional array containing
+                the \epsilon_j
+            constant(float): the constant
+        """
+        if self.conserves_particle_number:
+            hermitian_matrix = self.combined_hermitian_part
+            orbital_energies, diagonalizing_unitary = numpy.linalg.eigh(
+                    hermitian_matrix)
+            constant = self.constant
+        else:
+            majorana_matrix, majorana_constant = self.majorana_form()
+            canonical, orthogonal = antisymmetric_canonical_form(
+                    majorana_matrix)
+            orbital_energies = canonical[
+                    range(self.n_qubits),
+                    range(self.n_qubits, 2 * self.n_qubits)]
+            constant = -.5 * numpy.sum(orbital_energies) + majorana_constant
+
+        return orbital_energies, constant
+
     def majorana_form(self):
         """Return the Majorana represention of the Hamiltonian.
 
@@ -201,3 +241,145 @@ def majorana_operator(term=None, coefficient=1.):
     # Invalid input.
     else:
         raise ValueError('Operator specified incorrectly.')
+
+
+def diagonalizing_fermionic_unitary(antisymmetric_matrix):
+    """Compute the unitary that diagonalizes a quadratic Hamiltonian.
+
+    The input matrix represents a quadratic Hamiltonian in the Majorana basis.
+    The output matrix is a unitary that represents a transformation (mixing)
+    of the fermionic ladder operators. We use the convention that the
+    creation operators are listed before the annihilation operators.
+    The returned unitary has additional structure which ensures
+    that the transformed ladder operators also satisfy the fermionic
+    anticommutation relations.
+
+    Args:
+        antisymmetric_matrix(ndarray): A (2 * n_qubits) x (2 * n_qubits)
+            antisymmetric matrix representing a quadratic Hamiltonian in the
+            Majorana basis.
+    Returns:
+        diagonalizing_unitary(ndarray): A (2 * n_qubits) x (2 * n_qubits)
+            unitary matrix representing a transformation of the fermionic
+            ladder operators.
+    """
+    m, n = antisymmetric_matrix.shape
+    n_qubits = n // 2
+
+    # Get the orthogonal transformation that puts antisymmetric_matrix
+    # into canonical form
+    canonical, orthogonal = antisymmetric_canonical_form(antisymmetric_matrix)
+
+    # Create the matrix that converts between fermionic ladder and
+    # Majorana bases
+    normalized_identity = numpy.eye(n_qubits, dtype=complex) / numpy.sqrt(2.)
+    majorana_basis_change = numpy.eye(
+            2 * n_qubits, dtype=complex) / numpy.sqrt(2.)
+    majorana_basis_change[n_qubits:, n_qubits:] *= -1.j
+    majorana_basis_change[:n_qubits, n_qubits:] = normalized_identity
+    majorana_basis_change[n_qubits:, :n_qubits] = 1.j * normalized_identity
+
+    # Compute the unitary and return
+    diagonalizing_unitary = majorana_basis_change.T.conj().dot(
+            orthogonal.dot(majorana_basis_change))
+
+    return diagonalizing_unitary
+
+
+def antisymmetric_canonical_form(antisymmetric_matrix):
+    """Compute the canonical form of an antisymmetric matrix.
+
+    The input is a real, antisymmetric n x n matrix A, where n is even.
+    Its canonical form is::
+
+        A = R^T C R
+
+    where R is a real, orthogonal matrix and C has the form::
+
+        [  0     D ]
+        [ -D     0 ]
+
+    where D is a diagonal matrix with nonnegative entries.
+
+    Args:
+        antisymmetric_matrix(ndarray): An antisymmetric matrix with even
+            dimension.
+
+    Returns:
+        canonical(ndarray): The canonical form C of antisymmetric_matrix
+        orthogonal(ndarray): The orthogonal transformation R.
+    """
+    m, p = antisymmetric_matrix.shape
+
+    if m != p or p % 2 != 0:
+        raise ValueError('The input matrix must be square with even '
+                         'dimension.')
+
+    # Check that input matrix is antisymmetric
+    matrix_plus_transpose = antisymmetric_matrix + antisymmetric_matrix.T
+    maxval = numpy.max(numpy.abs(matrix_plus_transpose))
+    if maxval > EQ_TOLERANCE:
+        raise ValueError('The input matrix must be antisymmetric.')
+
+    # Compute Schur decomposition
+    canonical, orthogonal = schur(antisymmetric_matrix, output='real')
+
+    # The returned form is block diagonal; we need to permute rows and columns
+    # to put it into the form we want
+    n = p // 2
+    for i in range(1, n, 2):
+        swap_rows(canonical, i, n + i - 1)
+        swap_columns(canonical, i, n + i - 1)
+        swap_columns(orthogonal, i, n + i - 1)
+        if n % 2 != 0:
+            swap_rows(canonical, n - 1, n + i)
+            swap_columns(canonical, n - 1, n + i)
+            swap_columns(orthogonal, n - 1, n + i)
+
+    # Now we permute so that the upper right block is non-negative
+    for i in range(n):
+        if canonical[i, n + i] < -EQ_TOLERANCE:
+            swap_rows(canonical, i, n + i)
+            swap_columns(canonical, i, n + i)
+            swap_columns(orthogonal, i, n + i)
+
+    # Now we permute so that the nonzero entries are ordered by magnitude
+    # We use insertion sort
+    diagonal = canonical[range(n), range(n, 2 * n)]
+    for i in range(n):
+        # Insert the smallest element from the unsorted part of the list into
+        # index i
+        arg_min = numpy.argmin(diagonal[i:]) + i
+        if arg_min != i:
+            # Permute the upper right block
+            swap_rows(canonical, i, arg_min)
+            swap_columns(canonical, n + i, n + arg_min)
+            swap_columns(orthogonal, n + i, n + arg_min)
+            # Permute the lower left block
+            swap_rows(canonical, n + i, n + arg_min)
+            swap_columns(canonical, i, arg_min)
+            swap_columns(orthogonal, i, arg_min)
+            # Update diagonal
+            swap_rows(diagonal, i, arg_min)
+
+    return canonical, orthogonal.T
+
+
+def swap_rows(M, i, j):
+    """Swap rows i and j of matrix M."""
+    if len(M.shape) == 1:
+        M[i], M[j] = M[j], M[i]
+    else:
+        row_i = M[i, :].copy()
+        row_j = M[j, :].copy()
+        M[i, :], M[j, :] = row_j, row_i
+
+
+def swap_columns(M, i, j):
+    """Swap columns i and j of matrix M."""
+    if len(M.shape) == 1:
+        M[i], M[j] = M[j], M[i]
+    else:
+        column_i = M[:, i].copy()
+        column_j = M[:, j].copy()
+        M[:, i], M[:, j] = column_j, column_i
diff --git a/src/openfermion/transforms/_conversion.py b/src/openfermion/transforms/_conversion.py
index 7beba8d..f51badf 100644
--- a/src/openfermion/transforms/_conversion.py
+++ b/src/openfermion/transforms/_conversion.py
@@ -120,6 +120,9 @@ def get_interaction_operator(fermion_operator, n_qubits=None):
     # Loop through terms and assign to matrix.
     for term in fermion_operator.terms:
         coefficient = fermion_operator.terms[term]
+        # Ignore this term if the coefficient is zero
+        if abs(coefficient) < EQ_TOLERANCE:
+            continue
 
         # Handle constant shift.
         if len(term) == 0:
@@ -190,6 +193,9 @@ def get_quadratic_hamiltonian(fermion_operator,
     # Loop through terms and assign to matrix.
     for term in fermion_operator.terms:
         coefficient = fermion_operator.terms[term]
+        # Ignore this term if the coefficient is zero
+        if abs(coefficient) < EQ_TOLERANCE:
+            continue
 
         if len(term) == 0:
             # Constant term
@@ -238,7 +244,7 @@ def get_quadratic_hamiltonian(fermion_operator,
         else:
             # Operator contains non-quadratic terms
             raise QuadraticHamiltonianError('FermionOperator does not map '
-                                            'to QuadraticHamiltonian'
+                                            'to QuadraticHamiltonian '
                                             '(contains non-quadratic terms).')
 
     # Compute Hermitian part
diff --git a/src/openfermion/utils/__init__.py b/src/openfermion/utils/__init__.py
index bb46d80..7092aa6 100644
--- a/src/openfermion/utils/__init__.py
+++ b/src/openfermion/utils/__init__.py
@@ -22,8 +22,8 @@ from ._operator_utils import (commutator, count_qubits,
 
 from ._slater_determinants import (fermionic_gaussian_decomposition,
                                    givens_decomposition,
-                                   ground_state_preparation_circuit,
-                                   jw_get_quadratic_hamiltonian_ground_state)
+                                   gaussian_state_preparation_circuit,
+                                   jw_get_gaussian_state)
 
 from ._sparse_tools import (expectation,
                             expectation_computational_basis_state,
diff --git a/src/openfermion/utils/_slater_determinants.py b/src/openfermion/utils/_slater_determinants.py
index 9d0dd1a..c91200b 100644
--- a/src/openfermion/utils/_slater_determinants.py
+++ b/src/openfermion/utils/_slater_determinants.py
@@ -15,23 +15,35 @@ Slater determinants and fermionic Gaussian states."""
 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,
+from openfermion.ops._quadratic_hamiltonian import (
+        antisymmetric_canonical_form, diagonalizing_fermionic_unitary,
+        swap_columns)
+from openfermion.utils._sparse_tools import (jw_slater_determinant,
                                              kronecker_operators,
                                              pauli_matrix_map)
 
 
-def ground_state_preparation_circuit(quadratic_hamiltonian):
-    """Obtain a description of a circuit which prepares the ground state
-    of a quadratic Hamiltonian.
+def gaussian_state_preparation_circuit(
+        quadratic_hamiltonian, occupied_orbitals=None):
+    """Obtain a description of a circuit which prepares a fermionic Gaussian
+    state.
+
+    Fermionic Gaussian states can be regarded as eigenstates of quadratic
+    Hamiltonians. If the Hamiltonian conserves particle number, then these are
+    just Slater determinants.
 
     Args:
         quadratic_hamiltonian(QuadraticHamiltonian):
-            The Hamiltonian whose ground state is desired.
+            The Hamiltonian whose eigenstate is desired.
+        occupied_orbitals(list):
+            A list of integers representing the indices of the occupied
+            orbitals in the desired Gaussian state. If this is None
+            (the default), then it is assumed that the ground state is
+            desired, i.e., the orbitals with negative energies are filled.
 
     Returns:
         circuit_description(list[tuple]):
@@ -42,11 +54,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
+        start_orbitals(list):
+            The occupied orbitals 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.
+            be a Slater determinant (in the computational basis) with these
+            orbitals filled.
     """
     if not isinstance(quadratic_hamiltonian, QuadraticHamiltonian):
         raise ValueError('Input must be an instance of QuadraticHamiltonian.')
@@ -55,78 +67,95 @@ def ground_state_preparation_circuit(quadratic_hamiltonian):
         # The Hamiltonian conserves particle number, so we don't need
         # to use the most general procedure.
         hermitian_matrix = quadratic_hamiltonian.combined_hermitian_part
-        # Get the unitary rows which represent the ground state
-        # Slater determinant
-        energies, diagonalizing_unitary = numpy.linalg.eigh(
-                hermitian_matrix)
-        # We get the ground state by filling the orbitals that have
-        # negative energy
-        num_negative_energies = numpy.count_nonzero(
-                energies < -EQ_TOLERANCE)
-        slater_determinant_matrix = diagonalizing_unitary.T[
-                :num_negative_energies]
+        energies, diagonalizing_unitary = numpy.linalg.eigh(hermitian_matrix)
+
+        if occupied_orbitals is None:
+            # The ground state is desired, so we fill the orbitals that have
+            # negative energy
+            num_negative_energies = numpy.count_nonzero(
+                    energies < -EQ_TOLERANCE)
+            occupied_orbitals = range(num_negative_energies)
+
+        # Get the unitary rows which represent the Slater determinant
+        slater_determinant_matrix = diagonalizing_unitary.T[occupied_orbitals]
+
         # Get the circuit description
-        left_unitary, decomposition, diagonal = givens_decomposition(
+        decomposition, left_unitary, diagonal = givens_decomposition(
                 slater_determinant_matrix)
         circuit_description = list(reversed(decomposition))
-        n_electrons = num_negative_energies
+        start_orbitals = range(len(occupied_orbitals))
     else:
         # The Hamiltonian does not conserve particle number, so we
         # need to use the most general procedure.
         majorana_matrix, majorana_constant = (
                 quadratic_hamiltonian.majorana_form())
-        # Get the unitary rows which represent the ground state
-        # Slater determinant
         diagonalizing_unitary = diagonalizing_fermionic_unitary(
                 majorana_matrix)
-        slater_determinant_matrix = diagonalizing_unitary[
+
+        # Get the unitary rows which represent the Gaussian unitary
+        gaussian_unitary_matrix = diagonalizing_unitary[
                 quadratic_hamiltonian.n_qubits:]
+
         # Get the circuit description
-        left_unitary, decomposition, diagonal = (
+        decomposition, left_decomposition, diagonal, left_diagonal = (
                 fermionic_gaussian_decomposition(
-                    slater_determinant_matrix))
-        circuit_description = list(reversed(decomposition))
-        n_electrons = 0
-    return circuit_description, n_electrons
+                    gaussian_unitary_matrix))
+        if occupied_orbitals is None:
+            # The ground state is desired, so the circuit should be applied
+            # to the vaccuum state
+            start_orbitals = []
+            circuit_description = list(reversed(decomposition))
+        else:
+            start_orbitals = occupied_orbitals
+            # The circuit won't be applied to the ground state, so we need to
+            # use left_decomposition
+            circuit_description = list(reversed(
+                                       decomposition + left_decomposition))
+
+    return circuit_description, start_orbitals
+
 
+def jw_get_gaussian_state(quadratic_hamiltonian, occupied_orbitals=None):
+    """Compute an eigenvalue and eigenstate of a quadratic Hamiltonian.
 
-def jw_get_quadratic_hamiltonian_ground_state(quadratic_hamiltonian):
-    """Compute the lowest eigenvalue and eigenstate of a quadratic Hamiltonian.
+    Eigenstates of a quadratic Hamiltonian are also known as fermionic
+    Gaussian states.
+
+    Args:
+        quadratic_hamiltonian(QuadraticHamiltonian):
+            The Hamiltonian whose eigenstate is desired.
+        occupied_orbitals(list):
+            A list of integers representing the indices of the occupied
+            orbitals in the desired Gaussian state. If this is None
+            (the default), then it is assumed that the ground state is
+            desired, i.e., the orbitals with negative energies are filled.
 
     Returns:
-        ground_energy(float): The lowest eigenvalue.
-        state(sparse): The lowest eigenstate in scipy.sparse csc format.
+        energy(float): The eigenvalue.
+        state(sparse): The 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)
+    # Compute the energy
+    orbital_energies, constant = quadratic_hamiltonian.orbital_energies()
+    if occupied_orbitals is None:
+        if quadratic_hamiltonian.conserves_particle_number:
+            num_negative_energies = numpy.count_nonzero(
+                    orbital_energies < -EQ_TOLERANCE)
+            occupied_orbitals = range(num_negative_energies)
+        else:
+            occupied_orbitals = []
+    energy = numpy.sum(orbital_energies[occupied_orbitals]) + constant
 
-    # Obtain the circuit that prepares the ground state
-    circuit_description, n_electrons = ground_state_preparation_circuit(
-            quadratic_hamiltonian)
+    # Obtain the circuit that prepares the Gaussian state
+    circuit_description, start_orbitals = gaussian_state_preparation_circuit(
+            quadratic_hamiltonian, occupied_orbitals)
 
     # Initialize the starting state
-    state = jw_hartree_fock_state(n_electrons, n_qubits)
+    state = jw_slater_determinant(start_orbitals, n_qubits)
 
     # Apply the circuit
     particle_hole_transformation = (
@@ -140,7 +169,7 @@ def jw_get_quadratic_hamiltonian_ground_state(quadratic_hamiltonian):
                 state = jw_sparse_givens_rotation(
                             i, j, theta, phi, n_qubits).dot(state)
 
-    return ground_energy, state
+    return energy, state
 
 
 def fermionic_gaussian_decomposition(unitary_rows):
@@ -175,14 +204,20 @@ def fermionic_gaussian_decomposition(unitary_rows):
     of the form (i, j, theta, phi), which indicates a Givens rotation
     of rows i and j by angles theta and phi.
 
+    The matrix V^T D^*, the transpose of V times the complex conjugate of D,
+    can also be decomposed as a sequence of Givens rotations. This
+    decomposition is needed for a circuit that prepares an excited state.
+
     Args:
         unitary_rows(ndarray): A matrix with orthonormal rows and
             additional structure described above.
 
     Returns:
-        left_unitary(ndarray): An n x n matrix representing V.
         decomposition(list[tuple]): The decomposition of U.
+        left_decomposition(list[tuple]): The decomposition of V^T D^*.
         diagonal(ndarray): A list of the nonzero entries of D.
+        left_diagonal(ndarray): A list of the nonzero entries left from
+            the decomposition of V^T D^*.
     """
     current_matrix = numpy.copy(unitary_rows)
     n, p = current_matrix.shape
@@ -222,12 +257,12 @@ def fermionic_gaussian_decomposition(unitary_rows):
                 givens_rotate(left_unitary, givens_rotation, l, l + 1)
 
     # Initialize list to store decomposition of current_matrix
-    decomposition = list()
+    decomposition = []
     # There are 2 * n - 1 iterations (that is the circuit depth)
     for k in range(2 * n - 1):
         # Initialize the list of parallel operations to perform
         # in this iteration
-        parallel_ops = list()
+        parallel_ops = []
 
         # Perform a particle-hole transformation if necessary
         if k % 2 == 0 and abs(current_matrix[k // 2, n - 1]) > EQ_TOLERANCE:
@@ -269,7 +304,56 @@ def fermionic_gaussian_decomposition(unitary_rows):
 
     # Get the diagonal entries
     diagonal = current_matrix[range(n), range(n, 2 * n)]
-    return left_unitary, decomposition, diagonal
+
+    # Compute the decomposition of left_unitary^T * diagonal^*
+    current_matrix = left_unitary.T
+    for k in range(n):
+        current_matrix[:, k] *= diagonal[k].conj()
+    left_decomposition = []
+
+    for k in range(2 * (n - 1) - 1):
+        # Initialize the list of parallel operations to perform
+        # in this iteration
+        parallel_ops = []
+
+        # Get the (row, column) indices of elements to zero out in parallel.
+        if k < n - 1:
+            start_row = 0
+            start_column = n - 1 - k
+        else:
+            start_row = k - (n - 2)
+            start_column = k - (n - 3)
+        column_indices = range(start_column, n, 2)
+        row_indices = range(start_row, start_row + len(column_indices))
+        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,
+            # if needed
+            right_element = current_matrix[i, j].conj()
+            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,
+                                                         which='right')
+
+                # Add the parameters to the list
+                theta = numpy.arcsin(numpy.real(givens_rotation[1, 0]))
+                phi = numpy.angle(givens_rotation[1, 1])
+                parallel_ops.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
+        left_decomposition.append(tuple(parallel_ops))
+
+    # Get the diagonal entries
+    left_diagonal = current_matrix[range(n), range(n)]
+
+    return decomposition, left_decomposition, diagonal, left_diagonal
 
 
 def givens_decomposition(unitary_rows):
@@ -299,15 +383,13 @@ def givens_decomposition(unitary_rows):
             representing the matrix Q.
 
     Returns:
-        left_unitary: An m x m numpy array representing the matrix V.
-
         givens_rotations: A list of tuples of objects describing Givens
             rotations. The list looks like [(G_1, ), (G_2, G_3), ... ].
             The Givens rotations within a tuple can be implemented in parallel.
             The description of a Givens rotation is itself a tuple of the
             form (i, j, theta, phi), which represents a Givens rotation of
             rows i and j by angles theta and phi.
-
+        left_unitary: An m x m numpy array representing the matrix V.
         diagonal: A list of the nonzero entries of D.
     """
     current_matrix = numpy.copy(unitary_rows)
@@ -331,7 +413,7 @@ def givens_decomposition(unitary_rows):
                 givens_rotate(left_unitary, givens_rotation, l, l + 1)
 
     # Compute the decomposition of current_matrix into Givens rotations
-    givens_rotations = list()
+    givens_rotations = []
     # If m = n (the matrix is square) then we don't need to perform any
     # Givens rotations!
     if m != n:
@@ -371,7 +453,7 @@ def givens_decomposition(unitary_rows):
             column_indices = range(start_column, end_column, 2)
             indices_to_zero_out = zip(row_indices, column_indices)
 
-            parallel_rotations = list()
+            parallel_rotations = []
             for i, j in indices_to_zero_out:
                 # Compute the Givens rotation to zero out the (i, j) element,
                 # if needed
@@ -397,110 +479,7 @@ def givens_decomposition(unitary_rows):
     # Get the diagonal entries
     diagonal = current_matrix.diagonal()
 
-    return left_unitary, givens_rotations, diagonal
-
-
-def diagonalizing_fermionic_unitary(antisymmetric_matrix):
-    """Compute the unitary that diagonalizes a quadratic Hamiltonian.
-
-    The input matrix represents a quadratic Hamiltonian in the Majorana basis.
-    The output matrix is a unitary that represents a transformation (mixing)
-    of the fermionic ladder operators. We use the convention that the
-    creation operators are listed before the annihilation operators.
-    The returned unitary has additional structure which ensures
-    that the transformed ladder operators also satisfy the fermionic
-    anticommutation relations.
-
-    Args:
-        antisymmetric_matrix(ndarray): A (2 * n_qubits) x (2 * n_qubits)
-            antisymmetric matrix representing a quadratic Hamiltonian in the
-            Majorana basis.
-    Returns:
-        diagonalizing_unitary(ndarray): A (2 * n_qubits) x (2 * n_qubits)
-            unitary matrix representing a transformation of the fermionic
-            ladder operators.
-    """
-    m, n = antisymmetric_matrix.shape
-    n_qubits = n // 2
-
-    # Get the orthogonal transformation that puts antisymmetric_matrix
-    # into canonical form
-    canonical, orthogonal = antisymmetric_canonical_form(antisymmetric_matrix)
-
-    # Create the matrix that converts between fermionic ladder and
-    # Majorana bases
-    normalized_identity = numpy.eye(n_qubits, dtype=complex) / numpy.sqrt(2.)
-    majorana_basis_change = numpy.eye(
-            2 * n_qubits, dtype=complex) / numpy.sqrt(2.)
-    majorana_basis_change[n_qubits:, n_qubits:] *= -1.j
-    majorana_basis_change[:n_qubits, n_qubits:] = normalized_identity
-    majorana_basis_change[n_qubits:, :n_qubits] = 1.j * normalized_identity
-
-    # Compute the unitary and return
-    diagonalizing_unitary = majorana_basis_change.T.conj().dot(
-            orthogonal.dot(majorana_basis_change))
-
-    return diagonalizing_unitary
-
-
-def antisymmetric_canonical_form(antisymmetric_matrix):
-    """Compute the canonical form of an antisymmetric matrix.
-
-    The input is a real, antisymmetric n x n matrix A, where n is even.
-    Its canonical form is::
-
-        A = R^T C R
-
-    where R is a real, orthogonal matrix and C has the form::
-
-        [  0     D ]
-        [ -D     0 ]
-
-    where D is a diagonal matrix with nonnegative entries.
-
-    Args:
-        antisymmetric_matrix(ndarray): An antisymmetric matrix with even
-            dimension.
-
-    Returns:
-        canonical(ndarray): The canonical form C of antisymmetric_matrix
-        orthogonal(ndarray): The orthogonal transformation R.
-    """
-    m, n = antisymmetric_matrix.shape
-
-    if m != n or n % 2 != 0:
-        raise ValueError('The input matrix must be square with even '
-                         'dimension.')
-
-    # Check that input matrix is antisymmetric
-    matrix_plus_transpose = antisymmetric_matrix + antisymmetric_matrix.T
-    maxval = numpy.max(numpy.abs(matrix_plus_transpose))
-    if maxval > EQ_TOLERANCE:
-        raise ValueError('The input matrix must be antisymmetric.')
-
-    # Compute Schur decomposition
-    canonical, orthogonal = schur(antisymmetric_matrix, output='real')
-
-    # The returned form is block diagonal; we need to permute rows and columns
-    # to put it into the form we want
-    num_blocks = n // 2
-    for i in range(1, num_blocks, 2):
-        swap_rows(canonical, i, num_blocks + i - 1)
-        swap_columns(canonical, i, num_blocks + i - 1)
-        swap_columns(orthogonal, i, num_blocks + i - 1)
-        if num_blocks % 2 != 0:
-            swap_rows(canonical, num_blocks - 1, num_blocks + i)
-            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
+    return givens_rotations, left_unitary, diagonal
 
 
 def givens_matrix_elements(a, b, which='left'):
@@ -660,17 +639,3 @@ def jw_sparse_particle_hole_transformation_last_mode(n_qubits):
     """
     left_eye = eye(2 ** (n_qubits - 1), format='csc')
     return kronecker_operators([left_eye, pauli_matrix_map['X']])
-
-
-def swap_rows(M, i, j):
-    """Swap rows i and j of matrix M."""
-    row_i = M[i, :].copy()
-    row_j = M[j, :].copy()
-    M[i, :], M[j, :] = row_j, row_i
-
-
-def swap_columns(M, i, j):
-    """Swap columns i and j of matrix M."""
-    column_i = M[:, i].copy()
-    column_j = M[:, j].copy()
-    M[:, i], M[:, j] = column_j, column_i
diff --git a/src/openfermion/utils/_sparse_tools.py b/src/openfermion/utils/_sparse_tools.py
index a1d8c2f..55cf9b9 100644
--- a/src/openfermion/utils/_sparse_tools.py
+++ b/src/openfermion/utils/_sparse_tools.py
@@ -395,9 +395,9 @@ def get_ground_state(operator):
     """
     if isinstance(operator, QuadraticHamiltonian):
         from openfermion.utils._slater_determinants import (
-                jw_get_quadratic_hamiltonian_ground_state)
+                jw_get_gaussian_state)
         eigenvalue, eigenstate = (
-                jw_get_quadratic_hamiltonian_ground_state(operator))
+                jw_get_gaussian_state(operator))
     elif isinstance(operator, scipy.sparse.spmatrix):
         if not is_hermitian(operator):
             raise ValueError('operator must be Hermitian.')


Test patch:
diff --git a/src/openfermion/ops/_quadratic_hamiltonian_test.py b/src/openfermion/ops/_quadratic_hamiltonian_test.py
index f34c68a..8c2079b 100644
--- a/src/openfermion/ops/_quadratic_hamiltonian_test.py
+++ b/src/openfermion/ops/_quadratic_hamiltonian_test.py
@@ -21,7 +21,8 @@ from openfermion.ops import (FermionOperator,
                              QuadraticHamiltonian,
                              normal_ordered)
 from openfermion.ops._quadratic_hamiltonian import majorana_operator
-from openfermion.transforms import get_fermion_operator
+from openfermion.transforms import get_fermion_operator, get_sparse_operator
+from openfermion.utils import get_ground_state
 
 
 class QuadraticHamiltoniansTest(unittest.TestCase):
@@ -32,14 +33,8 @@ class QuadraticHamiltoniansTest(unittest.TestCase):
         self.chemical_potential = 2.
 
         # 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.hermitian_mat = random_hermitian_matrix(self.n_qubits)
+        self.antisymmetric_mat = random_antisymmetric_matrix(self.n_qubits)
 
         self.combined_hermitian = (
                 self.hermitian_mat -
@@ -92,6 +87,46 @@ class QuadraticHamiltoniansTest(unittest.TestCase):
         self.assertTrue(self.quad_ham_pc.conserves_particle_number)
         self.assertFalse(self.quad_ham_npc.conserves_particle_number)
 
+    def test_add_chemical_potential(self):
+        """Test adding a chemical potential."""
+        self.quad_ham_npc.add_chemical_potential(2.4)
+
+        combined_hermitian_part = self.quad_ham_npc.combined_hermitian_part
+        hermitian_part = self.quad_ham_npc.hermitian_part
+
+        want_combined = (self.combined_hermitian -
+                         2.4 * numpy.eye(self.n_qubits))
+        want_hermitian = self.hermitian_mat
+
+        for i in numpy.ndindex(combined_hermitian_part.shape):
+            self.assertAlmostEqual(combined_hermitian_part[i],
+                                   want_combined[i])
+
+        for i in numpy.ndindex(hermitian_part.shape):
+            self.assertAlmostEqual(hermitian_part[i], want_hermitian[i])
+
+        self.assertAlmostEqual(2.4 + self.chemical_potential,
+                               self.quad_ham_npc.chemical_potential)
+
+    def test_orbital_energies(self):
+        """Test getting the orbital energies."""
+        # Test the particle-number-conserving case
+        orbital_energies, constant = self.quad_ham_pc.orbital_energies()
+        # Test the ground energy
+        energy = numpy.sum(
+                orbital_energies[orbital_energies < -EQ_TOLERANCE]) + constant
+        quad_ham_pc_sparse = get_sparse_operator(self.quad_ham_pc)
+        ground_energy, ground_state = get_ground_state(quad_ham_pc_sparse)
+        self.assertAlmostEqual(energy, ground_energy)
+
+        # Test the non-particle-number-conserving case
+        orbital_energies, constant = self.quad_ham_npc.orbital_energies()
+        # Test the ground energy
+        energy = constant
+        quad_ham_npc_sparse = get_sparse_operator(self.quad_ham_npc)
+        ground_energy, ground_state = get_ground_state(quad_ham_npc_sparse)
+        self.assertAlmostEqual(energy, ground_energy)
+
     def test_majorana_form(self):
         """Test getting the Majorana form."""
         majorana_matrix, majorana_constant = self.quad_ham_npc.majorana_form()
@@ -131,3 +166,23 @@ class MajoranaOperatorTest(unittest.TestCase):
             majorana_op = majorana_operator((2, 2))
         with self.assertRaises(ValueError):
             majorana_op = majorana_operator('a')
+
+
+def random_hermitian_matrix(n, real=False):
+    """Generate a random n x n Hermitian matrix."""
+    if real:
+        rand_mat = numpy.random.randn(n, n)
+    else:
+        rand_mat = numpy.random.randn(n, n) + 1.j * numpy.random.randn(n, n)
+    hermitian_mat = rand_mat + rand_mat.T.conj()
+    return hermitian_mat
+
+
+def random_antisymmetric_matrix(n, real=False):
+    """Generate a random n x n antisymmetric matrix."""
+    if real:
+        rand_mat = numpy.random.randn(n, n)
+    else:
+        rand_mat = numpy.random.randn(n, n) + 1.j * numpy.random.randn(n, n)
+    antisymmetric_mat = rand_mat - rand_mat.T
+    return antisymmetric_mat
diff --git a/src/openfermion/utils/_slater_determinants_test.py b/src/openfermion/utils/_slater_determinants_test.py
index 29bfa25..812f79e 100644
--- a/src/openfermion/utils/_slater_determinants_test.py
+++ b/src/openfermion/utils/_slater_determinants_test.py
@@ -19,60 +19,136 @@ 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.ops._quadratic_hamiltonian import (
+        antisymmetric_canonical_form,
+        diagonalizing_fermionic_unitary,
+        swap_rows)
+from openfermion.ops._quadratic_hamiltonian_test import (
+        random_hermitian_matrix, random_antisymmetric_matrix)
+from openfermion.transforms import get_sparse_operator
 from openfermion.utils import (fermionic_gaussian_decomposition,
                                get_ground_state,
                                givens_decomposition,
-                               ground_state_preparation_circuit,
-                               jordan_wigner_sparse,
-                               jw_get_quadratic_hamiltonian_ground_state)
+                               gaussian_state_preparation_circuit,
+                               jw_get_gaussian_state,
+                               jw_slater_determinant)
 from openfermion.utils._slater_determinants import (
-        antisymmetric_canonical_form,
-        diagonalizing_fermionic_unitary,
         double_givens_rotate,
         givens_rotate,
         givens_matrix_elements,
-        swap_rows)
+        jw_sparse_givens_rotation,
+        jw_sparse_particle_hole_transformation_last_mode)
 
 
-class GroundStatePreparationCircuitTest(unittest.TestCase):
+class GaussianStatePreparationCircuitTest(unittest.TestCase):
+
+    def setUp(self):
+        self.n_qubits_range = range(3, 6)
+
+    def test_ground_state_particle_conserving(self):
+        """Test getting the ground state preparation circuit for a Hamiltonian
+        that conserves particle number."""
+        for n_qubits in self.n_qubits_range:
+            # Initialize a particle-number-conserving Hamiltonian
+            quadratic_hamiltonian = random_quadratic_hamiltonian(
+                    n_qubits, True, True)
+
+            # Compute the true ground state
+            sparse_operator = get_sparse_operator(quadratic_hamiltonian)
+            ground_energy, ground_state = get_ground_state(sparse_operator)
+
+            # Obtain the circuit
+            circuit_description, start_orbitals = (
+                    gaussian_state_preparation_circuit(quadratic_hamiltonian))
+
+            # Initialize the starting state
+            state = jw_slater_determinant(start_orbitals, 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)
+
+            # Check that the state obtained using the circuit is a ground state
+            difference = sparse_operator * state - ground_energy * state
+            discrepancy = 0.
+            if difference.nnz:
+                discrepancy = max(abs(difference.data))
+
+            self.assertTrue(discrepancy < EQ_TOLERANCE)
+
+    def test_ground_state_particle_nonconserving(self):
+        """Test getting the ground state preparation circuit for a Hamiltonian
+        that does not conserve particle number."""
+        for n_qubits in self.n_qubits_range:
+            # Initialize a particle-number-conserving Hamiltonian
+            quadratic_hamiltonian = random_quadratic_hamiltonian(
+                    n_qubits, False, True)
+
+            # Compute the true ground state
+            sparse_operator = get_sparse_operator(quadratic_hamiltonian)
+            ground_energy, ground_state = get_ground_state(sparse_operator)
+
+            # Obtain the circuit
+            circuit_description, start_orbitals = (
+                    gaussian_state_preparation_circuit(quadratic_hamiltonian))
+
+            # Initialize the starting state
+            state = jw_slater_determinant(start_orbitals, 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)
+
+            # Check that the state obtained using the circuit is a ground state
+            difference = sparse_operator * state - ground_energy * 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):
-            description, n_electrons = ground_state_preparation_circuit('a')
+            description, n_electrons = gaussian_state_preparation_circuit('a')
 
 
-class GetQuadraticHamiltonianGroundStateTest(unittest.TestCase):
+class JWGetGaussianStateTest(unittest.TestCase):
 
     def setUp(self):
         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
+    def test_ground_state_particle_conserving(self):
+        """Test getting the ground state of a Hamiltonian that conserves
+        particle number."""
         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)
+            quadratic_hamiltonian = random_quadratic_hamiltonian(
+                    n_qubits, True)
 
             # Compute the true ground state
-            fermion_operator = get_fermion_operator(quadratic_hamiltonian)
-            sparse_operator = jordan_wigner_sparse(fermion_operator)
+            sparse_operator = get_sparse_operator(quadratic_hamiltonian)
             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))
+            circuit_energy, circuit_state = jw_get_gaussian_state(
+                    quadratic_hamiltonian)
 
             # Check that the energies match
             self.assertAlmostEqual(ground_energy, circuit_energy)
@@ -86,35 +162,21 @@ class GetQuadraticHamiltonianGroundStateTest(unittest.TestCase):
 
             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
+    def test_ground_state_particle_nonconserving(self):
+        """Test getting the ground state of a Hamiltonian that does not
+        conserve particle number."""
         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)
+            quadratic_hamiltonian = random_quadratic_hamiltonian(
+                    n_qubits, False)
 
             # Compute the true ground state
-            fermion_operator = get_fermion_operator(quadratic_hamiltonian)
-            sparse_operator = jordan_wigner_sparse(fermion_operator)
+            sparse_operator = get_sparse_operator(quadratic_hamiltonian)
             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(
+                    jw_get_gaussian_state(
                         quadratic_hamiltonian))
 
             # Check that the energies match
@@ -129,10 +191,82 @@ class GetQuadraticHamiltonianGroundStateTest(unittest.TestCase):
 
             self.assertTrue(discrepancy < EQ_TOLERANCE)
 
+    def test_excited_state_particle_conserving(self):
+        """Test getting an excited state of a Hamiltonian that conserves
+        particle number."""
+        for n_qubits in self.n_qubits_range:
+            # Initialize a particle-number-conserving Hamiltonian
+            quadratic_hamiltonian = random_quadratic_hamiltonian(
+                    n_qubits, True)
+
+            # Pick some orbitals to occupy
+            num_occupied_orbitals = numpy.random.randint(1, n_qubits + 1)
+            occupied_orbitals = numpy.random.choice(
+                    range(n_qubits), num_occupied_orbitals, False)
+
+            # Compute the Gaussian state
+            circuit_energy, gaussian_state = jw_get_gaussian_state(
+                    quadratic_hamiltonian, occupied_orbitals)
+
+            # Compute the true energy
+            orbital_energies, constant = (
+                    quadratic_hamiltonian.orbital_energies())
+            energy = numpy.sum(orbital_energies[occupied_orbitals]) + constant
+
+            # Check that the energies match
+            self.assertAlmostEqual(energy, circuit_energy)
+
+            # Check that the state obtained using the circuit is an eigenstate
+            # with the correct eigenvalue
+            sparse_operator = get_sparse_operator(quadratic_hamiltonian)
+            difference = (sparse_operator * gaussian_state -
+                          energy * gaussian_state)
+            discrepancy = 0.
+            if difference.nnz:
+                discrepancy = max(abs(difference.data))
+
+            self.assertTrue(discrepancy < EQ_TOLERANCE)
+
+    def test_excited_state_particle_nonconserving(self):
+        """Test getting an excited state of a Hamiltonian that conserves
+        particle number."""
+        for n_qubits in self.n_qubits_range:
+            # Initialize a non-particle-number-conserving Hamiltonian
+            quadratic_hamiltonian = random_quadratic_hamiltonian(
+                    n_qubits, False)
+
+            # Pick some orbitals to occupy
+            num_occupied_orbitals = numpy.random.randint(1, n_qubits + 1)
+            occupied_orbitals = numpy.random.choice(
+                    range(n_qubits), num_occupied_orbitals, False)
+
+            # Compute the Gaussian state
+            circuit_energy, gaussian_state = jw_get_gaussian_state(
+                    quadratic_hamiltonian, occupied_orbitals)
+
+            # Compute the true energy
+            orbital_energies, constant = (
+                    quadratic_hamiltonian.orbital_energies())
+            energy = numpy.sum(orbital_energies[occupied_orbitals]) + constant
+
+            # Check that the energies match
+            self.assertAlmostEqual(energy, circuit_energy)
+
+            # Check that the state obtained using the circuit is an eigenstate
+            # with the correct eigenvalue
+            sparse_operator = get_sparse_operator(quadratic_hamiltonian)
+            difference = (sparse_operator * gaussian_state -
+                          energy * gaussian_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')
+            energy, state = jw_get_gaussian_state('a')
 
 
 class GivensDecompositionTest(unittest.TestCase):
@@ -151,7 +285,7 @@ class GivensDecompositionTest(unittest.TestCase):
             Q = Q[:m, :]
 
             # Get Givens decomposition of Q
-            V, givens_rotations, diagonal = givens_decomposition(Q)
+            givens_rotations, V, diagonal = givens_decomposition(Q)
 
             # Compute U
             U = numpy.eye(n, dtype=complex)
@@ -186,7 +320,7 @@ class GivensDecompositionTest(unittest.TestCase):
             Q = Q[:m, :]
 
             # Get Givens decomposition of Q
-            V, givens_rotations, diagonal = givens_decomposition(Q)
+            givens_rotations, V, diagonal = givens_decomposition(Q)
 
             # Compute U
             U = numpy.eye(n, dtype=complex)
@@ -224,12 +358,12 @@ class GivensDecompositionTest(unittest.TestCase):
         Q = Q[:m, :n]
 
         with self.assertRaises(ValueError):
-            V, givens_rotations, diagonal = givens_decomposition(Q)
+            givens_rotations, V, diagonal = givens_decomposition(Q)
 
     def test_identity(self):
         n = 3
         Q = numpy.eye(n, dtype=complex)
-        V, givens_rotations, diagonal = givens_decomposition(Q)
+        givens_rotations, V, diagonal = givens_decomposition(Q)
 
         # V should be the identity
         I = numpy.eye(n, dtype=complex)
@@ -250,7 +384,7 @@ class GivensDecompositionTest(unittest.TestCase):
         Q[0, 2] = 1.
         Q[1, 1] = 1.
         Q[2, 0] = 1.
-        V, givens_rotations, diagonal = givens_decomposition(Q)
+        givens_rotations, V, diagonal = givens_decomposition(Q)
 
         # There should be no Givens rotations
         self.assertEqual(givens_rotations, list())
@@ -273,7 +407,7 @@ class GivensDecompositionTest(unittest.TestCase):
         Q = Q[:m, :]
 
         # Get Givens decomposition of Q
-        V, givens_rotations, diagonal = givens_decomposition(Q)
+        givens_rotations, V, diagonal = givens_decomposition(Q)
 
         # There should be no Givens rotations
         self.assertEqual(givens_rotations, list())
@@ -299,18 +433,37 @@ class FermionicGaussianDecompositionTest(unittest.TestCase):
     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
+            antisymmetric_mat = random_antisymmetric_matrix(2 * n, real=True)
 
             # Get the diagonalizing fermionic unitary
             ferm_unitary = diagonalizing_fermionic_unitary(
-                    antisymmetric_matrix)
+                    antisymmetric_mat)
             lower_unitary = ferm_unitary[n:]
 
             # Get fermionic Gaussian decomposition of lower_unitary
-            left_unitary, decomposition, diagonal = (
+            decomposition, left_decomposition, diagonal, left_diagonal = (
                     fermionic_gaussian_decomposition(lower_unitary))
 
+            # Compute left_unitary
+            left_unitary = numpy.eye(n, dtype=complex)
+            for parallel_set in left_decomposition:
+                combined_op = numpy.eye(n, dtype=complex)
+                for op in reversed(parallel_set):
+                    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)
+                    givens_rotate(combined_op, givens_rotation, i, j)
+                left_unitary = combined_op.dot(left_unitary)
+            for i in range(n):
+                left_unitary[i] *= left_diagonal[i]
+            left_unitary = left_unitary.T
+            for i in range(n):
+                left_unitary[i] *= diagonal[i]
+
             # Check that left_unitary zeroes out the correct entries of
             # lower_unitary
             product = left_unitary.dot(lower_unitary)
@@ -353,14 +506,14 @@ class FermionicGaussianDecompositionTest(unittest.TestCase):
         n, p = (3, 7)
         rand_mat = numpy.random.randn(n, p)
         with self.assertRaises(ValueError):
-            left_unitary, decomposition, antidiagonal = (
+            decomposition, left_unitary, antidiagonal = (
                     fermionic_gaussian_decomposition(rand_mat))
 
     def test_bad_constraints(self):
         n = 3
         ones_mat = numpy.ones((n, 2 * n))
         with self.assertRaises(ValueError):
-            left_unitary, decomposition, antidiagonal = (
+            decomposition, left_unitary, antidiagonal = (
                     fermionic_gaussian_decomposition(ones_mat))
 
 
@@ -372,14 +525,8 @@ class DiagonalizingFermionicUnitaryTest(unittest.TestCase):
         self.chemical_potential = 2.
 
         # 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.hermitian_mat = random_hermitian_matrix(self.n_qubits)
+        self.antisymmetric_mat = random_antisymmetric_matrix(self.n_qubits)
 
         # Initialize a non-particle-number-conserving Hamiltonian
         self.quad_ham_npc = QuadraticHamiltonian(
@@ -475,6 +622,10 @@ class AntisymmetricCanonicalFormTest(unittest.TestCase):
                 else:
                     self.assertAlmostEqual(canonical[i, j], 0.)
 
+        diagonal = canonical[range(n), range(n, 2 * n)]
+        for i in range(n - 1):
+            self.assertTrue(diagonal[i] <= diagonal[i + 1])
+
 
 class GivensMatrixElementsTest(unittest.TestCase):
 
@@ -537,3 +688,37 @@ class DoubleGivensRotateTest(unittest.TestCase):
         G = givens_matrix_elements(v[0], v[1])
         with self.assertRaises(ValueError):
             double_givens_rotate(A, G, 0, 1, which='a')
+
+
+class JWSparseGivensRotationTest(unittest.TestCase):
+
+    def test_bad_input(self):
+        with self.assertRaises(ValueError):
+            givens_matrix = jw_sparse_givens_rotation(0, 2, 1., 1., 5)
+        with self.assertRaises(ValueError):
+            givens_matrix = jw_sparse_givens_rotation(4, 5, 1., 1., 5)
+
+
+def random_quadratic_hamiltonian(n_qubits,
+                                 conserves_particle_number=False,
+                                 real=False):
+    """Generate a random instance of QuadraticHamiltonian
+
+    Args:
+        n_qubits(int): the number of qubits
+        conserves_particle_number(bool): whether the returned Hamiltonian
+            should conserve particle number
+        real(bool): whether to use only real numbers
+
+    Returns:
+        QuadraticHamiltonian
+    """
+    constant = numpy.random.randn()
+    chemical_potential = numpy.random.randn()
+    hermitian_mat = random_hermitian_matrix(n_qubits, real)
+    if conserves_particle_number:
+        antisymmetric_mat = None
+    else:
+        antisymmetric_mat = random_antisymmetric_matrix(n_qubits, real)
+    return QuadraticHamiltonian(constant, hermitian_mat,
+                                antisymmetric_mat, chemical_potential)