Source code for pennylane.utils

# Copyright 2018-2021 Xanadu Quantum Technologies Inc.

# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at

# Unless required by applicable law or agreed to in writing, software
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
"""
This module contains utilities and auxiliary functions which are shared
across the PennyLane submodules.
"""
# pylint: disable=protected-access,too-many-branches
from collections.abc import Iterable
import functools
import inspect
import itertools
import numbers
from operator import matmul
import warnings

import numpy as np
import scipy

import pennylane as qml

[docs]def decompose_hamiltonian(H, hide_identity=False):
r"""Decomposes a Hermitian matrix into a linear combination of Pauli operators.

Args:
H (array[complex]): a Hermitian matrix of dimension :math:2^n\times 2^n
hide_identity (bool): does not include the :class:~.Identity observable within
the tensor products of the decomposition if True

Returns:
tuple[list[float], list[~.Observable]]: a list of coefficients and a list
of corresponding tensor products of Pauli observables that decompose the Hamiltonian.

**Example:**

We can use this function to compute the Pauli operator decomposition of an arbitrary Hermitian
matrix:

>>> A = np.array(
... [[-2, -2+1j, -2, -2], [-2-1j,  0,  0, -1], [-2,  0, -2, -1], [-2, -1, -1,  0]])
>>> coeffs, obs_list = decompose_hamiltonian(A)
>>> coeffs
[-1.0, -1.5, -0.5, -1.0, -1.5, -1.0, -0.5, 1.0, -0.5, -0.5]

We can use the output coefficients and tensor Pauli terms to construct a :class:~.Hamiltonian:

>>> H = qml.Hamiltonian(coeffs, obs_list)
>>> print(H)
(-1.0) [I0 I1]
+ (-1.5) [X1]
+ (-0.5) [Y1]
+ (-1.0) [Z1]
+ (-1.5) [X0]
+ (-1.0) [X0 X1]
+ (-0.5) [X0 Z1]
+ (1.0) [Y0 Y1]
+ (-0.5) [Z0 X1]
+ (-0.5) [Z0 Y1]

This Hamiltonian can then be used in defining VQE problems using :class:~.ExpvalCost.
"""
n = int(np.log2(len(H)))
N = 2 ** n

if H.shape != (N, N):
raise ValueError(
"The Hamiltonian should have shape (2**n, 2**n), for any qubit number n>=1"
)

if not np.allclose(H, H.conj().T):
raise ValueError("The Hamiltonian is not Hermitian")

paulis = [qml.Identity, qml.PauliX, qml.PauliY, qml.PauliZ]
obs = []
coeffs = []

for term in itertools.product(paulis, repeat=n):
matrices = [i._matrix() for i in term]
coeff = np.trace(functools.reduce(np.kron, matrices) @ H) / N
coeff = np.real_if_close(coeff).item()

if not np.allclose(coeff, 0):
coeffs.append(coeff)

if not all(t is qml.Identity for t in term) and hide_identity:
obs.append(
functools.reduce(
matmul,
[t(i) for i, t in enumerate(term) if t is not qml.Identity],
)
)
else:
obs.append(functools.reduce(matmul, [t(i) for i, t in enumerate(term)]))

return coeffs, obs

[docs]def sparse_hamiltonian(H):
r"""Computes the sparse matrix representation a Hamiltonian in the computational basis.

Args:
H (~.Hamiltonian): Hamiltonian operator for which the matrix representation should be
computed

Returns:
coo_matrix: a sparse matrix in scipy coordinate list (COO) format with dimension
:math:(2^n, 2^n), where :math:n is the number of wires

**Example:**

This function can be used by passing a qml.Hamiltonian object as:

>>> coeffs = [1, -0.45]
>>> obs = [qml.PauliZ(0) @ qml.PauliZ(1), qml.PauliY(0) @ qml.PauliZ(1)]
>>> H = qml.Hamiltonian(coeffs, obs)
>>> H_sparse = sparse_hamiltonian(H)
>>> H_sparse
<4x4 sparse matrix of type '<class 'numpy.complex128'>'
with 2 stored elements in COOrdinate format>

The resulting sparse matrix can be either used directly or transformed into a numpy array:

>>> H_sparse.toarray()
array([[ 1.+0.j  ,  0.+0.j  ,  0.+0.45j,  0.+0.j  ],
[ 0.+0.j  , -1.+0.j  ,  0.+0.j  ,  0.-0.45j],
[ 0.-0.45j,  0.+0.j  , -1.+0.j  ,  0.+0.j  ],
[ 0.+0.j  ,  0.+0.45j,  0.+0.j  ,  1.+0.j  ]])
"""
if not isinstance(H, qml.Hamiltonian):
raise TypeError("Passed Hamiltonian must be of type qml.Hamiltonian")

