Hyper-efficient unitary diamond distance

qiskit/quantum_info/__init__.py

@@ -64,6 +64,7 @@
 .. autofunction:: process_fidelity
 .. autofunction:: gate_error
 .. autofunction:: diamond_norm
+.. autofunction:: unitary_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,
+    unitary_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,
+    unitary_diamond_distance,
+    gate_error,
+    process_fidelity,
+)
 from .operator import Operator
 from .scalar_op import ScalarOp
 from .symplectic import (

qiskit/quantum_info/operators/measures.py

@@ -350,6 +350,94 @@ def cvx_bmat(mat_r, mat_i):
     return sol
 
 
+def unitary_diamond_distance(channel1: Operator, channel2: Operator) -> 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 the difference of two unitary channels
+    (or unitary operators). This is often used as a distance metric in quantum
+    information This method is a more specialised implementation of
+    :meth:`~qiskit.quantum_info.diamond_norm`.
+
+    Args:
+        channel1 (Operator): a unitary operator.
+        channel2 (Operator): a unitary operator.
+
+    Returns:
+        float: The completely-bounded trace norm :math:`\|U - V\|_{\diamond}`.
+
+    Raises:
+        ValueError: if the input operators do not have the same dimensions or aren't unitary.
+
+    Additional Information:
+        The implementation uses an unproved result in Aharonov et al. (1998). Geometrically,
+        we compute the distance :math:`d` between the origin and the convex hull
+        of the eigenvalues which is then plugged into :math:`\sqrt{1 - d^2}`.
+
+    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.
+    """
+    op1 = _input_formatter(channel1, Operator, "unitary_diamond_distance", "channel1")
+    op2 = _input_formatter(channel2, Operator, "unitary_diamond_distance", "channel2")
+
+    # Check operators are unitary and have same dimension
+    if not op1.is_unitary():
+        raise ValueError(
+            "Invalid operator supplied to channel1 of unitary_diamond_distance"
+            "operators must be unitary."
+        )
+    if not op2.is_unitary():
+        raise ValueError(
+            "Invalid operator supplied to channel2 of unitary_diamond_distance"
+            "operators must be unitary."
+        )
+    if op1.dim != op2.dim:
+        raise ValueError(
+            "Input quantum channel and target unitary must have the same "
+            f"dimensions ({op1.dim} != {op2.dim})."
+        )
+
+    # Compute the diamond norm
+    op1 = op1.data
+    op2 = op2.data
+    pre_diag = np.conj(op1).T @ op2
+    eigenvals = np.linalg.eigvals(pre_diag)
+    d = _find_poly_distance(eigenvals)
+    return 2 * np.sqrt(1 - d**2)
+
+
+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)
+    pos_max = phases.max()
+    neg_min = phases.min()
+    pos_min = np.where(phases > 0, phases, np.inf).min()
+    neg_max = np.where(phases <= 0, phases, -np.inf).max()
+
+    if neg_min > 0:
+        # all eigenvalues have positive phase, so hull is above x axis
+        return np.cos((pos_max - pos_min) / 2)
+
+    if pos_max <= 0:
+        # all eigenvalues have negative phase, so hull is below x axis
+        return np.cos((np.abs(neg_min) - np.abs(neg_max)) / 2)
+
+    big_angle = pos_max - neg_min
+    small_angle = pos_min - neg_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

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

@@ -0,0 +1,9 @@
+---
+features_quantum_info:
+  - |
+    Added :meth:`qiskit.quantum_info.unitary_diamond_distance` function for computing the
+    distance between two unitary channels. This is equivalent to using
+    :meth:`qiskit.quantum_info.unitary_diamond_distance` to calculate the distance of two 
+    channels: ``diamond_norm(chan1 - chan2)``. In the case of unitary channels this 
+    implementation is significantly more efficient. 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 unitary_diamond_distance
+from qiskit.quantum_info.random import random_unitary
+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,33 @@ 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_unitary_diamond_distance(self):
+        """Test the unitary_diamond_distance function for RZGates
+        with a specific set of angles."""
+        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(unitary_diamond_distance(op1, op2), target, places=7)
+
+    @combine(num_qubits=[1, 2, 3])
+    def test_unitary_diamond_distance_random(self, num_qubits):
+        """Tests the unitary_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)
+        op1 = random_unitary(2**num_qubits, seed=rng)
+        op2 = random_unitary(2**num_qubits, seed=765720)
+        target = diamond_norm(Choi(op1) - Choi(op2))
+        self.assertAlmostEqual(unitary_diamond_distance(op1, op2), target, places=4)
 
 
 if __name__ == "__main__":