Automatically detect spin symmetry in Bogoliubov transformation matrix (

#463)

* index negative energies

* block diag

* increase numpy requirement for block

* fix initial state

* deprecate orbital_energies

* fix ground_energy

requirements.txt

@@ -1,4 +1,4 @@
-numpy>=1.11.0
+numpy>=1.13.0
 scipy>=1.1.0
 h5py>=2.8
 future

src/openfermion/ops/_quadratic_hamiltonian.py

@@ -12,6 +12,8 @@
 
 """Hamiltonians that are quadratic in the fermionic ladder operators."""
 
+import warnings
+
 import numpy
 from scipy.linalg import schur
 
@@ -125,52 +127,12 @@ def add_chemical_potential(self, chemical_potential):
                                       numpy.eye(self.n_qubits))
         self.chemical_potential += chemical_potential
 
-    def orbital_energies(self, non_negative=False):
-        r"""Return the orbital energies.
-
-        Any quadratic Hamiltonian is unitarily equivalent to a Hamiltonian
-        of the form
-
-        .. math::
-
-            \sum_{j} \varepsilon_j b^\dagger_j b_j + \text{constant}.
-
-        We call the :math:`\varepsilon_j` the orbital energies.
-        The eigenvalues of the Hamiltonian are sums of subsets of the
-        orbital energies (up to the additive constant).
-
-        Args:
-            non_negative(bool): If True, always return a list of orbital
-                energies that are non-negative. This option is ignored if
-                the Hamiltonian does not conserve particle number, in which
-                case the returned orbital energies are always non-negative.
-
-        Returns
-        -------
-        orbital_energies(ndarray)
-            A one-dimensional array containing the :math:`\varepsilon_j`
-        constant(float)
-            The constant
-        """
-        if self.conserves_particle_number and not non_negative:
-            hermitian_matrix = self.combined_hermitian_part
-            orbital_energies, _ = numpy.linalg.eigh(
-                hermitian_matrix)
-            constant = self.constant
-        else:
-            majorana_matrix, majorana_constant = self.majorana_form()
-            canonical, _ = antisymmetric_canonical_form(
-                majorana_matrix)
-            orbital_energies = canonical[
-                range(self.n_qubits), range(self.n_qubits, 2 * self.n_qubits)]
-            constant = -0.5 * numpy.sum(orbital_energies) + majorana_constant
-
-        return orbital_energies, constant
-
     def ground_energy(self):
         """Return the ground energy."""
-        _, constant = self.orbital_energies(non_negative=True)
-        return constant
+        orbital_energies, _, constant = (
+                self.diagonalizing_bogoliubov_transform())
+        return numpy.sum(orbital_energies[
+            numpy.where(orbital_energies < 0.0)[0]]) + constant
 
     def majorana_form(self):
         r"""Return the Majorana represention of the Hamiltonian.
@@ -324,10 +286,30 @@ def index_map(i):
             spin_indices = [index_map(i) for i in range(n_sites)]
             matrix = self.combined_hermitian_part[
                     numpy.ix_(spin_indices, spin_indices)]
+            orbital_energies, diagonalizing_unitary_T = numpy.linalg.eigh(
+                    matrix)
         else:
             matrix = self.combined_hermitian_part
 
-        orbital_energies, diagonalizing_unitary_T = numpy.linalg.eigh(matrix)
+            if _is_spin_block_diagonal(matrix):
+                up_block = matrix[:n_modes//2, :n_modes//2]
+                down_block = matrix[n_modes//2:, n_modes//2:]
+
+                up_orbital_energies, up_diagonalizing_unitary_T = (
+                        numpy.linalg.eigh(up_block))
+                down_orbital_energies, down_diagonalizing_unitary_T = (
+                        numpy.linalg.eigh(down_block))
+
+                orbital_energies = numpy.concatenate(
+                    (up_orbital_energies, down_orbital_energies))
+                diagonalizing_unitary_T = numpy.block([
+                    [up_diagonalizing_unitary_T,
+                        numpy.zeros((n_modes//2, n_modes//2))],
+                    [numpy.zeros((n_modes//2, n_modes//2)),
+                        down_diagonalizing_unitary_T]])
+            else:
+                orbital_energies, diagonalizing_unitary_T = numpy.linalg.eigh(
+                        matrix)
 
         return orbital_energies, diagonalizing_unitary_T.T, self.constant
 
@@ -417,6 +399,44 @@ def diagonalizing_circuit(self):
 
         return circuit_description
 
+    # DEPRECATED FUNCTIONS
+    # ====================
+    def orbital_energies(self, non_negative=False):
+        r"""Return the orbital energies.
+
+        Any quadratic Hamiltonian is unitarily equivalent to a Hamiltonian
+        of the form
+
+        .. math::
+
+            \sum_{j} \varepsilon_j b^\dagger_j b_j + \text{constant}.
+
+        We call the :math:`\varepsilon_j` the orbital energies.
+        The eigenvalues of the Hamiltonian are sums of subsets of the
+        orbital energies (up to the additive constant).
+
+        Args:
+            non_negative(bool): If True, always return a list of orbital
+                energies that are non-negative. This option is ignored if
+                the Hamiltonian does not conserve particle number, in which
+                case the returned orbital energies are always non-negative.
+
+        Returns
+        -------
+        orbital_energies(ndarray)
+            A one-dimensional array containing the :math:`\varepsilon_j`
+        constant(float)
+            The constant
+        """
+        warnings.warn('The method `orbital_energies` is deprecated. '
+                      'Use the method `diagonalizing_bogoliubov_transform` '
+                      'instead.', DeprecationWarning)
+
+        orbital_energies, _, constant = (
+                self.diagonalizing_bogoliubov_transform())
+
+        return orbital_energies, constant
+
 
 def antisymmetric_canonical_form(antisymmetric_matrix):
     """Compute the canonical form of an antisymmetric matrix.
