diff --git a/bqskit/ir/gates/parameterized/__init__.py b/bqskit/ir/gates/parameterized/__init__.py index 34b51aee..44627a35 100644 --- a/bqskit/ir/gates/parameterized/__init__.py +++ b/bqskit/ir/gates/parameterized/__init__.py @@ -10,6 +10,7 @@ from bqskit.ir.gates.parameterized.cry import CRYGate from bqskit.ir.gates.parameterized.crz import CRZGate from bqskit.ir.gates.parameterized.cu import CUGate +from bqskit.ir.gates.parameterized.diagonal import DiagonalGate from bqskit.ir.gates.parameterized.fsim import FSIMGate from bqskit.ir.gates.parameterized.mpry import MPRYGate from bqskit.ir.gates.parameterized.mprz import MPRZGate @@ -42,6 +43,7 @@ 'CRYGate', 'CRZGate', 'CUGate', + 'DiagonalGate', 'FSIMGate', 'MPRYGate', 'MPRZGate', diff --git a/bqskit/ir/gates/parameterized/diagonal.py b/bqskit/ir/gates/parameterized/diagonal.py new file mode 100644 index 00000000..e83f3006 --- /dev/null +++ b/bqskit/ir/gates/parameterized/diagonal.py @@ -0,0 +1,84 @@ +"""This module implements a general Diagonal Gate.""" +from __future__ import annotations + +import numpy as np +import numpy.typing as npt + +from bqskit.ir.gates.qubitgate import QubitGate +from bqskit.qis.unitary.optimizable import LocallyOptimizableUnitary +from bqskit.qis.unitary.unitary import RealVector +from bqskit.qis.unitary.unitarymatrix import UnitaryMatrix +from bqskit.utils.cachedclass import CachedClass + + +class DiagonalGate( + QubitGate, + CachedClass, + LocallyOptimizableUnitary, +): + """ + A gate representing a general diagonal unitary. The top-left element is + fixed to 1, and the rest are set to exp(i * theta). + + This gate is used to optimize the Block ZXZ decomposition of a unitary. + """ + _qasm_name = 'diag' + + def __init__( + self, + num_qudits: int = 2, + ): + self._num_qudits = num_qudits + # 1 parameter per diagonal element, removing one for global phase + self._num_params = 2 ** num_qudits - 1 + + def get_unitary(self, params: RealVector = []) -> UnitaryMatrix: + """Return the unitary for this gate, see :class:`Unitary` for more.""" + self.check_parameters(params) + + mat = np.eye(2 ** self.num_qudits, dtype=np.complex128) + + for i in range(1, 2 ** self.num_qudits): + mat[i][i] = np.exp(1j * params[i - 1]) + + return UnitaryMatrix(mat) + + def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]: + """ + Return the gradient for this gate. + + See :class:`DifferentiableUnitary` for more info. + """ + self.check_parameters(params) + + grad = np.zeros( + ( + len(params), 2 ** self.num_qudits, + 2 ** self.num_qudits, + ), dtype=np.complex128, + ) + + for i, ind in enumerate(range(1, 2 ** self.num_qudits)): + grad[i][ind][ind] = 1j * np.exp(1j * params[i]) + + return grad + + def optimize(self, env_matrix: npt.NDArray[np.complex128]) -> list[float]: + """ + Return the optimal parameters with respect to an environment matrix. + + See :class:`LocallyOptimizableUnitary` for more info. + """ + self.check_env_matrix(env_matrix) + thetas = [0.0] * self.num_params + + base = env_matrix[0, 0] + if base == 0: + base = np.max(env_matrix[0, :]) + + for i in range(1, 2 ** self.num_qudits): + # Optimize each angle independently + a = np.angle(env_matrix[i, i] / base) + thetas[i - 1] = -1 * a + + return thetas diff --git a/bqskit/passes/__init__.py b/bqskit/passes/__init__.py index b05f5324..ee0015a8 100644 --- a/bqskit/passes/__init__.py +++ b/bqskit/passes/__init__.py @@ -291,11 +291,16 @@ from bqskit.passes.search.heuristics.astar import AStarHeuristic from bqskit.passes.search.heuristics.dijkstra import DijkstraHeuristic from bqskit.passes.search.heuristics.greedy import GreedyHeuristic +from bqskit.passes.synthesis.bzxz import BlockZXZPass +from bqskit.passes.synthesis.bzxz import FullBlockZXZPass from bqskit.passes.synthesis.diagonal import WalshDiagonalSynthesisPass from bqskit.passes.synthesis.leap import LEAPSynthesisPass from bqskit.passes.synthesis.pas import PermutationAwareSynthesisPass from bqskit.passes.synthesis.qfast import QFASTDecompositionPass from bqskit.passes.synthesis.qpredict import QPredictDecompositionPass +from bqskit.passes.synthesis.qsd import FullQSDPass +from bqskit.passes.synthesis.qsd import MGDPass +from bqskit.passes.synthesis.qsd import QSDPass from bqskit.passes.synthesis.qsearch import QSearchSynthesisPass from bqskit.passes.synthesis.synthesis import SynthesisPass from bqskit.passes.synthesis.target import SetTargetPass @@ -332,6 +337,11 @@ 'WalshDiagonalSynthesisPass', 'LEAPSynthesisPass', 'QSearchSynthesisPass', + 'FullQSDPass', + 'QSDPass', + 'MGDPass', + 'BlockZXZPass', + 'FullBlockZXZPass', 'QFASTDecompositionPass', 'QPredictDecompositionPass', 'CompressPass', diff --git a/bqskit/passes/processing/extract_diagonal.py b/bqskit/passes/processing/extract_diagonal.py new file mode 100644 index 00000000..898509b1 --- /dev/null +++ b/bqskit/passes/processing/extract_diagonal.py @@ -0,0 +1,188 @@ +"""This module implements the ExtractDiagonalPass.""" +from __future__ import annotations + +from typing import Any + +from bqskit.compiler.basepass import BasePass +from bqskit.compiler.passdata import PassData +from bqskit.ir.circuit import Circuit +from bqskit.ir.gates import DiagonalGate +from bqskit.ir.gates import VariableUnitaryGate +from bqskit.ir.gates.constant import CNOTGate +from bqskit.ir.operation import Operation +from bqskit.ir.opt.cost.functions import HilbertSchmidtResidualsGenerator +from bqskit.ir.opt.cost.generator import CostFunctionGenerator +from bqskit.qis.unitary.unitarymatrix import UnitaryMatrix + + +theorized_bounds = [0, 0, 3, 14, 61, 252] + + +def construct_linear_ansatz(num_qudits: int) -> Circuit: + """ + Generate a linear ansatz for extracting the diagonal of a unitary. + + This ansatz consists of a `num_qudits` width Diagonal Gate followed + by a ladder of CNOTs and single qubit gates. Right now, we try to use + one fewer CNOT than the theorized minimum number of CNOTs to represent + the unitary. + + This ansatz is simply a heuristic and does not have theoretical + backing. However, we see that for unitaries up to 5 qubits, this + ansatz does succeed most of the time with a threshold of 1e-8. + """ + theorized_num = theorized_bounds[num_qudits] + circuit = Circuit(num_qudits) + circuit.append_gate(DiagonalGate(num_qudits), tuple(range(num_qudits))) + for i in range(num_qudits): + circuit.append_gate(VariableUnitaryGate(1), (i,)) + for _ in range(theorized_num // (num_qudits - 1)): + # Apply n - 1 linear CNOTs + for i in range(num_qudits - 1): + circuit.append_gate(CNOTGate(), (i, i + 1)) + circuit.append_gate(VariableUnitaryGate(1), (i,)) + circuit.append_gate(VariableUnitaryGate(1), (i + 1,)) + return circuit + + +class ExtractDiagonalPass(BasePass): + """ + A pass that attempts to extract a diagonal matrix from a unitary matrix. + + https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=1269020 + + While there is a known algorithm for 2-qubit gates, we utilize + synthesis methods instead to scale to wider qubit gates. + + As a heuristic, we attempt to extrac the diagonal using a linear chain + ansatz of CNOT gates. We have found that up to 5 qubits, this ansatz + does succeed for most unitaries with fewer CNOTs than the theoretical + minimum number of CNOTs (utilizing the power of the Diagonal Gate in front). + """ + + def __init__( + self, + qudit_size: int = 2, + success_threshold: float = 1e-8, + cost: CostFunctionGenerator = HilbertSchmidtResidualsGenerator(), + instantiate_options: dict[str, Any] = {}, + ) -> None: + # We only support diagonal extraction for 2-5 qubits + assert qudit_size >= 2 and qudit_size <= 5 + self.qudit_size = qudit_size + self.success_threshold = success_threshold + self.cost = cost + self.instantiate_options: dict[str, Any] = { + 'cost_fn_gen': self.cost, + 'min_iters': 0, + 'diff_tol_r': 1e-4, + 'multistarts': 16, + 'method': 'qfactor', + } + self.instantiate_options.update(instantiate_options) + super().__init__() + + async def decompose( + self, + op: Operation, + target: UnitaryMatrix, + cost: CostFunctionGenerator = HilbertSchmidtResidualsGenerator(), + success_threshold: float = 1e-14, + instantiate_options: dict[str, Any] = {}, + ) -> tuple[ + Operation | None, + Circuit, + ]: + """ + Return the circuit that is generated from one levl of QSD. + + Args: + op (Operation): The VariableUnitaryGate Operation to decompose. + target (UnitaryMatrix): The target unitary. + cost (CostFunctionGenerator): The cost function generator to + determine if we have succeeded in decomposing the gate. + success_threshold (float): The threshold for the cost function. + instantiate_options (dict[str, Any]): The options to pass to the + instantiate method. + """ + + circ = Circuit(op.gate.num_qudits) + + if op.gate.num_qudits == 2: + # For now just try for 2 qubit + circ.append_gate(DiagonalGate(op.gate.num_qudits), (0, 1)) + circ.append_gate(VariableUnitaryGate(op.gate.num_qudits - 1), (0,)) + circ.append_gate(VariableUnitaryGate(op.gate.num_qudits - 1), (1,)) + circ.append_gate(CNOTGate(), (0, 1)) + circ.append_gate(VariableUnitaryGate(op.gate.num_qudits - 1), (0,)) + circ.append_gate(VariableUnitaryGate(op.gate.num_qudits - 1), (1,)) + circ.append_gate(CNOTGate(), (0, 1)) + circ.append_gate(VariableUnitaryGate(op.gate.num_qudits - 1), (0,)) + circ.append_gate(VariableUnitaryGate(op.gate.num_qudits - 1), (1,)) + elif op.gate.num_qudits == 3: + circ = construct_linear_ansatz(op.gate.num_qudits) + else: + circ = construct_linear_ansatz(op.gate.num_qudits) + + instantiated_circ = circ.instantiate( + target=target, + **instantiate_options, + ) + + if cost.calc_cost(instantiated_circ, target) < success_threshold: + diag_op = instantiated_circ.pop((0, 0)) + return diag_op, instantiated_circ + + default_circ = Circuit(op.gate.num_qudits) + default_circ.append_gate( + op.gate, + tuple(range(op.gate.num_qudits)), op.params, + ) + return None, default_circ + + async def run(self, circuit: Circuit, data: PassData) -> None: + """Synthesize `utry`, see :class:`SynthesisPass` for more.""" + num_ops = 0 + + num_gates_to_consider = circuit.count( + VariableUnitaryGate(self.qudit_size), + ) + + while num_gates_to_consider > 1: + # Find last Unitary + all_ops = list(circuit.operations_with_cycles(reverse=True)) + found = False + for cyc, op in all_ops: + if ( + isinstance(op.gate, VariableUnitaryGate) + and op.gate.num_qudits in [2, 3, 4] + ): + if found: + merge_op = op + merge_pt = (cyc, op.location[0]) + merge_location = op.location + break + else: + num_ops += 1 + gate = op + pt = (cyc, op.location[0]) + found = True + diag_op, circ = await self.decompose( + gate, + cost=self.cost, + target=gate.get_unitary(), + success_threshold=self.success_threshold, + instantiate_options=self.instantiate_options, + ) + + circuit.replace_with_circuit(pt, circ, as_circuit_gate=True) + num_gates_to_consider -= 1 + # Commute Diagonal into next op + if diag_op: + new_mat = diag_op.get_unitary() @ merge_op.get_unitary() + circuit.replace_gate( + merge_pt, merge_op.gate, merge_location, + VariableUnitaryGate.get_params(new_mat), + ) + + circuit.unfold_all() diff --git a/bqskit/passes/processing/treescan.py b/bqskit/passes/processing/treescan.py index 12376e59..9093a46c 100644 --- a/bqskit/passes/processing/treescan.py +++ b/bqskit/passes/processing/treescan.py @@ -5,12 +5,12 @@ from typing import Any from typing import Callable -from bqskit.compiler.basepass import BasePass from bqskit.compiler.passdata import PassData from bqskit.ir.circuit import Circuit from bqskit.ir.operation import Operation from bqskit.ir.opt.cost.functions import HilbertSchmidtResidualsGenerator from bqskit.ir.opt.cost.generator import CostFunctionGenerator +from bqskit.passes.processing.scan import ScanningGateRemovalPass from bqskit.runtime import get_runtime from bqskit.utils.typing import is_integer from bqskit.utils.typing import is_real_number @@ -18,7 +18,7 @@ _logger = logging.getLogger(__name__) -class TreeScanningGateRemovalPass(BasePass): +class TreeScanningGateRemovalPass(ScanningGateRemovalPass): """ The TreeScanningGateRemovalPass class. diff --git a/bqskit/passes/synthesis/__init__.py b/bqskit/passes/synthesis/__init__.py index eea09723..4246d49f 100644 --- a/bqskit/passes/synthesis/__init__.py +++ b/bqskit/passes/synthesis/__init__.py @@ -1,11 +1,16 @@ """This package implements synthesis passes and synthesis related classes.""" from __future__ import annotations +from bqskit.passes.synthesis.bzxz import BlockZXZPass +from bqskit.passes.synthesis.bzxz import FullBlockZXZPass from bqskit.passes.synthesis.diagonal import WalshDiagonalSynthesisPass from bqskit.passes.synthesis.leap import LEAPSynthesisPass from bqskit.passes.synthesis.pas import PermutationAwareSynthesisPass from bqskit.passes.synthesis.qfast import QFASTDecompositionPass from bqskit.passes.synthesis.qpredict import QPredictDecompositionPass +from bqskit.passes.synthesis.qsd import FullQSDPass +from bqskit.passes.synthesis.qsd import MGDPass +from bqskit.passes.synthesis.qsd import QSDPass from bqskit.passes.synthesis.qsearch import QSearchSynthesisPass from bqskit.passes.synthesis.synthesis import SynthesisPass from bqskit.passes.synthesis.target import SetTargetPass @@ -19,4 +24,9 @@ 'SetTargetPass', 'PermutationAwareSynthesisPass', 'WalshDiagonalSynthesisPass', + 'FullQSDPass', + 'MGDPass', + 'QSDPass', + 'BlockZXZPass', + 'FullBlockZXZPass', ] diff --git a/bqskit/passes/synthesis/bzxz.py b/bqskit/passes/synthesis/bzxz.py new file mode 100644 index 00000000..e5a33160 --- /dev/null +++ b/bqskit/passes/synthesis/bzxz.py @@ -0,0 +1,364 @@ +""" +This module implements the Full Block ZXZ Decomposition. + +Additionally, it defines the Block ZXZ Decomposition for a single pass. +""" +from __future__ import annotations + +import logging +from typing import Any + +import numpy as np +from scipy.linalg import eig +from scipy.linalg import svd + +from bqskit.compiler.basepass import BasePass +from bqskit.compiler.passdata import PassData +from bqskit.compiler.workflow import Workflow +from bqskit.ir.circuit import Circuit +from bqskit.ir.gates import CircuitGate +from bqskit.ir.gates.constant import HGate +from bqskit.ir.gates.constant import IdentityGate +from bqskit.ir.gates.constant import ZGate +from bqskit.ir.gates.parameterized.mprz import MPRZGate +from bqskit.ir.location import CircuitLocation +from bqskit.ir.operation import Operation +from bqskit.passes.processing.extract_diagonal import ExtractDiagonalPass +from bqskit.passes.processing.scan import ScanningGateRemovalPass +from bqskit.passes.processing.treescan import TreeScanningGateRemovalPass +from bqskit.passes.synthesis.qsd import MGDPass +from bqskit.passes.synthesis.qsd import QSDPass +from bqskit.qis.unitary.unitary import RealVector +from bqskit.qis.unitary.unitarymatrix import UnitaryMatrix +from bqskit.runtime import get_runtime + + +_logger = logging.getLogger(__name__) + + +class FullBlockZXZPass(BasePass): + """ + A pass performing a full Block ZXZ decomposition. + + References: + C.C. Paige, M. Wei, + History and generality of the CS decomposition, + Linear Algebra and its Applications, + Volumes 208–209, + 1994, + Pages 303-326, + ISSN 0024-3795, + https://doi.org/10.1016/0024-3795(94)90446-4. + """ + + def __init__( + self, + min_qudit_size: int = 2, + perform_scan: bool = False, + start_from_left: bool = True, + tree_depth: int = 0, + perform_extract: bool = True, + instantiate_options: dict[str, Any] = {}, + ) -> None: + """ + Perform the full Block ZXZ decomposition. + + Args: + min_qudit_size (int): Performs QSD until the circuit only has + VariableUnitaryGates with a number of qudits less than or equal + to this value. (Default: 2) + perform_scan (bool): Whether or not to perform the scanning + gate removal pass. (Default: False) + start_from_left (bool): Determines where the scan starts + attempting to remove gates from. If True, scan goes left + to right, otherwise right to left. (Default: True) + tree_depth (int): The depth of the tree to use in the + TreeScanningGateRemovalPass. If set to 0, we will instead + use the ScanningGateRemovalPass. (Default: 0) + perform_extract (bool): Whether or not to perform the diagonal + extraction pass. (Default: True) + instantiate_options (dict): The options to pass to the + scanning gate removal pass. (Default: {}) + """ + + self.start_from_left = start_from_left + self.min_qudit_size = min_qudit_size + instantiation_options = {'method': 'qfactor'} + instantiation_options.update(instantiate_options) + self.perform_scan = perform_scan + self.perform_extract = perform_extract + self.scan = ScanningGateRemovalPass( + start_from_left=start_from_left, + instantiate_options=instantiation_options, + ) + if tree_depth > 0: + self.scan = TreeScanningGateRemovalPass( + start_from_left=start_from_left, + instantiate_options=instantiation_options, + tree_depth=tree_depth, + ) + self.bzxz = BlockZXZPass(min_qudit_size=min_qudit_size) + self.mgd = MGDPass() + self.diag = ExtractDiagonalPass(qudit_size=min_qudit_size) + + async def run(self, circuit: Circuit, data: PassData) -> None: + """ + Perform succesive rounds of Block ZXZ decomposition until the circuit is + fully decomposed with no VariableUnitaryGates larger then + `min_qudit_size`. + + At the end, attempt to extrac the diagonal gate and commute through the + circuit to find optimal CNOT counts. + """ + passes: list[BasePass] = [] + start_num = max(x.num_qudits for x in circuit.operations()) + + for _ in range(self.min_qudit_size, start_num): + passes.append(self.bzxz) + if self.perform_scan: + passes.append(self.scan) + + if self.perform_extract: + passes.append(self.diag) + + # Once we have commuted the diagonal gate, we can break down the + # multiplexed gates + for _ in range(1, start_num): + passes.append(self.mgd) + if self.perform_scan: + passes.append(self.scan) + + await Workflow(passes).run(circuit, data) + + +class BlockZXZPass(BasePass): + """ + A pass performing one round of decomposition from the Block ZXZ algorithm. + + References: + Krol, Anna M., and Zaid Al-Ars. + "Highly Efficient Decomposition of n-Qubit Quantum Gates + Based on Block-ZXZ Decomposition." + arXiv preprint arXiv:2403.13692 (2024). + https://arxiv.org/html/2403.13692v1 + """ + + def __init__( + self, + min_qudit_size: int = 4, + ) -> None: + """ + Construct a single level of the QSDPass. + + Args: + min_qudit_size (int): Performs a decomposition on all gates + with width > min_qudit_size + """ + self.min_qudit_size = min_qudit_size + + @staticmethod + def initial_decompose(U: UnitaryMatrix) -> tuple[ + UnitaryMatrix, + UnitaryMatrix, + UnitaryMatrix, + UnitaryMatrix, + ]: + """ + This function decomposes the a given unitary U into 2 controlled + unitaries (B, and C), 1 multiplexed unitary (A) and 2 Hadamard gates. + + We return the A1 and A2 (sub-unitaries of the multipelexed unitary A), + as well as the sub unitary B and C. + + By sub-unitary, we mean the unitary that will be applied to the bottom + n - 1 qubits given the value of the top qubit 0. For a controlled + unitary, that sub-unitary is either the Identity (if qubit 0 is 0) or a + unitary on n-1 qubits if qubit 0 is 1. + + This follows equations 5-9 in the paper. + """ + # To solve for A1, A2, and C, we must divide U into 4 submatrices. + # The upper left block is X, the upper right block is Y, the lower left + # block is U_21, and the lower right block is U_22. + + X = U[0:len(U) // 2, 0:len(U) // 2] + Y = U[0:len(U) // 2, len(U) // 2:] + U_21 = U[len(U) // 2:, 0:len(U) // 2] + U_22 = U[len(U) // 2:, len(U) // 2:] + + # We let X = V_X @ s_X @ W_X where V_X and W_X are unitary and s_X is + # a diagonal matrix of singular values. We use the SVD to find these + # matrices. We do the same decomposition for Y. + V_X, s_X, W_XH = svd(X) + V_Y, s_Y, W_YH = svd(Y) + + # We can then define S_X and S_Y as V_X @ s_X @ V_X^H and + # V_Y @ s_Y @ V_Y^H + S_X = V_X @ np.diag(s_X) @ V_X.conj().T + S_Y = V_Y @ np.diag(s_Y) @ V_Y.conj().T + + # We define U_X and U_Y as V_X^H @ U @ W_X and V_Y^H @ U @ W_Y + U_X = V_X @ W_XH + U_Y = V_Y @ W_YH + + # Now, to perform the decomposition as defined in Section 4.1 + # We can set C dagger as i(U_Y^H @ U_X) + CH = 1j * U_Y.conj().T @ U_X + + # We then set A_1 as S_X @ U_X and A_2 as U_21 + U_22 @ (iU_Y^H @ U_X) + A_1 = UnitaryMatrix((S_X + 1j * S_Y) @ U_X) + A_2 = UnitaryMatrix(U_21 + U_22 @ (1j * U_Y.conj().T @ U_X)) + + # Finally, we can set B as 2(A_1^H @ U) - I + I = np.eye(len(U) // 2) + B = UnitaryMatrix((2 * (A_1.conj().T @ X)) - I) + + return A_1, A_2, B, UnitaryMatrix(CH.conj().T) + + @staticmethod + def demultiplex( + U_1: UnitaryMatrix, + U_2: UnitaryMatrix, + ) -> tuple[ + UnitaryMatrix, + RealVector, + UnitaryMatrix, + ]: + """ + Demultiplex U that is block_diag(U_1, U_2) to 2 Unitary Matrices V and W + and a diagonal matrix D. + + The Diagonal Matrix D corresponds to a Multiplexed Rz gate acting on the + most significant qubit. + + We return the Unitary Matrices V and W and the parameters for the + corresponding MPRZ gate. + """ + # U can be decomposed into U = (I otimes V )(D otimes D†)(I otimes W ) + + # We can find V,D^2 by performing an eigen decomposition of + # U_1 @ U_2† + d2, V = eig(U_1 @ U_2.conj().T) + d = np.sqrt(d2) + D = np.diag(d) + + # We can then multiply to solve for W + W = D @ V.conj().T @ U_2 + + # We can use d to find the parameters for the MPRZ gate. + # Note than because and Rz does exp(-i * theta / 2), we must + # multiply by -2. + d_params: list[float] = list(np.angle(d) * -2) + + return UnitaryMatrix(V), d_params, UnitaryMatrix(W) + + @staticmethod + def zxz(orig_u: UnitaryMatrix) -> Circuit: + """Return the circuit that is generated from one levl of Block ZXZ + decomposition.""" + + # First calculate the A, B, and C matrices for the initial decomp + A_1, A_2, B, C = BlockZXZPass.initial_decompose(orig_u) + + # Now decompose thee multiplexed A gate and the controlled C gate + I = IdentityGate(orig_u.num_qudits - 1).get_unitary() + VA, AZ_params, WA = BlockZXZPass.demultiplex(A_1, A_2) + VC, CZ_params, WC = BlockZXZPass.demultiplex(I, C) + + # Now calculate optimized B_tilde (Section 5.2) and decompose + # We merge in WA and VC into B and then merge in the additional + # CZ gates from the optimization + small_I = IdentityGate(orig_u.num_qudits - 2).get_unitary() + Z = ZGate(1).get_unitary() + B_tilde_1 = UnitaryMatrix(WA @ VC) + B_tilde_2 = UnitaryMatrix( + np.kron(Z, small_I) @ + WA @ B @ VC @ np.kron(Z, small_I), + ) + VB, BZ_params, WB = BlockZXZPass.demultiplex(B_tilde_1, B_tilde_2) + + # Define circuit locations + # We let the target qubit be qubit 0 (top qubit in diagram) + # The select qubits are the remaining qubits + controlled_qubit = 0 + select_qubits = list(range(1, orig_u.num_qudits)) + all_qubits = list(range(0, orig_u.num_qudits)) + + # Construct Circuit + circ = Circuit(orig_u.num_qudits) + + # Add WC gate + wc_gate, wc_params = QSDPass.create_unitary_gate(WC) + circ.append_gate(wc_gate, CircuitLocation(select_qubits), wc_params) + + # Add decomposed MPRZ gate circuit. + # Since the MPRZ gate circuit sets the target qubit as qubit + # num_qudits - 1, we must shift the qubits to the left + shifted_qubits = all_qubits[1:] + all_qubits[0:1] + circ.append_circuit( + MGDPass.decompose_mpx_two_levels( + decompose_ry=False, + params=CZ_params, + num_qudits=orig_u.num_qudits, + drop_last_cnot=True, + ), + shifted_qubits, + ) + + # Add the decomposed B-tilde gates WB and a Hadamard + combo_1_gate, combo_1_params = QSDPass.create_unitary_gate(WB) + circ.append_gate( + combo_1_gate, CircuitLocation(select_qubits), + combo_1_params, + ) + circ.append_gate(HGate(), CircuitLocation((controlled_qubit,))) + + # The central MPRZ_gate. We set the target to the controlled qubit, + # so there is no need to shift + z_gate = MPRZGate(len(all_qubits), 0) + circ.append_gate(z_gate, CircuitLocation(all_qubits), BZ_params) + + # Now add the decomposed B-tilde gates VB and a Hadamard + combo_2_gate, combo_2_params = QSDPass.create_unitary_gate(VB) + circ.append_gate( + combo_2_gate, CircuitLocation(select_qubits), + combo_2_params, + ) + circ.append_gate(HGate(), CircuitLocation((controlled_qubit,))) + + # Add the decomposed MPRZ gate circuit again on shifted qubits + circ.append_circuit( + MGDPass.decompose_mpx_two_levels( + decompose_ry=False, + params=AZ_params, + num_qudits=orig_u.num_qudits, + reverse=True, + drop_last_cnot=True, + ), + shifted_qubits, + ) + + va_gate, va_params = QSDPass.create_unitary_gate(VA) + circ.append_gate(va_gate, CircuitLocation(select_qubits), va_params) + + # assert np.allclose(orig_u, circ.get_unitary()) + return circ + + async def run(self, circuit: Circuit, data: PassData) -> None: + """Perform a single level of the Block ZXZ decomposition.""" + unitaries, pts, locations = QSDPass.get_variable_unitary_pts( + circuit, self.min_qudit_size, + ) + + if len(unitaries) > 0: + circs = await get_runtime().map(BlockZXZPass.zxz, unitaries) + # Do bulk replace + circ_gates = [CircuitGate(x) for x in circs] + circ_ops = [ + Operation(x, locations[i], x._circuit.params) + for i, x in enumerate(circ_gates) + ] + circuit.batch_replace(pts, circ_ops) + circuit.unfold_all() + + circuit.unfold_all() diff --git a/bqskit/passes/synthesis/qsd.py b/bqskit/passes/synthesis/qsd.py new file mode 100644 index 00000000..55aeba11 --- /dev/null +++ b/bqskit/passes/synthesis/qsd.py @@ -0,0 +1,552 @@ +"""This module implements the Quantum Shannon Decomposition.""" +from __future__ import annotations + +import logging +from typing import Any + +import numpy as np +from scipy.linalg import cossin +from scipy.linalg import schur + +from bqskit.compiler.basepass import BasePass +from bqskit.compiler.passdata import PassData +from bqskit.compiler.workflow import Workflow +from bqskit.ir.circuit import Circuit +from bqskit.ir.circuit import CircuitLocation +from bqskit.ir.circuit import CircuitPoint +from bqskit.ir.gates import CircuitGate +from bqskit.ir.gates.constant import CNOTGate +from bqskit.ir.gates.parameterized import RYGate +from bqskit.ir.gates.parameterized import RZGate +from bqskit.ir.gates.parameterized import VariableUnitaryGate +from bqskit.ir.gates.parameterized.mpry import MPRYGate +from bqskit.ir.gates.parameterized.mprz import MPRZGate +from bqskit.ir.operation import Operation +from bqskit.passes.processing import ScanningGateRemovalPass +from bqskit.passes.processing import TreeScanningGateRemovalPass +from bqskit.qis.permutation import PermutationMatrix +from bqskit.qis.unitary.unitary import RealVector +from bqskit.qis.unitary.unitarymatrix import UnitaryMatrix +from bqskit.runtime import get_runtime + +_logger = logging.getLogger(__name__) + + +class FullQSDPass(BasePass): + """ + A pass performing one round of decomposition from the QSD algorithm. + + Important: This pass runs on VariableUnitaryGates only. Make sure to convert + any gates you want to decompose to VariableUnitaryGates before running this + pass. + + Additionally, ScanningGateRemovalPass will operate on the context of the + entire circuit. If your circuit is large, it is best to set `perform_scan` + to False. + + References: + C.C. Paige, M. Wei, + History and generality of the CS decomposition, + Linear Algebra and its Applications, + Volumes 208–209, + 1994, + Pages 303-326, + ISSN 0024-3795, + https://doi.org/10.1016/0024-3795(94)90446-4. + """ + + def __init__( + self, + min_qudit_size: int = 2, + perform_scan: bool = False, + start_from_left: bool = True, + tree_depth: int = 0, + instantiate_options: dict[str, Any] = {}, + ) -> None: + """ + Construct a single level of the QSDPass. + + Args: + min_qudit_size (int): Performs QSD until the circuit only has + VariableUnitaryGates with a number of qudits less than or equal + to this value. (Default: 2) + perform_scan (bool): Whether or not to perform the scanning + gate removal pass. (Default: False) + start_from_left (bool): Determines where the scan starts + attempting to remove gates from. If True, scan goes left + to right, otherwise right to left. (Default: True) + tree_depth (int): The depth of the tree to use in the + TreeScanningGateRemovalPass. If set to 0, we will instead + use the ScanningGateRemovalPass. (Default: 0) + instantiate_options (dict): The options to pass to the + scanning gate removal pass. (Default: {}) + """ + self.start_from_left = start_from_left + self.min_qudit_size = min_qudit_size + instantiation_options = {'method': 'qfactor'} + instantiation_options.update(instantiate_options) + self.scan = ScanningGateRemovalPass( + start_from_left=start_from_left, + instantiate_options=instantiation_options, + ) + if tree_depth > 0: + self.scan = TreeScanningGateRemovalPass( + start_from_left=start_from_left, + instantiate_options=instantiation_options, + tree_depth=tree_depth, + ) + # Instantiate the helper QSD pass + self.qsd = QSDPass(min_qudit_size=min_qudit_size) + # Instantiate the helper Multiplex Gate Decomposition pass + self.mgd = MGDPass(decompose_twice=False) + self.perform_scan = perform_scan + + async def run(self, circuit: Circuit, data: PassData) -> None: + """Run a round of QSD, Multiplex Gate Decomposition, and Scanning Gate + Removal (optionally) until you reach the desired qudit size gates.""" + passes: list[BasePass] = [] + start_num = max(x.num_qudits for x in circuit.operations()) + for _ in range(self.min_qudit_size, start_num): + passes.append(self.qsd) + if self.perform_scan: + passes.append(self.scan) + passes.append(self.mgd) + await Workflow(passes).run(circuit, data) + + +class MGDPass(BasePass): + """ + A pass performing one round of decomposition of the MPRY and MPRZ gates in a + circuit. + + References: + C.C. Paige, M. Wei, + History and generality of the CS decomposition, + Linear Algebra and its Applications, + Volumes 208–209, + 1994, + Pages 303-326, + ISSN 0024-3795, + https://arxiv.org/pdf/quant-ph/0406176.pdf + """ + + def __init__(self, decompose_twice: bool = True) -> None: + """ + The MGDPass decomposes all MPRY and MPRZ gates in a circuit. + + Args: + decompose_twice (bool): Whether to decompose the MPRZ gate twice. + This will save 2 CNOT gates in the decomposition. If false, + the pass will only decompose one level. (Default: True) + """ + self.decompose_twice = decompose_twice + + @staticmethod + def decompose_mpx_one_level( + decompose_ry: bool, + params: RealVector, + num_qudits: int, + reverse: bool = False, + drop_last_cnot: bool = False, + ) -> Circuit: + """ + Decompose Multiplexed Gate one level. + + Args: + params (RealVector): The parameters for the original MPRZ gate + num_qudits (int): The number of qudits in the MPRZ gate + reverse (bool): Whether to reverse the order of the gates (you can + decompose the gates in either order to get the same result) + drop_last_cnot (bool): Whether to drop the last CNOT gate. This + should be set if you are doing a 2 level decomposition to save 2 + CNOT gates. + + Returns: + Circuit: The circuit that decomposes the MPRZ gate + """ + + new_gate: MPRZGate | MPRYGate | RZGate | RYGate = RZGate() + if decompose_ry: + new_gate = RYGate() + + if (num_qudits >= 3): + if decompose_ry: + new_gate = MPRYGate(num_qudits - 1, num_qudits - 2) + else: + # Remove 1 qubit, last qubit is controlled + new_gate = MPRZGate(num_qudits - 1, num_qudits - 2) + + left_params, right_params = MPRYGate.get_decomposition(params) + circ = Circuit(num_qudits) + new_gate_location = list(range(1, num_qudits)) + cx_location = (0, num_qudits - 1) + + ops = [ + Operation(new_gate, new_gate_location, left_params), + Operation(CNOTGate(), cx_location), + Operation(new_gate, new_gate_location, right_params), + Operation(CNOTGate(), cx_location), + ] + + if drop_last_cnot: + ops.pop() + + if reverse: + ops.reverse() + + for op in ops: + circ.append(op) + + return circ + + @staticmethod + def decompose_mpx_two_levels( + decompose_ry: bool, + params: RealVector, + num_qudits: int, + reverse: bool = False, + drop_last_cnot: bool = False, + ) -> Circuit: + """ + We decompose a multiplexed RZ gate 2 levels deep. This allows you to + remove 2 CNOTs as per Figure 2 in + https://arxiv.org/pdf/quant-ph/0406176.pdf. + + Furthermore, in the context of the Block ZXZ decomposition, you can + set `drop_last_cnot` to True. This CNOT gets merged into a central gate, + which saves another 2 CNOTs. This is shown in section 5.2 of + https://arxiv.org/pdf/2403.13692v1.pdf. + + Args: + decompose_ry (bool): Whether to decompose the MPRY gate + params (RealVector): The parameters for the original MPR gate + num_qudits (int): The number of qudits in the MPR gate + reverse (bool): Whether to reverse the order of the gates (you can + decompose the gates in either order to get the same result) + drop_last_cnot (bool): Whether to drop the last CNOT gate (only + should be set to True if you are doing section 5.2 optimization) + + Returns: + Circuit: The circuit that decomposes the MPR gate + """ + + if num_qudits <= 2: + # If you have less than 3 qubits, just decompose one level + return MGDPass.decompose_mpx_one_level( + decompose_ry, + params, + num_qudits, + reverse, + ) + + # Get params for first decomposition of the MPRZ gate + left_params, right_params = MPRYGate.get_decomposition(params) + + # Decompose the MPRZ gate into 2 MPRZ gates, dropping the last CNOT + # Also Reverse the circuit for the right side in order to do the + # optimization in section 5.2 + circ_left = MGDPass.decompose_mpx_one_level( + decompose_ry, + left_params, + num_qudits - 1, + reverse, + drop_last_cnot=True, + ) + + circ_right = MGDPass.decompose_mpx_one_level( + decompose_ry, + right_params, + num_qudits - 1, + not reverse, + drop_last_cnot=True, + ) + + # Now, construct the circuit. + # This will generate the original MPRZ gate with the target qubit + # set as qubit num_qudits - 1 + circ = Circuit(num_qudits) + cx_location_big = (0, num_qudits - 1) + + ops: list[Circuit | Operation] = [ + circ_left, + Operation(CNOTGate(), cx_location_big), + circ_right, + Operation(CNOTGate(), cx_location_big), + ] + + # We can remove the last CNOT gate as per section 5.2 + if drop_last_cnot: + ops.pop() + + if reverse: + ops.reverse() + + for op in ops: + if isinstance(op, Operation): + circ.append(op) + else: + circ.append_circuit(op, list(range(1, num_qudits))) + + return circ + + async def run(self, circuit: Circuit, data: PassData) -> None: + """Decompose all MPRY and MPRZ gates in the circuit one level.""" + ops: list[Operation] = [] + pts: list[CircuitPoint] = [] + locations: list[CircuitLocation] = [] + all_ops = list(circuit.operations_with_cycles(reverse=True)) + + # Gather all of the multiplexed operations + for cyc, op in all_ops: + if isinstance(op.gate, MPRYGate) or isinstance(op.gate, MPRZGate): + ops.append(op) + pts.append(CircuitPoint((cyc, op.location[0]))) + # Adjust location based on current target, move target to last + # qudit + loc = list(op.location) + loc = ( + loc[0:op.gate.target_qubit] + + loc[(op.gate.target_qubit + 1):] + + [loc[op.gate.target_qubit]] + ) + locations.append(CircuitLocation(loc)) + + if len(ops) > 0: + # Do a bulk QSDs -> circs + if self.decompose_twice: + circs = [ + MGDPass.decompose_mpx_two_levels( + isinstance(op.gate, MPRYGate), + op.params, + op.num_qudits, + ) for op in ops + ] + else: + circs = [ + MGDPass.decompose_mpx_one_level( + isinstance(op.gate, MPRYGate), + op.params, + op.num_qudits, + ) for op in ops + ] + circ_gates = [CircuitGate(x) for x in circs] + circ_ops = [ + Operation(x, locations[i], x._circuit.params) + for i, x in enumerate(circ_gates) + ] + circuit.batch_replace(pts, circ_ops) + circuit.unfold_all() + + circuit.unfold_all() + + +class QSDPass(BasePass): + """ + A pass performing one round of decomposition from the QSD algorithm. + + This decomposition takes each unitary of size n and decomposes it + into a circuit with 4 VariableUnitaryGates of size n - 1 and 3 multiplexed + rotation gates. + + Important: This pass runs on VariableUnitaryGates only. + + References: + https://arxiv.org/pdf/quant-ph/0406176 + """ + + def __init__( + self, + min_qudit_size: int = 4, + ) -> None: + """ + Construct a single level of the QSDPass. + + Args: + min_qudit_size (int): Performs a decomposition on all gates + with width > min_qudit_size + """ + self.min_qudit_size = min_qudit_size + + @staticmethod + def shift_down_unitary( + num_qudits: int, + end_qubits: int, + ) -> PermutationMatrix: + """Return the Permutation Matrix that shifts the qubits down by 1 + qubit.""" + top_qubits = num_qudits - end_qubits + now_bottom_qubits = list(reversed(range(top_qubits))) + now_top_qubits = list(range(num_qudits - end_qubits, num_qudits)) + final_qudits = now_top_qubits + now_bottom_qubits + return PermutationMatrix.from_qubit_location(num_qudits, final_qudits) + + @staticmethod + def shift_up_unitary( + num_qudits: int, + end_qubits: int, + ) -> PermutationMatrix: + """Return the Permutation Matrix that shifts the qubits down by 1 + qubit.""" + bottom_qubits = list(range(end_qubits)) + top_qubits = list(reversed(range(end_qubits, num_qudits))) + final_qudits = top_qubits + bottom_qubits + return PermutationMatrix.from_qubit_location(num_qudits, final_qudits) + + @staticmethod + def create_unitary_gate(u: UnitaryMatrix) -> tuple[ + VariableUnitaryGate, + RealVector, + ]: + """Create a VariableUnitaryGate from a UnitaryMatrix.""" + gate = VariableUnitaryGate(u.num_qudits) + params = np.concatenate((np.real(u).flatten(), np.imag(u).flatten())) + return gate, params + + @staticmethod + def create_multiplexed_circ( + us: list[UnitaryMatrix], + select_qubits: list[int], + ) -> Circuit: + """ + Takes a list of 2 unitaries of size n. Returns a circuit that decomposes + the unitaries into a circuit with 2 unitaries of size n-1 and a + multiplexed controlled gate. + + Args: + us (list[UnitaryMatrix]): The unitaries to decompose + select_qubits (list[int]): The qubits to use as select qubits + controlled_qubit (int): The qubit to use as the controlled qubit + + Returns: + Circuit: The circuit that decomposes the unitaries + + Using this paper: https://arxiv.org/pdf/quant-ph/0406176.pdf. Thm 12 + """ + u1 = us[0] + u2 = us[1] + assert (u1.num_qudits == u2.num_qudits) + all_qubits = list(range(len(select_qubits) + 1)) + # Use Schur Decomposition to split Us into V, D, and W matrices + D_2, V = schur(u1._utry @ u2.dagger._utry) + # D^2 will be diagonal since u1u2h is unitary + D = np.sqrt(np.diag(np.diag(D_2))) + # Calculate W @ U1 + left_mat = D @ V.conj().T @ u2._utry + left_gate, left_params = QSDPass.create_unitary_gate( + UnitaryMatrix(left_mat), + ) + + # Create Multi Controlled Z Gate + z_params: RealVector = np.array(-2 * np.angle(np.diag(D)).flatten()) + z_gate = MPRZGate(len(all_qubits), u1.num_qudits) + + # Create right gate + right_gate, right_params = QSDPass.create_unitary_gate(UnitaryMatrix(V)) + + circ = Circuit(u1.num_qudits + 1) + circ.append_gate(left_gate, CircuitLocation(select_qubits), left_params) + circ.append_gate(z_gate, CircuitLocation(all_qubits), z_params) + circ.append_gate( + right_gate, CircuitLocation( + select_qubits, + ), right_params, + ) + return circ + + @staticmethod + def mod_unitaries(u: UnitaryMatrix) -> UnitaryMatrix: + """Apply a permutation transform to the unitaries to the rest of the + circuit.""" + shift_up = QSDPass.shift_up_unitary(u.num_qudits, u.num_qudits - 1) + shift_down = QSDPass.shift_down_unitary(u.num_qudits, u.num_qudits - 1) + return shift_up @ u @ shift_down + + @staticmethod + def qsd(orig_u: UnitaryMatrix) -> Circuit: + ''' + Perform the Quantum Shannon Decomposition on a unitary matrix. + Args: + orig_u (UnitaryMatrix): The unitary matrix to decompose + + Returns: + Circuit: The circuit that decomposes the unitary + ''' + + # Shift the unitary qubits down by one + u = QSDPass.mod_unitaries(orig_u) + + # Perform CS Decomp to solve for multiplexed unitaries and theta_y + (u1, u2), theta_y, (v1h, v2h) = cossin( + u._utry, p=u.shape[0] / 2, q=u.shape[1] / 2, separate=True, + ) + assert (len(theta_y) == u.shape[0] / 2) + + # Create the multiplexed circuit + # This generates 2 circuits that multipex U,V with an MPRY gate + controlled_qubit = u.num_qudits - 1 + select_qubits = list(range(0, u.num_qudits - 1)) + all_qubits = list(range(u.num_qudits)) + circ_1 = QSDPass.create_multiplexed_circ( + [ + UnitaryMatrix(v1h), UnitaryMatrix(v2h), + ], + select_qubits, + ) + circ_2 = QSDPass.create_multiplexed_circ( + [ + UnitaryMatrix(u1), UnitaryMatrix(u2), + ], + select_qubits, + ) + gate_2 = MPRYGate(u.num_qudits, controlled_qubit) + + circ_1.append_gate(gate_2, CircuitLocation(all_qubits), 2 * theta_y) + circ_1.append_circuit( + circ_2, CircuitLocation(list(range(u.num_qudits))), + ) + return circ_1 + + @staticmethod + def get_variable_unitary_pts( + circuit: Circuit, + min_qudit_size: int, + ) -> tuple[list[UnitaryMatrix], list[CircuitPoint], list[CircuitLocation]]: + """Get all VariableUnitary Gates in the circuit wider than + `min_qudit_size` and return their unitaries, points, and locations.""" + unitaries: list[UnitaryMatrix] = [] + pts: list[CircuitPoint] = [] + locations: list[CircuitLocation] = [] + num_ops = 0 + all_ops = list(circuit.operations_with_cycles(reverse=True)) + + # Gather all of the VariableUnitary unitaries + for cyc, op in all_ops: + if ( + op.num_qudits > min_qudit_size + and isinstance(op.gate, VariableUnitaryGate) + ): + num_ops += 1 + unitaries.append(op.get_unitary()) + pts.append(CircuitPoint((cyc, op.location[0]))) + locations.append(op.location) + + return unitaries, pts, locations + + async def run(self, circuit: Circuit, data: PassData) -> None: + """Perform a single pass of Quantum Shannon Decomposition on the + circuit.""" + unitaries, pts, locations = QSDPass.get_variable_unitary_pts( + circuit, self.min_qudit_size, + ) + + if len(unitaries) > 0: + circs = await get_runtime().map(QSDPass.qsd, unitaries) + circ_gates = [CircuitGate(x) for x in circs] + circ_ops = [ + Operation(x, locations[i], x._circuit.params) + for i, x in enumerate(circ_gates) + ] + circuit.batch_replace(pts, circ_ops) + circuit.unfold_all() + + circuit.unfold_all() diff --git a/tests/passes/synthesis/test_bzxz.py b/tests/passes/synthesis/test_bzxz.py new file mode 100644 index 00000000..02993682 --- /dev/null +++ b/tests/passes/synthesis/test_bzxz.py @@ -0,0 +1,49 @@ +from __future__ import annotations + +import numpy as np + +from bqskit.compiler import Compiler +from bqskit.ir.circuit import Circuit +from bqskit.ir.gates.parameterized import VariableUnitaryGate +from bqskit.passes import BlockZXZPass +from bqskit.passes import FullBlockZXZPass +from bqskit.qis import UnitaryMatrix + + +def create_random_unitary_circ(num_qudits: int) -> Circuit: + """Create a Circuit with a random VariableUnitaryGate.""" + circuit = Circuit(num_qudits) + utry = UnitaryMatrix.random(num_qudits) + utry_params = np.concatenate(( + np.real(utry._utry).flatten(), + np.imag(utry._utry).flatten(), + )) + circuit.append_gate( + VariableUnitaryGate(num_qudits), + list(range(num_qudits)), + utry_params, + ) + return circuit + + +class TestBZXZ: + def test_small_qubit_bzxz(self, compiler: Compiler) -> None: + circuit = create_random_unitary_circ(4) + utry = circuit.get_unitary() + bzxz = BlockZXZPass(min_qudit_size=2) + circuit = compiler.compile(circuit, [bzxz]) + dist = circuit.get_unitary().get_distance_from(utry) + assert dist <= 1e-5 + + def test_full_bzxz_no_extract(self, compiler: Compiler) -> None: + circuit = create_random_unitary_circ(5) + utry = circuit.get_unitary() + bzxz = FullBlockZXZPass( + min_qudit_size=2, perform_scan=False, + perform_extract=False, + ) + circuit = compiler.compile(circuit, [bzxz]) + dist = circuit.get_unitary().get_distance_from(utry) + print(dist) + print(circuit.gate_counts) + assert dist <= 1e-5 diff --git a/tests/passes/synthesis/test_qsd.py b/tests/passes/synthesis/test_qsd.py new file mode 100644 index 00000000..6e500269 --- /dev/null +++ b/tests/passes/synthesis/test_qsd.py @@ -0,0 +1,62 @@ +from __future__ import annotations + +import numpy as np + +from bqskit.compiler import Compiler +from bqskit.ir.circuit import Circuit +from bqskit.ir.gates.parameterized import VariableUnitaryGate +from bqskit.passes import FullQSDPass +from bqskit.qis import UnitaryMatrix + + +def create_random_unitary_circ(num_qudits: int) -> Circuit: + """Create a Circuit with a random VariableUnitaryGate.""" + circuit = Circuit(num_qudits) + utry = UnitaryMatrix.random(num_qudits) + utry_params = np.concatenate(( + np.real(utry._utry).flatten(), + np.imag(utry._utry).flatten(), + )) + circuit.append_gate( + VariableUnitaryGate(num_qudits), + list(range(num_qudits)), + utry_params, + ) + return circuit + + +class TestQSD: + def test_three_qubit_qsd(self, compiler: Compiler) -> None: + circuit = create_random_unitary_circ(3) + utry = circuit.get_unitary() + # Run one pass of QSD + qsd = FullQSDPass(min_qudit_size=2, perform_scan=False) + circuit = compiler.compile(circuit, [qsd]) + dist = circuit.get_unitary().get_distance_from(utry) + assert circuit.count(VariableUnitaryGate(2)) == 4 + assert circuit.count(VariableUnitaryGate(3)) == 0 + assert dist <= 1e-5 + + def test_four_qubit_qubit(self, compiler: Compiler) -> None: + circuit = create_random_unitary_circ(4) + utry = circuit.get_unitary() + # Run two passes of QSD + qsd = FullQSDPass(min_qudit_size=2, perform_scan=False) + circuit = compiler.compile(circuit, [qsd]) + dist = circuit.get_unitary().get_distance_from(utry) + assert circuit.count(VariableUnitaryGate(2)) == 16 + assert circuit.count(VariableUnitaryGate(3)) == 0 + assert circuit.count(VariableUnitaryGate(4)) == 0 + assert dist <= 1e-5 + + def test_five_qubit_qsd(self, compiler: Compiler) -> None: + circuit = create_random_unitary_circ(5) + utry = circuit.get_unitary() + # Run two passes of QSD + qsd = FullQSDPass(min_qudit_size=3, perform_scan=False) + circuit = compiler.compile(circuit, [qsd]) + dist = circuit.get_unitary().get_distance_from(utry) + assert circuit.count(VariableUnitaryGate(3)) == 16 + assert circuit.count(VariableUnitaryGate(4)) == 0 + assert circuit.count(VariableUnitaryGate(5)) == 0 + assert dist <= 1e-5