Merge pull request #1487 from qiboteam/qfim

Add `quantum_fisher_information_matrix` to `quantum_info.metrics`

doc/source/api-reference/qibo.rst

@@ -1951,6 +1951,12 @@ Frame Potential
 .. autofunction:: qibo.quantum_info.frame_potential
 
 
+Quantum Fisher Information Matrix
+"""""""""""""""""""""""""""""""""
+
+.. autofunction:: qibo.quantum_info.quantum_fisher_information_matrix
+
+
 Linear Algebra Operations
 ^^^^^^^^^^^^^^^^^^^^^^^^^
 

src/qibo/backends/abstract.py

@@ -376,6 +376,13 @@ def calculate_singular_value_decomposition(self, matrix):  # pragma: no cover
         """Calculate the Singular Value Decomposition of ``matrix``."""
         raise_error(NotImplementedError)
 
+    @abc.abstractmethod
+    def calculate_jacobian_matrix(
+        self, circuit, parameters, initial_state=None, return_complex: bool = True
+    ):  # pragma: no cover
+        """Calculate the Jacobian matrix of ``circuit`` with respect to varables ``params``."""
+        raise_error(NotImplementedError)
+
     @abc.abstractmethod
     def calculate_hamiltonian_matrix_product(
         self, matrix1, matrix2

src/qibo/backends/numpy.py

@@ -799,6 +799,15 @@ def calculate_matrix_power(
     def calculate_singular_value_decomposition(self, matrix):
         return self.np.linalg.svd(matrix)
 
+    def calculate_jacobian_matrix(
+        self, circuit, parameters=None, initial_state=None, return_complex: bool = True
+    ):
+        raise_error(
+            NotImplementedError,
+            "This method is only implemented in backends that allow automatic differentiation, "
+            + "e.g. ``PytorchBackend`` and ``TensorflowBackend``.",
+        )
+
     # TODO: remove this method
     def calculate_hamiltonian_matrix_product(self, matrix1, matrix2):
         return matrix1 @ matrix2

src/qibo/backends/pytorch.py

@@ -239,6 +239,21 @@ def calculate_matrix_power(
         copied = super().calculate_matrix_power(copied, power, precision_singularity)
         return self.cast(copied, dtype=copied.dtype)
 
+    def calculate_jacobian_matrix(
+        self, circuit, parameters=None, initial_state=None, return_complex: bool = True
+    ):
+        copied = circuit.copy(deep=True)
+
+        def func(parameters):
+            """torch requires object(s) to be wrapped in a function."""
+            copied.set_parameters(parameters)
+            state = self.execute_circuit(copied, initial_state=initial_state).state()
+            if return_complex:
+                return self.np.real(state), self.np.imag(state)
+            return self.np.real(state)
+
+        return self.np.autograd.functional.jacobian(func, parameters)
+
     def _test_regressions(self, name):
         if name == "test_measurementresult_apply_bitflips":
             return [

src/qibo/backends/tensorflow.py

@@ -220,6 +220,26 @@ def calculate_singular_value_decomposition(self, matrix):
         S, U, V = self.tf.linalg.svd(matrix)
         return U, S, self.np.conj(self.np.transpose(V))
 
+    def calculate_jacobian_matrix(
+        self, circuit, parameters=None, initial_state=None, return_complex: bool = True
+    ):
+        copied = circuit.copy(deep=True)
+
+        # necessary for the tape to properly watch the variables
+        parameters = self.tf.Variable(parameters)
+
+        with self.tf.GradientTape(persistent=return_complex) as tape:
+            copied.set_parameters(parameters)
+            state = self.execute_circuit(copied, initial_state=initial_state).state()
+            real = self.np.real(state)
+            if return_complex:
+                imag = self.np.imag(state)
+
+        if return_complex:
+            return tape.jacobian(real, parameters), tape.jacobian(imag, parameters)
+
+        return tape.jacobian(real, parameters)
+
     def calculate_hamiltonian_matrix_product(self, matrix1, matrix2):
         if self.is_sparse(matrix1) or self.is_sparse(matrix2):
             raise_error(

src/qibo/quantum_info/metrics.py

@@ -845,6 +845,76 @@ def frame_potential(
     return potential / samples**2
 
 
+def quantum_fisher_information_matrix(
+    circuit,
+    parameters=None,
+    initial_state=None,
+    return_complex: bool = True,
+    backend=None,
+):
+    """Calculate the Quantum Fisher Information Matrix (QFIM) of a parametrized ``circuit``.
+
+    Given a set of ``parameters`` :math:`\\theta = \\{\\theta_{k}\\}_{k\\in[M]}` and a
+    parameterized unitary ``circuit`` :math:`U(\\theta)` acting on an ``initial_state``
+    :math:`\\ket{\\phi}`, the QFIM is such that its elements can be calculated as
+
+    .. math::
+        \\mathbf{F}_{jk} = 4 \\, \\text{Re}\\left\\{ \\braket{\\partial_{j} \\psi | \\partial_{k}
+            \\psi} - \\braket{\\partial_{j} \\psi | \\psi}\\!\\braket{\\psi | \\partial_{k} \\psi}
+            \\right\\} \\, ,
+
+    where we have used the short notations :math:`\\ket{\\psi} \\equiv \\ket{\\psi(\\theta)}
+    = U(\\theta) \\ket{\\phi}`, and :math:`\\ket{\\partial_{k} \\psi} \\equiv \\frac{\\partial}
+    {\\partial\\theta_{k}} \\ket{\\psi(\\theta)}`.
+    If the ``initial_state`` :math:`\\ket{\\phi}` is not specified, it defaults to
+    :math:`\\ket{0}^{\\otimes n}`.
+
+    Args:
+        circuit (:class:`qibo.models.circuit.Circuit`): parametrized circuit :math:`U(\\theta)`.
+        parameters (ndarray, optional): parameters whose QFIM to calculate.
+            If ``None``, QFIM is calculated with the paremeters from ``circuit``, i.e.
+            ``parameters = circuit.get_parameters()``. Defaults to ``None``.
+        initial_state (ndarray, optional): Initial configuration. It can be specified
+            by the setting the state vector using an array or a circuit. If ``None``,
+            the initial state is :math:`\\ket{0}^{\\otimes n}`. Defaults to ``None``.
+        return_complex (bool, optional): If ``True``, calculates the Jacobian matrix
+            of real and imaginary parts of :math:`\\ket{\\psi(\\theta)}`. If ``False``,
+            calculates only the Jacobian matrix of the real part. Defaults to ``True``.
+        backend (:class:`qibo.backends.abstract.Backend`, optional): backend to be used
+            in the execution. If ``None``, it uses :class:`qibo.backends.GlobalBackend`.
+            Defaults to ``None``.
+
+    Returns:
+        ndarray: Quantum Fisher Information :math:`\\mathbf{F}`.
+    """
+    backend = _check_backend(backend)
+
+    if parameters is None:
+        parameters = circuit.get_parameters()
+        parameters = backend.cast(parameters, dtype=float).flatten()
+
+    jacobian = backend.calculate_jacobian_matrix(
+        circuit, parameters, initial_state, return_complex
+    )
+
+    if return_complex:
+        jacobian = jacobian[0] + 1j * jacobian[1]
+
+    jacobian = backend.cast(jacobian, dtype=np.complex128)
+
+    copied = circuit.copy(deep=True)
+    copied.set_parameters(parameters)
+
+    state = backend.execute_circuit(copied, initial_state=initial_state).state()
+
+    overlaps = jacobian.T @ state
+
+    qfim = jacobian.T @ jacobian
+    qfim = qfim - backend.np.outer(overlaps, backend.np.conj(overlaps.T))
+
+    return 4 * backend.np.real(qfim)
+
+
 def _check_hermitian(matrix, backend=None):
     """Checks if a given matrix is Hermitian.
 

