Adding Block ZXZ Decomposition to BQSKit

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',

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

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',

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()

bqskit/passes/processing/treescan.py

@@ -5,20 +5,20 @@
 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
 
 _logger = logging.getLogger(__name__)
 
 
-class TreeScanningGateRemovalPass(BasePass):
+class TreeScanningGateRemovalPass(ScanningGateRemovalPass):
     """
     The TreeScanningGateRemovalPass class.
 

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',
 ]