n = len(H.wires)
matrix = scipy.sparse.coo_matrix((2 ** n, 2 ** n), dtype="complex128")

for coeffs, ops in zip(H.coeffs, H.ops):

obs = [scipy.sparse.coo_matrix(o.matrix) for o in qml.operation.Tensor(ops).obs]
mat = [scipy.sparse.eye(2, format="coo")] * n

for i, j in enumerate(ops.wires):
mat[j] = obs[i]

matrix += functools.reduce(lambda i, j: scipy.sparse.kron(i, j, format="coo"), mat) * coeffs

return matrix.tocoo()

def _flatten(x):
"""Iterate recursively through an arbitrarily nested structure in depth-first order.

See also :func:_unflatten.

Args:
x (array, Iterable, Any): each element of an array or an Iterable may itself be any of these types

Yields:
Any: elements of x in depth-first order
"""
if isinstance(x, np.ndarray):
yield from _flatten(x.flat)  # should we allow object arrays? or just "yield from x.flat"?
elif isinstance(x, qml.wires.Wires):
# Reursive calls to flatten Wires will cause infinite recursion (Wires atoms are Wires).
# Since Wires are always flat, just yield.
for item in x:
yield item
elif isinstance(x, Iterable) and not isinstance(x, (str, bytes)):
for item in x:
yield from _flatten(item)
else:
yield x

def _unflatten(flat, model):
"""Restores an arbitrary nested structure to a flattened iterable.

See also :func:_flatten.

Args:
flat (array): 1D array of items
model (array, Iterable, Number): model nested structure

Raises:
TypeError: if model contains an object of unsupported type

Returns:
Union[array, list, Any], array: first elements of flat arranged into the nested
structure of model, unused elements of flat
"""
if isinstance(model, (numbers.Number, str)):
return flat[0], flat[1:]

if isinstance(model, np.ndarray):
idx = model.size
res = np.array(flat)[:idx].reshape(model.shape)
return res, flat[idx:]

if isinstance(model, Iterable):
res = []
for x in model:
val, flat = _unflatten(flat, x)
res.append(val)
return res, flat

raise TypeError("Unsupported type in the model: {}".format(type(model)))

[docs]def unflatten(flat, model):
"""Wrapper for :func:_unflatten.

Args:
flat (array): 1D array of items
model (array, Iterable, Number): model nested structure

Raises:
ValueError: if flat has more elements than model
"""
# pylint:disable=len-as-condition
res, tail = _unflatten(np.asarray(flat), model)
if len(tail) != 0:
raise ValueError("Flattened iterable has more elements than the model.")
return res

def _inv_dict(d):
"""Reverse a dictionary mapping.

Returns multimap where the keys are the former values,
and values are sets of the former keys.

Args:
d (dict[a->b]): mapping to reverse

Returns:
dict[b->set[a]]: reversed mapping
"""
ret = {}
for k, v in d.items():
return ret

def _get_default_args(func):
"""Get the default arguments of a function.

Args:
func (callable): a function

Returns:
dict[str, tuple]: mapping from argument name to (positional idx, default value)
"""
signature = inspect.signature(func)
return {
k: (idx, v.default)
for idx, (k, v) in enumerate(signature.parameters.items())
if v.default is not inspect.Parameter.empty
}

[docs]@functools.lru_cache()
def pauli_eigs(n):
r"""Eigenvalues for :math:A^{\otimes n}, where :math:A is
Pauli operator, or shares its eigenvalues.

As an example if n==2, then the eigenvalues of a tensor product consisting
of two matrices sharing the eigenvalues with Pauli matrices is returned.

Args:
n (int): the number of qubits the matrix acts on
Returns:
list: the eigenvalues of the specified observable
"""
if n == 1:
return np.array([1, -1])
return np.concatenate([pauli_eigs(n - 1), -pauli_eigs(n - 1)])

[docs]def inv(operation_list):
"""Invert a list of operations or a :doc:template </introduction/templates>.

If the inversion happens inside a QNode, the operations are removed and requeued
in the reversed order for proper inversion.

.. warning::
Use of :func:~.inv() is deprecated and should be replaced with
:func:~.adjoint().

**Example:**

The following example illuminates the inversion of a template:

.. code-block:: python3

@qml.template
def ansatz(weights, wires):
for idx, wire in enumerate(wires):
qml.RX(weights[idx], wires=[wire])

for idx in range(len(wires) - 1):
qml.CNOT(wires=[wires[idx], wires[idx + 1]])

dev = qml.device('default.qubit', wires=2)

@qml.qnode(dev)
def circuit(weights):
qml.inv(ansatz(weights, wires=[0, 1]))
return qml.expval(qml.PauliZ(0) @ qml.PauliZ(1))

We may also invert an operation sequence:

.. code-block:: python3

dev = qml.device('default.qubit', wires=2)

@qml.qnode(dev)
def circuit1():
qml.T(wires=[0]).inv()
qml.S(wires=[0]).inv()
return qml.expval(qml.PauliZ(0) @ qml.PauliZ(1))

@qml.qnode(dev)
def circuit2():
return qml.expval(qml.PauliZ(0) @ qml.PauliZ(1))

Double checking that both circuits produce the same output:

>>> ZZ1 = circuit1()
>>> ZZ2 = circuit2()
>>> assert ZZ1 == ZZ2
True

Args:
operation_list (Iterable[~.Operation]): An iterable of operations

Returns:
List[~.Operation]: The inverted list of operations
"""

warnings.warn(
"Use of qml.inv() is deprecated and should be replaced with qml.adjoint().",
UserWarning,
)
if isinstance(operation_list, qml.operation.Operation):
operation_list = [operation_list]
elif operation_list is None:
raise ValueError(
"None was passed as an argument to inv. "
"This could happen if inversion of a template without the template decorator is attempted."
)
elif callable(operation_list):
raise ValueError(
"A function was passed as an argument to inv. "
"This could happen if inversion of a template function is attempted. "
"Please use inv on the function including its arguments, as in inv(template(args))."
)
elif isinstance(operation_list, qml.tape.QuantumTape):
return new_tape

elif not isinstance(operation_list, Iterable):
raise ValueError("The provided operation_list is not iterable.")

non_ops = [
(idx, op)
for idx, op in enumerate(operation_list)
if not isinstance(op, qml.operation.Operation)
]

if non_ops:
string_reps = [" operation_list[{}] = {}".format(idx, op) for idx, op in non_ops]
raise ValueError(
"The given operation_list does not only contain Operations."
+ "The following elements of the iterable were not Operations:"
+ ",".join(string_reps)
)

for op in operation_list:
try:
# remove the queued operation to be inverted
# from the existing queuing context
qml.QueuingContext.remove(op)
except KeyError:
# operation to be inverted does not
# exist on the queuing context
pass

def qfunc():
for o in operation_list:
o.queue()

with qml.tape.QuantumTape() as tape:

return tape

[docs]def expand(matrix, original_wires, expanded_wires):
r"""Expand a an operator matrix to more wires.

Args:
matrix (array): :math:2^n \times 2^n matrix where n = len(original_wires).
original_wires (Sequence[int]): original wires of matrix
expanded_wires (Union[Sequence[int], int]): expanded wires of matrix, can be shuffled.
If a single int m is given, corresponds to list(range(m))

Returns:
array: :math:2^m \times 2^m matrix where m = len(expanded_wires).
"""
if isinstance(expanded_wires, numbers.Integral):
expanded_wires = list(range(expanded_wires))

