diff --git a/Cargo.lock b/Cargo.lock index 6859c622967d..dee434e870a4 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -589,6 +589,15 @@ dependencies = [ "either", ] +[[package]] +name = "itertools" +version = "0.12.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ba291022dbbd398a455acf126c1e341954079855bc60dfdda641363bd6922569" +dependencies = [ + "either", +] + [[package]] name = "jod-thread" version = "0.1.2" @@ -1079,6 +1088,7 @@ dependencies = [ "faer-ext", "hashbrown 0.14.3", "indexmap 2.2.6", + "itertools 0.12.1", "ndarray", "num-bigint", "num-complex", diff --git a/crates/accelerate/Cargo.toml b/crates/accelerate/Cargo.toml index 05ca3f5b6395..a43fdc6ff506 100644 --- a/crates/accelerate/Cargo.toml +++ b/crates/accelerate/Cargo.toml @@ -21,6 +21,7 @@ num-complex = "0.4" num-bigint = "0.4" rustworkx-core = "0.14" faer = "0.18.2" +itertools = "0.12.1" qiskit-circuit.workspace = true [dependencies.smallvec] diff --git a/crates/accelerate/src/isometry.rs b/crates/accelerate/src/isometry.rs new file mode 100644 index 000000000000..8d0761666bb6 --- /dev/null +++ b/crates/accelerate/src/isometry.rs @@ -0,0 +1,367 @@ +// This code is part of Qiskit. +// +// (C) Copyright IBM 2024 +// +// 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 std::ops::BitAnd; + +use approx::abs_diff_eq; +use num_complex::{Complex64, ComplexFloat}; +use pyo3::prelude::*; +use pyo3::wrap_pyfunction; +use pyo3::Python; + +use hashbrown::HashSet; +use itertools::Itertools; +use ndarray::prelude::*; +use numpy::{IntoPyArray, PyReadonlyArray1, PyReadonlyArray2}; + +use crate::two_qubit_decompose::ONE_QUBIT_IDENTITY; + +/// Find special unitary matrix that maps [c0,c1] to [r,0] or [0,r] if basis_state=0 or +/// basis_state=1 respectively +#[pyfunction] +pub fn reverse_qubit_state( + py: Python, + state: [Complex64; 2], + basis_state: usize, + epsilon: f64, +) -> PyObject { + reverse_qubit_state_inner(&state, basis_state, epsilon) + .into_pyarray_bound(py) + .into() +} + +#[inline(always)] +fn l2_norm(vec: &[Complex64]) -> f64 { + vec.iter() + .fold(0., |acc, elem| acc + elem.norm_sqr()) + .sqrt() +} + +fn reverse_qubit_state_inner( + state: &[Complex64; 2], + basis_state: usize, + epsilon: f64, +) -> Array2 { + let r = l2_norm(state); + let r_inv = 1. / r; + if r < epsilon { + Array2::eye(2) + } else if basis_state == 0 { + array![ + [state[0].conj() * r_inv, state[1].conj() * r_inv], + [-state[1] * r_inv, state[0] * r_inv], + ] + } else { + array![ + [-state[1] * r_inv, state[0] * r_inv], + [state[0].conj() * r_inv, state[1].conj() * r_inv], + ] + } +} + +/// This method finds the single-qubit gates for a UCGate to disentangle a qubit: +/// we consider the n-qubit state v[:,0] starting with k zeros (in the computational basis). +/// The qubit with label n-s-1 is disentangled into the basis state k_s(k,s). + +#[pyfunction] +pub fn find_squs_for_disentangling( + py: Python, + v: PyReadonlyArray2, + k: usize, + s: usize, + epsilon: f64, + n: usize, +) -> Vec { + let v = v.as_array(); + let k_prime = 0; + let i_start = if b(k, s + 1) == 0 { + a(k, s + 1) + } else { + a(k, s + 1) + 1 + }; + let mut output: Vec> = (0..i_start).map(|_| Array2::eye(2)).collect(); + let mut squs: Vec> = (i_start..2_usize.pow((n - s - 1) as u32)) + .map(|i| { + reverse_qubit_state_inner( + &[ + v[[2 * i * 2_usize.pow(s as u32) + b(k, s), k_prime]], + v[[(2 * i + 1) * 2_usize.pow(s as u32) + b(k, s), k_prime]], + ], + k_s(k, s), + epsilon, + ) + }) + .collect(); + output.append(&mut squs); + output + .into_iter() + .map(|x| x.into_pyarray_bound(py).into()) + .collect() +} + +#[pyfunction] +pub fn apply_ucg( + py: Python, + m: PyReadonlyArray2, + k: usize, + single_qubit_gates: Vec>, +) -> PyObject { + let mut m = m.as_array().to_owned(); + let shape = m.shape(); + let num_qubits = shape[0].ilog2(); + let num_col = shape[1]; + let spacing: usize = 2_usize.pow(num_qubits - k as u32 - 1); + for j in 0..2_usize.pow(num_qubits - 1) { + let i = (j / spacing) * spacing + j; + let gate_index = i / (2_usize.pow(num_qubits - k as u32)); + for col in 0..num_col { + let gate = single_qubit_gates[gate_index].as_array(); + let a = m[[i, col]]; + let b = m[[i + spacing, col]]; + m[[i, col]] = gate[[0, 0]] * a + gate[[0, 1]] * b; + m[[i + spacing, col]] = gate[[1, 0]] * a + gate[[1, 1]] * b; + } + } + m.into_pyarray_bound(py).into() +} + +#[inline(always)] +fn bin_to_int(bin: &[u8]) -> usize { + bin.iter() + .fold(0_usize, |acc, digit| (acc << 1) + *digit as usize) +} + +#[pyfunction] +pub fn apply_diagonal_gate( + py: Python, + m: PyReadonlyArray2, + action_qubit_labels: Vec, + diag: PyReadonlyArray1, +) -> PyResult { + let diag = diag.as_slice()?; + let mut m = m.as_array().to_owned(); + let shape = m.shape(); + let num_qubits = shape[0].ilog2(); + let num_col = shape[1]; + for state in std::iter::repeat([0_u8, 1_u8]) + .take(num_qubits as usize) + .multi_cartesian_product() + { + let diag_index = action_qubit_labels + .iter() + .fold(0_usize, |acc, i| (acc << 1) + state[*i] as usize); + let i = bin_to_int(&state); + for j in 0..num_col { + m[[i, j]] = diag[diag_index] * m[[i, j]] + } + } + Ok(m.into_pyarray_bound(py).into()) +} + +#[pyfunction] +pub fn apply_diagonal_gate_to_diag( + mut m_diagonal: Vec, + action_qubit_labels: Vec, + diag: PyReadonlyArray1, + num_qubits: usize, +) -> PyResult> { + let diag = diag.as_slice()?; + if m_diagonal.is_empty() { + return Ok(m_diagonal); + } + for state in std::iter::repeat([0_u8, 1_u8]) + .take(num_qubits) + .multi_cartesian_product() + .take(m_diagonal.len()) + { + let diag_index = action_qubit_labels + .iter() + .fold(0_usize, |acc, i| (acc << 1) + state[*i] as usize); + let i = bin_to_int(&state); + m_diagonal[i] *= diag[diag_index] + } + Ok(m_diagonal) +} + +/// Helper method for _apply_multi_controlled_gate. This constructs the basis states the MG gate +/// is acting on for a specific state state_free of the qubits we neither control nor act on +fn construct_basis_states( + state_free: &[u8], + control_set: &HashSet, + target_label: usize, +) -> [usize; 2] { + let size = state_free.len() + control_set.len() + 1; + let mut e1: usize = 0; + let mut e2: usize = 0; + let mut j = 0; + for i in 0..size { + e1 <<= 1; + e2 <<= 1; + if control_set.contains(&i) { + e1 += 1; + e2 += 1; + } else if i == target_label { + e2 += 1; + } else { + assert!(j <= 1); + e1 += state_free[j] as usize; + e2 += state_free[j] as usize; + j += 1 + } + } + [e1, e2] +} + +#[pyfunction] +pub fn apply_multi_controlled_gate( + py: Python, + m: PyReadonlyArray2, + control_labels: Vec, + target_label: usize, + gate: PyReadonlyArray2, +) -> PyObject { + let mut m = m.as_array().to_owned(); + let gate = gate.as_array(); + let shape = m.shape(); + let num_qubits = shape[0].ilog2(); + let num_col = shape[1]; + let free_qubits = num_qubits as usize - control_labels.len() - 1; + let control_set: HashSet = control_labels.into_iter().collect(); + if free_qubits == 0 { + let [e1, e2] = construct_basis_states(&[], &control_set, target_label); + for i in 0..num_col { + let temp: Vec<_> = gate + .dot(&aview2(&[[m[[e1, i]]], [m[[e2, i]]]])) + .into_iter() + .take(2) + .collect(); + m[[e1, i]] = temp[0]; + m[[e2, i]] = temp[1]; + } + return m.into_pyarray_bound(py).into(); + } + for state_free in std::iter::repeat([0_u8, 1_u8]) + .take(free_qubits) + .multi_cartesian_product() + { + let [e1, e2] = construct_basis_states(&state_free, &control_set, target_label); + for i in 0..num_col { + let temp: Vec<_> = gate + .dot(&aview2(&[[m[[e1, i]]], [m[[e2, i]]]])) + .into_iter() + .take(2) + .collect(); + m[[e1, i]] = temp[0]; + m[[e2, i]] = temp[1]; + } + } + m.into_pyarray_bound(py).into() +} + +#[pyfunction] +pub fn ucg_is_identity_up_to_global_phase( + single_qubit_gates: Vec>, + epsilon: f64, +) -> bool { + let global_phase: Complex64 = if single_qubit_gates[0].as_array()[[0, 0]].abs() >= epsilon { + single_qubit_gates[0].as_array()[[0, 0]].finv() + } else { + return false; + }; + for raw_gate in single_qubit_gates { + let gate = raw_gate.as_array(); + if !abs_diff_eq!( + gate.mapv(|x| x * global_phase), + aview2(&ONE_QUBIT_IDENTITY), + epsilon = 1e-8 // Default tolerance from numpy for allclose() + ) { + return false; + } + } + true +} + +#[pyfunction] +fn diag_is_identity_up_to_global_phase(diag: Vec, epsilon: f64) -> bool { + let global_phase: Complex64 = if diag[0].abs() >= epsilon { + diag[0].finv() + } else { + return false; + }; + for d in diag { + if (global_phase * d - 1.0).abs() >= epsilon { + return false; + } + } + true +} + +#[pyfunction] +pub fn merge_ucgate_and_diag( + py: Python, + single_qubit_gates: Vec>, + diag: Vec, +) -> Vec { + single_qubit_gates + .iter() + .enumerate() + .map(|(i, raw_gate)| { + let gate = raw_gate.as_array(); + let res = aview2(&[ + [diag[2 * i], Complex64::new(0., 0.)], + [Complex64::new(0., 0.), diag[2 * i + 1]], + ]) + .dot(&gate); + res.into_pyarray_bound(py).into() + }) + .collect() +} + +#[inline(always)] +#[pyfunction] +fn k_s(k: usize, s: usize) -> usize { + if k == 0 { + 0 + } else { + let filter = 1 << s; + k.bitand(filter) >> s + } +} + +#[inline(always)] +#[pyfunction] +fn a(k: usize, s: usize) -> usize { + k / 2_usize.pow(s as u32) +} + +#[inline(always)] +#[pyfunction] +fn b(k: usize, s: usize) -> usize { + k - (a(k, s) * 2_usize.pow(s as u32)) +} + +#[pymodule] +pub fn isometry(m: &Bound) -> PyResult<()> { + m.add_wrapped(wrap_pyfunction!(diag_is_identity_up_to_global_phase))?; + m.add_wrapped(wrap_pyfunction!(find_squs_for_disentangling))?; + m.add_wrapped(wrap_pyfunction!(reverse_qubit_state))?; + m.add_wrapped(wrap_pyfunction!(apply_ucg))?; + m.add_wrapped(wrap_pyfunction!(apply_diagonal_gate))?; + m.add_wrapped(wrap_pyfunction!(apply_diagonal_gate_to_diag))?; + m.add_wrapped(wrap_pyfunction!(apply_multi_controlled_gate))?; + m.add_wrapped(wrap_pyfunction!(ucg_is_identity_up_to_global_phase))?; + m.add_wrapped(wrap_pyfunction!(merge_ucgate_and_diag))?; + m.add_wrapped(wrap_pyfunction!(a))?; + m.add_wrapped(wrap_pyfunction!(b))?; + m.add_wrapped(wrap_pyfunction!(k_s))?; + Ok(()) +} diff --git a/crates/accelerate/src/lib.rs b/crates/accelerate/src/lib.rs index bb7621dce346..0af8ea6a0fce 100644 --- a/crates/accelerate/src/lib.rs +++ b/crates/accelerate/src/lib.rs @@ -19,6 +19,7 @@ pub mod dense_layout; pub mod edge_collections; pub mod error_map; pub mod euler_one_qubit_decomposer; +pub mod isometry; pub mod nlayout; pub mod optimize_1q_gates; pub mod pauli_exp_val; @@ -28,6 +29,7 @@ pub mod sampled_exp_val; pub mod sparse_pauli_op; pub mod stochastic_swap; pub mod two_qubit_decompose; +pub mod uc_gate; pub mod utils; pub mod vf2_layout; diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 7dcb273ac163..5e833bd86fda 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -60,7 +60,7 @@ const TWO_PI: f64 = 2.0 * PI; const C1: c64 = c64 { re: 1.0, im: 0.0 }; -static ONE_QUBIT_IDENTITY: [[Complex64; 2]; 2] = [ +pub static ONE_QUBIT_IDENTITY: [[Complex64; 2]; 2] = [ [Complex64::new(1., 0.), Complex64::new(0., 0.)], [Complex64::new(0., 0.), Complex64::new(1., 0.)], ]; diff --git a/crates/accelerate/src/uc_gate.rs b/crates/accelerate/src/uc_gate.rs new file mode 100644 index 000000000000..3a5f74a6f0b1 --- /dev/null +++ b/crates/accelerate/src/uc_gate.rs @@ -0,0 +1,163 @@ +// This code is part of Qiskit. +// +// (C) Copyright IBM 2024 +// +// 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::{FRAC_1_SQRT_2, PI}; + +use faer_ext::{IntoFaerComplex, IntoNdarrayComplex}; +use ndarray::prelude::*; +use numpy::{IntoPyArray, PyReadonlyArray2}; + +use crate::euler_one_qubit_decomposer::det_one_qubit; + +const PI2: f64 = PI / 2.; +const EPS: f64 = 1e-10; + +// These constants are the non-zero elements of an RZ gate's unitary with an +// angle of pi / 2 +const RZ_PI2_11: Complex64 = Complex64::new(FRAC_1_SQRT_2, -FRAC_1_SQRT_2); +const RZ_PI2_00: Complex64 = Complex64::new(FRAC_1_SQRT_2, FRAC_1_SQRT_2); + +/// This method implements the decomposition given in equation (3) in +/// https://arxiv.org/pdf/quant-ph/0410066.pdf. +/// +/// The decomposition is used recursively to decompose uniformly controlled gates. +/// +/// a,b = single qubit unitaries +/// v,u,r = outcome of the decomposition given in the reference mentioned above +/// +/// (see there for the details). +fn demultiplex_single_uc( + a: ArrayView2, + b: ArrayView2, +) -> [Array2; 3] { + let x = a.dot(&b.mapv(|x| x.conj()).t()); + let det_x = det_one_qubit(x.view()); + let x11 = x[[0, 0]] / det_x.sqrt(); + let phi = det_x.arg(); + + let r1 = (Complex64::new(0., 1.) / 2. * (PI2 - phi / 2. - x11.arg())).exp(); + let r2 = (Complex64::new(0., 1.) / 2. * (PI2 - phi / 2. + x11.arg() + PI)).exp(); + + let r = array![[r1, Complex64::new(0., 0.)], [Complex64::new(0., 0.), r2],]; + + let decomp = r + .dot(&x) + .dot(&r) + .view() + .into_faer_complex() + .complex_eigendecomposition(); + let mut u: Array2 = decomp.u().into_ndarray_complex().to_owned(); + let s = decomp.s().column_vector(); + let mut diag: Array1 = + Array1::from_shape_fn(u.shape()[0], |i| s[i].to_num_complex()); + + // If d is not equal to diag(i,-i), then we put it into this "standard" form + // (see eq. (13) in https://arxiv.org/pdf/quant-ph/0410066.pdf) by interchanging + // the eigenvalues and eigenvectors + if (diag[0] + Complex64::new(0., 1.)).abs() < EPS { + diag = diag.slice(s![..;-1]).to_owned(); + u = u.slice(s![.., ..;-1]).to_owned(); + } + diag.mapv_inplace(|x| x.sqrt()); + let d = Array2::from_diag(&diag); + let v = d + .dot(&u.mapv(|x| x.conj()).t()) + .dot(&r.mapv(|x| x.conj()).t()) + .dot(&b); + [v, u, r] +} + +#[pyfunction] +pub fn dec_ucg_help( + py: Python, + sq_gates: Vec>, + num_qubits: u32, +) -> (Vec, PyObject) { + let mut single_qubit_gates: Vec> = sq_gates + .into_iter() + .map(|x| x.as_array().to_owned()) + .collect(); + let mut diag: Array1 = Array1::ones(2_usize.pow(num_qubits)); + let num_controls = num_qubits - 1; + for dec_step in 0..num_controls { + let num_ucgs = 2_usize.pow(dec_step); + // The decomposition works recursively and the followign loop goes over the different + // UCGates that arise in the decomposition + for ucg_index in 0..num_ucgs { + let len_ucg = 2_usize.pow(num_controls - dec_step); + for i in 0..len_ucg / 2 { + let shift = ucg_index * len_ucg; + let a = single_qubit_gates[shift + i].view(); + let b = single_qubit_gates[shift + len_ucg / 2 + i].view(); + // Apply the decomposition for UCGates given in equation (3) in + // https://arxiv.org/pdf/quant-ph/0410066.pdf + // to demultiplex one control of all the num_ucgs uniformly-controlled gates + // with log2(len_ucg) uniform controls + let [v, u, r] = demultiplex_single_uc(a, b); + // replace the single-qubit gates with v,u (the already existing ones + // are not needed any more) + single_qubit_gates[shift + i] = v; + single_qubit_gates[shift + len_ucg / 2 + i] = u; + // Now we decompose the gates D as described in Figure 4 in + // https://arxiv.org/pdf/quant-ph/0410066.pdf and merge some of the gates + // into the UCGates and the diagonal at the end of the circuit + + // Remark: The Rz(pi/2) rotation acting on the target qubit and the Hadamard + // gates arising in the decomposition of D are ignored for the moment (they will + // be added together with the C-NOT gates at the end of the decomposition + // (in the method dec_ucg())) + let r_conj_t = r.mapv(|x| x.conj()).t().to_owned(); + if ucg_index < num_ucgs - 1 { + // Absorb the Rz(pi/2) rotation on the control into the UC-Rz gate and + // merge the UC-Rz rotation with the following UCGate, + // which hasn't been decomposed yet + let k = shift + len_ucg + i; + + single_qubit_gates[k] = single_qubit_gates[k].dot(&r_conj_t); + single_qubit_gates[k].mapv_inplace(|x| x * RZ_PI2_00); + let k = k + len_ucg / 2; + single_qubit_gates[k] = single_qubit_gates[k].dot(&r); + single_qubit_gates[k].mapv_inplace(|x| x * RZ_PI2_11); + } else { + // Absorb the Rz(pi/2) rotation on the control into the UC-Rz gate and merge + // the trailing UC-Rz rotation into a diagonal gate at the end of the circuit + for ucg_index_2 in 0..num_ucgs { + let shift_2 = ucg_index_2 * len_ucg; + let k = 2 * (i + shift_2); + diag[k] *= r_conj_t[[0, 0]] * RZ_PI2_00; + diag[k + 1] *= r_conj_t[[1, 1]] * RZ_PI2_00; + let k = len_ucg + k; + diag[k] *= r[[0, 0]] * RZ_PI2_11; + diag[k + 1] *= r[[1, 1]] * RZ_PI2_11; + } + } + } + } + } + ( + single_qubit_gates + .into_iter() + .map(|x| x.into_pyarray_bound(py).into()) + .collect(), + diag.into_pyarray_bound(py).into(), + ) +} + +#[pymodule] +pub fn uc_gate(m: &Bound) -> PyResult<()> { + m.add_wrapped(wrap_pyfunction!(dec_ucg_help))?; + Ok(()) +} diff --git a/crates/pyext/src/lib.rs b/crates/pyext/src/lib.rs index b7c89872bf8d..a21b1307a88f 100644 --- a/crates/pyext/src/lib.rs +++ b/crates/pyext/src/lib.rs @@ -15,10 +15,11 @@ use pyo3::wrap_pymodule; use qiskit_accelerate::{ convert_2q_block_matrix::convert_2q_block_matrix, dense_layout::dense_layout, - error_map::error_map, euler_one_qubit_decomposer::euler_one_qubit_decomposer, nlayout::nlayout, - optimize_1q_gates::optimize_1q_gates, pauli_exp_val::pauli_expval, results::results, - sabre::sabre, sampled_exp_val::sampled_exp_val, sparse_pauli_op::sparse_pauli_op, - stochastic_swap::stochastic_swap, two_qubit_decompose::two_qubit_decompose, utils::utils, + error_map::error_map, euler_one_qubit_decomposer::euler_one_qubit_decomposer, + isometry::isometry, nlayout::nlayout, optimize_1q_gates::optimize_1q_gates, + pauli_exp_val::pauli_expval, results::results, sabre::sabre, sampled_exp_val::sampled_exp_val, + sparse_pauli_op::sparse_pauli_op, stochastic_swap::stochastic_swap, + two_qubit_decompose::two_qubit_decompose, uc_gate::uc_gate, utils::utils, vf2_layout::vf2_layout, }; @@ -31,6 +32,7 @@ fn _accelerate(m: &Bound) -> PyResult<()> { m.add_wrapped(wrap_pymodule!(dense_layout))?; m.add_wrapped(wrap_pymodule!(error_map))?; m.add_wrapped(wrap_pymodule!(euler_one_qubit_decomposer))?; + m.add_wrapped(wrap_pymodule!(isometry))?; m.add_wrapped(wrap_pymodule!(nlayout))?; m.add_wrapped(wrap_pymodule!(optimize_1q_gates))?; m.add_wrapped(wrap_pymodule!(pauli_expval))?; @@ -40,6 +42,7 @@ fn _accelerate(m: &Bound) -> PyResult<()> { m.add_wrapped(wrap_pymodule!(sparse_pauli_op))?; m.add_wrapped(wrap_pymodule!(stochastic_swap))?; m.add_wrapped(wrap_pymodule!(two_qubit_decompose))?; + m.add_wrapped(wrap_pymodule!(uc_gate))?; m.add_wrapped(wrap_pymodule!(utils))?; m.add_wrapped(wrap_pymodule!(vf2_layout))?; Ok(()) diff --git a/qiskit/__init__.py b/qiskit/__init__.py index 590a5698a77d..e4fbc1729e53 100644 --- a/qiskit/__init__.py +++ b/qiskit/__init__.py @@ -65,6 +65,8 @@ ) sys.modules["qiskit._accelerate.dense_layout"] = qiskit._accelerate.dense_layout sys.modules["qiskit._accelerate.error_map"] = qiskit._accelerate.error_map +sys.modules["qiskit._accelerate.isometry"] = qiskit._accelerate.isometry +sys.modules["qiskit._accelerate.uc_gate"] = qiskit._accelerate.uc_gate sys.modules["qiskit._accelerate.euler_one_qubit_decomposer"] = ( qiskit._accelerate.euler_one_qubit_decomposer ) diff --git a/qiskit/circuit/library/generalized_gates/isometry.py b/qiskit/circuit/library/generalized_gates/isometry.py index 1294feb26342..c180e7a14484 100644 --- a/qiskit/circuit/library/generalized_gates/isometry.py +++ b/qiskit/circuit/library/generalized_gates/isometry.py @@ -21,7 +21,6 @@ from __future__ import annotations -import itertools import math import numpy as np from qiskit.circuit.exceptions import CircuitError @@ -30,6 +29,7 @@ from qiskit.circuit.quantumregister import QuantumRegister from qiskit.exceptions import QiskitError from qiskit.quantum_info.operators.predicates import is_isometry +from qiskit._accelerate import isometry as isometry_rs from .diagonal import Diagonal from .uc import UCGate @@ -157,12 +157,16 @@ def _gates_to_uncompute(self): # correspond to the firstfew columns of the identity matrix up to diag, and hence we only # have to save a list containing them. for column_index in range(2**m): - self._decompose_column(circuit, q, diag, remaining_isometry, column_index) + remaining_isometry, diag = self._decompose_column( + circuit, q, diag, remaining_isometry, column_index + ) # extract phase of the state that was sent to the basis state ket(column_index) diag.append(remaining_isometry[column_index, 0]) # remove first column (which is now stored in diag) remaining_isometry = remaining_isometry[:, 1:] - if len(diag) > 1 and not _diag_is_identity_up_to_global_phase(diag, self._epsilon): + if len(diag) > 1 and not isometry_rs.diag_is_identity_up_to_global_phase( + diag, self._epsilon + ): diagonal = Diagonal(np.conj(diag)) circuit.append(diagonal, q_input) return circuit @@ -173,7 +177,10 @@ def _decompose_column(self, circuit, q, diag, remaining_isometry, column_index): """ n = int(math.log2(self.iso_data.shape[0])) for s in range(n): - self._disentangle(circuit, q, diag, remaining_isometry, column_index, s) + remaining_isometry, diag = self._disentangle( + circuit, q, diag, remaining_isometry, column_index, s + ) + return remaining_isometry, diag def _disentangle(self, circuit, q, diag, remaining_isometry, column_index, s): """ @@ -189,13 +196,19 @@ def _disentangle(self, circuit, q, diag, remaining_isometry, column_index, s): n = int(math.log2(self.iso_data.shape[0])) # MCG to set one entry to zero (preparation for disentangling with UCGate): - index1 = 2 * _a(k, s + 1) * 2**s + _b(k, s + 1) - index2 = (2 * _a(k, s + 1) + 1) * 2**s + _b(k, s + 1) + index1 = 2 * isometry_rs.a(k, s + 1) * 2**s + isometry_rs.b(k, s + 1) + index2 = (2 * isometry_rs.a(k, s + 1) + 1) * 2**s + isometry_rs.b(k, s + 1) target_label = n - s - 1 # Check if a MCG is required - if _k_s(k, s) == 0 and _b(k, s + 1) != 0 and np.abs(v[index2, k_prime]) > self._epsilon: + if ( + isometry_rs.k_s(k, s) == 0 + and isometry_rs.b(k, s + 1) != 0 + and np.abs(v[index2, k_prime]) > self._epsilon + ): # Find the MCG, decompose it and apply it to the remaining isometry - gate = _reverse_qubit_state([v[index1, k_prime], v[index2, k_prime]], 0, self._epsilon) + gate = isometry_rs.reverse_qubit_state( + [v[index1, k_prime], v[index2, k_prime]], 0, self._epsilon + ) control_labels = [ i for i, x in enumerate(_get_binary_rep_as_list(k, n)) @@ -205,57 +218,49 @@ def _disentangle(self, circuit, q, diag, remaining_isometry, column_index, s): circuit, q, gate, control_labels, target_label ) # apply the MCG to the remaining isometry - _apply_multi_controlled_gate(v, control_labels, target_label, gate) + v = isometry_rs.apply_multi_controlled_gate(v, control_labels, target_label, gate) # correct for the implementation "up to diagonal" - diag_mcg_inverse = np.conj(diagonal_mcg).tolist() - _apply_diagonal_gate(v, control_labels + [target_label], diag_mcg_inverse) + diag_mcg_inverse = np.conj(diagonal_mcg).astype(complex, copy=False) + v = isometry_rs.apply_diagonal_gate( + v, control_labels + [target_label], diag_mcg_inverse + ) # update the diag according to the applied diagonal gate - _apply_diagonal_gate_to_diag(diag, control_labels + [target_label], diag_mcg_inverse, n) + diag = isometry_rs.apply_diagonal_gate_to_diag( + diag, control_labels + [target_label], diag_mcg_inverse, n + ) # UCGate to disentangle a qubit: # Find the UCGate, decompose it and apply it to the remaining isometry single_qubit_gates = self._find_squs_for_disentangling(v, k, s) - if not _ucg_is_identity_up_to_global_phase(single_qubit_gates, self._epsilon): + if not isometry_rs.ucg_is_identity_up_to_global_phase(single_qubit_gates, self._epsilon): control_labels = list(range(target_label)) diagonal_ucg = self._append_ucg_up_to_diagonal( circuit, q, single_qubit_gates, control_labels, target_label ) # merge the diagonal into the UCGate for efficient application of both together - diagonal_ucg_inverse = np.conj(diagonal_ucg).tolist() - single_qubit_gates = _merge_UCGate_and_diag(single_qubit_gates, diagonal_ucg_inverse) + diagonal_ucg_inverse = np.conj(diagonal_ucg).astype(complex, copy=False) + single_qubit_gates = isometry_rs.merge_ucgate_and_diag( + single_qubit_gates, diagonal_ucg_inverse + ) # apply the UCGate (with the merged diagonal gate) to the remaining isometry - _apply_ucg(v, len(control_labels), single_qubit_gates) + v = isometry_rs.apply_ucg(v, len(control_labels), single_qubit_gates) # update the diag according to the applied diagonal gate - _apply_diagonal_gate_to_diag( + diag = isometry_rs.apply_diagonal_gate_to_diag( diag, control_labels + [target_label], diagonal_ucg_inverse, n ) # # correct for the implementation "up to diagonal" # diag_inv = np.conj(diag).tolist() # _apply_diagonal_gate(v, control_labels + [target_label], diag_inv) + return v, diag # This method finds the single-qubit gates for a UCGate to disentangle a qubit: # we consider the n-qubit state v[:,0] starting with k zeros (in the computational basis). # The qubit with label n-s-1 is disentangled into the basis state k_s(k,s). def _find_squs_for_disentangling(self, v, k, s): - k_prime = 0 - n = int(math.log2(self.iso_data.shape[0])) - if _b(k, s + 1) == 0: - i_start = _a(k, s + 1) - else: - i_start = _a(k, s + 1) + 1 - id_list = [np.eye(2, 2) for _ in range(i_start)] - squs = [ - _reverse_qubit_state( - [ - v[2 * i * 2**s + _b(k, s), k_prime], - v[(2 * i + 1) * 2**s + _b(k, s), k_prime], - ], - _k_s(k, s), - self._epsilon, - ) - for i in range(i_start, 2 ** (n - s - 1)) - ] - return id_list + squs + res = isometry_rs.find_squs_for_disentangling( + v, k, s, self._epsilon, n=int(math.log2(self.iso_data.shape[0])) + ) + return res # Append a UCGate up to diagonal to the circuit circ. def _append_ucg_up_to_diagonal(self, circ, q, single_qubit_gates, control_labels, target_label): @@ -338,146 +343,6 @@ def inv_gate(self): return self._inverse -# Find special unitary matrix that maps [c0,c1] to [r,0] or [0,r] if basis_state=0 or -# basis_state=1 respectively -def _reverse_qubit_state(state, basis_state, epsilon): - state = np.array(state) - r = np.linalg.norm(state) - if r < epsilon: - return np.eye(2, 2) - if basis_state == 0: - m = np.array([[np.conj(state[0]), np.conj(state[1])], [-state[1], state[0]]]) / r - else: - m = np.array([[-state[1], state[0]], [np.conj(state[0]), np.conj(state[1])]]) / r - return m - - -# Methods for applying gates to matrices (should be moved to Qiskit AER) - -# Input: matrix m with 2^n rows (and arbitrary many columns). Think of the columns as states -# on n qubits. The method applies a uniformly controlled gate (UCGate) to all the columns, where -# the UCGate is specified by the inputs k and single_qubit_gates: - -# k = number of controls. We assume that the controls are on the k most significant qubits -# (and the target is on the (k+1)th significant qubit) -# single_qubit_gates = [u_0,...,u_{2^k-1}], where the u_i's are 2*2 unitaries -# (provided as numpy arrays) - -# The order of the single-qubit unitaries is such that the first unitary u_0 is applied to the -# (k+1)th significant qubit if the control qubits are in the state ket(0...00), the gate u_1 is -# applied if the control qubits are in the state ket(0...01), and so on. - -# The input matrix m and the single-qubit gates have to be of dtype=complex. - - -def _apply_ucg(m, k, single_qubit_gates): - # ToDo: Improve efficiency by parallelizing the gate application. A generalized version of - # ToDo: this method should be implemented by the state vector simulator in Qiskit AER. - num_qubits = int(math.log2(m.shape[0])) - num_col = m.shape[1] - spacing = 2 ** (num_qubits - k - 1) - for j in range(2 ** (num_qubits - 1)): - i = (j // spacing) * spacing + j - gate_index = i // (2 ** (num_qubits - k)) - for col in range(num_col): - m[np.array([i, i + spacing]), np.array([col, col])] = np.ndarray.flatten( - single_qubit_gates[gate_index].dot(np.array([[m[i, col]], [m[i + spacing, col]]])) - ).tolist() - return m - - -# Apply a diagonal gate with diagonal entries liste in diag and acting on qubits with labels -# action_qubit_labels to a matrix m. -# The input matrix m has to be of dtype=complex -# The qubit labels are such that label 0 corresponds to the most significant qubit, label 1 to -# the second most significant qubit, and so on ... - - -def _apply_diagonal_gate(m, action_qubit_labels, diag): - # ToDo: Improve efficiency by parallelizing the gate application. A generalized version of - # ToDo: this method should be implemented by the state vector simulator in Qiskit AER. - num_qubits = int(math.log2(m.shape[0])) - num_cols = m.shape[1] - basis_states = list(itertools.product([0, 1], repeat=num_qubits)) - for state in basis_states: - state_on_action_qubits = [state[i] for i in action_qubit_labels] - diag_index = _bin_to_int(state_on_action_qubits) - i = _bin_to_int(state) - for j in range(num_cols): - m[i, j] = diag[diag_index] * m[i, j] - return m - - -# Special case of the method _apply_diagonal_gate, where the input m is a diagonal matrix on the -# log2(len(m_diagonal)) least significant qubits (this method is more efficient in this case -# than _apply_diagonal_gate). The input m_diagonal is provided as a list of diagonal entries. -# The diagonal diag is applied on the qubits with labels listed in action_qubit_labels. The input -# num_qubits gives the total number of considered qubits (this input is required to interpret the -# action_qubit_labels in relation to the least significant qubits). - - -def _apply_diagonal_gate_to_diag(m_diagonal, action_qubit_labels, diag, num_qubits): - if not m_diagonal: - return m_diagonal - basis_states = list(itertools.product([0, 1], repeat=num_qubits)) - for state in basis_states[: len(m_diagonal)]: - state_on_action_qubits = [state[i] for i in action_qubit_labels] - diag_index = _bin_to_int(state_on_action_qubits) - i = _bin_to_int(state) - m_diagonal[i] *= diag[diag_index] - return m_diagonal - - -# Apply a MC single-qubit gate (given by the 2*2 unitary input: gate) with controlling on -# the qubits with label control_labels and acting on the qubit with label target_label -# to a matrix m. The input matrix m and the gate have to be of dtype=complex. The qubit labels are -# such that label 0 corresponds to the most significant qubit, label 1 to the second most -# significant qubit, and so on ... - - -def _apply_multi_controlled_gate(m, control_labels, target_label, gate): - # ToDo: This method should be integrated into the state vector simulator in Qiskit AER. - num_qubits = int(math.log2(m.shape[0])) - num_cols = m.shape[1] - control_labels.sort() - free_qubits = num_qubits - len(control_labels) - 1 - basis_states_free = list(itertools.product([0, 1], repeat=free_qubits)) - for state_free in basis_states_free: - (e1, e2) = _construct_basis_states(state_free, control_labels, target_label) - for i in range(num_cols): - m[np.array([e1, e2]), np.array([i, i])] = np.ndarray.flatten( - gate.dot(np.array([[m[e1, i]], [m[e2, i]]])) - ).tolist() - return m - - -# Helper method for _apply_multi_controlled_gate. This constructs the basis states the MG gate -# is acting on for a specific state state_free of the qubits we neither control nor act on. - - -def _construct_basis_states(state_free, control_labels, target_label): - e1 = [] - e2 = [] - j = 0 - for i in range(len(state_free) + len(control_labels) + 1): - if i in control_labels: - e1.append(1) - e2.append(1) - elif i == target_label: - e1.append(0) - e2.append(1) - else: - e1.append(state_free[j]) - e2.append(state_free[j]) - j += 1 - out1 = _bin_to_int(e1) - out2 = _bin_to_int(e2) - return out1, out2 - - -# Some helper methods: - - # Get the qubits in the list qubits corresponding to the labels listed in labels. The total number # of qubits is given by num_qubits (and determines the convention for the qubit labeling) @@ -496,14 +361,6 @@ def _reverse_qubit_oder(qubits): # Convert list of binary digits to integer -def _bin_to_int(binary_digits_as_list): - return int("".join(str(x) for x in binary_digits_as_list), 2) - - -def _ct(m): - return np.transpose(np.conjugate(m)) - - def _get_binary_rep_as_list(n, num_digits): binary_string = np.binary_repr(n).zfill(num_digits) binary = [] @@ -511,64 +368,3 @@ def _get_binary_rep_as_list(n, num_digits): for c in line: binary.append(int(c)) return binary[-num_digits:] - - -# absorb a diagonal gate into a UCGate - - -def _merge_UCGate_and_diag(single_qubit_gates, diag): - for i, gate in enumerate(single_qubit_gates): - single_qubit_gates[i] = np.array([[diag[2 * i], 0.0], [0.0, diag[2 * i + 1]]]).dot(gate) - return single_qubit_gates - - -# Helper variables/functions for the column-by-column decomposition - - -# a(k,s) and b(k,s) are positive integers such that k = a(k,s)2^s + b(k,s) -# (with the maximal choice of a(k,s)) - - -def _a(k, s): - return k // 2**s - - -def _b(k, s): - return k - (_a(k, s) * 2**s) - - -# given a binary representation of k with binary digits [k_{n-1},..,k_1,k_0], -# the method k_s(k, s) returns k_s - - -def _k_s(k, s): - if k == 0: - return 0 - else: - num_digits = s + 1 - return _get_binary_rep_as_list(k, num_digits)[0] - - -# Check if a gate of a special form is equal to the identity gate up to global phase - - -def _ucg_is_identity_up_to_global_phase(single_qubit_gates, epsilon): - if not np.abs(single_qubit_gates[0][0, 0]) < epsilon: - global_phase = 1.0 / (single_qubit_gates[0][0, 0]) - else: - return False - for gate in single_qubit_gates: - if not np.allclose(global_phase * gate, np.eye(2, 2)): - return False - return True - - -def _diag_is_identity_up_to_global_phase(diag, epsilon): - if not np.abs(diag[0]) < epsilon: - global_phase = 1.0 / (diag[0]) - else: - return False - for d in diag: - if not np.abs(global_phase * d - 1) < epsilon: - return False - return True diff --git a/qiskit/circuit/library/generalized_gates/uc.py b/qiskit/circuit/library/generalized_gates/uc.py index 2d650e98466a..f54567123e02 100644 --- a/qiskit/circuit/library/generalized_gates/uc.py +++ b/qiskit/circuit/library/generalized_gates/uc.py @@ -21,7 +21,6 @@ from __future__ import annotations -import cmath import math import numpy as np @@ -33,14 +32,11 @@ from qiskit.circuit.quantumcircuit import QuantumCircuit from qiskit.circuit.exceptions import CircuitError from qiskit.exceptions import QiskitError - -# pylint: disable=cyclic-import -from qiskit.synthesis.one_qubit.one_qubit_decompose import OneQubitEulerDecomposer +from qiskit._accelerate import uc_gate from .diagonal import Diagonal _EPS = 1e-10 # global variable used to chop very small numbers to zero -_DECOMPOSER1Q = OneQubitEulerDecomposer("U3") class UCGate(Gate): @@ -203,99 +199,7 @@ def _dec_ucg_help(self): https://arxiv.org/pdf/quant-ph/0410066.pdf. """ single_qubit_gates = [gate.astype(complex) for gate in self.params] - diag = np.ones(2**self.num_qubits, dtype=complex) - num_contr = self.num_qubits - 1 - for dec_step in range(num_contr): - num_ucgs = 2**dec_step - # The decomposition works recursively and the following loop goes over the different - # UCGates that arise in the decomposition - for ucg_index in range(num_ucgs): - len_ucg = 2 ** (num_contr - dec_step) - for i in range(int(len_ucg / 2)): - shift = ucg_index * len_ucg - a = single_qubit_gates[shift + i] - b = single_qubit_gates[shift + len_ucg // 2 + i] - # Apply the decomposition for UCGates given in equation (3) in - # https://arxiv.org/pdf/quant-ph/0410066.pdf - # to demultiplex one control of all the num_ucgs uniformly-controlled gates - # with log2(len_ucg) uniform controls - v, u, r = self._demultiplex_single_uc(a, b) - # replace the single-qubit gates with v,u (the already existing ones - # are not needed any more) - single_qubit_gates[shift + i] = v - single_qubit_gates[shift + len_ucg // 2 + i] = u - # Now we decompose the gates D as described in Figure 4 in - # https://arxiv.org/pdf/quant-ph/0410066.pdf and merge some of the gates - # into the UCGates and the diagonal at the end of the circuit - - # Remark: The Rz(pi/2) rotation acting on the target qubit and the Hadamard - # gates arising in the decomposition of D are ignored for the moment (they will - # be added together with the C-NOT gates at the end of the decomposition - # (in the method dec_ucg())) - if ucg_index < num_ucgs - 1: - # Absorb the Rz(pi/2) rotation on the control into the UC-Rz gate and - # merge the UC-Rz rotation with the following UCGate, - # which hasn't been decomposed yet. - k = shift + len_ucg + i - single_qubit_gates[k] = single_qubit_gates[k].dot( - UCGate._ct(r) - ) * UCGate._rz(np.pi / 2).item((0, 0)) - k = k + len_ucg // 2 - single_qubit_gates[k] = single_qubit_gates[k].dot(r) * UCGate._rz( - np.pi / 2 - ).item((1, 1)) - else: - # Absorb the Rz(pi/2) rotation on the control into the UC-Rz gate and merge - # the trailing UC-Rz rotation into a diagonal gate at the end of the circuit - for ucg_index_2 in range(num_ucgs): - shift_2 = ucg_index_2 * len_ucg - k = 2 * (i + shift_2) - diag[k] = ( - diag[k] - * UCGate._ct(r).item((0, 0)) - * UCGate._rz(np.pi / 2).item((0, 0)) - ) - diag[k + 1] = ( - diag[k + 1] - * UCGate._ct(r).item((1, 1)) - * UCGate._rz(np.pi / 2).item((0, 0)) - ) - k = len_ucg + k - diag[k] *= r.item((0, 0)) * UCGate._rz(np.pi / 2).item((1, 1)) - diag[k + 1] *= r.item((1, 1)) * UCGate._rz(np.pi / 2).item((1, 1)) - return single_qubit_gates, diag - - def _demultiplex_single_uc(self, a, b): - """ - This method implements the decomposition given in equation (3) in - https://arxiv.org/pdf/quant-ph/0410066.pdf. - The decomposition is used recursively to decompose uniformly controlled gates. - a,b = single qubit unitaries - v,u,r = outcome of the decomposition given in the reference mentioned above - (see there for the details). - """ - # The notation is chosen as in https://arxiv.org/pdf/quant-ph/0410066.pdf. - x = a.dot(UCGate._ct(b)) - det_x = np.linalg.det(x) - x11 = x.item((0, 0)) / cmath.sqrt(det_x) - phi = cmath.phase(det_x) - r1 = cmath.exp(1j / 2 * (np.pi / 2 - phi / 2 - cmath.phase(x11))) - r2 = cmath.exp(1j / 2 * (np.pi / 2 - phi / 2 + cmath.phase(x11) + np.pi)) - r = np.array([[r1, 0], [0, r2]], dtype=complex) - d, u = np.linalg.eig(r.dot(x).dot(r)) - # If d is not equal to diag(i,-i), then we put it into this "standard" form - # (see eq. (13) in https://arxiv.org/pdf/quant-ph/0410066.pdf) by interchanging - # the eigenvalues and eigenvectors. - if abs(d[0] + 1j) < _EPS: - d = np.flip(d, 0) - u = np.flip(u, 1) - d = np.diag(np.sqrt(d)) - v = d.dot(UCGate._ct(u)).dot(UCGate._ct(r)).dot(b) - return v, u, r - - @staticmethod - def _ct(m): - return np.transpose(np.conjugate(m)) + return uc_gate.dec_ucg_help(single_qubit_gates, self.num_qubits) @staticmethod def _rz(alpha):