Source code for Inelastica.math.spectral

from __future__ import print_function

import numpy as N
import numpy.linalg as LA


[docs] def mm(* args): """ Matrix multiplication with arbitrary number of arguments and the SpectralMatrix type. """ args = list(args) # Look for SpectralMatrices where = N.where(N.array([isinstance(ii, SpectralMatrix) for ii in args]))[0] if len(where) > 0: res = __mmSpectralMatrix(args, where) else: res = __mm(args) return res
def __mmSpectralMatrix(args, where): # find smallest size = [args[where[ii]].L.shape[1] for ii in range(len(where))] smallest = N.where(N.array(size) == N.min(size))[0][0] # if more than one SpectralMatrix, expand others into normal matrices args = [[ii] for ii in args] # To get flatten to work for ii in where: if ii != where[smallest]: args[ii] = [args[ii][0].L, args[ii][0].R] args = sum(args, []) # Flatten list # Split at Spectral matrix where = N.where(N.array([isinstance(ii, SpectralMatrix) for ii in args]))[0] res = SpectralMatrix() res.L = __mm(args[:where[0]]+[args[where[0]].L]) res.R = __mm([args[where[0]].R]+args[where[0]+1:]) return res def __mm(args): # Normal matrix mult, order important for speed if matrices has different sizes if len(args) == 1: return args[0] # Find smallest matrix Lsize = N.array([ii.shape[0] for ii in args[:-1]]) if len(args[-1].shape) == 1: # Most rightmost is vector Rsize = N.array([ii.shape[1] for ii in args[1:-1]]+[1]) else: Rsize = N.array([ii.shape[1] for ii in args[1:]]) Lwhere = N.where(Lsize == N.min(Lsize))[0] Rwhere = N.where(Rsize == N.min(Rsize))[0] if N.min(Lsize) > N.min(Rsize): where = [ii+1 for ii in Rwhere] left = False else: where = Lwhere left = True if len(where) == len(args)-1 or len(args) == 2: # Order does not matter res = N.dot(args[0], args[1]) for ii in range(len(args)-2): res = N.dot(res, args[ii+2]) else: # Make sure the matrix mult is done in order of smaller # matrices first ii = where[0] if ii == 0: res = __mm([__mm([args[0], args[1]])]+args[2:]) elif ii == len(args)-1: res = __mm(args[:ii-1]+[__mm([args[ii-1], args[ii]])]) else: if left: res = __mm(args[:ii]+[__mm([args[ii]]+args[ii+1:])]) else: res = __mm([__mm(args[:ii]+[args[ii]])]+args[ii+1:]) return res
[docs] class SpectralMatrix(object): r""" Matrix class for spectral matrices. Idea: Split a spectral matrix :math:`\mathbf{A}^{(mxm)}` into two smaller matrices :math:`\mathbf{A} = \mathbf{L}.\mathbf{R}` where :math:`\mathbf{L}^{(mxn)} = [\lambda_1 \mathbf{v}_1, ..., \lambda_n \mathbf{v}_n]` and :math:`\mathbf{R}^{(nxm)} = [\mathbf{v}_1, ..., \mathbf{v}_n]^\dagger` are constructed from a call to ``LA.eigh(A)`` keeping only :math:`n` eigensolutions with eigenvalues above the specified cutoff. The magic takes place in Matrix multiply function. Todo: Adding two gives full matrix! Should be easy to fix. Attributes ---------- L : ndarray R : ndarray Parameters ---------- A : ndarray """ # self.L/R : Left / right matrices
[docs] def __init__(self, A=None, cutoff=1e-8): if isinstance(A, N.ndarray): # Initialize ... only Hermitian matrices ev, evec = LA.eigh(A) # Drop eigenvalues indx = N.where(N.abs(ev) > cutoff)[0] ev, evec = ev[indx], evec[:, indx] print("SpectralMatrix: Fraction of eigenvalues above cutoff (%.1e) is %i/%i"%(cutoff, len(indx), len(A))) self.L = N.dot(evec, N.diag(ev)) self.R = dagger(evec)
# print N.allclose(A, N.dot(self.L, self.R))
[docs] def full(self): r""" Returns the dense ndarray via matrix multiplication :math:`\mathbf{A}^{(mxm)} = \mathbf{L}.\mathbf{R}`. """ return mm(self.L, self.R)
def __add__(self, b, subtract=False): # Could be improved for addition of two spectral matrices if not isinstance(b, SpectralMatrix): return self.full()+b else: NN = self.L.shape[0] Na, Nb = self.L.shape[1], b.L.shape[1] if Na+Nb > 0.3*NN: # Too many eigenvalues print("Too many eigenvalues") if not subtract: return self.full()+b.full() else: return self.full()-b.full() else: res = SpectralMatrix() res.L = N.zeros((NN, Na+Nb), N.complex) # Always assume complex !? res.R = N.zeros((Na+Nb, NN), N.complex) # Always assume complex !? res.L[:, 0:Na], res.R[0:Na, :] = self.L, self.R if not subtract: res.L[:, Na:Na+Nb], res.R[Na:Na+Nb, :] = b.L, b.R else: res.L[:, Na:Na+Nb], res.R[Na:Na+Nb, :] = -b.L, b.R return res def __radd__(self, b): return self+b def __sub__(self, b): if not isinstance(b, SpectralMatrix): return self.full()-b else: return self.__add__(b, subtract=True) def __rsub__(self, b): if not isinstance(b, SpectralMatrix): return b-self.full() else: return b.__add__(self, subtract=True) def __mul__(self, b): tmp = SpectralMatrix() tmp.L, tmp.R = self.L*b, self.R return tmp def __rmul__(self, b): return self*b def __dagger__(self): tmp = SpectralMatrix() tmp.L = dagger(self.R) tmp.R = dagger(self.L) return tmp
[docs] def trace(a): """ Returns the trace of a normal or spectral matrix. """ if isinstance(a, SpectralMatrix): return N.trace(mm(a.R, a.L)) # Switch sum around to make faster! else: return N.trace(a)
[docs] def dagger(x): """ Returns the hermitian conjugation of a normal or spectral matrix. """ if isinstance(x, SpectralMatrix): return x.__dagger__() else: return N.transpose(N.conjugate(x))