Hyper-efficient unitary diamond distance

qiskit/quantum_info/__init__.py

@@ -64,6 +64,7 @@
 .. autofunction:: process_fidelity
 .. autofunction:: gate_error
 .. autofunction:: diamond_norm
+.. autofunction:: diamond_distance
 .. autofunction:: state_fidelity
 .. autofunction:: purity
 .. autofunction:: concurrence
@@ -128,7 +129,13 @@
 )
 from .operators.channel import PTM, Chi, Choi, Kraus, Stinespring, SuperOp
 from .operators.dihedral import CNOTDihedral
-from .operators.measures import average_gate_fidelity, diamond_norm, gate_error, process_fidelity
+from .operators.measures import (
+    average_gate_fidelity,
+    diamond_norm,
+    diamond_distance,
+    gate_error,
+    process_fidelity,
+)
 from .random import (
     random_clifford,
     random_cnotdihedral,

qiskit/quantum_info/operators/__init__.py

@@ -15,7 +15,13 @@
 from __future__ import annotations
 from .channel import PTM, Chi, Choi, Kraus, Stinespring, SuperOp
 from .dihedral import CNOTDihedral
-from .measures import average_gate_fidelity, diamond_norm, gate_error, process_fidelity
+from .measures import (
+    average_gate_fidelity,
+    diamond_norm,
+    diamond_distance,
+    gate_error,
+    process_fidelity,
+)
 from .operator import Operator
 from .scalar_op import ScalarOp
 from .symplectic import (

qiskit/quantum_info/operators/measures.py

@@ -22,6 +22,7 @@
 from qiskit.circuit.gate import Gate
 from qiskit.quantum_info.operators.base_operator import BaseOperator
 from qiskit.quantum_info.operators.operator import Operator
+from qiskit.quantum_info.operators.symplectic.pauli import Pauli
 from qiskit.quantum_info.operators.channel.quantum_channel import QuantumChannel
 from qiskit.quantum_info.operators.channel import Choi, SuperOp
 from qiskit.quantum_info.states.densitymatrix import DensityMatrix
@@ -350,6 +351,110 @@ def cvx_bmat(mat_r, mat_i):
     return sol
 
 
+def diamond_distance(op1: BaseOperator, op2: BaseOperator) -> float:
+    r"""Return the diamond distance between two unitary operators.
+
+    This function computes the completely-bounded trace-norm (often
+    referred to as the diamond norm) of a difference of two quantum maps.
+    It is often used as a distance metric in quantum information where
+
+    :math:`d(E, F) = ||E-F||_diamond`.
+
+
+    Args:
+        op1 (Operator): a unitary operator.
+        op2 (Operator): a unitary operator.
+
+    Returns:
+        float: The completely-bounded trace norm of `op1 - op2`.
+
+    Raises:
+        ValueError: if the input operators do not have the same dimensions.
+
+    Additional Information:
+        If both operators are unitary, the implementation uses a result in Aharonov
+        et al. (1998) resulting in a significant optimisation. Geometrically, we
+        compute the distance :math:`d` between the origin and the convex hull of
+        the eigenvalues of :math:`U^\dag V` which is plugged into :math:`\sqrt{1 - d^2}`.
+
+        If the operators are `Pauli` objects, we use an optimisation to find eigenvalues.
+
+    Reference:
+        D. Aharonov, A. Kitaev, and N. Nisan. “Quantum circuits with
+        mixed states” in Proceedings of the thirtieth annual ACM symposium
+        on Theory of computing, pp. 20-30, 1998.
+
+    .. note::
+
+        This function requires the optional CVXPY package to be installed.
+        Any additional kwargs will be passed to the ``cvxpy.solve``
+        function. See the CVXPY documentation for information on available
+        SDP solvers.
+    """
+    op1 = _input_formatter(op1, BaseOperator, "diamond_distance", "op1")
+    op2 = _input_formatter(op2, BaseOperator, "diamond_distance", "op2")
+
+    # Check base operators have same dimension
+    if op1.dim != op2.dim:
+        raise ValueError(
+            "Input quantum channel and target unitary must have the same "
+            f"dimensions ({op1.dim} != {op2.dim})."
+        )
+
+    # Attempt to run unitary optimisation
+    if (
+        isinstance(op1, Operator)
+        and isinstance(op2, Operator)
+        and op1.is_unitary()
+        and op2.is_unitary()
+    ):
+        # Compute the diamond norm
+        mat1 = op1.data
+        mat2 = op2.data
+        pre_diag = np.conj(mat1).T @ mat2
+        eigenvals = np.linalg.eigvals(pre_diag)
+        d = _find_poly_distance(eigenvals)
+        return 2 * np.sqrt(1 - d**2)
+    elif isinstance(op1, Pauli) and isinstance(op2, Pauli):
+        # Check Pauli equality up to phase
+        if (op1.z == op2.z).all() and (op1.x == op2.x).all():
+            return 0.0
+        else:
+            return 2.0
+    else:
+        # TODO: Implement special case for pauli channels (Benenti and Strini 2010)
+        # as well as a potential optimisation for clifford circuits
+
+        # Compute the diamond norm
+        return diamond_norm(Choi(op1) - Choi(op2))
+
+
+def _find_poly_distance(eigenvals: np.ndarray) -> float:
+    """Function to find the distance between origin and the convex hull of eigenvalues."""
+    phases = np.angle(eigenvals)
+    phase_max = phases.max()
+    phase_min = phases.min()
+
+    if phase_min > 0:  # all eigenvals have pos phase: hull is above x axis
+        return np.cos((phase_max - phase_min) / 2)
+
+    if phase_max <= 0:  # all eigenvals have neg phase: hull is below x axis
+        return np.cos((np.abs(phase_min) - np.abs(phase_max)) / 2)
+
+    pos_phase_min = np.where(phases > 0, phases, np.inf).min()
+    neg_phase_max = np.where(phases <= 0, phases, -np.inf).max()
+
+    big_angle = phase_max - phase_min
+    small_angle = pos_phase_min - neg_phase_max
+    if big_angle >= np.pi:
+        if small_angle <= np.pi:  # hull contains the origin
+            return 0
+        else:  # hull is left of y axis
+            return np.cos((2 * np.pi - small_angle) / 2)
+    else:  # hull is right of y axis
+        return np.cos(big_angle / 2)
+
+
 def _cvxpy_check(name):
     """Check that a supported CVXPY version is installed"""
     # Check if CVXPY package is installed
@@ -382,6 +487,10 @@ def _input_formatter(obj, fallback_class, func_name, arg_name):
     if hasattr(obj, "to_channel"):
         return obj.to_channel()
 
+    # Pauli-like input
+    if isinstance(obj, Pauli):
+        return obj
+
     # Unitary-like input
     if isinstance(obj, (Gate, BaseOperator)):
         return Operator(obj)

releasenotes/notes/unitary-diamond-feature-fd953e0d0bbbb073.yaml

@@ -0,0 +1,9 @@
+---
+features_quantum_info:
+  - |
+    Added :meth:`qiskit.quantum_info.diamond_distance` function for computing the
+    distance between two quantum channels. This is equivalent to using
+    :meth:`qiskit.quantum_info.diamond_norm` to calculate the distance of two 
+    channels: ``diamond_norm(chan1 - chan2)``. If both channels are unitary, the 
+    implementation runs an optimisation. Refer to
+    `#12341 <https://github.com/Qiskit/qiskit/issues/12341>` for more details.

test/python/quantum_info/operators/test_measures.py

@@ -12,6 +12,7 @@
 """Tests for operator measures."""
 
 import unittest
+import importlib.util
 from ddt import ddt
 
 import numpy as np
@@ -21,6 +22,9 @@
 from qiskit.quantum_info import average_gate_fidelity
 from qiskit.quantum_info import gate_error
 from qiskit.quantum_info import diamond_norm
+from qiskit.quantum_info import diamond_distance
+from qiskit.quantum_info.random import random_unitary, random_pauli, random_clifford
+from qiskit.circuit.library import RZGate
 from test import combine  # pylint: disable=wrong-import-order
 from test import QiskitTestCase  # pylint: disable=wrong-import-order
 
@@ -163,9 +167,7 @@ def test_channel_gate_error(self):
     @combine(num_qubits=[1, 2, 3])
     def test_diamond_norm(self, num_qubits):
         """Test the diamond_norm for {num_qubits}-qubit pauli channel."""
-        try:
-            import cvxpy
-        except ImportError:
+        if importlib.util.find_spec("cvxpy") is None:
             # Skip test if CVXPY not installed
             self.skipTest("CVXPY not installed.")
 
@@ -178,11 +180,65 @@ def test_diamond_norm(self, num_qubits):
             op = op + coeff * Choi(Operator.from_label(label))
         target = np.sum(np.abs(coeffs))
 
-        try:
-            value = diamond_norm(op)
-            self.assertAlmostEqual(value, target, places=4)
-        except cvxpy.SolverError:
-            self.skipTest("CVXPY solver failed.")
+        value = diamond_norm(op)
+        self.assertAlmostEqual(value, target, places=4)
+
+    def test_diamond_distance(self):
+        """Test the diamond_distance function for RZGates
+        with a specific set of angles."""
+        if importlib.util.find_spec("cvxpy") is None:
+            # Skip test if CVXPY not installed
+            self.skipTest("CVXPY not installed.")
+        angles = np.linspace(0, 2 * np.pi, 10, endpoint=False)
+        for angle in angles:
+            op1 = Operator(RZGate(angle))
+            op2 = Operator.from_label("I")
+            d2 = np.cos(angle / 2) ** 2  # analytical formula for hull distance
+            target = np.sqrt(1 - d2) * 2
+            self.assertAlmostEqual(diamond_distance(op1, op2), target, places=7)
+
+    @combine(num_qubits=[1, 2])
+    def test_diamond_distance_random(self, num_qubits):
+        """Tests the diamond_distance for random unitaries.
+        Compares results with semi-definite program."""
+        if importlib.util.find_spec("cvxpy") is None:
+            # Skip test if CVXPY not installed
+            self.skipTest("CVXPY not installed.")
+
+        rng = np.random.default_rng(1234)
+        for _ in range(5):
+            op1 = random_unitary(2**num_qubits, seed=rng)
+            op2 = random_unitary(2**num_qubits, seed=rng)
+            target = diamond_norm(Choi(op1) - Choi(op2))
+            self.assertAlmostEqual(diamond_distance(op1, op2), target, places=4)
+
+    @combine(num_qubits=[1, 2])
+    def test_diamond_distance_random_pauli(self, num_qubits):
+        """Test diamond_distance for non-CP channel"""
+        if importlib.util.find_spec("cvxpy") is None:
+            # Skip test if CVXPY not installed
+            self.skipTest("CVXPY not installed.")
+
+        rng = np.random.default_rng(1234)
+        for _ in range(5):
+            op1 = random_pauli(2**num_qubits, seed=rng)
+            op2 = random_pauli(2**num_qubits, seed=rng)
+            target = diamond_norm(Choi(op1) - Choi(op2))
+            self.assertAlmostEqual(diamond_distance(op1, op2), target, places=4)
+
+    @combine(num_qubits=[1, 2])
+    def test_diamond_distance_random_clifford(self, num_qubits):
+        """Test diamond_distance for non-CP channel"""
+        if importlib.util.find_spec("cvxpy") is None:
+            # Skip test if CVXPY not installed
+            self.skipTest("CVXPY not installed.")
+
+        rng = np.random.default_rng(1234)
+        for _ in range(5):
+            op1 = random_clifford(2**num_qubits, seed=rng)
+            op2 = random_clifford(2**num_qubits, seed=rng)
+            target = diamond_norm(Choi(op1) - Choi(op2))
+            self.assertAlmostEqual(diamond_distance(op1, op2), target, places=4)
 
 
 if __name__ == "__main__":