"""This module regroups functions used for decomposition of arbitrary
unitary operator into elementary gates regrouped in a quantum circuit."""
from __future__ import annotations
import math
from typing import Union
import numpy as np
from mpqp.core.circuit import QCircuit
from mpqp.gates import CNOT, Ry, Rz
from mpqp.tools import Matrix
from mpqp.tools.maths import is_power_of_two
from scipy.linalg import cossin
PRECISION = 1e-9
def _gray_code(n: int):
return n ^ (n >> 1)
def _unitary_SVD(U: Matrix) -> tuple[Matrix, Matrix, Matrix]:
r"""
Returns the decomposition of the unitary using eigenvalues decomposition.
This decomposition is used to decompose matrices that are the result of a cosine-sine decomposition.
Hence we have $$U = \begin{pmatrix}U_0 & 0\\ 0 & U_1 \end{pmatrix}$$ with $U_0$ and $U_1$ being unitary matrices themselves.
The goal of this SVD is to have following result : $$U = (I \otimes V) \times D \times (I \otimes W)$$.
D is in the form of a multiplexed Rz operator.
The whole equation should be in the form :
`U = \begin{pmatrix}V & 0 \\ 0 & V end{pmatrix} \times \begin{pmatrix}D & 0 \\ 0 & D^{\dagger} \end{pmatrix} \times \begin{pmatrix}W & 0 \\ 0 & W \end{pmatrix}`
In the context of the quantum Shannon decomposition U is the result of a cosine sine decomposition, hence U is always in the form : `U = \begin{pmatrix}U_0 & 0\\ 0 & U_1 \end{pmatrix}`
Args:
U : The unitary matrix to decompose.
Returns:
A tuple of 3 matrices (V,D,W) so that U = $$(I \otimes V) \times D \times (I \otimes W)$$
Reference:
.. [1] V. Shende, S. S. Bullock, and I. Markov, “Synthesis of quantum logic
circuits,” Computer-Aided Design of Integrated Circuits and Systems,
IEEE Transactions on, vol. 25, pp. 1000 – 1010, July 2006.
"""
# Decompose U into it's 2 submatrices G0 and G1
length = len(U)
# U is a block diagonal matrix. Here we the two blocks in G0 and G1
SU = U[0 : (length // 2)]
SU2 = U[(length // 2) : length]
g0 = []
g1 = []
for i in range(length // 2):
g0.append(list(SU[i][0 : length // 2]))
g1.append(list(SU2[i][length // 2 : length]))
G0 = np.array(g0)
G1 = np.array(g1)
# Build G as G = V @ D² @ V†
G = G0 @ G1.conj().T
eigvals, V = np.linalg.eig(G)
D = np.diag(np.sqrt(eigvals.astype(complex)))
# W = D @ V† @ G1
W = np.asarray(D @ V.conj().T @ G1, dtype=np.complex128)
D_dagg = D.conj().T
# Reconstruct the whole D matrix
padding = np.zeros(length // 2)
D_result = []
for i in range(length // 2):
D_result.append(list(D[0 : length // 2][i]) + list(padding))
for i in range(length // 2):
D_result.append(list(padding) + list(D_dagg[0 : length // 2][i]))
return V, np.array(D_result), W
def _gray_code_decomposition(
thetas: Matrix,
circuit: QCircuit,
position: int,
rotation: Union[type[Rz], type[Ry]],
) -> QCircuit:
"""
Returns the decomposition of a multiplexed Rz or Ry gate.
The circuit is composed of a succession of CNOTs and rotations according to the original operator.
This function is originally from cirq's implementation.
Args:
thetas : A list of floats that are the rotations to be applied on each qubits.
circuit: The circuit in which the decomposition is stocked.
position : On which qubit is the rotation is taking place.
rotation : Either a Ry or a Rz rotation.
Returns:
The circuit containing the decomposed operator.
Reference:
.. [1] Mikko Möttönen, Juha J. Vartiainen, Ville Bergholm, and Martti M. Salomaa. 2004. Quantum circuits for general multi-qubit gates. American Physical Society (APS) : 93-13.
"""
for i in range(len(thetas)):
angle = sum(
-angle if bin(_gray_code(i) & j).count('1') % 2 else angle
for j, angle in enumerate(thetas)
)
angle *= 2 / len(thetas)
# CNOT's control is the changed bit of two consecutive numbers in gray code
changed = _gray_code(i) ^ _gray_code(i + 1)
control = next(i for i in range(len(thetas)) if (changed >> i & 1))
control = max(-control - position - 1 + circuit.nb_qubits, 1)
if np.abs(angle) > PRECISION: # Dodge unnecessary rotations
circuit.add(rotation(angle, position))
circuit.add(CNOT(control + position, position))
return circuit
def _decompose(U: Matrix, circuit: QCircuit, position: int = 0) -> QCircuit:
"""
This function recursively decompose the matrix U into the circuit then returns it.
For 1 qubit operators it executes a ZYZ decomposition.
For higher dimensions it does a Quantum Shannon decomposition.
Reference:
.. [1] Mikko Möttönen, Juha J. Vartiainen, Ville Bergholm, and Martti M. Salomaa. 2004. Quantum circuits for general multi-qubit gates. American Physical Society (APS) : 93-13.
"""
if len(U) == 2: # Decompose a 1 qubit operator
delta = np.angle(np.linalg.det(np.asarray(U, dtype=np.complex128))) / len(U)
# extract the global phase so that SU is a special unitary, note that if U is a special unitary then delta = 0
SU = U / np.exp(1j * delta)
beta = 2 * math.acos(np.abs(SU[0][0]))
alpha = -np.angle(SU[0][0]) - np.angle(SU[1][0])
gamma = -np.angle(SU[0][0]) + np.angle(SU[1][0])
circuit.add(Rz(alpha, position))
circuit.add(Ry(beta, position))
circuit.add(Rz(gamma, position))
circuit.input_g_phase += delta # Stores the gphase in the circuit
return circuit
else: # 2 qubits or more
length = len(U)
U12, MuxRy, V12 = cossin(U, p=length // 2, q=length // 2, separate=False)
# Extracts the rotations of the multiplexed Ry for later decomposition
thetas = []
for i in range(MuxRy.shape[0] // 2):
thetas.append(np.arccos(MuxRy[i][i]))
thetas = np.array(thetas)
assert isinstance(U12, np.ndarray)
assert isinstance(V12, np.ndarray)
Vu, MuxRzu, Wu = _unitary_SVD(U12)
Vv, MuxRzv, Wv = _unitary_SVD(V12)
# Extracts the rotations of both multiplexed Rz for later decomposition
du = np.angle(
MuxRzu.diagonal() # pyright: ignore[reportCallIssue, reportArgumentType]
)
for i in range(len(du) // 2):
du[i] *= -1
dv = np.angle(
MuxRzv.diagonal() # pyright: ignore[reportCallIssue, reportArgumentType]
)
for i in range(len(dv) // 2):
dv[i] *= -1
# Now recursively decompose every obtained matrices.
circuit = _decompose(Wv, circuit, position + 1)
circuit = _gray_code_decomposition(
dv, circuit, position, Rz # pyright: ignore[reportArgumentType]
)
circuit = _decompose(Vv, circuit, position + 1)
circuit = _gray_code_decomposition(
thetas,
circuit,
position,
Ry,
)
circuit = _decompose(Wu, circuit, position + 1)
circuit = _gray_code_decomposition(
du, circuit, position, Rz # pyright: ignore[reportArgumentType]
)
circuit = _decompose(Vu, circuit, position + 1)
return circuit
def _optimize_circuit(circuit: QCircuit) -> QCircuit:
"""
Optimize the circuit of the decomposition by removing useless CNOT gates.
"""
length = len(circuit.instructions)
i = 0
while i < length - 2:
if isinstance(circuit.instructions[i], CNOT):
j = i + 1
while isinstance(circuit.instructions[j], CNOT):
if circuit.instructions[i] == circuit.instructions[j]:
circuit.instructions.pop(j)
circuit.instructions.pop(i)
length -= 2
i -= 1
break
j += 1
i += 1
return circuit
[docs]
def quantum_shannon_decomposition(U: Matrix) -> QCircuit:
"""
Returns a circuit containing the decomposition of a unitary.
The resulting circuit is composed of gates CNOT, Ry and Rz.
The Quantum Shannon Decomposition works by splitting an unitary into 3 matrices using the cosine sine decomposition.
Which decompose an unitary matrix into 2 n-sized matrices around a multiplexed Ry gate.
The two other matrices are then decomposed into 2 n-1 sized unitaries surrounding a multiplexed Rz gate.
This process repeats until the matrices are reduced to 1 qubit gates which can be easily split with the ZYZ decomposition.
Args:
U: The unitary matrix to be decomposed
Returns:
A quantum circuit equivalent to U.
References:
.. [1] Mikko Möttönen, Juha J. Vartiainen, Ville Bergholm, and Martti M. Salomaa. 2004. Quantum circuits for general multi-qubit gates. American Physical Society (APS) : 93-13.
Examples:
>>> U = np.array([[1,0],[0,1]])
>>> circuit = quantum_shannon_decomposition(U)
>>> print(matrix_eq(U, circuit.to_matrix()))
True
"""
if not is_power_of_two(len(U)):
raise ValueError(
f"The size of the unitary matrix should be of the form 2**n : got {len(U)}"
)
circuit = QCircuit(int(np.log2(len(U))))
circuit = _decompose(U, circuit, 0)
return _optimize_circuit(circuit)