@@ -495,3 +515,13 @@ def antisymmetric_canonical_form(antisymmetric_matrix):
             swap_rows(diagonal, i, arg_min)
 
     return canonical, orthogonal.T
+
+
+def _is_spin_block_diagonal(matrix) -> bool:
+    n = matrix.shape[0]
+    if n % 2:
+        return False
+    max_upper_right = numpy.max(numpy.abs(matrix[:n//2, n//2:]))
+    max_lower_left = numpy.max(numpy.abs(matrix[n//2:, :n//2]))
+    return (numpy.isclose(max_upper_right, 0.0)
+            and numpy.isclose(max_lower_left, 0.0))

src/openfermion/ops/_quadratic_hamiltonian_test.py

@@ -239,6 +239,29 @@ def test_diagonalizing_bogoliubov_transform_particle_conserving(self):
         quad_ham = get_quadratic_hamiltonian(
                 reorder(get_fermion_operator(quad_ham), up_then_down))
 
+        orbital_energies, transformation_matrix, _ = (
+                quad_ham.diagonalizing_bogoliubov_transform()
+        )
+        max_upper_right = numpy.max(numpy.abs(transformation_matrix[:5, 5:]))
+        max_lower_left = numpy.max(numpy.abs(transformation_matrix[5:, :5]))
+
+        numpy.testing.assert_allclose(
+                orbital_energies[:5], orbital_energies[5:])
+        numpy.testing.assert_allclose(
+                transformation_matrix.dot(
+                    quad_ham.combined_hermitian_part.T.dot(
+                        transformation_matrix.T.conj())),
+                numpy.diag(orbital_energies),
+                atol=1e-7)
+        numpy.testing.assert_allclose(max_upper_right, 0.0)
+        numpy.testing.assert_allclose(max_lower_left, 0.0)
+
+        # Specific spin sector
+        quad_ham = random_quadratic_hamiltonian(
+                5, conserves_particle_number=True, expand_spin=True)
+        quad_ham = get_quadratic_hamiltonian(
+                reorder(get_fermion_operator(quad_ham), up_then_down))
+
         for spin_sector in range(2):
             orbital_energies, transformation_matrix, _ = (
                     quad_ham.diagonalizing_bogoliubov_transform(

src/openfermion/utils/_slater_determinants.py

@@ -105,8 +105,7 @@ def gaussian_state_preparation_circuit(
         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(orbital_energies < 0.0)
-            occupied_orbitals = range(num_negative_energies)
+            occupied_orbitals = numpy.where(orbital_energies < 0.0)[0]
 
         # Get the unitary rows which represent the Slater determinant
         slater_determinant_matrix = transformation_matrix[occupied_orbitals]