Merge branch 'main' into rot_vis

Author: NoureldinYosri

qualtran/rotation_synthesis/__init__.py

@@ -14,9 +14,10 @@
 
 from qualtran.rotation_synthesis.math_config import NumpyConfig, with_dps
 from qualtran.rotation_synthesis.matrix import to_cirq, to_quirk, to_sequence
-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,10 +15,11 @@
 from __future__ import annotations
 
 import abc
-from typing import Optional, Sequence
+from typing import Optional, Sequence, Union
 
 import attrs
 import cirq
+import numpy as np
 
 import qualtran.rotation_synthesis._typing as rst
 import qualtran.rotation_synthesis.math_config as mc
@@ -35,7 +36,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
@@ -147,6 +148,44 @@ def to_quirk(self, fmt: str = "xz") -> str:
         cols = '[' + ','.join(f'[{g}]' for g in gates) + ']'
         return "https://algassert.com/quirk#circuit={\"cols\":%s}" % cols
 
+    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

@@ -48,6 +48,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

@@ -29,7 +29,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,
         )
 
 

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())
@@ -194,16 +210,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
 
 
 def _gate_from_name(gate_name: str) -> SU2CliffordT: