Add magnitude approximation protocol (#1750)

- adds the magnitude approximation protocol to approximate general
unitaries
- fixes the docstring for the Z-rotation protocols to be $e^{i\theta Z}$
( $= Rz(-2\theta)$ ) instead of $Rz(2\theta)$

Author: NoureldinYosri

qualtran/rotation_synthesis/__init__.py

@@ -13,9 +13,10 @@
 #  limitations under the License.
 
 from qualtran.rotation_synthesis.math_config import NumpyConfig, with_dps
-from qualtran.rotation_synthesis.protocols.clifford_t_synthesis import (
+from qualtran.rotation_synthesis.protocols import (
     diagonal_unitary_approx,
     fallback_protocol,
+    magnitude_approx,
     mixed_diagonal_protocol,
     mixed_fallback_protocol,
 )

qualtran/rotation_synthesis/channels/channel.py

@@ -15,9 +15,10 @@
 from __future__ import annotations
 
 import abc
-from typing import Optional, Sequence
+from typing import Optional, Sequence, Union
 
 import attrs
+import numpy as np
 
 import qualtran.rotation_synthesis._typing as rst
 import qualtran.rotation_synthesis.math_config as mc
@@ -33,7 +34,7 @@ def expected_num_ts(self, config: mc.MathConfig) -> rst.Real: ...
 
     @abc.abstractmethod
     def diamond_norm_distance_to_rz(self, theta: rst.Real, config: mc.MathConfig) -> rst.Real:
-        r"""Returns the diamond norm distance to $Rz(2\theta)$."""
+        r"""Returns the diamond norm distance to $e^{i\theta Z}$."""
 
 
 @attrs.frozen
@@ -110,6 +111,44 @@ def from_sequence(cls, seq: Sequence[str], twirl: bool = False) -> UnitaryChanne
         n = sum(g.startswith("T") for g in seq)
         return UnitaryChannel(u.matrix[0, 0], u.matrix[1, 0], n, twirl)
 
+    def diamond_norm_distance_to_unitary(
+        self, unitary: np.ndarray, config: mc.MathConfig
+    ) -> rst.Real:
+        r"""Returns the diamond norm distance between self and the given unitary.
+
+        From Theorem B.1 of arxiv:2203.10064, the diamond norm distance between two untiaries
+        $U, V$ is $|v_1 - v_0|$ where $v_i$ are the eigen values of $V^\dagger U$. Geometrically
+        this is the diameter of the smallest disc in the complex plane that contains both
+        eigenvalues.
+        """
+        # W = V^\dagger U
+        w = self.to_matrix().adjoint().numpy(config) @ unitary
+        # Compute the eigen values of W
+        a = config.one
+        b = -w[0, 0] - w[1, 1]
+        c = w[0, 0] * w[1, 1] - w[0, 1] * w[1, 0]
+        d = config.sqrt(b**2 - 4 * a * c)
+        eigv0, eigv1 = [(-b - d) / (2 * a), (-b + d) / (2 * a)]
+        # Compute the norm of the difference.
+        diameter_vec = eigv1 - eigv0
+        return config.sqrt(diameter_vec.real**2 + diameter_vec.imag**2)
+
+    @classmethod
+    def from_unitaries(
+        cls, *unitaries: Union[UnitaryChannel, su2_ct.SU2CliffordT]
+    ) -> UnitaryChannel:
+        if not unitaries:
+            raise ValueError('at least one unitary should be provided')
+
+        unitary = su2_ct.ISqrt2
+        for u in unitaries:
+            if isinstance(u, UnitaryChannel):
+                unitary = unitary @ u.to_matrix()
+            else:
+                unitary = unitary @ u
+        unitary = unitary.rescale()
+        return UnitaryChannel(unitary.matrix[0, 0], unitary.matrix[1, 0], unitary.num_t_gates())
+
 
 @attrs.frozen
 class ProjectiveChannel(Channel):

qualtran/rotation_synthesis/channels/channel_test.py

@@ -46,6 +46,10 @@ def test_unitary_from_sequence(gates):
 def test_diamond_distance_for_unitary(gates, theta, distance):
     c = ch.UnitaryChannel.from_sequence(gates)
     np.testing.assert_allclose(c.diamond_norm_distance_to_rz(theta, mc.NumpyConfig), distance)
+    u = np.zeros((2, 2), complex)
+    u[1, 1] = np.exp(-1j * theta)
+    u[0, 0] = u[1, 1].conjugate()
+    np.testing.assert_allclose(c.diamond_norm_distance_to_unitary(u, mc.NumpyConfig), distance)
 
 
 @pytest.mark.parametrize(

qualtran/rotation_synthesis/lattice/geometry.py

@@ -213,7 +213,7 @@ def plot(self, ax: Optional[plt.Axes] = None, add_label: bool = True, **patch_ar
             fig, ax = plt.subplots(1)
 
         theta = float(self.tilt(mc.NumpyConfig))
-        e = self.rotate(-theta, mc.NumpyConfig)
+        e = self.rotate(theta, mc.NumpyConfig)
         w = 2 / np.sqrt(float(e.D[0, 0]))
         h = 2 / np.sqrt(float(e.D[1, 1]))
         c = self.center.astype(float).tolist()

qualtran/rotation_synthesis/math_config.py

@@ -40,6 +40,7 @@ class MathConfig:
     _ceil: Callable[[Real], int] = math.ceil
     _arctan2: Callable[[Real, Real], Real] = math.atan2
     _arcsin: Callable[[Real], Real] = math.asin
+    _arccos: Callable[[Real], Real] = math.acos
     _number: Callable[[Real], Real] = float
 
     def number(self, x) -> Real:
@@ -75,6 +76,9 @@ def __hash__(self) -> int:
     def arcsin(self, x: Real) -> Real:
         return self._arcsin(x)
 
+    def arccos(self, x: Real) -> Real:
+        return self._arccos(x)
+
 
 NumpyConfig = MathConfig(
     "numpy",
@@ -91,6 +95,7 @@ def arcsin(self, x: Real) -> Real:
     lambda x: int(np.ceil(x)),
     np.arctan2,
     np.arcsin,
+    np.arccos,
     np.float128,
 )
 
@@ -124,5 +129,6 @@ def with_dps(dps: int) -> MathConfig:
         lambda x: int(mpmath.ceil(x)),
         mpmath.atan2,
         mpmath.asin,
+        mpmath.acos,
         mpmath.mpf,
     )

qualtran/rotation_synthesis/matrix/analytical_decomposition.py

@@ -0,0 +1,56 @@
+#  Copyright 2025 Google LLC
+#
+#  Licensed under the Apache License, Version 2.0 (the "License");
+#  you may not use this file except in compliance with the License.
+#  You may obtain a copy of the License at
+#
+#      https://www.apache.org/licenses/LICENSE-2.0
+#
+#  Unless required by applicable law or agreed to in writing, software
+#  distributed under the License is distributed on an "AS IS" BASIS,
+#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#  See the License for the specific language governing permissions and
+#  limitations under the License.
+
+import numpy as np
+
+import qualtran.rotation_synthesis._typing as rst
+import qualtran.rotation_synthesis.math_config as mc
+
+
+def su_unitary_to_zxz_angles(
+    u: np.ndarray, config: mc.MathConfig
+) -> tuple[rst.Real, rst.Real, rst.Real]:
+    det = u[0, 0] * u[1, 1] - u[0, 1] * u[1, 0]
+    # print(det)
+    # assert config.isclose(det, 1)
+
+    cos_phi = (u[1, 1] * u[0, 0] + u[1, 0] * u[0, 1]).real
+    # clip to the range [-1, 1] to ensure numerical stability.
+    cos_phi = min(1, max(-1, cos_phi))
+    phi = config.arccos(cos_phi)
+
+    sum_half_theta = config.arctan2(u[1, 1].imag, u[1, 1].real)
+    diff_half_theta = config.arctan2(u[1, 0].imag, u[1, 0].real) - 1.5 * config.pi
+
+    theta1 = sum_half_theta + diff_half_theta
+    theta2 = sum_half_theta - diff_half_theta
+
+    return theta1, phi, theta2
+
+
+def rx(phi: rst.Real, config: mc.MathConfig) -> np.ndarray:
+    c = config.cos(phi / 2)
+    s = config.sin(phi / 2)
+    return np.array([[c, -1j * s], [-1j * s, c]])
+
+
+def rz(theta: rst.Real, config: mc.MathConfig) -> np.ndarray:
+    v = config.cos(theta / 2) + 1j * config.sin(theta / 2)
+    return np.diag([v.conjugate(), v])
+
+
+def unitary_from_zxz(
+    theta1: rst.Real, phi: rst.Real, theta2: rst.Real, config: mc.MathConfig
+) -> np.ndarray:
+    return rz(theta1, config) @ rx(phi, config) @ rz(theta2, config)

qualtran/rotation_synthesis/matrix/analytical_decomposition_test.py

@@ -0,0 +1,46 @@
+#  Copyright 2025 Google LLC
+#
+#  Licensed under the Apache License, Version 2.0 (the "License");
+#  you may not use this file except in compliance with the License.
+#  You may obtain a copy of the License at
+#
+#      https://www.apache.org/licenses/LICENSE-2.0
+#
+#  Unless required by applicable law or agreed to in writing, software
+#  distributed under the License is distributed on an "AS IS" BASIS,
+#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#  See the License for the specific language governing permissions and
+#  limitations under the License.
+
+import numpy as np
+import pytest
+from scipy import stats
+
+import qualtran.rotation_synthesis.math_config as mc
+from qualtran.rotation_synthesis.matrix import analytical_decomposition as ad
+
+
+def _random_angles(n, seed):
+    rng = np.random.default_rng(seed)
+    for _ in range(n):
+        yield rng.random(3) * 2 * np.pi
+
+
+def _random_su2(n, seed):
+    for u in stats.unitary_group(2).rvs(n, seed):
+        yield u / np.linalg.det(u) ** 0.5
+
+
+@pytest.mark.parametrize(["theta1", "phi", "theta2"], _random_angles(100, seed=0))
+def test_decomposition_round_trip(theta1, phi, theta2):
+    u = ad.unitary_from_zxz(theta1, phi, theta2, mc.NumpyConfig)
+    angles = ad.su_unitary_to_zxz_angles(u, config=mc.NumpyConfig)
+    v = ad.unitary_from_zxz(*angles, config=mc.NumpyConfig)
+    np.testing.assert_allclose(actual=u, desired=v)
+
+
+@pytest.mark.parametrize("u", _random_su2(100, seed=0))
+def test_decomposition_round_trip_start_with_unitary(u):
+    angles = ad.su_unitary_to_zxz_angles(u, config=mc.NumpyConfig)
+    v = ad.unitary_from_zxz(*angles, config=mc.NumpyConfig)
+    np.testing.assert_allclose(actual=u, desired=v)

qualtran/rotation_synthesis/matrix/generation_test.py

@@ -28,7 +28,8 @@ def test_generated_rotations_determinant():
     all_rotations = generate_rotations(5)
     for n in range(len(all_rotations)):
         assert np.allclose(
-            [abs(np.linalg.det(r.numpy())) for r in all_rotations[n]], 2 * (2 + np.sqrt(2)) ** n
+            [abs(np.linalg.det(r.matrix.astype(complex))) for r in all_rotations[n]],
+            2 * (2 + np.sqrt(2)) ** n,
         )
 
 
@@ -58,4 +59,4 @@ def test_generated_rotations_unitary():
             assert gates is not None
             for p in gates:
                 u = _numpy_matrix_for_symbol[p] @ u
-            assert _are_close_up_to_global_phase(u, r.numpy())
+            assert _are_close_up_to_global_phase(u, r.matrix.astype(complex))

qualtran/rotation_synthesis/matrix/su2_ct.py

@@ -19,6 +19,7 @@
 import attrs
 import numpy as np
 
+import qualtran.rotation_synthesis.math_config as mc
 from qualtran.rotation_synthesis.rings import zsqrt2, zw
 
 
@@ -81,8 +82,23 @@ def __hash__(self):
     def __eq__(self, other):
         return np.all(self.matrix == other.matrix)
 
-    def numpy(self) -> np.ndarray:
-        return self.matrix.astype(complex)
+    def numpy(self, config: Optional[mc.MathConfig] = None) -> np.ndarray:
+        """Returns the numpy representation of the unitary.
+        Args:
+            config: An optional MathConfig used to convert the matrix entries to complex
+                numbers and normalize the result. If not given numpy methods are used.
+        """
+        if config is None:
+            result = self.matrix.astype(complex)
+            result = result / np.linalg.det(result) ** 0.5
+            return result
+        result = np.zeros((2, 2)) + 1j * config.zero
+        sqrt_det = config.sqrt(self.det().value(config.sqrt2))
+        for i in range(2):
+            for j in range(2):
+                result[i, j] = self.matrix[i, j].value(config.sqrt2)
+        result = result / sqrt_det
+        return result
 
     def adjoint(self) -> "SU2CliffordT":
         return SU2CliffordT(self.matrix.T.conj())
@@ -206,16 +222,32 @@ def is_valid(self) -> bool:
         _, _, n1 = self.matrix[1, 0].to_zsqrt2()
         return n0 == n1
 
+    def rescale(self) -> 'SU2CliffordT':
+        r"""Rescales the matrix such that its determinant is minimized.
 
-_KEY_MAP = {
-    (0, 0, 0, 0): "Tx",
-    (0, 0, 1, 0): "Tz",
-    (0, 1, 0, 0): "Ty",
-    (0, 1, 1, 0): "Tx",
-    (1, 0, 1, 0): "Tx",
-    (1, 1, 0, 0): "Tx",
-    (1, 1, 1, 0): "Tz",
-}
+        The determinant of the unitary can be written as $2\lambda^n$ where $\lambda=2+\sqrt{l}$
+        and $n$ is the number of $T$ gates needed to synthesize the matrix. When all entries of
+        the matrix are divisible by $\lambda$ then we can divide through by $\lambda$ to reduce $n$
+        """
+        u = self
+        while u.det() > 2 * zsqrt2.LAMBDA_KLIUCHNIKOV:
+            if not all(a.is_divisible_by(zw.LAMBDA_KLIUCHNIKOV) for a in u.matrix.flat):
+                break
+            u = SU2CliffordT(
+                [[x // zw.LAMBDA_KLIUCHNIKOV for x in row] for row in u.matrix], u.gates
+            )
+        return u
+
+    def num_t_gates(self) -> int:
+        """Returns the number of T gates needed to synthesize the matrix."""
+        det = self.det()
+        x = zsqrt2.ZSqrt2(2)
+        n = 0
+        while x < det:
+            x = x * zsqrt2.LAMBDA_KLIUCHNIKOV
+            n += 1
+        assert x == det
+        return n
 
 
 @functools.cache

qualtran/rotation_synthesis/matrix/su2_ct_test.py

@@ -18,6 +18,7 @@
 import numpy as np
 import pytest
 
+import qualtran.rotation_synthesis.math_config as mc
 from qualtran.rotation_synthesis.matrix import su2_ct
 from qualtran.rotation_synthesis.rings import zsqrt2, zw
 
@@ -54,9 +55,9 @@ def test_multiply(seq):
     for i in seq:
         g = g @ _GATES[i]
         k += i >= 3  # is a T gate.
-        g_numpy = (g_numpy @ _GATES[i].numpy()) / _SQRT2
+        g_numpy = (g_numpy @ _GATES[i].matrix.astype(complex)) / _SQRT2
     assert g.det() == 2 * _LAMBDA**k
-    np.testing.assert_allclose(g.numpy() / _SQRT2, g_numpy, atol=1e-9)
+    np.testing.assert_allclose(g.matrix.astype(complex) / _SQRT2, g_numpy, atol=1e-9)
 
 
 @pytest.mark.parametrize("g", _make_random_su(10, 3, random_cliffords=False, seed=0))
@@ -76,9 +77,9 @@ def test_adjoint(g):
     ["g", "g_numpy"], [[su2_ct.Tx, _TX_numpy], [su2_ct.Ty, _TY_numpy], [su2_ct.Tz, _TZ_numpy]]
 )
 def test_t_gates(g, g_numpy):
-    np.testing.assert_allclose(g.numpy() / _SQRT2 / np.sqrt(2 + _SQRT2), g_numpy)
+    np.testing.assert_allclose(g.matrix.astype(complex) / _SQRT2 / np.sqrt(2 + _SQRT2), g_numpy)
     np.testing.assert_allclose(
-        g.adjoint().numpy() / _SQRT2 / np.sqrt(2 + _SQRT2), g_numpy.T.conjugate()
+        g.adjoint().matrix.astype(complex) / _SQRT2 / np.sqrt(2 + _SQRT2), g_numpy.T.conjugate()
     )
 
 
@@ -101,8 +102,33 @@ def test_generate_cliffords():
     cirq_cliffords = [
         cirq.unitary(c) for c in cirq.SingleQubitCliffordGate.all_single_qubit_cliffords
     ]
-    assert np.allclose(np.abs([np.linalg.det(c.numpy()) for c in cliffords]), 2)
+    assert np.allclose(np.abs([np.linalg.det(c.matrix.astype(complex)) for c in cliffords]), 2)
     sqrt2 = np.sqrt(2)
     for c in cliffords:
-        u = c.numpy() / sqrt2
+        u = c.matrix.astype(complex) / sqrt2
         assert np.any([are_close_up_to_global_phase(u, c) for c in cirq_cliffords])
+
+
+@pytest.mark.parametrize("g", _make_random_su(50, 5, random_cliffords=True, seed=0))
+@pytest.mark.parametrize("config", [None, mc.NumpyConfig])
+def test_rescale(g: su2_ct.SU2CliffordT, config):
+    np.testing.assert_allclose(g.numpy(config), g.rescale().numpy(config))
+
+
+def test_num_t_gates():
+    for clifford in su2_ct.generate_cliffords():
+        assert clifford.num_t_gates() == 0
+
+        for t in su2_ct.Ts:
+            assert (clifford @ t).num_t_gates() == 1
+
+    for t1 in su2_ct.Ts:
+        for t2 in su2_ct.Ts:
+            assert (t1 @ t2).num_t_gates() == 2
+
+    for t in su2_ct.Ts:
+        # still two T gates
+        assert (t @ t.adjoint()).num_t_gates() == 2
+
+        # We need to call .rescale to remove the common factor and reduce the T count.
+        assert (t @ t.adjoint()).rescale().num_t_gates() == 0