Source code for pennylane.pauli.pauli_arithmetic
# Copyright 2018-2022 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.
"""The Pauli arithmetic abstract reduced representation classes"""
# pylint:disable=protected-access
from copy import copy
from functools import lru_cache, reduce
import numpy as np
from scipy import sparse
import pennylane as qml
from pennylane import math
from pennylane.ops import Identity, PauliX, PauliY, PauliZ, Prod, SProd, Sum
from pennylane.typing import TensorLike
from pennylane.wires import Wires, WiresLike
I = "I"
X = "X"
Y = "Y"
Z = "Z"
op_map = {
I: Identity,
X: PauliX,
Y: PauliY,
Z: PauliZ,
}
op_to_str_map = {
Identity: I,
PauliX: X,
PauliY: Y,
PauliZ: Z,
}
matI = np.eye(2)
matX = np.array([[0, 1], [1, 0]])
matY = np.array([[0, -1j], [1j, 0]])
matZ = np.array([[1, 0], [0, -1]])
mat_map = {
I: matI,
X: matX,
Y: matY,
Z: matZ,
}
anticom_map = {
I: {I: 0, X: 0, Y: 0, Z: 0},
X: {I: 0, X: 0, Y: 1, Z: 1},
Y: {I: 0, X: 1, Y: 0, Z: 1},
Z: {I: 0, X: 1, Y: 1, Z: 0},
}
@lru_cache
def _make_operation(op, wire):
return op_map[op](wire)
@lru_cache
def _cached_sparse_data(op):
"""Returns the sparse data and indices of a Pauli operator."""
if op == "I":
data = np.array([1.0, 1.0], dtype=np.complex128)
indices = np.array([0, 1], dtype=np.int64)
elif op == "X":
data = np.array([1.0, 1.0], dtype=np.complex128)
indices = np.array([1, 0], dtype=np.int64)
elif op == "Y":
data = np.array([-1.0j, 1.0j], dtype=np.complex128)
indices = np.array([1, 0], dtype=np.int64)
else: # op == "Z"
data = np.array([1.0, -1.0], dtype=np.complex128)
indices = np.array([0, 1], dtype=np.int64)
return data, indices
@lru_cache(maxsize=2)
def _cached_arange(n):
"Caches `np.arange` output to speed up sparse calculations."
return np.arange(n)
pauli_to_sparse_int = {I: 0, X: 1, Y: 1, Z: 0} # (I, Z) and (X, Y) have the same sparsity
def _ps_to_sparse_index(pauli_words, wires):
"""Represent the Pauli words sparse structure in a matrix of shape n_words x n_wires."""
indices = np.zeros((len(pauli_words), len(wires)))
for i, pw in enumerate(pauli_words):
if not pw.wires:
continue
wire_indices = np.array(wires.indices(pw.wires))
indices[i, wire_indices] = [pauli_to_sparse_int[pw[w]] for w in pw.wires]
return indices
_map_I = {
I: (1, I),
X: (1, X),
Y: (1, Y),
Z: (1, Z),
}
_map_X = {
I: (1, X),
X: (1, I),
Y: (1.0j, Z),
Z: (-1.0j, Y),
}
_map_Y = {
I: (1, Y),
X: (-1.0j, Z),
Y: (1, I),
Z: (1j, X),
}
_map_Z = {
I: (1, Z),
X: (1j, Y),
Y: (-1.0j, X),
Z: (1, I),
}
mul_map = {I: _map_I, X: _map_X, Y: _map_Y, Z: _map_Z}
[docs]class PauliWord(dict):
r"""
Immutable dictionary used to represent a Pauli Word,
associating wires with their respective operators.
Can be constructed from a standard dictionary.
.. note::
An empty :class:`~.PauliWord` will be treated as the multiplicative
identity (i.e identity on all wires). Its matrix is the identity matrix
(trivially the :math:`1\times 1` one matrix when no ``wire_order`` is passed to
``PauliWord({}).to_mat()``).
**Examples**
Initializing a Pauli word:
>>> w = PauliWord({"a": 'X', 2: 'Y', 3: 'Z'})
>>> w
X(a) @ Y(2) @ Z(3)
When multiplying Pauli words together, we obtain a :class:`~PauliSentence` with the resulting ``PauliWord`` as a key and the corresponding coefficient as its value.
>>> w1 = PauliWord({0:"X", 1:"Y"})
>>> w2 = PauliWord({1:"X", 2:"Z"})
>>> w1 @ w2
-1j * Z(1) @ Z(2) @ X(0)
We can multiply scalars to Pauli words or add/subtract them, resulting in a :class:`~PauliSentence` instance.
>>> 0.5 * w1 - 1.5 * w2 + 2
0.5 * X(0) @ Y(1)
+ -1.5 * X(1) @ Z(2)
+ 2 * I
"""
# this allows scalar multiplication from left with numpy arrays np.array(0.5) * pw1
# taken from [stackexchange](https://stackoverflow.com/questions/40694380/forcing-multiplication-to-use-rmul-instead-of-numpy-array-mul-or-byp/44634634#44634634)
__array_priority__ = 1000
def __missing__(self, key):
"""If the wire is not in the Pauli word,
then no operator acts on it, so return the Identity."""
return I
def __init__(self, mapping):
"""Strip identities from PauliWord on init!"""
for wire, op in mapping.copy().items():
if op == I:
del mapping[wire]
super().__init__(mapping)
@property
def pauli_rep(self):
"""Trivial pauli_rep"""
return PauliSentence({self: 1.0})
def __reduce__(self):
"""Defines how to pickle and unpickle a PauliWord. Otherwise, un-pickling
would cause __setitem__ to be called, which is forbidden on PauliWord.
For more information, see: https://docs.python.org/3/library/pickle.html#object.__reduce__
"""
return (PauliWord, (dict(self),))
def __copy__(self):
"""Copy the PauliWord instance."""
return PauliWord(dict(self.items()))
def __deepcopy__(self, memo):
res = self.__copy__()
memo[id(self)] = res
return res
def __setitem__(self, key, item):
"""Restrict setting items after instantiation."""
raise TypeError("PauliWord object does not support assignment")
[docs] def update(self, __m, **kwargs) -> None:
"""Restrict updating PW after instantiation."""
raise TypeError("PauliWord object does not support assignment")
def __hash__(self):
return hash(frozenset(self.items()))
def _matmul(self, other):
"""Private matrix multiplication that returns (pauli_word, coeff) tuple for more lightweight processing"""
base, iterator, swapped = (
(self, other, False) if len(self) >= len(other) else (other, self, True)
)
result = copy(dict(base))
coeff = 1
for wire, term in iterator.items():
if wire in base:
factor, new_op = mul_map[term][base[wire]] if swapped else mul_map[base[wire]][term]
if new_op == I:
del result[wire]
else:
coeff *= factor
result[wire] = new_op
elif term != I:
result[wire] = term
return PauliWord(result), coeff
def __matmul__(self, other):
"""Multiply two Pauli words together using the matrix product if wires overlap
and the tensor product otherwise.
Empty Pauli words are treated as the Identity operator on all wires.
Args:
other (PauliWord): The Pauli word to multiply with
Returns:
PauliSentence: coeff * new_word
"""
if isinstance(other, PauliSentence):
return PauliSentence({self: 1.0}) @ other
new_word, coeff = self._matmul(other)
return PauliSentence({new_word: coeff})
def __mul__(self, other):
"""Multiply a PauliWord by a scalar
Args:
other (Scalar): The scalar to multiply the PauliWord with
Returns:
PauliSentence
"""
if isinstance(other, TensorLike):
if not qml.math.ndim(other) == 0:
raise ValueError(
f"Attempting to multiply a PauliWord with an array of dimension {qml.math.ndim(other)}"
)
return PauliSentence({self: other})
raise TypeError(
f"PauliWord can only be multiplied by numerical data. Attempting to multiply by {other} of type {type(other)}"
)
__rmul__ = __mul__
def __add__(self, other):
"""Add PauliWord instances and scalars to PauliWord.
Returns a PauliSentence."""
# Note that the case of PauliWord + PauliSentence is covered in PauliSentence
if isinstance(other, PauliWord):
if other == self:
return PauliSentence({self: 2.0})
return PauliSentence({self: 1.0, other: 1.0})
if isinstance(other, TensorLike):
# Scalars are interepreted as scalar * Identity
IdWord = PauliWord({})
if IdWord == self:
return PauliSentence({self: 1.0 + other})
return PauliSentence({self: 1.0, IdWord: other})
return NotImplemented
__radd__ = __add__
def __iadd__(self, other):
"""Inplace addition"""
return self + other
def __sub__(self, other):
"""Subtract other PauliSentence, PauliWord, or scalar"""
return self + -1 * other
def __rsub__(self, other):
"""Subtract other PauliSentence, PauliWord, or scalar"""
return -1 * self + other
def __truediv__(self, other):
"""Divide a PauliWord by a scalar"""
if isinstance(other, TensorLike):
return self * (1 / other)
raise TypeError(
f"PauliWord can only be divided by numerical data. Attempting to divide by {other} of type {type(other)}"
)
[docs] def commutes_with(self, other):
"""Fast check if two PauliWords commute with each other"""
wires = set(self) & set(other)
if not wires:
return True
anticom_count = sum(anticom_map[self[wire]][other[wire]] for wire in wires)
return (anticom_count % 2) == 0
def _commutator(self, other):
"""comm between two PauliWords, returns tuple (new_word, coeff) for faster arithmetic"""
# This may be helpful to developers that need a more lightweight comm between pauli words
# without creating PauliSentence classes
if self.commutes_with(other):
return PauliWord({}), 0.0
new_word, coeff = self._matmul(other)
return new_word, 2 * coeff
[docs] def commutator(self, other):
"""
Compute commutator between a ``PauliWord`` :math:`P` and other operator :math:`O`
.. math:: [P, O] = P O - O P
When the other operator is a :class:`~PauliWord` or :class:`~PauliSentence`,
this method is faster than computing ``P @ O - O @ P``. It is what is being used
in :func:`~commutator` when setting ``pauli=True``.
Args:
other (Union[Operator, PauliWord, PauliSentence]): Second operator
Returns:
~PauliSentence: The commutator result in form of a :class:`~PauliSentence` instances.
**Examples**
You can compute commutators between :class:`~PauliWord` instances.
>>> pw = PauliWord({0:"X"})
>>> pw.commutator(PauliWord({0:"Y"}))
2j * Z(0)
You can also compute the commutator with other operator types if they have a Pauli representation.
>>> pw.commutator(qml.Y(0))
2j * Z(0)
"""
if isinstance(other, PauliWord):
new_word, coeff = self._commutator(other)
if coeff == 0:
return PauliSentence({})
return PauliSentence({new_word: coeff})
if isinstance(other, qml.operation.Operator):
op_self = PauliSentence({self: 1.0})
return op_self.commutator(other)
if isinstance(other, PauliSentence):
# for infix method, this would be handled by __ror__
return -1.0 * other.commutator(self)
raise NotImplementedError(
f"Cannot compute natively a commutator between PauliWord and {other} of type {type(other)}"
)
def __str__(self):
"""String representation of a PauliWord."""
if len(self) == 0:
return "I"
return " @ ".join(f"{op}({w})" for w, op in self.items())
def __repr__(self):
"""Terminal representation for PauliWord"""
return str(self)
@property
def wires(self):
"""Track wires in a PauliWord."""
return Wires(self)
[docs] def to_mat(self, wire_order=None, format="dense", coeff=1.0):
"""Returns the matrix representation.
Keyword Args:
wire_order (iterable or None): The order of qubits in the tensor product.
format (str): The format of the matrix. It is "dense" by default. Use "csr" for sparse.
coeff (float): Coefficient multiplying the resulting matrix.
Returns:
(Union[NumpyArray, ScipySparseArray]): Matrix representation of the Pauli word.
Raises:
ValueError: Can't get the matrix of an empty PauliWord.
"""
wire_order = self.wires if wire_order is None else Wires(wire_order)
if not wire_order.contains_wires(self.wires):
raise ValueError(
"Can't get the matrix for the specified wire order because it "
f"does not contain all the Pauli word's wires {self.wires}"
)
if len(self) == 0:
n = len(wire_order) if wire_order is not None else 0
return (
np.diag([coeff] * 2**n)
if format == "dense"
else coeff * sparse.eye(2**n, format=format, dtype="complex128")
)
if format == "dense":
return coeff * reduce(math.kron, (mat_map[self[w]] for w in wire_order))
return self._to_sparse_mat(wire_order, coeff)
def _to_sparse_mat(self, wire_order, coeff):
"""Compute the sparse matrix of the Pauli word times a coefficient, given a wire order.
See pauli_sparse_matrices.md for the technical details of the implementation."""
matrix_size = 2 ** len(wire_order)
matrix = sparse.csr_matrix((matrix_size, matrix_size), dtype="complex128")
# Avoid checks and copies in __init__ by directly setting the attributes of an empty matrix
matrix.data = self._get_csr_data(wire_order, coeff)
matrix.indices = self._get_csr_indices(wire_order)
matrix.indptr = _cached_arange(matrix_size + 1) # Non-zero entries by row (starting from 0)
return matrix
def _get_csr_data(self, wire_order, coeff):
"""Computes the sparse matrix data of the Pauli word times a coefficient, given a wire order."""
full_word = [self[wire] for wire in wire_order]
matrix_size = 2 ** len(wire_order)
if len(self) == 0:
return np.full(matrix_size, coeff, dtype=np.complex128)
data = np.empty(matrix_size, dtype=np.complex128) # Non-zero values
current_size = 2
data[:current_size], _ = _cached_sparse_data(full_word[-1])
data[:current_size] *= coeff # Multiply initial term better than the full matrix
for s in full_word[-2::-1]:
if s == "I":
data[current_size : 2 * current_size] = data[:current_size]
elif s == "X":
data[current_size : 2 * current_size] = data[:current_size]
elif s == "Y":
data[current_size : 2 * current_size] = 1j * data[:current_size]
data[:current_size] *= -1j
elif s == "Z":
data[current_size : 2 * current_size] = -data[:current_size]
current_size *= 2
return data
def _get_csr_data_2(self, wire_order, coeff):
"""Computes the sparse matrix data of the Pauli word times a coefficient, given a wire order."""
full_word = [self[wire] for wire in wire_order]
nwords = len(full_word)
if nwords < 2:
return np.array([1.0]), self._get_csr_data(wire_order, coeff)
outer = self._get_csr_data(wire_order[: nwords // 2], 1.0)
inner = self._get_csr_data(wire_order[nwords // 2 :], coeff)
return outer, inner
def _get_csr_indices(self, wire_order):
"""Computes the sparse matrix indices of the Pauli word times a coefficient, given a wire order."""
full_word = [self[wire] for wire in wire_order]
matrix_size = 2 ** len(wire_order)
if len(self) == 0:
return _cached_arange(matrix_size)
indices = np.empty(matrix_size, dtype=np.int64) # Column index of non-zero values
current_size = 2
_, indices[:current_size] = _cached_sparse_data(full_word[-1])
for s in full_word[-2::-1]:
if s == "I":
indices[current_size : 2 * current_size] = indices[:current_size] + current_size
elif s == "X":
indices[current_size : 2 * current_size] = indices[:current_size]
indices[:current_size] += current_size
elif s == "Y":
indices[current_size : 2 * current_size] = indices[:current_size]
indices[:current_size] += current_size
elif s == "Z":
indices[current_size : 2 * current_size] = indices[:current_size] + current_size
current_size *= 2
return indices
[docs] def operation(self, wire_order: WiresLike = ()):
"""Returns a native PennyLane :class:`~pennylane.operation.Operation` representing the PauliWord."""
if len(self) == 0:
return Identity(wires=wire_order)
if qml.capture.enabled():
# cant use lru_cache with program capture
factors = [op_map[op](wire) for wire, op in self.items()]
else:
factors = [_make_operation(op, wire) for wire, op in self.items()]
pauli_rep = PauliSentence({self: 1})
return factors[0] if len(factors) == 1 else Prod(*factors, _pauli_rep=pauli_rep)
[docs] def map_wires(self, wire_map: dict) -> "PauliWord":
"""Return a new PauliWord with the wires mapped."""
return self.__class__({wire_map.get(w, w): op for w, op in self.items()})
pw_id = PauliWord({}) # empty pauli word to be re-used
[docs]class PauliSentence(dict):
r"""Dictionary representing a linear combination of Pauli words, with the keys
as :class:`~pennylane.pauli.PauliWord` instances and the values correspond to coefficients.
.. note::
An empty :class:`~.PauliSentence` will be treated as the additive
identity (i.e ``0 * Identity()``). Its matrix is the all-zero matrix
(trivially the :math:`1\times 1` zero matrix when no ``wire_order`` is passed to
``PauliSentence({}).to_mat()``).
**Examples**
>>> from pennylane import PauliSentence, PauliWord
>>> ps = PauliSentence({
... PauliWord({0:'X', 1:'Y'}): 1.23,
... PauliWord({2:'Z', 0:'Y'}): -0.45j
... })
>>> ps
1.23 * X(0) @ Y(1)
+ (-0-0.45j) * Z(2) @ Y(0)
Combining Pauli words automatically results in Pauli sentences that can be used to construct more complicated operators.
>>> w1 = PauliWord({0:"X", 1:"Y"})
>>> w2 = PauliWord({1:"X", 2:"Z"})
>>> ps = 0.5 * w1 - 1.5 * w2 + 2
>>> ps + PauliWord({3:"Z"}) - 1
0.5 * X(0) @ Y(1)
+ -1.5 * X(1) @ Z(2)
+ 1 * I
+ 1.0 * Z(3)
Note that while the empty :class:`~PauliWord` ``PauliWord({})`` respresents the identity, the empty ``PauliSentence`` represents 0
>>> PauliSentence({})
0 * I
We can compute commutators using the ``PauliSentence.commutator()`` method
>>> op1 = PauliWord({0:"X", 1:"X"})
>>> op2 = PauliWord({0:"Y"}) + PauliWord({1:"Y"})
>>> op1.commutator(op2)
2j * Z(0) @ X(1)
+ 2j * X(0) @ Z(1)
Or, alternatively, use :func:`~commutator`.
>>> qml.commutator(op1, op2, pauli=True)
Note that we need to specify ``pauli=True`` as :func:`~.commutator` returns PennyLane operators by default.
"""
# this allows scalar multiplication from left with numpy arrays np.array(0.5) * ps1
# taken from [stackexchange](https://stackoverflow.com/questions/40694380/forcing-multiplication-to-use-rmul-instead-of-numpy-array-mul-or-byp/44634634#44634634)
__array_priority__ = 1000
@property
def pauli_rep(self):
"""Trivial pauli_rep"""
return self
def __missing__(self, key):
"""If the PauliWord is not in the sentence then the coefficient
associated with it should be 0."""
return 0.0
[docs] def trace(self):
r"""Return the normalized trace of the ``PauliSentence`` instance
.. math:: \frac{1}{2^n} \text{tr}\left( P \right).
The normalized trace does not scale with the number of qubits :math:`n`.
>>> PauliSentence({PauliWord({0:"I", 1:"I"}): 0.5}).trace()
0.5
>>> PauliSentence({PauliWord({}): 0.5}).trace()
0.5
"""
return self.get(pw_id, 0.0)
def __add__(self, other):
"""Add a PauliWord, scalar or other PauliSentence to a PauliSentence.
Empty Pauli sentences are treated as the additive identity
(i.e 0 * Identity on all wires). The non-empty Pauli sentence is returned.
"""
if isinstance(other, PauliSentence):
smaller_ps, larger_ps = (
(self, copy(other)) if len(self) < len(other) else (other, copy(self))
)
for key in smaller_ps:
larger_ps[key] += smaller_ps[key]
return larger_ps
if isinstance(other, PauliWord):
res = copy(self)
if other in res:
res[other] += 1.0
else:
res[other] = 1.0
return res
if isinstance(other, TensorLike):
# Scalars are interepreted as scalar * Identity
res = copy(self)
IdWord = PauliWord({})
if IdWord in res:
res[IdWord] += other
else:
res[IdWord] = other
return res
raise TypeError(f"Cannot add {other} of type {type(other)} to PauliSentence")
__radd__ = __add__
def __iadd__(self, other):
"""Inplace addition of two Pauli sentence together by adding terms of other to self"""
if isinstance(other, PauliSentence):
for key in other:
if key in self:
self[key] += other[key]
else:
self[key] = other[key]
return self
if isinstance(other, PauliWord):
if other in self:
self[other] += 1.0
else:
self[other] = 1.0
return self
if isinstance(other, TensorLike):
IdWord = PauliWord({})
if IdWord in self:
self[IdWord] += other
else:
self[IdWord] = other
return self
raise TypeError(f"Cannot add {other} of type {type(other)} to PauliSentence")
def __sub__(self, other):
"""Subtract other PauliSentence, PauliWord, or scalar"""
return self + -1 * other
def __rsub__(self, other):
"""Subtract other PauliSentence, PauliWord, or scalar"""
return -1 * self + other
def __copy__(self):
"""Copy the PauliSentence instance."""
copied_ps = {}
for pw, coeff in self.items():
copied_ps[copy(pw)] = coeff
return PauliSentence(copied_ps)
def __deepcopy__(self, memo):
res = self.__copy__()
memo[id(self)] = res
return res
def __matmul__(self, other):
"""Matrix / tensor product between two PauliSentences by iterating over each sentence and multiplying
the Pauli words pair-wise"""
if isinstance(other, PauliWord):
other = PauliSentence({other: 1.0})
final_ps = PauliSentence()
if len(self) == 0 or len(other) == 0:
return final_ps
for pw1 in self:
for pw2 in other:
prod_pw, coeff = pw1._matmul(pw2)
final_ps[prod_pw] = final_ps[prod_pw] + coeff * self[pw1] * other[pw2]
return final_ps
def __mul__(self, other):
"""Multiply a PauliWord by a scalar
Args:
other (Scalar): The scalar to multiply the PauliWord with
Returns:
PauliSentence
"""
if isinstance(other, TensorLike):
if not qml.math.ndim(other) == 0:
raise ValueError(
f"Attempting to multiply a PauliSentence with an array of dimension {qml.math.ndim(other)}"
)
return PauliSentence({key: other * value for key, value in self.items()})
raise TypeError(
f"PauliSentence can only be multiplied by numerical data. Attempting to multiply by {other} of type {type(other)}"
)
__rmul__ = __mul__
def __truediv__(self, other):
"""Divide a PauliSentence by a scalar"""
if isinstance(other, TensorLike):
return self * (1 / other)
raise TypeError(
f"PauliSentence can only be divided by numerical data. Attempting to divide by {other} of type {type(other)}"
)
[docs] def commutator(self, other):
"""
Compute commutator between a ``PauliSentence`` :math:`P` and other operator :math:`O`
.. math:: [P, O] = P O - O P
When the other operator is a :class:`~PauliWord` or :class:`~PauliSentence`,
this method is faster than computing ``P @ O - O @ P``. It is what is being used
in :func:`~commutator` when setting ``pauli=True``.
Args:
other (Union[Operator, PauliWord, PauliSentence]): Second operator
Returns:
~PauliSentence: The commutator result in form of a :class:`~PauliSentence` instances.
**Examples**
You can compute commutators between :class:`~PauliSentence` instances.
>>> pw1 = PauliWord({0:"X"})
>>> pw2 = PauliWord({1:"X"})
>>> ps1 = PauliSentence({pw1: 1., pw2: 2.})
>>> ps2 = PauliSentence({pw1: 0.5j, pw2: 1j})
>>> ps1.commutator(ps2)
0 * I
You can also compute the commutator with other operator types if they have a Pauli representation.
>>> ps1.commutator(qml.Y(0))
2j * Z(0)"""
final_ps = PauliSentence()
if isinstance(other, PauliWord):
for pw1 in self:
comm_pw, coeff = pw1._commutator(other)
if len(comm_pw) != 0:
final_ps[comm_pw] += coeff * self[pw1]
return final_ps
if not isinstance(other, PauliSentence):
if other.pauli_rep is None:
raise NotImplementedError(
f"Cannot compute a native commutator of a Pauli word or sentence with the operator {other} of type {type(other)}."
f"You can try to use qml.commutator(op1, op2, pauli=False) instead."
)
other = qml.pauli.pauli_sentence(other)
for pw1 in self:
for pw2 in other:
comm_pw, coeff = pw1._commutator(pw2)
if len(comm_pw) != 0:
final_ps[comm_pw] += coeff * self[pw1] * other[pw2]
return final_ps
def __str__(self):
"""String representation of the PauliSentence."""
if len(self) == 0:
return "0 * I"
return "\n+ ".join(f"{coeff} * {str(pw)}" for pw, coeff in self.items())
def __repr__(self):
"""Terminal representation for PauliSentence"""
return str(self)
@property
def wires(self):
"""Track wires of the PauliSentence."""
return Wires.all_wires((pw.wires for pw in self.keys()))
[docs] def to_mat(self, wire_order=None, format="dense", buffer_size=None):
"""Returns the matrix representation.
Keyword Args:
wire_order (iterable or None): The order of qubits in the tensor product.
format (str): The format of the matrix. It is "dense" by default. Use "csr" for sparse.
buffer_size (int or None): The maximum allowed memory in bytes to store intermediate results
in the calculation of sparse matrices. It defaults to ``2 ** 30`` bytes that make
1GB of memory. In general, larger buffers allow faster computations.
Returns:
(Union[NumpyArray, ScipySparseArray]): Matrix representation of the Pauli sentence.
Raises:
ValueError: Can't get the matrix of an empty PauliSentence.
"""
wire_order = self.wires if wire_order is None else Wires(wire_order)
if len(self) == 0:
n = len(wire_order) if wire_order is not None else 0
if format == "dense":
return np.zeros((2**n, 2**n))
return sparse.csr_matrix((2**n, 2**n), dtype="complex128").asformat(format)
if format == "dense":
return self._to_dense_mat(wire_order)
return self._to_sparse_mat(wire_order, buffer_size=buffer_size).asformat(format)
def _to_sparse_mat(self, wire_order, buffer_size=None):
"""Compute the sparse matrix of the Pauli sentence by efficiently adding the Pauli words
that it is composed of. See pauli_sparse_matrices.md for the technical details."""
pauli_words = list(self) # Ensure consistent ordering
n_wires = len(wire_order)
matrix_size = 2**n_wires
matrix = sparse.csr_matrix((matrix_size, matrix_size), dtype="complex128")
op_sparse_idx = _ps_to_sparse_index(pauli_words, wire_order)
_, unique_sparse_structures, unique_invs = np.unique(
op_sparse_idx, axis=0, return_index=True, return_inverse=True
)
pw_sparse_structures = unique_sparse_structures[unique_invs]
buffer_size = buffer_size or 2**30 # Default to 1GB of memory
# Convert bytes to number of matrices:
# complex128 (16) for each data entry and int64 (8) for each indices entry
buffer_size = max(1, buffer_size // ((16 + 8) * matrix_size))
mat_data = np.empty((matrix_size, buffer_size), dtype=np.complex128)
mat_indices = np.empty((matrix_size, buffer_size), dtype=np.int64)
n_matrices_in_buffer = 0
for sparse_structure in unique_sparse_structures:
indices, *_ = np.nonzero(pw_sparse_structures == sparse_structure)
mat = self._sum_same_structure_pws([pauli_words[i] for i in indices], wire_order)
mat_data[:, n_matrices_in_buffer] = mat.data
mat_indices[:, n_matrices_in_buffer] = mat.indices
n_matrices_in_buffer += 1
if n_matrices_in_buffer == buffer_size:
# Add partial results in batches to control the memory usage
matrix += self._sum_different_structure_pws(mat_indices, mat_data)
n_matrices_in_buffer = 0
matrix += self._sum_different_structure_pws(
mat_indices[:, :n_matrices_in_buffer], mat_data[:, :n_matrices_in_buffer]
)
matrix.eliminate_zeros()
return matrix
def _to_dense_mat(self, wire_order):
"""Compute the dense matrix of the Pauli sentence by efficiently adding the Pauli words
that it is composed of. See pauli_sparse_matrices.md for the technical details."""
pauli_words = list(self) # Ensure consistent ordering
try:
op_sparse_idx = _ps_to_sparse_index(pauli_words, wire_order)
except qml.wires.WireError as e:
raise ValueError(
"Can't get the matrix for the specified wire order because it "
f"does not contain all the Pauli sentence's wires {self.wires}"
) from e
_, unique_sparse_structures, unique_invs = np.unique(
op_sparse_idx, axis=0, return_index=True, return_inverse=True
)
pw_sparse_structures = unique_sparse_structures[unique_invs]
full_matrix = None
for sparse_structure in unique_sparse_structures:
indices, *_ = np.nonzero(pw_sparse_structures == sparse_structure)
mat = self._sum_same_structure_pws_dense([pauli_words[i] for i in indices], wire_order)
full_matrix = mat if full_matrix is None else qml.math.add(full_matrix, mat)
return full_matrix
[docs] def dot(self, vector, wire_order=None):
"""Computes the matrix-vector product of the Pauli sentence with a state vector.
See pauli_sparse_matrices.md for the technical details."""
wire_order = self.wires if wire_order is None else Wires(wire_order)
if not wire_order.contains_wires(self.wires):
raise ValueError(
"Can't get the matrix for the specified wire order because it "
f"does not contain all the Pauli sentence's wires {self.wires}"
)
pauli_words = list(self) # Ensure consistent ordering
op_sparse_idx = _ps_to_sparse_index(pauli_words, wire_order)
_, unique_sparse_structures, unique_invs = np.unique(
op_sparse_idx, axis=0, return_index=True, return_inverse=True
)
pw_sparse_structures = unique_sparse_structures[unique_invs]
dtype = np.complex64 if vector.dtype in (np.float32, np.complex64) else np.complex128
if vector.ndim == 1:
vector = vector.reshape(1, -1)
mv = np.zeros_like(vector, dtype=dtype)
for sparse_structure in unique_sparse_structures:
indices, *_ = np.nonzero(pw_sparse_structures == sparse_structure)
entries, data = self._get_same_structure_csr(
[pauli_words[i] for i in indices], wire_order
)
mv += vector[:, entries] * data.reshape(1, -1)
return mv.reshape(vector.shape)
def _get_same_structure_csr(self, pauli_words, wire_order):
"""Returns the CSR indices and data for Pauli words with the same sparse structure."""
indices = pauli_words[0]._get_csr_indices(wire_order)
nwires = len(wire_order)
nwords = len(pauli_words)
inner = np.empty((nwords, 2 ** (nwires - nwires // 2)), dtype=np.complex128)
outer = np.empty((nwords, 2 ** (nwires // 2)), dtype=np.complex128)
for i, word in enumerate(pauli_words):
outer[i, :], inner[i, :] = word._get_csr_data_2(
wire_order, coeff=qml.math.to_numpy(self[word])
)
data = outer.T @ inner
return indices, data.ravel()
def _sum_same_structure_pws_dense(self, pauli_words, wire_order):
matrix_size = 2 ** (len(wire_order))
base_matrix = sparse.csr_matrix((matrix_size, matrix_size), dtype="complex128")
data0 = pauli_words[0]._get_csr_data(wire_order, 1)
base_matrix.data = np.ones_like(data0)
base_matrix.indices = pauli_words[0]._get_csr_indices(wire_order)
base_matrix.indptr = _cached_arange(
matrix_size + 1
) # Non-zero entries by row (starting from 0)
base_matrix = base_matrix.toarray()
coeff = self[pauli_words[0]]
ml_interface = qml.math.get_interface(coeff)
if ml_interface == "torch":
data0 = qml.math.convert_like(data0, coeff)
data = coeff * data0 if qml.math.ndim(coeff) == 0 else qml.math.outer(coeff, data0)
for pw in pauli_words[1:]:
coeff = self[pw]
csr_data = pw._get_csr_data(wire_order, 1)
ml_interface = qml.math.get_interface(coeff)
if ml_interface == "torch":
csr_data = qml.math.convert_like(csr_data, coeff)
data += (
coeff * csr_data if qml.math.ndim(coeff) == 0 else qml.math.outer(coeff, csr_data)
)
return qml.math.einsum("ij,...i->...ij", base_matrix, data)
def _sum_same_structure_pws(self, pauli_words, wire_order):
"""Sums Pauli words with the same sparse structure."""
mat = pauli_words[0].to_mat(
wire_order, coeff=qml.math.to_numpy(self[pauli_words[0]]), format="csr"
)
for word in pauli_words[1:]:
mat.data += word.to_mat(
wire_order, coeff=qml.math.to_numpy(self[word]), format="csr"
).data
return mat
@staticmethod
def _sum_different_structure_pws(indices, data):
"""Sums Pauli words with different parse structures."""
size = indices.shape[0]
idx = np.argsort(indices, axis=1)
matrix = sparse.csr_matrix((size, size), dtype="complex128")
matrix.indices = np.take_along_axis(indices, idx, axis=1).ravel()
matrix.data = np.take_along_axis(data, idx, axis=1).ravel()
num_entries_per_row = indices.shape[1]
matrix.indptr = _cached_arange(size + 1) * num_entries_per_row
# remove zeros and things sufficiently close to zero
matrix.data[np.abs(matrix.data) < 1e-16] = 0 # Faster than np.isclose(matrix.data, 0)
matrix.eliminate_zeros()
return matrix
[docs] def operation(self, wire_order: WiresLike = ()):
"""Returns a native PennyLane :class:`~pennylane.operation.Operation` representing the PauliSentence."""
if len(self) == 0:
return qml.s_prod(0, Identity(wires=wire_order))
summands = []
wire_order = wire_order or self.wires
for pw, coeff in self.items():
pw_op = pw.operation(wire_order=list(wire_order))
rep = PauliSentence({pw: coeff})
summands.append(
pw_op if qml.math.all(coeff == 1) else SProd(coeff, pw_op, _pauli_rep=rep)
)
return summands[0] if len(summands) == 1 else Sum(*summands, _pauli_rep=self)
[docs] def simplify(self, tol=1e-8):
"""Remove any PauliWords in the PauliSentence with coefficients less than the threshold tolerance."""
items = list(self.items())
for pw, coeff in items:
if abs(coeff) <= tol:
del self[pw]
if len(self) == 0:
self = PauliSentence({}) # pylint: disable=self-cls-assignment
[docs] def map_wires(self, wire_map: dict) -> "PauliSentence":
"""Return a new PauliSentence with the wires mapped."""
return self.__class__({pw.map_wires(wire_map): coeff for pw, coeff in self.items()})
_modules/pennylane/pauli/pauli_arithmetic
Download Python script
Download Notebook
View on GitHub