diff --git a/qiskit/quantum_info/synthesis/xx_decompose/__init__.py b/qiskit/quantum_info/synthesis/xx_decompose/__init__.py new file mode 100644 index 000000000000..0b199e8946b3 --- /dev/null +++ b/qiskit/quantum_info/synthesis/xx_decompose/__init__.py @@ -0,0 +1,17 @@ +# 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. + +""" +Submodule symbol exports for XX decomposition. +""" + +from .decomposer import XXDecomposer diff --git a/qiskit/quantum_info/synthesis/xx_decompose/circuits.py b/qiskit/quantum_info/synthesis/xx_decompose/circuits.py new file mode 100644 index 000000000000..8abf8dfe3489 --- /dev/null +++ b/qiskit/quantum_info/synthesis/xx_decompose/circuits.py @@ -0,0 +1,298 @@ +# 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. + +""" +Tools for building optimal circuits out of XX interactions. + +Inputs: + + A set of native XX operations, described as strengths. + + A right-angled path, computed using the methods in `paths.py`. + +Output: + + A circuit which implements the target operation (expressed exactly as the exponential of + `a XX + b YY + c ZZ`) using the native operations and local gates. +""" + +from functools import reduce +import math +from operator import itemgetter + +import numpy as np + +from qiskit.circuit.quantumcircuit import QuantumCircuit +from qiskit.circuit.library.standard_gates import RXXGate, RYYGate, RZGate +from qiskit.exceptions import QiskitError + +from .paths import decomposition_hop +from .utilities import EPSILON, safe_arccos +from .weyl import ( + apply_reflection, + apply_shift, + canonical_rotation_circuit, + reflection_options, + shift_options, +) + + +# pylint:disable=invalid-name +def decompose_xxyy_into_xxyy_xx(a_target, b_target, a_source, b_source, interaction): + """ + Consumes a target canonical interaction CAN(a_target, b_target) and source interactions + CAN(a1, b1), CAN(a2), then manufactures a circuit identity of the form + + CAN(a_target, b_target) = (Zr, Zs) CAN(a_source, b_source) (Zu, Zv) CAN(interaction) (Zx, Zy). + + Returns the 6-tuple (r, s, u, v, x, y). + """ + + cplus, cminus = np.cos(a_source + b_source), np.cos(a_source - b_source) + splus, sminus = np.sin(a_source + b_source), np.sin(a_source - b_source) + ca, sa = np.cos(interaction), np.sin(interaction) + + uplusv = ( + 1 + / 2 + * safe_arccos( + cminus ** 2 * ca ** 2 + sminus ** 2 * sa ** 2 - np.cos(a_target - b_target) ** 2, + 2 * cminus * ca * sminus * sa, + ) + ) + uminusv = ( + 1 + / 2 + * safe_arccos( + cplus ** 2 * ca ** 2 + splus ** 2 * sa ** 2 - np.cos(a_target + b_target) ** 2, + 2 * cplus * ca * splus * sa, + ) + ) + + u, v = (uplusv + uminusv) / 2, (uplusv - uminusv) / 2 + + # NOTE: the target matrix is phase-free + middle_matrix = reduce( + np.dot, + [ + RXXGate(2 * a_source).to_matrix() @ RYYGate(2 * b_source).to_matrix(), + np.kron(RZGate(2 * u).to_matrix(), RZGate(2 * v).to_matrix()), + RXXGate(2 * interaction).to_matrix(), + ], + ) + + phase_solver = np.array( + [ + [ + 1 / 4, + 1 / 4, + 1 / 4, + 1 / 4, + ], + [ + 1 / 4, + -1 / 4, + -1 / 4, + 1 / 4, + ], + [ + 1 / 4, + 1 / 4, + -1 / 4, + -1 / 4, + ], + [ + 1 / 4, + -1 / 4, + 1 / 4, + -1 / 4, + ], + ] + ) + inner_phases = [ + np.angle(middle_matrix[0, 0]), + np.angle(middle_matrix[1, 1]), + np.angle(middle_matrix[1, 2]) + np.pi / 2, + np.angle(middle_matrix[0, 3]) + np.pi / 2, + ] + r, s, x, y = np.dot(phase_solver, inner_phases) + + # If there's a phase discrepancy, need to conjugate by an extra Z/2 (x) Z/2. + generated_matrix = reduce( + np.dot, + [ + np.kron(RZGate(2 * r).to_matrix(), RZGate(2 * s).to_matrix()), + middle_matrix, + np.kron(RZGate(2 * x).to_matrix(), RZGate(2 * y).to_matrix()), + ], + ) + if (abs(np.angle(generated_matrix[3, 0]) - np.pi / 2) < 0.01 and a_target > b_target) or ( + abs(np.angle(generated_matrix[3, 0]) + np.pi / 2) < 0.01 and a_target < b_target + ): + x += np.pi / 4 + y += np.pi / 4 + r -= np.pi / 4 + s -= np.pi / 4 + + return r, s, u, v, x, y + + +def xx_circuit_step(source, strength, target, embodiment): + """ + Builds a single step in an XX-based circuit. + + `source` and `target` are positive canonical coordinates; `strength` is the interaction strength + at this step in the circuit as a canonical coordinate (so that CX = RZX(pi/2) corresponds to + pi/4); and `embodiment` is a Qiskit circuit which enacts the canonical gate of the prescribed + interaction `strength`. + """ + + permute_source_for_overlap, permute_target_for_overlap = None, None + + # apply all possible reflections, shifts to the source + for source_reflection_name in reflection_options: + reflected_source_coord, source_reflection, reflection_phase_shift = apply_reflection( + source_reflection_name, source + ) + for source_shift_name in shift_options: + shifted_source_coord, source_shift, shift_phase_shift = apply_shift( + source_shift_name, reflected_source_coord + ) + + # check for overlap, back out permutation + source_shared, target_shared = None, None + for i, j in [(0, 0), (0, 1), (0, 2), (1, 0), (1, 1), (1, 2), (2, 0), (2, 1), (2, 2)]: + + if ( + abs(np.mod(abs(shifted_source_coord[i] - target[j]), np.pi)) < EPSILON + or abs(np.mod(abs(shifted_source_coord[i] - target[j]), np.pi) - np.pi) + < EPSILON + ): + source_shared, target_shared = i, j + break + if source_shared is None: + continue + + # pick out the other coordinates + source_first, source_second = [x for x in [0, 1, 2] if x != source_shared] + target_first, target_second = [x for x in [0, 1, 2] if x != target_shared] + + # check for arccos validity + r, s, u, v, x, y = decompose_xxyy_into_xxyy_xx( + float(target[target_first]), + float(target[target_second]), + float(shifted_source_coord[source_first]), + float(shifted_source_coord[source_second]), + float(strength), + ) + if any(math.isnan(val) for val in (r, s, u, v, x, y)): + continue + + # OK: this combination of things works. + # save the permutation which rotates the shared coordinate into ZZ. + permute_source_for_overlap = canonical_rotation_circuit(source_first, source_second) + permute_target_for_overlap = canonical_rotation_circuit(target_first, target_second) + break + + if permute_source_for_overlap is not None: + break + + if permute_source_for_overlap is None: + raise QiskitError( + f"Error during RZX decomposition: Could not find a suitable Weyl " + f"reflection to match {source} to {target} along {strength}." + ) + + prefix_circuit, affix_circuit = QuantumCircuit(2), QuantumCircuit(2) + + # the basic formula we're trying to work with is: + # target^p_t_f_o = + # rs * (source^s_reflection * s_shift)^p_s_f_o * uv * operation * xy + # but we're rearranging it into the form + # target = affix source prefix + # and computing just the prefix / affix circuits. + + # the outermost prefix layer comes from the (inverse) target permutation. + prefix_circuit.compose(permute_target_for_overlap.inverse(), inplace=True) + # the middle prefix layer comes from the local Z rolls. + prefix_circuit.rz(2 * x, [0]) + prefix_circuit.rz(2 * y, [1]) + prefix_circuit.compose(embodiment, inplace=True) + prefix_circuit.rz(2 * u, [0]) + prefix_circuit.rz(2 * v, [1]) + # the innermost prefix layer is source_reflection, shifted by source_shift, + # finally conjugated by p_s_f_o. + prefix_circuit.compose(permute_source_for_overlap, inplace=True) + prefix_circuit.compose(source_reflection, inplace=True) + prefix_circuit.global_phase += -np.log(reflection_phase_shift).imag + prefix_circuit.global_phase += -np.log(shift_phase_shift).imag + + # the affix circuit is constructed in reverse. + # first (i.e., innermost), we install the other half of the source transformations and p_s_f_o. + affix_circuit.compose(source_reflection.inverse(), inplace=True) + affix_circuit.compose(source_shift, inplace=True) + affix_circuit.compose(permute_source_for_overlap.inverse(), inplace=True) + # then, the other local rolls in the middle. + affix_circuit.rz(2 * r, [0]) + affix_circuit.rz(2 * s, [1]) + # finally, the other half of the p_t_f_o conjugation. + affix_circuit.compose(permute_target_for_overlap, inplace=True) + + return {"prefix_circuit": prefix_circuit, "affix_circuit": affix_circuit} + + +def canonical_xx_circuit(target, strength_sequence, basis_embodiments): + """ + Assembles a Qiskit circuit from a specified `strength_sequence` of XX-type interactions which + emulates the canonical gate at canonical coordinate `target`. The circuits supplied by + `basis_embodiments` are used to instantiate the individual XX actions. + + NOTE: The elements of `strength_sequence` are expected to be normalized so that np.pi/2 + corresponds to RZX(np.pi/2) = CX; `target` is taken to be a positive canonical coordinate; + and `basis_embodiments` maps `strength_sequence` elements to circuits which instantiate + these gates. + """ + # empty decompositions are easy! + if len(strength_sequence) == 0: + return QuantumCircuit(2) + + # assemble the prefix / affix circuits + prefix_circuit, affix_circuit = QuantumCircuit(2), QuantumCircuit(2) + while len(strength_sequence) > 1: + source = decomposition_hop(target, strength_sequence) + strength = strength_sequence[-1] + + preceding_prefix_circuit, preceding_affix_circuit = itemgetter( + "prefix_circuit", "affix_circuit" + )(xx_circuit_step(source, strength / 2, target, basis_embodiments[strength])) + + prefix_circuit.compose(preceding_prefix_circuit, inplace=True) + affix_circuit.compose(preceding_affix_circuit, inplace=True, front=True) + + target, strength_sequence = source, strength_sequence[:-1] + + circuit = prefix_circuit + + # lastly, deal with the "leading" gate. + if target[0] <= np.pi / 4: + circuit.compose(basis_embodiments[strength_sequence[0]], inplace=True) + else: + _, source_reflection, reflection_phase_shift = apply_reflection("reflect XX, YY", [0, 0, 0]) + _, source_shift, shift_phase_shift = apply_shift("X shift", [0, 0, 0]) + + circuit.compose(source_reflection, inplace=True) + circuit.compose(basis_embodiments[strength_sequence[0]], inplace=True) + circuit.compose(source_reflection.inverse(), inplace=True) + circuit.compose(source_shift, inplace=True) + circuit.global_phase += -np.log(shift_phase_shift).imag + circuit.global_phase += -np.log(reflection_phase_shift).imag + + circuit.compose(affix_circuit, inplace=True) + + return circuit diff --git a/qiskit/quantum_info/synthesis/xx_decompose/decomposer.py b/qiskit/quantum_info/synthesis/xx_decompose/decomposer.py new file mode 100644 index 000000000000..0ebfc5d58b7f --- /dev/null +++ b/qiskit/quantum_info/synthesis/xx_decompose/decomposer.py @@ -0,0 +1,278 @@ +# 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. + +""" +Driver for a synthesis routine which emits optimal XX-based circuits. +""" + +import heapq +import math +from operator import itemgetter +from typing import Callable, Optional + +import numpy as np + +from qiskit.circuit.quantumcircuit import QuantumCircuit +from qiskit.circuit.library.standard_gates import RXXGate, RZXGate +from qiskit.exceptions import QiskitError +from qiskit.extensions import UnitaryGate +from qiskit.quantum_info.operators import Operator, average_gate_fidelity +from qiskit.quantum_info.synthesis.one_qubit_decompose import ONE_QUBIT_EULER_BASIS_GATES +from qiskit.quantum_info.synthesis.two_qubit_decompose import TwoQubitWeylDecomposition +from qiskit.transpiler.passes.optimization.optimize_1q_decomposition import ( + Optimize1qGatesDecomposition, +) + +from .circuits import apply_reflection, apply_shift, canonical_xx_circuit +from .utilities import EPSILON +from .polytopes import XXPolytope + + +def _average_infidelity(p, q): + """ + Computes the infidelity distance between two points p, q expressed in positive canonical + coordinates. + """ + + a0, b0, c0 = p + a1, b1, c1 = q + + return 1 - 1 / 20 * ( + 4 + + 16 + * ( + math.cos(a0 - a1) ** 2 * math.cos(b0 - b1) ** 2 * math.cos(c0 - c1) ** 2 + + math.sin(a0 - a1) ** 2 * math.sin(b0 - b1) ** 2 * math.sin(c0 - c1) ** 2 + ) + ) + + +class XXDecomposer: + """ + A class for optimal decomposition of 2-qubit unitaries into 2-qubit basis gates of XX type + (i.e., each locally equivalent to CAN(alpha, 0, 0) for a possibly varying alpha). + + Args: + euler_basis: Basis string provided to OneQubitEulerDecomposer for 1Q synthesis. + Defaults to "U". + embodiments: A dictionary mapping interaction strengths alpha to native circuits which + embody the gate CAN(alpha, 0, 0). Strengths are taken so that pi/2 represents the class + of a full CX. + backup_optimizer: If supplied, defers synthesis to this callable when XXDecomposer + has no efficient decomposition of its own. Useful for special cases involving 2 or 3 + applications of XX(pi/2), in which case standard synthesis methods provide lower + 1Q gate count. + + NOTE: If embodiments is not passed, or if an entry is missing, it will be populated as needed + using the method _default_embodiment. + """ + + def __init__( + self, + euler_basis: str = "U", + embodiments: Optional[dict] = None, + backup_optimizer: Optional[Callable] = None, + ): + self._decomposer1q = Optimize1qGatesDecomposition(ONE_QUBIT_EULER_BASIS_GATES[euler_basis]) + self.gate = RZXGate(np.pi / 2) + self.embodiments = embodiments if embodiments is not None else {} + self.backup_optimizer = backup_optimizer + + self._check_embodiments() + + @staticmethod + def _default_embodiment(strength): + """ + If the user does not provide a custom implementation of XX(strength), then this routine + defines a default implementation using RZX. + """ + xx_circuit = QuantumCircuit(2) + + # NOTE: One could branch here on `strength == np.pi / 2` and decide to use a CX-based + # circuit in this one case where it's available. + xx_circuit.h(0) + xx_circuit.rzx(strength, 0, 1) + xx_circuit.h(0) + + return xx_circuit + + def _check_embodiments(self): + """ + Checks that `self.embodiments` is populated with legal circuit embodiments: the key-value + pair (angle, circuit) satisfies Operator(circuit) approx RXX(angle).to_matrix(). + """ + + for angle, embodiment in self.embodiments.items(): + actual = Operator(RXXGate(angle)) + purported = Operator(embodiment) + if average_gate_fidelity(actual, purported) < 1 - EPSILON: + raise QiskitError( + f"RXX embodiment provided for angle {angle} disagrees with RXXGate({angle})" + ) + + @staticmethod + def _best_decomposition(canonical_coordinate, available_strengths): + """ + Finds the cheapest sequence of `available_strengths` which supports the best approximation + to `canonical_coordinate`. Returns a dictionary with keys "cost", "point", and "operations". + + NOTE: `canonical_coordinate` is a positive canonical coordinate. `strengths` is a dictionary + mapping the available strengths to their (infidelity) costs, with the strengths + themselves normalized so that pi/2 represents CX = RZX(pi/2). + """ + best_point, best_cost, best_sequence = [0, 0, 0], 1.0, [] + priority_queue = [] + heapq.heappush(priority_queue, (0, [])) + canonical_coordinate = np.array(canonical_coordinate) + + while True: + sequence_cost, sequence = heapq.heappop(priority_queue) + + strength_polytope = XXPolytope.from_strengths(*[x / 2 for x in sequence]) + candidate_point = strength_polytope.nearest(canonical_coordinate) + candidate_cost = sequence_cost + _average_infidelity( + canonical_coordinate, candidate_point + ) + + if candidate_cost < best_cost: + best_point, best_cost, best_sequence = candidate_point, candidate_cost, sequence + + if strength_polytope.member(canonical_coordinate): + break + + for strength, extra_cost in available_strengths.items(): + if len(sequence) == 0 or strength <= sequence[-1]: + heapq.heappush( + priority_queue, (sequence_cost + extra_cost, sequence + [strength]) + ) + + return {"point": best_point, "cost": best_cost, "sequence": best_sequence} + + def num_basis_gates(self, unitary): + """ + Counts the number of gates that would be emitted during re-synthesis. + + NOTE: Used by ConsolidateBlocks. + """ + strengths = self._strength_to_infidelity(1.0) + + # get the associated _positive_ canonical coordinate + weyl_decomposition = TwoQubitWeylDecomposition(unitary) + target = [getattr(weyl_decomposition, x) for x in ("a", "b", "c")] + if target[-1] < -EPSILON: + target = [np.pi / 2 - target[0], target[1], -target[2]] + + best_sequence = self._best_decomposition(target, strengths)["sequence"] + return len(best_sequence) + + @staticmethod + def _strength_to_infidelity(basis_fidelity, approximate=False): + """ + Converts a dictionary mapping ZX strengths to fidelities to a dictionary mapping ZX + strengths to infidelities. Also supports one of the other formats Qiskit uses: if only a + lone float is supplied, it extends it from CX over CX/2 and CX/3 by linear decay. + """ + + if isinstance(basis_fidelity, float): + if not approximate: + slope, offset = 1e-10, 1e-12 + else: + slope, offset = (1 - basis_fidelity) / 2, (1 - basis_fidelity) / 2 + return { + strength: slope * strength / (np.pi / 2) + offset + for strength in [np.pi / 2, np.pi / 4, np.pi / 6] + } + elif isinstance(basis_fidelity, dict): + return { + strength: (1 - fidelity if approximate else 1e-12 + 1e-10 * strength / (np.pi / 2)) + for (strength, fidelity) in basis_fidelity.items() + } + + raise TypeError("Unknown basis_fidelity payload.") + + def __call__(self, unitary, basis_fidelity=1.0, approximate=True): + """ + Fashions a circuit which (perhaps `approximate`ly) models the special unitary operation + `unitary`, using the circuit templates supplied at initialization as `embodiments`. The + routine uses `basis_fidelity` to select the optimal circuit template, including when + performing exact synthesis; the contents of `basis_fidelity` is a dictionary mapping + interaction strengths (scaled so that CX = RZX(pi/2) corresponds to pi/2) to circuit + fidelities. + + Args: + unitary (Operator or ndarray): 4x4 unitary to synthesize. + basis_fidelity (dict or float): Fidelity of basis gates. Can be either (1) a dictionary + mapping XX angle values to fidelity at that angle; or (2) a single float f, + interpreted as {pi: f, pi/2: f/2, pi/3: f/3} . + approximate (bool): Approximates if basis fidelities are less than 1.0 . + Returns: + QuantumCircuit: Synthesized circuit. + """ + strength_to_infidelity = self._strength_to_infidelity( + basis_fidelity, approximate=approximate + ) + + # get the associated _positive_ canonical coordinate + weyl_decomposition = TwoQubitWeylDecomposition(unitary) + target = [getattr(weyl_decomposition, x) for x in ("a", "b", "c")] + if target[-1] < -EPSILON: + target = [np.pi / 2 - target[0], target[1], -target[2]] + + # scan for the best point + best_point, best_sequence = itemgetter("point", "sequence")( + self._best_decomposition(target, strength_to_infidelity) + ) + # build the circuit building this canonical gate + embodiments = { + k: self.embodiments.get(k, self._default_embodiment(k)) + for k, v in strength_to_infidelity.items() + } + circuit = canonical_xx_circuit(best_point, best_sequence, embodiments) + + if best_sequence == [np.pi / 2, np.pi / 2, np.pi / 2] and self.backup_optimizer is not None: + return self.backup_optimizer(unitary, basis_fidelity=basis_fidelity) + + # change to positive canonical coordinates + if weyl_decomposition.c >= -EPSILON: + # if they're the same... + corrected_circuit = QuantumCircuit(2) + corrected_circuit.rz(np.pi, [0]) + corrected_circuit.compose(circuit, [0, 1], inplace=True) + corrected_circuit.rz(-np.pi, [0]) + circuit = corrected_circuit + else: + # else they're in the "positive" scissors part... + corrected_circuit = QuantumCircuit(2) + _, source_reflection, _ = apply_reflection("reflect XX, ZZ", [0, 0, 0]) + _, source_shift, _ = apply_shift("X shift", [0, 0, 0]) + + corrected_circuit.compose(source_reflection.inverse(), inplace=True) + corrected_circuit.rz(np.pi, [0]) + corrected_circuit.compose(circuit, [0, 1], inplace=True) + corrected_circuit.rz(-np.pi, [0]) + corrected_circuit.compose(source_shift.inverse(), inplace=True) + corrected_circuit.compose(source_reflection, inplace=True) + corrected_circuit.global_phase += np.pi / 2 + + circuit = corrected_circuit + + circ = QuantumCircuit(2, global_phase=weyl_decomposition.global_phase) + + circ.append(UnitaryGate(weyl_decomposition.K2r), [0]) + circ.append(UnitaryGate(weyl_decomposition.K2l), [1]) + circ.compose(circuit, [0, 1], inplace=True) + circ.append(UnitaryGate(weyl_decomposition.K1r), [0]) + circ.append(UnitaryGate(weyl_decomposition.K1l), [1]) + + circ = self._decomposer1q(circ) + + return circ diff --git a/qiskit/quantum_info/synthesis/xx_decompose/paths.py b/qiskit/quantum_info/synthesis/xx_decompose/paths.py new file mode 100644 index 000000000000..739202f8801a --- /dev/null +++ b/qiskit/quantum_info/synthesis/xx_decompose/paths.py @@ -0,0 +1,395 @@ +# 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. + +""" +Routines for producing right-angled paths through the Weyl alcove. Consider a set of native +interactions with an associated minimal covering set of minimum-cost circuit polytopes, as well as a +target coordinate. The coverage set associates to the target coordinate a circuit type +C = (O1 ... On) consisting of a sequence of native interactions Oj. A _path_ is a sequence +(I P1 ... Pn) of intermediate Weyl points, where Pj is accessible from P(j-1) by Oj. A path is said +to be _right-angled_ when at each step one coordinate is fixed (up to possible Weyl reflection) when +expressed in canonical coordinates. + +The key inputs to our method are: + ++ A family of "b coordinates" which describe the target canonical coordinate. ++ A family of "a coordinates" which describe the source canonical coordinate. ++ A sequence of interaction strengths for which the b-coordinate can be modeled, with one selected + to be stripped from the sequence ("beta"). The others are bundled as the sum of the + sequence (s+), its maximum value (s1), and its second maximum value (s2). + +Given the b-coordinate and a set of intersection strengths, the procedure for backsolving for the +a-coordinates is then extracted from the monodromy polytope. + +NOTE: The constants in this file are auto-generated and are not meant to be edited by hand / read. +""" + +import numpy as np + +from .polytopes import ConvexPolytopeData, PolytopeData, manual_get_vertex, polytope_has_element + + +def get_augmented_coordinate(target_coordinate, strengths): + """ + Assembles a coordinate in the system used by `xx_region_polytope`. + """ + *strengths, beta = strengths + strengths = sorted(strengths + [0, 0]) + interaction_coordinate = [sum(strengths), strengths[-1], strengths[-2], beta] + return [*target_coordinate, *interaction_coordinate] + + +def decomposition_hop(target_coordinate, strengths): + """ + Given a `target_coordinate` and a list of interaction `strengths`, produces a new canonical + coordinate which is one step back along `strengths`. + + `target_coordinate` is taken to be in positive canonical coordinates, and the entries of + strengths are taken to be [0, pi], so that (sj / 2, 0, 0) is a positive canonical coordinate. + """ + + target_coordinate = [x / (np.pi / 2) for x in target_coordinate] + strengths = [x / np.pi for x in strengths] + + augmented_coordinate = get_augmented_coordinate(target_coordinate, strengths) + specialized_polytope = None + for cp in xx_region_polytope.convex_subpolytopes: + if not polytope_has_element(cp, augmented_coordinate): + continue + if "AF=B1" in cp.name: + af, _, _ = target_coordinate + elif "AF=B2" in cp.name: + _, af, _ = target_coordinate + elif "AF=B3" in cp.name: + _, _, af = target_coordinate + else: + raise ValueError("Couldn't find a coordinate to fix.") + + raw_convex_polytope = next( + (cpp for cpp in xx_lift_polytope.convex_subpolytopes if cpp.name == cp.name), None + ) + + coefficient_dict = {} + for inequality in raw_convex_polytope.inequalities: + if inequality[1] == 0 and inequality[2] == 0: + continue + offset = ( + inequality[0] # old constant term + + inequality[3] * af + + inequality[4] * augmented_coordinate[0] # b1 + + inequality[5] * augmented_coordinate[1] # b2 + + inequality[6] * augmented_coordinate[2] # b3 + + inequality[7] * augmented_coordinate[3] # s+ + + inequality[8] * augmented_coordinate[4] # s1 + + inequality[9] * augmented_coordinate[5] # s2 + + inequality[10] * augmented_coordinate[6] # beta + ) + + if offset <= coefficient_dict.get((inequality[1], inequality[2]), offset): + coefficient_dict[(inequality[1], inequality[2])] = offset + + specialized_polytope = PolytopeData( + convex_subpolytopes=[ + ConvexPolytopeData( + inequalities=[[v, h, l] for ((h, l), v) in coefficient_dict.items()] + ) + ] + ) + + break + + if specialized_polytope is None: + raise ValueError("Failed to match a constrained_polytope summand.") + + ah, al = manual_get_vertex(specialized_polytope) + return [x * (np.pi / 2) for x in sorted([ah, al, af], reverse=True)] + + +xx_region_polytope = PolytopeData( + convex_subpolytopes=[ + ConvexPolytopeData( + inequalities=[ + [0, 0, 0, 0, 0, 1, -1, 0], + [0, 0, 0, 0, 0, 0, 1, 0], + [1, 0, 0, 0, 0, -2, 0, 0], + [1, 0, 0, 0, 0, 0, 0, -2], + [0, 0, 1, -1, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, 0, 1, -1, -1, 0], + [1, -1, -1, 0, 0, 0, 0, 0], + [0, -1, -1, -1, 1, 0, 0, 1], + [0, 1, -1, 0, 0, 0, 0, 0], + [0, 1, -1, -1, 1, -2, 0, 1], + [0, 1, -1, -1, 1, 0, 0, -1], + [0, 0, 0, -1, 1, -1, 0, 0], + [0, 0, -1, 0, 1, -1, -1, 1], + [0, 0, 0, 0, 0, 0, 0, 1], + ], + equalities=[], + name="I ∩ A alcove ∩ " + "A unreflected ∩ ah slant ∩ al frustrum ∩ B alcove ∩ B unreflected ∩ AF=B3", + ), + ConvexPolytopeData( + inequalities=[ + [0, 0, 0, 0, 0, 1, -1, 0], + [0, 0, 0, 0, 0, 0, 1, 0], + [1, 0, 0, 0, 0, -2, 0, 0], + [1, 0, 0, 0, 0, 0, 0, -2], + [0, 0, 1, -1, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, 0, 1, -1, -1, 0], + [1, -1, -1, 0, 0, 0, 0, 0], + [1, -1, -1, -1, 1, -2, 0, 1], + [0, 1, -1, 0, 0, 0, 0, 0], + [-1, 1, -1, -1, 1, 0, 0, 1], + [1, -1, -1, -1, 1, 0, 0, -1], + [0, 0, 0, -1, 1, -1, 0, 0], + [0, 0, -1, 0, 1, -1, -1, 1], + [0, 0, 0, 0, 0, 0, 0, 1], + ], + equalities=[], + name="I ∩ A alcove ∩ " + "A reflected ∩ ah strength ∩ al frustrum ∩ B alcove ∩ B reflected ∩ AF=B3", + ), + ConvexPolytopeData( + inequalities=[ + [0, 0, 0, 0, 0, 1, -1, 0], + [0, 0, 0, 0, 0, 0, 1, 0], + [1, 0, 0, 0, 0, -2, 0, 0], + [1, 0, 0, 0, 0, 0, 0, -2], + [0, 1, -1, 0, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0], + [1, -1, -1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 1, -1, -1, 0], + [0, 1, -1, -1, 1, -2, 0, 1], + [0, -1, -1, -1, 1, 0, 0, 1], + [0, 0, 1, -1, 0, 0, 0, 0], + [1, -1, 1, -1, 0, 0, 0, -1], + [0, 1, 1, -1, 1, -2, 0, -1], + [0, -1, 1, -1, 1, 0, 0, -1], + [0, 0, 0, -1, 1, -1, -1, 1], + [0, 0, 0, 0, 0, 0, 0, 1], + ], + equalities=[], + name="I ∩ A alcove ∩ " + "A unreflected ∩ af slant ∩ al frustrum ∩ B alcove ∩ B unreflected ∩ AF=B1", + ), + ConvexPolytopeData( + inequalities=[ + [0, 0, 0, 0, 0, 1, -1, 0], + [0, 0, 0, 0, 0, 0, 1, 0], + [1, 0, 0, 0, 0, -2, 0, 0], + [1, 0, 0, 0, 0, 0, 0, -2], + [0, 1, -1, 0, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0], + [1, -1, -1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 1, -1, -1, 0], + [-1, 1, -1, -1, 1, 0, 0, 1], + [1, -1, -1, -1, 1, -2, 0, 1], + [0, 0, 1, -1, 0, 0, 0, 0], + [1, -1, 1, -1, 0, 0, 0, -1], + [-1, 1, 1, -1, 1, 0, 0, -1], + [1, -1, 1, -1, 1, -2, 0, -1], + [0, 0, 0, -1, 1, -1, -1, 1], + [0, 0, 0, 0, 0, 0, 0, 1], + ], + equalities=[], + name="I ∩ A alcove ∩ " + "A reflected ∩ af strength ∩ al frustrum ∩ B alcove ∩ B reflected ∩ AF=B1", + ), + ] +) +""" +For any choice of sequence of strengths, the theory of the monodromy polytope yields a polytope P +so that (b1, b2, b3) belongs to P if and only if (b1, b2, b3) is the positive canonical coordinate +of a program expressible by a circuit whose 2Q interactions are the prescribed sequence of +fractional CX strengths, interleaved with arbitrary single-qubit gates. The polytope above has the +same property as P, but (1) it is parametrized over the sequence of strengths, and (2) it is broken +into a certain sum of convex regions. Each region is tagged with a name, and it is the projection of +the region of `xx_lift_polytope` (see below) with the corresponding name. + +The coordinates are [k, b1, b2, b3, s+, s1, s2, beta]. +""" + + +xx_lift_polytope = PolytopeData( + convex_subpolytopes=[ + ConvexPolytopeData( + inequalities=[ + [0, 0, 0, 0, 0, 0, 0, 1, -1, -1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0], + [1, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0], + [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1], + [0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0], + [1, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0], + [1, -1, 0, -1, 0, 0, 0, 0, 0, 0, 0], + [0, -1, -1, -1, 0, 0, 0, 1, 0, 0, 0], + [0, 1, -1, -1, 0, 0, 0, 1, -2, 0, 0], + [0, 0, -1, 0, 0, 0, 0, 1, -1, -1, 0], + [0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0], + [1, 0, 0, 0, -1, -1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, -1, -1, -1, 1, 0, 0, 1], + [0, 0, 0, 0, 1, -1, -1, 1, -2, 0, 1], + [0, 0, 0, 0, 1, -1, -1, 1, 0, 0, -1], + [0, 0, 0, 0, 0, 0, -1, 1, -1, -1, 1], + [0, 0, 0, 0, 0, 0, -1, 1, -1, 0, 0], + [0, -1, -1, 0, 1, 1, 0, 0, 0, 0, 1], + [2, -1, -1, 0, -1, -1, 0, 0, 0, 0, -1], + [0, 1, 1, 0, -1, -1, 0, 0, 0, 0, 1], + [0, -1, 1, 0, 1, -1, 0, 0, 0, 0, 1], + [0, 1, -1, 0, -1, 1, 0, 0, 0, 0, 1], + [0, 1, -1, 0, 1, -1, 0, 0, 0, 0, -1], + ], + equalities=[[0, 0, 0, 1, 0, 0, -1, 0, 0, 0, 0]], + name="I ∩ A alcove ∩ " + "A unreflected ∩ ah slant ∩ al frustrum ∩ B alcove ∩ B unreflected ∩ AF=B3", + ), + ConvexPolytopeData( + inequalities=[ + [0, 0, 0, 0, 0, 0, 0, 1, -1, -1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0], + [1, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0], + [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1], + [0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0], + [1, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0], + [1, -1, 0, -1, 0, 0, 0, 0, 0, 0, 0], + [0, -1, -1, -1, 0, 0, 0, 1, 0, 0, 0], + [0, -1, -1, 1, 0, 0, 0, 1, -2, 0, 0], + [0, 0, -1, 0, 0, 0, 0, 1, -1, -1, 0], + [0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0], + [1, 0, 0, 0, -1, -1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, -1, -1, -1, 1, 0, 0, 1], + [0, 0, 0, 0, 1, -1, -1, 1, -2, 0, 1], + [0, 0, 0, 0, 1, -1, -1, 1, 0, 0, -1], + [0, 0, 0, 0, 0, 0, -1, 1, -1, -1, 1], + [0, 0, 0, 0, 0, 0, -1, 1, -1, 0, 0], + [0, -1, -1, 0, 0, 1, 1, 0, 0, 0, 1], + [2, -1, -1, 0, 0, -1, -1, 0, 0, 0, -1], + [0, 1, 1, 0, 0, -1, -1, 0, 0, 0, 1], + [0, -1, 1, 0, 0, 1, -1, 0, 0, 0, 1], + [0, 1, -1, 0, 0, -1, 1, 0, 0, 0, 1], + [0, 1, -1, 0, 0, 1, -1, 0, 0, 0, -1], + ], + equalities=[[0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0]], + name="I ∩ A alcove ∩ " + "A unreflected ∩ af slant ∩ al frustrum ∩ B alcove ∩ B unreflected ∩ AF=B1", + ), + ConvexPolytopeData( + inequalities=[ + [0, 0, 0, 0, 0, 0, 0, 1, -1, -1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0], + [1, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0], + [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1], + [0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0], + [1, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0], + [1, -1, 0, -1, 0, 0, 0, 0, 0, 0, 0], + [1, -1, -1, -1, 0, 0, 0, 1, -2, 0, 0], + [-1, -1, -1, 1, 0, 0, 0, 1, 0, 0, 0], + [0, 0, -1, 0, 0, 0, 0, 1, -1, -1, 0], + [0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0], + [1, 0, 0, 0, -1, -1, 0, 0, 0, 0, 0], + [-1, 0, 0, 0, 1, -1, -1, 1, 0, 0, 1], + [1, 0, 0, 0, -1, -1, -1, 1, 0, 0, -1], + [1, 0, 0, 0, -1, -1, -1, 1, -2, 0, 1], + [0, 0, 0, 0, 0, 0, -1, 1, -1, -1, 1], + [0, 0, 0, 0, 0, 0, -1, 1, -1, 0, 0], + [0, -1, -1, 0, 0, 1, 1, 0, 0, 0, 1], + [2, -1, -1, 0, 0, -1, -1, 0, 0, 0, -1], + [0, 1, 1, 0, 0, -1, -1, 0, 0, 0, 1], + [0, -1, 1, 0, 0, 1, -1, 0, 0, 0, 1], + [0, 1, -1, 0, 0, -1, 1, 0, 0, 0, 1], + [0, 1, -1, 0, 0, 1, -1, 0, 0, 0, -1], + ], + equalities=[[0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0]], + name="I ∩ A alcove ∩ " + "A reflected ∩ af strength ∩ al frustrum ∩ B alcove ∩ B reflected ∩ AF=B1", + ), + ConvexPolytopeData( + inequalities=[ + [0, 0, 0, 0, 0, 0, 0, 1, -1, -1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0], + [1, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0], + [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1], + [0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0], + [1, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0], + [1, -1, 0, -1, 0, 0, 0, 0, 0, 0, 0], + [1, -1, -1, -1, 0, 0, 0, 1, -2, 0, 0], + [-1, 1, -1, -1, 0, 0, 0, 1, 0, 0, 0], + [0, 0, -1, 0, 0, 0, 0, 1, -1, -1, 0], + [0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0], + [1, 0, 0, 0, -1, -1, 0, 0, 0, 0, 0], + [-1, 0, 0, 0, 1, -1, -1, 1, 0, 0, 1], + [1, 0, 0, 0, -1, -1, -1, 1, 0, 0, -1], + [1, 0, 0, 0, -1, -1, -1, 1, -2, 0, 1], + [0, 0, 0, 0, 0, 0, -1, 1, -1, -1, 1], + [0, 0, 0, 0, 0, 0, -1, 1, -1, 0, 0], + [0, -1, -1, 0, 1, 1, 0, 0, 0, 0, 1], + [2, -1, -1, 0, -1, -1, 0, 0, 0, 0, -1], + [0, 1, 1, 0, -1, -1, 0, 0, 0, 0, 1], + [0, -1, 1, 0, 1, -1, 0, 0, 0, 0, 1], + [0, 1, -1, 0, -1, 1, 0, 0, 0, 0, 1], + [0, 1, -1, 0, 1, -1, 0, 0, 0, 0, -1], + ], + equalities=[[0, 0, 0, 1, 0, 0, -1, 0, 0, 0, 0]], + name="I ∩ A alcove ∩ " + "A reflected ∩ ah strength ∩ al frustrum ∩ B alcove ∩ B reflected ∩ AF=B3", + ), + ] +) +""" +For any choice of sequence of strengths, the theory of the monodromy polytope yields a polytope Q +so that (a1, a2, a3, b1, b2, b3) belongs to Q if and only if: + + * (b1, b2, b3) is the positive canonical coordinate of a program expressible by a circuit whose + 2Q interactions are the prescribed sequence of fractional CX strengths, interleaved with + arbitrary single-qubit gates. + * (a1, a2, a3) is the positive canonical coordinate of a program expressible by a circuit whose + 2Q interactions are the prescribed sequence of fractional CX strengths _with the last + interaction omitted_, interleaved with arbitrary single-qubit gates. + * (b1, b2, b3) is accessible from (a1, a2, a3) using single-qubit gates and a single application + of the final fractional CX strength. + +The polytope above is a slight variation on Q: + + (1) It is parametrized over the sequence of strengths. + (2) It is broken into a certain sum of convex regions. If (b1, b2, b3) belongs to the projection + of one of these convex regions, then that same region's projection to (a1, a2, a3) is + guaranteed to be nonempty. (To test this condition, see `xx_region_polytope`.) + +Points in this region can thus be used to calculate circuits using `decompose_xxyy_into_xxyy_xx`. + +The coordinates are [k, ah, al, af, b1, b2, b3, s+, s1, s2, beta]. +""" diff --git a/qiskit/quantum_info/synthesis/xx_decompose/polytopes.py b/qiskit/quantum_info/synthesis/xx_decompose/polytopes.py new file mode 100644 index 000000000000..d069b806712e --- /dev/null +++ b/qiskit/quantum_info/synthesis/xx_decompose/polytopes.py @@ -0,0 +1,265 @@ +# 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. + +""" +Defines bare dataclasses which house polytope information, as well as a specialized data structure +which describes those two-qubit programs accessible to a given sequence of XX-type interactions. +""" + +from copy import copy +from dataclasses import dataclass, field +from itertools import combinations +from typing import List + +import numpy as np + +from qiskit.exceptions import QiskitError + +from .utilities import EPSILON + + +@dataclass +class ConvexPolytopeData: + """ + The raw data underlying a ConvexPolytope. Describes a single convex + polytope, specified by families of `inequalities` and `equalities`, each + entry of which respectively corresponds to + + inequalities[j][0] + sum_i inequalities[j][i] * xi >= 0 + + and + + equalities[j][0] + sum_i equalities[j][i] * xi == 0. + """ + + inequalities: List[List[int]] + equalities: List[List[int]] = field(default_factory=list) + name: str = "" + + +@dataclass +class PolytopeData: + """ + The raw data of a union of convex polytopes. + """ + + convex_subpolytopes: List[ConvexPolytopeData] + + +def polytope_has_element(polytope, point): + """ + Tests whether `polytope` contains `point. + """ + return all( + -EPSILON <= inequality[0] + sum(x * y for x, y in zip(point, inequality[1:])) + for inequality in polytope.inequalities + ) and all( + abs(equality[0] + sum(x * y for x, y in zip(point, equality[1:]))) <= EPSILON + for equality in polytope.equalities + ) + + +def manual_get_vertex(polytope, seed=42): + """ + Returns a single random vertex from `polytope`. + """ + rng = np.random.default_rng(seed) + + if isinstance(polytope, PolytopeData): + paragraphs = copy(polytope.convex_subpolytopes) + elif isinstance(polytope, ConvexPolytopeData): + paragraphs = [polytope] + else: + raise TypeError(f"{type(polytope)} is not polytope-like.") + + rng.shuffle(paragraphs) + for convex_subpolytope in paragraphs: + sentences = convex_subpolytope.inequalities + convex_subpolytope.equalities + if len(sentences) == 0: + continue + dimension = len(sentences[0]) - 1 + rng.shuffle(sentences) + for inequalities in combinations(sentences, dimension): + matrix = np.array([x[1:] for x in inequalities]) + b = np.array([x[0] for x in inequalities]) + try: + vertex = np.linalg.inv(-matrix) @ b + if polytope_has_element(convex_subpolytope, vertex): + return vertex + except np.linalg.LinAlgError: + pass + + raise QiskitError(f"Polytope has no feasible solutions:\n{polytope}") + + +@dataclass +class XXPolytope: + """ + Describes those two-qubit programs accessible to a given sequence of XX-type interactions. + + NOTE: Strengths are normalized so that CX corresponds to pi / 4, which differs from Qiskit's + conventions around RZX elsewhere. + """ + + # NOTE: This is _not_ a subclass of PolytopeData, because we're never going to call slow, + # generic PolytopeData functions on it. + + total_strength: float = 0.0 + max_strength: float = 0.0 + place_strength: float = 0.0 + + @classmethod + def from_strengths(cls, *strengths): + """ + Constructs an XXPolytope from a sequence of strengths. + """ + total_strength, max_strength, place_strength = 0, 0, 0 + for strength in strengths: + total_strength += strength + if strength >= max_strength: + max_strength, place_strength = strength, max_strength + elif strength >= place_strength: + place_strength = strength + + return XXPolytope( + total_strength=total_strength, max_strength=max_strength, place_strength=place_strength + ) + + def add_strength(self, new_strength: float = 0.0): + """ + Returns a new XXPolytope with one new XX interaction appended. + """ + return XXPolytope( + total_strength=self.total_strength + new_strength, + max_strength=max(self.max_strength, new_strength), + place_strength=( + self.max_strength + if new_strength > self.max_strength + else new_strength + if new_strength > self.place_strength + else self.place_strength + ), + ) + + @property + def _offsets(self): + """ + Returns b with A*x + b ≥ 0 iff x belongs to the XXPolytope. + """ + return np.array( + [ + 0, + 0, + 0, + np.pi / 2, + self.total_strength, + self.total_strength - 2 * self.max_strength, + self.total_strength - self.max_strength - self.place_strength, + ] + ) + + def member(self, point): + """ + Returns True when `point` is a member of `self`. + """ + global A # pylint:disable=global-statement + + reflected_point = point.copy().reshape(-1, 3) + rows = reflected_point[:, 0] >= np.pi / 4 + EPSILON + reflected_point[rows, 0] = np.pi / 2 - reflected_point[rows, 0] + reflected_point = reflected_point.reshape(point.shape) + + return np.all( + self._offsets + np.einsum("ij,...j->...i", A, reflected_point) >= -EPSILON, axis=-1 + ) + + def nearest(self, point): + """ + Finds the nearest point (in Euclidean or infidelity distance) to `self`. + """ + # pylint:disable=invalid-name + + # NOTE: A CAS says that there are no degenerate double intersections, and the only + # degenerate triple intersections are + # + # (1, -1, 0), (0, 0, 1), (-1, 1, 1) and (1, 1, 0), (0, 0, 1), (1, 1, 1). + # + # Skipping this pair won't save much work, so we don't bother. + + global A1, A1inv, A2, A2inv, A3, A3inv # pylint:disable=global-statement + # These global variables contain projection matrices, computed once-and-for-all, which + # produce the Euclidean-nearest projection. + + if isinstance(point, np.ndarray) and len(point.shape) == 1: + y0 = point.copy() + elif isinstance(point, list): + y0 = np.array(point) + else: + raise TypeError(f"Can't handle type of point: {point} ({type(point)})") + + reflected_p = y0[0] > np.pi / 4 + EPSILON + if reflected_p: + y0[0] = np.pi / 2 - y0[0] + + # short circuit in codimension 0 + if self.member(y0): + if reflected_p: + y0[0] = np.pi / 2 - y0[0] + return y0 + + # codimension 1 + b1 = self._offsets.reshape(7, 1) + A1y0 = np.einsum("ijk,k->ij", A1, y0) + nearest1 = np.einsum("ijk,ik->ij", A1inv, b1 + A1y0) - y0 + + # codimension 2 + b2 = np.array([*combinations(self._offsets, 2)]) + A2y0 = np.einsum("ijk,k->ij", A2, y0) + nearest2 = np.einsum("ijk,ik->ij", A2inv, b2 + A2y0) - y0 + + # codimension 3 + b3 = np.array([*combinations(self._offsets, 3)]) + nearest3 = np.einsum("ijk,ik->ij", A3inv, b3) + + # pick the nearest + nearest = -np.concatenate([nearest1, nearest2, nearest3]) + nearest = nearest[self.member(nearest)] + smallest_index = np.argmin(np.linalg.norm(nearest - y0, axis=1)) + + if reflected_p: + nearest[smallest_index][0] = np.pi / 2 - nearest[smallest_index][0] + return nearest[smallest_index] + + +A = np.array( + [ + [1, -1, 0], # a ≥ b + [0, 1, -1], # b ≥ c + [0, 0, 1], # c ≥ 0 + [-1, -1, 0], # pi/2 ≥ a + b + [-1, -1, -1], # strength + [1, -1, -1], # slant + [0, 0, -1], # frustrum + ] +) +A1 = A.reshape(-1, 1, 3) # pylint:disable=too-many-function-args +A1inv = np.linalg.pinv(A1) +A2 = np.array([np.array([x, y], dtype=float) for (x, y) in combinations(A, 2)]) +A2inv = np.linalg.pinv(A2) +A3 = np.array([np.array([x, y, z], dtype=float) for (x, y, z) in combinations(A, 3)]) +A3inv = np.linalg.pinv(A3) +""" +These globals house matrices, computed once-and-for-all, which project onto the Euclidean-nearest +planes spanned by the planes/points/lines defined by the faces of any XX polytope. + +See `XXPolytope.nearest`. +""" diff --git a/qiskit/quantum_info/synthesis/xx_decompose/utilities.py b/qiskit/quantum_info/synthesis/xx_decompose/utilities.py new file mode 100644 index 000000000000..62cd9ebbb7d0 --- /dev/null +++ b/qiskit/quantum_info/synthesis/xx_decompose/utilities.py @@ -0,0 +1,39 @@ +# 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. + +""" +Depository for generic utility snippets. +""" + +import warnings + +import numpy as np + + +EPSILON = 1e-6 # Fraction(1, 1_000_000) + + +# TODO: THIS IS A STOPGAP!!! +def safe_arccos(numerator, denominator): + """ + Computes arccos(n/d) with different (better?) numerical stability. + """ + threshold = 0.005 + + if abs(numerator) > abs(denominator) and abs(numerator - denominator) < threshold: + return 0.0 + elif abs(numerator) > abs(denominator) and abs(numerator + denominator) < threshold: + return np.pi + else: + with warnings.catch_warnings(): + warnings.filterwarnings("ignore", category=RuntimeWarning) + return np.arccos(numerator / denominator) diff --git a/qiskit/quantum_info/synthesis/xx_decompose/weyl.py b/qiskit/quantum_info/synthesis/xx_decompose/weyl.py new file mode 100644 index 000000000000..c0c30c443576 --- /dev/null +++ b/qiskit/quantum_info/synthesis/xx_decompose/weyl.py @@ -0,0 +1,132 @@ +# 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. + +""" +Simple circuit constructors for Weyl reflections. +""" + +import numpy as np + +from qiskit.circuit.quantumcircuit import QuantumCircuit +from qiskit.circuit.library.standard_gates import RXGate, RYGate, RZGate + + +reflection_options = { + "no reflection": ([1, 1, 1], 1, []), + "reflect XX, YY": ([-1, -1, 1], 1, [RZGate]), + "reflect XX, ZZ": ([-1, 1, -1], 1, [RYGate]), + "reflect YY, ZZ": ([1, -1, -1], 1, [RXGate]), +} +""" +A table of available reflection transformations on canonical coordinates. +Entries take the form + + readable_name: (reflection scalars, global phase, [gate constructors]), + +where reflection scalars (a, b, c) model the map (x, y, z) |-> (ax, by, cz), +global phase is a complex unit, and gate constructors are applied in sequence +and by conjugation to the first qubit and are passed pi as a parameter. +""" + +shift_options = { + "no shift": ([0, 0, 0], 1, []), + "Z shift": ([0, 0, 1], 1j, [RZGate]), + "Y shift": ([0, 1, 0], -1j, [RYGate]), + "Y,Z shift": ([0, 1, 1], 1, [RYGate, RZGate]), + "X shift": ([1, 0, 0], -1j, [RXGate]), + "X,Z shift": ([1, 0, 1], 1, [RXGate, RZGate]), + "X,Y shift": ([1, 1, 0], -1, [RXGate, RYGate]), + "X,Y,Z shift": ([1, 1, 1], -1j, [RXGate, RYGate, RZGate]), +} +""" +A table of available shift transformations on canonical coordinates. Entries +take the form + + readable name: (shift scalars, global phase, [gate constructors]), + +where shift scalars model the map + + (x, y, z) |-> (x + a pi / 2, y + b pi / 2, z + c pi / 2) , + +global phase is a complex unit, and gate constructors are applied to the first +and second qubits and are passed pi as a parameter. +""" + + +def apply_reflection(reflection_name, coordinate): + """ + Given a reflection type and a canonical coordinate, applies the reflection + and describes a circuit which enacts the reflection + a global phase shift. + """ + reflection_scalars, reflection_phase_shift, source_reflection_gates = reflection_options[ + reflection_name + ] + reflected_coord = [x * y for x, y in zip(reflection_scalars, coordinate)] + source_reflection = QuantumCircuit(2) + for gate in source_reflection_gates: + source_reflection.append(gate(np.pi), [0]) + + return reflected_coord, source_reflection, reflection_phase_shift + + +def apply_shift(shift_name, coordinate): + """ + Given a shift type and a canonical coordinate, applies the shift and + describes a circuit which enacts the shift + a global phase shift. + """ + shift_scalars, shift_phase_shift, source_shift_gates = shift_options[shift_name] + shifted_coord = [np.pi / 2 * x + y for x, y in zip(shift_scalars, coordinate)] + + source_shift = QuantumCircuit(2) + for gate in source_shift_gates: + source_shift.append(gate(np.pi), [0]) + source_shift.append(gate(np.pi), [1]) + + return shifted_coord, source_shift, shift_phase_shift + + +def canonical_rotation_circuit(first_index, second_index): + """ + Given a pair of distinct indices 0 ≤ (first_index, second_index) ≤ 2, + produces a two-qubit circuit which rotates a canonical gate + + a0 XX + a1 YY + a2 ZZ + + into + + a[first] XX + a[second] YY + a[other] ZZ . + """ + conj = QuantumCircuit(2) + + if (0, 1) == (first_index, second_index): + pass # no need to do anything + elif (0, 2) == (first_index, second_index): + conj.rx(-np.pi / 2, [0]) + conj.rx(np.pi / 2, [1]) + elif (1, 0) == (first_index, second_index): + conj.rz(-np.pi / 2, [0]) + conj.rz(-np.pi / 2, [1]) + elif (1, 2) == (first_index, second_index): + conj.rz(np.pi / 2, [0]) + conj.rz(np.pi / 2, [1]) + conj.ry(np.pi / 2, [0]) + conj.ry(-np.pi / 2, [1]) + elif (2, 0) == (first_index, second_index): + conj.rz(np.pi / 2, [0]) + conj.rz(np.pi / 2, [1]) + conj.rx(np.pi / 2, [0]) + conj.rx(-np.pi / 2, [1]) + elif (2, 1) == (first_index, second_index): + conj.ry(np.pi / 2, [0]) + conj.ry(-np.pi / 2, [1]) + + return conj diff --git a/qiskit/transpiler/layout.py b/qiskit/transpiler/layout.py index b823dda90096..0337e0afd360 100644 --- a/qiskit/transpiler/layout.py +++ b/qiskit/transpiler/layout.py @@ -308,7 +308,10 @@ def from_intlist(int_list, *qregs): num_qubits = sum(reg.size for reg in qregs) # Check if list is too short to cover all qubits if len(int_list) != num_qubits: - raise LayoutError("Integer list length must equal number of qubits in circuit.") + raise LayoutError( + f"Integer list length ({len(int_list)}) must equal number of qubits " + f"in circuit ({num_qubits}): {int_list}." + ) out = Layout() main_idx = 0 for qreg in qregs: diff --git a/qiskit/transpiler/passes/optimization/consolidate_blocks.py b/qiskit/transpiler/passes/optimization/consolidate_blocks.py index 1ecf6787ec71..290e3821ddaf 100644 --- a/qiskit/transpiler/passes/optimization/consolidate_blocks.py +++ b/qiskit/transpiler/passes/optimization/consolidate_blocks.py @@ -56,11 +56,7 @@ def __init__(self, kak_basis_gate=None, force_consolidate=False, basis_gates=Non if kak_basis_gate is not None: self.decomposer = TwoQubitBasisDecomposer(kak_basis_gate) elif basis_gates is not None: - kak_basis_gate = unitary_synthesis._choose_kak_gate(basis_gates) - if kak_basis_gate is not None: - self.decomposer = TwoQubitBasisDecomposer(kak_basis_gate) - else: - self.decomposer = None + self.decomposer = unitary_synthesis._basis_gates_to_decomposer_2q(basis_gates) else: self.decomposer = TwoQubitBasisDecomposer(CXGate()) diff --git a/qiskit/transpiler/passes/synthesis/unitary_synthesis.py b/qiskit/transpiler/passes/synthesis/unitary_synthesis.py index b6bc795bc18c..bf470068b5cd 100644 --- a/qiskit/transpiler/passes/synthesis/unitary_synthesis.py +++ b/qiskit/transpiler/passes/synthesis/unitary_synthesis.py @@ -23,8 +23,16 @@ from qiskit.dagcircuit.dagcircuit import DAGCircuit from qiskit.extensions.quantum_initializer import isometry from qiskit.quantum_info.synthesis import one_qubit_decompose +from qiskit.quantum_info.synthesis.xx_decompose import XXDecomposer from qiskit.quantum_info.synthesis.two_qubit_decompose import TwoQubitBasisDecomposer -from qiskit.circuit.library.standard_gates import iSwapGate, CXGate, CZGate, RXXGate, ECRGate +from qiskit.circuit.library.standard_gates import ( + iSwapGate, + CXGate, + CZGate, + RXXGate, + RZXGate, + ECRGate, +) from qiskit.transpiler.passes.synthesis import plugin from qiskit.providers.models import BackendProperties @@ -38,6 +46,7 @@ def _choose_kak_gate(basis_gates): "iswap": iSwapGate(), "rxx": RXXGate(pi / 2), "ecr": ECRGate(), + "rzx": RZXGate(pi / 4), # typically pi/6 is also available } kak_gate = None @@ -78,6 +87,23 @@ def _choose_bases(basis_gates, basis_dict=None): return out_basis +def _basis_gates_to_decomposer_2q(basis_gates, pulse_optimize=None): + kak_gate = _choose_kak_gate(basis_gates) + euler_basis = _choose_euler_basis(basis_gates) + + if isinstance(kak_gate, RZXGate): + backup_optimizer = TwoQubitBasisDecomposer( + CXGate(), euler_basis=euler_basis, pulse_optimize=pulse_optimize + ) + return XXDecomposer(euler_basis=euler_basis, backup_optimizer=backup_optimizer) + elif kak_gate is not None: + return TwoQubitBasisDecomposer( + kak_gate, euler_basis=euler_basis, pulse_optimize=pulse_optimize + ) + else: + return None + + class UnitarySynthesis(TransformationPass): """Synthesize gates according to their basis gates.""" @@ -332,15 +358,12 @@ def run(self, unitary, **options): qubits = options["coupling_map"][1] euler_basis = _choose_euler_basis(basis_gates) - kak_gate = _choose_kak_gate(basis_gates) - - decomposer1q, decomposer2q = None, None if euler_basis is not None: decomposer1q = one_qubit_decompose.OneQubitEulerDecomposer(euler_basis) - if kak_gate is not None: - decomposer2q = TwoQubitBasisDecomposer( - kak_gate, euler_basis=euler_basis, pulse_optimize=pulse_optimize - ) + else: + decomposer1q = None + + decomposer2q = _basis_gates_to_decomposer_2q(basis_gates, pulse_optimize=pulse_optimize) synth_dag = None wires = None diff --git a/qiskit/transpiler/passes/utils/gate_direction.py b/qiskit/transpiler/passes/utils/gate_direction.py index f617537e30cf..c6fe65de25b8 100644 --- a/qiskit/transpiler/passes/utils/gate_direction.py +++ b/qiskit/transpiler/passes/utils/gate_direction.py @@ -20,7 +20,7 @@ from qiskit.circuit import QuantumRegister from qiskit.dagcircuit import DAGCircuit -from qiskit.circuit.library.standard_gates import RYGate, HGate, CXGate, ECRGate +from qiskit.circuit.library.standard_gates import RYGate, HGate, CXGate, ECRGate, RZXGate class GateDirection(TransformationPass): @@ -39,6 +39,12 @@ class GateDirection(TransformationPass): │ ECR │ = └┬──────────┤│ ECR │├───┤ q_1: ┤1 ├ q_1: ─┤ RY(pi/2) ├┤0 ├┤ H ├ └──────┘ └──────────┘└──────┘└───┘ + + ┌──────┐ ┌───┐┌──────┐┌───┐ + q_0: ┤0 ├ q_0: ┤ H ├┤1 ├┤ H ├ + │ RZX │ = ├───┤│ RZX │├───┤ + q_1: ┤1 ├ q_1: ┤ H ├┤0 ├┤ H ├ + └──────┘ └───┘└──────┘└───┘ """ def __init__(self, coupling_map): @@ -69,6 +75,18 @@ def __init__(self, coupling_map): self._ecr_dag.apply_operation_back(HGate(), [qr[0]], []) self._ecr_dag.apply_operation_back(HGate(), [qr[1]], []) + @staticmethod + def _rzx_dag(parameter): + _rzx_dag = DAGCircuit() + qr = QuantumRegister(2) + _rzx_dag.add_qreg(qr) + _rzx_dag.apply_operation_back(HGate(), [qr[0]], []) + _rzx_dag.apply_operation_back(HGate(), [qr[1]], []) + _rzx_dag.apply_operation_back(RZXGate(parameter), [qr[1], qr[0]], []) + _rzx_dag.apply_operation_back(HGate(), [qr[0]], []) + _rzx_dag.apply_operation_back(HGate(), [qr[1]], []) + return _rzx_dag + def run(self, dag): """Run the GateDirection pass on `dag`. @@ -113,10 +131,12 @@ def run(self, dag): dag.substitute_node_with_dag(node, self._cx_dag) elif node.name == "ecr": dag.substitute_node_with_dag(node, self._ecr_dag) + elif node.name == "rzx": + dag.substitute_node_with_dag(node, self._rzx_dag(*node.op.params)) else: raise TranspilerError( - "Flipping of gate direction is only supported " - "for CX and ECR at this time." + f"Flipping of gate direction is only supported " + f"for CX, ECR, and RZX at this time, not {node.name}." ) return dag diff --git a/releasenotes/notes/feature-rzx-decomposition-c3b5a36b88303c1f.yaml b/releasenotes/notes/feature-rzx-decomposition-c3b5a36b88303c1f.yaml new file mode 100644 index 000000000000..e0e9cbde4b89 --- /dev/null +++ b/releasenotes/notes/feature-rzx-decomposition-c3b5a36b88303c1f.yaml @@ -0,0 +1,8 @@ +--- +features: + - | + This release introduces a decomposition method for two-qubit gates which + targets user-defined sets of RZX gates. Transpiler users can enable + decomposition for {RZX(pi/2), RZX(pi/4), and RZX(pi/6)} specifically by including + `'rzx'` in their `basis_gates` list. Quantum information package users can + find the method itself under the `XXDecomposer` class. diff --git a/requirements.txt b/requirements.txt index 866111a812b3..b73061f466a1 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,4 +1,3 @@ -contextvars>=2.4;python_version<'3.7' retworkx>=0.10.1 numpy>=1.17 ply>=3.10 @@ -11,3 +10,7 @@ python-dateutil>=2.8.0 stevedore>=3.0.0 symengine>=0.8 ; platform_machine == 'x86_64' or platform_machine == 'aarch64' or platform_machine == 'ppc64le' or platform_machine == 'amd64' or platform_machine == 'arm64' tweedledum>=1.1,<2.0 + +# remove these once support for python 3.6 is dropped +dataclasses>=0.8;python_version<'3.7' +contextvars>=2.4;python_version<'3.7' diff --git a/test/python/quantum_info/xx_decompose/__init__.py b/test/python/quantum_info/xx_decompose/__init__.py new file mode 100644 index 000000000000..a447510cec4a --- /dev/null +++ b/test/python/quantum_info/xx_decompose/__init__.py @@ -0,0 +1,13 @@ +# 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. + +"""Qiskit integration tests for XX-based synthesis.""" diff --git a/test/python/quantum_info/xx_decompose/test_circuits.py b/test/python/quantum_info/xx_decompose/test_circuits.py new file mode 100644 index 000000000000..2ee86fde94d2 --- /dev/null +++ b/test/python/quantum_info/xx_decompose/test_circuits.py @@ -0,0 +1,140 @@ +# 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. + +""" +Tests for qiskit-terra/qiskit/quantum_info/synthesis/xx_decompose/circuits.py . +""" + +from operator import itemgetter +import unittest + +import ddt +import numpy as np + +from qiskit.circuit import QuantumCircuit +from qiskit.circuit.library import RZGate +import qiskit.quantum_info.operators +from qiskit.quantum_info.synthesis.weyl import weyl_coordinates +from qiskit.quantum_info.synthesis.xx_decompose.circuits import ( + decompose_xxyy_into_xxyy_xx, + xx_circuit_step, +) + +from .utilities import canonical_matrix + +EPSILON = 0.001 + + +@ddt.ddt +class TestMonodromyCircuits(unittest.TestCase): + """Check circuit synthesis step routines.""" + + def __init__(self, *args, seed=42, **kwargs): + super().__init__(*args, **kwargs) + self.rng = np.random.default_rng(seed) + + def _generate_xxyy_test_case(self): + """ + Generates a random pentuple of values (a_source, b_source), beta, (a_target, b_target) s.t. + + CAN(a_source, b_source) * exp(-i(s ZI + t IZ)) * CAN(beta) + = + exp(-i(u ZI + v IZ)) * CAN(a_target, b_target) * exp(-i(x ZI + y IZ)) + + admits a solution in (s, t, u, v, x, y). + + Returns (source_coordinate, interaction, target_coordinate). + """ + source_coordinate = [self.rng.random(), self.rng.random(), 0.0] + source_coordinate = [ + source_coordinate[0] * np.pi / 8, + source_coordinate[1] * source_coordinate[0] * np.pi / 8, + 0.0, + ] + interaction = [self.rng.random() * np.pi / 8] + z_angles = [self.rng.random() * np.pi / 8, self.rng.random() * np.pi / 8] + prod = ( + canonical_matrix(*source_coordinate) + @ np.kron(RZGate(2 * z_angles[0]).to_matrix(), RZGate(2 * z_angles[1]).to_matrix()) + @ canonical_matrix(interaction[0], 0.0, 0.0) + ) + target_coordinate = weyl_coordinates(prod) + + self.assertAlmostEqual(target_coordinate[-1], 0.0, delta=EPSILON) + + return source_coordinate, interaction, target_coordinate + + # pylint:disable=invalid-name + def test_decompose_xxyy(self): + """ + Test that decompose_xxyy_into_xxyy_xx correctly recovers decompositions. + """ + + for _ in range(100): + source_coordinate, interaction, target_coordinate = self._generate_xxyy_test_case() + + r, s, u, v, x, y = decompose_xxyy_into_xxyy_xx( + target_coordinate[0], + target_coordinate[1], + source_coordinate[0], + source_coordinate[1], + interaction[0], + ) + + prod = ( + np.kron(RZGate(2 * r).to_matrix(), RZGate(2 * s).to_matrix()) + @ canonical_matrix(*source_coordinate) + @ np.kron(RZGate(2 * u).to_matrix(), RZGate(2 * v).to_matrix()) + @ canonical_matrix(interaction[0], 0.0, 0.0) + @ np.kron(RZGate(2 * x).to_matrix(), RZGate(2 * y).to_matrix()) + ) + expected = canonical_matrix(*target_coordinate) + self.assertTrue(np.all(np.abs(prod - expected) < EPSILON)) + + def test_xx_circuit_step(self): + """ + Test that `xx_circuit_step` correctly generates prefix/affix circuits relating source + canonical coordinates to target canonical coordinates along prescribed interactions, all + randomly selected. + """ + + for _ in range(100): + source_coordinate, interaction, target_coordinate = self._generate_xxyy_test_case() + + source_embodiment = qiskit.QuantumCircuit(2) + source_embodiment.append( + qiskit.extensions.UnitaryGate(canonical_matrix(*source_coordinate)), [0, 1] + ) + interaction_embodiment = qiskit.QuantumCircuit(2) + interaction_embodiment.append( + qiskit.extensions.UnitaryGate(canonical_matrix(*interaction)), [0, 1] + ) + + prefix_circuit, affix_circuit = itemgetter("prefix_circuit", "affix_circuit")( + xx_circuit_step( + source_coordinate, interaction[0], target_coordinate, interaction_embodiment + ) + ) + + target_embodiment = QuantumCircuit(2) + target_embodiment.compose(prefix_circuit, inplace=True) + target_embodiment.compose(source_embodiment, inplace=True) + target_embodiment.compose(affix_circuit, inplace=True) + self.assertTrue( + np.all( + np.abs( + qiskit.quantum_info.operators.Operator(target_embodiment).data + - canonical_matrix(*target_coordinate) + ) + < EPSILON + ) + ) diff --git a/test/python/quantum_info/xx_decompose/test_decomposer.py b/test/python/quantum_info/xx_decompose/test_decomposer.py new file mode 100644 index 000000000000..16ef3e9eab89 --- /dev/null +++ b/test/python/quantum_info/xx_decompose/test_decomposer.py @@ -0,0 +1,134 @@ +# 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. + +""" +Tests for qiskit-terra/qiskit/quantum_info/synthesis/xx_decompose/qiskit.py . +""" + +from statistics import mean +import unittest + +import ddt +import numpy as np +from scipy.stats import unitary_group + +import qiskit +from qiskit.quantum_info.operators import Operator +from qiskit.quantum_info.synthesis.xx_decompose.decomposer import ( + XXDecomposer, + TwoQubitWeylDecomposition, +) + +from .utilities import canonical_matrix + +EPSILON = 1e-8 + + +@ddt.ddt +class TestMonodromyQISKit(unittest.TestCase): + """Check QISKit routines.""" + + decomposer = XXDecomposer(euler_basis="PSX") + + def __init__(self, *args, seed=42, **kwargs): + super().__init__(*args, **kwargs) + self._random_state = np.random.Generator(np.random.PCG64(seed)) + + def test_random_compilation(self): + """Test that compilation gives correct results.""" + for _ in range(100): + unitary = unitary_group.rvs(4, random_state=self._random_state) + unitary /= np.linalg.det(unitary) ** (1 / 4) + + # decompose into CX, CX/2, and CX/3 + circuit = self.decomposer(unitary, approximate=False) + decomposed_unitary = Operator(circuit).data + + self.assertTrue(np.all(unitary - decomposed_unitary < EPSILON)) + + def test_compilation_determinism(self): + """Test that compilation is stable under multiple calls.""" + for _ in range(10): + unitary = unitary_group.rvs(4, random_state=self._random_state) + unitary /= np.linalg.det(unitary) ** (1 / 4) + + # decompose into CX, CX/2, and CX/3 + circuit1 = self.decomposer(unitary, approximate=False) + circuit2 = self.decomposer(unitary, approximate=False) + + self.assertEqual(circuit1, circuit2) + + @ddt.data(np.pi / 3, np.pi / 5, np.pi / 2) + def test_default_embodiment(self, angle): + """Test that _default_embodiment actually does yield XX gates.""" + embodiment = self.decomposer._default_embodiment(angle) + embodiment_matrix = Operator(embodiment).data + self.assertTrue(np.all(canonical_matrix(angle, 0, 0) - embodiment_matrix < EPSILON)) + + def test_check_embodiment(self): + """Test that XXDecomposer._check_embodiments correctly diagnoses il/legal embodiments.""" + + # build the member of the XX family corresponding to a single CX + good_angle = np.pi / 2 + good_embodiment = qiskit.QuantumCircuit(2) + good_embodiment.h(0) + good_embodiment.cx(0, 1) + good_embodiment.h(1) + good_embodiment.rz(np.pi / 2, 0) + good_embodiment.rz(np.pi / 2, 1) + good_embodiment.h(1) + good_embodiment.h(0) + good_embodiment.global_phase += np.pi / 4 + + # mismatch two members of the XX family + bad_angle = np.pi / 10 + bad_embodiment = qiskit.QuantumCircuit(2) + + # "self.assertDoesNotRaise" + XXDecomposer(embodiments={good_angle: good_embodiment}) + + self.assertRaises( + qiskit.exceptions.QiskitError, XXDecomposer, embodiments={bad_angle: bad_embodiment} + ) + + def test_compilation_improvement(self): + """Test that compilation to CX, CX/2, CX/3 improves over CX alone.""" + slope, offset = (64 * 90) / 1000000, 909 / 1000000 + 1 / 1000 + strength_table = self.decomposer._strength_to_infidelity( + basis_fidelity={ + strength: 1 - (slope * strength / (np.pi / 2) + offset) + for strength in [np.pi / 2, np.pi / 4, np.pi / 6] + }, + approximate=True, + ) + limited_strength_table = {np.pi / 2: strength_table[np.pi / 2]} + + clever_costs = [] + naive_costs = [] + for _ in range(200): + unitary = unitary_group.rvs(4, random_state=self._random_state) + unitary /= np.linalg.det(unitary) ** (1 / 4) + + weyl_decomposition = TwoQubitWeylDecomposition(unitary) + target = [getattr(weyl_decomposition, x) for x in ("a", "b", "c")] + if target[-1] < -EPSILON: + target = [np.pi / 2 - target[0], target[1], -target[2]] + + # decompose into CX, CX/2, and CX/3 + clever_costs.append(self.decomposer._best_decomposition(target, strength_table)["cost"]) + naive_costs.append( + self.decomposer._best_decomposition(target, limited_strength_table)["cost"] + ) + + # the following are taken from Fig 14 of the XX synthesis paper + self.assertAlmostEqual(mean(clever_costs), 1.445e-2, delta=5e-3) + self.assertAlmostEqual(mean(naive_costs), 2.058e-2, delta=5e-3) diff --git a/test/python/quantum_info/xx_decompose/test_polytopes.py b/test/python/quantum_info/xx_decompose/test_polytopes.py new file mode 100644 index 000000000000..20f3bc6ea65f --- /dev/null +++ b/test/python/quantum_info/xx_decompose/test_polytopes.py @@ -0,0 +1,61 @@ +# 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. + +""" +Tests for qiskit-terra/qiskit/quantum_info/synthesis/xx_decompose/xx_polytope.py . +""" + +import random +import unittest + +import ddt +import numpy as np +from numpy import pi + +from qiskit.quantum_info.synthesis.xx_decompose.polytopes import XXPolytope + +EPSILON = 0.001 + + +@ddt.ddt +class TestMonodromyXXPolytope(unittest.TestCase): + """Check specialized XX polytope routines.""" + + def __init__(self, *args, seed=42, **kwargs): + super().__init__(*args, **kwargs) + random.seed(seed) + np.random.seed(seed) + + @ddt.data( + ([pi / 4, pi / 4, pi / 4], [0.52359878, 0.39269908, 0.31415927]), + ([pi / 6, pi / 8, pi / 10], [pi / 6, pi / 8, pi / 10]), + ([pi / 4, pi / 4, 0], [0.615228561, 0.615228561, 0]), + ([pi / 4, pi / 8, 0], [pi / 4, pi / 8, 0]), + ) + @ddt.unpack + def test_nearest(self, offbody, expected): + """Check that the nearest point calculator recovers some known cases.""" + polytope = XXPolytope.from_strengths(pi / 6, pi / 8, pi / 10) + result = polytope.nearest(np.array(offbody)) + print(offbody, result) + self.assertTrue(np.all(np.abs(np.array(expected) - result) < EPSILON)) + + def test_add_strengths(self): + """ + Tests that adding new strengths to an existing XXPolytope is equivalent to forming the + appropriate XXPolytope from scratch. + """ + for _ in range(100): + strengths = [random.random() for _ in range(4)] + small_polytope = XXPolytope.from_strengths(*strengths[:-1]) + large_polytope = XXPolytope.from_strengths(*strengths) + self.assertEqual(small_polytope.add_strength(strengths[-1]), large_polytope) diff --git a/test/python/quantum_info/xx_decompose/test_weyl.py b/test/python/quantum_info/xx_decompose/test_weyl.py new file mode 100644 index 000000000000..3d8a915d82f4 --- /dev/null +++ b/test/python/quantum_info/xx_decompose/test_weyl.py @@ -0,0 +1,111 @@ +# 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. + +""" +Tests for qiskit-terra/qiskit/quantum_info/synthesis/xx_decompose/weyl.py . +""" + +from itertools import permutations +import unittest + +import ddt +import numpy as np + +from qiskit.quantum_info.operators import Operator +from qiskit.quantum_info.synthesis.xx_decompose.weyl import ( + apply_reflection, + apply_shift, + canonical_rotation_circuit, + reflection_options, + shift_options, +) + +from .utilities import canonical_matrix + +EPSILON = 0.001 + + +@ddt.ddt +class TestMonodromyWeyl(unittest.TestCase): + """Check Weyl action routines.""" + + def __init__(self, *args, seed=42, **kwargs): + super().__init__(*args, **kwargs) + self.rng = np.random.default_rng(seed) + + def test_reflections(self): + """Check that reflection circuits behave as expected.""" + for name in reflection_options: + coordinate = [self.rng.random() for _ in range(3)] + reflected_coordinate, reflection_circuit, reflection_phase = apply_reflection( + name, coordinate + ) + original_matrix = canonical_matrix(*coordinate) + reflected_matrix = canonical_matrix(*reflected_coordinate) + reflect_matrix = Operator(reflection_circuit).data + with np.printoptions(precision=3, suppress=True): + print(name) + print(reflect_matrix.conjugate().transpose(1, 0) @ original_matrix @ reflect_matrix) + print(reflected_matrix * reflection_phase) + self.assertTrue( + np.all( + np.abs( + reflect_matrix.conjugate().transpose(1, 0) + @ original_matrix + @ reflect_matrix + - reflected_matrix * reflection_phase + ) + < EPSILON + ) + ) + + def test_shifts(self): + """Check that shift circuits behave as expected.""" + for name in shift_options: + coordinate = [self.rng.random() for _ in range(3)] + shifted_coordinate, shift_circuit, shift_phase = apply_shift(name, coordinate) + original_matrix = canonical_matrix(*coordinate) + shifted_matrix = canonical_matrix(*shifted_coordinate) + shift_matrix = Operator(shift_circuit).data + with np.printoptions(precision=3, suppress=True): + print(name) + print(original_matrix @ shift_matrix) + print(shifted_matrix * shift_phase) + self.assertTrue( + np.all( + np.abs(original_matrix @ shift_matrix - shifted_matrix * shift_phase) < EPSILON + ) + ) + + def test_rotations(self): + """Check that rotation circuits behave as expected.""" + for permutation in permutations([0, 1, 2]): + coordinate = [self.rng.random() for _ in range(3)] + rotation_circuit = canonical_rotation_circuit(permutation[0], permutation[1]) + original_matrix = canonical_matrix(*coordinate) + rotation_matrix = Operator(rotation_circuit).data + rotated_matrix = canonical_matrix( + coordinate[permutation[0]], + coordinate[permutation[1]], + coordinate[permutation[2]], + ) + self.assertTrue( + np.all( + np.abs( + rotation_matrix.conjugate().transpose(1, 0) + @ original_matrix + @ rotation_matrix + - rotated_matrix + ) + < EPSILON + ) + ) diff --git a/test/python/quantum_info/xx_decompose/utilities.py b/test/python/quantum_info/xx_decompose/utilities.py new file mode 100644 index 000000000000..f4d97a2b9bf3 --- /dev/null +++ b/test/python/quantum_info/xx_decompose/utilities.py @@ -0,0 +1,26 @@ +# 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. + +""" +Helper functions for the XXDecomposer test suite. +""" + +from qiskit.circuit.library import RXXGate, RYYGate, RZZGate + + +def canonical_matrix(a=0.0, b=0.0, c=0.0): + """ + Produces the matrix form of a "canonical operator" + + exp(-i (a XX + b YY + c ZZ)) . + """ + return RXXGate(2 * a).to_matrix() @ RYYGate(2 * b).to_matrix() @ RZZGate(-2 * c).to_matrix()