Source code for pennylane.templates.state_preparations.mottonen
# Copyright 2018-2021 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.
r"""
Contains the MottonenStatePreparation template.
"""
from typing import Optional
import numpy as np
import pennylane as qml
from pennylane.operation import Operation
from pennylane.typing import TensorLike
def gray_code(rank):
"""Generates the
`Gray code <https://en.wikipedia.org/wiki/Gray_code>`__
of given rank, as numeric output.
Args:
rank (int): rank of the Gray code (i.e. number of bits)
Returns:
np.ndarray[int]: Array of ``2**rank`` integers that make up the Gray code.
"""
g = np.array([0, 1])
for i in range(1, rank):
g = np.concatenate([g, g[::-1] + 2**i])
return g
_walsh_hadamard_matrix = np.array([[1, 1], [1, -1]]) / 2
_cnot_matrix = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]]).reshape((2,) * 4)
def compute_theta(alpha: TensorLike, num_qubits: Optional[int] = None):
r"""Maps the input angles ``alpha`` of the multi-controlled rotations decomposition of a
uniformly controlled rotation to the rotation angles used in the
`Gray code <https://en.wikipedia.org/wiki/Gray_code>`__ implementation.
This function uses the fact that the transformation given by Eq. (3) in
`Möttönen et al. (2004) <https://arxiv.org/abs/quant-ph/0407010>`_ is equal to a Walsh-Hadamard
transform followed by some permutations, which can be expressed as a ladder of CNOT gates
applied to the angles, when interpreting them as a quantum state.
Args:
alpha (tensor_like): The array or tensor to be transformed. Must have a length that
is a power of two.
num_qubits (int): Number of qubits. If not given, it will be computed from ``alpha``.
If given, it should match the trailing dimension of ``alpha``.
Returns:
tensor_like: The transformed tensor with the same shape as the input ``alpha``.
Due to the execution of the transform as a sequence of tensor multiplications
with shapes ``(2, 2), (2, 2,... 2)->(2, 2,... 2)``, the theoretical scaling of this
method is the same as the one for the
`Fast Walsh-Hadamard transform <https://en.wikipedia.org/wiki/Fast_Walsh-Hadamard_transform>`__:
On :math:`n` qubits, there are :math:`n` calls to ``tensordot``, each multiplying a ``(2, 2)``
matrix to a ``(2,)*num_qubits`` vector, with a single axis being contracted. This means
that there are :math:`n` operations with a floating point operation count of
``4 * 2**(num_qubits-1)``, where ``4`` is the cost of a single ``(2, 2) @ (2,)`` contraction
and ``2**(n-1)`` is the number of copies due to the non-contracted :math:`n-1` axes.
Due to the large internal speedups of compiled matrix multiplication and compatibility
with autodifferentiation frameworks, the approach taken here is favourable over a manual
realization of the FWHT unless memory limitations restrict the creation of intermediate
arrays, which would make in-place techniques favourable.
Similarly, the permutation can be applied by contracting the angles with the reshaped CNOT
matrix.
"""
orig_shape = qml.math.shape(alpha)
num_qubits = num_qubits or int(qml.math.log2(orig_shape[-1]))
if num_qubits == 0:
# No processing occurs for num_qubits=0
return alpha
# Reshape the array so that we may apply the Hadamard transform to each axis individually
if broadcasted := len(orig_shape) > 1:
new_shape = (orig_shape[0],) + (2,) * num_qubits
else:
new_shape = (2,) * num_qubits
alpha = qml.math.reshape(alpha, new_shape)
# Apply Hadamard transform to each axis, shifted by one for broadcasting
for i in range(broadcasted, num_qubits + broadcasted):
alpha = qml.math.tensordot(_walsh_hadamard_matrix, alpha, axes=[[1], [i]])
# The axes are now in the ordering [qubit n-1, qubit n-2, ..., qubit 1, qubit 0, batch]
if num_qubits > 1:
# If there is more than one qubit, we need to reorder the angles, according to applying
# the CNOT ladder [CNOT([i, i+1]) for i in range(num_qubits-1)]
# The first CNOT thus targets the zeroth and first qubit, axes n-1 and n-2 (see above)
alpha = qml.math.tensordot(
_cnot_matrix, alpha, axes=[[2, 3], [num_qubits - 1, num_qubits - 2]]
)
# The axes are now ordered as [qubit 0, qubit 1, qubit n-1, qubit n-2, ..., qubit 2, batch]
# Following CNOTs use the same axes: the next control qubit (previous target qubit) always
# is in position ``1`` and the next target qubit always is the last qubit axis
# (``num_qubits-1``). For example, the first loop iteration moves the axes into positions
# [qubit 1, qubit 2, qubit 0, qubit n-1, qubit n-2, ... ,qubit 3, batch]
# and the iteration after that moves them to
# [qubit 2, qubit 3, qubit 1, qubit 0, qubit n-1, qubit n-2, ... ,qubit 4, batch]
for i in range(broadcasted + 1, num_qubits + broadcasted - 1):
alpha = qml.math.tensordot(_cnot_matrix, alpha, axes=[[2, 3], [1, num_qubits - 1]])
# In the end, we exchange the first two axes because we have the axes ordering
# [qubit n-2, qubit n-1, qubit n-3, qubit n-4, ... qubit 1, qubit 0, batch]
alpha = qml.math.moveaxis(alpha, 0, 1)
# Finally, the axis ordering has to be flipped entirely, moving the batch to the front
# and the qubits into the right ordering, [batch, qubit 0, qubit 1, ..., qubit n-1]
# For num_qubits=1 we just exchange the single qubit axis and the batching axis
return qml.math.reshape(qml.math.transpose(alpha), orig_shape)
def _uniform_rotation_dagger_ops(gate, alpha, control_wires, target_wire):
r"""Returns a list of operators that applies a uniformly-controlled rotation to the target qubit.
Args:
gate (.Operation): gate to be applied, needs to have exactly one parameter
alpha (tensor_like): angles to decompose the uniformly-controlled rotation into multi-controlled rotations
control_wires (array[int]): wires that act as control
target_wire (int): wire that acts as target
Returns:
list[.Operator]: sequence of operators defined by this function
"""
with qml.queuing.AnnotatedQueue() as q:
_apply_uniform_rotation_dagger(gate, alpha, control_wires, target_wire)
if qml.queuing.QueuingManager.recording():
for op in q.queue:
qml.apply(op)
return q.queue
def _apply_uniform_rotation_dagger(gate, alpha, control_wires, target_wire):
r"""Applies a uniformly-controlled rotation to the target qubit.
A uniformly-controlled rotation is a sequence of multi-controlled
rotations, each of which is conditioned on the control qubits being in a different state.
For example, a uniformly-controlled rotation with two control qubits describes a sequence of
four multi-controlled rotations, each applying the rotation only if the control qubits
are in states :math:`|00\rangle`, :math:`|01\rangle`, :math:`|10\rangle`, and :math:`|11\rangle`, respectively.
To implement a uniformly-controlled rotation using single qubit rotations and CNOT gates,
a decomposition based on Gray codes is used. For this purpose, the multi-controlled rotation
angles alpha have to be converted into a set of non-controlled rotation angles theta.
For more details, see `Möttönen and Vartiainen (2005), Fig 7a<https://arxiv.org/abs/quant-ph/0504100>`_.
Args:
gate (.Operation): gate to be applied, needs to have exactly one parameter
alpha (tensor_like): angles to decompose the uniformly-controlled rotation into multi-controlled rotations
control_wires (array[int]): wires that act as control
target_wire (int): wire that acts as target
"""
gray_code_rank = len(control_wires)
theta = compute_theta(alpha, num_qubits=gray_code_rank)
if gray_code_rank == 0:
if (
qml.math.is_abstract(theta)
or qml.math.requires_grad(theta)
or qml.math.all(theta[..., 0] != 0.0)
):
gate(theta[..., 0], wires=[target_wire])
return
code = gray_code(gray_code_rank)
control_indices = np.log2(code ^ np.roll(code, -1)).astype(int)
# For abstract or differentiated theta we will never skip a rotation. Likewise, if there
# is at least one non-zero angle (per batch if batched) for all rotations.
skip_none = qml.math.is_abstract(theta) or qml.math.requires_grad(theta)
if not skip_none:
nonzero = (
(theta != 0.0) if qml.math.ndim(theta) == 1 else qml.math.any(theta != 0.0, axis=0)
)
skip_none = qml.math.all(nonzero)
for i, control_index in enumerate(control_indices):
# If we do not _never_ skip, we might skip _some_ rotation
if skip_none or qml.math.all(theta[..., i] != 0.0):
gate(theta[..., i], wires=[target_wire])
qml.CNOT(wires=[control_wires[control_index], target_wire])
def _uniform_rotation_dagger_ops(gate, alpha, control_wires, target_wire):
r"""Returns a list of operators that applies a uniformly-controlled rotation to the target qubit.
Args:
gate (.Operation): gate to be applied, needs to have exactly one parameter
alpha (tensor_like): angles to decompose the uniformly-controlled rotation into multi-controlled rotations
control_wires (array[int]): wires that act as control
target_wire (int): wire that acts as target
Returns:
list[.Operator]: sequence of operators defined by this function
"""
with qml.queuing.AnnotatedQueue() as q:
_apply_uniform_rotation_dagger(gate, alpha, control_wires, target_wire)
if qml.queuing.QueuingManager.recording():
for op in q.queue:
qml.apply(op)
return q.queue
def _get_alpha_z(omega, n, k):
r"""Computes the rotation angles required to implement the uniformly-controlled Z rotation
applied to the :math:`k`th qubit.
The :math:`j`th angle is related to the phases omega of the desired amplitudes via:
.. math:: \alpha^{z,k}_j = \sum_{l=1}^{2^{k-1}} \frac{\omega_{(2j-1) 2^{k-1}+l} - \omega_{(2j-2) 2^{k-1}+l}}{2^{k-1}}
Args:
omega (tensor_like): phases of the state to prepare
n (int): total number of qubits for the uniformly-controlled rotation
k (int): index of current qubit
Returns:
array representing :math:`\alpha^{z,k}`
"""
indices1 = [
[(2 * j - 1) * 2 ** (k - 1) + l - 1 for l in range(1, 2 ** (k - 1) + 1)]
for j in range(1, 2 ** (n - k) + 1)
]
indices2 = [
[(2 * j - 2) * 2 ** (k - 1) + l - 1 for l in range(1, 2 ** (k - 1) + 1)]
for j in range(1, 2 ** (n - k) + 1)
]
term1 = qml.math.take(omega, indices=indices1, axis=-1)
term2 = qml.math.take(omega, indices=indices2, axis=-1)
diff = (term1 - term2) / 2 ** (k - 1)
return qml.math.sum(diff, axis=-1)
def _get_alpha_y(a, n, k):
r"""Computes the rotation angles required to implement the uniformly controlled Y rotation
applied to the :math:`k`th qubit.
The :math:`j`-th angle is related to the absolute values, a, of the desired amplitudes via:
.. math:: \alpha^{y,k}_j = 2 \arcsin \sqrt{ \frac{ \sum_{l=1}^{2^{k-1}} a_{(2j-1)2^{k-1} +l}^2 }{ \sum_{l=1}^{2^{k}} a_{(j-1)2^{k} +l}^2 } }
Args:
a (tensor_like): absolute values of the state to prepare
n (int): total number of qubits for the uniformly-controlled rotation
k (int): index of current qubit
Returns:
array representing :math:`\alpha^{y,k}`
"""
indices_numerator = [
[(2 * (j + 1) - 1) * 2 ** (k - 1) + l for l in range(2 ** (k - 1))]
for j in range(2 ** (n - k))
]
numerator = qml.math.take(a, indices=indices_numerator, axis=-1)
numerator = qml.math.sum(qml.math.abs(numerator) ** 2, axis=-1)
indices_denominator = [[j * 2**k + l for l in range(2**k)] for j in range(2 ** (n - k))]
denominator = qml.math.take(a, indices=indices_denominator, axis=-1)
denominator = qml.math.sum(qml.math.abs(denominator) ** 2, axis=-1)
# Divide only where denominator is zero, else leave initial value of zero.
# The equation guarantees that the numerator is also zero in the corresponding entries.
with np.errstate(divide="ignore", invalid="ignore"):
division = numerator / denominator
# Cast the numerator and denominator to ensure compatibility with interfaces
division = qml.math.cast(division, np.float64)
denominator = qml.math.cast(denominator, np.float64)
division = qml.math.where(denominator != 0.0, division, 0.0)
return 2 * qml.math.arcsin(qml.math.sqrt(division))
[docs]
class MottonenStatePreparation(Operation):
r"""
Prepares an arbitrary state on the given wires using a decomposition into gates developed
by `Möttönen et al. (2004) <https://arxiv.org/abs/quant-ph/0407010>`_.
The state is prepared via a sequence
of uniformly controlled rotations. A uniformly controlled rotation on a target qubit is
composed from all possible controlled rotations on the qubit and can be used to address individual
elements of the state vector.
In the work of Möttönen et al., inverse state preparation
is executed by first equalizing the phases of the state vector via uniformly controlled Z rotations,
and then rotating the now real state vector into the direction of the state :math:`|0\rangle` via
uniformly controlled Y rotations.
This code is adapted from code written by Carsten Blank for PennyLane-Qiskit.
.. warning::
Due to non-trivial classical processing of the state vector,
this template is not always fully differentiable.
Args:
state_vector (tensor_like): Input array of shape ``(2^n,)``, where ``n`` is the number of wires
the state preparation acts on. The input array must be normalized.
wires (Iterable): wires that the template acts on
Example:
``MottonenStatePreparation`` creates any arbitrary state on the given wires depending on the input state vector.
.. code-block:: python
dev = qml.device('default.qubit', wires=3)
@qml.qnode(dev)
def circuit(state):
qml.MottonenStatePreparation(state_vector=state, wires=range(3))
return qml.state()
state = np.array([1, 2j, 3, 4j, 5, 6j, 7, 8j])
state = state / np.linalg.norm(state)
print(qml.draw(circuit, level="device", max_length=80)(state))
.. code-block::
0: ──RY(2.35)─╭●───────────╭●──────────────╭●────────────────────────╭●
1: ──RY(2.09)─╰X──RY(0.21)─╰X─╭●───────────│────────────╭●───────────│─
2: ──RY(1.88)─────────────────╰X──RY(0.10)─╰X──RY(0.08)─╰X──RY(0.15)─╰X
──╭●────────╭●────╭●────╭●─╭GlobalPhase(-0.79)─┤ ╭State
──╰X────────╰X─╭●─│──╭●─│──├GlobalPhase(-0.79)─┤ ├State
───RZ(1.57)────╰X─╰X─╰X─╰X─╰GlobalPhase(-0.79)─┤ ╰State
The state preparation can be checked by running:
>>> print(np.allclose(state, circuit(state)))
True
"""
grad_method = None
ndim_params = (1,)
def __init__(self, state_vector, wires, id=None):
# check if the `state_vector` param is batched
batched = len(qml.math.shape(state_vector)) > 1
state_batch = state_vector if batched else [state_vector]
# apply checks to each state vector in the batch
for i, state in enumerate(state_batch):
shape = qml.math.shape(state)
if len(shape) != 1:
raise ValueError(
f"State vectors must be one-dimensional; vector {i} has shape {shape}."
)
n_amplitudes = shape[0]
if n_amplitudes != 2 ** len(qml.wires.Wires(wires)):
raise ValueError(
f"State vectors must be of length {2 ** len(wires)} or less; vector {i} has length {n_amplitudes}."
)
if not qml.math.is_abstract(state):
norm = qml.math.sum(qml.math.abs(state) ** 2)
if not qml.math.allclose(norm, 1.0, atol=1e-3):
raise ValueError(
f"State vectors have to be of norm 1.0, vector {i} has squared norm {norm}"
)
super().__init__(state_vector, wires=wires, id=id)
@property
def num_params(self):
return 1
[docs]
@staticmethod
def compute_decomposition(state_vector, wires): # pylint: disable=arguments-differ
r"""Representation of the operator as a product of other operators.
.. math:: O = O_1 O_2 \dots O_n.
.. seealso:: :meth:`~.MottonenStatePreparation.decomposition`.
Args:
state_vector (tensor_like): Normalized state vector of shape ``(2^len(wires),)``
wires (Any or Iterable[Any]): wires that the operator acts on
Returns:
list[.Operator]: decomposition of the operator
**Example**
>>> state_vector = torch.tensor([0.5, 0.5, 0.5, 0.5])
>>> qml.MottonenStatePreparation.compute_decomposition(state_vector, wires=["a", "b"])
[RY(array(1.57079633), wires=['a']),
RY(array(1.57079633), wires=['b']),
CNOT(wires=['a', 'b']),
CNOT(wires=['a', 'b'])]
"""
if len(qml.math.shape(state_vector)) > 1:
raise ValueError(
"Broadcasting with MottonenStatePreparation is not supported. Please use the "
"qml.transforms.broadcast_expand transform to use broadcasting with "
"MottonenStatePreparation."
)
a = qml.math.abs(state_vector)
omega = qml.math.angle(state_vector)
# change ordering of wires, since original code
# was written for IBM machines
wires_reverse = wires[::-1]
op_list = []
# Apply inverse y rotation cascade to prepare correct absolute values of amplitudes
for k in range(len(wires_reverse), 0, -1):
alpha_y_k = _get_alpha_y(a, len(wires_reverse), k)
control = wires_reverse[k:]
target = wires_reverse[k - 1]
op_list.extend(_uniform_rotation_dagger_ops(qml.RY, alpha_y_k, control, target))
# If necessary, apply inverse z rotation cascade to prepare correct phases of amplitudes
if (
qml.math.is_abstract(omega)
or qml.math.requires_grad(omega)
or not qml.math.allclose(omega, 0)
):
for k in range(len(wires_reverse), 0, -1):
alpha_z_k = _get_alpha_z(omega, len(wires_reverse), k)
control = wires_reverse[k:]
target = wires_reverse[k - 1]
if len(alpha_z_k) > 0:
op_list.extend(_uniform_rotation_dagger_ops(qml.RZ, alpha_z_k, control, target))
global_phase = qml.math.sum(-1 * qml.math.angle(state_vector) / len(state_vector))
op_list.extend([qml.GlobalPhase(global_phase, wires=wires)])
return op_list
_modules/pennylane/templates/state_preparations/mottonen
Download Python script
Download Notebook
View on GitHub