Merge pull request #1287 from qiboteam/dbi_fix

Expose DBI operators

src/qibo/hamiltonians/hamiltonians.py

@@ -180,7 +180,7 @@ def energy_fluctuation(self, state):
         energy = self.expectation(state)
         h = self.matrix
         h2 = Hamiltonian(nqubits=self.nqubits, matrix=h @ h, backend=self.backend)
-        average_h2 = self.backend.calculate_expectation_state(h2, state, normalize=True)
+        average_h2 = h2.expectation(state, normalize=True)
         return np.sqrt(np.abs(average_h2 - energy**2))
 
     def __add__(self, o):

src/qibo/models/dbi/double_bracket.py

@@ -56,34 +56,55 @@ def __init__(
     def __call__(
         self, step: float, mode: DoubleBracketGeneratorType = None, d: np.array = None
     ):
+        r"""We use convention that $H' = U^\dagger H U$ where $U=e^{-sW}$ with $W=[D,H]$ (or depending on `mode` an approximation, see `eval_dbr_unitary`). If $s>0$ then for $D = \Delta(H)$ the GWW DBR will give a $\sigma$-decrease, see https://arxiv.org/abs/2206.11772."""
+
+        operator = self.eval_dbr_unitary(step, mode, d)
+        operator_dagger = self.backend.cast(
+            np.matrix(self.backend.to_numpy(operator)).getH()
+        )
+        self.h.matrix = operator_dagger @ self.h.matrix @ operator
+        return operator
+
+    def eval_dbr_unitary(
+        self, step: float, mode: DoubleBracketGeneratorType = None, d: np.array = None
+    ):
+        """In call we will are working in the convention that $H' = U^\\dagger H U$ where $U=e^{-sW}$ with $W=[D,H]$ or an approximation of that by a group commutator. That is handy because if we switch from the DBI in the Heisenberg picture for the Hamiltonian, we get that the transformation of the state is $|\\psi'\rangle = U |\\psi\rangle$ so that $\\langle H\rangle_{\\psi'} = \\langle H' \rangle_\\psi$ (i.e. when writing the unitary acting on the state dagger notation is avoided).
+
+        The group commutator must approximate $U=e^{-s[D,H]}$. This is achieved by setting $r = \\sqrt{s}$ so that
+        $$V = e^{-irH}e^{irD}e^{irH}e^{-irD}$$
+        because
+        $$e^{-irH}De^{irH} = D+ir[D,H]+O(r^2)$$
+        so
+        $$V\approx e^{irD +i^2 r^2[D,H] + O(r^2) -irD} \approx U\\ .$$
+        See the app in https://arxiv.org/abs/2206.11772 for a derivation.
+        """
         if mode is None:
             mode = self.mode
 
         if mode is DoubleBracketGeneratorType.canonical:
             operator = self.backend.calculate_matrix_exp(
-                1.0j * step,
+                -1.0j * step,
                 self.commutator(self.diagonal_h_matrix, self.h.matrix),
             )
         elif mode is DoubleBracketGeneratorType.single_commutator:
             if d is None:
                 d = self.diagonal_h_matrix
             operator = self.backend.calculate_matrix_exp(
-                1.0j * step,
+                -1.0j * step,
                 self.commutator(d, self.h.matrix),
             )
         elif mode is DoubleBracketGeneratorType.group_commutator:
             if d is None:
                 d = self.diagonal_h_matrix
+
+            sqrt_step = np.sqrt(step)
             operator = (
-                self.h.exp(-step)
-                @ self.backend.calculate_matrix_exp(-step, d)
-                @ self.h.exp(step)
-                @ self.backend.calculate_matrix_exp(step, d)
+                self.h.exp(sqrt_step)
+                @ self.backend.calculate_matrix_exp(-sqrt_step, d)
+                @ self.h.exp(-sqrt_step)
+                @ self.backend.calculate_matrix_exp(sqrt_step, d)
             )
-        operator_dagger = self.backend.cast(
-            np.matrix(self.backend.to_numpy(operator)).getH()
-        )
-        self.h.matrix = operator @ self.h.matrix @ operator_dagger
+        return operator
 
     @staticmethod
     def commutator(a, b):

src/qibo/models/dbi/utils.py

@@ -133,11 +133,13 @@ def select_best_dbr_generator(
             )
         else:
             step_best = step
-        dbi_eval(step=step_best)
-        optimal_steps.append(step_best)
-        norms_off_diagonal_restriction.append(dbi_eval.off_diagonal_norm)
+        dbi_object(step=step)
+        optimal_steps.append(step)
+        norms_off_diagonal_restriction.append(dbi_object.off_diagonal_norm)
     # find best d
-    idx_max_loss = np.argmin(norms_off_diagonal_restriction)
+    idx_max_loss = norms_off_diagonal_restriction.index(
+        min(norms_off_diagonal_restriction)
+    )
     flip = flip_list[idx_max_loss]
     step_optimal = optimal_steps[idx_max_loss]
     dbi_eval = deepcopy(dbi_object)

src/qibo/models/utils.py

@@ -207,3 +207,10 @@ def fourier_coefficients(
     # Shift the filtered coefficients back
     filtered_coeffs = np.fft.ifftshift(shifted_filtered_coeffs)
     return filtered_coeffs
+
+
+def vqe_loss(params, circuit, hamiltonian):
+    circuit.set_parameters(params)
+    result = hamiltonian.backend.execute_circuit(circuit)
+    final_state = result.state()
+    return hamiltonian.expectation(final_state)

src/qibo/models/variational.py

@@ -1,8 +1,8 @@
 import numpy as np
 
 from qibo.config import raise_error
-from qibo.hamiltonians import Hamiltonian
 from qibo.models.evolution import StateEvolution
+from qibo.models.utils import vqe_loss
 
 
 class VQE:
@@ -42,6 +42,7 @@ def minimize(
         self,
         initial_state,
         method="Powell",
+        loss_func=None,
         jac=None,
         hess=None,
         hessp=None,
@@ -61,6 +62,7 @@ def minimize(
             method (str): the desired minimization method.
                 See :meth:`qibo.optimizers.optimize` for available optimization
                 methods.
+            loss (callable): loss function, the default one is :func:`qibo.models.utils.vqe_loss`.
             jac (dict): Method for computing the gradient vector for scipy optimizers.
             hess (dict): Method for computing the hessian matrix for scipy optimizers.
             hessp (callable): Hessian of objective function times an arbitrary
@@ -80,25 +82,20 @@ def minimize(
             the ``OptimizeResult``, for ``'cma'`` the ``CMAEvolutionStrategy.result``,
             and for ``'sgd'`` the options used during the optimization.
         """
-
-        def _loss(params, circuit, hamiltonian):
-            circuit.set_parameters(params)
-            result = hamiltonian.backend.execute_circuit(circuit)
-            final_state = result.state()
-            return hamiltonian.expectation(final_state)
-
+        if loss_func is None:
+            loss_func = vqe_loss
         if compile:
-            loss = self.hamiltonian.backend.compile(_loss)
+            loss = self.hamiltonian.backend.compile(loss_func)
         else:
-            loss = _loss
+            loss = loss_func
 
         if method == "cma":
             # TODO: check if we can use this shortcut
             # dtype = getattr(self.hamiltonian.backend.np, self.hamiltonian.backend._dtypes.get('DTYPE'))
             dtype = self.hamiltonian.backend.np.float64
-            loss = lambda p, c, h: dtype(_loss(p, c, h))
+            loss = lambda p, c, h: dtype(loss_func(p, c, h))
         elif method != "sgd":
-            loss = lambda p, c, h: self.hamiltonian.backend.to_numpy(_loss(p, c, h))
+            loss = lambda p, c, h: self.hamiltonian.backend.to_numpy(loss_func(p, c, h))
 
         result, parameters, extra = self.optimizers.optimize(
             loss,
@@ -636,8 +633,6 @@ def minimize(
             The corresponding best parameters.
             extra: variable with historical data for the energy and callbacks.
         """
-        import numpy as np
-
         parameters = np.array([delta_t, 0])
 
         def _loss(params, falqon, hamiltonian):
@@ -647,7 +642,7 @@ def _loss(params, falqon, hamiltonian):
 
         energy = [np.inf]
         callback_result = []
-        for it in range(1, max_layers + 1):
+        for _ in range(1, max_layers + 1):
             beta = self.hamiltonian.backend.to_numpy(
                 _loss(parameters, self, self.evol_hamiltonian)
             )

tests/test_models_dbi.py

@@ -10,7 +10,8 @@
 )
 from qibo.quantum_info import random_hermitian
 
-NSTEPS = 10
+NSTEPS = 50
+SEED = 10
 """Number of steps for evolution."""
 
 
@@ -30,7 +31,7 @@ def test_double_bracket_iteration_canonical(backend, nqubits):
 
 @pytest.mark.parametrize("nqubits", [1, 2])
 def test_double_bracket_iteration_group_commutator(backend, nqubits):
-    h0 = random_hermitian(2**nqubits, backend=backend)
+    h0 = random_hermitian(2**nqubits, backend=backend, seed=SEED)
     d = backend.cast(np.diag(np.diag(backend.to_numpy(h0))))
     dbi = DoubleBracketIteration(
         Hamiltonian(nqubits, h0, backend=backend),
@@ -47,6 +48,29 @@ def test_double_bracket_iteration_group_commutator(backend, nqubits):
     assert initial_off_diagonal_norm > dbi.off_diagonal_norm
 
 
+@pytest.mark.parametrize("nqubits", [3])
+def test_double_bracket_iteration_eval_dbr_unitary(backend, nqubits):
+    r"""The bound is $$||e^{-[D,H]}-GC||\le s^{3/2}(||[H,[D,H]||+||[D,[D,H]]||$$"""
+    h0 = random_hermitian(2**nqubits, backend=backend)
+    d = backend.cast(np.diag(np.diag(backend.to_numpy(h0))))
+    dbi = DoubleBracketIteration(
+        Hamiltonian(nqubits, h0, backend=backend),
+        mode=DoubleBracketGeneratorType.group_commutator,
+    )
+
+    for s in np.linspace(0.001, 0.01, NSTEPS):
+        u = dbi.eval_dbr_unitary(
+            s, d=d, mode=DoubleBracketGeneratorType.single_commutator
+        )
+        v = dbi.eval_dbr_unitary(
+            s, d=d, mode=DoubleBracketGeneratorType.group_commutator
+        )
+
+        assert np.linalg.norm(u - v) < 10 * s**1.49 * (
+            np.linalg.norm(h0) + np.linalg.norm(d)
+        ) * np.linalg.norm(h0) * np.linalg.norm(d)
+
+
 @pytest.mark.parametrize("nqubits", [1, 2])
 def test_double_bracket_iteration_single_commutator(backend, nqubits):
     h0 = random_hermitian(2**nqubits, backend=backend)

tests/test_models_dbi_utils.py

@@ -32,7 +32,7 @@ def test_generate_Z_operators(backend, nqubits):
 
 
 @pytest.mark.parametrize("nqubits", [1, 2])
-@pytest.mark.parametrize("step", [0.1, None])
+@pytest.mark.parametrize("step", [0.1, 0.2])
 def test_select_best_dbr_generator(backend, nqubits, step):
     h0 = random_hermitian(2**nqubits, seed=1, backend=backend)
     dbi = DoubleBracketIteration(