Add the Solovay-Kitaev algorithm to the transpiler passes (#5657)

* add sk pass

* add utils

* add sk skeleton

* add example test file

* Add implementations initial version

* make the test run

* fix lint

* fix phase, add test

* fix the phase calculation

* add accuracy test

* cleanup of unused methods, append more efficient

* Add unittests for Solovay Kitaev Utils

* Add unitttests for commutator_decompose

* Remove unittests for non-public functions

* Fix pylint issues

* refactor simplify

* fix the test imports

* improve and add docstrings

* fix lint in testfile

* fix lint in utils

* fix lint in sk

* Add unittests

* fix phase, lint and allow basis gates as strings

* rm print

* fix adjoint and dot

* Fix unittests and add docstring to them

* move to circuit to GateSequence

* Fix assertion

* allow caching the basic approxs

* Fix bug in append from GateSequence

* Remove unnecesssary code

* Fix bug in test setup

* Add unittests for multiqubit QFT-circuit

* Add unittests for QFT with more qubits

* make compatible with current tests

* Remove unnecessary testcases

* Remove unnecessary unittests

* Fix bug in unittests

* Add release notes

* import scipy; scipy.optimize.fsolve() does not work

scipy.optimize has to be imported, since it is not imported when scipy is imported

* numpy

* Update qiskit/transpiler/passes/synthesis/solovay_kitaev_utils.py

Co-authored-by: Luciano Bello <bel@zurich.ibm.com>

* document ValueError

* rm unused methods, add GateSequence tests

* fix lint

* try fixing lint and docs

* cleanup

* improve readability of math
* move commutator_decompose into utils
* rm tests that are never run
* rm trailing todos and commented code

* add default dataset, cleanup tests

* rm unused imports

* replace diamond norm by trace distance

* rm unused imports

* fix deprecated use of DAGNode.type

* black

* really black!

* fail nicely

* test

* attempt to fix default approx path

* Move decompositions into py file dict

* speed up the basic approx generation

* add candidate checking by string

further reduces the runtime by a factor of about 1.8!

* use a tree for basic approx generation

Co-authored-by: Jake Lishman <jake.lishman@ibm.com>

* fix tests

* more speedups!

* use kdtree!

* refactor: generation as sep. function

* remove this masssssive file!

* start plugin work

* first attempt at synthesis plugin

* plugin itself works

-- but not when called via transpile or the UnitarySynthesis

* unitary synth works!

* use class-level method for SolovayKitaev

* docstrings

* add tests which cover still missing feature

* rename to SKSynth and better reno

* re-organize files

- algorithm to qiskit.synthesis
- plugin to qiskit.transpiler.passes.synthesis

* cover sklearn-free case

* missing future import

* move sk pass to transpiler subdirectory

* reno & try fixing slow optional module-level imports

* attempt 2 to fix sklearn import

* comments from Matthew's review

* directly construct dag

* (previous commit incomplete: missed this here)

* add SK to plugin toctree

* try fixing sphinx

* Apply Jake's comments

Removed the SK plugin from the plugin.py docs for now

* Fix doc example

* Update docs

* Fix tests

* Fix docs error

Co-authored-by: Julien Gacon <gaconju@gmail.com>
Co-authored-by: Luciano Bello <bel@zurich.ibm.com>
Co-authored-by: Jake Lishman <jake.lishman@ibm.com>
Co-authored-by: Matthew Treinish <mtreinish@kortar.org>
This commit is contained in:
LNoorl 2022-12-09 23:04:59 +01:00 committed by GitHub
parent b4f4c342c2
commit 9c1c18cea4
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
14 changed files with 1734 additions and 3 deletions

View File

@ -15,6 +15,8 @@
Circuit Synthesis (:mod:`qiskit.synthesis`)
===========================================
.. automodule:: qiskit.synthesis.discrete_basis
.. currentmodule:: qiskit.synthesis
Evolution Synthesis

View File

@ -0,0 +1,91 @@
# This code is part of Qiskit.
#
# (C) Copyright IBM 2017, 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.
r"""
=======================================================================================
Approximately decompose 1q gates to a discrete basis using the Solovay-Kitaev algorithm
=======================================================================================
.. currentmodule:: qiskit.synthesis.discrete_basis
Approximately decompose 1q gates to a discrete basis using the Solovay-Kitaev algorithm.
The Solovay-Kitaev theorem [1] states that any single qubit gate can be approximated to
arbitrary precision by a set of fixed single-qubit gates, if the set generates a dense
subset in :math:`SU(2)`. This is an important result, since it means that any single-qubit
gate can be expressed in terms of a discrete, universal gate set that we know how to implement
fault-tolerantly. Therefore, the Solovay-Kitaev algorithm allows us to take any
non-fault tolerant circuit and rephrase it in a fault-tolerant manner.
This implementation of the Solovay-Kitaev algorithm is based on [2].
For example, the following circuit
.. parsed-literal::
q_0: RX(0.8)
can be decomposed into
.. parsed-literal::
global phase: 7π/8
q_0: H T H
with an L2-error of approximately 0.01.
Examples:
.. jupyter-execute::
import numpy as np
from qiskit.circuit import QuantumCircuit
from qiskit.circuit.library import TGate, HGate, TdgGate
from qiskit.converters import circuit_to_dag, dag_to_circuit
from qiskit.transpiler.passes.synthesis import SolovayKitaev
from qiskit.quantum_info import Operator
circuit = QuantumCircuit(1)
circuit.rx(0.8, 0)
dag = circuit_to_dag(circuit)
print('Original circuit:')
print(circuit.draw())
basis_gates = [TGate(), TdgGate(), HGate()]
skd = SolovayKitaev(recursion_degree=2)
discretized = dag_to_circuit(skd.run(dag))
print('Discretized circuit:')
print(discretized.draw())
print('Error:', np.linalg.norm(Operator(circuit).data - Operator(discretized).data))
References:
[1]: Kitaev, A Yu (1997). Quantum computations: algorithms and error correction.
Russian Mathematical Surveys. 52 (6): 11911249.
`Online <https://iopscience.iop.org/article/10.1070/RM1997v052n06ABEH002155>`_.
[2]: Dawson, Christopher M.; Nielsen, Michael A. (2005) The Solovay-Kitaev Algorithm.
`arXiv:quant-ph/0505030 <https://arxiv.org/abs/quant-ph/0505030>`_
"""
from .solovay_kitaev import SolovayKitaevDecomposition

View File

@ -0,0 +1,245 @@
# This code is part of Qiskit.
#
# (C) Copyright IBM 2017, 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.
"""Functions to compute the decomposition of an SO(3) matrix as balanced commutator."""
from __future__ import annotations
import math
import numpy as np
from .gate_sequence import _check_is_so3, GateSequence
def _compute_trace_so3(matrix: np.ndarray) -> float:
"""Computes trace of an SO(3)-matrix.
Args:
matrix: an SO(3)-matrix
Returns:
Trace of ``matrix``.
Raises:
ValueError: if ``matrix`` is not an SO(3)-matrix.
"""
_check_is_so3(matrix)
trace = np.matrix.trace(matrix)
trace_rounded = min(trace, 3)
return trace_rounded
def _compute_rotation_axis(matrix: np.ndarray) -> np.ndarray:
"""Computes rotation axis of SO(3)-matrix.
Args:
matrix: The SO(3)-matrix for which rotation angle needs to be computed.
Returns:
The rotation axis of the SO(3)-matrix ``matrix``.
Raises:
ValueError: if ``matrix`` is not an SO(3)-matrix.
"""
_check_is_so3(matrix)
trace = _compute_trace_so3(matrix)
theta = math.acos(0.5 * (trace - 1))
if math.sin(theta) > 1e-10:
x = 1 / (2 * math.sin(theta)) * (matrix[2][1] - matrix[1][2])
y = 1 / (2 * math.sin(theta)) * (matrix[0][2] - matrix[2][0])
z = 1 / (2 * math.sin(theta)) * (matrix[1][0] - matrix[0][1])
else:
x = 1.0
y = 0.0
z = 0.0
return np.array([x, y, z])
def _solve_decomposition_angle(matrix: np.ndarray) -> float:
"""Computes angle for balanced commutator of SO(3)-matrix ``matrix``.
Computes angle a so that the SO(3)-matrix ``matrix`` can be decomposed
as commutator [v,w] where v and w are both rotations of a about some axis.
The computation is done by solving a trigonometric equation using scipy.optimize.fsolve.
Args:
matrix: The SO(3)-matrix for which the decomposition angle needs to be computed.
Returns:
Angle a so that matrix = [v,w] with v and w rotations of a about some axis.
Raises:
ValueError: if ``matrix`` is not an SO(3)-matrix.
"""
from scipy.optimize import fsolve
_check_is_so3(matrix)
trace = _compute_trace_so3(matrix)
angle = math.acos((1 / 2) * (trace - 1))
def objective(phi):
rhs = 2 * math.sin(phi / 2) ** 2
rhs *= math.sqrt(1 - math.sin(phi / 2) ** 4)
lhs = math.sin(angle / 2)
return rhs - lhs
decomposition_angle = fsolve(objective, angle)[0]
return decomposition_angle
def _compute_rotation_from_angle_and_axis( # pylint: disable=invalid-name
angle: float, axis: np.ndarray
) -> np.ndarray:
"""Computes the SO(3)-matrix corresponding to the rotation of ``angle`` about ``axis``.
Args:
angle: The angle of the rotation.
axis: The axis of the rotation.
Returns:
SO(3)-matrix that represents a rotation of ``angle`` about ``axis``.
Raises:
ValueError: if ``axis`` is not a 3-dim unit vector.
"""
if axis.shape != (3,):
raise ValueError(f"Axis must be a 1d array of length 3, but has shape {axis.shape}.")
if abs(np.linalg.norm(axis) - 1.0) > 1e-4:
raise ValueError(f"Axis must have a norm of 1, but has {np.linalg.norm(axis)}.")
res = math.cos(angle) * np.identity(3) + math.sin(angle) * _cross_product_matrix(axis)
res += (1 - math.cos(angle)) * np.outer(axis, axis)
return res
def _compute_rotation_between(from_vector: np.ndarray, to_vector: np.ndarray) -> np.ndarray:
"""Computes the SO(3)-matrix for rotating ``from_vector`` to ``to_vector``.
Args:
from_vector: unit vector of size 3
to_vector: unit vector of size 3
Returns:
SO(3)-matrix that brings ``from_vector`` to ``to_vector``.
Raises:
ValueError: if at least one of ``from_vector`` of ``to_vector`` is not a 3-dim unit vector.
"""
from_vector = from_vector / np.linalg.norm(from_vector)
to_vector = to_vector / np.linalg.norm(to_vector)
dot = np.dot(from_vector, to_vector)
cross = _cross_product_matrix(np.cross(from_vector, to_vector))
rotation_matrix = np.identity(3) + cross + np.dot(cross, cross) / (1 + dot)
return rotation_matrix
def _cross_product_matrix(v: np.ndarray) -> np.ndarray:
"""Computes cross product matrix from vector.
Args:
v: Vector for which cross product matrix needs to be computed.
Returns:
The cross product matrix corresponding to vector ``v``.
"""
return np.array([[0, -v[2], v[1]], [v[2], 0, -v[0]], [-v[1], v[0], 0]])
def _compute_commutator_so3(a: np.ndarray, b: np.ndarray) -> np.ndarray:
"""Computes the commutator of the SO(3)-matrices ``a`` and ``b``.
The computation uses the fact that the inverse of an SO(3)-matrix is equal to its transpose.
Args:
a: SO(3)-matrix
b: SO(3)-matrix
Returns:
The commutator [a,b] of ``a`` and ``b`` w
Raises:
ValueError: if at least one of ``a`` or ``b`` is not an SO(3)-matrix.
"""
_check_is_so3(a)
_check_is_so3(b)
a_dagger = np.conj(a).T
b_dagger = np.conj(b).T
return np.dot(np.dot(np.dot(a, b), a_dagger), b_dagger)
def commutator_decompose(
u_so3: np.ndarray, check_input: bool = True
) -> tuple[GateSequence, GateSequence]:
r"""Decompose an :math:`SO(3)`-matrix, :math:`U` as a balanced commutator.
This function finds two :math:`SO(3)` matrices :math:`V, W` such that the input matrix
equals
.. math::
U = V^\dagger W^\dagger V W.
For this decomposition, the following statement holds
.. math::
||V - I||_F, ||W - I||_F \leq \frac{\sqrt{||U - I||_F}}{2},
where :math:`I` is the identity and :math:`||\cdot ||_F` is the Frobenius norm.
Args:
u_so3: SO(3)-matrix that needs to be decomposed as balanced commutator.
check_input: If True, checks whether the input matrix is actually SO(3).
Returns:
Tuple of GateSequences from SO(3)-matrices :math:`V, W`.
Raises:
ValueError: if ``u_so3`` is not an SO(3)-matrix.
"""
if check_input:
# assert that the input matrix is really SO(3)
_check_is_so3(u_so3)
identity = np.identity(3)
if not (
np.allclose(u_so3.dot(u_so3.T), identity) and np.allclose(u_so3.T.dot(u_so3), identity)
):
raise ValueError("Input matrix is not orthogonal.")
angle = _solve_decomposition_angle(u_so3)
# Compute rotation about x-axis with angle 'angle'
vx = _compute_rotation_from_angle_and_axis(angle, np.array([1, 0, 0]))
# Compute rotation about y-axis with angle 'angle'
wy = _compute_rotation_from_angle_and_axis(angle, np.array([0, 1, 0]))
commutator = _compute_commutator_so3(vx, wy)
u_so3_axis = _compute_rotation_axis(u_so3)
commutator_axis = _compute_rotation_axis(commutator)
sim_matrix = _compute_rotation_between(commutator_axis, u_so3_axis)
sim_matrix_dagger = np.conj(sim_matrix).T
v = np.dot(np.dot(sim_matrix, vx), sim_matrix_dagger)
w = np.dot(np.dot(sim_matrix, wy), sim_matrix_dagger)
return GateSequence.from_matrix(v), GateSequence.from_matrix(w)

View File

@ -0,0 +1,414 @@
# This code is part of Qiskit.
#
# (C) Copyright IBM 2017, 2020.
#
# 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.
"""Algebra utilities and the ``GateSequence`` class."""
from __future__ import annotations
from collections.abc import Sequence
import math
import numpy as np
from qiskit.circuit import Gate, QuantumCircuit, Qubit
from qiskit.dagcircuit import DAGCircuit
from qiskit.extensions import UnitaryGate
class GateSequence:
"""A class implementing a sequence of gates.
This class stores the sequence of gates along with the unitary they implement.
"""
def __init__(self, gates: Sequence[Gate] = ()) -> None:
"""Create a new sequence of gates.
Args:
gates: The gates in the sequence. The default is [].
"""
self.gates = list(gates)
self.matrices = [np.asarray(gate, dtype=np.complex128) for gate in gates]
self.labels = [gate.name for gate in gates]
# get U(2) representation of the gate sequence
u2_matrix = np.identity(2)
for matrix in self.matrices:
# idea: could this be optimized by a specific numpy operation?
u2_matrix = matrix.dot(u2_matrix)
# convert to SU(2)
su2_matrix, global_phase = _convert_u2_to_su2(u2_matrix)
# convert to SO(3), that's what the Solovay Kitaev algorithm uses
so3_matrix = _convert_su2_to_so3(su2_matrix)
# store the matrix and the global phase
self._eulers = None
self.name = " ".join(self.labels)
self.global_phase = global_phase
self.product = so3_matrix
self.product_su2 = su2_matrix
def remove_cancelling_pair(self, indices: Sequence[int]) -> None:
"""Remove a pair of indices that cancel each other and *do not* change the matrices."""
for index in list(indices[::-1]):
self.gates.pop(index)
self.labels.pop(index)
# restore name
self.name = " ".join(self.labels)
def __eq__(self, other: "GateSequence") -> bool:
"""Check if this GateSequence is the same as the other GateSequence.
Args:
other: The GateSequence that will be compared to ``self``.
Returns:
True if ``other`` is equivalent to ``self``, false otherwise.
"""
if not len(self.gates) == len(other.gates):
return False
for gate1, gate2 in zip(self.gates, other.gates):
if gate1 != gate2:
return False
if self.global_phase != other.global_phase:
return False
return True
def to_circuit(self):
"""Convert to a circuit.
If no gates set but the product is not the identity, returns a circuit with a
unitary operation to implement the matrix.
"""
if len(self.gates) == 0 and not np.allclose(self.product, np.identity(3)):
circuit = QuantumCircuit(1, global_phase=self.global_phase)
su2 = _convert_so3_to_su2(self.product)
circuit.unitary(su2, [0])
return circuit
circuit = QuantumCircuit(1, global_phase=self.global_phase)
for gate in self.gates:
circuit.append(gate, [0])
return circuit
def to_dag(self):
"""Convert to a :class:`.DAGCircuit`.
If no gates set but the product is not the identity, returns a circuit with a
unitary operation to implement the matrix.
"""
qreg = [Qubit()]
dag = DAGCircuit()
dag.add_qubits(qreg)
if len(self.gates) == 0 and not np.allclose(self.product, np.identity(3)):
su2 = _convert_so3_to_su2(self.product)
dag.apply_operation_back(UnitaryGate(su2), qreg)
return dag
dag.global_phase = self.global_phase
for gate in self.gates:
dag.apply_operation_back(gate, qreg)
return dag
def append(self, gate: Gate) -> "GateSequence":
"""Append gate to the sequence of gates.
Args:
gate: The gate to be appended.
Returns:
GateSequence with ``gate`` appended.
"""
# invalidate euler angles and name
self._eulers = None
# TODO: this recomputes the product whenever we append something, which could be more
# efficient by storing the current matrix and just multiplying the input gate to it
# self.product = convert_su2_to_so3(self._compute_product(self.gates))
matrix = np.array(gate, dtype=np.complex128)
su2, phase = _convert_u2_to_su2(matrix)
so3 = _convert_su2_to_so3(su2)
self.product = so3.dot(self.product)
self.product_su2 = su2.dot(self.product_su2)
self.global_phase = self.global_phase + phase
self.gates.append(gate)
if len(self.labels) > 0:
self.name += f" {gate.name}"
else:
self.name = gate.name
self.labels.append(gate.name)
self.matrices.append(matrix)
return self
def adjoint(self) -> "GateSequence":
"""Get the complex conjugate."""
# We're initializing an empty GateSequence and set the state manually, as we can
# efficiently infer the adjoint values from the current value instead of recomputing them.
adjoint = GateSequence()
adjoint.gates = [gate.inverse() for gate in reversed(self.gates)]
adjoint.labels = [inv.name for inv in adjoint.gates]
adjoint.name = " ".join(adjoint.labels)
adjoint.product = np.conj(self.product).T
adjoint.product_su2 = np.conj(self.product_su2).T
adjoint.global_phase = -self.global_phase
return adjoint
def copy(self) -> "GateSequence":
"""Create copy of the sequence of gates.
Returns:
A new ``GateSequence`` containing copy of list of gates.
"""
out = type(self).__new__(type(self))
out.labels = self.labels.copy()
out.gates = self.gates.copy()
out.matrices = self.matrices.copy()
out.global_phase = self.global_phase
out.product = self.product.copy()
out.product_su2 = self.product_su2.copy()
out.name = self.name
out._eulers = self._eulers
return out
def __len__(self) -> int:
"""Return length of sequence of gates.
Returns:
Length of list containing gates.
"""
return len(self.gates)
def __getitem__(self, index: int) -> Gate:
"""Returns the gate at ``index`` from the list of gates.
Args
index: Index of gate in list that will be returned.
Returns:
The gate at ``index`` in the list of gates.
"""
return self.gates[index]
def __repr__(self) -> str:
"""Return string representation of this object.
Returns:
Representation of this sequence of gates.
"""
out = "["
for gate in self.gates:
out += gate.name
out += ", "
out += "]"
out += ", product: "
out += str(self.product)
return out
def __str__(self) -> str:
"""Return string representation of this object.
Returns:
Representation of this sequence of gates.
"""
out = "["
for gate in self.gates:
out += gate.name
out += ", "
out += "]"
out += ", product: \n"
out += str(self.product)
return out
@classmethod
def from_matrix(cls, matrix: np.ndarray) -> "GateSequence":
"""Initialize the gate sequence from a matrix, without a gate sequence.
Args:
matrix: The matrix, can be SU(2) or SO(3).
Returns:
A ``GateSequence`` initialized from the input matrix.
Raises:
ValueError: If the matrix has an invalid shape.
"""
instance = cls()
if matrix.shape == (2, 2):
instance.product = _convert_su2_to_so3(matrix)
elif matrix.shape == (3, 3):
instance.product = matrix
else:
raise ValueError(f"Matrix must have shape (3, 3) or (2, 2) but has {matrix.shape}.")
instance.gates = []
return instance
def dot(self, other: "GateSequence") -> "GateSequence":
"""Compute the dot-product with another gate sequence.
Args:
other: The other gate sequence.
Returns:
The dot-product as gate sequence.
"""
# We're initializing an empty GateSequence and set the state manually, as we can more
# efficiently compute the multiplied values from the already constructed matrices.
composed = GateSequence()
composed.gates = other.gates + self.gates
composed.labels = other.labels + self.labels
composed.name = " ".join(composed.labels)
composed.product = np.dot(self.product, other.product)
composed.global_phase = self.global_phase + other.global_phase
return composed
def _convert_u2_to_su2(u2_matrix: np.ndarray) -> tuple[np.ndarray, float]:
"""Convert a U(2) matrix to SU(2) by adding a global phase."""
z = 1 / np.sqrt(np.linalg.det(u2_matrix))
su2_matrix = z * u2_matrix
phase = np.arctan2(np.imag(z), np.real(z))
return su2_matrix, phase
def _compute_euler_angles_from_so3(matrix: np.ndarray) -> tuple[float, float, float]:
"""Computes the Euler angles from the SO(3)-matrix u.
Uses the algorithm from Gregory Slabaugh,
see `here <https://www.gregslabaugh.net/publications/euler.pdf>`_.
Args:
matrix: The SO(3)-matrix for which the Euler angles need to be computed.
Returns:
Tuple (phi, theta, psi), where phi is rotation about z-axis, theta rotation about y-axis
and psi rotation about x-axis.
"""
matrix = np.round(matrix, decimals=10)
if matrix[2][0] != 1 and matrix[2][1] != -1:
theta = -math.asin(matrix[2][0])
psi = math.atan2(matrix[2][1] / math.cos(theta), matrix[2][2] / math.cos(theta))
phi = math.atan2(matrix[1][0] / math.cos(theta), matrix[0][0] / math.cos(theta))
return phi, theta, psi
else:
phi = 0
if matrix[2][0] == 1:
theta = math.pi / 2
psi = phi + math.atan2(matrix[0][1], matrix[0][2])
else:
theta = -math.pi / 2
psi = -phi + math.atan2(-matrix[0][1], -matrix[0][2])
return phi, theta, psi
def _compute_su2_from_euler_angles(angles: tuple[float, float, float]) -> np.ndarray:
"""Computes SU(2)-matrix from Euler angles.
Args:
angles: The tuple containing the Euler angles for which the corresponding SU(2)-matrix
needs to be computed.
Returns:
The SU(2)-matrix corresponding to the Euler angles in angles.
"""
phi, theta, psi = angles
uz_phi = np.array([[np.exp(-0.5j * phi), 0], [0, np.exp(0.5j * phi)]], dtype=complex)
uy_theta = np.array(
[[math.cos(theta / 2), math.sin(theta / 2)], [-math.sin(theta / 2), math.cos(theta / 2)]],
dtype=complex,
)
ux_psi = np.array(
[[math.cos(psi / 2), math.sin(psi / 2) * 1j], [math.sin(psi / 2) * 1j, math.cos(psi / 2)]],
dtype=complex,
)
return np.dot(uz_phi, np.dot(uy_theta, ux_psi))
def _convert_su2_to_so3(matrix: np.ndarray) -> np.ndarray:
"""Computes SO(3)-matrix from input SU(2)-matrix.
Args:
matrix: The SU(2)-matrix for which a corresponding SO(3)-matrix needs to be computed.
Returns:
The SO(3)-matrix corresponding to ``matrix``.
Raises:
ValueError: if ``matrix`` is not an SU(2)-matrix.
"""
_check_is_su2(matrix)
matrix = matrix.astype(complex)
a = np.real(matrix[0][0])
b = np.imag(matrix[0][0])
c = -np.real(matrix[0][1])
d = -np.imag(matrix[0][1])
rotation = np.array(
[
[a**2 - b**2 - c**2 + d**2, 2 * a * b + 2 * c * d, -2 * a * c + 2 * b * d],
[-2 * a * b + 2 * c * d, a**2 - b**2 + c**2 - d**2, 2 * a * d + 2 * b * c],
[2 * a * c + 2 * b * d, 2 * b * c - 2 * a * d, a**2 + b**2 - c**2 - d**2],
],
dtype=float,
)
return rotation
def _convert_so3_to_su2(matrix: np.ndarray) -> np.ndarray:
"""Converts an SO(3)-matrix to a corresponding SU(2)-matrix.
Args:
matrix: SO(3)-matrix to convert.
Returns:
SU(2)-matrix corresponding to SO(3)-matrix ``matrix``.
Raises:
ValueError: if ``matrix`` is not an SO(3)-matrix.
"""
_check_is_so3(matrix)
return _compute_su2_from_euler_angles(_compute_euler_angles_from_so3(matrix))
def _check_is_su2(matrix: np.ndarray) -> None:
"""Check whether ``matrix`` is SU(2), otherwise raise an error."""
if matrix.shape != (2, 2):
raise ValueError(f"Matrix must have shape (2, 2) but has {matrix.shape}.")
if abs(np.linalg.det(matrix) - 1) > 1e-4:
raise ValueError(f"Determinant of matrix must be 1, but is {np.linalg.det(matrix)}.")
def _check_is_so3(matrix: np.ndarray) -> None:
"""Check whether ``matrix`` is SO(3), otherwise raise an error."""
if matrix.shape != (3, 3):
raise ValueError(f"Matrix must have shape (3, 3) but has {matrix.shape}.")
if abs(np.linalg.det(matrix) - 1) > 1e-4:
raise ValueError(f"Determinant of matrix must be 1, but is {np.linalg.det(matrix)}.")

View File

@ -0,0 +1,163 @@
# This code is part of Qiskit.
#
# (C) Copyright IBM 2017, 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.
"""Functions to generate the basic approximations of single qubit gates for Solovay-Kitaev."""
from __future__ import annotations
import warnings
import collections
import numpy as np
import qiskit.circuit.library.standard_gates as gates
from qiskit.circuit import Gate
from qiskit.quantum_info.operators.predicates import matrix_equal
from qiskit.utils import optionals
from .gate_sequence import GateSequence
Node = collections.namedtuple("Node", ("labels", "sequence", "children"))
_1q_inverses = {
"i": "i",
"x": "x",
"y": "y",
"z": "z",
"h": "h",
"t": "tdg",
"tdg": "t",
"s": "sdg",
"sgd": "s",
}
_1q_gates = {
"i": gates.IGate(),
"x": gates.XGate(),
"y": gates.YGate(),
"z": gates.ZGate(),
"h": gates.HGate(),
"t": gates.TGate(),
"tdg": gates.TdgGate(),
"s": gates.SGate(),
"sdg": gates.SdgGate(),
"sx": gates.SXGate(),
"sxdg": gates.SXdgGate(),
}
def _check_candidate(candidate, existing_sequences, tol=1e-10):
if optionals.HAS_SKLEARN:
return _check_candidate_kdtree(candidate, existing_sequences, tol)
warnings.warn(
"The SolovayKitaev algorithm relies on scikit-learn's KDTree for a "
"fast search over the basis approximations. Without this, we fallback onto a "
"greedy search with is significantly slower. We highly suggest to install "
"scikit-learn to use this feature.",
category=RuntimeWarning,
)
return _check_candidate_greedy(candidate, existing_sequences, tol)
def _check_candidate_greedy(candidate, existing_sequences, tol=1e-10):
# do a quick, string-based check if the same sequence already exists
if any(candidate.name == existing.name for existing in existing_sequences):
return False
for existing in existing_sequences:
if matrix_equal(existing.product_su2, candidate.product_su2, ignore_phase=True, atol=tol):
# is the new sequence less or more efficient?
return len(candidate.gates) < len(existing.gates)
return True
@optionals.HAS_SKLEARN.require_in_call
def _check_candidate_kdtree(candidate, existing_sequences, tol=1e-10):
"""Check if there's a candidate implementing the same matrix up to ``tol``.
This uses a k-d tree search and is much faster than the greedy, list-based search.
"""
from sklearn.neighbors import KDTree
# do a quick, string-based check if the same sequence already exists
if any(candidate.name == existing.name for existing in existing_sequences):
return False
points = np.array([sequence.product.flatten() for sequence in existing_sequences])
candidate = np.array([candidate.product.flatten()])
kdtree = KDTree(points)
dist, _ = kdtree.query(candidate)
return dist[0][0] > tol
def _process_node(node: Node, basis: list[str], sequences: list[GateSequence]):
inverse_last = _1q_inverses[node.labels[-1]] if node.labels else None
for label in basis:
if label == inverse_last:
continue
sequence = node.sequence.copy()
sequence.append(_1q_gates[label])
if _check_candidate(sequence, sequences):
sequences.append(sequence)
node.children.append(Node(node.labels + (label,), sequence, []))
return node.children
def generate_basic_approximations(
basis_gates: list[str | Gate], depth: int, filename: str | None = None
) -> list[GateSequence]:
"""Generates a list of ``GateSequence``s with the gates in ``basic_gates``.
Args:
basis_gates: The gates from which to create the sequences of gates.
depth: The maximum depth of the approximations.
filename: If provided, the basic approximations are stored in this file.
Returns:
List of ``GateSequences`` using the gates in ``basic_gates``.
Raises:
ValueError: If ``basis_gates`` contains an invalid gate identifier.
"""
basis = []
for gate in basis_gates:
if isinstance(gate, str):
if gate not in _1q_gates.keys():
raise ValueError(f"Invalid gate identifier: {gate}")
basis.append(gate)
else: # gate is a qiskit.circuit.Gate
basis.append(gate.name)
tree = Node((), GateSequence(), [])
cur_level = [tree]
sequences = [tree.sequence]
for _ in [None] * depth:
next_level = []
for node in cur_level:
next_level.extend(_process_node(node, basis, sequences))
cur_level = next_level
if filename is not None:
data = {}
for sequence in sequences:
gatestring = sequence.name
data[gatestring] = sequence.product
np.save(filename, data)
return sequences

View File

@ -0,0 +1,198 @@
# This code is part of Qiskit.
#
# (C) Copyright IBM 2017, 2020.
#
# 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.
"""Synthesize a single qubit gate to a discrete basis set."""
from __future__ import annotations
import numpy as np
from qiskit.circuit import QuantumCircuit, Gate
from qiskit.dagcircuit import DAGCircuit
from .gate_sequence import GateSequence
from .commutator_decompose import commutator_decompose
from .generate_basis_approximations import generate_basic_approximations, _1q_gates, _1q_inverses
class SolovayKitaevDecomposition:
"""The Solovay Kitaev discrete decomposition algorithm.
This class is called recursively by the transpiler pass, which is why it is separeted.
See :mod:`qiskit.synthesis.discrete_basis` for more information.
"""
def __init__(
self, basic_approximations: str | dict[str, np.ndarray] | list[GateSequence] | None = None
) -> None:
"""
Args:
basic_approximations: A specification of the basic SU(2) approximations in terms
of discrete gates. At each iteration this algorithm, the remaining error is
approximated with the closest sequence of gates in this set.
If a ``str``, this specifies a ``.npy`` filename from which to load the
approximation. If a ``dict``, then this contains
``{gates: effective_SO3_matrix}`` pairs,
e.g. ``{"h t": np.array([[0, 0.7071, -0.7071], [0, -0.7071, -0.7071], [-1, 0, 0]]}``.
If a list, this contains the same information as the dict, but already converted to
:class:`.GateSequence` objects, which contain the SO(3) matrix and gates.
"""
if basic_approximations is None:
# generate a default basic approximation
basic_approximations = generate_basic_approximations(
basis_gates=["h", "t", "tdg"], depth=10
)
self.basic_approximations = self.load_basic_approximations(basic_approximations)
def load_basic_approximations(self, data: list | str | dict) -> list[GateSequence]:
"""Load basic approximations.
Args:
data: If a string, specifies the path to the file from where to load the data.
If a dictionary, directly specifies the decompositions as ``{gates: matrix}``.
There ``gates`` are the names of the gates producing the SO(3) matrix ``matrix``,
e.g. ``{"h t": np.array([[0, 0.7071, -0.7071], [0, -0.7071, -0.7071], [-1, 0, 0]]}``.
Returns:
A list of basic approximations as type ``GateSequence``.
Raises:
ValueError: If the number of gate combinations and associated matrices does not match.
"""
# is already a list of GateSequences
if isinstance(data, list):
return data
# if a file, load the dictionary
if isinstance(data, str):
data = np.load(data, allow_pickle=True)
sequences = []
for gatestring, matrix in data.items():
sequence = GateSequence()
sequence.gates = [_1q_gates[element] for element in gatestring.split()]
sequence.product = np.asarray(matrix)
sequences.append(sequence)
return sequences
def run(
self, gate_matrix: np.ndarray, recursion_degree: int, return_dag: bool = False
) -> QuantumCircuit | DAGCircuit:
r"""Run the algorithm.
Args:
gate_matrix: The 2x2 matrix representing the gate. Does not need to be SU(2).
recursion_degree: The recursion degree, called :math:`n` in the paper.
return_dag: If ``True`` return a :class:`.DAGCircuit`, else a :class:`.QuantumCircuit`.
Returns:
A one-qubit circuit approximating the ``gate_matrix`` in the specified discrete basis.
"""
# make input matrix SU(2) and get the according global phase
z = 1 / np.sqrt(np.linalg.det(gate_matrix))
gate_matrix_su2 = GateSequence.from_matrix(z * gate_matrix)
global_phase = np.arctan2(np.imag(z), np.real(z))
# get the decompositon as GateSequence type
decomposition = self._recurse(gate_matrix_su2, recursion_degree)
# simplify
_remove_identities(decomposition)
_remove_inverse_follows_gate(decomposition)
# convert to a circuit and attach the right phases
if return_dag:
out = decomposition.to_dag()
else:
out = decomposition.to_circuit()
out.global_phase = decomposition.global_phase - global_phase
return out
def _recurse(self, sequence: GateSequence, n: int) -> GateSequence:
"""Performs ``n`` iterations of the Solovay-Kitaev algorithm on ``sequence``.
Args:
sequence: GateSequence to which the Solovay-Kitaev algorithm is applied.
n: number of iterations that the algorithm needs to run.
Returns:
GateSequence that approximates ``sequence``.
Raises:
ValueError: if ``u`` does not represent an SO(3)-matrix.
"""
if sequence.product.shape != (3, 3):
raise ValueError("Shape of U must be (3, 3) but is", sequence.shape)
if n == 0:
return self.find_basic_approximation(sequence)
u_n1 = self._recurse(sequence, n - 1)
v_n, w_n = commutator_decompose(sequence.dot(u_n1.adjoint()).product)
v_n1 = self._recurse(v_n, n - 1)
w_n1 = self._recurse(w_n, n - 1)
return v_n1.dot(w_n1).dot(v_n1.adjoint()).dot(w_n1.adjoint()).dot(u_n1)
def find_basic_approximation(self, sequence: GateSequence) -> Gate:
"""Finds gate in ``self._basic_approximations`` that best represents ``sequence``.
Args:
sequence: The gate to find the approximation to.
Returns:
Gate in basic approximations that is closest to ``sequence``.
"""
# TODO explore using a k-d tree here
def key(x):
return np.linalg.norm(np.subtract(x.product, sequence.product))
best = min(self.basic_approximations, key=key)
return best
def _remove_inverse_follows_gate(sequence):
index = 0
while index < len(sequence.gates) - 1:
curr_gate = sequence.gates[index]
next_gate = sequence.gates[index + 1]
if curr_gate.name in _1q_inverses.keys():
remove = _1q_inverses[curr_gate.name] == next_gate.name
else:
remove = curr_gate.inverse() == next_gate
if remove:
# remove gates at index and index + 1
sequence.remove_cancelling_pair([index, index + 1])
# take a step back to see if we have uncovered a new pair, e.g.
# [h, s, sdg, h] at index = 1 removes s, sdg but if we continue at index 1
# we miss the uncovered [h, h] pair at indices 0 and 1
if index > 0:
index -= 1
else:
# next index
index += 1
def _remove_identities(sequence):
index = 0
while index < len(sequence.gates):
if sequence.gates[index].name == "id":
sequence.gates.pop(index)
else:
index += 1

View File

@ -134,7 +134,7 @@ Circuit Analysis
DAGLongestPath
Synthesis
=============
=========
.. autosummary::
:toctree: ../stubs/
@ -143,6 +143,8 @@ Synthesis
LinearFunctionsSynthesis
LinearFunctionsToPermutations
HighLevelSynthesis
SolovayKitaev
SolovayKitaevSynthesis
Post Layout (Post transpile qubit selection)
============================================
@ -245,6 +247,8 @@ from .synthesis import unitary_synthesis_plugin_names
from .synthesis import LinearFunctionsSynthesis
from .synthesis import LinearFunctionsToPermutations
from .synthesis import HighLevelSynthesis
from .synthesis import SolovayKitaev
from .synthesis import SolovayKitaevSynthesis
# calibration
from .calibration import PulseGates

View File

@ -24,7 +24,7 @@ class Collect1qRuns(AnalysisPass):
The blocks contain "op" nodes in topological order such that all gates
in a block act on the same qubits and are adjacent in the circuit.
After the execution, ``property_set['block_list']`` is set to a list of
After the execution, ``property_set['run_list']`` is set to a list of
tuples of "op" node.
"""
self.property_set["run_list"] = dag.collect_1q_runs()

View File

@ -73,7 +73,7 @@ class ConsolidateBlocks(TransformationPass):
# compute ordered indices for the global circuit wires
global_index_map = {wire: idx for idx, wire in enumerate(dag.qubits)}
blocks = self.property_set["block_list"]
blocks = self.property_set["block_list"] or []
basis_gate_name = self.decomposer.gate.name
all_block_gates = set()
for block in blocks:

View File

@ -16,3 +16,4 @@ from .unitary_synthesis import UnitarySynthesis
from .plugin import unitary_synthesis_plugin_names
from .linear_functions_synthesis import LinearFunctionsSynthesis, LinearFunctionsToPermutations
from .high_level_synthesis import HighLevelSynthesis, HLSConfig
from .solovay_kitaev_synthesis import SolovayKitaev, SolovayKitaevSynthesis

View File

@ -0,0 +1,196 @@
# This code is part of Qiskit.
#
# (C) Copyright IBM 2021.
#
# 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 Solovay-Kitaev synthesis plugin to Qiskit's transpiler.
"""
from __future__ import annotations
import numpy as np
from qiskit.converters import circuit_to_dag
from qiskit.dagcircuit import DAGCircuit
from qiskit.synthesis.discrete_basis.solovay_kitaev import SolovayKitaevDecomposition
from qiskit.synthesis.discrete_basis.generate_basis_approximations import (
generate_basic_approximations,
)
from qiskit.transpiler.basepasses import TransformationPass
from qiskit.transpiler.exceptions import TranspilerError
from .plugin import UnitarySynthesisPlugin
class SolovayKitaev(TransformationPass):
r"""Approximately decompose 1q gates to a discrete basis using the Solovay-Kitaev algorithm.
See :mod:`qiskit.synthesis.discrete_basis` for more information.
"""
def __init__(
self,
recursion_degree: int = 3,
basic_approximations: str | dict[str, np.ndarray] | None = None,
) -> None:
"""
Args:
recursion_degree: The recursion depth for the Solovay-Kitaev algorithm.
A larger recursion depth increases the accuracy and length of the
decomposition.
basic_approximations: The basic approximations for the finding the best discrete
decomposition at the root of the recursion. If a string, it specifies the ``.npy``
file to load the approximations from. If a dictionary, it contains
``{label: SO(3)-matrix}`` pairs. If None, a default based on the H, T and Tdg gates
up to combinations of depth 10 is generated.
"""
super().__init__()
self.recursion_degree = recursion_degree
self._sk = SolovayKitaevDecomposition(basic_approximations)
def run(self, dag: DAGCircuit) -> DAGCircuit:
"""Run the ``SolovayKitaev`` pass on `dag`.
Args:
dag: The input dag.
Returns:
Output dag with 1q gates synthesized in the discrete target basis.
Raises:
TranspilerError: if a gates does not have to_matrix
"""
for node in dag.op_nodes():
if not node.op.num_qubits == 1:
continue # ignore all non-single qubit gates
if not hasattr(node.op, "to_matrix"):
raise TranspilerError(
"SolovayKitaev does not support gate without "
f"to_matrix method: {node.op.name}"
)
matrix = node.op.to_matrix()
# call solovay kitaev
approximation = self._sk.run(matrix, self.recursion_degree, return_dag=True)
# convert to a dag and replace the gate by the approximation
dag.substitute_node_with_dag(node, approximation)
return dag
class SolovayKitaevSynthesis(UnitarySynthesisPlugin):
"""A Solovay-Kitaev Qiskit unitary synthesis plugin.
This plugin is invoked by :func:`~.compiler.transpile` when the ``unitary_synthesis_method``
parameter is set to ``"sk"``.
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:
basis_approximations (str | dict):
The basic approximations for the finding the best discrete decomposition at the root of the
recursion. If a string, it specifies the ``.npy`` file to load the approximations from.
If a dictionary, it contains ``{label: SO(3)-matrix}`` pairs. If None, a default based on
the specified ``basis_gates`` and ``depth`` is generated.
basis_gates (list):
A list of strings specifying the discrete basis gates to decompose to. If None,
defaults to ``["h", "t", "tdg"]``.
depth (int):
The gate-depth of the the basic approximations. All possible, unique combinations of the
basis gates up to length ``depth`` are considered. If None, defaults to 10.
recursion_degree (int):
The number of times the decomposition is recursively improved. If None, defaults to 3.
"""
# we cache an instance of the Solovay-Kitaev class to generate the
# computationally expensive basis approximation of single qubit gates only once
_sk = None
@property
def max_qubits(self):
"""Maximum number of supported qubits is ``1``."""
return 1
@property
def min_qubits(self):
"""Minimum number of supported qubits is ``1``."""
return 1
@property
def supports_natural_direction(self):
"""The plugin does not support natural direction, it does not assume
bidirectional two qubit gates."""
return True
@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. By default it synthesis to the
``["h", "t", "tdg"]`` gate basis."""
return True
@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
config = options.get("config") or {}
recursion_degree = config.get("recursion_degree", 3)
# if we didn't yet construct the Solovay-Kitaev instance, which contains
# the basic approximations, do it now
if SolovayKitaevSynthesis._sk is None:
basic_approximations = config.get("basic_approximations", None)
basis_gates = options.get("basis_gates", ["h", "t", "tdg"])
# if the basic approximations are not generated and not given,
# try to generate them if the basis set is specified
if basic_approximations is None:
depth = config.get("depth", 10)
basic_approximations = generate_basic_approximations(basis_gates, depth)
SolovayKitaevSynthesis._sk = SolovayKitaevDecomposition(basic_approximations)
approximate_circuit = SolovayKitaevSynthesis._sk.run(unitary, recursion_degree)
dag_circuit = circuit_to_dag(approximate_circuit)
return dag_circuit

View File

@ -0,0 +1,46 @@
---
features:
- |
Added the :class:`.SolovayKitaev` transpiler pass to run the Solovay-Kitaev algorithm for
approximating single-qubit unitaries using a discrete gate set. In combination with the basis
translator, this allows to convert any unitary circuit to a universal discrete gate set,
which could be implemented fault-tolerantly.
This pass can e.g. be used after compiling to U and CX gates:
.. jupyter-execute::
from qiskit import transpile
from qiskit.circuit.library import QFT
from qiskit.transpiler.passes.synthesis import SolovayKitaev
qft = QFT(3)
# optimize to general 1-qubit unitaries and CX
transpiled = transpile(qft, basis_gates=["u", "cx"], optimization_level=1)
skd = SolovayKitaev() # uses T Tdg and H as default basis
discretized = skd(transpiled)
print(discretized.count_ops())
The decomposition can also be used with the unitary synthesis plugin, as
the "sk" method on the :class:`~.UnitarySynthesis` transpiler pass:
.. jupyter-execute::
from qiskit import QuantumCircuit
from qiskit.quantum_info import Operator
from qiskit.transpiler.passes import UnitarySynthesis
circuit = QuantumCircuit(1)
circuit.rx(0.8, 0)
unitary = Operator(circuit).data
unitary_circ = QuantumCircuit(1)
unitary_circ.unitary(unitary, [0])
synth = UnitarySynthesis(basis_gates=["h", "s"], method="sk")
out = synth(unitary_circ)
out.draw('mpl')

View File

@ -101,6 +101,7 @@ setup(
"qiskit.unitary_synthesis": [
"default = qiskit.transpiler.passes.synthesis.unitary_synthesis:DefaultUnitarySynthesis",
"aqc = qiskit.transpiler.synthesis.aqc.aqc_plugin:AQCSynthesisPlugin",
"sk = qiskit.transpiler.passes.synthesis.solovay_kitaev_synthesis:SolovayKitaevSynthesis",
],
"qiskit.synthesis": [
"clifford.default = qiskit.transpiler.passes.synthesis.high_level_synthesis:DefaultSynthesisClifford",

View File

@ -0,0 +1,370 @@
# This code is part of Qiskit.
#
# (C) Copyright IBM 2017, 2020.
#
# 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.
"""Test the Solovay Kitaev transpilation pass."""
import unittest
import math
import numpy as np
import scipy
from ddt import ddt, data
from qiskit.test import QiskitTestCase
from qiskit import transpile
from qiskit.circuit import QuantumCircuit
from qiskit.circuit.library import TGate, TdgGate, HGate, SGate, SdgGate, IGate, QFT
from qiskit.converters import circuit_to_dag, dag_to_circuit
from qiskit.quantum_info import Operator
from qiskit.synthesis.discrete_basis.generate_basis_approximations import (
generate_basic_approximations,
)
from qiskit.synthesis.discrete_basis.commutator_decompose import commutator_decompose
from qiskit.synthesis.discrete_basis.gate_sequence import GateSequence
from qiskit.transpiler import PassManager
from qiskit.transpiler.exceptions import TranspilerError
from qiskit.transpiler.passes import UnitarySynthesis, Collect1qRuns, ConsolidateBlocks
from qiskit.transpiler.passes.synthesis import SolovayKitaev, SolovayKitaevSynthesis
def _trace_distance(circuit1, circuit2):
"""Return the trace distance of the two input circuits."""
op1, op2 = Operator(circuit1), Operator(circuit2)
return 0.5 * np.trace(scipy.linalg.sqrtm(np.conj(op1 - op2).T.dot(op1 - op2))).real
def _generate_x_rotation(angle: float) -> np.ndarray:
return np.array(
[[1, 0, 0], [0, math.cos(angle), -math.sin(angle)], [0, math.sin(angle), math.cos(angle)]]
)
def _generate_y_rotation(angle: float) -> np.ndarray:
return np.array(
[[math.cos(angle), 0, math.sin(angle)], [0, 1, 0], [-math.sin(angle), 0, math.cos(angle)]]
)
def _generate_z_rotation(angle: float) -> np.ndarray:
return np.array(
[[math.cos(angle), -math.sin(angle), 0], [math.sin(angle), math.cos(angle), 0], [0, 0, 1]]
)
def is_so3_matrix(array: np.ndarray) -> bool:
"""Check if the input array is a SO(3) matrix."""
if array.shape != (3, 3):
return False
if abs(np.linalg.det(array) - 1.0) > 1e-10:
return False
if False in np.isreal(array):
return False
return True
@ddt
class TestSolovayKitaev(QiskitTestCase):
"""Test the Solovay Kitaev algorithm and transformation pass."""
def setUp(self):
super().setUp()
self.basic_approx = generate_basic_approximations([HGate(), TGate(), TdgGate()], 3)
def test_unitary_synthesis(self):
"""Test the unitary synthesis transpiler pass with Solovay-Kitaev."""
circuit = QuantumCircuit(2)
circuit.rx(0.8, 0)
circuit.cx(0, 1)
circuit.x(1)
_1q = Collect1qRuns()
_cons = ConsolidateBlocks()
_synth = UnitarySynthesis(["h", "s"], method="sk")
passes = PassManager([_1q, _cons, _synth])
compiled = passes.run(circuit)
diff = np.linalg.norm(Operator(compiled) - Operator(circuit))
self.assertLess(diff, 1)
self.assertEqual(set(compiled.count_ops().keys()), {"h", "s", "cx"})
def test_plugin(self):
"""Test calling the plugin directly."""
circuit = QuantumCircuit(1)
circuit.rx(0.8, 0)
unitary = Operator(circuit).data
plugin = SolovayKitaevSynthesis()
out = plugin.run(unitary, basis_gates=["h", "s"])
reference = QuantumCircuit(1, global_phase=3 * np.pi / 4)
reference.h(0)
reference.s(0)
reference.h(0)
self.assertEqual(dag_to_circuit(out), reference)
def test_generating_default_approximation(self):
"""Test the approximation set is generated by default."""
skd = SolovayKitaev()
circuit = QuantumCircuit(1)
dummy = skd(circuit)
self.assertIsNotNone(skd._sk.basic_approximations)
def test_i_returns_empty_circuit(self):
"""Test that ``SolovayKitaev`` returns an empty circuit when
it approximates the I-gate."""
circuit = QuantumCircuit(1)
circuit.i(0)
skd = SolovayKitaev(3, self.basic_approx)
decomposed_circuit = skd(circuit)
self.assertEqual(QuantumCircuit(1), decomposed_circuit)
def test_exact_decomposition_acts_trivially(self):
"""Test that the a circuit that can be represented exactly is represented exactly."""
circuit = QuantumCircuit(1)
circuit.t(0)
circuit.h(0)
circuit.tdg(0)
synth = SolovayKitaev(3, self.basic_approx)
dag = circuit_to_dag(circuit)
decomposed_dag = synth.run(dag)
decomposed_circuit = dag_to_circuit(decomposed_dag)
self.assertEqual(circuit, decomposed_circuit)
def test_fails_with_no_to_matrix(self):
"""Test failer if gate does not have to_matrix."""
circuit = QuantumCircuit(1)
circuit.initialize("0")
synth = SolovayKitaev(3, self.basic_approx)
dag = circuit_to_dag(circuit)
with self.assertRaises(TranspilerError) as cm:
_ = synth.run(dag)
self.assertEqual(
"SolovayKitaev does not support gate without to_matrix method: initialize",
cm.exception.message,
)
def test_str_basis_gates(self):
"""Test specifying the basis gates by string works."""
circuit = QuantumCircuit(1)
circuit.rx(0.8, 0)
basic_approx = generate_basic_approximations(["h", "t", "s"], 3)
synth = SolovayKitaev(2, basic_approx)
dag = circuit_to_dag(circuit)
discretized = dag_to_circuit(synth.run(dag))
reference = QuantumCircuit(1, global_phase=7 * np.pi / 8)
reference.h(0)
reference.t(0)
reference.h(0)
self.assertEqual(discretized, reference)
def test_approximation_on_qft(self):
"""Test the Solovay-Kitaev decomposition on the QFT circuit."""
qft = QFT(3)
transpiled = transpile(qft, basis_gates=["u", "cx"], optimization_level=1)
skd = SolovayKitaev(1)
with self.subTest("1 recursion"):
discretized = skd(transpiled)
self.assertLess(_trace_distance(transpiled, discretized), 15)
skd.recursion_degree = 2
with self.subTest("2 recursions"):
discretized = skd(transpiled)
self.assertLess(_trace_distance(transpiled, discretized), 7)
@ddt
class TestGateSequence(QiskitTestCase):
"""Test the ``GateSequence`` class."""
def test_append(self):
"""Test append."""
seq = GateSequence([IGate()])
seq.append(HGate())
ref = GateSequence([IGate(), HGate()])
self.assertEqual(seq, ref)
def test_eq(self):
"""Test equality."""
base = GateSequence([HGate(), HGate()])
seq1 = GateSequence([HGate(), HGate()])
seq2 = GateSequence([IGate()])
seq3 = GateSequence([HGate(), HGate()])
seq3.global_phase = 0.12
seq4 = GateSequence([IGate(), HGate()])
with self.subTest("equal"):
self.assertEqual(base, seq1)
with self.subTest("same product, but different repr (-> false)"):
self.assertNotEqual(base, seq2)
with self.subTest("differing global phase (-> false)"):
self.assertNotEqual(base, seq3)
with self.subTest("same num gates, but different gates (-> false)"):
self.assertNotEqual(base, seq4)
def test_to_circuit(self):
"""Test converting a gate sequence to a circuit."""
seq = GateSequence([HGate(), HGate(), TGate(), SGate(), SdgGate()])
ref = QuantumCircuit(1)
ref.h(0)
ref.h(0)
ref.t(0)
ref.s(0)
ref.sdg(0)
# a GateSequence is SU(2), so add the right phase
z = 1 / np.sqrt(np.linalg.det(Operator(ref)))
ref.global_phase = np.arctan2(np.imag(z), np.real(z))
self.assertEqual(seq.to_circuit(), ref)
def test_adjoint(self):
"""Test adjoint."""
seq = GateSequence([TGate(), SGate(), HGate(), IGate()])
inv = GateSequence([IGate(), HGate(), SdgGate(), TdgGate()])
self.assertEqual(seq.adjoint(), inv)
def test_copy(self):
"""Test copy."""
seq = GateSequence([IGate()])
copied = seq.copy()
seq.gates.append(HGate())
self.assertEqual(len(seq.gates), 2)
self.assertEqual(len(copied.gates), 1)
@data(0, 1, 10)
def test_len(self, n):
"""Test __len__."""
seq = GateSequence([IGate()] * n)
self.assertEqual(len(seq), n)
def test_getitem(self):
"""Test __getitem__."""
seq = GateSequence([IGate(), HGate(), IGate()])
self.assertEqual(seq[0], IGate())
self.assertEqual(seq[1], HGate())
self.assertEqual(seq[2], IGate())
self.assertEqual(seq[-2], HGate())
def test_from_su2_matrix(self):
"""Test from_matrix with an SU2 matrix."""
matrix = np.array([[1, 1], [1, -1]], dtype=complex) / np.sqrt(2)
matrix /= np.sqrt(np.linalg.det(matrix))
seq = GateSequence.from_matrix(matrix)
ref = GateSequence([HGate()])
self.assertEqual(seq.gates, list())
self.assertTrue(np.allclose(seq.product, ref.product))
self.assertEqual(seq.global_phase, 0)
def test_from_so3_matrix(self):
"""Test from_matrix with an SO3 matrix."""
matrix = np.array([[0, 0, -1], [0, -1, 0], [-1, 0, 0]])
seq = GateSequence.from_matrix(matrix)
ref = GateSequence([HGate()])
self.assertEqual(seq.gates, list())
self.assertTrue(np.allclose(seq.product, ref.product))
self.assertEqual(seq.global_phase, 0)
def test_from_invalid_matrix(self):
"""Test from_matrix with invalid matrices."""
with self.subTest("2x2 but not SU2"):
matrix = np.array([[1, 1], [1, -1]], dtype=complex) / np.sqrt(2)
with self.assertRaises(ValueError):
_ = GateSequence.from_matrix(matrix)
with self.subTest("not 2x2 or 3x3"):
with self.assertRaises(ValueError):
_ = GateSequence.from_matrix(np.array([[1]]))
def test_dot(self):
"""Test dot."""
seq1 = GateSequence([HGate()])
seq2 = GateSequence([TGate(), SGate()])
composed = seq1.dot(seq2)
ref = GateSequence([TGate(), SGate(), HGate()])
# check the product matches
self.assertTrue(np.allclose(ref.product, composed.product))
# check the circuit & phases are equivalent
self.assertTrue(Operator(ref.to_circuit()).equiv(composed.to_circuit()))
@ddt
class TestSolovayKitaevUtils(QiskitTestCase):
"""Test the public functions in the Solovay Kitaev utils."""
@data(
_generate_x_rotation(0.1),
_generate_y_rotation(0.2),
_generate_z_rotation(0.3),
np.dot(_generate_z_rotation(0.5), _generate_y_rotation(0.4)),
np.dot(_generate_y_rotation(0.5), _generate_x_rotation(0.4)),
)
def test_commutator_decompose_return_type(self, u_so3: np.ndarray):
"""Test that ``commutator_decompose`` returns two SO(3) gate sequences."""
v, w = commutator_decompose(u_so3)
self.assertTrue(is_so3_matrix(v.product))
self.assertTrue(is_so3_matrix(w.product))
self.assertIsInstance(v, GateSequence)
self.assertIsInstance(w, GateSequence)
@data(
_generate_x_rotation(0.1),
_generate_y_rotation(0.2),
_generate_z_rotation(0.3),
np.dot(_generate_z_rotation(0.5), _generate_y_rotation(0.4)),
np.dot(_generate_y_rotation(0.5), _generate_x_rotation(0.4)),
)
def test_commutator_decompose_decomposes_correctly(self, u_so3):
"""Test that ``commutator_decompose`` exactly decomposes the input."""
v, w = commutator_decompose(u_so3)
v_so3 = v.product
w_so3 = w.product
actual_commutator = np.dot(v_so3, np.dot(w_so3, np.dot(np.conj(v_so3).T, np.conj(w_so3).T)))
self.assertTrue(np.allclose(actual_commutator, u_so3))
if __name__ == "__main__":
unittest.main()