tests/test_quantum_info_metrics.py

@@ -3,6 +3,7 @@
 
 from qibo import Circuit, gates
 from qibo.config import PRECISION_TOL
+from qibo.models.encodings import _generate_rbs_angles, unary_encoder
 from qibo.quantum_info.metrics import (
     average_gate_fidelity,
     bures_angle,
@@ -18,6 +19,7 @@
     process_fidelity,
     process_infidelity,
     purity,
+    quantum_fisher_information_matrix,
     trace_distance,
 )
 from qibo.quantum_info.random_ensembles import (
@@ -386,3 +388,42 @@ def test_frame_potential(backend, nqubits, power_t, samples):
     )
 
     backend.assert_allclose(potential, potential_haar, rtol=1e-2, atol=1e-2)
+
+
+@pytest.mark.parametrize("params_flag", [None, True])
+@pytest.mark.parametrize("return_complex", [False, True])
+@pytest.mark.parametrize("nqubits", [4, 8])
+def test_qfim(backend, nqubits, return_complex, params_flag):
+    if backend.name not in ["pytorch", "tensorflow"]:
+        circuit = Circuit(nqubits)
+        params = np.random.rand(3)
+        params = backend.cast(params, dtype=params.dtype)
+        with pytest.raises(NotImplementedError):
+            test = quantum_fisher_information_matrix(circuit, params, backend=backend)
+    else:
+        # QFIM from https://arxiv.org/abs/2405.20408 is known analytically
+        data = np.random.rand(nqubits)
+        data = backend.cast(data, dtype=data.dtype)
+
+        params = _generate_rbs_angles(data, nqubits, "diagonal")
+        params = backend.cast(params, dtype=np.float64)
+
+        target = [1]
+        for param in params[:-1]:
+            elem = float(target[-1] * backend.np.sin(param) ** 2)
+            target.append(elem)
+        target = 4 * backend.np.diag(backend.cast(target, dtype=np.float64))
+
+        # numerical qfim from quantum_info
+        circuit = unary_encoder(data, "diagonal")
+
+        if params_flag is not None:
+            circuit.set_parameters(params)
+        else:
+            params = params_flag
+
+        qfim = quantum_fisher_information_matrix(
+            circuit, params, return_complex=return_complex, backend=backend
+        )
+
+        backend.assert_allclose(qfim, target, atol=1e-6)