Source code for pennylane.templates.subroutines.qsvt
# Copyright 2018-2023 Xanadu Quantum Technologies Inc.
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
"""
Contains the QSVT template and qsvt wrapper function.
"""
# pylint: disable=too-many-arguments
import copy
from typing import Literal
import numpy as np
from numpy.polynomial import Polynomial, chebyshev
import pennylane as qml
from pennylane.operation import AnyWires, Operation
from pennylane.ops.op_math import adjoint
from pennylane.queuing import QueuingManager
from pennylane.wires import Wires
# pylint: disable=too-many-branches, unused-argument
[docs]def qsvt(A, poly, encoding_wires=None, block_encoding=None, **kwargs):
r"""
Implements the Quantum Singular Value Transformation (QSVT) for a matrix or Hamiltonian ``A``,
using a polynomial defined by ``poly`` and a block encoding specified by ``block_encoding``.
.. math::
\begin{pmatrix}
A & * \\
* & * \\
\end{pmatrix}
\xrightarrow{QSVT}
\begin{pmatrix}
\text{poly}(A) + i \dots & * \\
* & * \\
\end{pmatrix}
The polynomial transformation is encoded as the real part of the top left term after applying the operator.
This function calculates the required phase angles from the polynomial using :func:`~.poly_to_angles`.
.. note::
The function ``poly_to_angles``, used within ``qsvt``, is not JIT-compatible, which
prevents ``poly`` from being traceable in ``qsvt``. However, ``A`` is traceable
and can be optimized by JIT within this function.
Args:
A (Union[tensor_like, Operator]): The matrix on which the QSVT will be applied.
This can be an array or an object that has a Pauli representation. See :func:`~.pauli_decompose`.
poly (tensor_like): coefficients of the polynomial ordered from lowest to highest power
encoding_wires (Sequence[int]): The qubit wires used for the block encoding. See Usage Details below for
more information on ``encoding_wires`` depending on the block encoding used.
block_encoding (str): Specifies the type of block encoding to use. Options include:
- ``"prepselprep"``: Embeds the Hamiltonian ``A`` using :class:`~pennylane.PrepSelPrep`.
Default encoding for Hamiltonians.
- ``"qubitization"``: Embeds the Hamiltonian ``A`` using :class:`~pennylane.Qubitization`.
- ``"embedding"``: Embeds the matrix ``A`` using :class:`~pennylane.BlockEncode`.
Template not hardware compatible. Default encoding for matrices.
- ``"fable"``: Embeds the matrix ``A`` using :class:`~pennylane.FABLE`. Template hardware compatible.
Returns:
(Operator): A quantum operator implementing QSVT on the matrix ``A`` with the
specified encoding and projector phases.
.. seealso:: :class:`~.QSVT`
Example:
.. code-block:: python
# P(x) = -x + 0.5 x^3 + 0.5 x^5
poly = np.array([0, -1, 0, 0.5, 0, 0.5])
hamiltonian = qml.dot([0.3, 0.7], [qml.Z(1), qml.X(1) @ qml.Z(2)])
dev = qml.device("default.qubit")
@qml.qnode(dev)
def circuit():
qml.qsvt(hamiltonian, poly, encoding_wires=[0], block_encoding="prepselprep")
return qml.state()
matrix = qml.matrix(circuit, wire_order=[0, 1, 2])()
.. code-block:: pycon
>>> print(matrix[:4, :4].real)
[[-0.1625 0. -0.3793 0. ]
[ 0. -0.1625 0. 0.3793]
[-0.3793 0. 0.1625 0. ]
[ 0. 0.3793 0. 0.1625]]
.. details::
:title: Usage Details
If the input ``A`` is a Hamiltonian, the valid ``block_encoding`` values are
``"prepselprep"`` and ``"qubitization"``. In this case, ``encoding_wires`` refers to the
``control`` parameter in the templates :class:`~pennylane.PrepSelPrep` and :class:`~pennylane.Qubitization`,
respectively. These wires represent the auxiliary qubits necessary for the block encoding of
the Hamiltonian. The number of ``encoding_wires`` required must be :math:`\lceil \log_2(m) \rceil`,
where :math:`m` is the number of terms in the Hamiltonian.
.. code-block:: python
# P(x) = -1 + 0.2 x^2 + 0.5 x^4
poly = np.array([-1, 0, 0.2, 0, 0.5])
hamiltonian = qml.dot([0.3, 0.4, 0.3], [qml.Z(2), qml.X(2) @ qml.Z(3), qml.X(2)])
dev = qml.device("default.qubit")
@qml.qnode(dev)
def circuit():
qml.qsvt(hamiltonian, poly, encoding_wires=[0, 1], block_encoding="prepselprep")
return qml.state()
matrix = qml.matrix(circuit, wire_order=[0, 1, 2, 3])()
.. code-block:: pycon
>>> print(np.round(matrix[:4, :4], 4).real)
[[-0.7158 0. 0. 0. ]
[ 0. -0.975 0. 0. ]
[ 0. 0. -0.7158 0. ]
[ 0. -0. 0. -0.975 ]]
Alternatively, if the input ``A`` is a matrix, the valid values for ``block_encoding`` are
``"embedding"`` and ``"fable"``. In this case, the ``encoding_wires`` parameter corresponds to
the ``wires`` attribute in the templates :class:`~pennylane.BlockEncode` and :class:`~pennylane.FABLE`,
respectively. Note that for QSVT to work, the input matrix must be Hermitian.
.. code-block:: python
# P(x) = -1 + 0.2 x^2 + 0.5 x^4
poly = np.array([-0.1, 0, 0.2, 0, 0.5])
A = np.array([[-0.1, 0, 0, 0.1], [0, 0.2, 0, 0], [0, 0, -0.2, -0.2], [0.1, 0, -0.2, -0.1]])
dev = qml.device("default.qubit")
@qml.qnode(dev)
def circuit():
qml.qsvt(A, poly, encoding_wires=[0, 1, 2, 3, 4], block_encoding="fable")
return qml.state()
matrix = qml.matrix(circuit, wire_order=[0, 1, 2, 3, 4])()
.. code-block:: pycon
>>> print(np.round(matrix[:4, :4], 4).real)
[[-0.0954 0. -0.0056 -0.0054]
[ 0. -0.0912 -0. -0. ]
[-0.0056 0. -0.0788 0.0164]
[-0.0054 -0. 0.0164 -0.0842]]
Note that for the FABLE block encoding to function correctly, it must comply with the following:
.. math::
d \|A\|^2 \leq 1,
where :math:`d` is the maximum dimension of :math:`A` and :math:`\|A\|` is the 2-norm of :math:`A`.
In the previous example this is satisfied since :math:`d = 4` and :math:`\|A\|^2 = 0.2`:
.. code-block:: pycon
>>> print(4* np.linalg.norm(A, ord='fro')**2)
0.8000000000000004
"""
projectors = []
# If the input A is a Hamiltonian
if hasattr(A, "pauli_rep"):
if block_encoding not in ["prepselprep", "qubitization", None]:
raise ValueError(
"block_encoding = {block_encoding} not supported for A of type {type(A)}. When A is a Hamiltonian or has a Pauli decomposition, block_encoding should take the value 'prepselprep' or 'qubitization'. Otherwise, please provide the matrix of the Hamiltonian as input. For more details, see the 'qml.matrix' function."
)
if any(wire in qml.wires.Wires(encoding_wires) for wire in A.wires):
raise ValueError(
f"Control wires in '{block_encoding}' should be different from the hamiltonian wires"
)
encoding = (
qml.Qubitization(A, control=encoding_wires)
if block_encoding == "qubitization"
else qml.PrepSelPrep(A, control=encoding_wires)
)
angles = qml.poly_to_angles(poly, "QSVT")
projectors = [
qml.PCPhase(angle, dim=2 ** len(A.wires), wires=encoding_wires + A.wires)
for angle in angles
]
else:
if block_encoding not in ["embedding", "fable", None]:
raise ValueError(
"block_encoding = {block_encoding} not supported for A of type {type(A)}. When A is a matrix block_encoding should take the value 'embedding' or 'fable'. Otherwise, please provide an input with a Pauli decomposition. For more details, see the 'qml.pauli_decompose' function."
)
A = qml.math.atleast_2d(A)
max_dimension = 1 if len(qml.math.array(A).shape) == 0 else max(A.shape)
if block_encoding == "fable":
if qml.math.linalg.norm(max_dimension * qml.math.ravel(A), np.inf) > 1:
raise ValueError(
"The subnormalization factor should be lower than 1. Ensure that the product of the maximum dimension of A and its square norm is less than 1."
)
# FABLE encodes A / 2^n, need to rescale to obtain desired block-encoding
fable_norm = int(np.ceil(np.log2(max_dimension)))
encoding = qml.FABLE(2**fable_norm * A, wires=encoding_wires)
angles = qml.poly_to_angles(poly, "QSVT")
projectors = [qml.PCPhase(angle, dim=len(A), wires=encoding_wires) for angle in angles]
else:
c, r = qml.math.shape(A)
angles = qml.poly_to_angles(poly, "QSVT")
for idx, phi in enumerate(angles):
dim = c if idx % 2 else r
projectors.append(qml.PCPhase(phi, dim=dim, wires=encoding_wires))
encoding = qml.BlockEncode(A, wires=encoding_wires)
return QSVT(encoding, projectors)
[docs]class QSVT(Operation):
r"""QSVT(UA,projectors)
Implements the
`quantum singular value transformation <https://arxiv.org/abs/1806.01838>`__ (QSVT) circuit.
.. note ::
This template allows users to define hardware-compatible block encoding and
projector-controlled phase shift circuits. For a QSVT implementation that is
tailored to work directly with an input matrix and a transformation polynomial
see :func:`~.qsvt`.
Given an :class:`~.Operator` :math:`U`, which block encodes the matrix :math:`A`, and a list of
projector-controlled phase shift operations :math:`\vec{\Pi}_\phi`, this template applies a
circuit for the quantum singular value transformation as follows.
When the number of projector-controlled phase shifts is even (:math:`d` is odd), the QSVT
circuit is defined as:
.. math::
U_{QSVT} = \tilde{\Pi}_{\phi_1}U\left[\prod^{(d-1)/2}_{k=1}\Pi_{\phi_{2k}}U^\dagger
\tilde{\Pi}_{\phi_{2k+1}}U\right]\Pi_{\phi_{d+1}}.
And when the number of projector-controlled phase shifts is odd (:math:`d` is even):
.. math::
U_{QSVT} = \left[\prod^{d/2}_{k=1}\Pi_{\phi_{2k-1}}U^\dagger\tilde{\Pi}_{\phi_{2k}}U\right]
\Pi_{\phi_{d+1}}.
This circuit applies a polynomial transformation (:math:`Poly^{SV}`) to the singular values of
the block encoded matrix:
.. math::
\begin{align}
U_{QSVT}(A, \vec{\phi}) &=
\begin{bmatrix}
Poly^{SV}(A) & \cdot \\
\cdot & \cdot
\end{bmatrix}.
\end{align}
.. seealso::
:func:`~.qsvt` and `A Grand Unification of Quantum Algorithms <https://arxiv.org/pdf/2105.02859.pdf>`_.
Args:
UA (Operator): the block encoding circuit, specified as an :class:`~.Operator`,
like :class:`~.BlockEncode`
projectors (Sequence[Operator]): a list of projector-controlled phase
shifts that implement the desired polynomial
Raises:
ValueError: if the input block encoding is not an operator
**Example**
To implement QSVT in a circuit, we can use the following method:
>>> dev = qml.device("default.qubit", wires=[0])
>>> block_encoding = qml.Hadamard(wires=0) # note H is a block encoding of 1/sqrt(2)
>>> phase_shifts = [qml.RZ(-2 * theta, wires=0) for theta in (1.23, -0.5, 4)] # -2*theta to match convention
>>>
>>> @qml.qnode(dev)
>>> def example_circuit():
... qml.QSVT(block_encoding, phase_shifts)
... return qml.expval(qml.Z(0))
>>>
>>> example_circuit()
0.5403023058681395
We can visualize the circuit as follows:
>>> print(qml.draw(example_circuit)())
0: ──QSVT─┤ <Z>
To see the implementation details, we can expand the circuit:
>>> q_script = qml.tape.QuantumScript(ops=[qml.QSVT(block_encoding, phase_shifts)])
>>> print(q_script.expand().draw(decimals=2))
0: ──RZ(-2.46)──H──RZ(1.00)──H†──RZ(-8.00)─┤
See the Usage Details section for more examples on implementing QSVT with different block
encoding methods.
.. details::
:title: Usage Details
The QSVT operation can be used with different block encoding methods, depending on the
initial operator for which the singular value transformation is applied and the desired
backend device. Examples are provided below.
If we want to transform the singular values of a matrix,
the matrix can be block-encoded with either the :class:`~.BlockEncode` or :class:`~.FABLE`
operations. Note that :class:`~.BlockEncode` is more efficient on simulator devices but
it cannot be used with hardware backends because it currently has no gate decomposition.
The :class:`~.FABLE` operation is less efficient on simulator devices but is hardware
compatible.
The following example applies the polynomial :math:`p(x) = -x + 0.5x^3 + 0.5x^5` to an
arbitrary hermitian matrix using :class:`~.BlockEncode` for block encoding.
.. code-block::
poly = np.array([0, -1, 0, 0.5, 0, 0.5])
angles = qml.poly_to_angles(poly, "QSVT")
input_matrix = np.array([[0.2, 0.1], [0.1, -0.1]])
wires = [0, 1]
block_encode = qml.BlockEncode(input_matrix, wires=wires)
projectors = [
qml.PCPhase(angles[i], dim=len(input_matrix), wires=wires)
for i in range(len(angles))
]
dev = qml.device("default.qubit")
@qml.qnode(dev)
def circuit():
qml.QSVT(block_encode, projectors)
return qml.state()
.. code-block:: pycon
>>> circuit()
array([-0.194205 +0.66654551j, -0.097905 +0.35831418j,
0.3319832 -0.51047262j, -0.09551437+0.01043668j])
If we want to transform the singular values of a linear
combination of unitaries, e.g., a Hamiltonian, it can be block-encoded with operations
such as :class:`~.PrepSelPrep` or :class:`~.Qubitization`. Note that both of these operations
have a gate decomposition and can be implemented on hardware. The following example applies the polynomial
:math:`p(x) = -x + 0.5x^3 + 0.5x^5` to the Hamiltonian :math:`H = 0.1X_3 - 0.7X_3Z_4 - 0.2Z_3Y_4`,
block-encoded with :class:`~.PrepSelPrep`.
.. code-block::
poly = np.array([0, -1, 0, 0.5, 0, 0.5])
H = 0.1 * qml.X(2) - 0.7 * qml.X(2) @ qml.Z(3) - 0.2 * qml.Z(2)
control_wires = [0, 1]
block_encode = qml.PrepSelPrep(H, control=control_wires)
angles = qml.poly_to_angles(poly, "QSVT")
projectors = [
qml.PCPhase(angles[i], dim=2 ** len(H.wires), wires=control_wires + H.wires)
for i in range(len(angles))
]
dev = qml.device("default.qubit")
@qml.qnode(dev)
def circuit():
qml.QSVT(block_encode, projectors)
return qml.state()
.. code-block:: pycon
>>> circuit()
array([ 1.44000000e-01+1.01511390e-01j, 0.00000000e+00+0.00000000e+00j,
4.32000000e-01+3.04534169e-01j, 0.00000000e+00+0.00000000e+00j,
1.92998954e-17+5.00377363e-17j, 0.00000000e+00+0.00000000e+00j,
5.59003542e-01+9.65699229e-02j, 0.00000000e+00+0.00000000e+00j,
4.22566958e-01+7.30000000e-02j, 0.00000000e+00+0.00000000e+00j,
-3.16925218e-01-5.47500000e-02j, 0.00000000e+00+0.00000000e+00j,
-2.98448441e-17-3.10878188e-17j, 0.00000000e+00+0.00000000e+00j,
-2.79501771e-01-4.82849614e-02j, 0.00000000e+00+0.00000000e+00j])
"""
num_wires = AnyWires
"""int: Number of wires that the operator acts on."""
grad_method = None
"""Gradient computation method."""
def _flatten(self):
data = (self.hyperparameters["UA"], self.hyperparameters["projectors"])
return data, tuple()
@classmethod
def _primitive_bind_call(cls, *args, **kwargs):
return cls._primitive.bind(*args, **kwargs)
@classmethod
def _unflatten(cls, data, _) -> "QSVT":
return cls(*data)
def __init__(self, UA, projectors, id=None):
if not isinstance(UA, qml.operation.Operator):
raise ValueError("Input block encoding must be an Operator")
self._hyperparameters = {
"UA": UA,
"projectors": projectors,
}
total_wires = qml.wires.Wires.all_wires(
[proj.wires for proj in projectors]
) + qml.wires.Wires(UA.wires)
super().__init__(wires=total_wires, id=id)
[docs] def map_wires(self, wire_map: dict):
# pylint: disable=protected-access
new_op = copy.deepcopy(self)
new_op._wires = Wires([wire_map.get(wire, wire) for wire in self.wires])
new_op._hyperparameters["UA"] = qml.map_wires(new_op._hyperparameters["UA"], wire_map)
new_op._hyperparameters["projectors"] = [
qml.map_wires(proj, wire_map) for proj in new_op._hyperparameters["projectors"]
]
return new_op
@property
def data(self):
r"""Flattened list of operator data in this QSVT operation.
This ensures that the backend of a ``QuantumScript`` which contains a
``QSVT`` operation can be inferred with respect to the types of the
``QSVT`` block encoding and projector-controlled phase shift data.
"""
return tuple(datum for op in self._operators for datum in op.data)
@data.setter
def data(self, new_data):
# We need to check if ``new_data`` is empty because ``Operator.__init__()`` will attempt to
# assign the QSVT data to an empty tuple (since no positional arguments are provided).
if new_data:
for op in self._operators:
if op.num_params > 0:
op.data = new_data[: op.num_params]
new_data = new_data[op.num_params :]
def __copy__(self):
# Override Operator.__copy__() to avoid setting the "data" property before the new instance
# is assigned hyper-parameters since QSVT data is derived from the hyper-parameters.
clone = QSVT.__new__(QSVT)
# Ensure the operators in the hyper-parameters are copied instead of aliased.
clone._hyperparameters = {
"UA": copy.copy(self._hyperparameters["UA"]),
"projectors": list(map(copy.copy, self._hyperparameters["projectors"])),
}
for attr, value in vars(self).items():
if attr != "_hyperparameters":
setattr(clone, attr, value)
return clone
@property
def _operators(self) -> list[qml.operation.Operator]:
"""Flattened list of operators that compose this QSVT operation."""
return [self._hyperparameters["UA"], *self._hyperparameters["projectors"]]
[docs] @staticmethod
def compute_decomposition(
*_data, UA, projectors, **_kwargs
): # pylint: disable=arguments-differ
r"""Representation of the operator as a product of other operators.
The :class:`~.QSVT` is decomposed into alternating block encoding
and projector-controlled phase shift operators. This is defined by the following
equations, where :math:`U` is the block encoding operation and both :math:`\Pi_\phi` and
:math:`\tilde{\Pi}_\phi` are projector-controlled phase shifts with angle :math:`\phi`.
When the number of projector-controlled phase shifts is even (:math:`d` is odd), the QSVT
circuit is defined as:
.. math::
U_{QSVT} = \Pi_{\phi_1}U\left[\prod^{(d-1)/2}_{k=1}\Pi_{\phi_{2k}}U^\dagger
\tilde{\Pi}_{\phi_{2k+1}}U\right]\Pi_{\phi_{d+1}}.
And when the number of projector-controlled phase shifts is odd (:math:`d` is even):
.. math::
U_{QSVT} = \left[\prod^{d/2}_{k=1}\Pi_{\phi_{2k-1}}U^\dagger\tilde{\Pi}_{\phi_{2k}}U\right]
\Pi_{\phi_{d+1}}.
.. seealso:: :meth:`~.QSVT.decomposition`.
Args:
UA (Operator): the block encoding circuit, specified as a :class:`~.Operator`
projectors (list[Operator]): a list of projector-controlled phase
shift circuits that implement the desired polynomial
Returns:
list[.Operator]: decomposition of the operator
"""
op_list = []
UA_adj = copy.copy(UA)
for idx, op in enumerate(projectors[:-1]):
if qml.QueuingManager.recording():
qml.apply(op)
op_list.append(op)
if idx % 2 == 0:
if qml.QueuingManager.recording():
qml.apply(UA)
op_list.append(UA)
else:
op_list.append(adjoint(UA_adj))
if qml.QueuingManager.recording():
qml.apply(projectors[-1])
op_list.append(projectors[-1])
return op_list
[docs] def label(self, decimals=None, base_label=None, cache=None):
op_label = base_label or self.__class__.__name__
return op_label
[docs] def queue(self, context=QueuingManager):
context.remove(self._hyperparameters["UA"])
for op in self._hyperparameters["projectors"]:
context.remove(op)
context.append(self)
return self
[docs] @staticmethod
def compute_matrix(*args, **kwargs):
r"""Representation of the operator as a canonical matrix in the computational basis (static method).
The canonical matrix is the textbook matrix representation that does not consider wires.
Implicitly, this assumes that the wires of the operator correspond to the global wire order.
.. seealso:: :meth:`~.Operator.matrix` and :func:`~.matrix`
Args:
*params (list): trainable parameters of the operator, as stored in the ``parameters`` attribute
**hyperparams (dict): non-trainable hyperparameters of the operator, as stored in the ``hyperparameters`` attribute
Returns:
tensor_like: matrix representation
"""
# pylint: disable=unused-argument
op_list = []
UA = kwargs["UA"]
projectors = kwargs["projectors"]
with (
QueuingManager.stop_recording()
): # incase this method is called in a queue context, this prevents
UA_copy = copy.copy(UA) # us from queuing operators unnecessarily
for idx, op in enumerate(projectors[:-1]):
op_list.append(op)
if idx % 2 == 0:
op_list.append(UA)
else:
op_list.append(adjoint(UA_copy))
op_list.append(projectors[-1])
mat = qml.matrix(qml.prod(*tuple(op_list[::-1])))
return mat
def _complementary_poly(poly_coeffs):
r"""
Computes the complementary polynomial Q given a polynomial P.
The polynomial Q is complementary to P if it satisfies the following equation:
.. math:
|P(e^{i\theta})|^2 + |Q(e^{i\theta})|^2 = 1, \quad \forall \quad \theta \in \left[0, 2\pi\right]
The method is based on computing an auxiliary polynomial R, finding its roots,
and reconstructing Q by using information extracted from the roots.
For more details see `arXiv:2308.01501 <https://arxiv.org/abs/2308.01501>`_.
Args:
poly_coeffs (tensor-like): coefficients of the complex polynomial P
Returns:
tensor-like: coefficients of the complementary polynomial Q
"""
poly_degree = len(poly_coeffs) - 1
# Build the polynomial R(z) = z^degree * (1 - conj(P(1/z)) * P(z)), deduced from (eq.33) and (eq.34) of
# `arXiv:2308.01501 <https://arxiv.org/abs/2308.01501>`_.
# Note that conj(P(1/z)) * P(z) could be expressed as z^-degree * conj(P(z)[::-1]) * P(z)
R = Polynomial.basis(poly_degree) - Polynomial(poly_coeffs) * Polynomial(
np.conj(poly_coeffs[::-1])
)
r_roots = R.roots()
inside_circle = [root for root in r_roots if np.abs(root) <= 1]
outside_circle = [root for root in r_roots if np.abs(root) > 1]
scale_factor = np.sqrt(np.abs(R.coef[-1] * np.prod(outside_circle)))
Q_poly = scale_factor * Polynomial.fromroots(inside_circle)
return Q_poly.coef
def _compute_qsp_angle(poly_coeffs):
r"""
Computes the Quantum Signal Processing (QSP) angles given the coefficients of a polynomial F.
The method for computing the QSP angles is adapted from the approach described in [`arXiv:2406.04246 <https://arxiv.org/abs/2406.04246>`_] for Generalized-QSP.
Args:
poly_coeffs (tensor-like): coefficients of the input polynomial F
Returns:
(tensor-like): QSP angles corresponding to the input polynomial F
.. details::
:title: Implementation details
Based on the appendix A in `arXiv:2406.04246 <https://arxiv.org/abs/2406.04246>`_, the target polynomial :math:`F`
is transformed into a new polynomial :math:`P` by following the steps below:
0. The input to the function are the coefficients of :math:`F`, e.g. :math:`[c_0, 0, c_1, 0, c_2]`.
1. We express the polynomial in the Chebyshev basis by applying the Chebyshev transform. This generates
a new representation :math:`[a_0, 0, a_1, 0, a_2]`.
2. We generate :math:`P` by reordering the array, moving the zeros to the initial positions.
:math:`P = [0, 0, a_0, a_1, a_2]`.
The polynomial :math:`P` can now be used in Algorithm 1 of [`arXiv:2308.01501 <https://arxiv.org/abs/2308.01501>`_]
in order to find the desired angles.
The above algorithm is specific to Generalized-QSP so an adaptation has been made to return the required angles:
- The :math:`R(\theta, \phi, \lambda)` gate, is simplified into a :math:`R_Y(\theta)` gate.
- The calculation of :math:`\theta_d` is updated to :math:`\theta_d = \tan^{-1}(\frac{a_d}{b_d})`.
In this way, the sign introduced by :math:`\phi_d` and :math:`\lambda_d` is absorbed
in the :math:`\theta_d` value.
"""
parity = (len(poly_coeffs) - 1) % 2
P = np.concatenate(
[np.zeros(len(poly_coeffs) // 2), chebyshev.poly2cheb(poly_coeffs)[parity::2]]
) * (1 - 1e-12)
complementary = _complementary_poly(P)
polynomial_matrix = np.array([P, complementary])
num_terms = polynomial_matrix.shape[1]
rotation_angles = np.zeros(num_terms)
# Adaptation of Algorithm 1 of [arXiv:2308.01501]
with qml.QueuingManager.stop_recording():
for idx in range(num_terms - 1, -1, -1):
poly_a, poly_b = polynomial_matrix[:, idx]
rotation_angles[idx] = np.arctan2(poly_b.real, poly_a.real)
rotation_op = qml.matrix(qml.RY(-2 * rotation_angles[idx], wires=0))
updated_poly_matrix = rotation_op @ polynomial_matrix
polynomial_matrix = np.array(
[updated_poly_matrix[0][1 : idx + 1], updated_poly_matrix[1][0:idx]]
)
return rotation_angles
def _gqsp_u3_gate(theta, phi, lambd):
r"""
Computes the U3 gate matrix for Generalized Quantum Signal Processing (GQSP) as described
in [`arXiv:2406.04246 <https://arxiv.org/abs/2406.04246>`_]
"""
exp_phi = qml.math.exp(1j * phi)
exp_lambda = qml.math.exp(1j * lambd)
exp_lambda_phi = qml.math.exp(1j * (lambd + phi))
matrix = np.array(
[
[exp_lambda_phi * qml.math.cos(theta), exp_phi * qml.math.sin(theta)],
[exp_lambda * qml.math.sin(theta), -qml.math.cos(theta)],
],
dtype=complex,
)
return matrix
def _compute_gqsp_angles(poly_coeffs):
r"""
Computes the Generalized Quantum Signal Processing (GQSP) angles given the coefficients of a polynomial P.
The method for computing the GQSP angles is based on the algorithm described in [`arXiv:2406.04246 <https://arxiv.org/abs/2406.04246>`_].
The complementary polynomial is calculated using root-finding methods.
Args:
poly_coeffs (tensor-like): Coefficients of the input polynomial P.
Returns:
(tensor-like): QSP angles corresponding to the input polynomial P. The shape is (3, P-degree)
"""
complementary = _complementary_poly(poly_coeffs)
# Algorithm 1 in [arXiv:2308.01501]
input_data = qml.math.array([poly_coeffs, complementary])
num_elements = input_data.shape[1]
angles_theta, angles_phi, angles_lambda = qml.math.zeros([3, num_elements])
for idx in range(num_elements - 1, -1, -1):
component_a, component_b = input_data[:, idx]
angles_theta[idx] = qml.math.arctan2(np.abs(component_b), qml.math.abs(component_a))
angles_phi[idx] = (
0
if qml.math.isclose(qml.math.abs(component_b), 0, atol=1e-10)
else qml.math.angle(component_a * qml.math.conj(component_b))
)
if idx == 0:
angles_lambda[0] = qml.math.angle(component_b)
else:
updated_matrix = (
_gqsp_u3_gate(angles_theta[idx], angles_phi[idx], 0).conj().T @ input_data
)
input_data = qml.math.array([updated_matrix[0][1 : idx + 1], updated_matrix[1][0:idx]])
return angles_theta, angles_phi, angles_lambda
[docs]def transform_angles(angles, routine1, routine2):
r"""
Converts angles for quantum signal processing (QSP) and quantum singular value transformation (QSVT) routines.
The transformation is based on Appendix A.2 of `arXiv:2105.02859 <https://arxiv.org/abs/2105.02859>`_. Note that QSVT is equivalent to taking the reflection convention of QSP.
Args:
angles (tensor-like): angles to be transformed
routine1 (str): the current routine for which the angles are obtained, must be either ``"QSP"`` or ``"QSVT"``
routine2 (str): the target routine for which the angles should be transformed,
must be either ``"QSP"`` or ``"QSVT"``
Returns:
tensor-like: the transformed angles as an array
**Example**
.. code-block::
>>> qsp_angles = np.array([0.2, 0.3, 0.5])
>>> qsvt_angles = qml.transform_angles(qsp_angles, "QSP", "QSVT")
>>> print(qsvt_angles)
[-6.86858347 1.87079633 -0.28539816]
.. details::
:title: Usage Details
This example applies the polynomial :math:`P(x) = x - \frac{x^3}{2} + \frac{x^5}{3}` to a block-encoding
of :math:`x = 0.2`.
.. code-block::
poly = np.array([0, 1.0, 0, -1/2, 0, 1/3])
qsp_angles = qml.poly_to_angles(poly, "QSP")
qsvt_angles = qml.transform_angles(qsp_angles, "QSP", "QSVT")
x = 0.2
# Encodes x in the top left of the matrix
block_encoding = qml.RX(2 * np.arccos(x), wires=0)
projectors = [qml.PCPhase(angle, dim=1, wires=0) for angle in qsvt_angles]
@qml.qnode(qml.device("default.qubit"))
def circuit_qsvt():
qml.QSVT(block_encoding, projectors)
return qml.state()
output = qml.matrix(circuit_qsvt, wire_order=[0])()[0, 0]
expected = sum(coef * (x**i) for i, coef in enumerate(poly))
print("output qsvt: ", output.real)
print("P(x) = ", expected)
.. code-block:: pycon
output qsvt: 0.19610666666647059
P(x) = 0.19610666666666668
"""
if routine1 == "QSP" and routine2 == "QSVT":
num_angles = len(angles)
update_vals = np.empty(num_angles)
update_vals[0] = 3 * np.pi / 4 - (3 + num_angles % 4) * np.pi / 2
update_vals[1:-1] = np.pi / 2
update_vals[-1] = -np.pi / 4
update_vals = qml.math.convert_like(update_vals, angles)
return angles + update_vals
if routine1 == "QSVT" and routine2 == "QSP":
num_angles = len(angles)
update_vals = np.empty(num_angles)
update_vals[0] = 3 * np.pi / 4 - (3 + num_angles % 4) * np.pi / 2
update_vals[1:-1] = np.pi / 2
update_vals[-1] = -np.pi / 4
update_vals = qml.math.convert_like(update_vals, angles)
return angles - update_vals
raise AssertionError(
f"Invalid conversion. The conversion between {routine1} --> {routine2} is not defined."
)
[docs]def poly_to_angles(poly, routine, angle_solver: Literal["root-finding"] = "root-finding"):
r"""
Computes the angles needed to implement a polynomial transformation with quantum signal processing (QSP)
or quantum singular value transformation (QSVT).
The polynomial :math:`P(x) = \sum_n a_n x^n` must satisfy :math:`|P(x)| \leq 1` for :math:`x \in [-1, 1]`, the
coefficients :math:`a_n` must be real and the exponents :math:`n` must be either all even or all odd.
For more details see `arXiv:2105.02859 <https://arxiv.org/abs/2105.02859>`_.
Args:
poly (tensor-like): coefficients of the polynomial ordered from lowest to highest power
routine (str): the routine for which the angles are computed. Must be either ``"QSP"`` or ``"QSVT"``.
angle_solver (str): the method used to calculate the angles. Default is ``"root-finding"``.
``"root-finding"`` is currently the only supported method.
Returns:
(tensor-like): computed angles for the specified routine
Raises:
AssertionError: if ``poly`` is not valid
AssertionError: if ``routine`` or ``angle_solver`` is not supported
**Example**
This example generates the ``QSVT`` angles for the polynomial :math:`P(x) = x - \frac{x^3}{2} + \frac{x^5}{3}`.
.. code-block::
>>> poly = np.array([0, 1.0, 0, -1/2, 0, 1/3])
>>> qsvt_angles = qml.poly_to_angles(poly, "QSVT")
>>> print(qsvt_angles)
[-5.49778714 1.57079633 1.57079633 0.5833829 1.61095884 0.74753829]
.. details::
:title: Usage Details
This example applies the polynomial :math:`P(x) = x - \frac{x^3}{2} + \frac{x^5}{3}` to a block-encoding
of :math:`x = 0.2`.
.. code-block::
poly = np.array([0, 1.0, 0, -1/2, 0, 1/3])
qsvt_angles = qml.poly_to_angles(poly, "QSVT")
x = 0.2
# Encode x in the top left of the matrix
block_encoding = qml.RX(2 * np.arccos(x), wires=0)
projectors = [qml.PCPhase(angle, dim=1, wires=0) for angle in qsvt_angles]
@qml.qnode(qml.device("default.qubit"))
def circuit_qsvt():
qml.QSVT(block_encoding, projectors)
return qml.state()
output = qml.matrix(circuit_qsvt, wire_order=[0])()[0, 0]
expected = sum(coef * (x**i) for i, coef in enumerate(poly))
print("output qsvt: ", output.real)
print("P(x) = ", expected)
.. code-block:: pycon
output qsvt: 0.19610666666647059
P(x) = 0.19610666666666668
"""
poly = qml.math.trim_zeros(qml.math.array(poly, like="numpy"), trim="b")
if len(poly) == 1:
raise AssertionError("The polynomial must have at least degree 1.")
for x in [-1, 0, 1]:
if qml.math.abs(qml.math.sum(coeff * x**i for i, coeff in enumerate(poly))) > 1:
# Check that |P(x)| ≤ 1. Only points -1, 0, 1 will be checked.
raise AssertionError("The polynomial must satisfy that |P(x)| ≤ 1 for all x in [-1, 1]")
if routine in ["QSVT", "QSP"]:
if not (
np.isclose(qml.math.sum(qml.math.abs(poly[::2])), 0.0)
or np.isclose(qml.math.sum(qml.math.abs(poly[1::2])), 0.0)
):
raise AssertionError(
"The polynomial has no definite parity. All odd or even entries in the array must take a value of zero."
)
assert np.allclose(
np.array(poly, dtype=np.complex128).imag, 0
), "Array must not have an imaginary part"
if routine == "QSVT":
if angle_solver == "root-finding":
return transform_angles(_compute_qsp_angle(poly), "QSP", "QSVT")
raise AssertionError("Invalid angle solver method. We currently support 'root-finding'")
if routine == "QSP":
if angle_solver == "root-finding":
return _compute_qsp_angle(poly)
raise AssertionError("Invalid angle solver method. Valid value is 'root-finding'")
if routine == "GQSP":
return _compute_gqsp_angles(poly)
raise AssertionError("Invalid routine. Valid values are 'QSP' and 'QSVT'")
_modules/pennylane/templates/subroutines/qsvt
Download Python script
Download Notebook
View on GitHub