Skip to content

Commit

Permalink
Oxidize parameter calculation for OneQubitEulerDecomposer (Qiskit#9185)
Browse files Browse the repository at this point in the history
* Oxidize parameter calculation for OneQubitEulerDecomposer

This commit ports the per basis parameter calculation from python to
rust. In looking at the runtime performance regression caused by Qiskit#8917
the majority is just that we're doing more work synthesizing to more
available basis to potentially produce better quality results. Profiling
the transpiler pass shows we're spending a non-trivial amount of time in
numpy/scipy (depending on whether it's before or after Qiskit#9179) computing
the determinant of the unitary. This is likely because those determinant
functions are designed to work with an arbitrarily large square matrix
while for the 1 qubit decomposer we're only ever working with a 2x2.
To remove this overhead this commit writes a dedicated rust function to
compute the determinant of a 2x2 complex matrix and then also adds
dedicated functions to calculate the angles for given basis to rust
as we can easily just return the end result from rust.

Related Qiskit#8774

* Eliminate python function for staticmethod definition

This commit removes one layer of function calls for the
OneQubitEulerDecomposer's staticmethods for calculating parameters.
Previously, there was a python function which called an inner rust
function, this eliminates one layer and attaches the rust function
directly as a static method to the class definition.
  • Loading branch information
mtreinish authored and Cryoris committed Jan 12, 2023
1 parent d77e804 commit 09aeb2d
Show file tree
Hide file tree
Showing 4 changed files with 130 additions and 69 deletions.
3 changes: 3 additions & 0 deletions qiskit/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,9 @@
sys.modules["qiskit._accelerate.sampled_exp_val"] = qiskit._accelerate.sampled_exp_val
sys.modules["qiskit._accelerate.vf2_layout"] = qiskit._accelerate.vf2_layout
sys.modules["qiskit._accelerate.error_map"] = qiskit._accelerate.error_map
sys.modules[
"qiskit._accelerate.euler_one_qubit_decomposer"
] = qiskit._accelerate.euler_one_qubit_decomposer


# Extend namespace for backwards compat
Expand Down
74 changes: 5 additions & 69 deletions qiskit/quantum_info/synthesis/one_qubit_decompose.py
Original file line number Diff line number Diff line change
Expand Up @@ -14,8 +14,6 @@
Decompose a single-qubit unitary via Euler angles.
"""

import math
import cmath
import numpy as np

from qiskit.circuit.quantumcircuit import QuantumCircuit
Expand All @@ -35,6 +33,7 @@
)
from qiskit.exceptions import QiskitError
from qiskit.quantum_info.operators.predicates import is_unitary_matrix
from qiskit._accelerate import euler_one_qubit_decomposer

DEFAULT_ATOL = 1e-12

Expand Down Expand Up @@ -217,62 +216,10 @@ def angles_and_phase(self, unitary):
"""
return self._params(unitary)

@staticmethod
def _params_zyz(mat):
"""Return the Euler angles and phase for the ZYZ basis."""
# We rescale the input matrix to be special unitary (det(U) = 1)
# This ensures that the quaternion representation is real
coeff = np.linalg.det(mat) ** (-0.5)
phase = -cmath.phase(coeff)
su_mat = coeff * mat # U in SU(2)
# OpenQASM SU(2) parameterization:
# U[0, 0] = exp(-i(phi+lambda)/2) * cos(theta/2)
# U[0, 1] = -exp(-i(phi-lambda)/2) * sin(theta/2)
# U[1, 0] = exp(i(phi-lambda)/2) * sin(theta/2)
# U[1, 1] = exp(i(phi+lambda)/2) * cos(theta/2)
theta = 2 * math.atan2(abs(su_mat[1, 0]), abs(su_mat[0, 0]))
phiplambda2 = cmath.phase(su_mat[1, 1])
phimlambda2 = cmath.phase(su_mat[1, 0])
phi = phiplambda2 + phimlambda2
lam = phiplambda2 - phimlambda2
return theta, phi, lam, phase

@staticmethod
def _params_zxz(mat):
"""Return the Euler angles and phase for the ZXZ basis."""
theta, phi, lam, phase = OneQubitEulerDecomposer._params_zyz(mat)
return theta, phi + np.pi / 2, lam - np.pi / 2, phase

@staticmethod
def _params_xyx(mat):
"""Return the Euler angles and phase for the XYX basis."""
# We use the fact that
# Rx(a).Ry(b).Rx(c) = H.Rz(a).Ry(-b).Rz(c).H
mat_zyz = 0.5 * np.array(
[
[
mat[0, 0] + mat[0, 1] + mat[1, 0] + mat[1, 1],
mat[0, 0] - mat[0, 1] + mat[1, 0] - mat[1, 1],
],
[
mat[0, 0] + mat[0, 1] - mat[1, 0] - mat[1, 1],
mat[0, 0] - mat[0, 1] - mat[1, 0] + mat[1, 1],
],
],
dtype=complex,
)
theta, phi, lam, phase = OneQubitEulerDecomposer._params_zyz(mat_zyz)
newphi, newlam = _mod_2pi(phi + np.pi), _mod_2pi(lam + np.pi)
return theta, newphi, newlam, phase + (newphi + newlam - phi - lam) / 2

@staticmethod
def _params_xzx(umat):
det = np.linalg.det(umat)
phase = (-1j * np.log(det)).real / 2
mat = umat / np.sqrt(det)
mat_zxz = _h_conjugate(mat)
theta, phi, lam, phase_zxz = OneQubitEulerDecomposer._params_zxz(mat_zxz)
return theta, phi, lam, phase + phase_zxz
_params_zyz = staticmethod(euler_one_qubit_decomposer.params_zyz)
_params_zxz = staticmethod(euler_one_qubit_decomposer.params_zxz)
_params_xyx = staticmethod(euler_one_qubit_decomposer.params_xyx)
_params_xzx = staticmethod(euler_one_qubit_decomposer.params_xzx)

@staticmethod
def _params_u3(mat):
Expand Down Expand Up @@ -593,14 +540,3 @@ def _mod_2pi(angle: float, atol: float = 0):
if abs(wrapped - np.pi) < atol:
wrapped = -np.pi
return wrapped


def _h_conjugate(su2):
"""Return su2 conjugated by Hadamard gate. No warning if input matrix is not in su2."""
return np.array(
[
[su2[0, 0].real + 1j * su2[1, 0].imag, 1j * su2[0, 0].imag + su2[1, 0].real],
[1j * su2[0, 0].imag - su2[1, 0].real, su2[0, 0].real - 1j * su2[1, 0].imag],
],
dtype=complex,
)
118 changes: 118 additions & 0 deletions src/euler_one_qubit_decomposer.rs
Original file line number Diff line number Diff line change
@@ -0,0 +1,118 @@
// This code is part of Qiskit.
//
// (C) Copyright IBM 2022
//
// This code is licensed under the Apache License, Version 2.0. You may
// obtain a copy of this license in the LICENSE.txt file in the root directory
// of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
//
// Any modifications or derivative works of this code must retain this
// copyright notice, and modified files need to carry a notice indicating
// that they have been altered from the originals.

use num_complex::{Complex64, ComplexFloat};
use pyo3::prelude::*;
use pyo3::wrap_pyfunction;
use pyo3::Python;
use std::f64::consts::PI;

use ndarray::prelude::*;
use numpy::PyReadonlyArray2;

#[inline]
fn det_one_qubit(mat: ArrayView2<Complex64>) -> Complex64 {
mat[[0, 0]] * mat[[1, 1]] - mat[[0, 1]] * mat[[1, 0]]
}

#[inline]
fn complex_phase(x: Complex64) -> f64 {
x.im.atan2(x.re)
}

#[inline]
fn mod_2pi(angle: f64) -> f64 {
(angle + PI) % (2. * PI) - PI
}

fn params_zyz_inner(mat: ArrayView2<Complex64>) -> [f64; 4] {
let coeff: Complex64 = 1. / det_one_qubit(mat).sqrt();
let phase = -complex_phase(coeff);
let tmp_1_0 = (coeff * mat[[1, 0]]).abs();
let tmp_0_0 = (coeff * mat[[0, 0]]).abs();
let theta = 2. * tmp_1_0.atan2(tmp_0_0);
let phiplambda2 = complex_phase(coeff * mat[[1, 1]]);
let phimlambda2 = complex_phase(coeff * mat[[1, 0]]);
let phi = phiplambda2 + phimlambda2;
let lam = phiplambda2 - phimlambda2;
[theta, phi, lam, phase]
}

fn params_zxz_inner(mat: ArrayView2<Complex64>) -> [f64; 4] {
let [theta, phi, lam, phase] = params_zyz_inner(mat);
[theta, phi + PI / 2., lam - PI / 2., phase]
}

#[pyfunction]
pub fn params_zyz(unitary: PyReadonlyArray2<Complex64>) -> [f64; 4] {
let mat = unitary.as_array();
params_zyz_inner(mat)
}

#[pyfunction]
pub fn params_xyx(unitary: PyReadonlyArray2<Complex64>) -> [f64; 4] {
let mat = unitary.as_array();
let mat_zyz = arr2(&[
[
0.5 * (mat[[0, 0]] + mat[[0, 1]] + mat[[1, 0]] + mat[[1, 1]]),
0.5 * (mat[[0, 0]] - mat[[0, 1]] + mat[[1, 0]] - mat[[1, 1]]),
],
[
0.5 * (mat[[0, 0]] + mat[[0, 1]] - mat[[1, 0]] - mat[[1, 1]]),
0.5 * (mat[[0, 0]] - mat[[0, 1]] - mat[[1, 0]] + mat[[1, 1]]),
],
]);
let [theta, phi, lam, phase] = params_zyz_inner(mat_zyz.view());
let new_phi = mod_2pi(phi + PI);
let new_lam = mod_2pi(lam + PI);
[
theta,
new_phi,
new_lam,
phase + (new_phi + new_lam - phi - lam) / 2.,
]
}

#[pyfunction]
pub fn params_xzx(unitary: PyReadonlyArray2<Complex64>) -> [f64; 4] {
let umat = unitary.as_array();
let det = det_one_qubit(umat);
let phase = (Complex64::new(0., -1.) * det.ln()).re / 2.;
let sqrt_det = det.sqrt();
let mat_zyz = arr2(&[
[
Complex64::new((umat[[0, 0]] / sqrt_det).re, (umat[[1, 0]] / sqrt_det).im),
Complex64::new((umat[[1, 0]] / sqrt_det).re, (umat[[0, 0]] / sqrt_det).im),
],
[
Complex64::new(-(umat[[1, 0]] / sqrt_det).re, (umat[[0, 0]] / sqrt_det).im),
Complex64::new((umat[[0, 0]] / sqrt_det).re, -(umat[[1, 0]] / sqrt_det).im),
],
]);
let [theta, phi, lam, phase_zxz] = params_zxz_inner(mat_zyz.view());
[theta, phi, lam, phase + phase_zxz]
}

#[pyfunction]
pub fn params_zxz(unitary: PyReadonlyArray2<Complex64>) -> [f64; 4] {
let mat = unitary.as_array();
params_zxz_inner(mat)
}

#[pymodule]
pub fn euler_one_qubit_decomposer(_py: Python, m: &PyModule) -> PyResult<()> {
m.add_wrapped(wrap_pyfunction!(params_zyz))?;
m.add_wrapped(wrap_pyfunction!(params_xyx))?;
m.add_wrapped(wrap_pyfunction!(params_xzx))?;
m.add_wrapped(wrap_pyfunction!(params_zxz))?;
Ok(())
}
4 changes: 4 additions & 0 deletions src/lib.rs
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@ use pyo3::Python;
mod dense_layout;
mod edge_collections;
mod error_map;
mod euler_one_qubit_decomposer;
mod nlayout;
mod optimize_1q_gates;
mod pauli_exp_val;
Expand Down Expand Up @@ -55,5 +56,8 @@ fn _accelerate(_py: Python<'_>, m: &PyModule) -> PyResult<()> {
m.add_wrapped(wrap_pymodule!(optimize_1q_gates::optimize_1q_gates))?;
m.add_wrapped(wrap_pymodule!(sampled_exp_val::sampled_exp_val))?;
m.add_wrapped(wrap_pymodule!(vf2_layout::vf2_layout))?;
m.add_wrapped(wrap_pymodule!(
euler_one_qubit_decomposer::euler_one_qubit_decomposer
))?;
Ok(())
}

0 comments on commit 09aeb2d

Please sign in to comment.