Fix Quantum Shannon Decomposition

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,30 @@ 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
+            )
+        )
+        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))
+        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 +150,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,24 +207,37 @@ 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
 
     # 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><i| )
     # Yield ops from multiplexed Rz part
     yield from _multiplexed_cossin(demux_qubits, -np.angle(d), ops.rz)
 
     # Yield operations for QSD on V
-    yield from quantum_shannon_decomposition(demux_qubits[1:], V)
+    yield from quantum_shannon_decomposition(demux_qubits[1:], V, atol=1e-6)
 
 
 def _nth_gray(n: int) -> int:
@@ -223,8 +297,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)

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('global_phase', np.exp(1j * np.linspace(0.1, 2 * np.pi, 10)))
+@pytest.mark.parametrize('gate', [cirq.X, cirq.Y, cirq.Z, cirq.H, cirq.S])
+def test_cliffords_with_global_phase(gate, global_phase):
+    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)
+
+
+@pytest.mark.parametrize('global_phase', np.exp(1j * np.linspace(0.1, 2 * np.pi, 10)))
+def test_global_phase(global_phase):
+    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])
+@pytest.mark.parametrize('global_phase', np.exp(1j * np.linspace(0, 2 * np.pi, 10)))
+def test_two_qubit_gate(gate, global_phase):
+    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))])
+@pytest.mark.parametrize('global_phase', np.exp(1j * np.linspace(0, 2 * np.pi, 10)))
+def test_three_qubit_gate(gate, global_phase):
+    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)
+
+
+@pytest.mark.parametrize('global_phase', np.exp(1j * np.linspace(0, 2 * np.pi, 4)))
+def test_qft5(global_phase):
+    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)