Repository: quantumlib/OpenFermion
Pull request #523 (2019-05-23 20:35:05)

Adds a first draft of tools to work with number preserving sparse operators
These are the tools (and the tests for them) that I wrote to create scipy sparse operators directly in some number (and possibly spin / excitation rank) subspace. Feedback about what to do to get them ready to pull into OpenFermion would be great.

Thanks,

Bill

ncrubin: @Wuggins This would be great to have but no immediate pressing need.  Let me know what you'd like to do with this PR.  I am happy to keep it open for the time being if you want to do more development.  
Wuggins: Hey @kevinsung, are there other changes I should make?
kevinsung: You missed a few of my comments; perhaps you need to load the "hidden conversations" to see them?
Wuggins: @kevinsung. Thank you, that was a silly mistake to make! I fixed those remaining things and have pushed the changes again.

Patch:
diff --git a/src/openfermion/transforms/_conversion.py b/src/openfermion/transforms/_conversion.py
index f1ffdf8..1d51171 100644
--- a/src/openfermion/transforms/_conversion.py
+++ b/src/openfermion/transforms/_conversion.py
@@ -16,6 +16,7 @@ from __future__ import absolute_import
 import itertools
 
 import numpy
+import scipy
 from future.utils import iteritems
 
 from openfermion.config import EQ_TOLERANCE
@@ -304,7 +305,7 @@ def get_quadratic_hamiltonian(fermion_operator,
 def get_diagonal_coulomb_hamiltonian(fermion_operator,
                                      n_qubits=None,
                                      ignore_incompatible_terms=False):
-    """Convert a FermionOperator to a DiagonalCoulombHamiltonian.
+    r"""Convert a FermionOperator to a DiagonalCoulombHamiltonian.
     
     Args:
         fermion_operator(FermionOperator): The operator to convert.
@@ -556,3 +557,331 @@ def get_boson_operator(operator, hbar=1.):
                         "supported for get_boson_operator.")
 
     return boson_operator
