Commit 2861162f authored by Tobias Hangleiter's avatar Tobias Hangleiter
Browse files

Add linalg and qi modules

parent 985d587e
This module contains utility functions for linear algebra manipulations
:meth:`abs2` :
Fast absolute value squared
Fast tensor product of an arbitrary number of inputs using
Multiple matrix product along an axis
Partial trace on a bipartite system
Maximum absolute difference of two arrays
Maximum relative difference of two arrays
Fast evaluation of a Pauli exponential
Set entries whose absolute value is below a certain threshold to zero
Determine if two vectors or operators are equal up to a global phase
Hilbert-Schmidt inner product
from functools import reduce
from math import copysign
from typing import List, Sequence, Tuple, Union
import numpy as np
import qutip as qt
from numpy import ndarray
from numpy.linalg import norm
from .qi import P_np
import numba as nb
numba = True
except ImportError:
numba = False
def abs2(arr):
Fast function to calculate the absolute value squared,
.. math::
|\cdot|^2 = \Re(\cdot)^2 + \Im(\cdot)^2.
Equivalent to::
arr : array_like
Array-like object to take the absolute value square of. Needs to
provide a ``real`` and ``imag`` attribute.
return arr.real**2 + arr.imag**2
# Vectorize if numba is available
if numba:
abs2 = nb.vectorize(
[nb.float64(nb.complex128), nb.float32(nb.complex64)]
def tensor(*args):
Fast tensor product using einsum. The arguments may be of arbitrary
dimension but are required to be square on its last two axes and the shapes
of the remaining axes must be the same. For example, the following shapes
are compatible::
(a, b, c, d, d), (a, b, c, e, e) -> (a, b, c, d*e, d*e)
if len(set(arg.shape[:-2] for arg in args)) != 1:
raise ValueError('Require all args to have the same shape except ' +
'on the last two axes.')
if not all(arg.shape[-2] == arg.shape[-1] for arg in args):
raise ValueError('Require all args to be square in its last two axes.')
def binary_tensor(A, B):
d = A.shape[-1]*B.shape[-1]
return np.einsum('...ij,...kl->...ikjl', A, B).reshape(*A.shape[:-2],
d, d)
return reduce(binary_tensor, args)
def mdot(arr: Sequence, axis: int = 0) -> ndarray:
"""Multiple matrix products along axis"""
return reduce(np.matmul, arr.swapaxes(0, axis))
def ptrace(arr: ndarray, dims: Sequence[int], system: int = 1) -> ndarray:
Partial trace of bipartite matrix.
arr : ndarray
Matrix of bipartite system on the last two axes. There are no
constraints on the shape of the other axes.
.. math::
\verb|arr|\equiv A\otimes B
dims : Sequence of ints
The dimensions of system A and B
system : int
The system to trace out, either 0 (A) or 1 (B)
>>> d_A, d_B = 3, 2
>>> A = np.arange(d_A**2).reshape(1, d_A, d_A)
>>> B = np.arange(d_B**2).reshape(1, d_B, d_B)
>>> C = tensor(A, B)
>>> ptrace(C, (d_A, d_B), 0)
array([[[ 0, 12],
[24, 36]]])
>>> ptrace(C, (d_A, d_B), 1)
array([[[ 0, 3, 6],
[ 9, 12, 15],
[18, 21, 24]]])
if system == 0:
einsum_str = '...ijik'
elif system == 1:
einsum_str = '...ikjk'
return np.einsum(einsum_str, arr.reshape(*arr.shape[:-2], *dims, *dims))
def pauli_expm(a: Sequence, nonzero_ind: List[bool]) -> ndarray:
Faster method of calculating a Pauli exponential
.. math::
\exp(-i(\vec{a}\cdot\vec{\sigma})) = I\cos|\vec{a}|
- i(\vec{a}\cdot\vec{\sigma})\sin|\vec{a}|
than :meth:`qutip.Qobj.expm()` or :meth:`scipy.linalg.expm`.
a : array_like, shape (n, ...)
Array with 0 < n <= 3 cartesian components on the first axis.
nonzero_ind : list of three bools
List of three Booleans indicating those cartesian components that are
nonzero, i.e., which of the three (x, y, z) are the first axis of *a*.
Accordingly, ``sum(nonzero_ind) == a.shape[0]``.
A Hadamard gate, i.e. a rotation by pi about X + Z (dividing by sqrt(2)
normalizes the rotation axis vector):
>>> a = np.full((2,), np.pi/2/np.sqrt(2))
>>> pauli_expm(a, [True, False, True]) # a[0]=a[2]=pi/2/sqrt(2), a[1]=0
if not sum(nonzero_ind) == a.shape[0]:
raise ValueError('nonzero_ind should have the same number of True '
f'entries as a.shape[0] = {a.shape[0]}!')
a = np.asarray(a)
P = np.asarray(P_np)
idx = np.asarray(nonzero_ind)
# Twice the rotation angle
theta = norm(a, axis=0)
# Rotation axis
n = a/theta
# Replace nans originating from divide by zero
n[np.isnan(n)] = 0
real = np.einsum('ij,...->ij...', P[0], np.cos(theta))
imag = np.einsum('i...,ijk,...->jk...', n, P[1:][idx], np.sin(theta))
result = real - 1j*imag
return result.transpose([2, 0, 1]) if result.ndim == 3 else result
def remove_float_errors(arr: ndarray, eps_scale: float = None):
Clean up arr by removing floating point numbers smaller than the dtype's
precision multiplied by eps_scale. Treats real and imaginary parts
if eps_scale is None:
atol = np.finfo(arr.dtype).eps*arr.shape[-1]
atol = np.finfo(arr.dtype).eps*eps_scale
# Hack around arr.imag sometimes not being writable
if arr.dtype == complex:
arr = arr.real + 1j*arr.imag
arr.real[np.abs(arr.real) <= atol] = 0
arr.imag[np.abs(arr.imag) <= atol] = 0
arr = arr.real
arr.real[np.abs(arr.real) <= atol] = 0
return arr
def check_phase_eq(psi: Union[qt.Qobj, ndarray],
phi: Union[qt.Qobj, ndarray],
eps: float = None,
normalized: bool = False) -> Tuple[bool, float]:
Checks whether psi and phi are equal up to a global phase, i.e.
.. math::
|\psi\rangle = e^{i\chi}|\phi\rangle \Leftrightarrow
\langle \phi|\psi\rangle = e^{i\chi},
and returns the phase. If the first return value is false, the second is
meaningless in this context. psi and phi can also be operators.
psi, phi : Qobj or ndarray
Vectors or operators to be compared
eps : float
The tolerance below which the two objects are treated as equal, i.e.,
the function returns ``True`` if ``abs(1 - modulus) <= eps``.
normalized : bool
Flag indicating if *psi* and *phi* are normalized with respect to the
Hilbert-Schmidt inner product :meth:`dot_HS`.
>>> psi = qt.sigmax()
>>> phi = qt.sigmax()*np.exp(1j*1.2345)
>>> check_phase_eq(psi, phi)
(True, 1.2345)
d = max(psi.shape)
psi, phi = [obj.full() if isinstance(obj, qt.Qobj) else obj
for obj in (psi, phi)]
if eps is None:
# Tolerance the floating point eps times the dimension squared
eps = max(np.finfo(psi.dtype).eps, np.finfo(phi.dtype).eps)*d**2
if psi.ndim - psi.shape.count(1) == 1:
# Vector
if normalized:
inner_product = (psi.T.conj() @ phi).squeeze()
inner_product = ((psi.T.conj() @ phi) /
elif psi.ndim == 2:
# Matrix
if normalized:
inner_product = dot_HS(psi, phi, eps)
inner_product = dot_HS(psi, phi, eps) /\
np.sqrt(dot_HS(psi, psi, eps)*dot_HS(phi, phi, eps))
raise ValueError('Invalid dimension')
phase = np.angle(inner_product)
modulus = abs(inner_product)
return abs(1 - modulus) <= eps, phase
def dot_HS(U: Union[ndarray, qt.Qobj],
V: Union[ndarray, qt.Qobj],
eps: float = None) -> float:
r"""Return the Hilbert-Schmidt inner product of U and V,
.. math::
\langle U, V\rangle_\mathrm{HS} := \mathrm{tr}(U^\dagger V).
U, V : Qobj or ndarray
Objects to compute the inner product of.
>>> U, V = qt.sigmax(), qt.sigmay()
>>> dot_HS(U, V)
>>> dot_HS(U, U)
if isinstance(U, qt.Qobj):
U = U.full()
if isinstance(V, qt.Qobj):
V = V.full()
res = np.einsum('ij,ij', U.conj(), V)
if eps is None:
# Tolerance is the dtype precision times the number of flops for the
# matrix multiplication
eps = np.finfo(U.dtype).eps**V.shape[-1]
except ValueError:
# dtype is int and therefore exact
eps = 0
# Deal with real and imaginary part separately
if abs(res.real) <= eps:
real = 0.0
elif abs(1 - abs(res.real)) <= eps:
real = copysign(1, res.real)
real = res.real
if abs(res.imag) <= eps:
imag = 0
elif abs(1 - abs(res.imag)) <= eps:
imag = copysign(1, res.imag)
imag = res.imag
return real + 1j*imag if imag != 0 else real
def sparsity(arr: ndarray, eps: float = None) -> float:
Return the sparsity of the array *arr*.
arr: array_like
The array
eps: float (default: dtype eps)
Entries smaller than this value will be treated as zero
sparsity: float
The sparsity of the array
eps = np.finfo(arr.dtype).eps if eps is None else eps
return (np.abs(arr) <= eps).sum()/arr.size
def density(arr: ndarray, eps: float = None) -> float:
Return the density of the array *arr*.
arr: array_like
The array
eps: float (default: dtype eps)
Entries smaller than this value will be treated as zero
density: float
The density of the array
eps = np.finfo(arr.dtype).eps if eps is None else eps
return (np.abs(arr) > eps).sum()/arr.size
def max_abs_diff(a, b):
"""Maximum absolute difference"""
return np.max(np.abs(a - b))
def max_rel_diff(a, b):
"""Maximum relative difference"""
return np.nanmax(np.abs((a - b)/np.abs(b)))
"""This module contains utilities pertaining to quantum information"""
import numpy as np
import qutip as qt
# Pauli matrices ...
# ... as QuTiP Qobj's
P_qt = [qt.qeye(2),
# ... as NumPy arrays
P_np = [P.full() for P in P_qt]
# ... as sparse matrices
P_sp = [ for P in P_qt]
# Clifford group generators
clifford_generators = [qt.snot(), qt.phasegate(np.pi/2).tidyup(), qt.cnot()]
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment