Add Clifford Tableau decomposition function

cirq-core/cirq/__init__.py

@@ -303,6 +303,7 @@
     AlignRight,
     compute_cphase_exponents_for_fsim_decomposition,
     ConvertToCzAndSingleGates,
+    decompose_clifford_tableau_to_operations,
     decompose_cphase_into_two_fsim,
     decompose_multi_controlled_x,
     decompose_multi_controlled_rotation,

cirq-core/cirq/ops/common_gates.py

@@ -1073,8 +1073,11 @@ def _act_on_(self, args: Any):
                     tableau.zs[:, q2].copy(),
                     tableau.xs[:, q2].copy(),
                 )
+                tableau.rs[:] ^= tableau.xs[:, q2] & tableau.zs[:, q2]
                 tableau.rs[:] ^= (
-                    tableau.xs[:, q1] & tableau.zs[:, q2] & (tableau.xs[:, q2] ^ tableau.zs[:, q1])
+                    tableau.xs[:, q1]
+                    & tableau.zs[:, q2]
+                    & (~(tableau.xs[:, q2] ^ tableau.zs[:, q1]))
                 )
                 tableau.xs[:, q2] ^= tableau.xs[:, q1]
                 tableau.zs[:, q1] ^= tableau.zs[:, q2]

cirq-core/cirq/optimizers/__init__.py

@@ -22,6 +22,8 @@
     AlignRight,
 )
 
