TETRIS-ADAPT-VQE

tangelo/algorithms/variational/__init__.py

@@ -18,3 +18,4 @@
 from .sa_oo_vqe_solver import SA_OO_Solver
 from .iqcc_solver import iQCC_solver
 from .iqcc_ilc_solver import iQCC_ILC_solver
+from .tetris_adapt_vqe_solver import TETRISADAPTSolver

tangelo/algorithms/variational/adapt_vqe_solver.py

@@ -166,6 +166,9 @@ def build(self):
                             }
 
         self.vqe_solver = VQESolver(self.vqe_options)
+
+        # Circuits without variational parameter raise a warning in VQESolver.
+        warnings.filterwarnings("ignore", message="No variational gate found in the circuit.")
         self.vqe_solver.build()
 
         # If applicable, give vqe_solver access to molecule object
@@ -229,21 +232,21 @@ def simulate(self):
 
             full_circuit = (self.vqe_solver.ansatz.circuit if self.ref_state is None else
                             self.vqe_solver.reference_circuit + self.vqe_solver.ansatz.circuit)
-            pool_select = self.rank_pool(self.pool_commutators, full_circuit,
-                                         backend=self.vqe_solver.backend, tolerance=self.tol)
+            gradients = self.compute_gradients(full_circuit, backend=self.vqe_solver.backend)
+            pool_select = self.choose_operator(gradients, tolerance=self.tol)
 
             # If pool selection returns an operator that changes the energy by
             # more than self.tol. Else, the loop is complete and the energy is
             # considered as converged.
-            if pool_select > -1:
-
+            if pool_select:
                 # Adding a new operator + initializing its parameters to 0.
                 # Previous parameters are kept as they were.
-                params += [0.]
-                if self.pool_type == 'fermion':
-                    self.vqe_solver.ansatz.add_operator(self.pool_operators[pool_select], self.fermionic_operators[pool_select])
-                else:
-                    self.vqe_solver.ansatz.add_operator(self.pool_operators[pool_select])
+                params += [0.] * len(pool_select)
+                for ps in pool_select:
+                    if self.pool_type == 'fermion':
+                        self.vqe_solver.ansatz.add_operator(self.pool_operators[ps], self.fermionic_operators[ps])
+                    else:
+                        self.vqe_solver.ansatz.add_operator(self.pool_operators[ps])
                 self.vqe_solver.initial_var_params = params
 
                 # Performs a VQE simulation and append the energy to a list.
@@ -258,24 +261,26 @@ def simulate(self):
                 self.converged = True
                 break
 
+        # Reconstructing the optimal circuit at the end of the ADAPT iterations
+        # or when the algorithm has converged.
+        self.ansatz.build_circuit(self.optimal_var_params)
+        self.optimal_circuit = (self.vqe_solver.ansatz.circuit if self.ref_state is None else
+                                self.vqe_solver.reference_circuit + self.vqe_solver.ansatz.circuit)
+
         return self.energies[-1]
 
-    def rank_pool(self, pool_commutators, circuit, backend, tolerance=1e-3):
-        """Rank pool of operators with a specific circuit.
+    def compute_gradients(self, circuit, backend):
+        """Compute gradients for the operator pool with a specific circuit.
 
         Args:
-            pool_commutators (QubitOperator): Commutator [H, operator] for each
-                generator.
             circuit (tangelo.linq.Circuit): Circuit for measuring each commutator.
             backend (tangelo.linq.Simulator): Backend to compute expectation values.
-            tolerance (float): Minimum value for gradient to be considered.
 