+
+
+def get_number_preserving_sparse_operator(
+        fermion_op,
+        num_qubits,
+        num_electrons,
+        spin_preserving=False,
+        reference_determinant=None,
+        excitation_level=None):
+    """Initialize a Scipy sparse matrix in a specific symmetry sector.
+
+    This method initializes a Scipy sparse matrix from a FermionOperator,
+    explicitly working in a particular particle number sector. Optionally, it
+    can also restrict the space to contain only states with a particular Sz.
+
+    Finally, the Hilbert space can also be restricted to only those states
+    which are reachable by excitations up to a fixed rank from an initial
+    reference determinant.
+
+    Args:
+        fermion_op(FermionOperator): An instance of the FermionOperator class.
+            It should not contain terms which do not preserve particle number.
+            If spin_preserving is set to True it should also not contain terms
+            which do not preserve the Sz (it is assumed that the ordering of
+            the indices goes alpha, beta, alpha, beta, ...).
+        num_qubits(int): The total number of qubits / spin-orbitals in the
+            system.
+        num_electrons(int): The number of particles in the desired Hilbert
+            space.
+        spin_preserving(bool): Whether or not the constructed operator should
+            be defined in a space which has support only on states with the
+            same Sz value as the reference_determinant.
+        reference_determinant(list(bool)): A list, whose length is equal to
+            num_qubits, which specifies which orbitals should be occupied in
+            the reference state. If spin_preserving is set to True then the Sz
+            value of this reference state determines the Sz value of the
+            symmetry sector in which the generated operator acts. If a value
+            for excitation_level is provided then the excitations are generated
+            with respect to the reference state. In any case, the ordering of
+            the states in the matrix representation of the operator depends on
+            reference_determinant and the state corresponding to
+            reference_determinant is the vector [1.0, 0.0, 0.0 ... 0.0]. Can be
+            set to None in order to take the first num_electrons orbitals to be
+            the occupied orbitals.
+        excitation_level(int): The number of excitations from the reference
+            state which should be included in the generated operator's matrix
+            representation. Can be set to None to include all levels of
+            excitation.
+
+    Returns:
+        sparse_op(scipy.sparse.csc_matrix): A sparse matrix representation of
+            fermion_op in the basis set by the arguments.
+    """
+
+    # We use the Hartree-Fock determinant as a reference if none is provided.
+    if reference_determinant is None:
+        reference_determinant = numpy.array([i < num_electrons for i in
+                                             range(num_qubits)])
+    else:
+        reference_determinant = numpy.asarray(reference_determinant)
+
+    if excitation_level is None:
+        excitation_level = num_electrons
+
+    state_array = numpy.asarray(list(_iterate_basis_(
+        reference_determinant, excitation_level, spin_preserving)))
+    # Create a 1d array with each determinant encoded as an integer for sorting purposes.
+    int_state_array = state_array.dot(1 << numpy.arange(state_array.shape[1])[::-1])
+    sorting_indices = numpy.argsort(int_state_array)
+
+    space_size = state_array.shape[0]
+
+    fermion_op = normal_ordered(fermion_op)
+
+    sparse_op = scipy.sparse.csc_matrix((space_size, space_size), dtype=float)
+
+    for term, coefficient in fermion_op.terms.items():
+        if len(term) == 0:
+            constant = coefficient * scipy.sparse.identity(
+                space_size, dtype=float, format='csc')
+
+            sparse_op += constant
+
+        else:
+            term_op = _build_term_op_(term,
+                                        state_array,
+                                        int_state_array,
+                                        sorting_indices)
+
+            sparse_op += coefficient * term_op
+
+    return sparse_op
+
+
+def _iterate_basis_(reference_determinant, excitation_level, spin_preserving):
+    """A helper method which iterates over the specified basis states.
+
+    Note that this method always yields the states in order of their excitation
+    rank from the reference_determinant.
+
+    Args:
+        reference_determinant(list(bool)): A list of bools which indicates
+            which orbitals are occupied and which are unoccupied in the
+            reference state.
+        excitation_level(int): The maximum excitation rank to iterate over.
+        spin_preserving(bool): A bool which, if set to True, constrains the
+            method to iterate over only those states which have the same Sz as
+            reference_determinant.
+
+    Yields:
+        Lists of bools which indicate which orbitals are occupied and which are
+            unoccupied in the current determinant.
+    """
+    if not spin_preserving:
+        for order in range(excitation_level + 1):
+            for determinant in _iterate_basis_order_(reference_determinant, order):
+                yield determinant
+
+    else:
+        alpha_excitation_level = min((numpy.sum(reference_determinant[::2]),
+                                      excitation_level))
+        beta_excitation_level = min((numpy.sum(reference_determinant[1::2]),
+                                     excitation_level))
+
+        for order in range(excitation_level + 1):
+            for alpha_order in range(alpha_excitation_level + 1):
+                beta_order = order - alpha_order
+                if (beta_order < 0 or beta_order > beta_excitation_level):
+                    continue
+
+                for determinant in _iterate_basis_spin_order_(
+                        reference_determinant, alpha_order, beta_order):
+                    yield determinant
+
+
+
+def _iterate_basis_order_(reference_determinant, order):
+    """A helper for iterating over determinants of a fixed excitation rank.
+
+    Args:
+        reference_determinant(list(bool)): The reference state with respect to
+            which we are iterating over excited determinants.
+        order(int): The number of excitations from the modes which are occupied
+            in the reference_determinant.
+
+    Yields:
+        Lists of bools which indicate which orbitals are occupied and which are
+            unoccupied in the current determinant.
+        """
+    occupied_indices = numpy.where(reference_determinant)[0]
+    unoccupied_indices = numpy.where(numpy.invert(reference_determinant))[0]
+
+    for occ_ind, unocc_ind in itertools.product(
+            itertools.combinations(occupied_indices, order),
+            itertools.combinations(unoccupied_indices, order)):
+        basis_state = reference_determinant.copy()
+
+        occ_ind = list(occ_ind)
+        unocc_ind = list(unocc_ind)
+
+        basis_state[occ_ind] = False
+        basis_state[unocc_ind] = True
+
+        yield basis_state
+
+
+def _iterate_basis_spin_order_(reference_determinant, alpha_order, beta_order):
+    """Iterates over states with a fixed excitation rank for each spin sector.
+
+    This helper method assumes that the two spin sectors are interleaved:
+    [1_alpha, 1_beta, 2_alpha, 2_beta, ...].
+
+    Args:
+        reference_determinant(list(bool)): The reference state with respect to
+            which we are iterating over excited determinants.
+        alpha_order(int): The number of excitations from the alpha spin sector
+            of the reference_determinant.
+        beta_order(int): The number of excitations from the beta spin sector of
+            the reference_determinant.
+
+    Yields:
+        Lists of bools which indicate which orbitals are occupied and which are
+            unoccupied in the current determinant.
+        """
+    occupied_alpha_indices = numpy.where(
+        reference_determinant[::2])[0] * 2
+    unoccupied_alpha_indices = numpy.where(
+        numpy.invert(reference_determinant[::2]))[0] * 2
+    occupied_beta_indices = numpy.where(
+        reference_determinant[1::2])[0] * 2 + 1
+    unoccupied_beta_indices = numpy.where(
+        numpy.invert(reference_determinant[1::2]))[0] * 2 + 1
+
+    for (alpha_occ_ind,
+            alpha_unocc_ind,
+            beta_occ_ind,
+            beta_unocc_ind) in itertools.product(
+            itertools.combinations(occupied_alpha_indices, alpha_order),
+            itertools.combinations(unoccupied_alpha_indices, alpha_order),
+            itertools.combinations(occupied_beta_indices, beta_order),
+            itertools.combinations(unoccupied_beta_indices, beta_order)):
+        basis_state = reference_determinant.copy()
+
+        alpha_occ_ind = list(alpha_occ_ind)
+        alpha_unocc_ind = list(alpha_unocc_ind)
+        beta_occ_ind = list(beta_occ_ind)
+        beta_unocc_ind = list(beta_unocc_ind)
+
+        basis_state[alpha_occ_ind] = False
+        basis_state[alpha_unocc_ind] = True
+        basis_state[beta_occ_ind] = False
+        basis_state[beta_unocc_ind] = True
+
+        yield basis_state
+
+
+def _build_term_op_(term, state_array, int_state_array, sorting_indices):
+    """Builds a scipy sparse representation of a term from a FermionOperator.
+
+    Args:
+        term(tuple of tuple(int, int)s): The argument is a tuple of tuples
+            representing a product of normal ordered fermionic creation and
+            annihilation operators, each of which is of the form (int, int)
+            where the first int indicates which site the operator acts on and
+            the second int indicates whether the operator is a creation
+            operator (1) or an annihilation operator (0). See the
+            implementation of FermionOperator for more details.
+        state_array(ndarray(bool)): A Numpy array which encodes each of the
+            determinants in the space we are working in a bools which indicate
+            the occupation of each mode. See the implementation of
+            get_number_preserving_sparse_operator for more details.
+        int_state_array(ndarray(int)): A one dimensional Numpy array which
+            encodes the integer representation of the binary number
+            corresponding to each determinant in state_array.
+        sorting_indices(ndarray.view): A Numpy view which sorts
+            int_state_array. This, together with int_state_array, allows for a
+            quick lookup of the position of a particular determinant in
+            state_array by converting it to its integer representation and
+            searching through the sorted int_state_array.
+
+    Raises:
+        ValueError: If term does not represent a particle number conserving
+            operator.
+
+    Returns:
+        A scipy.sparse.csc_matrix which corresponds to the operator specified
+            by term expressed in the basis corresponding to the other arguments
+            of the method."""
+
+    space_size = state_array.shape[0]
+
+    needs_to_be_occupied = []
+    needs_to_be_unoccupied = []
+
+    # We keep track of the number of creation and annihilation operators and
+    # ensure that there are an equal number of them in order to help detect
+    # invalid inputs.
+    delta = 0
+    for index, op_type in reversed(term):
+        if op_type == 0:
+            needs_to_be_occupied.append(index)
+            delta -= 1
+        else:
+            if index not in needs_to_be_occupied:
+                needs_to_be_unoccupied.append(index)
+            delta += 1
+
+    if delta != 0:
+        raise ValueError(
+            "The supplied operator doesn't preserve particle number")
+
+
+    # We search for every state which has the necessary orbitals occupied and
+    # unoccupied in order to not be immediately zeroed out based on the
+    # creation and annihilation operators specified in term.
+    maybe_valid_states = numpy.where(
+        numpy.logical_and(
+            numpy.all(state_array[:,needs_to_be_occupied], axis=1),
+            numpy.logical_not(numpy.any(state_array[:, needs_to_be_unoccupied], axis=1))))[0]
+
+    data = []
+    row_ind = []
+    col_ind = []
+    shape = (space_size, space_size)
+
+    # For each state that is not immediately zeroed out by the action of our
+    # operator we check to see if the determinant which this state gets mapped
+    # to is in the space we are considering.
+    # Note that a failure to find any state does not necessarily indicate that
+    # term specifies an invalid operator. For example, if we are restricting
+    # ourselves to double excitations from a fixed reference state then the
+    # action of term on some of our basis states may lead to determinants with
+    # more than two excitations from the reference. These more than double
+    # excited determinants are not included in the matrix representation (and
+    # hence, will not be present in state_array).
+    for k, state in enumerate(maybe_valid_states):
+        determinant = state_array[state, :]
+        target_determinant = determinant.copy()
+
+        parity = 1
+        for i, _ in reversed(term):
+            area_to_check = target_determinant[0:i]
+            parity *= (-1) ** numpy.sum(area_to_check)
+
+            target_determinant[i] = not target_determinant[i]
+
+        int_encoding = target_determinant.dot(
+            1 << numpy.arange(target_determinant.size)[::-1])
+
+        target_state_index_sorted = numpy.searchsorted(int_state_array,
+                                                       int_encoding,
+                                                       sorter=sorting_indices)
+
+        target_state = sorting_indices[target_state_index_sorted]
+
+        if int_state_array[target_state] == int_encoding:
+            # Then target state is in the space considered:
+            data.append(parity)
+            row_ind.append(target_state)
+            col_ind.append(state)
+
+    data = numpy.asarray(data)
+    row_ind = numpy.asarray(row_ind)
+    col_ind = numpy.asarray(col_ind)
+
+    term_op = scipy.sparse.csc_matrix((data, (row_ind, col_ind)), shape=shape)
+
+    return term_op


