Merge pull request #2416 from owenagnel/master

Improve dnorm efficiency in unitary case.

doc/biblio.rst

@@ -40,6 +40,10 @@ Bibliography
     C. Wood, J. Biamonte, D. G. Cory, *Tensor networks and graphical calculus for
     open quantum systems*. :arxiv:`1111.6950`
 
+.. [AKN98]
+    D. Aharonov, A. Kitaev, and N. Nisan, *Quantum circuits with mixed states*, in Proceedings of the 
+    thirtieth annual ACM symposium on Theory of computing, 20-30 (1998). :arxiv:`quant-ph/9806029`
+
 .. [dAless08]
     D. d’Alessandro, *Introduction to Quantum Control and Dynamics*, (Chapman & Hall/CRC, 2008).
 

doc/changes/2416.feature

@@ -0,0 +1,2 @@
+Updated `qutip.core.metrics.dnorm` to have an efficient speedup when finding the difference of two unitaries. We use a result on page 18 of 
+D. Aharonov, A. Kitaev, and N. Nisan, (1998).

qutip/core/metrics.py

@@ -434,12 +434,21 @@ def hellinger_dist(A, B, sparse=False, tol=0):
 
 def dnorm(A, B=None, solver="CVXOPT", verbose=False, force_solve=False,
           sparse=True):
-    """
+    r"""
     Calculates the diamond norm of the quantum map q_oper, using
     the simplified semidefinite program of [Wat13]_.
 
     The diamond norm SDP is solved by using `CVXPY <https://www.cvxpy.org/>`_.
 
+    If B is provided and both A and B are unitary, then the diamond norm
+    of the difference is calculated more efficiently using the following
+    geometric interpretation:
+    :math:`\|A - B\|_{\diamond}` equals :math:`2 \sqrt(1 - d^2)`, where
+    :math:`d`is the distance between the origin and the convex hull of the
+    eigenvalues of :math:`A B^{\dagger}`.
+    See [AKN98]_ page 18, in the paragraph immediately below the proof of 12.6,
+    as a reference.
+
     Parameters
     ----------
     A : Qobj
@@ -472,59 +481,40 @@ def dnorm(A, B=None, solver="CVXOPT", verbose=False, force_solve=False,
     if cvxpy is None:  # pragma: no cover
         raise ImportError("dnorm() requires CVXPY to be installed.")
 
+    if B is not None and A.dims != B.dims:
+        raise TypeError("A and B do not have the same dimensions.")
+
     # We follow the strategy of using Watrous' simpler semidefinite
     # program in its primal form. This is the same strategy used,
     # for instance, by both pyGSTi and SchattenNorms.jl. (By contrast,
     # QETLAB uses the dual problem.)
 
-    # Check if A and B are both unitaries. If so, then we can without
-    # loss of generality choose B to be the identity by using the
-    # unitary invariance of the diamond norm,
-    #     || A - B ||_♢ = || A B⁺ - I ||_♢.
-    # Then, using the technique mentioned by each of Johnston and
-    # da Silva,
-    #     || A B⁺ - I ||_♢ = max_{i, j} | \lambda_i(A B⁺) - \lambda_j(A B⁺) |,
-    # where \lambda_i(U) is the ith eigenvalue of U.
-
-    # There's a lot of conditions to check for this path.  Only check if they
-    # aren't superoperators.  The difference of unitaries optimization is
-    # currently only implemented for d == 2. Much of the code below is more
-    # general, though, in anticipation of generalizing the optimization.
+    # Check if A and B are both unitaries. If so we can use the geometric
+    # interpretation mentioned in D. Aharonov, A. Kitaev, and N. Nisan. (1998).
+    # We find the eigenvalues of AB⁺ and the distance d between the origin
+    # and the complex hull of these. Plugging this into 2√1-d² gives the
+    # diamond norm.
+
     if (
         not force_solve
+        and A.isunitary
         and B is not None
-        and A.isoper and B.isoper
-        and A.shape[0] == 2
-    ):
-        # Make an identity the same size as A and B to
-        # compare against.
-        I = qeye_like(A)
-        # Compare to B first, so that an error is raised
-        # as soon as possible.
-        Bd = B.dag()
-        if (
-            (B * Bd - I).norm() < 1e-6 and
-            (A * A.dag() - I).norm() < 1e-6
-        ):
-            # Now we are on the fast path, so let's compute the
-            # eigenvalues, then find the diameter of the smallest circle
-            # containing all of them.
-            #
-            # For now, this is only implemented for dim = 2, such that
-            # generalizing here will allow for generalizing the optimization.
-            # A reasonable approach would probably be to use Welzl's algorithm
-            # (https://en.wikipedia.org/wiki/Smallest-circle_problem).
-            U = A * B.dag()
-            eigs = U.eigenenergies()
-            eig_distances = np.abs(eigs[:, None] - eigs[None, :])
-            return np.max(eig_distances)
-
-    # Force the input superoperator to be a Choi matrix.
+        and B.isunitary
+    ):  # Special optimisation for a difference of unitaries.
+        U = A * B.dag()
+        eigs = U.eigenenergies()
+        d = _find_poly_distance(eigs)
+        return 2 * np.sqrt(1 - d**2)  # plug d into formula
+
     J = to_choi(A)
 