+from cirq.optimizers.clifford_decomposition import decompose_clifford_tableau_to_operations
+
 from cirq.optimizers.cphase_to_fsim import (
     compute_cphase_exponents_for_fsim_decomposition,
     decompose_cphase_into_two_fsim,

cirq-core/cirq/optimizers/clifford_decomposition.py

@@ -0,0 +1,151 @@
+# Copyright 2021 The Cirq Developers
+#
+# 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.
+
+"""Utility methods related to decompose clifford gate into circuits."""
+
+from typing import List, TYPE_CHECKING
+import functools
+
+import numpy as np
+from cirq import ops, protocols, qis, sim
+
+if TYPE_CHECKING:
+    import cirq
+
+
+def _X(q, args, operations, qubits):
+    args.axes = [q]
+    protocols.act_on(ops.X, args, allow_decompose=False)
+    operations.append(ops.X(qubits[q]))
+
+
+def _Z(q, args, operations, qubits):
+    args.axes = [q]
+    protocols.act_on(ops.Z, args, allow_decompose=False)
+    operations.append(ops.Z(qubits[q]))
+
+
+def _Sdg(q, args, operations, qubits):
+    # Apply the tableau with S^+, so the inverse of operation is S
+    args.axes = [q]
+    protocols.act_on(ops.ZPowGate() ** 1.5, args, allow_decompose=False)
+    operations.append(ops.S(qubits[q]))
+
+
+def _H(q, args, operations, qubits):
+    args.axes = [q]
+    protocols.act_on(ops.H, args, allow_decompose=False)
+    operations.append(ops.H(qubits[q]))
+
+
+def _CNOT(q1, q2, args, operations, qubits):
+    args.axes = [q1, q2]
+    protocols.act_on(ops.CNOT, args, allow_decompose=False)
+    operations.append(ops.CNOT(qubits[q1], qubits[q2]))
+
+
+def _SWAP(q1, q2, args, operations, qubits):
+    args.axes = [q1, q2]
+    protocols.act_on(ops.SWAP, args, allow_decompose=False)
+    operations.append(ops.SWAP(qubits[q1], qubits[q2]))
+
+
+def decompose_clifford_tableau_to_operations(
+    qubits: List['cirq.Qid'], clifford_tableau: qis.CliffordTableau
+) -> List[ops.Operation]:
+    """Decompose an n-qubit Clifford Tableau into a list of one/two qubit operations.
+
+    Args:
+        qubits: The list of qubits being operated on.
+        clifford_tableau: The Clifford Tableau for decomposition.
+
+    Returns:
+        A list of operations reconstructs the same Clifford tableau.
+    """
+    if len(qubits) != clifford_tableau.n:
+        raise ValueError(
+            f"The number of qubits must be the same as the number of Clifford Tableau."
+        )
+    assert (
+        clifford_tableau._validate()
+    ), "The provided clifford_tableau must satisfy the symplectic property."
+
+    t: qis.CliffordTableau = clifford_tableau.copy()
+    operations: List[ops.Operation] = []
+    args = sim.ActOnCliffordTableauArgs(
+        tableau=t, axes=[], prng=np.random.RandomState(), log_of_measurement_results={}
+    )
+
+    _X_with_ops = functools.partial(_X, args=args, operations=operations, qubits=qubits)
+    _Z_with_ops = functools.partial(_Z, args=args, operations=operations, qubits=qubits)
+    _H_with_ops = functools.partial(_H, args=args, operations=operations, qubits=qubits)
+    _S_with_ops = functools.partial(_Sdg, args=args, operations=operations, qubits=qubits)
+    _CNOT_with_ops = functools.partial(_CNOT, args=args, operations=operations, qubits=qubits)
+    _SWAP_with_ops = functools.partial(_SWAP, args=args, operations=operations, qubits=qubits)
+
+    # The procedure is based on Theorem 8 in
+    # [1] S. Aaronson, D. Gottesman, *Improved Simulation of Stabilizer Circuits*,
+    #     Phys. Rev. A 70, 052328 (2004). https://arxiv.org/abs/quant-ph/0406196
+    # with modification by doing it row-by-row instead.
+
+    # Suppose we have a Clifford Tableau:
+    #                   Xs  Zs
+    # Destabilizers:  [ A | B ]
+    # Stabilizers:    [ C | D ]
+    for i in range(t.n):
+        # Step 1a: Make the diagonal element of A as 1 by Hadamard gate if necessary.
+        if not t.xs[i, i] and t.zs[i, i]:
+            _H_with_ops(i)
+        # Step 1b: Make the diagonal element of A as 1 by swapping gate if necessary.
+        if not t.xs[i, i]:
+            for j in range(i + 1, t.n):
+                if t.xs[i, j]:
+                    _SWAP_with_ops(i, j)
+                    break
+        # Step 1c: We may still not be able to find non-zero element in whole Xs row. Then,
+        # apply swap + Hadamard from zs. It is guaranteed to find one by lemma 5 in [1].
+        if not t.xs[i, i]:
+            for j in range(i + 1, t.n):
+                if t.zs[i, j]:
+                    _H_with_ops(j)
+                    _SWAP_with_ops(i, j)
+                    break
+
+        # Step 2: Eliminate the elements in A By CNOT and phase gate (i-th row)
+        # first i rows of destabilizers: [ I  0 | 0  0 ]
+        _ = [_CNOT_with_ops(i, j) for j in range(i + 1, t.n) if t.xs[i, j]]
+        if np.any(t.zs[i, i:]):
+            if not t.zs[i, i]:
+                _S_with_ops(i)
+            _ = [_CNOT_with_ops(j, i) for j in range(i + 1, t.n) if t.zs[i, j]]
+            _S_with_ops(i)
+
+        # Step 3: Eliminate the elements in D By CNOT and phase gate (i-th row)
+        # first i rows of stabilizers: [ 0  0 | I  0 ]
+        _ = [_CNOT_with_ops(j, i) for j in range(i + 1, t.n) if t.zs[i + t.n, j]]
+        if np.any(t.xs[i + t.n, i:]):
+            # Swap xs and zs
+            _H_with_ops(i)
+            _ = [_CNOT_with_ops(i, j) for j in range(i + 1, t.n) if t.xs[i + t.n, j]]
+            if t.zs[i + t.n, i]:
+                _S_with_ops(i)
+            _H_with_ops(i)
+
+    # Step 4: Correct the phase of tableau
+    _ = [_Z_with_ops(i) for i, p in enumerate(t.rs[: t.n]) if p]
+    _ = [_X_with_ops(i) for i, p in enumerate(t.rs[t.n :]) if p]
+
+    # Step 5: invert the operations by reversing the orde: (AB)^{+} = B^{+} A^{+}.
+    # Note only S gate is not self-adjoint.
+    return operations[::-1]

cirq-core/cirq/optimizers/clifford_decomposition_test.py

@@ -0,0 +1,174 @@
+# Copyright 2021 The Cirq Developers
+#
+# 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 pytest
+import numpy as np
+
+import cirq
+from cirq.testing import assert_allclose_up_to_global_phase
+
+
+def test_misaligned_qubits():
+    qubits = cirq.LineQubit.range(1)
+    tableau = cirq.CliffordTableau(num_qubits=2)
+    with pytest.raises(ValueError):
+        cirq.decompose_clifford_tableau_to_operations(qubits, tableau)
+
+
+def test_clifford_decompose_one_qubit():
+    """Two random instance for one qubit decomposition."""
+    qubits = cirq.LineQubit.range(1)
+    args = cirq.ActOnCliffordTableauArgs(
+        tableau=cirq.CliffordTableau(num_qubits=1),
+        axes=[0],
+        prng=np.random.RandomState(),
+        log_of_measurement_results={},
+    )
+    cirq.act_on(cirq.X, args, allow_decompose=False)
+    cirq.act_on(cirq.H, args, allow_decompose=False)
+    cirq.act_on(cirq.S, args, allow_decompose=False)
+    expect_circ = cirq.Circuit(cirq.X(qubits[0]), cirq.H(qubits[0]), cirq.S(qubits[0]))
+    ops = cirq.decompose_clifford_tableau_to_operations(qubits, args.tableau)
+    circ = cirq.Circuit(ops)
+    assert_allclose_up_to_global_phase(cirq.unitary(expect_circ), cirq.unitary(circ), atol=1e-7)
+
+    qubits = cirq.LineQubit.range(1)
+    args = cirq.ActOnCliffordTableauArgs(
+        tableau=cirq.CliffordTableau(num_qubits=1),
+        axes=[0],
+        prng=np.random.RandomState(),
+        log_of_measurement_results={},
+    )
+    cirq.act_on(cirq.Z, args, allow_decompose=False)
+    cirq.act_on(cirq.H, args, allow_decompose=False)
+    cirq.act_on(cirq.S, args, allow_decompose=False)
+    cirq.act_on(cirq.H, args, allow_decompose=False)
+    cirq.act_on(cirq.X, args, allow_decompose=False)
+    expect_circ = cirq.Circuit(
+        cirq.Z(qubits[0]),
+        cirq.H(qubits[0]),
+        cirq.S(qubits[0]),
+        cirq.H(qubits[0]),
+        cirq.X(qubits[0]),
+    )
+    ops = cirq.decompose_clifford_tableau_to_operations(qubits, args.tableau)
+    circ = cirq.Circuit(ops)
+    assert_allclose_up_to_global_phase(cirq.unitary(expect_circ), cirq.unitary(circ), atol=1e-7)
+
+
+def test_clifford_decompose_two_qubits():
+    """Two random instance for two qubits decomposition."""
+    qubits = cirq.LineQubit.range(2)
+    args = cirq.ActOnCliffordTableauArgs(
+        tableau=cirq.CliffordTableau(num_qubits=2),
+        axes=[0],
+        prng=np.random.RandomState(),
+        log_of_measurement_results={},
+    )
+    cirq.act_on(cirq.H, args, allow_decompose=False)
+    args.axes = [0, 1]
+    cirq.act_on(cirq.CNOT, args, allow_decompose=False)
+    expect_circ = cirq.Circuit(cirq.H(qubits[0]), cirq.CNOT(qubits[0], qubits[1]))
+    ops = cirq.decompose_clifford_tableau_to_operations(qubits, args.tableau)
+    circ = cirq.Circuit(ops)
+    assert_allclose_up_to_global_phase(cirq.unitary(expect_circ), cirq.unitary(circ), atol=1e-7)
+
+    qubits = cirq.LineQubit.range(2)
+    args = cirq.ActOnCliffordTableauArgs(
+        tableau=cirq.CliffordTableau(num_qubits=2),
+        axes=[0],
+        prng=np.random.RandomState(),
+        log_of_measurement_results={},
+    )
+    cirq.act_on(cirq.H, args, allow_decompose=False)
+    args.axes = [0, 1]
+    cirq.act_on(cirq.CNOT, args, allow_decompose=False)
+    args.axes = [0]
+    cirq.act_on(cirq.H, args, allow_decompose=False)
+    cirq.act_on(cirq.S, args, allow_decompose=False)
+    args.axes = [1]
+    cirq.act_on(cirq.X, args, allow_decompose=False)
+    expect_circ = cirq.Circuit(
+        cirq.H(qubits[0]),
+        cirq.CNOT(qubits[0], qubits[1]),
+        cirq.H(qubits[0]),
+        cirq.S(qubits[0]),
+        cirq.X(qubits[1]),
+    )
+
+    ops = cirq.decompose_clifford_tableau_to_operations(qubits, args.tableau)
+    circ = cirq.Circuit(ops)
+    assert_allclose_up_to_global_phase(cirq.unitary(expect_circ), cirq.unitary(circ), atol=1e-7)
+
+
+def test_clifford_decompose_by_unitary():
+    """Validate the decomposition of random Clifford Tableau by unitary matrix.
+
+    Due to the exponential growth in dimension, it cannot validate very large number of qubits.
+    """
+    n, num_ops = 5, 20
+    gate_candidate = [cirq.X, cirq.Y, cirq.Z, cirq.H, cirq.S, cirq.CNOT, cirq.CZ]
+    for seed in range(100):
+        prng = np.random.RandomState(seed)
+        t = cirq.CliffordTableau(num_qubits=n)
+        qubits = cirq.LineQubit.range(n)
+        expect_circ = cirq.Circuit()
+        args = cirq.ActOnCliffordTableauArgs(
+            tableau=t, axes=[], prng=prng, log_of_measurement_results={}
+        )
+        for _ in range(num_ops):
+            g = prng.randint(len(gate_candidate))
+            indices = (prng.randint(n),) if g < 5 else prng.choice(n, 2, replace=False)
+            args.axes = indices
+            cirq.act_on(gate_candidate[g], args, allow_decompose=False)
+            expect_circ.append(gate_candidate[g].on(*[qubits[i] for i in indices]))
+        ops = cirq.decompose_clifford_tableau_to_operations(qubits, args.tableau)
+        circ = cirq.Circuit(ops)
+        circ.append(cirq.I.on_each(qubits))
+        expect_circ.append(cirq.I.on_each(qubits))
+        assert_allclose_up_to_global_phase(cirq.unitary(expect_circ), cirq.unitary(circ), atol=1e-7)
+
+
+def test_clifford_decompose_by_reconstruction():
+    """Validate the decomposition of random Clifford Tableau by reconstruction.
+
+    This approach can validate large number of qubits compared with the unitary one.
+    """
+    n, num_ops = 100, 500
+    gate_candidate = [cirq.X, cirq.Y, cirq.Z, cirq.H, cirq.S, cirq.CNOT, cirq.CZ]
+    for seed in range(10):
+        prng = np.random.RandomState(seed)
+        t = cirq.CliffordTableau(num_qubits=n)
+        qubits = cirq.LineQubit.range(n)
+        expect_circ = cirq.Circuit()
+        args = cirq.ActOnCliffordTableauArgs(
+            tableau=t, axes=[], prng=prng, log_of_measurement_results={}
+        )
+        for _ in range(num_ops):
+            g = prng.randint(len(gate_candidate))
+            indices = (prng.randint(n),) if g < 5 else prng.choice(n, 2, replace=False)
+            args.axes = indices
+            cirq.act_on(gate_candidate[g], args, allow_decompose=False)
+            expect_circ.append(gate_candidate[g].on(*[qubits[i] for i in indices]))
+        ops = cirq.decompose_clifford_tableau_to_operations(qubits, args.tableau)
+
+        reconstruct_t = cirq.CliffordTableau(num_qubits=n)
+        reconstruct_args = cirq.ActOnCliffordTableauArgs(
+            tableau=reconstruct_t, axes=[], prng=prng, log_of_measurement_results={}
+        )
+        for op in ops:
+            reconstruct_args.axes = [q.x for q in op.qubits]
+            cirq.act_on(op.gate, reconstruct_args, allow_decompose=False)
+
+        assert t == reconstruct_t