-        Returns:
-            int: Index of the operators with the highest gradient. If it is not
-                bigger than tolerance, returns -1.
+       Returns:
+            list of float: Operator gradients.
         """
 
-        gradient = [abs(backend.get_expectation_value(element, circuit)) for element in pool_commutators]
+        gradient = [abs(backend.get_expectation_value(element, circuit)) for element in self.pool_commutators]
         for deflate_circuit in self.deflation_circuits:
             for i, pool_op in enumerate(self.pool_operators):
                 op_circuit = Circuit([Gate(op[1], op[0]) for tuple in pool_op.terms for op in tuple])
@@ -285,12 +290,29 @@ def rank_pool(self, pool_commutators, circuit, backend, tolerance=1e-3):
                 pool_over = deflate_circuit.inverse() + circuit
                 f_dict, _ = backend.simulate(pool_over)
                 gradient[i] += self.deflation_coeff * grad * f_dict.get("0"*self.vqe_solver.ansatz.circuit.width, 0)
-        max_partial = max(gradient)
+
+        return gradient
+
+    def choose_operator(self, gradients, tolerance=1e-3):
+        """Choose next operator to add according to the ADAPT-VQE algorithm.
+
+        Args:
+            gradients (list of float): Operator gradients corresponding to
+                self.pool_operators.
+            tolerance (float): Minimum value for gradient to be considered.
+
+        Returns:
+            list of int: Index (list of length=1) of the operator with the
+                highest gradient. If it is not bigger than tolerance, returns
+                an empty list.
+        """
+        sorted_op_indices = sorted(range(len(gradients)), key=lambda k: gradients[k])
+        max_partial = gradients[sorted_op_indices[-1]]
 
         if self.verbose:
-            print(f"LARGEST PARTIAL DERIVATIVE: {max_partial :4E} \t[{gradient.index(max_partial)}]")
+            print(f"LARGEST PARTIAL DERIVATIVE: {max_partial :4E}")
 
-        return gradient.index(max_partial) if max_partial >= tolerance else -1
+        return [sorted_op_indices[-1]] if max_partial >= tolerance else []
 
     def get_resources(self):
         """Returns resources currently used in underlying VQE."""
@@ -306,13 +328,6 @@ def LBFGSB_optimizer(self, func, var_params):
         self.optimal_var_params = result.x
         self.optimal_energy = result.fun
 
-        # Reconstructing the optimal circuit at the end of the ADAPT iterations
-        # or when the algorithm has converged.
-        if self.converged or self.iteration == self.max_cycles:
-            self.ansatz.build_circuit(self.optimal_var_params)
-            self.optimal_circuit = (self.vqe_solver.ansatz.circuit if self.ref_state is None else
-                                    self.vqe_solver.reference_circuit + self.vqe_solver.ansatz.circuit)
-
         if self.verbose:
             print(f"VQESolver optimization results:")
             print(f"\tOptimal VQE energy: {result.fun}")

tangelo/algorithms/variational/tests/test_tetris_adapt_vqe_solver.py

@@ -0,0 +1,42 @@
+# Copyright 2021 Good Chemistry Company.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#   http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import unittest
+
+from tangelo.algorithms.variational import TETRISADAPTSolver
+from tangelo.molecule_library import mol_H2_sto3g
+
+
+class TETRISADAPTSolverTest(unittest.TestCase):
+
+    def test_build_tetris_adapt(self):
+        """Try instantiating TETRISADAPTSolver with basic input."""
+
+        opt_dict = {"molecule": mol_H2_sto3g, "max_cycles": 15}
+        adapt_solver = TETRISADAPTSolver(opt_dict)
+        adapt_solver.build()
+
+    def test_single_cycle_tetris_adapt(self):
+        """Try instantiating TETRISADAPTSolver with basic input."""
+
+        opt_dict = {"molecule": mol_H2_sto3g, "max_cycles": 1, "verbose": False}
+        adapt_solver = TETRISADAPTSolver(opt_dict)
+        adapt_solver.build()
+        adapt_solver.simulate()
+
+        self.assertAlmostEqual(adapt_solver.optimal_energy, -1.13727, places=4)
+
+
+if __name__ == "__main__":
+    unittest.main()

tangelo/algorithms/variational/tetris_adapt_vqe_solver.py

@@ -0,0 +1,74 @@
+# Copyright 2021 Good Chemistry Company.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#   http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""Module that defines the TETRIS-ADAPT-VQE algorithm framework. Differences vs
+ADAPT-VQE: more than one operators acting on different qubits can be added to
+the adaptive ansatz during the same ADAPT cycle. This algorithm creates denser
+circuits.
+
+Ref:
+    Panagiotis G. Anastasiou, Yanzhu Chen, Nicholas J. Mayhall, Edwin Barnes,
+    and Sophia E. Economou.
+    TETRIS-ADAPT-VQE: An adaptive algorithm that yields shallower, denser
+    circuit ansätze
+    arXiv:2209.10562 (2022)
+"""
+
+from tangelo.algorithms.variational import ADAPTSolver
+
+
+class TETRISADAPTSolver(ADAPTSolver):
+    """TETRIS-ADAPT-VQE class. This is an iterative algorithm that uses VQE. A
+    single method is redefined from ADAPTSolver to allow the addition of many
+    operators per ADAPT cycle.
+    """
+
+    def choose_operator(self, gradients, tolerance=1e-3):
+        """Choose the next operator(s) to add according to the TETRIS-ADAPT-VQE
+        algorithm.
+
+        Args:
+            gradients (list of float): Operator gradients (absolute values)
+                corresponding to self.pool_operators.
+            tolerance (float): Minimum value for gradient to be considered.
+
+        Returns:
+            list of int: Indice(s) of the operator(s) to be considered for this
+                ADAPT cycle.
+        """
+
+        qubit_indices = set(range(self.ansatz.circuit.width))
+
+        # Sorting the pool operators according to the gradients.
+        sorted_op_indices = sorted(range(len(gradients)), key=lambda k: gradients[k])
+
+        op_indices_to_add = list()
+        for i in sorted_op_indices[::-1]:
+
+            # If gradient is lower than the tolerance, all the remaining
+            # operators have a lower gradient also. If there is no "available"
+            # qubit anymore, no more operator can be added.
+            if gradients[i] < tolerance or len(qubit_indices) == 0:
+                break
+
+            op_indices = self.pool_operators[i].qubit_indices
+
+            # If the operator acts on a subset of "available" qubits, it can be
+            # considered for this ADAPT cycle. Those qubit indices are then
+            # removed from consideration for this ADAPT cycle.
+            if op_indices.issubset(qubit_indices):
+                qubit_indices -= op_indices
+                op_indices_to_add.append(i)
+
+        return op_indices_to_add