-    if B is not None:
+    if B is not None:  # If B is provided, calculate difference
         J -= to_choi(B)
 
+    if not force_solve and J.iscptp:
+        # diamond norm of a CPTP map is 1 (Prop 3.44 Watrous 2018)
+        return 1.0
+
     # Watrous 2012 also points out that the diamond norm of Lambda
     # is the same as the completely-bounded operator-norm (∞-norm)
     # of the dual map of Lambda. We can evaluate that norm much more
@@ -583,3 +573,34 @@ def unitarity(oper):
     """
     Eu = _to_superpauli(oper).full()[1:, 1:]
     return np.linalg.norm(Eu, 'fro')**2 / len(Eu)
+
+
+def _find_poly_distance(eigenvals: np.ndarray) -> float:
+    """
+    Returns the distance between the origin and the convex hull of eigenvalues.
+
+    The complex eigenvalues must have unit length (i.e. lie on the circle
+    about the origin).
+    """
+    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)

qutip/tests/core/test_metrics.py

@@ -439,22 +439,25 @@ def test_qubit_triangle(self, dimension):
         assert dnorm(A + B) <= dnorm(A) + dnorm(B) + 1e-7
 
     @pytest.mark.repeat(3)
-    @pytest.mark.parametrize("generator", [
-        pytest.param(rand_super_bcsz, id="super"),
-        pytest.param(rand_unitary, id="unitary"),
-    ])
-    def test_force_solve(self, dimension, generator):
-        """
-        Metrics: checks that special cases for dnorm agree with SDP solutions.
-        """
-        A, B = generator(dimension), generator(dimension)
+    def test_unitary_case(self, dimension):
+        """Check that the diamond norm is one for unitary maps."""
+        A, B = rand_unitary(dimension), rand_unitary(dimension)
         assert (
-            dnorm(A, B, force_solve=False)
+            dnorm(A, B)
             == pytest.approx(dnorm(A, B, force_solve=True), abs=1e-5)
         )
 
     @pytest.mark.repeat(3)
-    def test_cptp(self, dimension, sparse):
+    def test_cp_case(self, dimension):
+        """Check that the diamond norm is one for unitary maps."""
+        A = rand_super_bcsz(dimension, enforce_tp=False)
+        assert (
+            dnorm(A)
+            == pytest.approx(dnorm(A, force_solve=True), abs=1e-5)
+        )
+
+    @pytest.mark.repeat(3)
+    def test_cptp_case(self, dimension, sparse):
         """Check that the diamond norm is one for CPTP maps."""
         A = rand_super_bcsz(dimension)
         assert A.iscptp