DBI 3rd order GCI

src/qibo/models/dbi/double_bracket.py

@@ -23,7 +23,12 @@ class DoubleBracketGeneratorType(Enum):
     """Use single commutator."""
     group_commutator = auto()
     """Use group commutator approximation"""
-    # TODO: add double commutator (does it converge?)
+    group_commutator_3rd_order = auto()
+    """Implements: $e^{\frac{\\sqrt{5}-1}{2}isH}e^{\frac{\\sqrt{5}-1}{2}isD}e^{-isH}e^{isD}e^{\frac{3-\\sqrt{5}}{2}isH}e^{isD}
+    \approx e^{-s^2[H,D]} + O(s^4)$
+    which is equation (8) in https://arxiv.org/abs/2111.12177]
+    s must be taken as $\\sqrt{s}$ to approximate the flow using the commutator
+    """
 
 
 class DoubleBracketCostFunction(str, Enum):
@@ -125,6 +130,17 @@ def __call__(
                 @ self.h.exp(step)
                 @ self.backend.calculate_matrix_exp(step, d)
             )
+        elif mode is DoubleBracketGeneratorType.group_commutator_3rd_order:
+            if d is None:
+                d = self.diagonal_h_matrix
+            operator = (
+                self.h.exp(-step * (np.sqrt(5) - 1) / 2)
+                @ self.backend.calculate_matrix_exp(-step * (np.sqrt(5) - 1) / 2, d)
+                @ self.h.exp(step)
+                @ self.backend.calculate_matrix_exp(step * (np.sqrt(5) + 1) / 2, d)
+                @ self.h.exp(-step * (3 - np.sqrt(5)) / 2)
+                @ self.backend.calculate_matrix_exp(-step, d)
+            )
         operator_dagger = self.backend.cast(
             np.array(np.matrix(self.backend.to_numpy(operator)).getH())
         )

tests/test_models_dbi.py

@@ -58,6 +58,25 @@ def test_double_bracket_iteration_group_commutator(backend, nqubits):
     assert initial_off_diagonal_norm > dbi.off_diagonal_norm
 
 
+@pytest.mark.parametrize("nqubits", [1, 2])
+def test_double_bracket_iteration_group_commutator_3(backend, nqubits):
+    """Check 3rd order group commutator mode."""
+    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),
+        mode=DoubleBracketGeneratorType.group_commutator_3rd_order,
+    )
+    initial_off_diagonal_norm = dbi.off_diagonal_norm
+
+    # test first iteration with default d
+    dbi(mode=DoubleBracketGeneratorType.group_commutator_3, step=0.01)
+    for _ in range(NSTEPS):
+        dbi(step=0.01, d=d)
+
+    assert initial_off_diagonal_norm > dbi.off_diagonal_norm
+
+
 @pytest.mark.parametrize("nqubits", [1, 2])
 def test_double_bracket_iteration_single_commutator(backend, nqubits):
     """Check single commutator mode."""