Test patch:
diff --git a/src/openfermion/transforms/_conversion_test.py b/src/openfermion/transforms/_conversion_test.py
index ad7985f..b0674f8 100644
--- a/src/openfermion/transforms/_conversion_test.py
+++ b/src/openfermion/transforms/_conversion_test.py
@@ -14,8 +14,10 @@
 from __future__ import absolute_import
 
 import unittest
+import os
 
 import numpy
+import scipy
 import pytest
 
 from openfermion.hamiltonians import fermi_hubbard
@@ -26,13 +28,19 @@ from openfermion.ops import (BosonOperator,
                              QuadOperator,
                              QubitOperator)
 from openfermion.transforms import jordan_wigner
-from openfermion.utils import is_hermitian, normal_ordered
+from openfermion.utils import (is_hermitian,
+                               normal_ordered,
+                               jw_number_restrict_operator,
+                               jw_sz_restrict_operator,
+                               sparse_eigenspectrum)
 
+from openfermion.hamiltonians._molecular_data import MolecularData
 from openfermion.ops._interaction_operator import InteractionOperatorError
 from openfermion.ops._quadratic_hamiltonian import QuadraticHamiltonianError
 from openfermion.utils._testing_utils import (
         random_hermitian_matrix,
         random_quadratic_hamiltonian)
