Repository: quantumlib/OpenFermion
Pull request #111 (2017-11-17 05:40:21)

Minor restructuring and pretty docs
- I moved `jw_get_gaussian_state` from `utils/_slater_determinants.py` to `utils/_sparse_tools.py`. It seems to fit there better because it returns a sparse matrix. I reverted `get_ground_state` to the original version (without type-checking). You can get a ground state of a QuadraticHamiltonian as a sparse matrix by simply calling `jw_get_gaussian_state` on it.
- I made the Givens rotations docs prettier and added a reference to my arXiv paper where it is explained.


Patch:
diff --git a/src/openfermion/ops/_quadratic_hamiltonian.py b/src/openfermion/ops/_quadratic_hamiltonian.py
index d719ca7..72bc42d 100644
--- a/src/openfermion/ops/_quadratic_hamiltonian.py
+++ b/src/openfermion/ops/_quadratic_hamiltonian.py
@@ -32,15 +32,16 @@ class QuadraticHamiltonian(PolynomialTensor):
     .. math::
 
         \sum_{p, q} (M_{pq} - \mu \delta_{pq}) a^\dagger_p a_q
-        + \\frac12 \sum_{p, q} (A_{pq} a^\dagger_p a^\dagger_q + \\text{h.c.})
+        + \\frac12 \sum_{p, q}
+            (\\Delta_{pq} a^\dagger_p a^\dagger_q + \\text{h.c.})
         + \\text{constant}
 
     where
 
-        * :math:`M` is a Hermitian `n_qubits` x `n_qubits` matrix.
-        * :math:`A` is an antisymmetric `n_qubits` x `n_qubits` matrix.
-        * :math:`\mu` is a real number representing the chemical potential.
-        * :math:`\delta_{pq}` is the Kronecker delta symbol.
+        - :math:`M` is a Hermitian `n_qubits` x `n_qubits` matrix.
+        - :math:`\\Delta` is an antisymmetric `n_qubits` x `n_qubits` matrix.
+        - :math:`\mu` is a real number representing the chemical potential.
+        - :math:`\delta_{pq}` is the Kronecker delta symbol.
 
     We separate the chemical potential :math:`\mu` from :math:`M` so that
     we can use it to adjust the expectation value of the total number of
@@ -59,10 +60,13 @@ class QuadraticHamiltonian(PolynomialTensor):
             constant(float): A constant term in the operator.
             hermitian_part(ndarray): The matrix :math:`M`, which represents the
                 coefficients of the particle-number-conserving terms.
-                This is an `n_qubits` x `n_qubits` numpy array of complex numbers.
-            antisymmetric_part(ndarray): The matrix :math:`A`, which represents
-                the coefficients of the non-particle-number-conserving terms.
-                This is an `n_qubits` x `n_qubits` numpy array of complex numbers.
+                This is an `n_qubits` x `n_qubits` numpy array of complex
+                numbers.
+            antisymmetric_part(ndarray): The matrix :math:`\\Delta`,
+                which represents the coefficients of the
+                non-particle-number-conserving terms.
+                This is an `n_qubits` x `n_qubits` numpy array of complex
+                numbers.
             chemical_potential(float): The chemical potential :math:`\mu`.
         """
         n_qubits = hermitian_part.shape[0]
@@ -126,7 +130,7 @@ class QuadraticHamiltonian(PolynomialTensor):
 
         Any quadratic Hamiltonian is unitarily equivalent to a Hamiltonian
         of the form
-        
+
         .. math::
 
             \sum_{j} \\varepsilon_j a^\dagger_j a_j + \\text{constant}.
@@ -175,8 +179,9 @@ class QuadraticHamiltonian(PolynomialTensor):
 
             f_{j + N} = \\frac{i}{\\sqrt{2}} (a^\dagger_j - a_j)
 
-        and :math:`A` is a (2 * `n_qubits`) x (2 * `n_qubits`) real antisymmetric matrix.
-        This function returns the matrix :math:`A` and the constant.
+        and :math:`A` is a (2 * `n_qubits`) x (2 * `n_qubits`) real
+        antisymmetric matrix. This function returns the matrix
+        :math:`A` and the constant.
         """
         hermitian_part = self.combined_hermitian_part
         antisymmetric_part = self.antisymmetric_part
diff --git a/src/openfermion/utils/__init__.py b/src/openfermion/utils/__init__.py
index 7092aa6..3eaf6b9 100644
--- a/src/openfermion/utils/__init__.py
+++ b/src/openfermion/utils/__init__.py
@@ -22,8 +22,7 @@ from ._operator_utils import (commutator, count_qubits,
 
 from ._slater_determinants import (fermionic_gaussian_decomposition,
                                    givens_decomposition,
-                                   gaussian_state_preparation_circuit,
-                                   jw_get_gaussian_state)
+                                   gaussian_state_preparation_circuit)
 
 from ._sparse_tools import (expectation,
                             expectation_computational_basis_state,
@@ -33,6 +32,7 @@ from ._sparse_tools import (expectation,
                             is_hermitian,
                             jordan_wigner_sparse,
                             jw_hartree_fock_state,
+                            jw_get_gaussian_state,
                             jw_get_ground_states_by_particle_number,
                             jw_number_restrict_operator,
                             jw_slater_determinant,
diff --git a/src/openfermion/utils/_slater_determinants.py b/src/openfermion/utils/_slater_determinants.py
index af099bc..eac54d8 100644
--- a/src/openfermion/utils/_slater_determinants.py
+++ b/src/openfermion/utils/_slater_determinants.py
@@ -15,16 +15,11 @@ Slater determinants and fermionic Gaussian states."""
 from __future__ import absolute_import
 
 import numpy
-from scipy.sparse import csc_matrix, eye
 
 from openfermion.config import EQ_TOLERANCE
 from openfermion.ops import QuadraticHamiltonian
 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)
+        diagonalizing_fermionic_unitary, swap_columns)
 
 
 def gaussian_state_preparation_circuit(
@@ -34,7 +29,8 @@ def gaussian_state_preparation_circuit(
 
     Fermionic Gaussian states can be regarded as eigenstates of quadratic
     Hamiltonians. If the Hamiltonian conserves particle number, then these are
-    just Slater determinants.
+    just Slater determinants. See arXiv:1711.05395 for a detailed description
+    of how this procedure works.
 
     Args:
         quadratic_hamiltonian(QuadraticHamiltonian):
@@ -45,16 +41,19 @@ def gaussian_state_preparation_circuit(
             (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]):
+    Returns
+    -------
+        circuit_description (list[tuple]):
             A list of operations describing the circuit. Each operation
             is a tuple of objects describing elementary operations that
             can be performed in parallel. Each elementary operation
             is either the string 'pht', indicating a particle-hole
             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.
-        start_orbitals(list):
+            the form :math:`(i, j, \\theta, \\phi)`,
+            indicating a Givens rotation
+            of modes :math:`i` and :math:`j` by angles :math:`\\theta`
+            and :math:`\\phi`.
+        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 these
@@ -115,110 +114,74 @@ def gaussian_state_preparation_circuit(
     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 fermionic_gaussian_decomposition(unitary_rows):
+    """Decompose a matrix into a sequence of Givens rotations and
+    particle-hole transformations on the last fermionic mode.
 
-    Eigenstates of a quadratic Hamiltonian are also known as fermionic
-    Gaussian states.
+    The input is an :math:`N \\times 2N` matrix :math:`W` with orthonormal
+    rows. Furthermore, :math:`W` must have the block form
 
-    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.
+    .. math::
 
-    Returns:
-        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.')
+        W = ( W_1 \hspace{4pt} W_2 )
 
-    n_qubits = quadratic_hamiltonian.n_qubits
+    where :math:`W_1` and :math:`W_2` satisfy
 
-    # Compute the energy
-    orbital_energies, constant = quadratic_hamiltonian.orbital_energies()
-    if occupied_orbitals is None:
-        # The ground energy is desired
-        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 Gaussian state
-    circuit_description, start_orbitals = gaussian_state_preparation_circuit(
-            quadratic_hamiltonian, occupied_orbitals)
-
-    # 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)
+    .. math::
 
-    return energy, state
+        W_1  W_1^\dagger + W_2  W_2^\dagger &= I
 
+        W_1  W_2^T + W_2  W_1^T &= 0.
 
-def fermionic_gaussian_decomposition(unitary_rows):
-    """Decompose a matrix into a sequence of Givens rotations and
-    particle-hole transformations on the last fermionic mode.
+    Then :math:`W` can be decomposed as
 
-    The input is an n x (2 * n) matrix W with orthonormal rows.
-    Furthermore, W has the block form::
+    .. math::
 
-        W = [ W_1  |  W_2 ]
+        V  W  U^\dagger = ( 0 \hspace{6pt} D )
 
-    where W_1 and W_2 satisfy::
+    where :math:`V` and :math:`U` are unitary matrices and :math:`D`
+    is a diagonal unitary matrix. Furthermore, :math:`U` can be decomposed
+    as follows:
 
-        W_1 * W_1^\dagger + W_2 * W_2^\dagger = I
-        W_1 * W_2^T + W_2 * W_1^T = 0
+    .. math::
 
-    W can be decomposed as::
+        U = B G_{k} \cdots B G_3 G_2 B G_1 B,
 
-        V * W * U^\dagger = [ 0  |  D ]
+    where each :math:`G_i` is a Givens rotation, and :math:`B` represents
+    swapping the :math:`N`-th column with the :math:`2N`-th column,
+    which corresponds to a particle-hole transformation
+    on the last fermionic mode. This particle-hole transformation maps
+    :math:`a^\dagger_N` to :math:`a_N` and vice versa, while leaving the
+    other fermionic ladder operators invariant.
 
-    where V and U are unitary matrices and D is a diagonal unitary matrix.
-    Furthermore, we can decompose U as a sequence of Givens rotations
-    and particle-hole transformations on the last fermionic mode.
-    This particle-hole transformation maps a^\dagger_n to a^n and vice
-    versa, while leaving the other ladder operators invariant.
-
-    The decomposition of U is returned as a list of tuples of objects
+    The decomposition of :math:`U` is returned as a list of tuples of objects
     describing rotations and particle-hole transformations. The list looks
     something like [('pht', ), (G_1, ), ('pht', G_2), ... ].
     The objects within a tuple are either the string 'pht', which indicates
     a particle-hole transformation on the last fermionic mode, or a tuple
-    of the form (i, j, theta, phi), which indicates a Givens rotation
-    of rows i and j by angles theta and phi.
+    of the form :math:`(i, j, \\theta, \\phi)`, which indicates a
+    Givens rotation of rows :math:`i` and :math:`j` by angles
+    :math:`\\theta` and :math:`\\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.
+    The matrix :math:`V^T 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:
-        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^*.
+    Returns
+    -------
+        decomposition (list[tuple]):
+            The decomposition of :math:`U`.
+        left_decomposition (list[tuple]):
+            The decomposition of :math:`V^T D^*`.
+        diagonal (ndarray):
+            A list of the nonzero entries of :math:`D`.
+        left_diagonal (ndarray):
+            A list of the nonzero entries left from the decomposition
+            of :math:`V^T D^*`.
     """
     current_matrix = numpy.copy(unitary_rows)
     n, p = current_matrix.shape
@@ -360,38 +323,55 @@ def fermionic_gaussian_decomposition(unitary_rows):
 def givens_decomposition(unitary_rows):
     """Decompose a matrix into a sequence of Givens rotations.
 
-    The input is an m x n matrix Q with m <= n. The rows of Q are orthonormal.
-    Q can be decomposed as follows:
+    The input is an :math:`m \\times n` matrix :math:`Q` with :math:`m \leq n`.
+    The rows of :math:`Q` are orthonormal.
+    :math:`Q` can be decomposed as follows:
 
-        V * Q * U^\dagger = D
+    .. math::
 
-    where V and U are unitary matrices, and D is an m x n matrix with the
-    first m columns forming a diagonal matrix and the rest of the columns
-    being zero. Furthermore, we can decompose U as
+        V Q U^\dagger = D
 
-        U = G_k * ... * G_1
+    where :math:`V` and :math:`U` are unitary matrices, and :math:`D`
+    is an :math:`m \\times n` matrix with the
+    first :math:`m` columns forming a diagonal matrix and the rest of the
+    columns being zero. Furthermore, we can decompose :math:`U` as
 
-    where G_1, ..., G_k are complex Givens rotations, which are invertible
-    n x n matrices. We describe a complex Givens rotation by the column
-    indices (i, j) that it acts on, plus two angles (theta, phi) that
-    characterize the corresponding 2x2 unitary matrix
+    .. math::
 
-        [ cos(theta)    -e^{i phi} sin(theta) ]
-        [ sin(theta)     e^{i phi} cos(theta) ]
+        U = G_k ... G_1
+
+    where :math:`G_1, \\ldots, G_k` are complex Givens rotations.
+    A Givens rotation is a rotation within the two-dimensional subspace
+    spanned by two coordinate axes. Within the two relevant coordinate
+    axes, a Givens rotation has the form
+
+    .. math::
+
+        \\begin{pmatrix}
+            \\cos(\\theta) & -e^{i \\phi} \\sin(\\theta) \\\\
+            \\sin(\\theta) &     e^{i \\phi} \\cos(\\theta)
+        \\end{pmatrix}.
 
     Args:
         unitary_rows: A numpy array or matrix with orthonormal rows,
             representing the matrix Q.
 
-    Returns:
-        givens_rotations: A list of tuples of objects describing Givens
+    Returns
+    -------
+        givens_rotations (list[tuple]):
+            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.
+            form :math:`(i, j, \\theta, \\phi)`, which represents a
+            Givens rotation of coordinates
+            :math:`i` and :math:`j` by angles :math:`\\theta` and
+            :math:`\\phi`.
+        left_unitary (ndarray):
+            An :math:`m \\times m` numpy array representing the matrix
+            :math:`V`.
+        diagonal (ndarray):
+            A list of the nonzero entries of :math:`D`.
     """
     current_matrix = numpy.copy(unitary_rows)
     m, n = current_matrix.shape
@@ -603,40 +583,3 @@ def double_givens_rotate(operator, givens_rotation, i, j, which='row'):
                       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']])
diff --git a/src/openfermion/utils/_sparse_tools.py b/src/openfermion/utils/_sparse_tools.py
index f13c0ef..7a9414c 100644
--- a/src/openfermion/utils/_sparse_tools.py
+++ b/src/openfermion/utils/_sparse_tools.py
@@ -28,7 +28,8 @@ from openfermion.config import *
 from openfermion.ops import (FermionOperator, hermitian_conjugated,
                              normal_ordered, number_operator,
                              QuadraticHamiltonian, QubitOperator)
-from openfermion.utils import commutator, fourier_transform, Grid
+from openfermion.utils import (commutator, fourier_transform,
+                               gaussian_state_preparation_circuit, Grid)
 from openfermion.hamiltonians._jellium import (momentum_vector,
                                                position_vector,
                                                grid_indices)
@@ -367,6 +368,104 @@ def jw_get_ground_states_by_particle_number(sparse_operator, particle_number,
     return ground_energy, ground_states
 
 
+def jw_get_gaussian_state(quadratic_hamiltonian, occupied_orbitals=None):
+    """Compute an 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
+    -------
+        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 energy
+    orbital_energies, constant = quadratic_hamiltonian.orbital_energies()
+    if occupied_orbitals is None:
+        # The ground energy is desired
+        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 Gaussian state
+    circuit_description, start_orbitals = gaussian_state_preparation_circuit(
+            quadratic_hamiltonian, occupied_orbitals)
+
+    # 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)
+
+    return energy, state
+
+
+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 = scipy.sparse.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 = scipy.sparse.eye(2 ** i, format='csc')
+    right_eye = scipy.sparse.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 = scipy.sparse.eye(2 ** (n_qubits - 1), format='csc')
+    return kronecker_operators([left_eye, pauli_matrix_map['X']])
+
+
 def get_density_matrix(states, probabilities):
     n_qubits = states[0].shape[0]
     density_matrix = scipy.sparse.csc_matrix(
@@ -386,34 +485,24 @@ def is_hermitian(sparse_operator):
     return True
 
 
-def get_ground_state(operator):
+def get_ground_state(sparse_operator):
     """Compute lowest eigenvalue and eigenstate.
-
     Returns:
         eigenvalue: The lowest eigenvalue, a float.
         eigenstate: The lowest eigenstate in scipy.sparse csc format.
     """
-    if isinstance(operator, QuadraticHamiltonian):
-        from openfermion.utils._slater_determinants import (
-                jw_get_gaussian_state)
-        eigenvalue, eigenstate = jw_get_gaussian_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.')
+    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)
 
-    return eigenvalue, eigenstate
+    order = numpy.argsort(values)
+    values = values[order]
+    vectors = vectors[:, order]
+    eigenvalue = values[0]
+    eigenstate = scipy.sparse.csc_matrix(vectors[:, 0])
+    return eigenvalue, eigenstate.T
 
 
 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 f4e1d3a..5271f0b 100644
--- a/src/openfermion/utils/_slater_determinants_test.py
+++ b/src/openfermion/utils/_slater_determinants_test.py
@@ -20,8 +20,7 @@ from scipy.linalg import qr
 from openfermion.config import EQ_TOLERANCE
 from openfermion.ops import QuadraticHamiltonian
 from openfermion.ops._quadratic_hamiltonian import (
-        diagonalizing_fermionic_unitary,
-        swap_rows)
+        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
@@ -34,7 +33,8 @@ from openfermion.utils import (fermionic_gaussian_decomposition,
 from openfermion.utils._slater_determinants import (
         double_givens_rotate,
         givens_rotate,
-        givens_matrix_elements,
+        givens_matrix_elements)
+from openfermion.utils._sparse_tools import (
         jw_sparse_givens_rotation,
         jw_sparse_particle_hole_transformation_last_mode)
 
@@ -128,146 +128,6 @@ class GaussianStatePreparationCircuitTest(unittest.TestCase):
             description, n_electrons = gaussian_state_preparation_circuit('a')
 
 
-class JWGetGaussianStateTest(unittest.TestCase):
-
-    def setUp(self):
-        self.n_qubits_range = range(2, 10)
-
-    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:
-            # Initialize a particle-number-conserving Hamiltonian
-            quadratic_hamiltonian = random_quadratic_hamiltonian(
-                    n_qubits, True)
-
-            # Compute the true ground state
-            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_gaussian_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_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:
-            # Initialize a non-particle-number-conserving Hamiltonian
-            quadratic_hamiltonian = random_quadratic_hamiltonian(
-                    n_qubits, False)
-
-            # Compute the true ground state
-            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_gaussian_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_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_gaussian_state('a')
-
-
 class GivensDecompositionTest(unittest.TestCase):
 
     def setUp(self):
diff --git a/src/openfermion/utils/_sparse_tools_test.py b/src/openfermion/utils/_sparse_tools_test.py
index da4d3b9..6b374b6 100644
--- a/src/openfermion/utils/_sparse_tools_test.py
+++ b/src/openfermion/utils/_sparse_tools_test.py
@@ -26,6 +26,8 @@ from openfermion.transforms import (get_fermion_operator, get_sparse_operator,
 from openfermion.utils import fourier_transform, Grid
 from openfermion.utils._jellium_hf_state import (
     lowest_single_particle_energy_states)
+from openfermion.utils._slater_determinants_test import (
+        random_quadratic_hamiltonian)
 from openfermion.utils._sparse_tools import *
 
 
@@ -259,6 +261,146 @@ class JWGetGroundStatesByParticleNumberTest(unittest.TestCase):
         return
 
 
+class JWGetGaussianStateTest(unittest.TestCase):
+
+    def setUp(self):
+        self.n_qubits_range = range(2, 10)
+
+    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:
+            # Initialize a particle-number-conserving Hamiltonian
+            quadratic_hamiltonian = random_quadratic_hamiltonian(
+                    n_qubits, True)
+
+            # Compute the true ground state
+            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_gaussian_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_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:
+            # Initialize a non-particle-number-conserving Hamiltonian
+            quadratic_hamiltonian = random_quadratic_hamiltonian(
+                    n_qubits, False)
+
+            # Compute the true ground state
+            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_gaussian_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_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_gaussian_state('a')
+
+
 class GroundStateTest(unittest.TestCase):
     def test_get_ground_state_hermitian(self):
         ground = get_ground_state(get_sparse_operator(
@@ -276,36 +418,6 @@ class GroundStateTest(unittest.TestCase):
         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):