diff --git a/cirq-core/cirq/__init__.py b/cirq-core/cirq/__init__.py index 847b29ef163..39490bdc8a5 100644 --- a/cirq-core/cirq/__init__.py +++ b/cirq-core/cirq/__init__.py @@ -329,6 +329,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, diff --git a/cirq-core/cirq/optimizers/__init__.py b/cirq-core/cirq/optimizers/__init__.py index f3c739f90db..9a1ab286dca 100644 --- a/cirq-core/cirq/optimizers/__init__.py +++ b/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, diff --git a/cirq-core/cirq/optimizers/clifford_decomposition.py b/cirq-core/cirq/optimizers/clifford_decomposition.py new file mode 100644 index 00000000000..557676b2a3a --- /dev/null +++ b/cirq-core/cirq/optimizers/clifford_decomposition.py @@ -0,0 +1,184 @@ +# 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 to decompose Clifford gates 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: int, + args: sim.ActOnCliffordTableauArgs, + operations: List[ops.Operation], + qubits: List['cirq.Qid'], +): + protocols.act_on(ops.X, args, qubits=[qubits[q]], allow_decompose=False) + operations.append(ops.X(qubits[q])) + + +def _Z( + q: int, + args: sim.ActOnCliffordTableauArgs, + operations: List[ops.Operation], + qubits: List['cirq.Qid'], +): + protocols.act_on(ops.Z, args, qubits=[qubits[q]], allow_decompose=False) + operations.append(ops.Z(qubits[q])) + + +def _Sdg( + q: int, + args: sim.ActOnCliffordTableauArgs, + operations: List[ops.Operation], + qubits: List['cirq.Qid'], +): + # Apply the tableau with S^\{dagger} + protocols.act_on(ops.S ** -1, args, qubits=[qubits[q]], allow_decompose=False) + operations.append(ops.S(qubits[q])) + + +def _H( + q: int, + args: sim.ActOnCliffordTableauArgs, + operations: List[ops.Operation], + qubits: List['cirq.Qid'], +): + protocols.act_on(ops.H, args, qubits=[qubits[q]], allow_decompose=False) + operations.append(ops.H(qubits[q])) + + +def _CNOT( + q1: int, + q2: int, + args: sim.ActOnCliffordTableauArgs, + operations: List[ops.Operation], + qubits: List['cirq.Qid'], +): + protocols.act_on(ops.CNOT, args, qubits=[qubits[q1], qubits[q2]], allow_decompose=False) + operations.append(ops.CNOT(qubits[q1], qubits[q2])) + + +def _SWAP( + q1: int, + q2: int, + args: sim.ActOnCliffordTableauArgs, + operations: List[ops.Operation], + qubits: List['cirq.Qid'], +): + protocols.act_on(ops.SWAP, args, qubits=[qubits[q1], qubits[q2]], 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. + + The implementation is based on Theorem 8 in [1]. + [1] S. Aaronson, D. Gottesman, *Improved Simulation of Stabilizer Circuits*, + Phys. Rev. A 70, 052328 (2004). https://arxiv.org/abs/quant-ph/0406196 + + 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. + + Raises: + ValueError: The length of input qubit mismatch with the size of 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, qubits=qubits, 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 equal to 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 equal to 1 by SWAP 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 order: (AB)^{+} = B^{+} A^{+}. + # Note only S gate is not self-adjoint. + return operations[::-1] diff --git a/cirq-core/cirq/optimizers/clifford_decomposition_test.py b/cirq-core/cirq/optimizers/clifford_decomposition_test.py new file mode 100644 index 00000000000..64ba424ce51 --- /dev/null +++ b/cirq-core/cirq/optimizers/clifford_decomposition_test.py @@ -0,0 +1,171 @@ +# 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), + qubits=qubits, + prng=np.random.RandomState(), + log_of_measurement_results={}, + ) + cirq.act_on(cirq.X, args, qubits=[qubits[0]], allow_decompose=False) + cirq.act_on(cirq.H, args, qubits=[qubits[0]], allow_decompose=False) + cirq.act_on(cirq.S, args, qubits=[qubits[0]], 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), + qubits=qubits, + prng=np.random.RandomState(), + log_of_measurement_results={}, + ) + cirq.act_on(cirq.Z, args, qubits=[qubits[0]], allow_decompose=False) + cirq.act_on(cirq.H, args, qubits=[qubits[0]], allow_decompose=False) + cirq.act_on(cirq.S, args, qubits=[qubits[0]], allow_decompose=False) + cirq.act_on(cirq.H, args, qubits=[qubits[0]], allow_decompose=False) + cirq.act_on(cirq.X, args, qubits=[qubits[0]], 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), + qubits=qubits, + prng=np.random.RandomState(), + log_of_measurement_results={}, + ) + cirq.act_on(cirq.H, args, qubits=[qubits[0]], allow_decompose=False) + cirq.act_on(cirq.CNOT, args, qubits=[qubits[0], qubits[1]], 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), + qubits=qubits, + prng=np.random.RandomState(), + log_of_measurement_results={}, + ) + cirq.act_on(cirq.H, args, qubits=[qubits[0]], allow_decompose=False) + cirq.act_on(cirq.CNOT, args, qubits=[qubits[0], qubits[1]], allow_decompose=False) + cirq.act_on(cirq.H, args, qubits=[qubits[0]], allow_decompose=False) + cirq.act_on(cirq.S, args, qubits=[qubits[0]], allow_decompose=False) + cirq.act_on(cirq.X, args, qubits=[qubits[1]], 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, qubits=qubits, 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) + cirq.act_on( + gate_candidate[g], args, qubits=[qubits[i] for i in indices], 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, qubits=qubits, 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) + cirq.act_on( + gate_candidate[g], args, qubits=[qubits[i] for i in indices], 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, qubits=qubits, prng=prng, log_of_measurement_results={} + ) + for op in ops: + cirq.act_on(op.gate, reconstruct_args, qubits=op.qubits, allow_decompose=False) + + assert t == reconstruct_t