# -*- coding: utf-8 -*-
# =============================================================================
# filter_functions
# Copyright (C) 2019 Quantum Technology Group, RWTH Aachen University
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Contact email: tobias.hangleiter@rwth-aachen.de
# =============================================================================
"""
This module defines the Basis class, a subclass of NumPy's ndarray, to
represent operator bases.
Classes
-------
:class:`Basis`
The operator basis as an array of shape (d**2, d, d) with d the
dimension of the Hilbert space
Functions
---------
:func:`normalize`
Function to normalize a ``Basis`` instance
:func:`expand`
Function to expand an array of operators in a given basis
:func:`ggm_expand`
Fast function to expand an array of operators in a Generalized
Gell-Mann basis
"""
from itertools import product
from typing import Optional, Sequence, Union
from warnings import warn
import numpy as np
import opt_einsum as oe
from numpy import linalg as nla
from numpy.core import ndarray
from scipy import linalg as sla
from sparse import COO
from . import util
__all__ = ['Basis', 'expand', 'ggm_expand', 'normalize']
[docs]class Basis(ndarray):
r"""
Class for operator bases. There are several ways to instantiate a
Basis object:
- by just calling this constructor with a (possibly incomplete)
array of basis matrices. No checks regarding orthonormality or
hermiticity are performed.
- by calling one of the classes alternative constructors
(classmethods):
- :meth:`pauli`: Pauli operator basis
- :meth:`ggm`: Generalized Gell-Mann basis
- :meth:`from_partial` Generate an complete basis from
partial elements
These bases guarantee the following properties:
- hermitian
- orthonormal
- [traceless] (can be controlled by a flag)
Since Basis is a subclass of NumPy's ``ndarray``, it inherits all of
its attributes, e.g. ``shape``. The following attributes behave
slightly differently to a ndarray, however
- ``A == B`` is ``True`` if all elements evaluate almost equal,
i.e. equivalent to ``np.allclose(A, B)``.
- ``basis.T`` transposes the last two axes of ``basis``. For a
full basis, this corresponds to transposing each element
individually. For a basis element, it corresponds to normal
transposition.
Parameters
----------
basis_array: array_like, shape (n, d, d)
An array or list of square matrices that are elements of an
operator basis spanning :math:`\mathbb{C}^{d\times d}`. *n*
should be smaller than or equal to *d**2*.
traceless: bool, optional (default: auto)
Controls whether a traceless basis is forced. Here, traceless
means that the first element of the basis is the identity and
the remaining elements are matrices of trace zero. If an element
of ``basis_array`` is neither traceless nor the identity and
``traceless == True``, an exception will be raised. Defaults to
``True`` if basis_array is traceless and ``False`` if not.
btype: str, optional (default: ``'custom'``)
A string describing the basis type. For example, a basis created
by the factory method :meth:`pauli` has *btype* 'pauli'.
labels: sequence of str, optional
A list of labels for the individual basis elements. Defaults to
'C_0', 'C_1', ...
Attributes
----------
Other than the attributes inherited from ``ndarray``, a ``Basis``
instance has the following attributes:
btype: str
Basis type.
labels: sequence of str
The labels for the basis elements.
d: int
Dimension of the space spanned by the basis.
H: Basis
Hermitian conjugate.
isherm: bool
If the basis is hermitian.
isorthonorm: bool
If the basis is orthonormal.
istraceless: bool
If the basis is traceless except for an identity element
iscomplete: bool
If the basis is complete, ie spans the full space.
sparse: COO, shape (n, d, d)
Representation in the COO format supplied by the ``sparse``
package.
four_element_traces: COO, shape (n, n, n, n)
Traces over all possible combinations of four elements of self.
This is required for the calculation of the error transfer
matrix and thus cached in the Basis instance.
Most of the attributes above are properties which are lazily
evaluated and cached.
Methods
-------
Other than the methods inherited from ``ndarray``, a ``Basis``
instance has the following methods:
normalize(b)
Normalizes the basis (used internally when creating a basis from
elements)
tidyup(eps_scale=None)
Cleans up floating point errors in-place to make zeros actual
zeros. ``eps_scale`` is an optional argument multiplied to the
data type's ``eps`` to get the absolute tolerance.
"""
def __new__(cls, basis_array: Sequence, traceless: Optional[bool] = None,
btype: Optional[str] = None, labels: Optional[Sequence[str]] = None) -> 'Basis':
"""Constructor."""
if not hasattr(basis_array, '__getitem__'):
raise TypeError('Invalid data type. Must be array_like')
if isinstance(basis_array, cls):
basis = basis_array
else:
try:
# Allow single 2d element
if len(basis_array.shape) == 2:
basis_array = [basis_array]
except AttributeError:
pass
basis = util.parse_operators(basis_array, 'basis_array')
if basis.shape[0] > np.product(basis.shape[1:]):
raise ValueError('Given overcomplete set of basis matrices. '
'Not linearly independent.')
basis = basis.view(cls)
basis.btype = btype or 'Custom'
basis.d = basis.shape[-1]
if labels is not None:
if len(labels) != len(basis):
raise ValueError(f'Got {len(labels)} basis labels but expected {len(basis)}')
basis.labels = labels
else:
basis.labels = [f'$C_{{{i}}}$' for i in range(len(basis))]
return basis
def __array_finalize__(self, basis: 'Basis') -> None:
"""Required for subclassing ndarray."""
if basis is None:
return
self.btype = getattr(basis, 'btype', 'Custom')
self.labels = getattr(basis, 'labels', [f'$C_{{{i}}}$' for i in range(len(basis))])
self.d = getattr(basis, 'd', basis.shape[-1])
self._sparse = None
self._four_element_traces = None
self._isherm = None
self._isorthonorm = None
self._istraceless = None
self._iscomplete = None
self._eps = np.finfo(complex).eps
self._atol = self._eps*self.d**3
self._rtol = 0
def __eq__(self, other: object) -> bool:
"""Compare for equality."""
try:
if self.shape != other.shape:
return False
except AttributeError:
# Not ndarray
return np.equal(self, other)
return np.allclose(self.view(ndarray), other.view(ndarray),
atol=self._atol, rtol=self._rtol)
def __contains__(self, item: ndarray) -> bool:
"""Implement 'in' operator."""
return any(np.isclose(item.view(ndarray), self.view(ndarray),
rtol=self._rtol, atol=self._atol).all(axis=(1, 2)))
def __array_wrap__(self, out_arr, context=None):
"""
Fixes problem that ufuncs return 0-d arrays instead of scalars.
https://github.com/numpy/numpy/issues/5819#issue-72454838
"""
if out_arr.ndim:
return ndarray.__array_wrap__(self, out_arr, context)
def _print_checks(self) -> None:
"""Print checks for debug purposes."""
checks = ('isherm', 'istraceless', 'iscomplete', 'isorthonorm')
for check in checks:
print(check, ':\t', getattr(self, check))
@property
def isherm(self) -> bool:
"""Returns True if all basis elements are hermitian."""
if self._isherm is None:
self._isherm = (self.H == self)
return self._isherm
@property
def isorthonorm(self) -> bool:
"""Returns True if basis is orthonormal."""
if self._isorthonorm is None:
# All the basis is orthonormal iff the matrix consisting of all
# d**2 elements written as d**2-dimensional column vectors is
# unitary.
if self.ndim == 2 or len(self) == 1:
# Only one basis element
self._isorthonorm = True
else:
# Size of the result after multiplication
dim = self.shape[0]
U = self.reshape((dim, -1))
actual = U.conj() @ U.T
target = np.identity(dim)
atol = self._eps*(self.d**2)**3
self._isorthonorm = np.allclose(actual.view(ndarray), target,
atol=atol, rtol=self._rtol)
return self._isorthonorm
@property
def istraceless(self) -> bool:
"""
Returns True if basis is traceless except for possibly the identity.
"""
if self._istraceless is None:
trace = np.einsum('...jj', self)
trace = util.remove_float_errors(trace, self.d**2)
nonzero = trace.nonzero()
if nonzero[0].size == 0:
self._istraceless = True
elif nonzero[0].size == 1:
# Single element has nonzero trace, check if (proportional to)
# identity
elem = self[nonzero][0].view(ndarray) if self.ndim == 3 else self.view(ndarray)
offdiag_nonzero = elem[~np.eye(self.d, dtype=bool)].nonzero()
diag_equal = np.diag(elem) == elem[0, 0]
if diag_equal.all() and not offdiag_nonzero[0].any():
# Element is (proportional to) the identity, this we define
# as 'traceless' since a complete basis cannot have only
# traceless elems.
self._istraceless = True
else:
# Element not the identity, therefore not traceless
self._istraceless = False
else:
self._istraceless = False
return self._istraceless
@property
def iscomplete(self) -> bool:
"""Returns True if basis is complete."""
if self._iscomplete is None:
A = self.reshape(self.shape[0], -1)
rank = np.linalg.matrix_rank(A)
self._iscomplete = rank == self.d**2
return self._iscomplete
@property
def H(self) -> 'Basis':
"""Return the basis hermitian conjugated element-wise."""
return self.T.conj()
@property
def T(self) -> 'Basis':
"""Return the basis transposed element-wise."""
if self.ndim == 3:
return self.transpose(0, 2, 1)
if self.ndim == 2:
return self.transpose(1, 0)
return self
@property
def sparse(self) -> COO:
"""Return the basis as a sparse COO array"""
if self._sparse is None:
self._sparse = COO.from_numpy(self)
return self._sparse
@property
def four_element_traces(self) -> COO:
r"""
Return all traces of the form
:math:`\mathrm{tr}(C_i C_j C_k C_l)` as a sparse COO array for
:math:`i,j,k,l > 0` (i.e. excluding the identity).
"""
if self._four_element_traces is None:
# Most of the traces are zero, therefore store the result in a
# sparse array. For GGM bases, which are inherently sparse, it
# makes sense for any dimension to also calculate with sparse
# arrays. For Pauli bases, which are very dense, this is not so
# efficient but unavoidable for d > 12.
path = [(0, 1), (0, 1), (0, 1)]
if self.btype == 'Pauli' and self.d <= 12:
# For d == 12, the result is ~270 MB.
self._four_element_traces = COO.from_numpy(oe.contract('iab,jbc,kcd,lda->ijkl',
*(self,)*4, optimize=path))
else:
self._four_element_traces = oe.contract('iab,jbc,kcd,lda->ijkl', *(self.sparse,)*4,
backend='sparse', optimize=path)
return self._four_element_traces
@four_element_traces.setter
def four_element_traces(self, traces):
self._four_element_traces = traces
[docs] def normalize(self, copy: bool = False) -> Union[None, 'Basis']:
"""Normalize the basis."""
if copy:
return normalize(self)
self /= _norm(self)
[docs] def tidyup(self, eps_scale: Optional[float] = None) -> None:
"""Wraps util.remove_float_errors."""
if eps_scale is None:
atol = self._atol
else:
atol = self._eps*eps_scale
self.real[np.abs(self.real) <= atol] = 0
self.imag[np.abs(self.imag) <= atol] = 0
[docs] @classmethod
def pauli(cls, n: int) -> 'Basis':
r"""
Returns a Pauli basis for :math:`n` qubits, i.e. the basis spans
the space :math:`\mathbb{C}^{d\times d}` with :math:`d = 2^n`:
.. math::
\mathcal{P} = \{I, X, Y, Z\}^{\otimes n}.
The elements :math:`\sigma_i` are normalized with respect to the
Hilbert-Schmidt inner product,
.. math::
\langle\sigma_i,\sigma_j\rangle
&= \mathrm{Tr}\,\sigma_i^\dagger\sigma_j \\
&= \delta_{ij}.
Parameters
----------
n: int
The number of qubits.
Returns
-------
basis: Basis
The Basis object representing the Pauli basis.
"""
normalization = np.sqrt(2**n)
combinations = np.indices((4,)*n).reshape(n, 4**n)
sigma = util.tensor(*util.paulis[combinations], rank=2)
sigma /= normalization
return cls(sigma, btype='Pauli',
labels=[''.join(tup) for tup in product(['I', 'X', 'Y', 'Z'], repeat=n)])
[docs] @classmethod
def ggm(cls, d: int) -> 'Basis':
r"""
Returns a generalized Gell-Mann basis in :math:`d` dimensions
[Bert08]_ where the elements :math:`\Lambda_i` are normalized
with respect to the Hilbert-Schmidt inner product,
.. math::
\langle\Lambda_i,\Lambda_j\rangle
&= \mathrm{Tr}\,\Lambda_i^\dagger\Lambda_j \\
&= \delta_{ij}.
Parameters
----------
d: int
The dimensionality of the space spanned by the basis
Returns
-------
basis: Basis
The Basis object representing the GGM.
References
----------
.. [Bert08]
Bertlmann, R. A., & Krammer, P. (2008). Bloch vectors for
qudits. Journal of Physics A: Mathematical and Theoretical,
41(23). https://doi.org/10.1088/1751-8113/41/23/235303
"""
n_sym = int(d*(d - 1)/2)
sym_rng = np.arange(1, n_sym + 1)
diag_rng = np.arange(1, d)
# Indices for offdiagonal elements
j = np.repeat(np.arange(d - 1), np.arange(d - 1, 0, -1))
k = np.arange(1, n_sym+1) - (j*(2*d - j - 3)/2).astype(int)
j_offdiag = tuple(j)
k_offdiag = tuple(k)
# Indices for diagonal elements
j_diag = tuple(i for l in range(d) for i in range(l))
l_diag = tuple(i for i in range(1, d))
inv_sqrt2 = 1/np.sqrt(2)
Lambda = np.zeros((d**2, d, d), dtype=complex)
Lambda[0] = np.eye(d)/np.sqrt(d)
# First n matrices are symmetric
Lambda[sym_rng, j_offdiag, k_offdiag] = inv_sqrt2
Lambda[sym_rng, k_offdiag, j_offdiag] = inv_sqrt2
# Second n matrices are antisymmetric
Lambda[sym_rng+n_sym, j_offdiag, k_offdiag] = -1j*inv_sqrt2
Lambda[sym_rng+n_sym, k_offdiag, j_offdiag] = 1j*inv_sqrt2
# Remaining matrices have entries on the diagonal only
Lambda[np.repeat(diag_rng, diag_rng)+2*n_sym, j_diag, j_diag] = 1
Lambda[diag_rng + 2*n_sym, l_diag, l_diag] = -diag_rng
# Normalize
Lambda[2*n_sym + 1:, range(d), range(d)] /= np.tile(
np.sqrt(diag_rng*(diag_rng + 1))[:, None], (1, d)
)
return cls(Lambda, btype='GGM', labels=[rf'$\Lambda_{{{i}}}$' for i in range(len(Lambda))])
[docs] @classmethod
def from_partial(cls, partial_basis_array: Sequence,
traceless: Optional[bool] = None,
btype: Optional[str] = None,
labels: Optional[Sequence[str]] = None) -> 'Basis':
r"""Generate complete and orthonormal basis from a partial set.
The basis is completed using singular value decomposition to
determine the null space of the expansion coefficients of the
partial basis with respect to another complete basis.
Parameters
----------
partial_basis_array: array_like
A sequence of basis elements.
traceless: bool, optional
If a traceless basis should be generated (i.e. the first
element is the identity and all the others have trace zero).
btype: str, optional
A custom identifier.
labels: Sequence[str], optional
A list of custom labels for each element. If
`len(labels) == len(partial_basis_array)`, the newly created
elements get labels 'C_i'.
Returns
-------
basis: Basis, shape (d**2, d, d)
The orthonormal basis.
Raises
------
ValueError
If the given elements are not orthonormal.
ValueError
If the given elements are not traceless but
`traceless==True`.
ValueError
If not len(partial_basis_array) or d**2 labels were given.
"""
if btype is None:
btype = 'From partial'
if (labels is None
and hasattr(partial_basis_array, 'labels')
and len(partial_basis_array.labels) == len(partial_basis_array)):
# Need to check if labels and array are same length as indexing
# is unaware of our custom attributes
labels = partial_basis_array.labels
basis, labels = _full_from_partial(partial_basis_array, traceless, labels)
return cls(basis, btype=btype, labels=labels)
def _full_from_partial(elems: Sequence, traceless: bool, labels: Sequence[str]) -> Basis:
"""
Internal function to parse the basis elements *elems*. By default,
checks are performed for orthogonality and linear independence. If
either fails an exception is raised. Returns a full hermitian and
orthonormal basis.
"""
# Convert elems to basis to have access to its handy attributes
elems = Basis(elems)
elems.normalize(copy=False)
if not elems.isherm:
warn("(Some) elems not hermitian! The resulting basis also won't be.")
if not elems.isorthonorm:
raise ValueError("The basis elements are not orthonormal!")
if traceless is None:
traceless = elems.istraceless
else:
if traceless and not elems.istraceless:
raise ValueError("The basis elements are not traceless (up to an identity element) "
+ "but a traceless basis was requested!")
if labels is not None and len(labels) not in (len(elems), elems.d**2):
raise ValueError(f'Got {len(labels)} labels but expected {len(elems)} or {elems.d**2}')
# Get a Generalized Gell-Mann basis to expand in (fulfills the desired
# properties hermiticity and orthonormality, and therefore also linear
# combinations, ie basis expansions, of it will). Split off the identity so
# that for traceless bases we can put it in the front.
if traceless:
Id, ggm = np.split(Basis.ggm(elems.d), [1])
else:
ggm = Basis.ggm(elems.d)
coeffs = expand(elems, ggm, hermitian=elems.isherm, tidyup=True)
# Throw out coefficient vectors that are all zero (should only happen for
# the identity)
coeffs = coeffs[(coeffs != 0).any(axis=-1)]
if coeffs.size != 0:
# Get d**2 - len(coeffs) vectors spanning the nullspace of coeffs.
# Those together with coeffs span the whole space, and therefore also
# the linear combinations of GGMs weighted with the coefficients will
# span the whole matrix space
coeffs = np.concatenate((coeffs, sla.null_space(coeffs).T))
# Our new basis is given by linear combinations of GGMs with coeffs
basis = np.einsum('ij,jkl', coeffs, ggm)
else:
# Resulting array is of size zero, i.e. we can just return the GGMs
basis = ggm
# Add the identity again and normalize the new basis
if traceless:
basis = np.concatenate((Id, basis)).view(Basis)
else:
basis = basis.view(Basis)
# Clean up
basis.tidyup()
if labels is not None and len(labels) == len(elems):
# Fill up labels for newly generated elements
labels = list(labels)
if traceless:
# sort Identity label to the front, default to first if not found
# (should not happen since traceless checks that it is present)
id_idx = next((i for i, elem in enumerate(elems)
if np.allclose(Id.view(ndarray), elem.view(ndarray),
rtol=elems._rtol, atol=elems._atol)), 0)
labels.insert(0, labels.pop(id_idx))
labels.extend('$C_{{{}}}$'.format(i) for i in range(len(labels), len(basis)))
return basis, labels
def _norm(b: Sequence) -> ndarray:
"""Frobenius norm with two singleton dimensions inserted at the end."""
b = np.asanyarray(b)
norm = nla.norm(b, axis=(-1, -2))
return norm[..., None, None]
[docs]def normalize(b: Basis) -> Basis:
r"""
Return a copy of the basis *b* normalized with respect to the
Frobenius norm [Gol85]_:
:math:`||A||_F = \left[\sum_{i,j} |a_{i,j}|^2\right]^{1/2}`
or equivalently, with respect to the Hilbert-Schmidt inner product
as implemented by :func:`~filter_functions.util.dot_HS`.
References
----------
.. [Gol85]
G. H. Golub and C. F. Van Loan, *Matrix Computations*,
Baltimore, MD, Johns Hopkins University Press, 1985, pg. 15
"""
return (b/_norm(b)).squeeze().view(Basis)
[docs]def expand(M: Union[ndarray, Basis], basis: Union[ndarray, Basis],
normalized: bool = True, hermitian: bool = False, tidyup: bool = False) -> ndarray:
r"""
Expand the array *M* in the basis given by *basis*.
Parameters
----------
M: array_like
The square matrix (d, d) or array of square matrices (..., d, d)
to be expanded in *basis*
basis: array_like
The basis of shape (m, d, d) in which to expand.
normalized: bool {True}
Wether the basis is normalized.
hermitian: bool (default: False)
If M is hermitian along its last two axes, the result will be
real.
tidyup: bool {False}
Whether to set values below the floating point eps to zero.
Returns
-------
coefficients: ndarray
The coefficient array with shape (..., m) or (m,) if *M* was 2-d
Notes
-----
For an orthogonal matrix basis :math:`\mathcal{C} = \big\{C_k\in
\mathbb{C}^{d\times d}: \langle C_k,C_l\rangle_\mathrm{HS} \propto
\delta_{kl}\big\}_{k=0}^{d^2-1}` with the Hilbert-Schmidt inner
product as implemented by :func:`~filter_functions.util.dot_HS` and
:math:`M\in\mathbb{C}^{d\times d}`, the expansion of
:math:`M` in terms of :math:`\mathcal{C}` is given by
.. math::
M &= \sum_j c_j C_j, \\
c_j &= \frac{\mathrm{tr}\big(M C_j\big)}
{\mathrm{tr}\big(C_j^\dagger C_j\big)}.
"""
def cast(arr):
return arr.real if hermitian and basis.isherm else arr
coefficients = cast(np.tensordot(M, basis, axes=[(-2, -1), (-1, -2)]))
if not normalized:
coefficients /= cast(np.einsum('bij,bji->b', basis, basis))
return util.remove_float_errors(coefficients) if tidyup else coefficients
[docs]def ggm_expand(M: Union[ndarray, Basis], traceless: bool = False,
hermitian: bool = False) -> ndarray:
r"""
Expand the matrix *M* in a Generalized Gell-Mann basis [Bert08]_.
This function makes use of the explicit construction prescription of
the basis and thus makes do without computing the expansion
coefficients as the overlap between the matrix and each basis
element.
Parameters
----------
M: array_like
The square matrix (d, d) or array of square matrices (..., d, d)
to be expanded in a GGM basis.
traceless: bool (default: False)
Include the basis element proportional to the identity in the
expansion. If it is known beforehand that M is traceless, the
corresponding coefficient is zero and thus doesn't need to be
computed.
hermitian: bool (default: False)
If M is hermitian along its last two axes, the result will be
real.
Returns
-------
coefficients: ndarray
The real coefficient array with shape (d**2,) or (..., d**2)
References
----------
.. [Bert08]
Bertlmann, R. A., & Krammer, P. (2008). Bloch vectors for
qudits. Journal of Physics A: Mathematical and Theoretical,
41(23). https://doi.org/10.1088/1751-8113/41/23/235303
"""
if M.shape[-1] != M.shape[-2]:
raise ValueError('M should be square in its last two axes')
def cast(arr):
return arr.real if hermitian else arr
# Squeeze out an extra dimension to be shape agnostic
square = M.ndim < 3
if square:
M = M[None, ...]
d = M.shape[-1]
n_sym = int(d*(d - 1)/2)
sym_rng = np.arange(1, n_sym + 1)
diag_rng = np.arange(1, d)
# Map linear index to index tuple of upper triangle
j = np.repeat(np.arange(d - 1), np.arange(d - 1, 0, -1))
k = np.arange(1, n_sym+1) - (j*(2*d - j - 3)/2).astype(int)
offdiag_idx = [tuple(j), tuple(k)]
# Indices for upper triangular part of M
triu_idx = tuple([...] + offdiag_idx)
# Indices for lower triangular part of M
tril_idx = tuple([...] + offdiag_idx[::-1])
# Indices for diagonal elements of M up to first to last
diag_idx = tuple([...] + [tuple(i for i in range(d - 1))]*2)
# Indices of diagonal elements of M starting from first
diag_idx_shifted = tuple([...] + [tuple(i for i in range(1, d))]*2)
# Compute the coefficients
coeffs = np.zeros((*M.shape[:-2], d**2), dtype=float if hermitian else complex)
if not traceless:
# First element is proportional to the trace of M
coeffs[..., 0] = cast(np.einsum('...jj', M))/np.sqrt(d)
# Elements proportional to the symmetric GGMs
coeffs[..., sym_rng] = cast(M[triu_idx] + M[tril_idx])/np.sqrt(2)
# Elements proportional to the antisymmetric GGMs
coeffs[..., sym_rng + n_sym] = cast(1j*(M[triu_idx] - M[tril_idx]))/np.sqrt(2)
# Elements proportional to the diagonal GGMs
coeffs[..., diag_rng + 2*n_sym] = cast(M[diag_idx].cumsum(axis=-1)
- diag_rng*M[diag_idx_shifted])
coeffs[..., diag_rng + 2*n_sym] /= cast(np.sqrt(diag_rng*(diag_rng + 1)))
return coeffs.squeeze() if square else coeffs
def equivalent_pauli_basis_elements(idx: Union[Sequence[int], int], N: int) -> ndarray:
"""
Get the indices of the equivalent (up to identities tensored to it)
basis elements of Pauli bases of qubits at position idx in the total
Pauli basis for N qubits.
"""
idx = [idx] if isinstance(idx, int) else idx
multi_index = np.ix_(*[range(4) if i in idx else [0]
for i in range(N)])
elem_idx = np.ravel_multi_index(multi_index, [4]*N).ravel()
return elem_idx
def remap_pauli_basis_elements(order: Sequence[int], N: int) -> ndarray:
"""
For a N-qubit Pauli basis, transpose the order of the subsystems and
return the indices that permute the old basis to the new.
"""
# Index tuples for single qubit paulis that give the n-qubit paulis when
# tensored together
pauli_idx = np.indices((4,)*N).reshape(N, 4**N).T
# Indices of the N-qubit basis reordered according to order
linear_idx = [np.ravel_multi_index([idx_tup[i] for i in order], (4,)*N)
for idx_tup in pauli_idx]
return np.array(linear_idx)