diff --git a/qiskit/synthesis/unitary/__init__.py b/qiskit/synthesis/unitary/__init__.py new file mode 100644 index 000000000000..f19592ee8bdb --- /dev/null +++ b/qiskit/synthesis/unitary/__init__.py @@ -0,0 +1,13 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2017 - 2023. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""Module containing unitary synthesis methods.""" diff --git a/qiskit/synthesis/unitary/aqc/__init__.py b/qiskit/synthesis/unitary/aqc/__init__.py new file mode 100644 index 000000000000..086782087f72 --- /dev/null +++ b/qiskit/synthesis/unitary/aqc/__init__.py @@ -0,0 +1,178 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2017 - 2023. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +r""" +===================================================================== +Approximate Quantum Compiler (:mod:`qiskit.synthesis.unitary.aqc`) +===================================================================== + +.. currentmodule:: qiskit.synthesis.unitary.aqc + +Implementation of Approximate Quantum Compiler as described in the paper [1]. + +Interface +========= + +The main public interface of this module is reached by passing ``unitary_synthesis_method='aqc'`` to +:obj:`~.compiler.transpile`. This will swap the synthesis method to use :obj:`AQCSynthesisPlugin`. +The individual classes are: + +.. autosummary:: + :toctree: ../stubs + :template: autosummary/class_no_inherited_members.rst + + AQC + AQCSynthesisPlugin + ApproximateCircuit + ApproximatingObjective + CNOTUnitCircuit + CNOTUnitObjective + DefaultCNOTUnitObjective + FastCNOTUnitObjective + + +Mathematical Detail +=================== + +We are interested in compiling a quantum circuit, which we formalize as finding the best +circuit representation in terms of an ordered gate sequence of a target unitary matrix +:math:`U\in U(d)`, with some additional hardware constraints. In particular, we look at +representations that could be constrained in terms of hardware connectivity, as well +as gate depth, and we choose a gate basis in terms of CNOT and rotation gates. +We recall that the combination of CNOT and rotation gates is universal in :math:`SU(d)` and +therefore it does not limit compilation. + +To properly define what we mean by best circuit representation, we define the metric +as the Frobenius norm between the unitary matrix of the compiled circuit :math:`V` and +the target unitary matrix :math:`U`, i.e., :math:`\|V - U\|_{\mathrm{F}}`. This choice +is motivated by mathematical programming considerations, and it is related to other +formulations that appear in the literature. Let's take a look at the problem in more details. + +Let :math:`n` be the number of qubits and :math:`d=2^n`. Given a CNOT structure :math:`ct` +and a vector of rotation angles :math:`\theta`, the parametric circuit forms a matrix +:math:`Vct(\theta)\in SU(d)`. If we are given a target circuit forming a matrix +:math:`U\in SU(d)`, then we would like to compute + +.. math:: + + argmax_{\theta}\frac{1}{d}|\langle Vct(\theta),U\rangle| + +where the inner product is the Frobenius inner product. Note that +:math:`|\langle V,U\rangle|\leq d` for all unitaries :math:`U` and :math:`V`, so the objective +has range in :math:`[0,1]`. + +Our strategy is to maximize + +.. math:: + + \frac{1}{d}\Re \langle Vct(\theta),U\rangle + +using its gradient. We will now discuss the specifics by going through an example. + +While the range of :math:`Vct` is a subset of :math:`SU(d)` by construction, the target +circuit may form a general unitary matrix. However, for any :math:`U\in U(d)`, + +.. math:: + + \frac{\exp(2\pi i k/d)}{\det(U)^{1/d}}U\in SU(d)\text{ for all }k\in\{0,\ldots,d-1\}. + +Thus, we should normalize the target circuit by its global phase and then approximately +compile the normalized circuit. We can add the global phase back in afterwards. + +In the algorithm let :math:`U'` denote the un-normalized target matrix and :math:`U` +the normalized target matrix. Now that we have :math:`U`, we give the gradient function +to the Nesterov's method optimizer and compute :math:`\theta`. + +To add the global phase back in, we can form the control circuit as + +.. math:: + + \frac{\langle Vct(\theta),U'\rangle}{|\langle Vct(\theta),U'\rangle|}Vct(\theta). + +Note that while we optimized using Nesterov's method in the paper, this was for its convergence +guarantees, not its speed in practice. It is much faster to use L-BFGS which is used as a +default optimizer in this implementation. + +A basic usage of the AQC algorithm should consist of the following steps:: + + # Define a target circuit as a unitary matrix + unitary = ... + + # Define a number of qubits for the algorithm, at least 3 qubits + num_qubits = int(round(np.log2(unitary.shape[0]))) + + # Choose a layout of the CNOT structure for the approximate circuit, e.g. ``spin`` for + # a linear layout. + layout = options.get("layout") or "spin" + + # Choose a connectivity type, e.g. ``full`` for full connectivity between qubits. + connectivity = options.get("connectivity") or "full" + + # Define a targeted depth of the approximate circuit in the number of CNOT units. + depth = int(options.get("depth") or 0) + + # Generate a network made of CNOT units + cnots = make_cnot_network( + num_qubits=num_qubits, + network_layout=layout, + connectivity_type=connectivity, + depth=depth + ) + + # Create an optimizer to be used by AQC + optimizer = L_BFGS_B() + + # Create an instance + aqc = AQC(optimizer) + + # Create a template circuit that will approximate our target circuit + approximate_circuit = CNOTUnitCircuit(num_qubits=num_qubits, cnots=cnots) + + # Create an objective that defines our optimization problem + approximating_objective = DefaultCNOTUnitObjective(num_qubits=num_qubits, cnots=cnots) + + # Run optimization process to compile the unitary + aqc.compile_unitary( + target_matrix=unitary, + approximate_circuit=approximate_circuit, + approximating_objective=approximating_objective + ) + +Now ``approximate_circuit`` is a circuit that approximates the target unitary to a certain +degree and can be used instead of the original matrix. + +This uses a helper function, :obj:`make_cnot_network`. + +.. autofunction:: make_cnot_network + +One can take advantage of accelerated version of objective function. It implements the same +mathematical algorithm as the default one ``DefaultCNOTUnitObjective`` but runs several times +faster. Instantiation of accelerated objective function class is similar to the default case: + + # Create an objective that defines our optimization problem + approximating_objective = FastCNOTUnitObjective(num_qubits=num_qubits, cnots=cnots) + +The rest of the code in the above example does not change. + +References: + + [1]: Liam Madden, Andrea Simonetto, Best Approximate Quantum Compiling Problems. + `arXiv:2106.05649 `_ +""" + +from .approximate import ApproximateCircuit, ApproximatingObjective +from .aqc import AQC +from .aqc_plugin import AQCSynthesisPlugin +from .cnot_structures import make_cnot_network +from .cnot_unit_circuit import CNOTUnitCircuit +from .cnot_unit_objective import CNOTUnitObjective, DefaultCNOTUnitObjective +from .fast_gradient.fast_gradient import FastCNOTUnitObjective diff --git a/qiskit/synthesis/unitary/aqc/approximate.py b/qiskit/synthesis/unitary/aqc/approximate.py new file mode 100644 index 000000000000..1b647d36c8b5 --- /dev/null +++ b/qiskit/synthesis/unitary/aqc/approximate.py @@ -0,0 +1,116 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2023. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +"""Base classes for an approximate circuit definition.""" +from __future__ import annotations +from abc import ABC, abstractmethod +from typing import Optional, SupportsFloat +import numpy as np + +from qiskit import QuantumCircuit + + +class ApproximateCircuit(QuantumCircuit, ABC): + """A base class that represents an approximate circuit.""" + + def __init__(self, num_qubits: int, name: Optional[str] = None) -> None: + """ + Args: + num_qubits: number of qubit this circuit will span. + name: a name of the circuit. + """ + super().__init__(num_qubits, name=name) + + @property + @abstractmethod + def thetas(self) -> np.ndarray: + """ + The property is not implemented and raises a ``NotImplementedException`` exception. + + Returns: + a vector of parameters of this circuit. + """ + raise NotImplementedError + + @abstractmethod + def build(self, thetas: np.ndarray) -> None: + """ + Constructs this circuit out of the parameters(thetas). Parameter values must be set before + constructing the circuit. + + Args: + thetas: a vector of parameters to be set in this circuit. + """ + raise NotImplementedError + + +class ApproximatingObjective(ABC): + """ + A base class for an optimization problem definition. An implementing class must provide at least + an implementation of the ``objective`` method. In such case only gradient free optimizers can + be used. Both method, ``objective`` and ``gradient``, preferable to have in an implementation. + """ + + def __init__(self) -> None: + # must be set before optimization + self._target_matrix: np.ndarray | None = None + + @abstractmethod + def objective(self, param_values: np.ndarray) -> SupportsFloat: + """ + Computes a value of the objective function given a vector of parameter values. + + Args: + param_values: a vector of parameter values for the optimization problem. + + Returns: + a float value of the objective function. + """ + raise NotImplementedError + + @abstractmethod + def gradient(self, param_values: np.ndarray) -> np.ndarray: + """ + Computes a gradient with respect to parameters given a vector of parameter values. + + Args: + param_values: a vector of parameter values for the optimization problem. + + Returns: + an array of gradient values. + """ + raise NotImplementedError + + @property + def target_matrix(self) -> np.ndarray: + """ + Returns: + a matrix being approximated + """ + return self._target_matrix + + @target_matrix.setter + def target_matrix(self, target_matrix: np.ndarray) -> None: + """ + Args: + target_matrix: a matrix to approximate in the optimization procedure. + """ + self._target_matrix = target_matrix + + @property + @abstractmethod + def num_thetas(self) -> int: + """ + + Returns: + the number of parameters in this optimization problem. + """ + raise NotImplementedError diff --git a/qiskit/synthesis/unitary/aqc/aqc.py b/qiskit/synthesis/unitary/aqc/aqc.py new file mode 100644 index 000000000000..1b5f928eb71d --- /dev/null +++ b/qiskit/synthesis/unitary/aqc/aqc.py @@ -0,0 +1,117 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2023. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +"""A generic implementation of Approximate Quantum Compiler.""" +from typing import Optional + +import numpy as np + +from qiskit.algorithms.optimizers import L_BFGS_B, Optimizer +from qiskit.quantum_info import Operator +from .approximate import ApproximateCircuit, ApproximatingObjective + + +class AQC: + """ + A generic implementation of Approximate Quantum Compiler. This implementation is agnostic of + the underlying implementation of the approximate circuit, objective, and optimizer. Users may + pass corresponding implementations of the abstract classes: + + * Optimizer is an instance of :class:`~qiskit.algorithms.optimizers.Optimizer` and used to run + the optimization process. A choice of optimizer may affect overall convergence, required time + for the optimization process and achieved objective value. + + * Approximate circuit represents a template which parameters we want to optimize. Currently, + there's only one implementation based on 4-rotations CNOT unit blocks: + :class:`.CNOTUnitCircuit`. See the paper for more details. + + * Approximate objective is tightly coupled with the approximate circuit implementation and + provides two methods for computing objective function and gradient with respect to approximate + circuit parameters. This objective is passed to the optimizer. Currently, there are two + implementations based on 4-rotations CNOT unit blocks: :class:`.DefaultCNOTUnitObjective` and + its accelerated version :class:`.FastCNOTUnitObjective`. Both implementations share the same + idea of maximization the Hilbert-Schmidt product between the target matrix and its + approximation. The former implementation approach should be considered as a baseline one. It + may suffer from performance issues, and is mostly suitable for a small number of qubits + (up to 5 or 6), whereas the latter, accelerated one, can be applied to larger problems. + + * One should take into consideration the exponential growth of matrix size with the number of + qubits because the implementation not only creates a potentially large target matrix, but + also allocates a number of temporary memory buffers comparable in size to the target matrix. + """ + + def __init__( + self, + optimizer: Optional[Optimizer] = None, + seed: Optional[int] = None, + ): + """ + Args: + optimizer: an optimizer to be used in the optimization procedure of the search for + the best approximate circuit. By default, :obj:`.L_BFGS_B` is used with max + iterations set to 1000. + seed: a seed value to be user by a random number generator. + """ + super().__init__() + self._optimizer = optimizer + self._seed = seed + + def compile_unitary( + self, + target_matrix: np.ndarray, + approximate_circuit: ApproximateCircuit, + approximating_objective: ApproximatingObjective, + initial_point: Optional[np.ndarray] = None, + ) -> None: + """ + Approximately compiles a circuit represented as a unitary matrix by solving an optimization + problem defined by ``approximating_objective`` and using ``approximate_circuit`` as a + template for the approximate circuit. + + Args: + target_matrix: a unitary matrix to approximate. + approximate_circuit: a template circuit that will be filled with the parameter values + obtained in the optimization procedure. + approximating_objective: a definition of the optimization problem. + initial_point: initial values of angles/parameters to start optimization from. + """ + matrix_dim = target_matrix.shape[0] + # check if it is actually a special unitary matrix + target_det = np.linalg.det(target_matrix) + if not np.isclose(target_det, 1): + su_matrix = target_matrix / np.power(target_det, (1 / matrix_dim), dtype=complex) + global_phase_required = True + else: + su_matrix = target_matrix + global_phase_required = False + + # set the matrix to approximate in the algorithm + approximating_objective.target_matrix = su_matrix + + optimizer = self._optimizer or L_BFGS_B(maxiter=1000) + + if initial_point is None: + np.random.seed(self._seed) + initial_point = np.random.uniform(0, 2 * np.pi, approximating_objective.num_thetas) + + opt_result = optimizer.minimize( + fun=approximating_objective.objective, + x0=initial_point, + jac=approximating_objective.gradient, + ) + + approximate_circuit.build(opt_result.x) + + approx_matrix = Operator(approximate_circuit).data + + if global_phase_required: + alpha = np.angle(np.trace(np.dot(approx_matrix.conj().T, target_matrix))) + approximate_circuit.global_phase = alpha diff --git a/qiskit/synthesis/unitary/aqc/aqc_plugin.py b/qiskit/synthesis/unitary/aqc/aqc_plugin.py new file mode 100644 index 000000000000..990f9d35a1fd --- /dev/null +++ b/qiskit/synthesis/unitary/aqc/aqc_plugin.py @@ -0,0 +1,144 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2023. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +""" +An AQC synthesis plugin to Qiskit's transpiler. +""" +import numpy as np + +from qiskit.converters import circuit_to_dag +from qiskit.transpiler.passes.synthesis.plugin import UnitarySynthesisPlugin + + +class AQCSynthesisPlugin(UnitarySynthesisPlugin): + """ + An AQC-based Qiskit unitary synthesis plugin. + + This plugin is invoked by :func:`~.compiler.transpile` when the ``unitary_synthesis_method`` + parameter is set to ``"aqc"``. + + This plugin supports customization and additional parameters can be passed to the plugin + by passing a dictionary as the ``unitary_synthesis_plugin_config`` parameter of + the :func:`~qiskit.compiler.transpile` function. + + Supported parameters in the dictionary: + + network_layout (str) + Type of network geometry, one of {``"sequ"``, ``"spin"``, ``"cart"``, ``"cyclic_spin"``, + ``"cyclic_line"``}. Default value is ``"spin"``. + + connectivity_type (str) + type of inter-qubit connectivity, {``"full"``, ``"line"``, ``"star"``}. Default value + is ``"full"``. + + depth (int) + depth of the CNOT-network, i.e. the number of layers, where each layer consists of a + single CNOT-block. + + optimizer (:class:`~qiskit.algorithms.optimizers.Optimizer`) + An instance of optimizer to be used in the optimization process. + + seed (int) + A random seed. + + initial_point (:class:`~numpy.ndarray`) + Initial values of angles/parameters to start the optimization process from. + """ + + @property + def max_qubits(self): + """Maximum number of supported qubits is ``14``.""" + return 14 + + @property + def min_qubits(self): + """Minimum number of supported qubits is ``3``.""" + return 3 + + @property + def supports_natural_direction(self): + """The plugin does not support natural direction, + it assumes bidirectional two qubit gates.""" + return False + + @property + def supports_pulse_optimize(self): + """The plugin does not support optimization of pulses.""" + return False + + @property + def supports_gate_lengths(self): + """The plugin does not support gate lengths.""" + return False + + @property + def supports_gate_errors(self): + """The plugin does not support gate errors.""" + return False + + @property + def supported_bases(self): + """The plugin does not support bases for synthesis.""" + return None + + @property + def supports_basis_gates(self): + """The plugin does not support basis gates and by default it synthesizes a circuit using + ``["rx", "ry", "rz", "cx"]`` gate basis.""" + return False + + @property + def supports_coupling_map(self): + """The plugin does not support coupling maps.""" + return False + + def run(self, unitary, **options): + + # Runtime imports to avoid the overhead of these imports for + # plugin discovery and only use them if the plugin is run/used + from qiskit.algorithms.optimizers import L_BFGS_B + from qiskit.synthesis.unitary.aqc.aqc import AQC + from qiskit.synthesis.unitary.aqc.cnot_structures import make_cnot_network + from qiskit.synthesis.unitary.aqc.cnot_unit_circuit import CNOTUnitCircuit + from qiskit.synthesis.unitary.aqc.cnot_unit_objective import DefaultCNOTUnitObjective + + num_qubits = int(round(np.log2(unitary.shape[0]))) + + config = options.get("config") or {} + + network_layout = config.get("network_layout", "spin") + connectivity_type = config.get("connectivity_type", "full") + depth = config.get("depth", 0) + + cnots = make_cnot_network( + num_qubits=num_qubits, + network_layout=network_layout, + connectivity_type=connectivity_type, + depth=depth, + ) + + optimizer = config.get("optimizer", L_BFGS_B(maxiter=1000)) + seed = config.get("seed") + aqc = AQC(optimizer, seed) + + approximate_circuit = CNOTUnitCircuit(num_qubits=num_qubits, cnots=cnots) + approximating_objective = DefaultCNOTUnitObjective(num_qubits=num_qubits, cnots=cnots) + + initial_point = config.get("initial_point") + aqc.compile_unitary( + target_matrix=unitary, + approximate_circuit=approximate_circuit, + approximating_objective=approximating_objective, + initial_point=initial_point, + ) + + dag_circuit = circuit_to_dag(approximate_circuit) + return dag_circuit diff --git a/qiskit/synthesis/unitary/aqc/cnot_structures.py b/qiskit/synthesis/unitary/aqc/cnot_structures.py new file mode 100644 index 000000000000..91873c0f01d3 --- /dev/null +++ b/qiskit/synthesis/unitary/aqc/cnot_structures.py @@ -0,0 +1,299 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2023. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +""" +These are the CNOT structure methods: anything that you need for creating CNOT structures. +""" +import logging + +import numpy as np + +_NETWORK_LAYOUTS = ["sequ", "spin", "cart", "cyclic_spin", "cyclic_line"] +_CONNECTIVITY_TYPES = ["full", "line", "star"] + + +logger = logging.getLogger(__name__) + + +def _lower_limit(num_qubits: int) -> int: + """ + Returns lower limit on the number of CNOT units that guarantees exact representation of + a unitary operator by quantum gates. + + Args: + num_qubits: number of qubits. + + Returns: + lower limit on the number of CNOT units. + """ + num_cnots = round(np.ceil((4**num_qubits - 3 * num_qubits - 1) / 4.0)) + return num_cnots + + +def make_cnot_network( + num_qubits: int, + network_layout: str = "spin", + connectivity_type: str = "full", + depth: int = 0, +) -> np.ndarray: + """ + Generates a network consisting of building blocks each containing a CNOT gate and possibly some + single-qubit ones. This network models a quantum operator in question. Note, each building + block has 2 input and outputs corresponding to a pair of qubits. What we actually return here + is a chain of indices of qubit pairs shared by every building block in a row. + + Args: + num_qubits: number of qubits. + network_layout: type of network geometry, ``{"sequ", "spin", "cart", "cyclic_spin", + "cyclic_line"}``. + connectivity_type: type of inter-qubit connectivity, ``{"full", "line", "star"}``. + depth: depth of the CNOT-network, i.e. the number of layers, where each layer consists of + a single CNOT-block; default value will be selected, if ``L <= 0``. + + Returns: + A matrix of size ``(2, N)`` matrix that defines layers in cnot-network, where ``N`` + is either equal ``L``, or defined by a concrete type of the network. + + Raises: + ValueError: if unsupported type of CNOT-network layout or number of qubits or combination + of parameters are passed. + """ + if num_qubits < 2: + raise ValueError("Number of qubits must be greater or equal to 2") + + if depth <= 0: + new_depth = _lower_limit(num_qubits) + logger.debug( + "Number of CNOT units chosen as the lower limit: %d, got a non-positive value: %d", + new_depth, + depth, + ) + depth = new_depth + + if network_layout == "sequ": + links = _get_connectivity(num_qubits=num_qubits, connectivity=connectivity_type) + return _sequential_network(num_qubits=num_qubits, links=links, depth=depth) + + elif network_layout == "spin": + return _spin_network(num_qubits=num_qubits, depth=depth) + + elif network_layout == "cart": + cnots = _cartan_network(num_qubits=num_qubits) + logger.debug( + "Optimal lower bound: %d; Cartan CNOTs: %d", _lower_limit(num_qubits), cnots.shape[1] + ) + return cnots + + elif network_layout == "cyclic_spin": + if connectivity_type != "full": + raise ValueError(f"'{network_layout}' layout expects 'full' connectivity") + + return _cyclic_spin_network(num_qubits, depth) + + elif network_layout == "cyclic_line": + if connectivity_type != "line": + raise ValueError(f"'{network_layout}' layout expects 'line' connectivity") + + return _cyclic_line_network(num_qubits, depth) + else: + raise ValueError( + f"Unknown type of CNOT-network layout, expects one of {_NETWORK_LAYOUTS}, " + f"got {network_layout}" + ) + + +def _get_connectivity(num_qubits: int, connectivity: str) -> dict: + """ + Generates connectivity structure between qubits. + + Args: + num_qubits: number of qubits. + connectivity: type of connectivity structure, ``{"full", "line", "star"}``. + + Returns: + dictionary of allowed links between qubits. + + Raises: + ValueError: if unsupported type of CNOT-network layout is passed. + """ + if num_qubits == 1: + links = {0: [0]} + + elif connectivity == "full": + # Full connectivity between qubits. + links = {i: list(range(num_qubits)) for i in range(num_qubits)} + + elif connectivity == "line": + # Every qubit is connected to its immediate neighbours only. + links = {i: [i - 1, i, i + 1] for i in range(1, num_qubits - 1)} + + # first qubit + links[0] = [0, 1] + + # last qubit + links[num_qubits - 1] = [num_qubits - 2, num_qubits - 1] + + elif connectivity == "star": + # Every qubit is connected to the first one only. + links = {i: [0, i] for i in range(1, num_qubits)} + + # first qubit + links[0] = list(range(num_qubits)) + + else: + raise ValueError( + f"Unknown connectivity type, expects one of {_CONNECTIVITY_TYPES}, got {connectivity}" + ) + return links + + +def _sequential_network(num_qubits: int, links: dict, depth: int) -> np.ndarray: + """ + Generates a sequential network. + + Args: + num_qubits: number of qubits. + links: dictionary of connectivity links. + depth: depth of the network (number of layers of building blocks). + + Returns: + A matrix of ``(2, N)`` that defines layers in qubit network. + """ + layer = 0 + cnots = np.zeros((2, depth), dtype=int) + while True: + for i in range(0, num_qubits - 1): + for j in range(i + 1, num_qubits): + if j in links[i]: + cnots[0, layer] = i + cnots[1, layer] = j + layer += 1 + if layer >= depth: + return cnots + + +def _spin_network(num_qubits: int, depth: int) -> np.ndarray: + """ + Generates a spin-like network. + + Args: + num_qubits: number of qubits. + depth: depth of the network (number of layers of building blocks). + + Returns: + A matrix of size ``2 x L`` that defines layers in qubit network. + """ + layer = 0 + cnots = np.zeros((2, depth), dtype=int) + while True: + for i in range(0, num_qubits - 1, 2): + cnots[0, layer] = i + cnots[1, layer] = i + 1 + layer += 1 + if layer >= depth: + return cnots + + for i in range(1, num_qubits - 1, 2): + cnots[0, layer] = i + cnots[1, layer] = i + 1 + layer += 1 + if layer >= depth: + return cnots + + +def _cartan_network(num_qubits: int) -> np.ndarray: + """ + Cartan decomposition in a recursive way, starting from n = 3. + + Args: + num_qubits: number of qubits. + + Returns: + 2xN matrix that defines layers in qubit network, where N is the + depth of Cartan decomposition. + + Raises: + ValueError: if number of qubits is less than 3. + """ + n = num_qubits + if n > 3: + cnots = np.asarray([[0, 0, 0], [1, 1, 1]]) + mult = np.asarray([[n - 2, n - 3, n - 2, n - 3], [n - 1, n - 1, n - 1, n - 1]]) + for _ in range(n - 2): + cnots = np.hstack((np.tile(np.hstack((cnots, mult)), 3), cnots)) + mult[0, -1] -= 1 + mult = np.tile(mult, 2) + elif n == 3: + cnots = np.asarray( + [ + [0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0], + [1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1], + ] + ) + else: + raise ValueError(f"The number of qubits must be >= 3, got {n}.") + + return cnots + + +def _cyclic_spin_network(num_qubits: int, depth: int) -> np.ndarray: + """ + Same as in the spin-like network, but the first and the last qubits are also connected. + + Args: + num_qubits: number of qubits. + depth: depth of the network (number of layers of building blocks). + + Returns: + A matrix of size ``2 x L`` that defines layers in qubit network. + """ + + cnots = np.zeros((2, depth), dtype=int) + z = 0 + while True: + for i in range(0, num_qubits, 2): + if i + 1 <= num_qubits - 1: + cnots[0, z] = i + cnots[1, z] = i + 1 + z += 1 + if z >= depth: + return cnots + + for i in range(1, num_qubits, 2): + if i + 1 <= num_qubits - 1: + cnots[0, z] = i + cnots[1, z] = i + 1 + z += 1 + elif i == num_qubits - 1: + cnots[0, z] = i + cnots[1, z] = 0 + z += 1 + if z >= depth: + return cnots + + +def _cyclic_line_network(num_qubits: int, depth: int) -> np.ndarray: + """ + Generates a line based CNOT structure. + + Args: + num_qubits: number of qubits. + depth: depth of the network (number of layers of building blocks). + + Returns: + A matrix of size ``2 x L`` that defines layers in qubit network. + """ + + cnots = np.zeros((2, depth), dtype=int) + for i in range(depth): + cnots[0, i] = (i + 0) % num_qubits + cnots[1, i] = (i + 1) % num_qubits + return cnots diff --git a/qiskit/synthesis/unitary/aqc/cnot_unit_circuit.py b/qiskit/synthesis/unitary/aqc/cnot_unit_circuit.py new file mode 100644 index 000000000000..ed0de7e6bffc --- /dev/null +++ b/qiskit/synthesis/unitary/aqc/cnot_unit_circuit.py @@ -0,0 +1,103 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2023. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +""" +This is the Parametric Circuit class: anything that you need for a circuit +to be parametrized and used for approximate compiling optimization. +""" +from __future__ import annotations +from typing import Optional + +import numpy as np + +from .approximate import ApproximateCircuit + + +class CNOTUnitCircuit(ApproximateCircuit): + """A class that represents an approximate circuit based on CNOT unit blocks.""" + + def __init__( + self, + num_qubits: int, + cnots: np.ndarray, + tol: Optional[float] = 0.0, + name: Optional[str] = None, + ) -> None: + """ + Args: + num_qubits: the number of qubits in this circuit. + cnots: an array of dimensions ``(2, L)`` indicating where the CNOT units will be placed. + tol: angle parameter less or equal this (small) value is considered equal zero and + corresponding gate is not inserted into the output circuit (because it becomes + identity one in this case). + name: name of this circuit + + Raises: + ValueError: if an unsupported parameter is passed. + """ + super().__init__(num_qubits=num_qubits, name=name) + + if cnots.ndim != 2 or cnots.shape[0] != 2: + raise ValueError("CNOT structure must be defined as an array of the size (2, N)") + + self._cnots = cnots + self._num_cnots = cnots.shape[1] + self._tol = tol + + # Thetas to be optimized by the AQC algorithm + self._thetas: np.ndarray | None = None + + @property + def thetas(self) -> np.ndarray: + """ + Returns a vector of rotation angles used by CNOT units in this circuit. + + Returns: + Parameters of the rotation gates in this circuit. + """ + return self._thetas + + def build(self, thetas: np.ndarray) -> None: + """ + Constructs a Qiskit quantum circuit out of the parameters (angles) of this circuit. If a + parameter value is less in absolute value than the specified tolerance then the + corresponding rotation gate will be skipped in the circuit. + """ + n = self.num_qubits + self._thetas = thetas + cnots = self._cnots + + for k in range(n): + # add initial three rotation gates for each qubit + p = 4 * self._num_cnots + 3 * k + k = n - k - 1 + if np.abs(thetas[2 + p]) > self._tol: + self.rz(thetas[2 + p], k) + if np.abs(thetas[1 + p]) > self._tol: + self.ry(thetas[1 + p], k) + if np.abs(thetas[0 + p]) > self._tol: + self.rz(thetas[0 + p], k) + + for c in range(self._num_cnots): + p = 4 * c + # Extract where the CNOT goes + q1 = n - 1 - int(cnots[0, c]) + q2 = n - 1 - int(cnots[1, c]) + # Construct a CNOT unit + self.cx(q1, q2) + if np.abs(thetas[0 + p]) > self._tol: + self.ry(thetas[0 + p], q1) + if np.abs(thetas[1 + p]) > self._tol: + self.rz(thetas[1 + p], q1) + if np.abs(thetas[2 + p]) > self._tol: + self.ry(thetas[2 + p], q2) + if np.abs(thetas[3 + p]) > self._tol: + self.rx(thetas[3 + p], q2) diff --git a/qiskit/synthesis/unitary/aqc/cnot_unit_objective.py b/qiskit/synthesis/unitary/aqc/cnot_unit_objective.py new file mode 100644 index 000000000000..b260e9062677 --- /dev/null +++ b/qiskit/synthesis/unitary/aqc/cnot_unit_objective.py @@ -0,0 +1,299 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2023. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +""" +A definition of the approximate circuit compilation optimization problem based on CNOT unit +definition. +""" +from __future__ import annotations +import typing +from abc import ABC + +import numpy as np +from numpy import linalg as la + +from .approximate import ApproximatingObjective +from .elementary_operations import ry_matrix, rz_matrix, place_unitary, place_cnot, rx_matrix + + +class CNOTUnitObjective(ApproximatingObjective, ABC): + """ + A base class for a problem definition based on CNOT unit. This class may have different + subclasses for objective and gradient computations. + """ + + def __init__(self, num_qubits: int, cnots: np.ndarray) -> None: + """ + Args: + num_qubits: number of qubits. + cnots: a CNOT structure to be used in the optimization procedure. + """ + super().__init__() + self._num_qubits = num_qubits + self._cnots = cnots + self._num_cnots = cnots.shape[1] + + @property + def num_cnots(self): + """ + Returns: + A number of CNOT units to be used by the approximate circuit. + """ + return self._num_cnots + + @property + def num_thetas(self): + """ + Returns: + Number of parameters (angles) of rotation gates in this circuit. + """ + return 3 * self._num_qubits + 4 * self._num_cnots + + +class DefaultCNOTUnitObjective(CNOTUnitObjective): + """A naive implementation of the objective function based on CNOT units.""" + + def __init__(self, num_qubits: int, cnots: np.ndarray) -> None: + """ + Args: + num_qubits: number of qubits. + cnots: a CNOT structure to be used in the optimization procedure. + """ + super().__init__(num_qubits, cnots) + + # last objective computations to be re-used by gradient + self._last_thetas: np.ndarray | None = None + self._cnot_right_collection: np.ndarray | None = None + self._cnot_left_collection: np.ndarray | None = None + self._rotation_matrix: int | np.ndarray | None = None + self._cnot_matrix: np.ndarray | None = None + + def objective(self, param_values: np.ndarray) -> typing.SupportsFloat: + # rename parameters just to make shorter and make use of our dictionary + thetas = param_values + n = self._num_qubits + d = int(2**n) + cnots = self._cnots + num_cnots = self.num_cnots + + # to save intermediate computations we define the following matrices + # this is the collection of cnot unit matrices ordered from left to + # right as in the circuit, not matrix product + cnot_unit_collection = np.zeros((d, d * num_cnots), dtype=complex) + # this is the collection of matrix products of the cnot units up + # to the given position from the right of the circuit + cnot_right_collection = np.zeros((d, d * num_cnots), dtype=complex) + # this is the collection of matrix products of the cnot units up + # to the given position from the left of the circuit + cnot_left_collection = np.zeros((d, d * num_cnots), dtype=complex) + # first, we construct each cnot unit matrix + for cnot_index in range(num_cnots): + theta_index = 4 * cnot_index + + # cnot qubit indices for the cnot unit identified by cnot_index + q1 = int(cnots[0, cnot_index]) + q2 = int(cnots[1, cnot_index]) + + # rotations that are applied on the q1 qubit + ry1 = ry_matrix(thetas[0 + theta_index]) + rz1 = rz_matrix(thetas[1 + theta_index]) + + # rotations that are applied on the q2 qubit + ry2 = ry_matrix(thetas[2 + theta_index]) + rx2 = rx_matrix(thetas[3 + theta_index]) + + # combine the rotations on qubits q1 and q2 + single_q1 = np.dot(rz1, ry1) + single_q2 = np.dot(rx2, ry2) + + # we place single qubit matrices at the corresponding locations in the (2^n, 2^n) matrix + full_q1 = place_unitary(single_q1, n, q1) + full_q2 = place_unitary(single_q2, n, q2) + + # we place a cnot matrix at the qubits q1 and q2 in the full matrix + cnot_q1q2 = place_cnot(n, q1, q2) + + # compute the cnot unit matrix and store in cnot_unit_collection + cnot_unit_collection[:, d * cnot_index : d * (cnot_index + 1)] = la.multi_dot( + [full_q2, full_q1, cnot_q1q2] + ) + + # this is the matrix corresponding to the intermediate matrix products + # it will end up being the matrix product of all the cnot unit matrices + # first we multiply from the right-hand side of the circuit + cnot_matrix = np.eye(d) + for cnot_index in range(num_cnots - 1, -1, -1): + cnot_matrix = np.dot( + cnot_matrix, cnot_unit_collection[:, d * cnot_index : d * (cnot_index + 1)] + ) + cnot_right_collection[:, d * cnot_index : d * (cnot_index + 1)] = cnot_matrix + # now we multiply from the left-hand side of the circuit + cnot_matrix = np.eye(d) + for cnot_index in range(num_cnots): + cnot_matrix = np.dot( + cnot_unit_collection[:, d * cnot_index : d * (cnot_index + 1)], cnot_matrix + ) + cnot_left_collection[:, d * cnot_index : d * (cnot_index + 1)] = cnot_matrix + + # this is the matrix corresponding to the initial rotations + # we start with 1 and kronecker product each qubit's rotations + rotation_matrix: int | np.ndarray = 1 + for q in range(n): + theta_index = 4 * num_cnots + 3 * q + rz0 = rz_matrix(thetas[0 + theta_index]) + ry1 = ry_matrix(thetas[1 + theta_index]) + rz2 = rz_matrix(thetas[2 + theta_index]) + rotation_matrix = np.kron(rotation_matrix, la.multi_dot([rz0, ry1, rz2])) + + # the matrix corresponding to the full circuit is the cnot part and + # rotation part multiplied together + circuit_matrix = np.dot(cnot_matrix, rotation_matrix) + + # compute error + error = 0.5 * (la.norm(circuit_matrix - self._target_matrix, "fro") ** 2) + + # cache computations for gradient + self._last_thetas = thetas + self._cnot_left_collection = cnot_left_collection + self._cnot_right_collection = cnot_right_collection + self._rotation_matrix = rotation_matrix + self._cnot_matrix = cnot_matrix + + return error + + def gradient(self, param_values: np.ndarray) -> np.ndarray: + # just to make shorter + thetas = param_values + # if given thetas are the same as used at the previous objective computations, then + # we re-use computations, otherwise we have to re-compute objective + if not np.all(np.isclose(thetas, self._last_thetas)): + self.objective(thetas) + + # the partial derivative of the circuit with respect to an angle + # is the same circuit with the corresponding pauli gate, multiplied + # by a global phase of -1j / 2, next to the rotation gate (it commutes) + pauli_x = np.multiply(-1j / 2, np.asarray([[0, 1], [1, 0]])) + pauli_y = np.multiply(-1j / 2, np.asarray([[0, -1j], [1j, 0]])) + pauli_z = np.multiply(-1j / 2, np.asarray([[1, 0], [0, -1]])) + + n = self._num_qubits + d = int(2**n) + cnots = self._cnots + num_cnots = self.num_cnots + + # the partial derivative of the cost function is -Re + # where V' is the partial derivative of the circuit + # first we compute the partial derivatives in the cnot part + der = np.zeros(4 * num_cnots + 3 * n) + for cnot_index in range(num_cnots): + theta_index = 4 * cnot_index + + # cnot qubit indices for the cnot unit identified by cnot_index + q1 = int(cnots[0, cnot_index]) + q2 = int(cnots[1, cnot_index]) + + # rotations that are applied on the q1 qubit + ry1 = ry_matrix(thetas[0 + theta_index]) + rz1 = rz_matrix(thetas[1 + theta_index]) + + # rotations that are applied on the q2 qubit + ry2 = ry_matrix(thetas[2 + theta_index]) + rx2 = rx_matrix(thetas[3 + theta_index]) + + # combine the rotations on qubits q1 and q2 + # note we have to insert an extra pauli gate to take the derivative + # of the appropriate rotation gate + for i in range(4): + if i == 0: + single_q1 = la.multi_dot([rz1, pauli_y, ry1]) + single_q2 = np.dot(rx2, ry2) + elif i == 1: + single_q1 = la.multi_dot([pauli_z, rz1, ry1]) + single_q2 = np.dot(rx2, ry2) + elif i == 2: + single_q1 = np.dot(rz1, ry1) + single_q2 = la.multi_dot([rx2, pauli_y, ry2]) + else: + single_q1 = np.dot(rz1, ry1) + single_q2 = la.multi_dot([pauli_x, rx2, ry2]) + + # we place single qubit matrices at the corresponding locations in + # the (2^n, 2^n) matrix + full_q1 = place_unitary(single_q1, n, q1) + full_q2 = place_unitary(single_q2, n, q2) + + # we place a cnot matrix at the qubits q1 and q2 in the full matrix + cnot_q1q2 = place_cnot(n, q1, q2) + + # partial derivative of that particular cnot unit, size of (2^n, 2^n) + der_cnot_unit = la.multi_dot([full_q2, full_q1, cnot_q1q2]) + # der_cnot_unit is multiplied by the matrix product of cnot units to the left + # of it (if there are any) and to the right of it (if there are any) + if cnot_index == 0: + der_cnot_matrix = np.dot( + self._cnot_right_collection[:, d : 2 * d], + der_cnot_unit, + ) + elif num_cnots - 1 == cnot_index: + der_cnot_matrix = np.dot( + der_cnot_unit, + self._cnot_left_collection[:, d * (num_cnots - 2) : d * (num_cnots - 1)], + ) + else: + der_cnot_matrix = la.multi_dot( + [ + self._cnot_right_collection[ + :, d * (cnot_index + 1) : d * (cnot_index + 2) + ], + der_cnot_unit, + self._cnot_left_collection[:, d * (cnot_index - 1) : d * cnot_index], + ] + ) + + # the matrix corresponding to the full circuit partial derivative + # is the partial derivative of the cnot part multiplied by the usual + # rotation part + der_circuit_matrix = np.dot(der_cnot_matrix, self._rotation_matrix) + # we compute the partial derivative of the cost function + der[i + theta_index] = -np.real( + np.trace(np.dot(der_circuit_matrix.conj().T, self._target_matrix)) + ) + + # now we compute the partial derivatives in the rotation part + # we start with 1 and kronecker product each qubit's rotations + for i in range(3 * n): + der_rotation_matrix: int | np.ndarray = 1 + for q in range(n): + theta_index = 4 * num_cnots + 3 * q + rz0 = rz_matrix(thetas[0 + theta_index]) + ry1 = ry_matrix(thetas[1 + theta_index]) + rz2 = rz_matrix(thetas[2 + theta_index]) + # for the appropriate rotation gate that we are taking + # the partial derivative of, we have to insert the + # corresponding pauli matrix + if i - 3 * q == 0: + rz0 = np.dot(pauli_z, rz0) + elif i - 3 * q == 1: + ry1 = np.dot(pauli_y, ry1) + elif i - 3 * q == 2: + rz2 = np.dot(pauli_z, rz2) + der_rotation_matrix = np.kron(der_rotation_matrix, la.multi_dot([rz0, ry1, rz2])) + + # the matrix corresponding to the full circuit partial derivative + # is the usual cnot part multiplied by the partial derivative of + # the rotation part + der_circuit_matrix = np.dot(self._cnot_matrix, der_rotation_matrix) + # we compute the partial derivative of the cost function + der[4 * num_cnots + i] = -np.real( + np.trace(np.dot(der_circuit_matrix.conj().T, self._target_matrix)) + ) + + return der diff --git a/qiskit/synthesis/unitary/aqc/elementary_operations.py b/qiskit/synthesis/unitary/aqc/elementary_operations.py new file mode 100644 index 000000000000..65e6bea5f050 --- /dev/null +++ b/qiskit/synthesis/unitary/aqc/elementary_operations.py @@ -0,0 +1,108 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2023. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +""" +These are a number of elementary functions that are required for the AQC routines to work. +""" + +import numpy as np + +from qiskit.circuit.library import RXGate, RZGate, RYGate + + +def place_unitary(unitary: np.ndarray, n: int, j: int) -> np.ndarray: + """ + Computes I(j - 1) tensor product U tensor product I(n - j), where U is a unitary matrix + of size ``(2, 2)``. + + Args: + unitary: a unitary matrix of size ``(2, 2)``. + n: num qubits. + j: position where to place a unitary. + + Returns: + a unitary of n qubits with u in position j. + """ + return np.kron(np.kron(np.eye(2**j), unitary), np.eye(2 ** (n - 1 - j))) + + +def place_cnot(n: int, j: int, k: int) -> np.ndarray: + """ + Places a CNOT from j to k. + + Args: + n: number of qubits. + j: control qubit. + k: target qubit. + + Returns: + a unitary of n qubits with CNOT placed at ``j`` and ``k``. + """ + if j < k: + unitary = np.kron( + np.kron(np.eye(2**j), [[1, 0], [0, 0]]), np.eye(2 ** (n - 1 - j)) + ) + np.kron( + np.kron( + np.kron(np.kron(np.eye(2**j), [[0, 0], [0, 1]]), np.eye(2 ** (k - j - 1))), + [[0, 1], [1, 0]], + ), + np.eye(2 ** (n - 1 - k)), + ) + else: + unitary = np.kron( + np.kron(np.eye(2**j), [[1, 0], [0, 0]]), np.eye(2 ** (n - 1 - j)) + ) + np.kron( + np.kron( + np.kron(np.kron(np.eye(2**k), [[0, 1], [1, 0]]), np.eye(2 ** (j - k - 1))), + [[0, 0], [0, 1]], + ), + np.eye(2 ** (n - 1 - j)), + ) + return unitary + + +def rx_matrix(phi: float) -> np.ndarray: + """ + Computes an RX rotation by the angle of ``phi``. + + Args: + phi: rotation angle. + + Returns: + an RX rotation matrix. + """ + return RXGate(phi).to_matrix() + + +def ry_matrix(phi: float) -> np.ndarray: + """ + Computes an RY rotation by the angle of ``phi``. + + Args: + phi: rotation angle. + + Returns: + an RY rotation matrix. + """ + return RYGate(phi).to_matrix() + + +def rz_matrix(phi: float) -> np.ndarray: + """ + Computes an RZ rotation by the angle of ``phi``. + + Args: + phi: rotation angle. + + Returns: + an RZ rotation matrix. + """ + return RZGate(phi).to_matrix() diff --git a/qiskit/synthesis/unitary/aqc/fast_gradient/__init__.py b/qiskit/synthesis/unitary/aqc/fast_gradient/__init__.py new file mode 100644 index 000000000000..45496159f84d --- /dev/null +++ b/qiskit/synthesis/unitary/aqc/fast_gradient/__init__.py @@ -0,0 +1,164 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2023. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +r""" +================================================================================ +Fast implementation of objective function class +(:mod:`qiskit.synthesis.unitary.aqc.fast_gradient`) +================================================================================ + +.. currentmodule:: qiskit.synthesis.unitary.aqc.fast_gradient + +Extension to the implementation of Approximate Quantum Compiler as described in the paper [1]. + +Interface +========= + +The main public class of this module is FastCNOTUnitObjective. It replaces the default objective +function implementation :class:`.DefaultCNOTUnitObjective` for faster computation. +The individual classes include the public one (FastCNOTUnitObjective) and few +internal ones: + +.. autosummary:: + :toctree: ../stubs + :template: autosummary/class_no_inherited_members.rst + + FastCNOTUnitObjective + LayerBase + Layer1Q + Layer2Q + PMatrix + + +Mathematical Details +==================== + +In what follows we briefly outline the main ideas underlying the accelerated implementation +of objective function class. + +* The key ingredient of approximate compiling is the efficient optimization procedure + that minimizes :math:`\|V - U\|_{\mathrm{F}}` on a classical computer, where :math:`U` + is a given (target) unitary matrix and :math:`V` is a matrix of approximating quantum + circuit. Alternatively, we maximize the Hilbert-Schmidt product between :math:`U` and + :math:`V` as outlined in the main part of the documentation. + +* The circuit :math:`V` can be represented as a sequence of 2-qubit gates (layers) + applied one after another. The corresponding matrix takes the form: + :math:`V = C_0 C_1 \ldots C_{L-1} F`, where :math:`L` is the length of the sequence + (number of layers). If the total number of qubits :math:`n > 2`, every + :math:`C_i = C_i(\Theta_i)` is a sparse, :math:`2^n \times 2^n` matrix of 2-qubit gate + (CNOT unit block) parameterized by a sub-set of parameters :math:`\Theta_i` + (4 parameters per unit block), and :math:`F` is a matrix that comprises the action + of all 1-qubit gates in front of approximating circuit. See the paper [1] for details. + +* Over the course of optimization we compute the value of objective function and its + gradient, which implies computation of :math:`V` and its derivatives + :math:`{\partial V}/{\partial \Theta_i}` for all :math:`i`, given the current estimation + of all the parameters :math:`\Theta`. + +* A naive implementation of the product :math:`V = C_0 C_1 \ldots C_{L-1} F` and its + derivatives would include computation and memorization of forward and backward partial + products as required by the backtracking algorithm. This is wasteful in terms of + performance and resource allocation. + +* Minimization of :math:`\|V - U\|_{\mathrm{F}}^2` is equivalent to maximization of + :math:`\text{Re}\left(\text{Tr}\left(U^{\dagger} V\right)\right)`. By cyclic permutation + of the sequence of matrices under trace operation, we can avoid memorization of intermediate + partial products of gate matrices :math:`C_i`. Note, matrix size grows exponentially with + the number of qubits, quickly becoming prohibitively large. + +* Sparse structure of :math:`C_i` can be exploited to speed up matrix-matrix multiplication. + However, using sparse matrices as such does not give performance gain because sparse patterns + tend to violate data proximity inside the cache memory of modern CPUs. Instead, we make use + of special structure of gate matrices :math:`C_i` coupled with permutation ones. Although + permutation is not cache friendly either, its impact is seemingly less severe than that + of sparse matrix multiplication (at least in Python implementation). + +* On every optimization iteration we, first, compute :math:`V = C_0 C_1 \ldots C_{L-1} F` + given the current estimation of all the parameters :math:`\Theta`. + +* As for the gradient of objective function, it can be shown (by moving cyclically around + an individual matrices under trace operation) that: + +.. math:: + \text{Tr}\left( U^{\dagger} \frac{\partial V}{\partial \Theta_{l,k}} \right) = + \langle \text{vec}\left(E_l\right), \text{vec}\left( + \frac{\partial C_l}{\partial \Theta_{l,k}}\right) \rangle, + +where :math:`\Theta_{l,k}` is a :math:`k`-th parameter of :math:`l`-th CNOT unit block, +and :math:`E_l=C_{l-1}\left(C_{l-2}\left(\cdots\left(C_0\left(U^{\dagger}V +C_0^{\dagger}\right)C_1^{\dagger}\right) \cdots\right)C_{l-1}^{\dagger}\right)C_l^{\dagger}` +is an intermediate matrix. + +* For every :math:`l`-th gradient component, we compute the trace using the matrix + :math:`E_l`, then this matrix is updated by multiplication on left and on the right + by corresponding gate matrices :math:`C_l` and :math:`C_{l+1}^{\dagger}` respectively + and proceed to the next gradient component. + +* We save computations and resources by not storing intermediate partial products of + :math:`C_i`. Instead, incrementally updated matrix :math:`E_l` keeps all related + information. Also, vectorization of involved matrices (see the above formula) allows + us to replace matrix-matrix multiplication by "cheaper" vector-vector one under the + trace operation. + +* The matrices :math:`C_i` are sparse. However, even for relatively small matrices + (< 1M elements) sparse-dense multiplication can be very slow. Construction of sparse + matrices takes a time as well. We should update every gate matrix on each iteration + of optimization loop. + +* In fact, any gate matrix :math:`C_i` can be transformed to what we call a standard + form: :math:`C_i = P^T \widetilde{C}_i P`, where :math:`P` is an easily computable + permutation matrix and :math:`\widetilde{C}_i` has a block-diagonal layout: + +.. math:: + \widetilde{C}_i = \left( + \begin{array}{ccc} + G_{4 \times 4} & \ddots & 0 \\ + \ddots & \ddots & \ddots \\ + 0 & \ddots & G_{4 \times 4} + \end{array} + \right) + +* The 2-qubit gate matrix :math:`G_{4 \times 4}` is repeated along diagonal of the full + :math:`2^n \times 2^n` :math:`\widetilde{C}_i`. + +* We do not actually create neither matrix :math:`\widetilde{C}_i` nor :math:`P`. + In fact, only :math:`G_{4 \times 4}` and a permutation array (of size :math:`2^n`) + are kept in memory. + +* Consider left-hand side multiplication by some dense, :math:`2^n \times 2^n` matrix :math:`M`: + +.. math:: + C_i M = P^T \widetilde{C}_i P M = P^T \left( \widetilde{C}_i \left( P M \right) \right) + +* First, we permute rows of :math:`M`, which is equivalent to the product :math:`P M`, but + without expensive multiplication of two :math:`2^n \times 2^n` matrices. + +* Second, we compute :math:`\widetilde{C}_i P M` multiplying every block-diagonal sub-matrix + :math:`G_{4 \times 4}` by the corresponding rows of :math:`P M`. This is the dense-dense + matrix multiplication, which is very well optimized on modern CPUs. Important: the total + number of operations is :math:`O(2^{2 n})` in contrast to :math:`O(2^{3 n})` as in general + case. + +* Third, we permute rows of :math:`\widetilde{C}_i P M` by applying :math:`P^T`. + +* Right-hand side multiplication is done in a similar way. + +* In summary, we save computational resources by exploiting some properties of 2-qubit gate + matrices :math:`C_i` and using hardware optimized multiplication of dense matrices. There + is still a room for further improvement, of course. + +References: + + [1]: Liam Madden, Andrea Simonetto, Best Approximate Quantum Compiling Problems. + `arXiv:2106.05649 `_ +""" diff --git a/qiskit/synthesis/unitary/aqc/fast_gradient/fast_grad_utils.py b/qiskit/synthesis/unitary/aqc/fast_gradient/fast_grad_utils.py new file mode 100644 index 000000000000..1c6e480fb820 --- /dev/null +++ b/qiskit/synthesis/unitary/aqc/fast_gradient/fast_grad_utils.py @@ -0,0 +1,237 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2023. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +""" +Utility functions in the fast gradient implementation. +""" +from __future__ import annotations +from typing import Union +import numpy as np + + +def is_permutation(x: np.ndarray) -> bool: + """ + Checks if array is really an index permutation. + + Args: + 1D-array of integers that supposedly represents a permutation. + + Returns: + True, if array is really a permutation of indices. + """ + return ( + isinstance(x, np.ndarray) + and x.ndim == 1 + and x.dtype == np.int64 + and np.all(np.sort(x) == np.arange(x.size, dtype=np.int64)) + ) + + +def reverse_bits(x: Union[int, np.ndarray], nbits: int, enable: bool) -> Union[int, np.ndarray]: + """ + Reverses the bit order in a number of ``nbits`` length. + If ``x`` is an array, then operation is applied to every entry. + + Args: + x: either a single integer or an array of integers. + nbits: number of meaningful bits in the number x. + enable: apply reverse operation, if enabled, otherwise leave unchanged. + + Returns: + a number or array of numbers with reversed bits. + """ + + if not enable: + if isinstance(x, int): + pass + else: + x = x.copy() + return x + + if isinstance(x, int): + res: int | np.ndarray = int(0) + else: + x = x.copy() + res = np.full_like(x, fill_value=0) + + for _ in range(nbits): + res <<= 1 + res |= x & 1 + x >>= 1 + return res + + +def swap_bits(num: int, a: int, b: int) -> int: + """ + Swaps the bits at positions 'a' and 'b' in the number 'num'. + + Args: + num: an integer number where bits should be swapped. + a: index of the first bit to be swapped. + b: index of the second bit to be swapped. + + Returns: + the number with swapped bits. + """ + x = ((num >> a) ^ (num >> b)) & 1 + return num ^ ((x << a) | (x << b)) + + +def bit_permutation_1q(n: int, k: int) -> np.ndarray: + """ + Constructs index permutation that brings a circuit consisting of a single + 1-qubit gate to "standard form": ``kron(I(2^n/2), G)``, as we call it. Here n + is the number of qubits, ``G`` is a 2x2 gate matrix, ``I(2^n/2)`` is the identity + matrix of size ``(2^n/2)x(2^n/2)``, and the full size of the circuit matrix is + ``(2^n)x(2^n)``. Circuit matrix in standard form becomes block-diagonal (with + sub-matrices ``G`` on the main diagonal). Multiplication of such a matrix and + a dense one is much faster than generic dense-dense product. Moreover, + we do not need to keep the entire circuit matrix in memory but just 2x2 ``G`` + one. This saves a lot of memory when the number of qubits is large. + + Args: + n: number of qubits. + k: index of qubit where single 1-qubit gate is applied. + + Returns: + permutation that brings the whole layer to the standard form. + """ + perm = np.arange(2**n, dtype=np.int64) + if k != n - 1: + for v in range(2**n): + perm[v] = swap_bits(v, k, n - 1) + return perm + + +def bit_permutation_2q(n: int, j: int, k: int) -> np.ndarray: + """ + Constructs index permutation that brings a circuit consisting of a single + 2-qubit gate to "standard form": ``kron(I(2^n/4), G)``, as we call it. Here ``n`` + is the number of qubits, ``G`` is a 4x4 gate matrix, ``I(2^n/4)`` is the identity + matrix of size ``(2^n/4)x(2^n/4)``, and the full size of the circuit matrix is + ``(2^n)x(2^n)``. Circuit matrix in standard form becomes block-diagonal (with + sub-matrices ``G`` on the main diagonal). Multiplication of such a matrix and + a dense one is much faster than generic dense-dense product. Moreover, + we do not need to keep the entire circuit matrix in memory but just 4x4 ``G`` + one. This saves a lot of memory when the number of qubits is large. + + Args: + n: number of qubits. + j: index of control qubit where single 2-qubit gate is applied. + k: index of target qubit where single 2-qubit gate is applied. + + Returns: + permutation that brings the whole layer to the standard form. + """ + dim = 2**n + perm = np.arange(dim, dtype=np.int64) + if j < n - 2: + if k < n - 2: + for v in range(dim): + perm[v] = swap_bits(swap_bits(v, j, n - 2), k, n - 1) + elif k == n - 2: + for v in range(dim): + perm[v] = swap_bits(swap_bits(v, n - 2, n - 1), j, n - 2) + else: + for v in range(dim): + perm[v] = swap_bits(v, j, n - 2) + elif j == n - 2: + if k < n - 2: + for v in range(dim): + perm[v] = swap_bits(v, k, n - 1) + else: + pass + else: + if k < n - 2: + for v in range(dim): + perm[v] = swap_bits(swap_bits(v, n - 2, n - 1), k, n - 1) + else: + for v in range(dim): + perm[v] = swap_bits(v, n - 2, n - 1) + return perm + + +def inverse_permutation(perm: np.ndarray) -> np.ndarray: + """ + Returns inverse permutation. + + Args: + perm: permutation to be reversed. + + Returns: + inverse permutation. + """ + inv = np.zeros_like(perm) + inv[perm] = np.arange(perm.size, dtype=np.int64) + return inv + + +def make_rx(phi: float, out: np.ndarray) -> np.ndarray: + """ + Makes a 2x2 matrix that corresponds to X-rotation gate. + This is a fast implementation that does not allocate the output matrix. + + Args: + phi: rotation angle. + out: placeholder for the result (2x2, complex-valued matrix). + + Returns: + rotation gate, same object as referenced by "out". + """ + a = 0.5 * phi + cs, sn = np.cos(a).item(), -1j * np.sin(a).item() + out[0, 0] = cs + out[0, 1] = sn + out[1, 0] = sn + out[1, 1] = cs + return out + + +def make_ry(phi: float, out: np.ndarray) -> np.ndarray: + """ + Makes a 2x2 matrix that corresponds to Y-rotation gate. + This is a fast implementation that does not allocate the output matrix. + + Args: + phi: rotation angle. + out: placeholder for the result (2x2, complex-valued matrix). + + Returns: + rotation gate, same object as referenced by "out". + """ + a = 0.5 * phi + cs, sn = np.cos(a).item(), np.sin(a).item() + out[0, 0] = cs + out[0, 1] = -sn + out[1, 0] = sn + out[1, 1] = cs + return out + + +def make_rz(phi: float, out: np.ndarray) -> np.ndarray: + """ + Makes a 2x2 matrix that corresponds to Z-rotation gate. + This is a fast implementation that does not allocate the output matrix. + + Args: + phi: rotation angle. + out: placeholder for the result (2x2, complex-valued matrix). + + Returns: + rotation gate, same object as referenced by "out". + """ + exp = np.exp(0.5j * phi).item() + out[0, 0] = 1.0 / exp + out[0, 1] = 0 + out[1, 0] = 0 + out[1, 1] = exp + return out diff --git a/qiskit/synthesis/unitary/aqc/fast_gradient/fast_gradient.py b/qiskit/synthesis/unitary/aqc/fast_gradient/fast_gradient.py new file mode 100644 index 000000000000..3424735a782a --- /dev/null +++ b/qiskit/synthesis/unitary/aqc/fast_gradient/fast_gradient.py @@ -0,0 +1,225 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2023. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +""" +Implementation of the fast objective function class. +""" + +import warnings +import numpy as np + +from .layer import ( + LayerBase, + Layer2Q, + Layer1Q, + init_layer2q_matrices, + init_layer2q_deriv_matrices, + init_layer1q_matrices, + init_layer1q_deriv_matrices, +) +from .pmatrix import PMatrix +from ..cnot_unit_objective import CNOTUnitObjective + + +class FastCNOTUnitObjective(CNOTUnitObjective): + """ + Implementation of objective function and gradient calculator, which is + similar to + :class:`~qiskit.transpiler.aqc.DefaultCNOTUnitObjective` + but several times faster. + """ + + def __init__(self, num_qubits: int, cnots: np.ndarray): + super().__init__(num_qubits, cnots) + + if not 2 <= num_qubits <= 16: + raise ValueError("expects number of qubits in the range [2..16]") + + dim = 2**num_qubits + self._ucf_mat = PMatrix(num_qubits) # U^dagger @ C @ F + self._fuc_mat = PMatrix(num_qubits) # F @ U^dagger @ C + self._circ_thetas = np.zeros((self.num_thetas,)) # last thetas used + + # array of C-layers: + self._c_layers = np.asarray([object()] * self.num_cnots, dtype=LayerBase) + # array of F-layers: + self._f_layers = np.asarray([object()] * num_qubits, dtype=LayerBase) + # 4x4 C-gate matrices: + self._c_gates = np.full((self.num_cnots, 4, 4), fill_value=0, dtype=np.complex128) + # derivatives of 4x4 C-gate matrices: + self._c_dervs = np.full((self.num_cnots, 4, 4, 4), fill_value=0, dtype=np.complex128) + # 4x4 F-gate matrices: + self._f_gates = np.full((num_qubits, 2, 2), fill_value=0, dtype=np.complex128) + # derivatives of 4x4 F-gate matrices: + self._f_dervs = np.full((num_qubits, 3, 2, 2), fill_value=0, dtype=np.complex128) + # temporary NxN matrices: + self._tmp1 = np.full((dim, dim), fill_value=0, dtype=np.complex128) + self._tmp2 = np.full((dim, dim), fill_value=0, dtype=np.complex128) + + # Create layers of 2-qubit gates. + for q in range(self.num_cnots): + j = int(cnots[0, q]) + k = int(cnots[1, q]) + self._c_layers[q] = Layer2Q(num_qubits=num_qubits, j=j, k=k) + + # Create layers of 1-qubit gates. + for k in range(num_qubits): + self._f_layers[k] = Layer1Q(num_qubits=num_qubits, k=k) + + def objective(self, param_values: np.ndarray) -> float: + """ + Computes the objective function and some intermediate data for + the subsequent gradient computation. + See description of the base class method. + """ + depth, n = self.num_cnots, self._num_qubits + + # Memorize the last angular parameters used to compute the objective. + if self._circ_thetas.size == 0: + self._circ_thetas = np.zeros((self.num_thetas,)) + np.copyto(self._circ_thetas, param_values) + + thetas4d = param_values[: 4 * depth].reshape(depth, 4) + thetas3n = param_values[4 * depth :].reshape(n, 3) + + init_layer2q_matrices(thetas=thetas4d, dst=self._c_gates) + init_layer2q_deriv_matrices(thetas=thetas4d, dst=self._c_dervs) + init_layer1q_matrices(thetas=thetas3n, dst=self._f_gates) + init_layer1q_deriv_matrices(thetas=thetas3n, dst=self._f_dervs) + + self._init_layers() + self._calc_ucf_fuc() + objective_value = self._calc_objective_function() + return objective_value + + def gradient(self, param_values: np.ndarray) -> np.ndarray: + """ + Computes the gradient of objective function. + See description of the base class method. + """ + + # If thetas are the same as used for objective value calculation + # before calling this function, then we re-use the computations, + # otherwise we have to re-compute the objective. + tol = float(np.sqrt(np.finfo(float).eps)) + if not np.allclose(param_values, self._circ_thetas, atol=tol, rtol=tol): + self.objective(param_values) + warnings.warn("gradient is computed before the objective") + + grad = np.full((param_values.size,), fill_value=0, dtype=np.float64) + grad4d = grad[: 4 * self.num_cnots].reshape(self.num_cnots, 4) + grad3n = grad[4 * self.num_cnots :].reshape(self._num_qubits, 3) + self._calc_gradient4d(grad4d) + self._calc_gradient3n(grad3n) + return grad + + def _init_layers(self): + """ + Initializes C-layers and F-layers by corresponding gate matrices. + """ + c_gates = self._c_gates + c_layers = self._c_layers + for q in range(self.num_cnots): + c_layers[q].set_from_matrix(mat=c_gates[q]) + + f_gates = self._f_gates + f_layers = self._f_layers + for q in range(self._num_qubits): + f_layers[q].set_from_matrix(mat=f_gates[q]) + + def _calc_ucf_fuc(self): + """ + Computes matrices ``ucf_mat`` and ``fuc_mat``. Both remain non-finalized. + """ + ucf_mat = self._ucf_mat + fuc_mat = self._fuc_mat + tmp1 = self._tmp1 + c_layers = self._c_layers + f_layers = self._f_layers + depth, n = self.num_cnots, self._num_qubits + + # tmp1 = U^dagger. + np.conj(self.target_matrix.T, out=tmp1) + + # ucf_mat = fuc_mat = U^dagger @ C = U^dagger @ C_{depth-1} @ ... @ C_{0}. + self._ucf_mat.set_matrix(tmp1) + for q in range(depth - 1, -1, -1): + ucf_mat.mul_right_q2(c_layers[q], temp_mat=tmp1, dagger=False) + fuc_mat.set_matrix(ucf_mat.finalize(temp_mat=tmp1)) + + # fuc_mat = F @ U^dagger @ C = F_{n-1} @ ... @ F_{0} @ U^dagger @ C. + for q in range(n): + fuc_mat.mul_left_q1(f_layers[q], temp_mat=tmp1) + + # ucf_mat = U^dagger @ C @ F = U^dagger @ C @ F_{n-1} @ ... @ F_{0}. + for q in range(n - 1, -1, -1): + ucf_mat.mul_right_q1(f_layers[q], temp_mat=tmp1, dagger=False) + + def _calc_objective_function(self) -> float: + """ + Computes the value of objective function. + """ + ucf = self._ucf_mat.finalize(temp_mat=self._tmp1) + trace_ucf = np.trace(ucf) + fobj = abs((2**self._num_qubits) - float(np.real(trace_ucf))) + + return fobj + + def _calc_gradient4d(self, grad4d: np.ndarray): + """ + Calculates a part gradient contributed by 2-qubit gates. + """ + fuc = self._fuc_mat + tmp1, tmp2 = self._tmp1, self._tmp2 + c_gates = self._c_gates + c_dervs = self._c_dervs + c_layers = self._c_layers + for q in range(self.num_cnots): + # fuc[q] <-- C[q-1] @ fuc[q-1] @ C[q].conj.T. Note, c_layers[q] has + # been initialized in _init_layers(), however, c_layers[q-1] was + # reused on the previous step, see below, so we need to restore it. + if q > 0: + c_layers[q - 1].set_from_matrix(mat=c_gates[q - 1]) + fuc.mul_left_q2(c_layers[q - 1], temp_mat=tmp1) + fuc.mul_right_q2(c_layers[q], temp_mat=tmp1, dagger=True) + fuc.finalize(temp_mat=tmp1) + # Compute gradient components. We reuse c_layers[q] several times. + for i in range(4): + c_layers[q].set_from_matrix(mat=c_dervs[q, i]) + grad4d[q, i] = (-1) * np.real( + fuc.product_q2(layer=c_layers[q], tmp1=tmp1, tmp2=tmp2) + ) + + def _calc_gradient3n(self, grad3n: np.ndarray): + """ + Calculates a part gradient contributed by 1-qubit gates. + """ + ucf = self._ucf_mat + tmp1, tmp2 = self._tmp1, self._tmp2 + f_gates = self._f_gates + f_dervs = self._f_dervs + f_layers = self._f_layers + for q in range(self._num_qubits): + # ucf[q] <-- F[q-1] @ ucf[q-1] @ F[q].conj.T. Note, f_layers[q] has + # been initialized in _init_layers(), however, f_layers[q-1] was + # reused on the previous step, see below, so we need to restore it. + if q > 0: + f_layers[q - 1].set_from_matrix(mat=f_gates[q - 1]) + ucf.mul_left_q1(f_layers[q - 1], temp_mat=tmp1) + ucf.mul_right_q1(f_layers[q], temp_mat=tmp1, dagger=True) + ucf.finalize(temp_mat=tmp1) + # Compute gradient components. We reuse f_layers[q] several times. + for i in range(3): + f_layers[q].set_from_matrix(mat=f_dervs[q, i]) + grad3n[q, i] = (-1) * np.real( + ucf.product_q1(layer=f_layers[q], tmp1=tmp1, tmp2=tmp2) + ) diff --git a/qiskit/synthesis/unitary/aqc/fast_gradient/layer.py b/qiskit/synthesis/unitary/aqc/fast_gradient/layer.py new file mode 100644 index 000000000000..dab5e9ee09e7 --- /dev/null +++ b/qiskit/synthesis/unitary/aqc/fast_gradient/layer.py @@ -0,0 +1,370 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2023. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +""" +Layer classes for the fast gradient implementation. +""" +from __future__ import annotations +from abc import abstractmethod, ABC +from typing import Optional +import numpy as np +from .fast_grad_utils import ( + bit_permutation_1q, + reverse_bits, + inverse_permutation, + bit_permutation_2q, + make_rz, + make_ry, +) + + +class LayerBase(ABC): + """ + Base class for any layer implementation. Each layer here is represented + by a 2x2 or 4x4 gate matrix ``G`` (applied to 1 or 2 qubits respectively) + interleaved with the identity ones: + ``Layer = I kron I kron ... kron G kron ... kron I kron I`` + """ + + @abstractmethod + def set_from_matrix(self, mat: np.ndarray): + """ + Updates this layer from an external gate matrix. + + Args: + mat: external gate matrix that initializes this layer's one. + """ + raise NotImplementedError() + + @abstractmethod + def get_attr(self) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """ + Returns gate matrix, direct and inverse permutations. + + Returns: + (1) gate matrix; (2) direct permutation; (3) inverse permutations. + """ + raise NotImplementedError() + + +class Layer1Q(LayerBase): + """ + Layer represents a simple circuit where 1-qubit gate matrix (of size 2x2) + interleaves with the identity ones. + """ + + def __init__(self, num_qubits: int, k: int, g2x2: Optional[np.ndarray] = None): + """ + Args: + num_qubits: number of qubits. + k: index of the bit where gate is applied. + g2x2: 2x2 matrix that makes up this layer along with identity ones, + or None (should be set up later). + """ + super().__init__() + + # 2x2 gate matrix (1-qubit gate). + self._gmat = np.full((2, 2), fill_value=0, dtype=np.complex128) + if isinstance(g2x2, np.ndarray): + np.copyto(self._gmat, g2x2) + + bit_flip = True + dim = 2**num_qubits + row_perm = reverse_bits( + bit_permutation_1q(n=num_qubits, k=k), nbits=num_qubits, enable=bit_flip + ) + col_perm = reverse_bits(np.arange(dim, dtype=np.int64), nbits=num_qubits, enable=bit_flip) + self._perm = np.full((dim,), fill_value=0, dtype=np.int64) + self._perm[row_perm] = col_perm + self._inv_perm = inverse_permutation(self._perm) + + def set_from_matrix(self, mat: np.ndarray): + """See base class description.""" + np.copyto(self._gmat, mat) + + def get_attr(self) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """See base class description.""" + return self._gmat, self._perm, self._inv_perm + + +class Layer2Q(LayerBase): + """ + Layer represents a simple circuit where 2-qubit gate matrix (of size 4x4) + interleaves with the identity ones. + """ + + def __init__(self, num_qubits: int, j: int, k: int, g4x4: Optional[np.ndarray] = None): + """ + Args: + num_qubits: number of qubits. + j: index of the first (control) bit. + k: index of the second (target) bit. + g4x4: 4x4 matrix that makes up this layer along with identity ones, + or None (should be set up later). + """ + super().__init__() + + # 4x4 gate matrix (2-qubit gate). + self._gmat = np.full((4, 4), fill_value=0, dtype=np.complex128) + if isinstance(g4x4, np.ndarray): + np.copyto(self._gmat, g4x4) + + bit_flip = True + dim = 2**num_qubits + row_perm = reverse_bits( + bit_permutation_2q(n=num_qubits, j=j, k=k), nbits=num_qubits, enable=bit_flip + ) + col_perm = reverse_bits(np.arange(dim, dtype=np.int64), nbits=num_qubits, enable=bit_flip) + self._perm = np.full((dim,), fill_value=0, dtype=np.int64) + self._perm[row_perm] = col_perm + self._inv_perm = inverse_permutation(self._perm) + + def set_from_matrix(self, mat: np.ndarray): + """See base class description.""" + np.copyto(self._gmat, mat) + + def get_attr(self) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """See base class description.""" + return self._gmat, self._perm, self._inv_perm + + +def init_layer1q_matrices(thetas: np.ndarray, dst: np.ndarray) -> np.ndarray: + """ + Initializes 4x4 matrices of 2-qubit gates defined in the paper. + + Args: + thetas: depth x 4 matrix of gate parameters for every layer, where + "depth" is the number of layers. + dst: destination array of size depth x 4 x 4 that will receive gate + matrices of each layer. + + Returns: + Returns the "dst" array. + """ + n = thetas.shape[0] + tmp = np.full((4, 2, 2), fill_value=0, dtype=np.complex128) + for k in range(n): + th = thetas[k] + a = make_rz(th[0], out=tmp[0]) + b = make_ry(th[1], out=tmp[1]) + c = make_rz(th[2], out=tmp[2]) + np.dot(np.dot(a, b, out=tmp[3]), c, out=dst[k]) + return dst + + +def init_layer1q_deriv_matrices(thetas: np.ndarray, dst: np.ndarray) -> np.ndarray: + """ + Initializes 4x4 derivative matrices of 2-qubit gates defined in the paper. + + Args: + thetas: depth x 4 matrix of gate parameters for every layer, where + "depth" is the number of layers. + dst: destination array of size depth x 4 x 4 x 4 that will receive gate + derivative matrices of each layer; there are 4 parameters per gate, + hence, 4 derivative matrices per layer. + + Returns: + Returns the "dst" array. + """ + n = thetas.shape[0] + y = np.asarray([[0, -0.5], [0.5, 0]], dtype=np.complex128) + z = np.asarray([[-0.5j, 0], [0, 0.5j]], dtype=np.complex128) + tmp = np.full((5, 2, 2), fill_value=0, dtype=np.complex128) + for k in range(n): + th = thetas[k] + a = make_rz(th[0], out=tmp[0]) + b = make_ry(th[1], out=tmp[1]) + c = make_rz(th[2], out=tmp[2]) + + za = np.dot(z, a, out=tmp[3]) + np.dot(np.dot(za, b, out=tmp[4]), c, out=dst[k, 0]) + yb = np.dot(y, b, out=tmp[3]) + np.dot(a, np.dot(yb, c, out=tmp[4]), out=dst[k, 1]) + zc = np.dot(z, c, out=tmp[3]) + np.dot(a, np.dot(b, zc, out=tmp[4]), out=dst[k, 2]) + return dst + + +def init_layer2q_matrices(thetas: np.ndarray, dst: np.ndarray) -> np.ndarray: + """ + Initializes 4x4 matrices of 2-qubit gates defined in the paper. + + Args: + thetas: depth x 4 matrix of gate parameters for every layer, where + "depth" is the number of layers. + dst: destination array of size depth x 4 x 4 that will receive gate + matrices of each layer. + + Returns: + Returns the "dst" array. + """ + depth = thetas.shape[0] + for k in range(depth): + th = thetas[k] + cs0 = np.cos(0.5 * th[0]).item() + sn0 = np.sin(0.5 * th[0]).item() + ep1 = np.exp(0.5j * th[1]).item() + en1 = np.exp(-0.5j * th[1]).item() + cs2 = np.cos(0.5 * th[2]).item() + sn2 = np.sin(0.5 * th[2]).item() + cs3 = np.cos(0.5 * th[3]).item() + sn3 = np.sin(0.5 * th[3]).item() + ep1cs0 = ep1 * cs0 + ep1sn0 = ep1 * sn0 + en1cs0 = en1 * cs0 + en1sn0 = en1 * sn0 + sn2cs3 = sn2 * cs3 + sn2sn3 = sn2 * sn3 + cs2cs3 = cs2 * cs3 + sn3cs2j = 1j * sn3 * cs2 + sn2sn3j = 1j * sn2sn3 + + flat_dst = dst[k].ravel() + flat_dst[:] = [ + -(sn2sn3j - cs2cs3) * en1cs0, + -(sn2cs3 + sn3cs2j) * en1cs0, + (sn2cs3 + sn3cs2j) * en1sn0, + (sn2sn3j - cs2cs3) * en1sn0, + (sn2cs3 - sn3cs2j) * en1cs0, + (sn2sn3j + cs2cs3) * en1cs0, + -(sn2sn3j + cs2cs3) * en1sn0, + (-sn2cs3 + sn3cs2j) * en1sn0, + (-sn2sn3j + cs2cs3) * ep1sn0, + -(sn2cs3 + sn3cs2j) * ep1sn0, + -(sn2cs3 + sn3cs2j) * ep1cs0, + (-sn2sn3j + cs2cs3) * ep1cs0, + (sn2cs3 - sn3cs2j) * ep1sn0, + (sn2sn3j + cs2cs3) * ep1sn0, + (sn2sn3j + cs2cs3) * ep1cs0, + (sn2cs3 - sn3cs2j) * ep1cs0, + ] + return dst + + +def init_layer2q_deriv_matrices(thetas: np.ndarray, dst: np.ndarray) -> np.ndarray: + """ + Initializes 4 x 4 derivative matrices of 2-qubit gates defined in the paper. + + Args: + thetas: depth x 4 matrix of gate parameters for every layer, where + "depth" is the number of layers. + dst: destination array of size depth x 4 x 4 x 4 that will receive gate + derivative matrices of each layer; there are 4 parameters per gate, + hence, 4 derivative matrices per layer. + + Returns: + Returns the "dst" array. + """ + depth = thetas.shape[0] + for k in range(depth): + th = thetas[k] + cs0 = np.cos(0.5 * th[0]).item() + sn0 = np.sin(0.5 * th[0]).item() + ep1 = np.exp(0.5j * th[1]).item() * 0.5 + en1 = np.exp(-0.5j * th[1]).item() * 0.5 + cs2 = np.cos(0.5 * th[2]).item() + sn2 = np.sin(0.5 * th[2]).item() + cs3 = np.cos(0.5 * th[3]).item() + sn3 = np.sin(0.5 * th[3]).item() + ep1cs0 = ep1 * cs0 + ep1sn0 = ep1 * sn0 + en1cs0 = en1 * cs0 + en1sn0 = en1 * sn0 + sn2cs3 = sn2 * cs3 + sn2sn3 = sn2 * sn3 + sn3cs2 = sn3 * cs2 + cs2cs3 = cs2 * cs3 + sn2cs3j = 1j * sn2cs3 + sn2sn3j = 1j * sn2sn3 + sn3cs2j = 1j * sn3cs2 + cs2cs3j = 1j * cs2cs3 + + flat_dst = dst[k, 0].ravel() + flat_dst[:] = [ + (sn2sn3j - cs2cs3) * en1sn0, + (sn2cs3 + sn3cs2j) * en1sn0, + (sn2cs3 + sn3cs2j) * en1cs0, + (sn2sn3j - cs2cs3) * en1cs0, + (-sn2cs3 + sn3cs2j) * en1sn0, + -(sn2sn3j + cs2cs3) * en1sn0, + -(sn2sn3j + cs2cs3) * en1cs0, + (-sn2cs3 + sn3cs2j) * en1cs0, + (-sn2sn3j + cs2cs3) * ep1cs0, + -(sn2cs3 + sn3cs2j) * ep1cs0, + (sn2cs3 + sn3cs2j) * ep1sn0, + (sn2sn3j - cs2cs3) * ep1sn0, + (sn2cs3 - sn3cs2j) * ep1cs0, + (sn2sn3j + cs2cs3) * ep1cs0, + -(sn2sn3j + cs2cs3) * ep1sn0, + (-sn2cs3 + sn3cs2j) * ep1sn0, + ] + + flat_dst = dst[k, 1].ravel() + flat_dst[:] = [ + -(sn2sn3 + cs2cs3j) * en1cs0, + (sn2cs3j - sn3cs2) * en1cs0, + -(sn2cs3j - sn3cs2) * en1sn0, + (sn2sn3 + cs2cs3j) * en1sn0, + -(sn2cs3j + sn3cs2) * en1cs0, + (sn2sn3 - cs2cs3j) * en1cs0, + (-sn2sn3 + cs2cs3j) * en1sn0, + (sn2cs3j + sn3cs2) * en1sn0, + (sn2sn3 + cs2cs3j) * ep1sn0, + (-sn2cs3j + sn3cs2) * ep1sn0, + (-sn2cs3j + sn3cs2) * ep1cs0, + (sn2sn3 + cs2cs3j) * ep1cs0, + (sn2cs3j + sn3cs2) * ep1sn0, + (-sn2sn3 + cs2cs3j) * ep1sn0, + (-sn2sn3 + cs2cs3j) * ep1cs0, + (sn2cs3j + sn3cs2) * ep1cs0, + ] + + flat_dst = dst[k, 2].ravel() + flat_dst[:] = [ + -(sn2cs3 + sn3cs2j) * en1cs0, + (sn2sn3j - cs2cs3) * en1cs0, + -(sn2sn3j - cs2cs3) * en1sn0, + (sn2cs3 + sn3cs2j) * en1sn0, + (sn2sn3j + cs2cs3) * en1cs0, + (-sn2cs3 + sn3cs2j) * en1cs0, + (sn2cs3 - sn3cs2j) * en1sn0, + -(sn2sn3j + cs2cs3) * en1sn0, + -(sn2cs3 + sn3cs2j) * ep1sn0, + (sn2sn3j - cs2cs3) * ep1sn0, + (sn2sn3j - cs2cs3) * ep1cs0, + -(sn2cs3 + sn3cs2j) * ep1cs0, + (sn2sn3j + cs2cs3) * ep1sn0, + (-sn2cs3 + sn3cs2j) * ep1sn0, + (-sn2cs3 + sn3cs2j) * ep1cs0, + (sn2sn3j + cs2cs3) * ep1cs0, + ] + + flat_dst = dst[k, 3].ravel() + flat_dst[:] = [ + -(sn2cs3j + sn3cs2) * en1cs0, + (sn2sn3 - cs2cs3j) * en1cs0, + (-sn2sn3 + cs2cs3j) * en1sn0, + (sn2cs3j + sn3cs2) * en1sn0, + -(sn2sn3 + cs2cs3j) * en1cs0, + (sn2cs3j - sn3cs2) * en1cs0, + -(sn2cs3j - sn3cs2) * en1sn0, + (sn2sn3 + cs2cs3j) * en1sn0, + -(sn2cs3j + sn3cs2) * ep1sn0, + (sn2sn3 - cs2cs3j) * ep1sn0, + (sn2sn3 - cs2cs3j) * ep1cs0, + -(sn2cs3j + sn3cs2) * ep1cs0, + -(sn2sn3 + cs2cs3j) * ep1sn0, + (sn2cs3j - sn3cs2) * ep1sn0, + (sn2cs3j - sn3cs2) * ep1cs0, + -(sn2sn3 + cs2cs3j) * ep1cs0, + ] + return dst diff --git a/qiskit/synthesis/unitary/aqc/fast_gradient/pmatrix.py b/qiskit/synthesis/unitary/aqc/fast_gradient/pmatrix.py new file mode 100644 index 000000000000..b68064323698 --- /dev/null +++ b/qiskit/synthesis/unitary/aqc/fast_gradient/pmatrix.py @@ -0,0 +1,312 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2023. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +""" +Matrix designed for fast multiplication by permutation and block-diagonal ones. +""" + +import numpy as np +from .layer import Layer1Q, Layer2Q + + +class PMatrix: + """ + Wrapper around a matrix that enables fast multiplication by permutation + matrices and block-diagonal ones. + """ + + def __init__(self, num_qubits: int): + """ + Initializes the internal structures of this object but does not set + the matrix yet. + + Args: + num_qubits: number of qubits. + """ + dim = 2**num_qubits + self._mat = np.empty(0) + self._dim = dim + self._temp_g2x2 = np.zeros((2, 2), dtype=np.complex128) + self._temp_g4x4 = np.zeros((4, 4), dtype=np.complex128) + self._temp_2x2 = self._temp_g2x2.copy() + self._temp_4x4 = self._temp_g4x4.copy() + self._identity_perm = np.arange(dim, dtype=np.int64) + self._left_perm = self._identity_perm.copy() + self._right_perm = self._identity_perm.copy() + self._temp_perm = self._identity_perm.copy() + self._temp_slice_dim_x_2 = np.zeros((dim, 2), dtype=np.complex128) + self._temp_slice_dim_x_4 = np.zeros((dim, 4), dtype=np.complex128) + self._idx_mat = self._init_index_matrix(dim) + self._temp_block_diag = np.zeros(self._idx_mat.shape, dtype=np.complex128) + + def set_matrix(self, mat: np.ndarray): + """ + Copies specified matrix to internal storage. Once the matrix + is set, the object is ready for use. + + **Note**, the matrix will be copied, mind the size issues. + + Args: + mat: matrix we want to multiply on the left and on the right by + layer matrices. + """ + if self._mat.size == 0: + self._mat = mat.copy() + else: + np.copyto(self._mat, mat) + + @staticmethod + def _init_index_matrix(dim: int) -> np.ndarray: + """ + Fast multiplication can be implemented by picking up a subset of + entries in a sparse matrix. + + Args: + dim: problem dimensionality. + + Returns: + 2d-array of indices for the fast multiplication. + """ + all_idx = np.arange(dim * dim, dtype=np.int64).reshape(dim, dim) + idx = np.full((dim // 4, 4 * 4), fill_value=0, dtype=np.int64) + b = np.full((4, 4), fill_value=0, dtype=np.int64) + for i in range(0, dim, 4): + b[:, :] = all_idx[i : i + 4, i : i + 4] + idx[i // 4, :] = b.T.ravel() + return idx + + def mul_right_q1(self, layer: Layer1Q, temp_mat: np.ndarray, dagger: bool): + """ + Multiplies ``NxN`` matrix, wrapped by this object, by a 1-qubit layer + matrix of the right, where ``N`` is the actual size of matrices involved, + ``N = 2^{num. of qubits}``. + + Args: + layer: 1-qubit layer, i.e. the layer with just one non-trivial + 1-qubit gate and other gates are just identity operators. + temp_mat: a temporary NxN matrix used as a workspace. + dagger: if true, the right-hand side matrix will be taken as + conjugate transposed. + """ + + gmat, perm, inv_perm = layer.get_attr() + mat = self._mat + dim = perm.size + + # temp_mat <-- mat[:, right_perm[perm]] = mat[:, right_perm][:, perm]: + np.take(mat, np.take(self._right_perm, perm, out=self._temp_perm), axis=1, out=temp_mat) + + # mat <-- mat[:, right_perm][:, perm] @ kron(I(N/4), gmat)^dagger, where + # conjugate-transposition might be or might be not applied: + gmat_right = np.conj(gmat, out=self._temp_g2x2).T if dagger else gmat + for i in range(0, dim, 2): + mat[:, i : i + 2] = np.dot( + temp_mat[:, i : i + 2], gmat_right, out=self._temp_slice_dim_x_2 + ) + + # Update right permutation: + self._right_perm[:] = inv_perm + + def mul_right_q2(self, layer: Layer2Q, temp_mat: np.ndarray, dagger: bool = True): + """ + Multiplies ``NxN`` matrix, wrapped by this object, by a 2-qubit layer + matrix on the right, where ``N`` is the actual size of matrices involved, + ``N = 2^{num. of qubits}``. + + Args: + layer: 2-qubit layer, i.e. the layer with just one non-trivial + 2-qubit gate and other gates are just identity operators. + temp_mat: a temporary NxN matrix used as a workspace. + dagger: if true, the right-hand side matrix will be taken as + conjugate transposed. + """ + + gmat, perm, inv_perm = layer.get_attr() + mat = self._mat + dim = perm.size + + # temp_mat <-- mat[:, right_perm[perm]] = mat[:, right_perm][:, perm]: + np.take(mat, np.take(self._right_perm, perm, out=self._temp_perm), axis=1, out=temp_mat) + + # mat <-- mat[:, right_perm][:, perm] @ kron(I(N/4), gmat)^dagger, where + # conjugate-transposition might be or might be not applied: + gmat_right = np.conj(gmat, out=self._temp_g4x4).T if dagger else gmat + for i in range(0, dim, 4): + mat[:, i : i + 4] = np.dot( + temp_mat[:, i : i + 4], gmat_right, out=self._temp_slice_dim_x_4 + ) + + # Update right permutation: + self._right_perm[:] = inv_perm + + def mul_left_q1(self, layer: Layer1Q, temp_mat: np.ndarray): + """ + Multiplies ``NxN`` matrix, wrapped by this object, by a 1-qubit layer + matrix of the left, where ``dim`` is the actual size of matrices involved, + ``dim = 2^{num. of qubits}``. + + Args: + layer: 1-qubit layer, i.e. the layer with just one non-trivial + 1-qubit gate and other gates are just identity operators. + temp_mat: a temporary NxN matrix used as a workspace. + """ + mat = self._mat + gmat, perm, inv_perm = layer.get_attr() + dim = perm.size + + # temp_mat <-- mat[left_perm[perm]] = mat[left_perm][perm]: + np.take(mat, np.take(self._left_perm, perm, out=self._temp_perm), axis=0, out=temp_mat) + + # mat <-- kron(I(dim/4), gmat) @ mat[left_perm][perm]: + if dim > 512: + # Faster for large matrices. + for i in range(0, dim, 2): + np.dot(gmat, temp_mat[i : i + 2, :], out=mat[i : i + 2, :]) + else: + # Faster for small matrices. + half = dim // 2 + np.copyto( + mat.reshape((2, half, dim)), np.swapaxes(temp_mat.reshape((half, 2, dim)), 0, 1) + ) + np.dot(gmat, mat.reshape(2, -1), out=temp_mat.reshape(2, -1)) + np.copyto( + mat.reshape((half, 2, dim)), np.swapaxes(temp_mat.reshape((2, half, dim)), 0, 1) + ) + + # Update left permutation: + self._left_perm[:] = inv_perm + + def mul_left_q2(self, layer: Layer2Q, temp_mat: np.ndarray): + """ + Multiplies ``NxN`` matrix, wrapped by this object, by a 2-qubit layer + matrix on the left, where ``dim`` is the actual size of matrices involved, + ``dim = 2^{num. of qubits}``. + + Args: + layer: 2-qubit layer, i.e. the layer with just one non-trivial + 2-qubit gate and other gates are just identity operators. + temp_mat: a temporary NxN matrix used as a workspace. + """ + mat = self._mat + gmat, perm, inv_perm = layer.get_attr() + dim = perm.size + + # temp_mat <-- mat[left_perm[perm]] = mat[left_perm][perm]: + np.take(mat, np.take(self._left_perm, perm, out=self._temp_perm), axis=0, out=temp_mat) + + # mat <-- kron(I(dim/4), gmat) @ mat[left_perm][perm]: + if dim > 512: + # Faster for large matrices. + for i in range(0, dim, 4): + np.dot(gmat, temp_mat[i : i + 4, :], out=mat[i : i + 4, :]) + else: + # Faster for small matrices. + half = dim // 4 + np.copyto( + mat.reshape((4, half, dim)), np.swapaxes(temp_mat.reshape((half, 4, dim)), 0, 1) + ) + np.dot(gmat, mat.reshape(4, -1), out=temp_mat.reshape(4, -1)) + np.copyto( + mat.reshape((half, 4, dim)), np.swapaxes(temp_mat.reshape((4, half, dim)), 0, 1) + ) + + # Update left permutation: + self._left_perm[:] = inv_perm + + def product_q1(self, layer: Layer1Q, tmp1: np.ndarray, tmp2: np.ndarray) -> np.complex128: + """ + Computes and returns: ``Trace(mat @ C) = Trace(mat @ P^T @ gmat @ P) = + Trace((P @ mat @ P^T) @ gmat) = Trace(C @ (P @ mat @ P^T)) = + vec(gmat^T)^T @ vec(P @ mat @ P^T)``, where mat is ``NxN`` matrix wrapped + by this object, ``C`` is matrix representation of the layer ``L``, and gmat + is 2x2 matrix of underlying 1-qubit gate. + + **Note**: matrix of this class must be finalized beforehand. + + Args: + layer: 1-qubit layer. + tmp1: temporary, external matrix used as a workspace. + tmp2: temporary, external matrix used as a workspace. + + Returns: + trace of the matrix product. + """ + mat = self._mat + gmat, perm, _ = layer.get_attr() + + # tmp2 = P @ mat @ P^T: + np.take(np.take(mat, perm, axis=0, out=tmp1), perm, axis=1, out=tmp2) + + # matrix dot product = Tr(transposed(kron(I(dim/4), gmat)), (P @ mat @ P^T)): + gmat_t, tmp3 = self._temp_g2x2, self._temp_2x2 + np.copyto(gmat_t, gmat.T) + _sum = 0.0 + for i in range(0, mat.shape[0], 2): + tmp3[:, :] = tmp2[i : i + 2, i : i + 2] + _sum += np.dot(gmat_t.ravel(), tmp3.ravel()) + + return np.complex128(_sum) + + def product_q2(self, layer: Layer2Q, tmp1: np.ndarray, tmp2: np.ndarray) -> np.complex128: + """ + Computes and returns: ``Trace(mat @ C) = Trace(mat @ P^T @ gmat @ P) = + Trace((P @ mat @ P^T) @ gmat) = Trace(C @ (P @ mat @ P^T)) = + vec(gmat^T)^T @ vec(P @ mat @ P^T)``, where mat is ``NxN`` matrix wrapped + by this object, ``C`` is matrix representation of the layer ``L``, and gmat + is 4x4 matrix of underlying 2-qubit gate. + + **Note**: matrix of this class must be finalized beforehand. + + Args: + layer: 2-qubit layer. + tmp1: temporary, external matrix used as a workspace. + tmp2: temporary, external matrix used as a workspace. + + Returns: + trace of the matrix product. + """ + mat = self._mat + gmat, perm, _ = layer.get_attr() + + # Compute the matrix dot product: + # Tr(transposed(kron(I(dim/4), gmat)), (P @ mat @ P^T)): + + # The fastest version so far, but requires two NxN temp. matrices. + # tmp2 = P @ mat @ P^T: + np.take(np.take(mat, perm, axis=0, out=tmp1), perm, axis=1, out=tmp2) + + bldia = self._temp_block_diag + np.take(tmp2.ravel(), self._idx_mat.ravel(), axis=0, out=bldia.ravel()) + bldia *= gmat.reshape(-1, gmat.size) + return np.complex128(np.sum(bldia)) + + def finalize(self, temp_mat: np.ndarray) -> np.ndarray: + """ + Applies the left (row) and right (column) permutations to the matrix. + at the end of computation process. + + Args: + temp_mat: temporary, external matrix. + + Returns: + finalized matrix with all transformations applied. + """ + mat = self._mat + + # mat <-- mat[left_perm][:, right_perm] = P_left @ mat @ transposed(P_right) + np.take(mat, self._left_perm, axis=0, out=temp_mat) + np.take(temp_mat, self._right_perm, axis=1, out=mat) + + # Set both permutations to identity once they have been applied. + self._left_perm[:] = self._identity_perm + self._right_perm[:] = self._identity_perm + return self._mat diff --git a/qiskit/transpiler/synthesis/graysynth.py b/qiskit/transpiler/synthesis/graysynth.py index 7f0ab1c08ba0..495bdf935746 100644 --- a/qiskit/transpiler/synthesis/graysynth.py +++ b/qiskit/transpiler/synthesis/graysynth.py @@ -23,8 +23,14 @@ # TODO: Deprecate in 0.24.0 and remove in 0.26.0 from qiskit.synthesis.linear.cnot_synth import * from qiskit.synthesis.linear_phase.cnot_phase_synth import * +from qiskit.utils.deprecation import deprecate_func +@deprecate_func( + since="0.45.0", + additional_msg="Instead, use the function ``synth_cnot_count_full_pmh`` from the module" + + "``qiskit.synthesis.linear``.", +) def cnot_synth(state, section_size=2): """ Synthesize linear reversible circuits for all-to-all architecture @@ -55,6 +61,11 @@ def cnot_synth(state, section_size=2): return synth_cnot_count_full_pmh(state, section_size=section_size) +@deprecate_func( + since="0.45.0", + additional_msg="Instead, use the function ``synth_cnot_phase_aam`` from the module" + + "``qiskit.synthesis.linear_phase``.", +) def graysynth(cnots, angles, section_size=2): """This function is an implementation of the GraySynth algorithm of Amy, Azimadeh and Mosca. diff --git a/test/python/synthesis/aqc/__init__.py b/test/python/synthesis/aqc/__init__.py new file mode 100644 index 000000000000..26f7536d3514 --- /dev/null +++ b/test/python/synthesis/aqc/__init__.py @@ -0,0 +1,11 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2023. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. diff --git a/test/python/synthesis/aqc/fast_gradient/__init__.py b/test/python/synthesis/aqc/fast_gradient/__init__.py new file mode 100644 index 000000000000..fdb172d367f0 --- /dev/null +++ b/test/python/synthesis/aqc/fast_gradient/__init__.py @@ -0,0 +1,11 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2022. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. diff --git a/test/python/synthesis/aqc/fast_gradient/test_cmp_gradients.py b/test/python/synthesis/aqc/fast_gradient/test_cmp_gradients.py new file mode 100644 index 000000000000..11c7fb4a7484 --- /dev/null +++ b/test/python/synthesis/aqc/fast_gradient/test_cmp_gradients.py @@ -0,0 +1,110 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2023. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +""" +Tests equivalence of the default and fast gradient computation routines. +""" + +import unittest +from typing import Tuple +from time import perf_counter +from test.python.synthesis.aqc.fast_gradient.utils_for_testing import rand_circuit, rand_su_mat +import numpy as np +from qiskit.synthesis.unitary.aqc.fast_gradient.fast_gradient import FastCNOTUnitObjective +from qiskit.synthesis.unitary.aqc.cnot_unit_objective import DefaultCNOTUnitObjective +from qiskit.test import QiskitTestCase + + +class TestCompareGradientImpls(QiskitTestCase): + """ + Tests equivalence of the default and fast gradient implementations. + """ + + max_num_qubits = 3 # maximum number of qubits in tests + max_depth = 10 # maximum circuit depth in tests + + def setUp(self): + super().setUp() + np.random.seed(0x0696969) + + def _compare( + self, num_qubits: int, depth: int + ) -> Tuple[int, int, float, float, float, float, float]: + """ + Calculates gradient and objective function value for the original + and the fast implementations, and compares the outputs from both. + Returns relative errors. Also, accumulates performance metrics. + """ + + cnots = rand_circuit(num_qubits=num_qubits, depth=depth) + depth = cnots.shape[1] # might be less than initial depth + u_mat = rand_su_mat(2**num_qubits) + dflt_obj = DefaultCNOTUnitObjective(num_qubits=num_qubits, cnots=cnots) + fast_obj = FastCNOTUnitObjective(num_qubits=num_qubits, cnots=cnots) + thetas = np.random.rand(4 * depth + 3 * num_qubits) * 2 * np.pi + thetas = thetas.astype(np.float64) + + dflt_obj.target_matrix = u_mat + fast_obj.target_matrix = u_mat + + # Compute fast gradient. + start = perf_counter() + f1 = fast_obj.objective(param_values=thetas) + g1 = fast_obj.gradient(param_values=thetas) + t1 = perf_counter() - start + + # Compute default gradient. + start = perf_counter() + f2 = dflt_obj.objective(param_values=thetas) + g2 = dflt_obj.gradient(param_values=thetas) + t2 = perf_counter() - start + + fobj_rel_err = abs(f1 - f2) / f2 + grad_rel_err = np.linalg.norm(g1 - g2) / np.linalg.norm(g2) + speedup = t2 / max(t1, np.finfo(np.float64).eps ** 2) + + tol = float(np.sqrt(np.finfo(np.float64).eps)) + self.assertLess(fobj_rel_err, tol) + self.assertLess(grad_rel_err, tol) + self.assertTrue(np.allclose(g1, g2, atol=tol, rtol=tol)) + return int(num_qubits), int(depth), fobj_rel_err, grad_rel_err, speedup, t1, t2 + + def test_cmp_gradients(self): + """ + Tests equivalence of gradients. + """ + + # Configurations of the number of qubits and depths we want to try. + configs = [ + (n, depth) + for n in range(2, self.max_num_qubits + 1) + for depth in np.sort( + np.random.permutation( + np.arange(3 if n <= 3 else 7, 9 if n <= 3 else self.max_depth) + )[0:10] + ) + ] + + results = [] + + # Run the tests sequentially. + for nqubits, depth in configs: + results.append(self._compare(nqubits, depth)) + + tol = float(np.sqrt(np.finfo(float).eps)) + for _, _, ferr, gerr, _, _, _ in results: + self.assertTrue(ferr < tol) + self.assertTrue(gerr < tol) + + +if __name__ == "__main__": + unittest.main() diff --git a/test/python/synthesis/aqc/fast_gradient/test_layer1q.py b/test/python/synthesis/aqc/fast_gradient/test_layer1q.py new file mode 100644 index 000000000000..e83150261dc5 --- /dev/null +++ b/test/python/synthesis/aqc/fast_gradient/test_layer1q.py @@ -0,0 +1,141 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2023. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +""" +Tests for Layer1Q implementation. +""" + +import unittest +from random import randint +import test.python.synthesis.aqc.fast_gradient.utils_for_testing as tut +import numpy as np +import qiskit.synthesis.unitary.aqc.fast_gradient.layer as lr +from qiskit.synthesis.unitary.aqc.fast_gradient.pmatrix import PMatrix +from qiskit.test import QiskitTestCase + + +class TestLayer1q(QiskitTestCase): + """ + Tests for Layer1Q class. + """ + + max_num_qubits = 5 # maximum number of qubits in tests + num_repeats = 50 # number of repetitions in tests + + def setUp(self): + super().setUp() + np.random.seed(0x0696969) + + def test_layer1q_matrix(self): + """ + Tests: (1) the correctness of Layer2Q matrix construction; + (2) matrix multiplication interleaved with permutations. + """ + + mat_kind = "complex" + eps = 100.0 * np.finfo(float).eps + max_rel_err = 0.0 + for n in range(2, self.max_num_qubits + 1): + + dim = 2**n + iden = tut.eye_int(n) + for k in range(n): + m_mat = tut.rand_matrix(dim=dim, kind=mat_kind) + t_mat, g_mat = tut.make_test_matrices2x2(n=n, k=k, kind=mat_kind) + lmat = lr.Layer1Q(num_qubits=n, k=k, g2x2=g_mat) + g2, perm, inv_perm = lmat.get_attr() + self.assertTrue(m_mat.dtype == t_mat.dtype == g_mat.dtype == g2.dtype) + self.assertTrue(np.all(g_mat == g2)) + self.assertTrue(np.all(iden[perm].T == iden[inv_perm])) + + g_mat = np.kron(tut.eye_int(n - 1), g_mat) + + # T == P^t @ G @ P. + err = tut.relative_error(t_mat, iden[perm].T @ g_mat @ iden[perm]) + self.assertLess(err, eps, "err = {:0.16f}".format(err)) + max_rel_err = max(max_rel_err, err) + + # Multiplication by permutation matrix of the left can be + # replaced by row permutations. + tm = t_mat @ m_mat + err1 = tut.relative_error(iden[perm].T @ g_mat @ m_mat[perm], tm) + err2 = tut.relative_error((g_mat @ m_mat[perm])[inv_perm], tm) + + # Multiplication by permutation matrix of the right can be + # replaced by column permutations. + mt = m_mat @ t_mat + err3 = tut.relative_error(m_mat @ iden[perm].T @ g_mat @ iden[perm], mt) + err4 = tut.relative_error((m_mat[:, perm] @ g_mat)[:, inv_perm], mt) + + self.assertTrue( + err1 < eps and err2 < eps and err3 < eps and err4 < eps, + "err1 = {:f}, err2 = {:f}, " + "err3 = {:f}, err4 = {:f}".format(err1, err2, err3, err4), + ) + max_rel_err = max(max_rel_err, err1, err2, err3, err4) + + def test_pmatrix_class(self): + """ + Test the class PMatrix. + """ + + _eps = 100.0 * np.finfo(float).eps + mat_kind = "complex" + max_rel_err = 0.0 + for n in range(2, self.max_num_qubits + 1): + + dim = 2**n + tmp1 = np.ndarray((dim, dim), dtype=np.complex128) + tmp2 = tmp1.copy() + for _ in range(self.num_repeats): + k0 = randint(0, n - 1) + k1 = randint(0, n - 1) + k2 = randint(0, n - 1) + k3 = randint(0, n - 1) + k4 = randint(0, n - 1) + + t0, g0 = tut.make_test_matrices2x2(n=n, k=k0, kind=mat_kind) + t1, g1 = tut.make_test_matrices2x2(n=n, k=k1, kind=mat_kind) + t2, g2 = tut.make_test_matrices2x2(n=n, k=k2, kind=mat_kind) + t3, g3 = tut.make_test_matrices2x2(n=n, k=k3, kind=mat_kind) + t4, g4 = tut.make_test_matrices2x2(n=n, k=k4, kind=mat_kind) + + c0 = lr.Layer1Q(num_qubits=n, k=k0, g2x2=g0) + c1 = lr.Layer1Q(num_qubits=n, k=k1, g2x2=g1) + c2 = lr.Layer1Q(num_qubits=n, k=k2, g2x2=g2) + c3 = lr.Layer1Q(num_qubits=n, k=k3, g2x2=g3) + c4 = lr.Layer1Q(num_qubits=n, k=k4, g2x2=g4) + + m_mat = tut.rand_matrix(dim=dim, kind=mat_kind) + ttmtt = t0 @ t1 @ m_mat @ np.conj(t2).T @ np.conj(t3).T + + pmat = PMatrix(n) + pmat.set_matrix(m_mat) + pmat.mul_left_q1(layer=c1, temp_mat=tmp1) + pmat.mul_left_q1(layer=c0, temp_mat=tmp1) + pmat.mul_right_q1(layer=c2, temp_mat=tmp1, dagger=True) + pmat.mul_right_q1(layer=c3, temp_mat=tmp1, dagger=True) + alt_ttmtt = pmat.finalize(temp_mat=tmp1) + + err1 = tut.relative_error(alt_ttmtt, ttmtt) + self.assertLess(err1, _eps, "relative error: {:f}".format(err1)) + + prod = np.complex128(np.trace(ttmtt @ t4)) + alt_prod = pmat.product_q1(layer=c4, tmp1=tmp1, tmp2=tmp2) + err2 = abs(alt_prod - prod) / abs(prod) + self.assertLess(err2, _eps, "relative error: {:f}".format(err2)) + + max_rel_err = max(max_rel_err, err1, err2) + + +if __name__ == "__main__": + unittest.main() diff --git a/test/python/synthesis/aqc/fast_gradient/test_layer2q.py b/test/python/synthesis/aqc/fast_gradient/test_layer2q.py new file mode 100644 index 000000000000..b916d72bd6de --- /dev/null +++ b/test/python/synthesis/aqc/fast_gradient/test_layer2q.py @@ -0,0 +1,149 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2023. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +""" +Tests for Layer2Q implementation. +""" + +import unittest +from random import randint +import test.python.synthesis.aqc.fast_gradient.utils_for_testing as tut +import numpy as np +import qiskit.synthesis.unitary.aqc.fast_gradient.layer as lr +from qiskit.synthesis.unitary.aqc.fast_gradient.pmatrix import PMatrix +from qiskit.test import QiskitTestCase + + +class TestLayer2q(QiskitTestCase): + """ + Tests for Layer2Q class. + """ + + max_num_qubits = 5 # maximum number of qubits in tests + num_repeats = 50 # number of repetitions in tests + + def setUp(self): + super().setUp() + np.random.seed(0x0696969) + + def test_layer2q_matrix(self): + """ + Tests: (1) the correctness of Layer2Q matrix construction; + (2) matrix multiplication interleaved with permutations. + """ + + mat_kind = "complex" + _eps = 100.0 * np.finfo(float).eps + max_rel_err = 0.0 + for n in range(2, self.max_num_qubits + 1): + + dim = 2**n + iden = tut.eye_int(n) + for j in range(n): + for k in range(n): + if j == k: + continue + m_mat = tut.rand_matrix(dim=dim, kind=mat_kind) + t_mat, g_mat = tut.make_test_matrices4x4(n=n, j=j, k=k, kind=mat_kind) + lmat = lr.Layer2Q(num_qubits=n, j=j, k=k, g4x4=g_mat) + g2, perm, inv_perm = lmat.get_attr() + self.assertTrue(m_mat.dtype == t_mat.dtype == g_mat.dtype == g2.dtype) + self.assertTrue(np.all(g_mat == g2)) + self.assertTrue(np.all(iden[perm].T == iden[inv_perm])) + + g_mat = np.kron(tut.eye_int(n - 2), g_mat) + + # T == P^t @ G @ P. + err = tut.relative_error(t_mat, iden[perm].T @ g_mat @ iden[perm]) + self.assertLess(err, _eps, "err = {:0.16f}".format(err)) + max_rel_err = max(max_rel_err, err) + + # Multiplication by permutation matrix of the left can be + # replaced by row permutations. + tm = t_mat @ m_mat + err1 = tut.relative_error(iden[perm].T @ g_mat @ m_mat[perm], tm) + err2 = tut.relative_error((g_mat @ m_mat[perm])[inv_perm], tm) + + # Multiplication by permutation matrix of the right can be + # replaced by column permutations. + mt = m_mat @ t_mat + err3 = tut.relative_error(m_mat @ iden[perm].T @ g_mat @ iden[perm], mt) + err4 = tut.relative_error((m_mat[:, perm] @ g_mat)[:, inv_perm], mt) + + self.assertTrue( + err1 < _eps and err2 < _eps and err3 < _eps and err4 < _eps, + "err1 = {:f}, err2 = {:f}, " + "err3 = {:f}, err4 = {:f}".format(err1, err2, err3, err4), + ) + max_rel_err = max(max_rel_err, err1, err2, err3, err4) + + def test_pmatrix_class(self): + """ + Test the class PMatrix. + """ + + _eps = 100.0 * np.finfo(float).eps + mat_kind = "complex" + max_rel_err = 0.0 + for n in range(2, self.max_num_qubits + 1): + + dim = 2**n + tmp1 = np.ndarray((dim, dim), dtype=np.complex128) + tmp2 = tmp1.copy() + for _ in range(self.num_repeats): + j0 = randint(0, n - 1) + k0 = (j0 + randint(1, n - 1)) % n + j1 = randint(0, n - 1) + k1 = (j1 + randint(1, n - 1)) % n + j2 = randint(0, n - 1) + k2 = (j2 + randint(1, n - 1)) % n + j3 = randint(0, n - 1) + k3 = (j3 + randint(1, n - 1)) % n + j4 = randint(0, n - 1) + k4 = (j4 + randint(1, n - 1)) % n + + t0, g0 = tut.make_test_matrices4x4(n=n, j=j0, k=k0, kind=mat_kind) + t1, g1 = tut.make_test_matrices4x4(n=n, j=j1, k=k1, kind=mat_kind) + t2, g2 = tut.make_test_matrices4x4(n=n, j=j2, k=k2, kind=mat_kind) + t3, g3 = tut.make_test_matrices4x4(n=n, j=j3, k=k3, kind=mat_kind) + t4, g4 = tut.make_test_matrices4x4(n=n, j=j4, k=k4, kind=mat_kind) + + c0 = lr.Layer2Q(num_qubits=n, j=j0, k=k0, g4x4=g0) + c1 = lr.Layer2Q(num_qubits=n, j=j1, k=k1, g4x4=g1) + c2 = lr.Layer2Q(num_qubits=n, j=j2, k=k2, g4x4=g2) + c3 = lr.Layer2Q(num_qubits=n, j=j3, k=k3, g4x4=g3) + c4 = lr.Layer2Q(num_qubits=n, j=j4, k=k4, g4x4=g4) + + m_mat = tut.rand_matrix(dim=dim, kind=mat_kind) + ttmtt = t0 @ t1 @ m_mat @ np.conj(t2).T @ np.conj(t3).T + + pmat = PMatrix(n) + pmat.set_matrix(m_mat) + pmat.mul_left_q2(layer=c1, temp_mat=tmp1) + pmat.mul_left_q2(layer=c0, temp_mat=tmp1) + pmat.mul_right_q2(layer=c2, temp_mat=tmp1) + pmat.mul_right_q2(layer=c3, temp_mat=tmp1) + alt_ttmtt = pmat.finalize(temp_mat=tmp1) + + err1 = tut.relative_error(alt_ttmtt, ttmtt) + self.assertLess(err1, _eps, "relative error: {:f}".format(err1)) + + prod = np.complex128(np.trace(ttmtt @ t4)) + alt_prod = pmat.product_q2(layer=c4, tmp1=tmp1, tmp2=tmp2) + err2 = abs(alt_prod - prod) / abs(prod) + self.assertLess(err2, _eps, "relative error: {:f}".format(err2)) + + max_rel_err = max(max_rel_err, err1, err2) + + +if __name__ == "__main__": + unittest.main() diff --git a/test/python/synthesis/aqc/fast_gradient/test_utils.py b/test/python/synthesis/aqc/fast_gradient/test_utils.py new file mode 100644 index 000000000000..723af7d99a1c --- /dev/null +++ b/test/python/synthesis/aqc/fast_gradient/test_utils.py @@ -0,0 +1,176 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2023. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +""" +Tests for utility functions. +""" + +import unittest +import random +import test.python.synthesis.aqc.fast_gradient.utils_for_testing as tut +import numpy as np +import qiskit.synthesis.unitary.aqc.fast_gradient.fast_grad_utils as myu +from qiskit.test import QiskitTestCase +from qiskit.synthesis.unitary.aqc.elementary_operations import rx_matrix as _rx +from qiskit.synthesis.unitary.aqc.elementary_operations import ry_matrix as _ry +from qiskit.synthesis.unitary.aqc.elementary_operations import rz_matrix as _rz + + +class TestUtils(QiskitTestCase): + """ + Tests utility functions used by the fast gradient implementation. + """ + + max_nqubits_rev_bits = 7 # max. number of qubits in test_reverse_bits + max_nqubits_mk_unit = 5 # max. number of qubits in test_make_unit_vector + max_nqubits_swap = 7 # max. number of qubits in test_swap + max_nqubits_perms = 5 # max. number of qubits in test_permutations + num_repeats_gates2x2 = 100 # number of repetitions in test_gates2x2 + + def setUp(self): + super().setUp() + np.random.seed(0x0696969) + + def test_reverse_bits(self): + """ + Tests the function utils.reverse_bits(). + """ + + nbits = self.max_nqubits_rev_bits + for x in range(2**nbits): + y = myu.reverse_bits(x, nbits, enable=True) + self.assertTrue(x == myu.reverse_bits(y, nbits, enable=True)) + for i in range(nbits): + self.assertTrue(((x >> i) & 1) == ((y >> (nbits - i - 1)) & 1)) + + x = np.arange(2**nbits, dtype=np.int64) + y = myu.reverse_bits(x, nbits, enable=True) + self.assertTrue(np.all(myu.reverse_bits(y, nbits, enable=True) == x)) + + def test_make_unit_vector(self): + """ + Tests the function utils_for_testing.make_unit_vector(). + """ + + # Computes unit vector as Kronecker product of (1 0) or (0 1) sub-vectors. + def _make_unit_vector_simple(num: int, nbits: int) -> np.ndarray: + self.assertTrue(isinstance(num, int) and isinstance(nbits, int)) + self.assertTrue(0 <= num < 2**nbits) + vec = 1 + for i in range(nbits): + x = (num >> i) & 1 + vec = np.kron(vec, np.asarray([1 - x, x], dtype=np.int64)) + self.assertTrue(vec.dtype == np.int64 and vec.shape == (2**nbits,)) + return vec + + for n in range(1, self.max_nqubits_mk_unit + 1): + for k in range(2**n): + v1 = _make_unit_vector_simple(k, n) + v2 = tut.make_unit_vector(k, n) + self.assertTrue(np.all(v1 == v2)) + + def test_swap(self): + """ + Tests the function utils.swap_bits(). + """ + + nbits = self.max_nqubits_swap + for x in range(2**nbits): + for a in range(nbits): + for b in range(nbits): + y = myu.swap_bits(x, a, b) + self.assertTrue(x == myu.swap_bits(y, a, b)) + self.assertTrue(((x >> a) & 1) == ((y >> b) & 1)) + self.assertTrue(((x >> b) & 1) == ((y >> a) & 1)) + + def test_permutations(self): + """ + Tests various properties of permutations and their compositions. + """ + + def _perm2_mat(perm: np.ndarray) -> np.ndarray: + """Creates permutation matrix from a permutation.""" + return np.eye(perm.size, dtype=np.int64)[perm] + + def _permutation_and_matrix(size: int) -> (np.ndarray, np.ndarray, np.ndarray, np.ndarray): + """Creates direct and inverse permutations and their matrices.""" + perm = np.random.permutation(np.arange(size, dtype=np.int64)) + inv_perm = myu.inverse_permutation(perm) + p_mat, q_mat = _perm2_mat(perm), _perm2_mat(inv_perm) + return perm, inv_perm, p_mat, q_mat + + for n in range(1, self.max_nqubits_perms + 1): + + dim = 2**n + for _ in range(100): + perm1, inv_perm1, p1_mat, q1_mat = _permutation_and_matrix(dim) + perm2, inv_perm2, p2_mat, q2_mat = _permutation_and_matrix(dim) + + iden = tut.eye_int(n) + self.assertTrue(np.all(p1_mat == q1_mat.T) and np.all(p1_mat.T == q1_mat)) + self.assertTrue(np.all(p2_mat == q2_mat.T) and np.all(p2_mat.T == q2_mat)) + self.assertTrue(np.all(p1_mat @ q1_mat == iden) and np.all(q1_mat @ p1_mat == iden)) + self.assertTrue(np.all(p2_mat @ q2_mat == iden) and np.all(q2_mat @ p2_mat == iden)) + + self.assertTrue(np.all(perm2[perm1] == np.take(perm2, perm1))) + self.assertTrue(np.all(perm1[perm2] == np.take(perm1, perm2))) + + self.assertTrue(np.all(p1_mat @ p2_mat == _perm2_mat(perm2[perm1]))) + self.assertTrue(np.all(p2_mat @ p1_mat == _perm2_mat(perm1[perm2]))) + + self.assertTrue(np.all(p1_mat @ q2_mat == _perm2_mat(inv_perm2[perm1]))) + self.assertTrue(np.all(q2_mat @ p1_mat == _perm2_mat(perm1[inv_perm2]))) + + self.assertTrue(np.all(p2_mat @ q1_mat == _perm2_mat(inv_perm1[perm2]))) + self.assertTrue(np.all(q1_mat @ p2_mat == _perm2_mat(perm2[inv_perm1]))) + + _, inv_perm3, p3_mat, _ = _permutation_and_matrix(dim) + perm4, _, _, q4_mat = _permutation_and_matrix(dim) + + mat = tut.rand_matrix(dim=dim, kind="randint") + + pmq = p1_mat @ mat @ q2_mat + self.assertTrue(np.all(pmq == mat[perm1][:, perm2])) + self.assertTrue(np.all(pmq == np.take(np.take(mat, perm1, axis=0), perm2, axis=1))) + pmq = p1_mat @ mat @ p2_mat + self.assertTrue(np.all(pmq == mat[perm1][:, inv_perm2])) + self.assertTrue( + np.all(pmq == np.take(np.take(mat, perm1, axis=0), inv_perm2, axis=1)) + ) + + ppmpq = p1_mat @ p2_mat @ mat @ p3_mat @ q4_mat + self.assertTrue(np.all(ppmpq == mat[perm2][perm1][:, inv_perm3][:, perm4])) + self.assertTrue(np.all(ppmpq == mat[perm2[perm1]][:, inv_perm3[perm4]])) + self.assertTrue( + np.all( + ppmpq + == np.take(np.take(mat, perm2[perm1], axis=0), inv_perm3[perm4], axis=1) + ) + ) + + def test_gates2x2(self): + """ + Tests implementation of 2x2 gate matrices vs the ones that create and + return temporary array upon every call (but have more compact code). + """ + + tol = np.finfo(np.float64).eps * 2.0 + out = np.full((2, 2), fill_value=0, dtype=np.complex128) + for test in range(self.num_repeats_gates2x2): # pylint: disable=unused-variable + phi = random.random() * 2.0 * np.pi + self.assertTrue(np.allclose(myu.make_rx(phi, out=out), _rx(phi), atol=tol, rtol=tol)) + self.assertTrue(np.allclose(myu.make_ry(phi, out=out), _ry(phi), atol=tol, rtol=tol)) + self.assertTrue(np.allclose(myu.make_rz(phi, out=out), _rz(phi), atol=tol, rtol=tol)) + + +if __name__ == "__main__": + unittest.main() diff --git a/test/python/synthesis/aqc/fast_gradient/utils_for_testing.py b/test/python/synthesis/aqc/fast_gradient/utils_for_testing.py new file mode 100644 index 000000000000..85bb263d184a --- /dev/null +++ b/test/python/synthesis/aqc/fast_gradient/utils_for_testing.py @@ -0,0 +1,179 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2023. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +""" +Utility functions for debugging and testing. +""" + +from typing import Tuple +import numpy as np +from scipy.stats import unitary_group +import qiskit.synthesis.unitary.aqc.fast_gradient.fast_grad_utils as fgu + + +def relative_error(a_mat: np.ndarray, b_mat: np.ndarray) -> float: + """ + Computes relative residual between two matrices in Frobenius norm. + """ + return float(np.linalg.norm(a_mat - b_mat, "fro")) / float(np.linalg.norm(b_mat, "fro")) + + +def make_unit_vector(pos: int, nbits: int) -> np.ndarray: + """ + Makes a unit vector ``e = (0 ... 0 1 0 ... 0)`` of size ``2^n`` with + unit at the position ``num``. **Note**: result depends on bit ordering. + + Args: + pos: position of unit in vector. + nbits: number of meaningful bit in the number "pos". + + Returns: + unit vector of size ``2^n``. + """ + vec = np.zeros((2**nbits,), dtype=np.int64) + vec[fgu.reverse_bits(pos, nbits, enable=True)] = 1 + return vec + + +def eye_int(n: int) -> np.ndarray: + """ + Creates an identity matrix with integer entries. + + Args: + n: number of bits. + + Returns: + unit matrix of size ``2^n`` with integer entries. + """ + return np.eye(2**n, dtype=np.int64) + + +def kron3(a_mat: np.ndarray, b_mat: np.ndarray, c_mat: np.ndarray) -> np.ndarray: + """ + Computes Kronecker product of 3 matrices. + """ + return np.kron(a_mat, np.kron(b_mat, c_mat)) + + +def kron5( + a_mat: np.ndarray, b_mat: np.ndarray, c_mat: np.ndarray, d_mat: np.ndarray, e_mat: np.ndarray +) -> np.ndarray: + """ + Computes Kronecker product of 5 matrices. + """ + return np.kron(np.kron(np.kron(np.kron(a_mat, b_mat), c_mat), d_mat), e_mat) + + +def rand_matrix(dim: int, kind: str = "complex") -> np.ndarray: + """ + Generates a random complex or integer value matrix. + + Args: + dim: matrix size dim-x-dim. + kind: "complex" or "randint". + + Returns: + a random matrix. + """ + if kind == "complex": + return ( + np.random.rand(dim, dim).astype(np.complex128) + + np.random.rand(dim, dim).astype(np.complex128) * 1j + ) + else: + return np.random.randint(low=1, high=100, size=(dim, dim), dtype=np.int64) + + +def make_test_matrices2x2(n: int, k: int, kind: str = "complex") -> Tuple[np.ndarray, np.ndarray]: + """ + Creates a ``2^n x 2^n`` random matrix made as a Kronecker product of identity + ones and a single 1-qubit gate. This models a layer in quantum circuit with + an arbitrary 1-qubit gate somewhere in the middle. + + Args: + n: number of qubits. + k: index of qubit a 1-qubit gate is acting on. + kind: entries of the output matrix are defined as: + "complex", "primes" or "randint". + Returns: + A tuple of (1) ``2^n x 2^n`` random matrix; (2) ``2 x 2`` matrix of 1-qubit + gate used for matrix construction. + """ + if kind == "primes": + a_mat = np.asarray([[2, 3], [5, 7]], dtype=np.int64) + else: + a_mat = rand_matrix(dim=2, kind=kind) + m_mat = kron3(eye_int(k), a_mat, eye_int(n - k - 1)) + return m_mat, a_mat + + +def make_test_matrices4x4( + n: int, j: int, k: int, kind: str = "complex" +) -> Tuple[np.ndarray, np.ndarray]: + """ + Creates a ``2^n x 2^n`` random matrix made as a Kronecker product of identity + ones and a single 2-qubit gate. This models a layer in quantum circuit with + an arbitrary 2-qubit gate somewhere in the middle. + + Args: + n: number of qubits. + j: index of the first qubit the 2-qubit gate acting on. + k: index of the second qubit the 2-qubit gate acting on. + kind: entries of the output matrix are defined as: + "complex", "primes" or "randint". + + Returns: + A tuple of (1) ``2^n x 2^n`` random matrix; (2) ``4 x 4`` matrix of + 2-qubit gate used for matrix construction. + """ + if kind == "primes": + a_mat = np.asarray([[2, 3], [5, 7]], dtype=np.int64) + b_mat = np.asarray([[11, 13], [17, 19]], dtype=np.int64) + c_mat = np.asarray([[47, 53], [41, 43]], dtype=np.int64) + d_mat = np.asarray([[31, 37], [23, 29]], dtype=np.int64) + else: + a_mat, b_mat = rand_matrix(dim=2, kind=kind), rand_matrix(dim=2, kind=kind) + c_mat, d_mat = rand_matrix(dim=2, kind=kind), rand_matrix(dim=2, kind=kind) + + if j < k: + m_mat = kron5(eye_int(j), a_mat, eye_int(k - j - 1), b_mat, eye_int(n - k - 1)) + kron5( + eye_int(j), c_mat, eye_int(k - j - 1), d_mat, eye_int(n - k - 1) + ) + else: + m_mat = kron5(eye_int(k), b_mat, eye_int(j - k - 1), a_mat, eye_int(n - j - 1)) + kron5( + eye_int(k), d_mat, eye_int(j - k - 1), c_mat, eye_int(n - j - 1) + ) + g_mat = np.kron(a_mat, b_mat) + np.kron(c_mat, d_mat) + return m_mat, g_mat + + +def rand_circuit(num_qubits: int, depth: int) -> np.ndarray: + """ + Generates a random circuit of unit blocks for debugging and testing. + """ + blocks = np.tile(np.arange(num_qubits).reshape(num_qubits, 1), depth) + for i in range(depth): + np.random.shuffle(blocks[:, i]) + return blocks[0:2, :].copy() + + +def rand_su_mat(dim: int) -> np.ndarray: + """ + Generates a random SU matrix. + Args: + dim: matrix size ``dim-x-dim``. + Returns: + random SU matrix. + """ + u_mat = unitary_group.rvs(dim) + u_mat /= np.linalg.det(u_mat) ** (1.0 / float(dim)) + return u_mat diff --git a/test/python/synthesis/aqc/sample_data.py b/test/python/synthesis/aqc/sample_data.py new file mode 100644 index 000000000000..d2537a0a163f --- /dev/null +++ b/test/python/synthesis/aqc/sample_data.py @@ -0,0 +1,425 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2023. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +"""Sample data for the AQC tests.""" + +import numpy as np + +# Circuit to optimize +ORIGINAL_CIRCUIT = np.asarray( + [ + [ + -0.07783734 + 0.22716113j, + -0.00910959 + 0.04320976j, + -0.21030301 + 0.20425801j, + -0.59549525 + 0.23857045j, + 0.26052641 - 0.21291574j, + 0.04239446 - 0.51982577j, + 0.11510024 + 0.0862973j, + -0.14629519 - 0.1248565j, + ], + [ + 0.36992751 + 0.24267278j, + -0.24662827 - 0.1550994j, + 0.01884297 - 0.19588483j, + -0.12679802 - 0.03807062j, + -0.47228903 - 0.23051667j, + 0.25913925 - 0.07348255j, + -0.33918991 - 0.31963187j, + -0.06086797 - 0.30571598j, + ], + [ + 0.0222799 - 0.19736565j, + -0.11059812 + 0.005871j, + -0.09659148 + 0.37062522j, + 0.10305057 + 0.2136938j, + 0.43803384 - 0.18087079j, + 0.2915554 + 0.21651464j, + -0.5518552 - 0.2278622j, + -0.07433388 + 0.16384345j, + ], + [ + -0.10575111 - 0.20500221j, + -0.46919823 + 0.28035591j, + -0.00532403 + 0.18394858j, + -0.40488016 - 0.02266676j, + -0.31535024 - 0.20471488j, + 0.04892042 + 0.4144072j, + 0.26845334 + 0.03135037j, + 0.02306506 + 0.24632236j, + ], + [ + -0.08478682 + 0.17291421j, + -0.00175922 + 0.29513196j, + -0.4834236 + 0.33435023j, + 0.40144427 - 0.05581059j, + -0.12199095 - 0.10961794j, + 0.03892304 + 0.15557211j, + 0.13831459 + 0.12635286j, + -0.14552454 - 0.50710571j, + ], + [ + 0.45282021 - 0.33188642j, + 0.2698306 + 0.52711561j, + 0.10654936 - 0.22260092j, + -0.10493819 - 0.16484395j, + 0.24892473 - 0.32326511j, + -0.09647525 + 0.02038148j, + 0.0574283 - 0.07771336j, + 0.07001849 - 0.21125733j, + ], + [ + -0.27951269 - 0.31668187j, + -0.19034985 + 0.03974792j, + 0.07497032 + 0.31316668j, + 0.02401214 - 0.27077229j, + -0.0952781 + 0.06841821j, + -0.34445156 - 0.38097419j, + -0.11546146 - 0.43206826j, + 0.28988791 - 0.21115787j, + ], + [ + 0.32899805 + 0.12250135j, + 0.32760646 - 0.13874956j, + -0.17776306 + 0.39125565j, + -0.26523096 - 0.05930179j, + -0.157175 + 0.11327617j, + -0.12985701 + 0.18256439j, + -0.19442067 + 0.22433888j, + 0.56113453 + 0.03006748j, + ], + ] +) + + +INITIAL_THETAS = np.asarray( + [ + 5.840950175706398, + 1.9878462360990703, + 1.1555959752434493, + 1.2852901366370226, + 3.5671215612441567, + 3.741917727527864, + 6.0602234589642325, + 4.104032738049671, + 4.705519181402343, + 4.106500609750883, + 4.698030703573854, + 6.040068359785516, + 0.05270523037861366, + 0.6688097436989993, + 1.8768107855557548, + 4.124353101010888, + 5.088202331694798, + 5.480042886388346, + 6.061059610413401, + 4.547049138788613, + 4.036791540460537, + 4.507894048871352, + 2.9380112118120447, + 2.045708862158222, + 2.7623685280761947, + 4.584771723524366, + 6.245577841248841, + 4.252922960633062, + 4.9688844247279915, + 1.073885953348928, + 0.16869897515173574, + 5.028874556936197, + 5.678256173892097, + 0.15504520280603118, + 3.0897395261041516, + 3.3065587353199613, + 3.74707815436836, + 0.3264588839853597, + 5.624013371278708, + 4.575831363947143, + 5.141844767695358, + 3.1429922509263313, + 5.090570189425126, + 0.6029880309099467, + 1.3757036977838175, + 1.6255798065519669, + 2.941195195443305, + 2.8863269571102226, + 4.457981426310531, + 1.118740030466607, + 3.339198104910743, + 1.0539555073048037, + 4.830600316082553, + 5.831867556758204, + 3.829561596723921, + 0.9436307270844012, + 3.0764153106528846, + 2.3709282696632825, + 5.331919923390631, + 5.724592720350565, + 2.411792644955149, + 1.9823192246938757, + 3.571325789530997, + 1.1800955180886494, + 0.7906857392089286, + ] +) + + +CARTAN_3 = [ + [0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0], + [1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1], +] + + +CARTAN_4 = [ + [ + 0, + 0, + 0, + 2, + 1, + 2, + 1, + 0, + 0, + 0, + 2, + 1, + 2, + 1, + 0, + 0, + 0, + 2, + 1, + 2, + 1, + 0, + 0, + 0, + 2, + 1, + 2, + 0, + 2, + 1, + 2, + 0, + 0, + 0, + 0, + 2, + 1, + 2, + 1, + 0, + 0, + 0, + 2, + 1, + 2, + 1, + 0, + 0, + 0, + 2, + 1, + 2, + 1, + 0, + 0, + 0, + 2, + 1, + 2, + 0, + 2, + 1, + 2, + 0, + 0, + 0, + 0, + 2, + 1, + 2, + 1, + 0, + 0, + 0, + 2, + 1, + 2, + 1, + 0, + 0, + 0, + 2, + 1, + 2, + 1, + 0, + 0, + 0, + 2, + 1, + 2, + 0, + 2, + 1, + 2, + 0, + 0, + 0, + 0, + 2, + 1, + 2, + 1, + 0, + 0, + 0, + 2, + 1, + 2, + 1, + 0, + 0, + 0, + 2, + 1, + 2, + 1, + 0, + 0, + 0, + ], + [ + 1, + 1, + 1, + 3, + 3, + 3, + 3, + 1, + 1, + 1, + 3, + 3, + 3, + 3, + 1, + 1, + 1, + 3, + 3, + 3, + 3, + 1, + 1, + 1, + 3, + 3, + 3, + 3, + 3, + 3, + 3, + 3, + 1, + 1, + 1, + 3, + 3, + 3, + 3, + 1, + 1, + 1, + 3, + 3, + 3, + 3, + 1, + 1, + 1, + 3, + 3, + 3, + 3, + 1, + 1, + 1, + 3, + 3, + 3, + 3, + 3, + 3, + 3, + 3, + 1, + 1, + 1, + 3, + 3, + 3, + 3, + 1, + 1, + 1, + 3, + 3, + 3, + 3, + 1, + 1, + 1, + 3, + 3, + 3, + 3, + 1, + 1, + 1, + 3, + 3, + 3, + 3, + 3, + 3, + 3, + 3, + 1, + 1, + 1, + 3, + 3, + 3, + 3, + 1, + 1, + 1, + 3, + 3, + 3, + 3, + 1, + 1, + 1, + 3, + 3, + 3, + 3, + 1, + 1, + 1, + ], +] diff --git a/test/python/synthesis/aqc/test_aqc.py b/test/python/synthesis/aqc/test_aqc.py new file mode 100644 index 000000000000..d045281b6101 --- /dev/null +++ b/test/python/synthesis/aqc/test_aqc.py @@ -0,0 +1,130 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2022. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +""" +Tests AQC framework using hardcoded and randomly generated circuits. +""" +import unittest +from test.python.transpiler.aqc.sample_data import ORIGINAL_CIRCUIT, INITIAL_THETAS +import numpy as np +from qiskit.algorithms.optimizers import L_BFGS_B +from qiskit.quantum_info import Operator +from qiskit.test import QiskitTestCase +from qiskit.synthesis.unitary.aqc.aqc import AQC +from qiskit.synthesis.unitary.aqc.cnot_unit_circuit import CNOTUnitCircuit +from qiskit.synthesis.unitary.aqc.cnot_structures import make_cnot_network +from qiskit.synthesis.unitary.aqc.cnot_unit_objective import DefaultCNOTUnitObjective +from qiskit.synthesis.unitary.aqc.fast_gradient.fast_gradient import FastCNOTUnitObjective + + +class TestAqc(QiskitTestCase): + """Main tests of approximate quantum compiler.""" + + def test_aqc(self): + """Tests AQC on a hardcoded circuit/matrix.""" + + seed = 12345 + + num_qubits = int(round(np.log2(ORIGINAL_CIRCUIT.shape[0]))) + cnots = make_cnot_network( + num_qubits=num_qubits, network_layout="spin", connectivity_type="full", depth=0 + ) + + optimizer = L_BFGS_B(maxiter=200) + + aqc = AQC(optimizer=optimizer, seed=seed) + + target_matrix = ORIGINAL_CIRCUIT + approximate_circuit = CNOTUnitCircuit(num_qubits, cnots) + approximating_objective = DefaultCNOTUnitObjective(num_qubits, cnots) + + aqc.compile_unitary( + target_matrix=target_matrix, + approximate_circuit=approximate_circuit, + approximating_objective=approximating_objective, + initial_point=INITIAL_THETAS, + ) + + approx_matrix = Operator(approximate_circuit).data + error = 0.5 * (np.linalg.norm(approx_matrix - ORIGINAL_CIRCUIT, "fro") ** 2) + self.assertTrue(error < 1e-3) + + def test_aqc_fastgrad(self): + """ + Tests AQC on a MCX circuit/matrix with random initial guess using + the accelerated implementation. + """ + seed = 12345 + np.random.seed(seed) + + num_qubits = int(3) + cnots = make_cnot_network( + num_qubits=num_qubits, network_layout="spin", connectivity_type="full", depth=0 + ) + + optimizer = L_BFGS_B(maxiter=200) + aqc = AQC(optimizer=optimizer, seed=seed) + + # Make multi-control CNOT gate matrix. + # Another option: target_matrix = ORIGINAL_CIRCUIT + target_matrix = np.eye(2**num_qubits, dtype=np.complex128) + target_matrix[-2:, -2:] = [[0, 1], [1, 0]] + + circ = CNOTUnitCircuit(num_qubits, cnots) + objv = FastCNOTUnitObjective(num_qubits, cnots) + + aqc.compile_unitary( + target_matrix=target_matrix, + approximate_circuit=circ, + approximating_objective=objv, + initial_point=2 * np.pi * np.random.rand(objv.num_thetas), + ) + + approx_matrix = Operator(circ).data + error = 0.5 * (np.linalg.norm(approx_matrix - target_matrix, "fro") ** 2) + self.assertTrue(error < 1e-3) + + def test_aqc_determinant_minus_one(self): + """ + Tests AQC on a matrix with determinant −1. + """ + seed = 12345 + + num_qubits = int(3) + cnots = make_cnot_network( + num_qubits=num_qubits, network_layout="spin", connectivity_type="full", depth=0 + ) + + optimizer = L_BFGS_B(maxiter=200) + aqc = AQC(optimizer=optimizer, seed=seed) + + target_matrix = np.eye(2**num_qubits, dtype=int) + # Make a unitary with determinant −1 by swapping any two columns + target_matrix[:, 2], target_matrix[:, 3] = target_matrix[:, 3], target_matrix[:, 2].copy() + + approximate_circuit = CNOTUnitCircuit(num_qubits, cnots) + approximating_objective = DefaultCNOTUnitObjective(num_qubits, cnots) + + aqc.compile_unitary( + target_matrix=target_matrix, + approximate_circuit=approximate_circuit, + approximating_objective=approximating_objective, + initial_point=INITIAL_THETAS, + ) + + approx_matrix = Operator(approximate_circuit).data + + error = 0.5 * (np.linalg.norm(approx_matrix - target_matrix, "fro") ** 2) + self.assertTrue(error < 1e-3) + + +if __name__ == "__main__": + unittest.main() diff --git a/test/python/synthesis/aqc/test_aqc_plugin.py b/test/python/synthesis/aqc/test_aqc_plugin.py new file mode 100644 index 000000000000..1193201335b4 --- /dev/null +++ b/test/python/synthesis/aqc/test_aqc_plugin.py @@ -0,0 +1,102 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2023. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +""" +Tests AQC plugin. +""" + +import numpy as np + +from qiskit import QuantumCircuit +from qiskit.algorithms.optimizers import SLSQP +from qiskit.converters import dag_to_circuit, circuit_to_dag +from qiskit.quantum_info import Operator +from qiskit.test import QiskitTestCase +from qiskit.transpiler import PassManager +from qiskit.transpiler.passes import UnitarySynthesis +from qiskit.synthesis.unitary.aqc.aqc_plugin import AQCSynthesisPlugin + + +class TestAQCSynthesisPlugin(QiskitTestCase): + """Basic tests of the AQC synthesis plugin.""" + + def setUp(self): + super().setUp() + self._qc = QuantumCircuit(3) + self._qc.mcx( + [ + 0, + 1, + ], + 2, + ) + + self._target_unitary = Operator(self._qc).data + self._seed_config = {"seed": 12345} + + def test_aqc_plugin(self): + """Basic test of the plugin.""" + plugin = AQCSynthesisPlugin() + dag = plugin.run(self._target_unitary, config=self._seed_config) + + approx_circuit = dag_to_circuit(dag) + approx_unitary = Operator(approx_circuit).data + + np.testing.assert_array_almost_equal(self._target_unitary, approx_unitary, 3) + + def test_plugin_setup(self): + """Tests the plugin via unitary synthesis pass""" + transpiler_pass = UnitarySynthesis( + basis_gates=["rx", "ry", "rz", "cx"], method="aqc", plugin_config=self._seed_config + ) + + dag = circuit_to_dag(self._qc) + dag = transpiler_pass.run(dag) + + approx_circuit = dag_to_circuit(dag) + approx_unitary = Operator(approx_circuit).data + + np.testing.assert_array_almost_equal(self._target_unitary, approx_unitary, 3) + + def test_plugin_configuration(self): + """Tests plugin with a custom configuration.""" + config = { + "network_layout": "sequ", + "connectivity_type": "full", + "depth": 0, + "seed": 12345, + "optimizer": SLSQP(), + } + transpiler_pass = UnitarySynthesis( + basis_gates=["rx", "ry", "rz", "cx"], method="aqc", plugin_config=config + ) + + dag = circuit_to_dag(self._qc) + dag = transpiler_pass.run(dag) + + approx_circuit = dag_to_circuit(dag) + approx_unitary = Operator(approx_circuit).data + + np.testing.assert_array_almost_equal(self._target_unitary, approx_unitary, 3) + + def test_with_pass_manager(self): + """Tests the plugin via pass manager""" + qc = QuantumCircuit(3) + qc.unitary(np.eye(8), [0, 1, 2]) + aqc = PassManager( + [ + UnitarySynthesis( + basis_gates=["u", "cx"], method="aqc", plugin_config=self._seed_config + ) + ] + ).run(qc) + approx_unitary = Operator(aqc).data + np.testing.assert_array_almost_equal(np.eye(8), approx_unitary, 3) diff --git a/test/python/synthesis/aqc/test_cnot_networks.py b/test/python/synthesis/aqc/test_cnot_networks.py new file mode 100644 index 000000000000..b31f401017ed --- /dev/null +++ b/test/python/synthesis/aqc/test_cnot_networks.py @@ -0,0 +1,54 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2023. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +""" +Tests building up CNOT unit structures. +""" +from test.python.transpiler.aqc.sample_data import CARTAN_4, CARTAN_3 + +import numpy as np +from ddt import ddt, data, unpack + +from qiskit.test import QiskitTestCase +from qiskit.synthesis.unitary.aqc import make_cnot_network + + +@ddt +class TestCNOTNetworks(QiskitTestCase): + """Tests constructing various CNOT structures.""" + + @data( + # number of qubits, network layout, connectivity, depth, output + (3, "spin", "full", 3, [[0, 1, 0], [1, 2, 1]]), + (4, "spin", "full", 4, [[0, 2, 1, 0], [1, 3, 2, 1]]), + (3, "sequ", "full", 4, [[0, 0, 1, 0], [1, 2, 2, 1]]), + (4, "sequ", "full", 7, [[0, 0, 0, 1, 1, 2, 0], [1, 2, 3, 2, 3, 3, 1]]), + (3, "sequ", "line", 3, [[0, 1, 0], [1, 2, 1]]), + (4, "sequ", "line", 4, [[0, 1, 2, 0], [1, 2, 3, 1]]), + (3, "sequ", "star", 3, [[0, 0, 0], [1, 2, 1]]), + (4, "sequ", "star", 4, [[0, 0, 0, 0], [1, 2, 3, 1]]), + (3, "cart", "full", 0, CARTAN_3), + (4, "cart", "full", 0, CARTAN_4), + (3, "cyclic_spin", "full", 3, [[0, 1, 0], [1, 2, 1]]), + (4, "cyclic_spin", "full", 5, [[0, 2, 1, 3, 0], [1, 3, 2, 0, 1]]), + (3, "cyclic_line", "line", 4, [[0, 1, 2, 0], [1, 2, 0, 1]]), + (4, "cyclic_line", "line", 5, [[0, 1, 2, 3, 0], [1, 2, 3, 0, 1]]), + ) + @unpack + def test_basic_networks(self, num_qubits, network_layout, connectivity, depth, output): + """Tests basic CNOT networks.""" + cnots = make_cnot_network( + num_qubits=num_qubits, + network_layout=network_layout, + connectivity_type=connectivity, + depth=depth, + ) + np.testing.assert_array_equal(cnots, output) diff --git a/test/python/synthesis/aqc/test_gradient.py b/test/python/synthesis/aqc/test_gradient.py new file mode 100644 index 000000000000..132c85aa2cc6 --- /dev/null +++ b/test/python/synthesis/aqc/test_gradient.py @@ -0,0 +1,108 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2023. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +""" +Tests analytical gradient vs the one computed via finite differences. +""" + +import unittest +from test.python.synthesis.aqc.sample_data import ORIGINAL_CIRCUIT, INITIAL_THETAS +import numpy as np +from qiskit.test import QiskitTestCase +from qiskit.synthesis.unitary.aqc.cnot_structures import make_cnot_network +from qiskit.synthesis.unitary.aqc.cnot_unit_objective import DefaultCNOTUnitObjective + + +class TestGradientAgainstFiniteDiff(QiskitTestCase): + """ + Compares analytical gradient vs the one computed via finite difference + approximation. Also, the test demonstrates that the difference between + analytical and numerical gradients is up to quadratic term in Taylor + expansion for small deltas. + """ + + def setUp(self): + super().setUp() + np.random.seed(0x0696969) + + def test_gradient(self): + """ + Gradient test for specified number of qubits and circuit depth. + """ + num_qubits = 3 + num_cnots = 14 + + cnots = make_cnot_network( + num_qubits=num_qubits, network_layout="spin", connectivity_type="full", depth=num_cnots + ) + + # we pick a target matrix from the existing sample data + target_matrix = ORIGINAL_CIRCUIT + objective = DefaultCNOTUnitObjective(num_qubits, cnots) + objective.target_matrix = target_matrix + + # thetas = np.random.rand(objective.num_thetas) * (2.0 * np.pi) + thetas = INITIAL_THETAS + fobj0 = objective.objective(thetas) + grad0 = objective.gradient(thetas) + + grad0_dir = grad0 / np.linalg.norm(grad0) + numerical_grad = np.zeros(thetas.size) + thetas_delta = np.zeros(thetas.size) + + # Every angle has a magnitude between 0 and 2*pi. We choose the + # angle increment (tau) about the same order of magnitude, tau <= 1, + # and then gradually decrease it towards zero. + tau = 1.0 + diff_prev = 0.0 + orders = [] + errors = [] + steps = 9 + for step in range(steps): + # Estimate gradient approximation error. + for i in range(thetas.size): + np.copyto(thetas_delta, thetas) + thetas_delta[i] -= tau + fobj1 = objective.objective(thetas_delta) + np.copyto(thetas_delta, thetas) + thetas_delta[i] += tau + fobj2 = objective.objective(thetas_delta) + numerical_grad[i] = (fobj2 - fobj1) / (2.0 * tau) + errors.append(np.linalg.norm(grad0 - numerical_grad) / np.linalg.norm(grad0)) + + # Estimate approximation order (should be quadratic for small tau). + # Note, we take perturbation in gradient direction. More rigorous + # approach would take a random direction, although quadratic + # convergence is less pronounced in this case. + perturbation = grad0_dir * tau + # circuit.set_thetas(thetas + perturbation) + fobj = objective.objective(thetas + perturbation) + diff = abs(fobj - fobj0 - np.dot(grad0, perturbation)) + orders.append( + 0.0 if step == 0 else float((np.log(diff_prev) - np.log(diff)) / np.log(2.0)) + ) + + tau /= 2.0 + diff_prev = diff + + # check errors + prev_error = errors[0] + for error in errors[1:]: + self.assertLess(error, prev_error * 0.75) + prev_error = error + + # check orders, skipping first zero + self.assertTrue(np.count_nonzero(np.asarray(orders[1:]) > 1.8) >= 3) + self.assertTrue(np.count_nonzero(np.asarray(orders[1:]) < 3.0) >= 3) + + +if __name__ == "__main__": + unittest.main() diff --git a/test/python/synthesis/test_cnot_phase_synthesis.py b/test/python/synthesis/test_cnot_phase_synthesis.py index d96d7cbc8448..8bf6c28a8cc3 100644 --- a/test/python/synthesis/test_cnot_phase_synthesis.py +++ b/test/python/synthesis/test_cnot_phase_synthesis.py @@ -75,7 +75,11 @@ def test_gray_synth(self, synth_func): [0, 1, 0, 0, 1, 0], ] angles = ["s", "t", "z", "s", "t", "t"] - c_gray = synth_func(cnots, angles) + if synth_func.__name__ == "graysynth": + with self.assertWarns(DeprecationWarning): + c_gray = synth_func(cnots, angles) + else: + c_gray = synth_func(cnots, angles) unitary_gray = UnitaryGate(Operator(c_gray)) # Create the circuit displayed above: @@ -133,7 +137,11 @@ def test_paper_example(self, synth_func): """ cnots = [[0, 1, 1, 1, 1, 1], [1, 0, 0, 1, 1, 1], [1, 0, 0, 1, 0, 0], [0, 0, 1, 0, 1, 0]] angles = ["t"] * 6 - c_gray = synth_func(cnots, angles) + if synth_func.__name__ == "graysynth": + with self.assertWarns(DeprecationWarning): + c_gray = synth_func(cnots, angles) + else: + c_gray = synth_func(cnots, angles) unitary_gray = UnitaryGate(Operator(c_gray)) # Create the circuit displayed above: @@ -186,7 +194,11 @@ def test_ccz(self, synth_func): """ cnots = [[1, 0, 0, 1, 1, 0, 1], [0, 1, 0, 1, 0, 1, 1], [0, 0, 1, 0, 1, 1, 1]] angles = ["t", "t", "t", "tdg", "tdg", "tdg", "t"] - c_gray = synth_func(cnots, angles) + if synth_func.__name__ == "graysynth": + with self.assertWarns(DeprecationWarning): + c_gray = synth_func(cnots, angles) + else: + c_gray = synth_func(cnots, angles) unitary_gray = UnitaryGate(Operator(c_gray)) # Create the circuit displayed above: @@ -252,7 +264,11 @@ def test_patel_markov_hayes(self, synth_func): [1, 1, 0, 1, 1, 1], [0, 0, 1, 1, 1, 0], ] - c_patel = synth_func(state) + if synth_func.__name__ == "cnot_synth": + with self.assertWarns(DeprecationWarning): + c_patel = synth_func(state) + else: + c_patel = synth_func(state) unitary_patel = UnitaryGate(Operator(c_patel)) # Create the circuit displayed above: