|
| 1 | +# Copyright 2023 The Cirq Developers |
| 2 | +# |
| 3 | +# Licensed under the Apache License, Version 2.0 (the "License"); |
| 4 | +# you may not use this file except in compliance with the License. |
| 5 | +# You may obtain a copy of the License at |
| 6 | +# |
| 7 | +# https://www.apache.org/licenses/LICENSE-2.0 |
| 8 | +# |
| 9 | +# Unless required by applicable law or agreed to in writing, software |
| 10 | +# distributed under the License is distributed on an "AS IS" BASIS, |
| 11 | +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| 12 | +# See the License for the specific language governing permissions and |
| 13 | +# limitations under the License. |
| 14 | + |
| 15 | + |
| 16 | +"""Utility methods for decomposing arbitrary n-qubit (2^n x 2^n) unitary. |
| 17 | +
|
| 18 | +Based on: |
| 19 | +Synthesis of Quantum Logic Circuits. Tech. rep. 2006, |
| 20 | +https://arxiv.org/abs/quant-ph/0406176 |
| 21 | +""" |
| 22 | +from typing import List, Callable, TYPE_CHECKING |
| 23 | + |
| 24 | +from scipy.linalg import cossin |
| 25 | + |
| 26 | +import numpy as np |
| 27 | + |
| 28 | +from cirq import ops |
| 29 | +from cirq.linalg import decompositions, predicates |
| 30 | + |
| 31 | +if TYPE_CHECKING: |
| 32 | + import cirq |
| 33 | + from cirq.ops import op_tree |
| 34 | + |
| 35 | + |
| 36 | +def quantum_shannon_decomposition(qubits: 'List[cirq.Qid]', u: np.ndarray) -> 'op_tree.OpTree': |
| 37 | + """Decomposes n-qubit unitary into CX/YPow/ZPow/CNOT gates, preserving global phase. |
| 38 | +
|
| 39 | + The algorithm is described in Shende et al.: |
| 40 | + Synthesis of Quantum Logic Circuits. Tech. rep. 2006, |
| 41 | + https://arxiv.org/abs/quant-ph/0406176 |
| 42 | +
|
| 43 | + Args: |
| 44 | + qubits: List of qubits in order of significance |
| 45 | + u: Numpy array for unitary matrix representing gate to be decomposed |
| 46 | +
|
| 47 | + Calls: |
| 48 | + (Base Case) |
| 49 | + 1. _single_qubit_decomposition |
| 50 | + OR |
| 51 | + (Recursive Case) |
| 52 | + 1. _msb_demuxer |
| 53 | + 2. _multiplexed_cossin |
| 54 | + 3. _msb_demuxer |
| 55 | +
|
| 56 | + Yields: |
| 57 | + A single 2-qubit or 1-qubit operations from OP TREE |
| 58 | + composed from the set |
| 59 | + { CNOT, rz, ry, ZPowGate } |
| 60 | +
|
| 61 | + Raises: |
| 62 | + ValueError: If the u matrix is non-unitary |
| 63 | + ValueError: If the u matrix is not of shape (2^n,2^n) |
| 64 | + """ |
| 65 | + if not predicates.is_unitary(u): # Check that u is unitary |
| 66 | + raise ValueError( |
| 67 | + "Expected input matrix u to be unitary, \ |
| 68 | + but it fails cirq.is_unitary check" |
| 69 | + ) |
| 70 | + |
| 71 | + n = u.shape[0] |
| 72 | + if n & (n - 1): |
| 73 | + raise ValueError( |
| 74 | + f"Expected input matrix u to be a (2^n x 2^n) shaped numpy array, \ |
| 75 | + but instead got shape {u.shape}" |
| 76 | + ) |
| 77 | + |
| 78 | + if n == 2: |
| 79 | + # Yield a single-qubit decomp if u is 2x2 |
| 80 | + yield from _single_qubit_decomposition(qubits[0], u) |
| 81 | + return |
| 82 | + |
| 83 | + # Perform a cosine-sine (linalg) decomposition on u |
| 84 | + # X = [ u1 , 0 ] [ cos(theta) , -sin(theta) ] [ v1 , 0 ] |
| 85 | + # [ 0 , u2 ] [ sin(theta) , cos(theta) ] [ 0 , v2 ] |
| 86 | + (u1, u2), theta, (v1, v2) = cossin(u, n / 2, n / 2, separate=True) |
| 87 | + |
| 88 | + # Yield ops from decomposition of multiplexed v1/v2 part |
| 89 | + yield from _msb_demuxer(qubits, v1, v2) |
| 90 | + |
| 91 | + # Observe that middle part looks like Σ_i( Ry(theta_i)⊗|i><i| ) |
| 92 | + # Then most significant qubit is Ry multiplexed over all other qubits |
| 93 | + # Yield ops from multiplexed Ry part |
| 94 | + yield from _multiplexed_cossin(qubits, theta, ops.ry) |
| 95 | + |
| 96 | + # Yield ops from decomposition of multiplexed u1/u2 part |
| 97 | + yield from _msb_demuxer(qubits, u1, u2) |
| 98 | + |
| 99 | + |
| 100 | +def _single_qubit_decomposition(qubit: 'cirq.Qid', u: np.ndarray) -> 'op_tree.OpTree': |
| 101 | + """Decomposes single-qubit gate, and returns list of operations, keeping phase invariant. |
| 102 | +
|
| 103 | + Args: |
| 104 | + qubit: Qubit on which to apply operations |
| 105 | + u: (2 x 2) Numpy array for unitary representing 1-qubit gate to be decomposed |
| 106 | +
|
| 107 | + Yields: |
| 108 | + A single operation from OP TREE of 3 operations (rz,ry,ZPowGate) |
| 109 | + """ |
| 110 | + # Perform native ZYZ decomposition |
| 111 | + phi_0, phi_1, phi_2 = decompositions.deconstruct_single_qubit_matrix_into_angles(u) |
| 112 | + |
| 113 | + # Determine global phase picked up |
| 114 | + phase = np.angle(u[0, 0] / (np.exp(-1j * (phi_0) / 2) * np.cos(phi_1 / 2))) |
| 115 | + |
| 116 | + # Append first two operations operations |
| 117 | + yield ops.rz(phi_0).on(qubit) |
| 118 | + yield ops.ry(phi_1).on(qubit) |
| 119 | + |
| 120 | + # Append third operation with global phase added |
| 121 | + yield ops.ZPowGate(exponent=phi_2 / np.pi, global_shift=phase / phi_2).on(qubit) |
| 122 | + |
| 123 | + |
| 124 | +def _msb_demuxer( |
| 125 | + demux_qubits: 'List[cirq.Qid]', u1: np.ndarray, u2: np.ndarray |
| 126 | +) -> 'op_tree.OpTree': |
| 127 | + """Demultiplexes a unitary matrix that is multiplexed in its most-significant-qubit. |
| 128 | +
|
| 129 | + Decomposition structure: |
| 130 | + [ u1 , 0 ] = [ V , 0 ][ D , 0 ][ W , 0 ] |
| 131 | + [ 0 , u2 ] [ 0 , V ][ 0 , D* ][ 0 , W ] |
| 132 | +
|
| 133 | + Gives: ( u1 )( u2* ) = ( V )( D^2 )( V* ) |
| 134 | + and: W = ( D )( V* )( u2 ) |
| 135 | +
|
| 136 | + Args: |
| 137 | + demux_qubits: Subset of total qubits involved in this unitary gate |
| 138 | + u1: Upper-left quadrant of total unitary to be decomposed (see diagram) |
| 139 | + u2: Lower-right quadrant of total unitary to be decomposed (see diagram) |
| 140 | +
|
| 141 | + Calls: |
| 142 | + 1. quantum_shannon_decomposition |
| 143 | + 2. _multiplexed_cossin |
| 144 | + 3. quantum_shannon_decomposition |
| 145 | +
|
| 146 | + Yields: Single operation from OP TREE of 2-qubit and 1-qubit operations |
| 147 | + """ |
| 148 | + # Perform a diagonalization to find values |
| 149 | + u = u1 @ u2.T.conjugate() |
| 150 | + dsquared, V = np.linalg.eig(u) |
| 151 | + d = np.sqrt(dsquared) |
| 152 | + D = np.diag(d) |
| 153 | + W = D @ V.T.conjugate() @ u2 |
| 154 | + |
| 155 | + # Last term is given by ( I ⊗ W ), demultiplexed |
| 156 | + # Remove most-significant (demuxed) control-qubit |
| 157 | + # Yield operations for QSD on W |
| 158 | + yield from quantum_shannon_decomposition(demux_qubits[1:], W) |
| 159 | + |
| 160 | + # Use complex phase of d_i to give theta_i (so d_i* gives -theta_i) |
| 161 | + # Observe that middle part looks like Σ_i( Rz(theta_i)⊗|i><i| ) |
| 162 | + # Yield ops from multiplexed Rz part |
| 163 | + yield from _multiplexed_cossin(demux_qubits, -np.angle(d), ops.rz) |
| 164 | + |
| 165 | + # Yield operations for QSD on V |
| 166 | + yield from quantum_shannon_decomposition(demux_qubits[1:], V) |
| 167 | + |
| 168 | + |
| 169 | +def _nth_gray(n: int) -> int: |
| 170 | + # Return the nth Gray Code number |
| 171 | + return n ^ (n >> 1) |
| 172 | + |
| 173 | + |
| 174 | +def _multiplexed_cossin( |
| 175 | + cossin_qubits: 'List[cirq.Qid]', angles: List[float], rot_func: Callable = ops.ry |
| 176 | +) -> 'op_tree.OpTree': |
| 177 | + """Performs a multiplexed rotation over all qubits in this unitary matrix, |
| 178 | +
|
| 179 | + Uses ry and rz multiplexing for quantum shannon decomposition |
| 180 | +
|
| 181 | + Args: |
| 182 | + cossin_qubits: Subset of total qubits involved in this unitary gate |
| 183 | + angles: List of angles to be multiplexed over for the given type of rotation |
| 184 | + rot_func: Rotation function used for this multiplexing implementation |
| 185 | + (cirq.ry or cirq.rz) |
| 186 | +
|
| 187 | + Calls: |
| 188 | + No major calls |
| 189 | +
|
| 190 | + Yields: Single operation from OP TREE from set 1- and 2-qubit gates: {ry,rz,CNOT} |
| 191 | + """ |
| 192 | + # Most significant qubit is main qubit with rotation function applied |
| 193 | + main_qubit = cossin_qubits[0] |
| 194 | + |
| 195 | + # All other qubits are control qubits |
| 196 | + control_qubits = cossin_qubits[1:] |
| 197 | + |
| 198 | + for j in range(len(angles)): |
| 199 | + # The rotation includes a factor of (-1) for each bit in the Gray Code |
| 200 | + # if the position of that bit is also 1 |
| 201 | + # The number of factors of -1 is counted using the 1s in the |
| 202 | + # binary representation of the (gray(j) & i) |
| 203 | + # Here, i gives the index for the angle, and |
| 204 | + # j is the iteration of the decomposition |
| 205 | + rotation = sum( |
| 206 | + -angle if bin(_nth_gray(j) & i).count('1') % 2 else angle |
| 207 | + for i, angle in enumerate(angles) |
| 208 | + ) |
| 209 | + |
| 210 | + # Divide by a factor of 2 for each additional select qubit |
| 211 | + # This is due to the halving in the decomposition applied recursively |
| 212 | + rotation = rotation * 2 / len(angles) |
| 213 | + |
| 214 | + # The XOR of the this gray code with the next will give the 1 for the bit |
| 215 | + # corresponding to the CNOT select, else 0 |
| 216 | + select_string = _nth_gray(j) ^ _nth_gray(j + 1) |
| 217 | + |
| 218 | + # Find the index number where the bit is 1 |
| 219 | + select_qubit = next(i for i in range(len(angles)) if (select_string >> i & 1)) |
| 220 | + |
| 221 | + # Negate the value, since we must index starting at most significant qubit |
| 222 | + # Also the final value will overflow, and it should be the MSB, |
| 223 | + # so introduce max function |
| 224 | + select_qubit = max(-select_qubit - 1, -len(control_qubits)) |
| 225 | + |
| 226 | + # Add a rotation on the main qubit |
| 227 | + yield rot_func(rotation).on(main_qubit) |
| 228 | + |
| 229 | + # Add a CNOT from the select qubit to the main qubit |
| 230 | + yield ops.CNOT(control_qubits[select_qubit], main_qubit) |
0 commit comments