diff --git a/cirq-core/cirq/transformers/analytical_decompositions/quantum_shannon_decomposition.py b/cirq-core/cirq/transformers/analytical_decompositions/quantum_shannon_decomposition.py index 125b8aa5f72..085e59db297 100644 --- a/cirq-core/cirq/transformers/analytical_decompositions/quantum_shannon_decomposition.py +++ b/cirq-core/cirq/transformers/analytical_decompositions/quantum_shannon_decomposition.py @@ -27,22 +27,40 @@ from cirq import ops from cirq.linalg import decompositions, predicates +from cirq.transformers.analytical_decompositions.two_qubit_to_cz import ( + two_qubit_matrix_to_cz_operations, +) +from cirq.transformers.analytical_decompositions.three_qubit_decomposition import ( + three_qubit_matrix_to_operations, +) +from cirq.protocols import unitary_protocol +from cirq.circuits.frozen_circuit import FrozenCircuit if TYPE_CHECKING: import cirq from cirq.ops import op_tree -def quantum_shannon_decomposition(qubits: 'List[cirq.Qid]', u: np.ndarray) -> 'op_tree.OpTree': - """Decomposes n-qubit unitary into CX/YPow/ZPow/CNOT gates, preserving global phase. +def quantum_shannon_decomposition( + qubits: 'List[cirq.Qid]', u: np.ndarray, atol: float = 1e-8 +) -> 'op_tree.OpTree': + """Decomposes n-qubit unitary 1-q, 2-q and GlobalPhase gates, preserving global phase. + + The gates used are CX/YPow/ZPow/CNOT/GlobalPhase/CZ/PhasedXZGate/PhasedXPowGate. The algorithm is described in Shende et al.: Synthesis of Quantum Logic Circuits. Tech. rep. 2006, https://arxiv.org/abs/quant-ph/0406176 + Note: Shannon decomposition is sensitive to the numerical accuracy of doing eigendecomposition. + Eigendecomposition is obtained using `np.linalg.eig` and the resulting difference between + the input and output unitary is heavily affected by the accuracy of `np.linalg.eig`. + + Args: qubits: List of qubits in order of significance u: Numpy array for unitary matrix representing gate to be decomposed + atol: Absolute tolerance of floating point checks. Calls: (Base Case) @@ -62,7 +80,7 @@ def quantum_shannon_decomposition(qubits: 'List[cirq.Qid]', u: np.ndarray) -> 'o ValueError: If the u matrix is non-unitary ValueError: If the u matrix is not of shape (2^n,2^n) """ - if not predicates.is_unitary(u): # Check that u is unitary + if not predicates.is_unitary(u, atol=atol): # Check that u is unitary raise ValueError( "Expected input matrix u to be unitary, \ but it fails cirq.is_unitary check" @@ -80,6 +98,32 @@ def quantum_shannon_decomposition(qubits: 'List[cirq.Qid]', u: np.ndarray) -> 'o yield from _single_qubit_decomposition(qubits[0], u) return + if n == 4: + operations = tuple( + two_qubit_matrix_to_cz_operations( + qubits[0], qubits[1], u, allow_partial_czs=True, clean_operations=True, atol=atol + ) + ) + yield from operations + i, j = np.unravel_index(np.argmax(np.abs(u)), u.shape) + new_unitary = unitary_protocol.unitary(FrozenCircuit.from_moments(*operations)) + global_phase = np.angle(u[i, j]) - np.angle(new_unitary[i, j]) + if np.abs(global_phase) > 1e-9: + yield ops.global_phase_operation(np.exp(1j * global_phase)) + return + + if n == 8: + operations = tuple( + three_qubit_matrix_to_operations(qubits[0], qubits[1], qubits[2], u, atol=atol) + ) + yield from operations + i, j = np.unravel_index(np.argmax(np.abs(u)), u.shape) + new_unitary = unitary_protocol.unitary(FrozenCircuit.from_moments(*operations)) + global_phase = np.angle(u[i, j]) - np.angle(new_unitary[i, j]) + if np.abs(global_phase) > 1e-9: + yield ops.global_phase_operation(np.exp(1j * global_phase)) + return + # Perform a cosine-sine (linalg) decomposition on u # X = [ u1 , 0 ] [ cos(theta) , -sin(theta) ] [ v1 , 0 ] # [ 0 , u2 ] [ sin(theta) , cos(theta) ] [ 0 , v2 ] @@ -108,17 +152,36 @@ def _single_qubit_decomposition(qubit: 'cirq.Qid', u: np.ndarray) -> 'op_tree.Op A single operation from OP TREE of 3 operations (rz,ry,ZPowGate) """ # Perform native ZYZ decomposition - phi_0, phi_1, phi_2 = decompositions.deconstruct_single_qubit_matrix_into_angles(u) + phi_0, phi_1, phi_2 = np.array( + decompositions.deconstruct_single_qubit_matrix_into_angles(u) + ) % (2 * np.pi) # Determine global phase picked up - phase = np.angle(u[0, 0] / (np.exp(-1j * (phi_0) / 2) * np.cos(phi_1 / 2))) - - # Append first two operations operations - yield ops.rz(phi_0).on(qubit) - yield ops.ry(phi_1).on(qubit) - - # Append third operation with global phase added - yield ops.ZPowGate(exponent=phi_2 / np.pi, global_shift=phase / phi_2).on(qubit) + global_phase = np.angle(u[0, 0]) + phi_0 / 2 + phi_2 / 2 + if np.abs(u[0, 0]) < 1e-9: + global_phase = np.angle(u[1, 0]) + phi_0 / 2 - phi_2 / 2 + + if np.abs(phi_2) > 1e-18: + # Append first two operations operations + yield ops.rz(phi_0).on(qubit) + yield ops.ry(phi_1).on(qubit) + + # Append third operation with global phase added + yield ops.ZPowGate(exponent=phi_2 / np.pi, global_shift=global_phase / phi_2 - 0.5)(qubit) + elif np.abs(phi_1) > 1e-18: + # Just a Z -> Y rotation so we attach the global phase to the Y rotation. + if np.abs(phi_0) > 1e-18: + yield ops.rz(phi_0)(qubit) + yield ops.YPowGate(exponent=phi_1 / np.pi, global_shift=global_phase / phi_1 - 0.5)(qubit) + elif np.abs(phi_0) > 1e-18: + # Just an Rz with a potential global phase. + yield ops.ZPowGate(exponent=phi_0 / np.pi, global_shift=global_phase / phi_0 - 0.5)(qubit) + elif np.abs(global_phase) > 1e-18: + # Global Phase. + yield ops.global_phase_operation(np.exp(1j * global_phase)) + else: + # Identity. + return def _msb_demuxer( @@ -146,8 +209,21 @@ def _msb_demuxer( Yields: Single operation from OP TREE of 2-qubit and 1-qubit operations """ # Perform a diagonalization to find values + u1 = u1.astype(np.complex128) + u2 = u2.astype(np.complex128) u = u1 @ u2.T.conjugate() - dsquared, V = np.linalg.eig(u) + if predicates.is_hermitian(u): + # If `u` is hermitian, use the more accurate `eigh` method. + dsquared, V = np.linalg.eigh(u) + else: + dsquared, V = np.linalg.eig(u) + # Use Gram–Schmidt to obtain orthonormal eigenvectors for each of the subspaces. + for i in range(V.shape[0]): + for j in range(i): + if np.abs(dsquared[i] - dsquared[j]) < 1e-9: + V[:, i] -= np.dot(V[:, j].conj(), V[:, i]) * V[:, j] + V[:, i] /= np.linalg.norm(V[:, i]) # normalize. + dsquared = dsquared.astype(np.complex128) d = np.sqrt(dsquared) D = np.diag(d) W = D @ V.T.conjugate() @ u2 @@ -155,7 +231,7 @@ def _msb_demuxer( # Last term is given by ( I ⊗ W ), demultiplexed # Remove most-significant (demuxed) control-qubit # Yield operations for QSD on W - yield from quantum_shannon_decomposition(demux_qubits[1:], W) + yield from quantum_shannon_decomposition(demux_qubits[1:], W, atol=1e-6) # Use complex phase of d_i to give theta_i (so d_i* gives -theta_i) # Observe that middle part looks like Σ_i( Rz(theta_i)⊗|i> int: @@ -223,8 +299,9 @@ def _multiplexed_cossin( # so introduce max function select_qubit = max(-select_qubit - 1, -len(control_qubits)) - # Add a rotation on the main qubit - yield rot_func(rotation).on(main_qubit) + if np.abs(rotation) > 1e-9: + # Add a rotation on the main qubit + yield rot_func(rotation).on(main_qubit) # Add a CNOT from the select qubit to the main qubit yield ops.CNOT(control_qubits[select_qubit], main_qubit) diff --git a/cirq-core/cirq/transformers/analytical_decompositions/quantum_shannon_decomposition_test.py b/cirq-core/cirq/transformers/analytical_decompositions/quantum_shannon_decomposition_test.py index 3aab31ac478..740260852fb 100644 --- a/cirq-core/cirq/transformers/analytical_decompositions/quantum_shannon_decomposition_test.py +++ b/cirq-core/cirq/transformers/analytical_decompositions/quantum_shannon_decomposition_test.py @@ -34,9 +34,8 @@ def test_random_qsd_n_qubit(n_qubits): circuit = cirq.Circuit(quantum_shannon_decomposition(qubits, U)) # Test return is equal to inital unitary assert cirq.approx_eq(U, circuit.unitary(), atol=1e-9) - # Test all operations in gate set - gates = (common_gates.Rz, common_gates.Ry, common_gates.ZPowGate, common_gates.CXPowGate) - assert all(isinstance(op.gate, gates) for op in circuit.all_operations()) + # Test all operations have at most 2 qubits. + assert all(cirq.num_qubits(op) <= 2 for op in circuit.all_operations()) def test_qsd_n_qubit_errors(): @@ -53,9 +52,8 @@ def test_random_single_qubit_decomposition(): circuit = cirq.Circuit(_single_qubit_decomposition(qubit, U)) # Test return is equal to inital unitary assert cirq.approx_eq(U, circuit.unitary(), atol=1e-9) - # Test all operations in gate set - gates = (common_gates.Rz, common_gates.Ry, common_gates.ZPowGate, common_gates.CXPowGate) - assert all(isinstance(op.gate, gates) for op in circuit.all_operations()) + # Test all operations have at most 2 qubits. + assert all(cirq.num_qubits(op) <= 2 for op in circuit.all_operations()) def test_msb_demuxer(): @@ -66,9 +64,8 @@ def test_msb_demuxer(): circuit = cirq.Circuit(_msb_demuxer(qubits, U1, U2)) # Test return is equal to inital unitary assert cirq.approx_eq(U_full, circuit.unitary(), atol=1e-9) - # Test all operations in gate set - gates = (common_gates.Rz, common_gates.Ry, common_gates.ZPowGate, common_gates.CXPowGate) - assert all(isinstance(op.gate, gates) for op in circuit.all_operations()) + # Test all operations have at most 2 qubits. + assert all(cirq.num_qubits(op) <= 2 for op in circuit.all_operations()) def test_multiplexed_cossin(): @@ -110,3 +107,95 @@ def test_multiplexed_cossin(): ) def test_nth_gray(n, gray): assert _nth_gray(n) == gray + + +def test_ghz_circuit_decomposes(): + # Test case from #6725 + ghz_circuit = cirq.Circuit(cirq.H(cirq.q(0)), cirq.CNOT(cirq.q(0), cirq.q(1))) + ghz_unitary = cirq.unitary(ghz_circuit) + decomposed_circuit = cirq.Circuit( + quantum_shannon_decomposition(cirq.LineQubit.range(2), ghz_unitary) + ) + new_unitary = cirq.unitary(decomposed_circuit) + np.testing.assert_allclose(new_unitary, ghz_unitary, atol=1e-6) + + +def test_qft_decomposes(): + # Test case from #6666 + qs = cirq.LineQubit.range(4) + qft_circuit = cirq.Circuit(cirq.qft(*qs)) + qft_unitary = cirq.unitary(qft_circuit) + decomposed_circuit = cirq.Circuit(quantum_shannon_decomposition(qs, qft_unitary)) + new_unitary = cirq.unitary(decomposed_circuit) + np.testing.assert_allclose(new_unitary, qft_unitary, atol=1e-6) + + +# Cliffords test the different corner cases of the ZYZ decomposition. +@pytest.mark.parametrize( + ['gate', 'num_ops'], + [ + (cirq.I, 0), + (cirq.X, 2), # rz & ry + (cirq.Y, 1), # ry + (cirq.Z, 1), # rz + (cirq.H, 2), # rz & ry + (cirq.S, 1), # rz & ry + ], +) +def test_cliffords(gate, num_ops): + desired_unitary = cirq.unitary(gate) + shannon_circuit = cirq.Circuit(quantum_shannon_decomposition((cirq.q(0),), desired_unitary)) + new_unitary = cirq.unitary(shannon_circuit) + assert len([*shannon_circuit.all_operations()]) == num_ops + if num_ops: + np.testing.assert_allclose(new_unitary, desired_unitary) + + +@pytest.mark.parametrize('gate', [cirq.X, cirq.Y, cirq.Z, cirq.H, cirq.S]) +def test_cliffords_with_global_phase(gate): + global_phase = np.exp(1j * np.random.choice(np.linspace(0.1, 2 * np.pi, 10))) + desired_unitary = cirq.unitary(gate) * global_phase + shannon_circuit = cirq.Circuit(quantum_shannon_decomposition((cirq.q(0),), desired_unitary)) + new_unitary = cirq.unitary(shannon_circuit) + np.testing.assert_allclose(new_unitary, desired_unitary) + + +def test_global_phase(): + global_phase = np.exp(1j * np.random.choice(np.linspace(0, 2 * np.pi, 10))) + shannon_circuit = cirq.Circuit( + quantum_shannon_decomposition((cirq.q(0),), np.eye(2) * global_phase) + ) + new_unitary = cirq.unitary(shannon_circuit) + np.testing.assert_allclose(np.diag(new_unitary), global_phase) + + +@pytest.mark.parametrize('gate', [cirq.CZ, cirq.CNOT, cirq.XX, cirq.YY, cirq.ZZ]) +def test_two_qubit_gate(gate): + global_phase = np.exp(1j * np.random.choice(np.linspace(0, 2 * np.pi, 10))) + desired_unitary = cirq.unitary(gate) * global_phase + shannon_circuit = cirq.Circuit( + quantum_shannon_decomposition(cirq.LineQubit.range(2), desired_unitary) + ) + new_unitary = cirq.unitary(shannon_circuit) + np.testing.assert_allclose(new_unitary, desired_unitary, atol=1e-6) + + +@pytest.mark.parametrize('gate', [cirq.CCNOT, cirq.qft(*cirq.LineQubit.range(3))]) +def test_three_qubit_gate(gate): + global_phase = np.exp(1j * np.random.choice(np.linspace(0, 2 * np.pi, 10))) + desired_unitary = cirq.unitary(gate) * global_phase + shannon_circuit = cirq.Circuit( + quantum_shannon_decomposition(cirq.LineQubit.range(3), desired_unitary) + ) + new_unitary = cirq.unitary(shannon_circuit) + np.testing.assert_allclose(new_unitary, desired_unitary, atol=1e-6) + + +def test_qft5(): + global_phase = np.exp(1j * np.random.choice(np.linspace(0, 2 * np.pi, 10))) + desired_unitary = cirq.unitary(cirq.qft(*cirq.LineQubit.range(5))) * global_phase + shannon_circuit = cirq.Circuit( + quantum_shannon_decomposition(cirq.LineQubit.range(5), desired_unitary) + ) + new_unitary = cirq.unitary(shannon_circuit) + np.testing.assert_allclose(new_unitary, desired_unitary, atol=1e-6)