diff --git a/qiskit/circuit/library/arithmetic/piecewise_chebyshev.py b/qiskit/circuit/library/arithmetic/piecewise_chebyshev.py new file mode 100644 index 000000000000..65cdc4a444a5 --- /dev/null +++ b/qiskit/circuit/library/arithmetic/piecewise_chebyshev.py @@ -0,0 +1,327 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2017, 2020. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""Piecewise polynomial Chebyshev approximation to a given f(x).""" + +from typing import Callable, List, Optional +import numpy as np +from numpy.polynomial.chebyshev import Chebyshev + +from qiskit.circuit import QuantumRegister, AncillaRegister +from qiskit.circuit.library.blueprintcircuit import BlueprintCircuit +from qiskit.circuit.exceptions import CircuitError + +from .piecewise_polynomial_pauli_rotations import PiecewisePolynomialPauliRotations + + +class PiecewiseChebyshev(BlueprintCircuit): + r"""Piecewise Chebyshev approximation to an input function. + + For a given function :math:`f(x)` and degree :math:`d`, this class implements a piecewise + polynomial Chebyshev approximation on :math:`n` qubits to :math:`f(x)` on the given intervals. + All the polynomials in the approximation are of degree :math:`d`. + + The values of the parameters are calculated according to [1]. + + Examples: + + .. jupyer-execute:: + + import numpy as np + from qiskit import QuantumCircuit + from qiskit.circuit.library.arithmetic.piecewise_chebyshev import PiecewiseChebyshev + f_x, degree, breakpoints, num_state_qubits = lambda x: np.arcsin(1 / x), 2, [2, 4], 2 + pw_approximation = PiecewiseChebyshev(f_x, degree, breakpoints, num_state_qubits) + pw_approximation._build() + qc = QuantumCircuit(pw_approximation.num_qubits) + qc.h(list(range(num_state_qubits))) + qc.append(pw_approximation.to_instruction(), qc.qubits) + print(qc.draw(output='mpl')) + + References: + + [1]: Haener, T., Roetteler, M., & Svore, K. M. (2018). + Optimizing Quantum Circuits for Arithmetic. + `arXiv:1805.12445 `_ + """ + + def __init__(self, + f_x: Callable[[int], float], + degree: Optional[int] = None, + breakpoints: Optional[List[int]] = None, + num_state_qubits: Optional[int] = None, + name: str = 'pw_cheb') -> None: + r""" + Args: + f_x: the function to be approximated. + degree: the degree of the polynomials. + Defaults to ``1``. + breakpoints: the breakpoints to define the piecewise-linear function. + Defaults to the full interval. + num_state_qubits: number of qubits representing the state. + name: The name of the circuit object. + """ + super().__init__(name=name) + + # define internal parameters + self._num_state_qubits = None + + # Store parameters + self._f_x = f_x + self._degree = degree if degree is not None else 1 + self._breakpoints = breakpoints if breakpoints is not None else [0] + + self._polynomials = None + + self.num_state_qubits = num_state_qubits + + def _check_configuration(self, raise_on_failure: bool = True) -> bool: + valid = True + + if self._f_x is None: + valid = False + if raise_on_failure: + raise AttributeError('The function to be approximated has not been set.') + + if self._degree is None: + valid = False + if raise_on_failure: + raise AttributeError('The degree of the polynomials has not been set.') + + if self._breakpoints is None: + valid = False + if raise_on_failure: + raise AttributeError('The breakpoints have not been set.') + + if self.num_state_qubits is None: + valid = False + if raise_on_failure: + raise AttributeError('The number of qubits has not been set.') + + if self.num_qubits < self.num_state_qubits + 1: + valid = False + if raise_on_failure: + raise CircuitError('Not enough qubits in the circuit, need at least ' + '{}.'.format(self.num_state_qubits + 1)) + + return valid + + @property + def f_x(self) -> Callable[[int], float]: + """The function to be approximated. + + Returns: + The function to be approximated. + """ + return self._f_x + + @f_x.setter + def f_x(self, f_x: Optional[Callable[[int], float]]) -> None: + """Set the function to be approximated. + + Note that this may change the underlying quantum register, if the number of state qubits + changes. + + Args: + f_x: The new function to be approximated. + """ + if self._f_x is None or f_x != self._f_x: + self._invalidate() + self._f_x = f_x + + self._reset_registers(self.num_state_qubits) + + @property + def degree(self) -> int: + """The degree of the polynomials. + + Returns: + The degree of the polynomials. + """ + return self._degree + + @degree.setter + def degree(self, degree: Optional[int]) -> None: + """Set the error tolerance. + + Note that this may change the underlying quantum register, if the number of state qubits + changes. + + Args: + degree: The new degree. + """ + if self._degree is None or degree != self._degree: + self._invalidate() + self._degree = degree + + self._reset_registers(self.num_state_qubits) + + @property + def breakpoints(self) -> List[int]: + """The breakpoints for the piecewise approximation. + + Returns: + The breakpoints for the piecewise approximation. + """ + breakpoints = self._breakpoints + + # it the state qubits are set ensure that the breakpoints match beginning and end + if self.num_state_qubits is not None: + num_states = 2 ** self.num_state_qubits + + # If the last breakpoint is < num_states, add the identity polynomial + if breakpoints[-1] < num_states: + breakpoints = breakpoints + [num_states] + + # If the first breakpoint is > 0, add the identity polynomial + if breakpoints[0] > 0: + breakpoints = [0] + breakpoints + + return breakpoints + + @breakpoints.setter + def breakpoints(self, breakpoints: Optional[List[int]]) -> None: + """Set the breakpoints for the piecewise approximation. + + Note that this may change the underlying quantum register, if the number of state qubits + changes. + + Args: + breakpoints: The new breakpoints for the piecewise approximation. + """ + if self._breakpoints is None or breakpoints != self._breakpoints: + self._invalidate() + self._breakpoints = breakpoints + + self._reset_registers(self.num_state_qubits) + + @property + def polynomials(self) -> List[List[float]]: + """The polynomials for the piecewise approximation. + + Returns: + The polynomials for the piecewise approximation. + """ + if self.num_state_qubits is None: + return [[]] + + # note this must be the private attribute since we handle missing breakpoints at + # 0 and 2 ^ num_qubits here (e.g. if the function we approximate is not defined at 0 + # and the user takes that into account we just add an identity) + num_intervals = len(self._breakpoints) + + # Calculate the polynomials + polynomials = [] + for i in range(0, num_intervals - 1): + # Calculate the polynomial approximating the function on the current interval + poly = Chebyshev.interpolate(self._f_x, self._degree, + domain=[self._breakpoints[i], + self._breakpoints[i + 1]]) + # Convert polynomial to the standard basis and rescale it for the rotation gates + poly = 2 * poly.convert(kind=np.polynomial.Polynomial).coef + # Convert to list and append + polynomials.append(poly.tolist()) + + # If the last breakpoint is < 2 ** num_qubits, add the identity polynomial + if self._breakpoints[-1] < 2 ** self.num_state_qubits: + polynomials = polynomials + [[2 * np.arcsin(1)]] + + # If the first breakpoint is > 0, add the identity polynomial + if self._breakpoints[0] > 0: + polynomials = [[2 * np.arcsin(1)]] + polynomials + + return polynomials + + @polynomials.setter + def polynomials(self, polynomials: Optional[List[List[float]]]) -> None: + """Set the polynomials for the piecewise approximation. + + Note that this may change the underlying quantum register, if the number of state qubits + changes. + + Args: + polynomials: The new breakpoints for the piecewise approximation. + """ + if self._polynomials is None or polynomials != self._polynomials: + self._invalidate() + self._polynomials = polynomials + + self._reset_registers(self.num_state_qubits) + + @property + def num_state_qubits(self) -> int: + r"""The number of state qubits representing the state :math:`|x\rangle`. + + Returns: + The number of state qubits. + """ + return self._num_state_qubits + + @num_state_qubits.setter + def num_state_qubits(self, num_state_qubits: Optional[int]) -> None: + """Set the number of state qubits. + + Note that this may change the underlying quantum register, if the number of state qubits + changes. + + Args: + num_state_qubits: The new number of qubits. + """ + if self._num_state_qubits is None or num_state_qubits != self._num_state_qubits: + self._invalidate() + self._num_state_qubits = num_state_qubits + + # Set breakpoints if they haven't been set + if num_state_qubits is not None and self._breakpoints is None: + self.breakpoints = [0, 2 ** num_state_qubits] + + self._reset_registers(num_state_qubits) + + def _reset_registers(self, num_state_qubits: Optional[int]) -> None: + if num_state_qubits is not None: + qr_state = QuantumRegister(num_state_qubits, 'state') + qr_target = QuantumRegister(1, 'target') + self.qregs = [qr_state, qr_target] + self._ancillas = [] + self._qubits = qr_state[:] + qr_target[:] + + num_ancillas = num_state_qubits + if num_ancillas > 0: + qr_ancilla = AncillaRegister(num_ancillas) + self.add_register(qr_ancilla) + + else: + self.qregs = [] + self._qubits = [] + self._ancillas = [] + + def _build(self): + """Build the circuit. The operation is considered successful when q_objective is :math:`|1>` + """ + # do not build the circuit if _data is already populated + if self._data is not None: + return + + self._data = [] + + # check whether the configuration is valid + self._check_configuration() + + poly_r = PiecewisePolynomialPauliRotations(self.num_state_qubits, + self.breakpoints, self.polynomials) + + qr_state = self.qubits[:self.num_state_qubits] + qr_target = [self.qubits[self.num_state_qubits]] + qr_ancillas = self.qubits[self.num_state_qubits + 1:] + + # Apply polynomial approximation + self.append(poly_r.to_instruction(), qr_state[:] + qr_target + qr_ancillas[:]) diff --git a/qiskit/circuit/library/arithmetic/piecewise_polynomial_pauli_rotations.py b/qiskit/circuit/library/arithmetic/piecewise_polynomial_pauli_rotations.py index bfe8726b6a5b..ff576c4da80a 100644 --- a/qiskit/circuit/library/arithmetic/piecewise_polynomial_pauli_rotations.py +++ b/qiskit/circuit/library/arithmetic/piecewise_polynomial_pauli_rotations.py @@ -127,6 +127,10 @@ def breakpoints(self) -> List[int]: Returns: The list of breakpoints. """ + if self.num_state_qubits is not None and len(self._breakpoints) == len(self.coeffs) and\ + self._breakpoints[-1] < 2 ** self.num_state_qubits: + return self._breakpoints + [2 ** self.num_state_qubits] + return self._breakpoints @breakpoints.setter @@ -246,20 +250,27 @@ def _reset_registers(self, num_state_qubits: Optional[int]) -> None: if self.contains_zero_breakpoint: num_ancillas -= 1 if num_ancillas > 0: - self._ancillas = [] qr_ancilla = AncillaRegister(num_ancillas) self.add_register(qr_ancilla) + else: + qr_ancilla = [] + + self._qubits = qr_state[:] + qr_target[:] + qr_ancilla[:] + self._ancillas = qr_ancilla[:] else: self.qregs = [] + self._qubits = [] + self._ancillas = [] def _build(self): - # Add the last breakpoint if necessary - if self.num_state_qubits is not None and len(self._breakpoints) == len(self._coeffs) and\ - self._breakpoints[-1] < 2 ** self.num_state_qubits: - self.breakpoints = self._breakpoints + [2 ** self.num_state_qubits] + # do not build the circuit if _data is already populated + if self._data is not None: + return + + self._data = [] - # The number of ancilla might have changed, so reset registers - super()._build() + # check whether the configuration is valid + self._check_configuration() qr_state = self.qubits[:self.num_state_qubits] qr_target = [self.qubits[self.num_state_qubits]] diff --git a/qiskit/circuit/library/arithmetic/polynomial_pauli_rotations.py b/qiskit/circuit/library/arithmetic/polynomial_pauli_rotations.py index c43767a68b66..3c093a40e25e 100644 --- a/qiskit/circuit/library/arithmetic/polynomial_pauli_rotations.py +++ b/qiskit/circuit/library/arithmetic/polynomial_pauli_rotations.py @@ -303,7 +303,14 @@ def _get_rotation_coefficients(self) -> Dict[Sequence[int], float]: return rotation_coeffs def _build(self): - super()._build() + # do not build the circuit if _data is already populated + if self._data is not None: + return + + self._data = [] + + # check whether the configuration is valid + self._check_configuration() qr_state = self.qubits[:self.num_state_qubits] qr_target = self.qubits[self.num_state_qubits] diff --git a/qiskit/circuit/library/blueprintcircuit.py b/qiskit/circuit/library/blueprintcircuit.py index bf26cbcff327..2131e477b6e1 100644 --- a/qiskit/circuit/library/blueprintcircuit.py +++ b/qiskit/circuit/library/blueprintcircuit.py @@ -59,7 +59,7 @@ def _check_configuration(self, raise_on_failure: bool = True) -> bool: def _build(self) -> None: """Build the circuit.""" # do not build the circuit if _data is already populated - if self._data: + if self._data is not None: return self._data = [] diff --git a/releasenotes/notes/add-piecewise-polynomial-chebyshev-66c82ffe1d749dd6.yaml b/releasenotes/notes/add-piecewise-polynomial-chebyshev-66c82ffe1d749dd6.yaml new file mode 100644 index 000000000000..3bb21dc3f530 --- /dev/null +++ b/releasenotes/notes/add-piecewise-polynomial-chebyshev-66c82ffe1d749dd6.yaml @@ -0,0 +1,6 @@ +--- +features: + - | + For a given function :math:`f(x)` and degree :math:`d`, this class implements a piecewise + polynomial Chebyshev approximation on :math:`n` qubits to :math:`f(x)` on the given intervals. + All the polynomials in the approximation are of degree :math:`d`. diff --git a/test/python/circuit/library/test_piecewise_chebyshev.py b/test/python/circuit/library/test_piecewise_chebyshev.py new file mode 100644 index 000000000000..78940f91b7b0 --- /dev/null +++ b/test/python/circuit/library/test_piecewise_chebyshev.py @@ -0,0 +1,127 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2017, 2020. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""Test the piecewise Chebyshev approximation.""" + +import unittest +from collections import defaultdict +import numpy as np +from ddt import ddt, data, unpack + +from qiskit.test.base import QiskitTestCase +from qiskit import BasicAer, execute +from qiskit.circuit import QuantumCircuit +from qiskit.circuit.library.arithmetic.piecewise_chebyshev import PiecewiseChebyshev + + +@ddt +class TestPiecewiseChebyshev(QiskitTestCase): + """Test the piecewise Chebyshev approximation.""" + + def assertFunctionIsCorrect(self, function_circuit, reference): + """Assert that ``function_circuit`` implements the reference function ``reference``.""" + function_circuit._build() + num_state_qubits = function_circuit.num_state_qubits + num_ancilla_qubits = function_circuit.num_ancillas + circuit = QuantumCircuit(num_state_qubits + 1 + num_ancilla_qubits) + circuit.h(list(range(num_state_qubits))) + circuit.append(function_circuit.to_instruction(), list(range(circuit.num_qubits))) + + backend = BasicAer.get_backend('statevector_simulator') + statevector = execute(circuit, backend).result().get_statevector() + + probabilities = defaultdict(float) + for i, statevector_amplitude in enumerate(statevector): + i = bin(i)[2:].zfill(circuit.num_qubits)[num_ancilla_qubits:] + probabilities[i] += np.real(np.abs(statevector_amplitude) ** 2) + + unrolled_probabilities = [] + unrolled_expectations = [] + for i, probability in probabilities.items(): + x, last_qubit = int(i[1:], 2), i[0] + if last_qubit == '0': + expected_amplitude = np.cos(reference(x)) / np.sqrt(2 ** num_state_qubits) + else: + expected_amplitude = np.sin(reference(x)) / np.sqrt(2 ** num_state_qubits) + + unrolled_probabilities += [probability] + unrolled_expectations += [np.real(np.abs(expected_amplitude) ** 2)] + + np.testing.assert_array_almost_equal(unrolled_probabilities, unrolled_expectations, + decimal=3) + + @ data((lambda x: np.arcsin(1 / x), 2, [2, 4], 2), + (lambda x: x / 8, 1, [1, 8], 3) + ) + @ unpack + def test_piecewise_chebyshev(self, f_x, degree, breakpoints, num_state_qubits): + """Test the piecewise Chebyshev approximation.""" + + def pw_poly(x): + if breakpoints[0] <= x < breakpoints[-1]: + return f_x(x) + return np.arcsin(1) + + pw_approximation = PiecewiseChebyshev(f_x, degree, breakpoints, num_state_qubits) + + self.assertFunctionIsCorrect(pw_approximation, pw_poly) + + def test_piecewise_chebyshev_mutability(self): + """Test the mutability of the piecewise Chebyshev approximation.""" + + def pw_poly(x, f_x): + if breakpoints[0] <= x < breakpoints[-1]: + return f_x(x) + return np.arcsin(1) + + def f_x_1(x): + return x / 2 + + pw_approximation = PiecewiseChebyshev(f_x_1) + + with self.subTest(msg='missing number of state qubits'): + with self.assertRaises(AttributeError): # no state qubits set + print(pw_approximation.draw()) + + with self.subTest(msg='default setup, just setting number of state qubits'): + pw_approximation.num_state_qubits = 2 + pw_approximation.f_x = f_x_1 + # set to the default breakpoints for pw_poly + breakpoints = [0, 4] + pw_approximation.breakpoints = breakpoints + self.assertFunctionIsCorrect(pw_approximation, lambda x: pw_poly(x, f_x_1)) + + def f_x_2(x): + return x / 4 + with self.subTest(msg='setting non-default values'): + breakpoints = [0, 2] + degree = 2 + pw_approximation.breakpoints = breakpoints + pw_approximation.degree = degree + pw_approximation.f_x = f_x_2 + self.assertFunctionIsCorrect(pw_approximation, lambda x: pw_poly(x, f_x_2)) + + def f_x_3(x): + return x ** 2 + + with self.subTest(msg='changing all values'): + pw_approximation.num_state_qubits = 4 + breakpoints = [1, 3, 6] + degree = 3 + pw_approximation.breakpoints = breakpoints + pw_approximation.degree = degree + pw_approximation.f_x = f_x_3 + self.assertFunctionIsCorrect(pw_approximation, lambda x: pw_poly(x, f_x_3)) + + +if __name__ == '__main__': + unittest.main()