N = len(original_wires)
M = len(expanded_wires)
D = M - N

if not set(expanded_wires).issuperset(original_wires):
raise ValueError("Invalid target subsystems provided in 'original_wires' argument.")

if matrix.shape != (2 ** N, 2 ** N):
raise ValueError(
"Matrix parameter must be of size (2**len(original_wires), 2**len(original_wires))"
)

dims = [2] * (2 * N)
tensor = matrix.reshape(dims)

if D > 0:
extra_dims = [2] * (2 * D)
identity = np.eye(2 ** D).reshape(extra_dims)
expanded_tensor = np.tensordot(tensor, identity, axes=0)
# Fix order of tensor factors
expanded_tensor = np.moveaxis(expanded_tensor, range(2 * N, 2 * N + D), range(N, N + D))
else:
expanded_tensor = tensor

wire_indices = []
for wire in original_wires:
wire_indices.append(expanded_wires.index(wire))

wire_indices = np.array(wire_indices)

# Order tensor factors according to wires
original_indices = np.array(range(N))
expanded_tensor = np.moveaxis(expanded_tensor, original_indices, wire_indices)
expanded_tensor = np.moveaxis(expanded_tensor, original_indices + M, wire_indices + M)

return expanded_tensor.reshape((2 ** M, 2 ** M))

[docs]def expand_vector(vector, original_wires, expanded_wires):
r"""Expand a vector to more wires.

Args:
vector (array): :math:2^n vector where n = len(original_wires).
original_wires (Sequence[int]): original wires of vector
expanded_wires (Union[Sequence[int], int]): expanded wires of vector, can be shuffled
If a single int m is given, corresponds to list(range(m))

Returns:
array: :math:2^m vector where m = len(expanded_wires).
"""
if isinstance(expanded_wires, numbers.Integral):
expanded_wires = list(range(expanded_wires))

N = len(original_wires)
M = len(expanded_wires)
D = M - N

if not set(expanded_wires).issuperset(original_wires):
raise ValueError("Invalid target subsystems provided in 'original_wires' argument.")

if vector.shape != (2 ** N,):
raise ValueError("Vector parameter must be of length 2**len(original_wires)")

dims = [2] * N
tensor = vector.reshape(dims)

if D > 0:
extra_dims = [2] * D
ones = np.ones(2 ** D).reshape(extra_dims)
expanded_tensor = np.tensordot(tensor, ones, axes=0)
else:
expanded_tensor = tensor

wire_indices = []
for wire in original_wires:
wire_indices.append(expanded_wires.index(wire))

wire_indices = np.array(wire_indices)

# Order tensor factors according to wires
original_indices = np.array(range(N))
expanded_tensor = np.moveaxis(expanded_tensor, original_indices, wire_indices)

return expanded_tensor.reshape(2 ** M)

[docs]def frobenius_inner_product(A, B, normalize=False):
r"""Frobenius inner product between two matrices.

.. math::

\langle A, B \rangle_F = \sum_{i,j=1}^n A_{ij} B_{ij} = \operatorname{tr} (A^T B)

The Frobenius inner product is equivalent to the Hilbert-Schmidt inner product for
matrices with real-valued entries.

Args:
A (array[float]): First matrix, assumed to be a square array.
B (array[float]): Second matrix, assumed to be a square array.
normalize (bool): If True, divide the inner_product by the Frobenius norms of A and B.
Defaults to False.

Returns:
float: Frobenius inner product of A and B

**Example**

>>> A = np.random.random((3,3))
>>> B = np.random.random((3,3))
>>> qml.utils.frobenius_inner_product(A, B)
3.091948202943376
"""
inner_product = np.sum(A * B)

if normalize:
inner_product /= np.linalg.norm(A, "fro") * np.linalg.norm(B, "fro")

return inner_product