"""Different math functions needed at some point in the library."""
# Copyright 2018-2019 TeNPy Developers, GNU GPLv3
import numpy as np
import warnings
from . import misc
import scipy.linalg
import scipy.sparse.linalg
int_I3 = np.eye(3, dtype=int)
LeviCivita3 = np.array([[np.cross(b, a) for a in int_I3] for b in int_I3])
__all__ = [
'matvec_to_array', 'entropy', 'gcd', 'gcd_array', 'lcm', 'speigs', 'speigsh', 'perm_sign',
'qr_li', 'rq_li'
]
[docs]def matvec_to_array(H):
"""transform an linear operator with a `matvec` method into a dense numpy array.
Parameters
----------
H : linear operator
should have `shape`, `dtype` attributes and a `matvec` method.
Returns
-------
H_dense : ndarray, shape ``(H.dim, H.dim)``
a dense array version of `H`.
"""
dim, dim2 = H.shape
assert (dim == dim2)
X = np.zeros((dim, dim), H.dtype)
v = np.zeros((dim), H.dtype)
for i in range(dim):
v[i] = 1
X[:, i] = H.matvec(v)
v[i] = 0
return X
##########################################################################
##########################################################################
# Actual Math functions
[docs]def entropy(p, n=1):
r"""Calculate the entropy of a distribution.
Assumes that p is a normalized distribution (``np.sum(p)==1.``).
Parameters
----------
p : 1D array
A normalized distribution.
n : 1 | float | np.inf
Selects the entropy, see below.
Returns
-------
entropy : float
Shannon-entropy :math:`-\sum_i p_i \log(p_i)` (n=1) or
Renyi-entropy :math:`\frac{1}{1-n} \log(\sum_i p_i^n)` (n != 1)
of the distribution `p`.
"""
p = p[p > 1.e-30] # just for stability reasons / to avoid NaN in log
if n == 1:
return -np.inner(np.log(p), p)
elif n == np.inf:
return -np.log(np.max(p))
else: # general n != 1, inf
return np.log(np.sum(p**n)) / (1. - n)
[docs]def gcd(a, b):
"""Computes the greatest common divisor (GCD) of two numbers.
Return 0 if both a, b are zero, otherwise always return a non-negative number.
"""
a, b = abs(a), abs(b)
while b:
a, b = b, a % b
return a
[docs]def gcd_array(a):
"""Return the greatest common divisor of all of entries in `a`"""
a = np.array(a).reshape(-1)
if len(a) <= 0:
raise ValueError
t = a[0]
for x in a[1:]:
if t == 1:
break
t = gcd(t, x)
return t
[docs]def lcm(a, b):
"""Returns the least common multiple (LCM) of two positive numbers."""
a0, b0 = a, b
while b:
a, b = b, a % b
return a0 * (b0 // a)
[docs]def speigs(A, k, *args, **kwargs):
"""Wrapper around :func:`scipy.sparse.linalg.eigs`, lifting the restriction ``k < rank(A)-1``.
Parameters
----------
A : MxM ndarray or like :class:`scipy.sparse.linalg.LinearOperator`
the (square) linear operator for which the eigenvalues should be computed.
k : int
the number of eigenvalues to be computed.
*args, **kwargs :
further arguments are directly given to :func:`scipy.sparse.linalg.eigs`
Returns
-------
w : ndarray
array of min(`k`, A.shape[0]) eigenvalues
v : ndarray
array of min(`k`, A.shape[0]) eigenvectors, ``v[:, i]`` is the `i`-th eigenvector.
Only returned if ``kwargs['return_eigenvectors'] == True``.
"""
d = A.shape[0]
if A.shape != (d, d):
raise ValueError("A.shape not a square matrix: " + str(A.shape))
if k < d - 1:
return scipy.sparse.linalg.eigs(A, k, *args, **kwargs)
else:
if k > d:
warnings.warn("trimming speigs k to smaller matrix dimension d", stacklevel=2)
k = d
if isinstance(A, np.ndarray):
Amat = A
else:
Amat = matvec_to_array(A) # Constructs the matrix
ret_eigv = kwargs.get('return_eigenvectors', args[7] if len(args) > 7 else True)
which = kwargs.get('which', args[2] if len(args) > 2 else 'LM')
if ret_eigv:
W, V = np.linalg.eig(Amat)
keep = misc.argsort(W, which)[:k]
return W[keep], V[:, keep]
else:
W = np.linalg.eigvals(Amat)
keep = misc.argsort(W, which)[:k]
return W
[docs]def speigsh(A, k, *args, **kwargs):
"""Wrapper around :func:`scipy.sparse.linalg.eigsh`, lifting the restriction ``k < rank(A)-1``.
Parameters
----------
A : MxM ndarray or like :class:`scipy.sparse.linalg.LinearOperator`
The (square) hermitian linear operator for which the eigenvalues should be computed.
k : int
The number of eigenvalues to be computed.
*args, **kwargs :
Further arguments are directly given to :func:`scipy.sparse.linalg.eigs`.
Returns
-------
w : ndarray
Array of min(`k`, A.shape[0]) eigenvalues.
v : ndarray
Array of min(`k`, A.shape[0]) eigenvectors, ``v[:, i]`` is the `i`-th eigenvector.
Only returned if ``kwargs['return_eigenvectors'] == True``.
"""
d = A.shape[0]
if A.shape != (d, d):
raise ValueError("A.shape not a square matrix: " + str(A.shape))
if k < d - 1:
return scipy.sparse.linalg.eigsh(A, k, *args, **kwargs)
else:
if k > d:
warnings.warn("trimming speigsh k to smaller matrix dimension d", stacklevel=2)
k = d
if isinstance(A, np.ndarray):
Amat = A
else:
Amat = matvec_to_array(A) # Constructs the matrix
ret_eigv = kwargs.get('return_eigenvectors', args[7] if len(args) > 7 else True)
which = kwargs.get('which', args[2] if len(args) > 2 else 'LM')
if ret_eigv:
W, V = np.linalg.eigh(Amat)
keep = misc.argsort(W, which)[:k]
return W[keep], V[:, keep]
else:
W = np.linalg.eigvalsh(Amat)
keep = misc.argsort(W, which)[:k]
return W
[docs]def perm_sign(p):
"""Given a permutation `p` of numbers, returns its sign. (+1 or -1)
Assumes that all the elements are distinct, if not, you get crap.
Examples
--------
>>> for p in itertools.permutations(range(3))]):
... print('{p!s}: {sign!s}'.format(p=p, sign=perm_sign(p)))
(0, 1, 2): 1
(0, 2, 1): -1
(1, 0, 2): -1
(1, 2, 0): 1
(2, 0, 1): 1
(2, 1, 0): -1
"""
rp = np.argsort(p)
p = np.argsort(rp)
s = 1
for i, v in enumerate(p):
if i == v:
continue
# by the way we loop, i < v, so we find where i is.
# p[i] = p[rp[i]] # we don't have to do that becasue we never need p[i] again
p[rp[i]] = v
rp[v] = rp[i]
s = -s
return s
[docs]def qr_li(A, cutoff=1.e-15):
"""QR decomposition with cutoff to discard nearly linear dependent columns in `Q`.
Perform a QR decomposition with pivoting, discard columns where ``R[i,i] < cuttoff``,
reverse the permututation from pivoting and perform another QR decomposition to ensure that
`R` is upper right.
Parameters
----------
A : :class:`numpy.ndarray`
Matrix to be decomposed as ``A = Q.R``
Returns
-------
Q, R : :class:`numpy.ndarray`
Decomposition of `A` into isometry `Q^d Q = 1` and upper right `R` with diagonal entries
larger than `cutoff`.
"""
Q, R, P = scipy.linalg.qr(A, mode='economic', pivoting=True)
keep = np.abs(np.diag(R)) > cutoff
assert len(keep) == R.shape[0]
Q = Q[:, keep]
R = R[keep, :]
# here, A P = Q R, thus A = Q R inv(P)
R = R[:, misc.inverse_permutation(P)]
q, R = scipy.linalg.qr(R, mode='economic', pivoting=False)
return np.dot(Q, q), R
[docs]def rq_li(A, cutoff=1.e-15):
"""RQ decomposition with cutoff to discard nearly linear dependent columns in `Q`.
Uses :func:`qr_li` on tranpose of `A`.
Note that `R` is nonzero in the lowest left corner; `R` has entries below the diagonal
for non-square `R`.
Parameters
----------
A : :class:`numpy.ndarray`
Matrix to be decomposed as ``A = Q.R``
Returns
-------
R, Q : :class:`numpy.ndarray`
Decomposition of `A` into isometry `Q Q^d = 1` and upper right `R` with diagonal entries
larger than `cutoff`. If ``M, N = A.shape``, then ``R.shape = M, K`` and ``Q.shape = K, N``
with ``K <= min(M, N)``.
"""
q, r = qr_li(A.T[:, ::-1], cutoff)
return r.T[::-1, ::-1], q.T[::-1, :]