Implement double bracket flow

src/qibo/hamiltonians/double_bracket.py

@@ -0,0 +1,138 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import numpy.typing as npt
+
+import qibo
+from qibo.hamiltonians.hamiltonians import Hamiltonian
+
+FIGSIZE = (5, 5)
+
+
+class DoubleBracket(Hamiltonian):
+    def __init__(self, nqubits: int, matrix: npt.DTypeLike = None, backend=None):
+        """
+        Class implementing the Double Bracket flow algorithm.
+        For more details, see https://arxiv.org/pdf/2206.11772.pdf
+
+        Args:
+            nqubits(int): Number of qubits.
+            matrix(npt.DTypeLike, Optional): Hamiltonian matrix.
+            backend(Optional): Qibo backend.
+
+        Example:
+            .. code-block:: python
+
+            import numpy as np
+            from qibo.hamiltonians.double_bracket import DoubleBracket
+            from qibo.hamiltonians.hamiltonians import Hamiltonian
+
+            nqubits = 3
+
+            example_matrix = [[ 3.,  0.,  0.,  1.,  0.,  0.,  1.,  0.],
+                   [ 0.,  1.,  1.,  0.,  0.,  0.,  0.,  1.],
+                   [ 0.,  1.,  1.,  0.,  1.,  0.,  0.,  0.],
+                   [ 1.,  0.,  0., -1.,  0.,  1.,  0.,  0.],
+                   [ 0.,  0.,  1.,  0.,  1.,  0.,  0.,  1.],
+                   [ 0.,  0.,  0.,  1.,  0., -1.,  1.,  0.],
+                   [ 1.,  0.,  0.,  0.,  0.,  1., -1.,  0.],
+                   [ 0.,  1.,  0.,  0.,  1.,  0.,  0., -3.]]
+
+            example_matrix = np.array(example_matrix)
+            ham = Hamiltonian(nqubits, example_matrix) # implement an hamiltonian in hamiltonians
+            db = DoubleBracket(nqubits, matrix = ham.matrix)
+
+            for _ in range(1000):
+                db.evolve( db.group_commutator( 1, Hamiltonian(nqubits,db.delta)))
+
+
+            db.plot_abs_matrix()
+            db.plot_norm()
+        """
+        super().__init__(nqubits, matrix, backend)
+        self.history = []
+        self.norm = []
+
+    @property
+    def delta(self):
+        """
+        Hamiltonian diagonal restriction.
+        """
+        return self.backend.cast(np.diag(self.matrix.diagonal()))
+
+    @property
+    def sigma(self):
+        """
+        Hamiltonian off-diagonal restriction.
+        """
+        return self.matrix - self.delta
+
+    def db_generator(self, step):
+        """Generator of the GWW flow."""
+        w = Hamiltonian(
+            nqubits=self.nqubits,
+            matrix=self.delta @ self.sigma - self.sigma @ self.delta,
+        )
+        return w.exp(step)
+
+    def group_commutator(self, step: float, gen: Hamiltonian):
+        """
+        Group commutator.
+
+        Args:
+            step (float)
+            gen (Hamiltonian): Flow generator.
+        """
+        return gen.exp(-1 * step) @ self.exp(-1 * step) @ gen.exp(step) @ self.exp(step)
+
+    def evolve(self, evolution_op: Hamiltonian):
+        """
+        Evolution of the Hamiltonian according to the group commutator generated by `gen`.
+
+        Args:
+            gen (Hamiltonian): generator of the group commutator.
+        """
+        self.history.append(
+            self.matrix
+        )  # TODO: Delete it, if we don't need to save the matrix evolution
+        self.matrix = (
+            evolution_op
+            @ self.matrix
+            @ self.backend.cast(np.transpose(np.conjugate(evolution_op)))
+        )
+        self.norm.append(np.linalg.norm(self.sigma))
+
+    def plot_abs_matrix(self, path: str = None):
+        """
+        Plot the Hamiltonian absolute values.
+
+        Args:
+            path (str, Optional): if a path is provided, the plot is dumped.
+
+        """
+        fig, ax = plt.subplots(figsize=FIGSIZE)
+        plot = plt.imshow(np.absolute(self.matrix), cmap="RdBu")
+        plt.xticks([])
+        plt.yticks([])
+        fig.colorbar(plot, ax=ax, extend="both")
+        if path:
+            plt.savefig(path)
+        else:
+            plt.show()
+
+    def plot_norm(self, path: str = None):
+        """
+
+        Plot the evolution of the Hamiltonian off diagonal norm.
+
+        Args:
+            path (str, Optional): if a path is provided, the plot is dumped.
+
+        """
+
+        plt.plot(self.norm, "-o")
+        plt.title(r"Norm off-diagonal $\vert\vert\sigma(H_k)\vert\vert$")
+        plt.xlabel("epochs")
+        if path:
+            plt.savefig(path)
+        else:
+            plt.show()