Add support for the Matsumoto-Amano normal form. (#1801)

Fixes #1762

Author: codrut3

qualtran/rotation_synthesis/matrix/_clifford_t_repr.py

@@ -16,6 +16,7 @@
 from typing import cast, Mapping, Optional
 
 import cirq
+import numpy as np
 
 import qualtran.rotation_synthesis.matrix._generation as ctg
 import qualtran.rotation_synthesis.matrix._su2_ct as _su2_ct
@@ -79,16 +80,72 @@ def _xz_sequence(
     return None
 
 
+def _matsumoto_amano_syllable(matrix: _su2_ct.SU2CliffordT) -> list[str]:
+    """Computes the next syllable in the Matsumoto-Amano decomposition for matrix.
+
+    Returns:
+        the next syllable as a list of gates, representing either T, HT, or SHT.
+
+    Raises:
+        ValueError if the parity matrix doesn't match any of the forms in
+        Lemma 4.10, https://arxiv.org/abs/1312.6584.
+    """
+    parity = matrix.bloch_form_parity()
+    # Parity matrix must have a 0 column, see Lemma 4.10, https://arxiv.org/abs/1312.6584.
+    # We move it to be last.
+    for i in range(2):
+        if np.all(parity[:, i] == 0):
+            parity[:, [i, 2]] = parity[:, [2, i]]
+            break
+    if np.array_equal(parity, np.array([[1, 1, 0], [1, 1, 0], [0, 0, 0]])):
+        # Leftmost syllable is T
+        return ['T']
+    elif np.array_equal(parity, np.array([[0, 0, 0], [1, 1, 0], [1, 1, 0]])):
+        # Leftmost syllable is HT
+        return ['H', 'T']
+    elif np.array_equal(parity, np.array([[1, 1, 0], [0, 0, 0], [1, 1, 0]])):
+        # Leftmost syllable is SHT
+        return ['S', 'H', 'T']
+    else:
+        raise ValueError(f'Unexpected parity matrix:\n{parity}')
+
+
+def _matsumoto_amano_sequence(matrix: _su2_ct.SU2CliffordT) -> tuple[str, ...]:
+    r"""Represents the Clifford+T operator in the Matsumoto-Amano normal form.
+
+    Returns:
+        a list of gates matching the regular expression $(T|\eps)(HT|SHT)^*C$,
+        where C is a Clifford operator, itself represented as a list of H and S gates.
+        The gates are returned in the order they need to be applied to generate the
+        input matrix.
+
+    Raises:
+        ValueError if during the decomposition an invalid SU2CliffordT matrix is created.
+    """
+    gates = []
+    while matrix.det() != _TWO:
+        syllable = _matsumoto_amano_syllable(matrix)
+        gates += syllable
+        for gate in syllable:
+            matrix = _su2_ct.GATE_MAP[gate].adjoint() @ matrix
+        new = matrix.scale_down()
+        if new is None or not new.is_valid():
+            raise ValueError('Invalid SU2CliffordT matrix')
+        matrix = new
+    return clifford(matrix) + tuple(gates[::-1])
+
+
 def to_sequence(matrix: _su2_ct.SU2CliffordT, fmt: str = 'xyz') -> tuple[str, ...]:
     r"""Returns a sequence of Clifford+T that produces the given matrix.
 
     Allowable format strings are
      - 'xyz' uses Tx, Ty, Tz gates.
      - 'xz' uses $Tx, Tz, Tx^\dagger, Tz^\dagger$
+     - 't' uses only Tz gates, and returns the Matsumoto-Amano normal form.
 
     Args:
         matrix: The matrix to represent.
-        fmt: A string from the set {'xz', 'xyz'} representing the allowed T gates described above.
+        fmt: A string from the set {'xz', 'xyz', 't'} representing the allowed T gates described above.
 
     Returns:
         A tuple of strings representing the gates.
@@ -100,6 +157,8 @@ def to_sequence(matrix: _su2_ct.SU2CliffordT, fmt: str = 'xyz') -> tuple[str, ..
         return _xyz_sequence(matrix)
     if fmt == 'xz':
         return cast(tuple[str, ...], _xz_sequence(matrix))
+    if fmt == 't':
+        return _matsumoto_amano_sequence(matrix)
     raise ValueError(f'{type=} is not supported')
 
 
@@ -116,6 +175,7 @@ def to_sequence(matrix: _su2_ct.SU2CliffordT, fmt: str = 'xyz') -> tuple[str, ..
     "Tx*": cirq.rx(-math.pi / 4),
     "Ty*": cirq.ry(-math.pi / 4),
     "Tz*": cirq.rz(-math.pi / 4),
+    "T": cirq.rz(math.pi / 4),
 }
 
 
@@ -130,6 +190,8 @@ def _to_quirk_name(name: str, allow_global_phase: bool = False) -> str:
         if name == "S*":
             return "\"Z^-½\""
         if name.startswith("T"):
+            if name == "T":
+                return "\"Z^¼\""
             if name.endswith("*"):
                 return "\"" + name[1].upper() + "^-¼" + "\""
             return "\"" + name[1].upper() + "^¼" + "\""
@@ -145,6 +207,8 @@ def _to_quirk_name(name: str, allow_global_phase: bool = False) -> str:
         if name == "S*":
             return '{"id":"Rzft","arg":"-pi/2"}'
         if name.startswith("T"):
+            if name == "T":
+                return '{"id":"Rzft","arg":"pi/4"}'
             angle = ['pi/4', '-pi/4'][name.endswith('*')]
             return '{"id":"R%sft","arg":"%s"}' % (name[1].lower(), angle)
     raise ValueError(f"{name=} is not supported")

qualtran/rotation_synthesis/matrix/_su2_ct.py

@@ -242,6 +242,62 @@ def num_t_gates(self) -> int:
         assert x == det
         return n
 
+    def bloch_sphere_form(self) -> tuple[np.ndarray, int]:
+        r"""Represents the unitary operator as a scaled element of SO(3).
+
+        The entries of the returned matrix are in $\mathbb{Z}[\sqrt{2}]$,
+        scaled by $\frac{1}{2\sqrt{2}(2+\sqrt{2})^n}$. The scaling factor
+        is implicit. The additional $\sqrt{2}$ factor is needed to ensure
+        the entries are in $\mathbb{Z}[\sqrt{2}]$.
+
+        See Section 4 in https://arxiv.org/abs/1312.6584.
+
+        Returns:
+            the scaled SO(3) matrix and the value n.
+        """
+        u = self.matrix[0, 0]
+        v = self.matrix[1, 0]
+        return (
+            np.array(
+                [
+                    [
+                        (u**2 - v.conj() ** 2).real_zsqrt2(),
+                        (u**2 - v**2).imag_zsqrt2(),
+                        2 * (u * v.conj()).real_zsqrt2(),
+                    ],
+                    [
+                        -(u**2 - v.conj() ** 2).imag_zsqrt2(),
+                        (u**2 + v**2).real_zsqrt2(),
+                        -2 * (u * v.conj()).imag_zsqrt2(),
+                    ],
+                    [
+                        -2 * (u * v).real_zsqrt2(),
+                        -2 * (u * v).imag_zsqrt2(),
+                        (u * u.conj() - v * v.conj()).real_zsqrt2(),
+                    ],
+                ]
+            ),
+            self.num_t_gates(),
+        )
+
+    def bloch_form_parity(self) -> np.ndarray:
+        """Returns the n-parity of the SO(3) Bloch sphere representation.
+
+        See Definition 4.7, https://arxiv.org/abs/1312.6584 for the definition of k-parity.
+        From Lemma 4.10, the least denominator exponent for the Bloch form equals the T-count,
+        and so the n-parity is well-defined.
+        """
+        bf, n = self.bloch_sphere_form()
+        bf = bf * _zsqrt2.SQRT_2**n
+        scale_factor = _zsqrt2.ZSqrt2(2, 0) * _zsqrt2.SQRT_2 * _zsqrt2.LAMBDA_KLIUCHNIKOV**n
+        if not all(x.is_divisible_by(scale_factor) for x in bf.flat):
+            raise ValueError(
+                "Not all entries of the \\sqrt{2}^n * SO(3) matrix are divisible "
+                f"by 2\\sqrt{2}(2+\\sqrt{2})^n. The matrix is:\n{bf}"
+            )
+        scaled_bf = [[x // scale_factor for x in row] for row in bf]
+        return np.array([[x.a % 2 for x in row] for row in scaled_bf])
+
 
 def _gate_from_name(gate_name: str) -> SU2CliffordT:
     adjoint = gate_name.count('*') % 2 == 1
@@ -308,6 +364,7 @@ def _key_map():
 # Tx, Ty, Tz scaled by sqrt(2*(2+sqrt(2)))
 Tx = SU2CliffordT(_I * _zw.SQRT_2 + _I - _X * _zw.J, ("Tx",))
 Ty = SU2CliffordT(_I * _zw.SQRT_2 + _I - _Y * _zw.J, ("Ty",))
+# Tz = sqrt(2) * (1 + w^*) * T and has determinant = 2(2+sqrt(2))
 Tz = SU2CliffordT(_I * _zw.SQRT_2 + _I - _Z * _zw.J, ("Tz",))
 Ts = [Tx, Ty, Tz]
 
@@ -322,6 +379,7 @@ def _key_map():
     "X": XSqrt2,
     "Y": YSqrt2,
     "Z": ZSqrt2,
+    "T": Tz,
 }
 
 PARAMETRIC_FORM_BASES = [

qualtran/rotation_synthesis/matrix/clifford_t_repr_test.py

@@ -60,9 +60,47 @@ def test_to_xz_seq(g: _su2_ct.SU2CliffordT):
 
 
 @pytest.mark.parametrize("g", _make_random_su(50, 5, random_cliffords=True, seed=0))
-def test_to_cirq(g):
-    circuit = cirq.Circuit(ctr.to_cirq(g, 'xyz'))
+@pytest.mark.parametrize("fmt", ('xyz', 'xz', 't'))
+def test_to_cirq(g: _su2_ct.SU2CliffordT, fmt: str):
+    g = g.rescale()
+    circuit = cirq.Circuit(ctr.to_cirq(g, fmt))
     unitary = cirq.unitary(circuit)
     u = g.matrix.astype(complex)
     u = u / np.linalg.det(u) ** 0.5
     assert cirq.protocols.equal_up_to_global_phase(u, unitary)
+
+
+@pytest.mark.parametrize(
+    ["g", "fmt", "expected"],
+    [
+        [_su2_ct.HSqrt2, "t", ('"Z^½"', '"X"', '"Z^½"', '"X"', '"H"')],
+        [_su2_ct.SSqrt2, "t", ('{"id":"Rzft","arg":"pi/2"}',)],
+        [_su2_ct.Tz, "t", ('{"id":"Rzft","arg":"pi/4"}',)],
+    ],
+)
+def test_to_quirk(g: _su2_ct.SU2CliffordT, fmt: str, expected: tuple[str, ...]):
+    assert ctr.to_quirk(g, fmt) == expected
+
+
+@pytest.mark.parametrize("g", _make_random_su(50, 10, random_cliffords=True, seed=0))
+def test_to_matsumoto_amano_seq(g: _su2_ct.SU2CliffordT):
+    g = g.rescale()
+    seq = ctr.to_sequence(g, 't')
+    prev_t = -1
+    # Check that the reversed list matches the regular expression $(T|\eps)(HT|SHT)^*C$.
+    seq_r = list(reversed(seq))
+    for i in range(len(seq_r)):
+        assert seq_r[i] in ('T', 'H', 'S', 'Z', 'X')
+        if i == 0 and seq_r[0] == 'T':
+            prev_t = 0
+            continue
+        if seq_r[i] == 'T':
+            # Check that this is one of HT, SHT syllabes
+            assert prev_t in (i - 2, i - 3)
+            if prev_t == i - 2:
+                assert seq_r[i - 1] == 'H'
+            else:
+                assert seq_r[i - 2] + seq_r[i - 1] == 'SH'
+            prev_t = i
+    got = _su2_ct.SU2CliffordT.from_sequence(seq)
+    assert got == g

qualtran/rotation_synthesis/matrix/su2_ct_test.py

@@ -113,3 +113,55 @@ def test_num_t_gates():
 
         # We need to call .rescale to remove the common factor and reduce the T count.
         assert (t @ t.adjoint()).rescale().num_t_gates() == 0
+
+
+_H_bloch_form_numpy = np.array([[0, 0, 1], [0, -1, 0], [1, 0, 0]])
+_S_bloch_form_numpy = np.array([[0, -1, 0], [1, 0, 0], [0, 0, 1]])
+_T_bloch_form_numpy = np.array([[1, -1, 0], [1, 1, 0], [0, 0, _SQRT2]]) / _SQRT2
+
+
+@pytest.mark.parametrize(
+    ["g", "bloch_form", "n"],
+    [
+        [_su2_ct.HSqrt2, _H_bloch_form_numpy, 0],
+        [_su2_ct.SSqrt2, _S_bloch_form_numpy, 0],
+        [_su2_ct.Tz, _T_bloch_form_numpy, 1],
+    ],
+)
+def test_bloch_sphere_form_generators(g: _su2_ct.SU2CliffordT, bloch_form: np.ndarray, n: int):
+    g_so3, m = g.bloch_sphere_form()
+    assert m == n
+    np.testing.assert_allclose(g_so3.astype(float) / (2 * _SQRT2 * (2 + _SQRT2) ** m), bloch_form)
+
+
+def _matrix_from_pauli(coords: np.ndarray) -> np.ndarray:
+    return coords[0] * _X + coords[1] * _Y + coords[2] * _Z
+
+
+@pytest.mark.parametrize("g", _make_random_su(10, 5, random_cliffords=True, seed=0))
+@pytest.mark.parametrize("vector", np.random.choice(2, size=(10, 3)))
+def test_bloch_sphere_form_random(g: _su2_ct.SU2CliffordT, vector: np.ndarray):
+    g_so3, _ = g.bloch_sphere_form()
+    g_action = (
+        g.matrix.astype(complex) @ _matrix_from_pauli(vector) @ g.adjoint().matrix.astype(complex)
+    )
+    g_so3_action = g_so3.astype(float) @ vector.T / _SQRT2
+    np.testing.assert_allclose(_matrix_from_pauli(g_so3_action), g_action, atol=1e-7)
+
+
+_T_parity_numpy = np.array([[1, 1, 0], [1, 1, 0], [0, 0, 0]])
+_HT_parity_numpy = np.array([[0, 0, 0], [1, 1, 0], [1, 1, 0]])
+_SHT_parity_numpy = np.array([[1, 1, 0], [0, 0, 0], [1, 1, 0]])
+
+
+@pytest.mark.parametrize(
+    ["g", "parity"],
+    [
+        [_su2_ct.Tz, _T_parity_numpy],
+        [_su2_ct.HSqrt2 @ _su2_ct.Tz, _HT_parity_numpy],
+        [_su2_ct.SSqrt2 @ _su2_ct.HSqrt2 @ _su2_ct.Tz, _SHT_parity_numpy],
+    ],
+)
+def test_bloch_form_parity(g: _su2_ct.SU2CliffordT, parity: np.ndarray):
+    g_parity = g.bloch_form_parity()
+    np.testing.assert_equal(g_parity, parity)

qualtran/rotation_synthesis/rings/_zw.py

@@ -42,7 +42,7 @@ class ZW:
 
     Elements of $\mathbb{Z}[\omega]$ can be represented by 4 integers as
     $$
-        \sum_{i=0}^4 a_i \omega^i
+        \sum_{i=0}^3 a_i \omega^i
     $$
 
     Attributes:
@@ -206,6 +206,16 @@ def to_zsqrt2(self) -> tuple[_zsqrt2.ZSqrt2, _zsqrt2.ZSqrt2, bool]:
             r == 1,
         )
 
+    def real_zsqrt2(self) -> _zsqrt2.ZSqrt2:
+        r"""Returns \sqrt{2} * the real part of the element."""
+        a, _, need_w = self.to_zsqrt2()
+        return a * _zsqrt2.SQRT_2 + _zsqrt2.ZSqrt2(need_w, 0)
+
+    def imag_zsqrt2(self) -> _zsqrt2.ZSqrt2:
+        r"""Returns \sqrt{2} * the imaginary part of the element."""
+        _, b, need_w = self.to_zsqrt2()
+        return b * _zsqrt2.SQRT_2 + _zsqrt2.ZSqrt2(need_w, 0)
+
     def norm(self) -> int:
         res = self * self.conj() * self.sqrt2_conj() * self.sqrt2_conj().conj()
         assert all(c == 0 for c in res.coords[1:])

qualtran/rotation_synthesis/rings/zw_test.py

@@ -138,3 +138,13 @@ def test_factor_prime(p):
     else:
         assert r.coords[1:] == (0, 0, 0)
         assert r.coords[0] % p == 0
+
+
+@pytest.mark.parametrize("x", _create_random_elements(20))
+def test_real_zsqrt2(x: rings.ZW):
+    np.testing.assert_allclose((x.real_zsqrt2()).value(_SQRT_2), x.value(_SQRT_2).real * _SQRT_2)
+
+
+@pytest.mark.parametrize("x", _create_random_elements(20))
+def test_imag_zsqrt2(x: rings.ZW):
+    np.testing.assert_allclose((x.imag_zsqrt2()).value(_SQRT_2), x.value(_SQRT_2).imag * _SQRT_2)