+from openfermion.config import THIS_DIRECTORY
 
 from openfermion.transforms._conversion import (
         get_boson_operator,
@@ -41,7 +49,9 @@ from openfermion.transforms._conversion import (
         get_interaction_operator,
         get_quad_operator,
         get_quadratic_hamiltonian,
-        get_sparse_operator)
+        get_sparse_operator,
+        get_number_preserving_sparse_operator,
+        _iterate_basis_)
 
 
 class GetInteractionOperatorTest(unittest.TestCase):
@@ -507,3 +517,505 @@ def test_get_fermion_operator_majorana_operator():
 def test_get_fermion_operator_wrong_type():
     with pytest.raises(TypeError):
         _ = get_fermion_operator(QubitOperator())
+
+
+class GetNumberPreservingSparseOperatorIntegrationTestLiH(unittest.TestCase):
+
+    def setUp(self):
+        # Set up molecule.
+        geometry = [('Li', (0., 0., 0.)), ('H', (0., 0., 1.45))]
+        basis = 'sto-3g'
+        multiplicity = 1
+        filename = os.path.join(THIS_DIRECTORY, 'data',
+                                'H1-Li1_sto-3g_singlet_1.45')
+        self.molecule = MolecularData(
+            geometry, basis, multiplicity, filename=filename)
+        self.molecule.load()
+
+        # Get molecular Hamiltonian.
+        self.molecular_hamiltonian = self.molecule.get_molecular_hamiltonian()
+
+        self.hubbard_hamiltonian = fermi_hubbard(
+            2, 2, 1.0, 4.0,
+            chemical_potential=.2,
+            magnetic_field=0.0,
+            spinless=False)
+
+
+    def test_number_on_reference(self):
+        sum_n_op = FermionOperator()
+        sum_sparse_n_op = get_number_preserving_sparse_operator(
+            sum_n_op,
+            self.molecule.n_qubits,
+            self.molecule.n_electrons,
+            spin_preserving=False)
+
+        space_size = sum_sparse_n_op.shape[0]
+        reference = numpy.zeros((space_size))
+        reference[0] = 1.0
+
+        for i in range(self.molecule.n_qubits):
+            n_op = FermionOperator(((i, 1), (i, 0)))
+            sum_n_op += n_op
+
+            sparse_n_op = get_number_preserving_sparse_operator(
+                n_op,
+                self.molecule.n_qubits,
+                self.molecule.n_electrons,
+                spin_preserving=False)
+
+            sum_sparse_n_op += sparse_n_op
+
+            expectation = reference.dot(sparse_n_op.dot(reference))
+
+            if i < self.molecule.n_electrons:
+                assert expectation == 1.0
+            else:
+                assert expectation == 0.0
+
+        convert_after_adding = get_number_preserving_sparse_operator(
+            sum_n_op,
+            self.molecule.n_qubits,
+            self.molecule.n_electrons,
+            spin_preserving=False)
+
+        assert scipy.sparse.linalg.norm(convert_after_adding - sum_sparse_n_op) < 1E-9
+
+        assert reference.dot(sum_sparse_n_op.dot(reference)) - self.molecule.n_electrons < 1E-9
+
+
+    def test_space_size_correct(self):
+        hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian)
+
+        sparse_ham = get_number_preserving_sparse_operator(
+            hamiltonian_fop,
+            self.molecule.n_qubits,
+            self.molecule.n_electrons,
+            spin_preserving=True)
+
+        space_size = sparse_ham.shape[0]
+
+        # Naive Hilbert space size is 2**12, or 4096.
+        assert space_size == 225
+
+
+    def test_hf_energy(self):
+        hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian)
+
+        sparse_ham = get_number_preserving_sparse_operator(
+            hamiltonian_fop,
+            self.molecule.n_qubits,
+            self.molecule.n_electrons,
+            spin_preserving=True)
+
+        space_size = sparse_ham.shape[0]
+        reference = numpy.zeros((space_size))
+        reference[0] = 1.0
+
+        sparse_hf_energy = reference.dot(sparse_ham.dot(reference))
+
+        assert numpy.linalg.norm(sparse_hf_energy - self.molecule.hf_energy) < 1E-9
+
+
+    def test_one_body_hf_energy(self):
+        one_body_part = self.molecular_hamiltonian
+        one_body_part.two_body_tensor = numpy.zeros_like(one_body_part.two_body_tensor)
+
+        one_body_fop = get_fermion_operator(one_body_part)
+        one_body_regular_sparse_op = get_sparse_operator(one_body_fop)
+
+        make_hf_fop = FermionOperator(((3, 1), (2, 1), (1, 1), (0, 1)))
+        make_hf_sparse_op = get_sparse_operator(make_hf_fop, n_qubits=12)
+
+        hf_state = numpy.zeros((2**12))
+        hf_state[0] = 1.0
+        hf_state = make_hf_sparse_op.dot(hf_state)
+
+        regular_sparse_hf_energy = (hf_state.dot(one_body_regular_sparse_op.dot(hf_state))).real
+
+        one_body_sparse_op = get_number_preserving_sparse_operator(
+            one_body_fop,
+            self.molecule.n_qubits,
+            self.molecule.n_electrons,
+            spin_preserving=True)
+
+        space_size = one_body_sparse_op.shape[0]
+        reference = numpy.zeros((space_size))
+        reference[0] = 1.0
+
+        sparse_hf_energy = reference.dot(one_body_sparse_op.dot(reference))
+
+        assert numpy.linalg.norm(sparse_hf_energy - regular_sparse_hf_energy) < 1E-9
+
+
+    def test_ground_state_energy(self):
+        hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian)
+
+        sparse_ham = get_number_preserving_sparse_operator(
+            hamiltonian_fop,
+            self.molecule.n_qubits,
+            self.molecule.n_electrons,
+            spin_preserving=True)
+
+        eig_val, eig_vec = scipy.sparse.linalg.eigsh(sparse_ham, k=1, which='SA')
+
+        assert numpy.abs(eig_val[0] - self.molecule.fci_energy) < 1E-9
+
+
+    def test_doubles_are_subset(self):
+        reference_determinants = [
+            [True, True, True, True, False, False,
+             False, False, False, False, False, False],
+            [True, True, False, False, False, False,
+             True, True, False, False, False, False]
+        ]
+
+        for reference_determinant in reference_determinants:
+            reference_determinant = numpy.asarray(reference_determinant)
+            doubles_state_array = numpy.asarray(
+                list(_iterate_basis_(reference_determinant,
+                                       excitation_level=2,
+                                       spin_preserving=True)))
+            doubles_int_state_array = doubles_state_array.dot(
+                1 << numpy.arange(doubles_state_array.shape[1])[::-1])
+
+            all_state_array = numpy.asarray(
+                list(_iterate_basis_(reference_determinant,
+                                       excitation_level=4,
+                                       spin_preserving=True)))
+            all_int_state_array = all_state_array.dot(
+                1 << numpy.arange(all_state_array.shape[1])[::-1])
+
+            for item in doubles_int_state_array:
+                assert item in all_int_state_array
+
+        for reference_determinant in reference_determinants:
+            reference_determinant = numpy.asarray(reference_determinant)
+            doubles_state_array = numpy.asarray(
+                list(_iterate_basis_(reference_determinant,
+                                       excitation_level=2,
+                                       spin_preserving=True)))
+            doubles_int_state_array = doubles_state_array.dot(
+                1 << numpy.arange(doubles_state_array.shape[1])[::-1])
+
+            all_state_array = numpy.asarray(
+                list(_iterate_basis_(reference_determinant,
+                                       excitation_level=4,
+                                       spin_preserving=False)))
+            all_int_state_array = all_state_array.dot(
+                1 << numpy.arange(all_state_array.shape[1])[::-1])
+
+            for item in doubles_int_state_array:
+                assert item in all_int_state_array
+
+        for reference_determinant in reference_determinants:
+            reference_determinant = numpy.asarray(reference_determinant)
+            doubles_state_array = numpy.asarray(
+                list(_iterate_basis_(reference_determinant,
+                                       excitation_level=2,
+                                       spin_preserving=False)))
+            doubles_int_state_array = doubles_state_array.dot(
+                1 << numpy.arange(doubles_state_array.shape[1])[::-1])
+
+            all_state_array = numpy.asarray(
+                list(_iterate_basis_(reference_determinant,
+                                       excitation_level=4,
+                                       spin_preserving=False)))
+            all_int_state_array = all_state_array.dot(
+                1 << numpy.arange(all_state_array.shape[1])[::-1])
+
+            for item in doubles_int_state_array:
+                assert item in all_int_state_array
+
+
+    def test_full_ham_hermitian(self):
+        hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian)
+
+        sparse_ham = get_number_preserving_sparse_operator(
+            hamiltonian_fop,
+           self.molecule.n_qubits,
+           self.molecule.n_electrons,
+            spin_preserving=True)
+
+        assert scipy.sparse.linalg.norm(sparse_ham - sparse_ham.getH()) < 1E-9
+
+
+    def test_full_ham_hermitian_non_spin_preserving(self):
+        hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian)
+
+        sparse_ham = get_number_preserving_sparse_operator(
+            hamiltonian_fop,
+           self.molecule.n_qubits,
+           self.molecule.n_electrons,
+            spin_preserving=False)
+
+        assert scipy.sparse.linalg.norm(sparse_ham - sparse_ham.getH()) < 1E-9
+
+
+    def test_singles_simple_one_body_term_hermitian(self):
+        fop = FermionOperator(((3, 1), (1, 0)))
+        fop_conj = FermionOperator(((1, 1), (3, 0)))
+
+        sparse_op = get_number_preserving_sparse_operator(
+            fop,
+            self.molecule.n_qubits,
+            self.molecule.n_electrons,
+            spin_preserving=True,
+            excitation_level=1)
+
+        sparse_op_conj = get_number_preserving_sparse_operator(
+            fop_conj,
+            self.molecule.n_qubits,
+            self.molecule.n_electrons,
+            spin_preserving=True,
+            excitation_level=1)
+
+        assert scipy.sparse.linalg.norm(sparse_op - sparse_op_conj.getH()) < 1E-9
+
+
+    def test_singles_simple_two_body_term_hermitian(self):
+        fop = FermionOperator(((3, 1), (8, 1), (1, 0), (4, 0)))
+        fop_conj = FermionOperator(((4, 1), (1, 1), (8, 0), (3, 0)))
+
+        sparse_op = get_number_preserving_sparse_operator(
+            fop,
+            self.molecule.n_qubits,
+            self.molecule.n_electrons,
+            spin_preserving=True,
+            excitation_level=1)
+
+        sparse_op_conj = get_number_preserving_sparse_operator(
+            fop_conj,
+            self.molecule.n_qubits,
+            self.molecule.n_electrons,
+            spin_preserving=True,
+            excitation_level=1)
+
+        assert scipy.sparse.linalg.norm(sparse_op - sparse_op_conj.getH()) < 1E-9
+
+
+    def test_singles_repeating_two_body_term_hermitian(self):
+        fop = FermionOperator(((3, 1), (1, 1), (5, 0), (1, 0)))
+        fop_conj = FermionOperator(((5, 1), (1, 1), (3, 0), (1, 0)))
+
+        sparse_op = get_number_preserving_sparse_operator(
+            fop,
+            self.molecule.n_qubits,
+            self.molecule.n_electrons,
+            spin_preserving=True,
+            excitation_level=1)
+
+        sparse_op_conj = get_number_preserving_sparse_operator(
+            fop_conj,
+            self.molecule.n_qubits,
+            self.molecule.n_electrons,
+            spin_preserving=True,
+            excitation_level=1)
+
+        assert scipy.sparse.linalg.norm(sparse_op - sparse_op_conj.getH()) < 1E-9
+
+
+    def test_singles_ham_hermitian(self):
+        hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian)
+
+        sparse_ham = get_number_preserving_sparse_operator(
+            hamiltonian_fop,
+            self.molecule.n_qubits,
+            self.molecule.n_electrons,
+            spin_preserving=True,
+            excitation_level=1)
+
+        assert scipy.sparse.linalg.norm(sparse_ham - sparse_ham.getH()) < 1E-9
+
+
+    def test_singles_ham_hermitian_non_spin_preserving(self):
+        hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian)
+
+        sparse_ham = get_number_preserving_sparse_operator(
+            hamiltonian_fop,
+            self.molecule.n_qubits,
+            self.molecule.n_electrons,
+            spin_preserving=False,
+            excitation_level=1)
+
+        assert scipy.sparse.linalg.norm(sparse_ham - sparse_ham.getH()) < 1E-9
+
+
+    def test_cisd_energy(self):
+        hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian)
+
+        sparse_ham = get_number_preserving_sparse_operator(
+            hamiltonian_fop,
+            self.molecule.n_qubits,
+            self.molecule.n_electrons,
+            spin_preserving=True,
+            excitation_level=2)
+
+        eig_val, eig_vec = scipy.sparse.linalg.eigsh(sparse_ham, k=1, which='SA')
+
+        assert numpy.abs(eig_val[0] - self.molecule.cisd_energy) < 1E-9
+
+
+    def test_cisd_energy_non_spin_preserving(self):
+        hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian)
+
+        sparse_ham = get_number_preserving_sparse_operator(
+            hamiltonian_fop,
+            self.molecule.n_qubits,
+            self.molecule.n_electrons,
+            spin_preserving=False,
+            excitation_level=2)
+
+        eig_val, eig_vec = scipy.sparse.linalg.eigsh(sparse_ham, k=1, which='SA')
+
+        assert numpy.abs(eig_val[0] - self.molecule.cisd_energy) < 1E-9
+
+    def test_cisd_matches_fci_energy_two_electron_hubbard(self):
+        hamiltonian_fop = self.hubbard_hamiltonian
+
+        sparse_ham_cisd = get_number_preserving_sparse_operator(
+            hamiltonian_fop,
+            8,
+            2,
+            spin_preserving=True,
+            excitation_level=2)
+
+        sparse_ham_fci = get_sparse_operator(
+            hamiltonian_fop,
+            n_qubits=8)
+
+        eig_val_cisd, eig_vec_cisd = scipy.sparse.linalg.eigsh(sparse_ham_cisd, k=1, which='SA')
+        eig_val_fci, eig_vec_fci = scipy.sparse.linalg.eigsh(sparse_ham_fci, k=1, which='SA')
+
+        assert numpy.abs(eig_val_cisd[0] - eig_val_fci[0]) < 1E-9
+
+    def test_weird_determinant_matches_fci_energy_two_electron_hubbard(self):
+        hamiltonian_fop = self.hubbard_hamiltonian
+
+        sparse_ham_cisd = get_number_preserving_sparse_operator(
+            hamiltonian_fop,
+            8,
+            2,
+            spin_preserving=True,
+            excitation_level=2,
+            reference_determinant = numpy.asarray(
+                [False, False, True, True, False, False, False, False]))
+
+        sparse_ham_fci = get_sparse_operator(
+            hamiltonian_fop,
+            n_qubits=8)
+
+        eig_val_cisd, eig_vec_cisd = scipy.sparse.linalg.eigsh(sparse_ham_cisd, k=1, which='SA')
+        eig_val_fci, eig_vec_fci = scipy.sparse.linalg.eigsh(sparse_ham_fci, k=1, which='SA')
+
+        assert numpy.abs(eig_val_cisd[0] - eig_val_fci[0]) < 1E-9
+
+    def test_number_restricted_spectra_match_molecule(self):
+        hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian)
+
+        sparse_ham_number_preserving = get_number_preserving_sparse_operator(
+            hamiltonian_fop,
+            self.molecule.n_qubits,
+            self.molecule.n_electrons,
+            spin_preserving=False)
+
+        sparse_ham = get_sparse_operator(hamiltonian_fop,
+                                         self.molecule.n_qubits)
+
+        sparse_ham_restricted_number_preserving = jw_number_restrict_operator(
+            sparse_ham,
+            n_electrons=self.molecule.n_electrons,
+            n_qubits=self.molecule.n_qubits)
+
+        spectrum_from_new_sparse_method = sparse_eigenspectrum(
+            sparse_ham_number_preserving)
+
+        spectrum_from_old_sparse_method = sparse_eigenspectrum(
+            sparse_ham_restricted_number_preserving)
+
+        spectral_deviation = numpy.amax(numpy.absolute(
+            spectrum_from_new_sparse_method - spectrum_from_old_sparse_method))
+        self.assertAlmostEqual(spectral_deviation, 0.)
+
+    def test_number_restricted_spectra_match_hubbard(self):
+        hamiltonian_fop = self.hubbard_hamiltonian
+
+        sparse_ham_number_preserving = get_number_preserving_sparse_operator(
+            hamiltonian_fop,
+            8,
+            4,
+            spin_preserving=False)
+
+        sparse_ham = get_sparse_operator(hamiltonian_fop,
+                                         8)
+
+        sparse_ham_restricted_number_preserving = jw_number_restrict_operator(
+            sparse_ham,
+            n_electrons=4,
+            n_qubits=8)
+
+        spectrum_from_new_sparse_method = sparse_eigenspectrum(
+            sparse_ham_number_preserving)
+
+        spectrum_from_old_sparse_method = sparse_eigenspectrum(
+            sparse_ham_restricted_number_preserving)
+
+        spectral_deviation = numpy.amax(numpy.absolute(
+            spectrum_from_new_sparse_method - spectrum_from_old_sparse_method))
+        self.assertAlmostEqual(spectral_deviation, 0.)
+
+    def test_number_sz_restricted_spectra_match_molecule(self):
+        hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian)
+
+        sparse_ham_number_sz_preserving = get_number_preserving_sparse_operator(
+            hamiltonian_fop,
+            self.molecule.n_qubits,
+            self.molecule.n_electrons,
+            spin_preserving=True)
+
+        sparse_ham = get_sparse_operator(hamiltonian_fop,
+                                         self.molecule.n_qubits)
+
+        sparse_ham_restricted_number_sz_preserving = jw_sz_restrict_operator(
+            sparse_ham,
+            0,
+            n_electrons=self.molecule.n_electrons,
+            n_qubits=self.molecule.n_qubits)
+
+        spectrum_from_new_sparse_method = sparse_eigenspectrum(
+            sparse_ham_number_sz_preserving)
+
+        spectrum_from_old_sparse_method = sparse_eigenspectrum(
+            sparse_ham_restricted_number_sz_preserving)
+
+        spectral_deviation = numpy.amax(numpy.absolute(
+            spectrum_from_new_sparse_method - spectrum_from_old_sparse_method))
+        self.assertAlmostEqual(spectral_deviation, 0.)
+
+    def test_number_sz_restricted_spectra_match_hubbard(self):
+        hamiltonian_fop = self.hubbard_hamiltonian
+
+        sparse_ham_number_sz_preserving = get_number_preserving_sparse_operator(
+            hamiltonian_fop,
+            8,
+            4,
+            spin_preserving=True)
+
+        sparse_ham = get_sparse_operator(hamiltonian_fop,
+                                         8)
+
+        sparse_ham_restricted_number_sz_preserving = jw_sz_restrict_operator(
+            sparse_ham,
+            0,
+            n_electrons=4,
+            n_qubits=8)
+
+        spectrum_from_new_sparse_method = sparse_eigenspectrum(
+            sparse_ham_number_sz_preserving)
+
+        spectrum_from_old_sparse_method = sparse_eigenspectrum(
+            sparse_ham_restricted_number_sz_preserving)
+
+        spectral_deviation = numpy.amax(numpy.absolute(
+            spectrum_from_new_sparse_method - spectrum_from_old_sparse_method))
+        self.assertAlmostEqual(spectral_deviation, 0.)