tangelo/toolboxes/operators/operators.py

@@ -124,6 +124,23 @@ def get_max_number_hamiltonian_terms(self, n_qubits):
 
         return sum([comb(n_qubits, 2*i, exact=True) * 3**(n_qubits-2*i) for i in range(n_qubits//2)])
 
+    @property
+    def qubit_indices(self):
+        """Return a set of integers corresponding to qubit indices the qubit
+        operator acts on.
+
+        Returns:
+            set: Set of qubit indices.
+        """
+
+        qubit_indices = set()
+        for term in self.terms:
+            if term:
+                indices = list(zip(*term))[0]
+                qubit_indices.update(indices)
+
+        return qubit_indices
+
 
 class QubitHamiltonian(QubitOperator):
     """QubitHamiltonian objects are essentially openfermion.QubitOperator

tangelo/toolboxes/operators/tests/test_operators.py

@@ -16,7 +16,8 @@
 
 import numpy as np
 
-from tangelo.toolboxes.operators import QubitHamiltonian, FermionOperator, QubitOperator, count_qubits, qubitop_to_qubitham
+from tangelo.toolboxes.operators import QubitHamiltonian, FermionOperator, \
+    QubitOperator, count_qubits, qubitop_to_qubitham
 
 
 class OperatorsUtilitiesTest(unittest.TestCase):
@@ -74,6 +75,21 @@ def test_get_coeff(self):
         np.testing.assert_array_almost_equal(ref_two_body, two_body)
 
 
+class QubitOperatorTest(unittest.TestCase):
+
+    def test_qubit_indices(self):
+        """Test the qubit_indices property of the QubitOperator class."""
+
+        q1 = QubitOperator()
+        self.assertEqual(q1.qubit_indices, set())
+
+        q2 = QubitOperator("Z0 Z1")
+        self.assertEqual(q2.qubit_indices, {0, 1})
+
+        q3 = QubitOperator("Z0 Z1") + QubitOperator("Z1 Z2")
+        self.assertEqual(q3.qubit_indices, {0, 1, 2})
+
+
 class QubitHamiltonianTest(unittest.TestCase):
 
     def test_instantiate_QubitHamiltonian(self):