Hi, I am trying to find the eigenvalues and eigenvectors as well as
the inverse for a large number of small matrices. The matrix size
(MxM) will typically range from 2x2 to 8x8 at most. The number of
matrices (N) can be from 100 up to a million or more. My current
solution is to define "eig" and "inv" to be,
def inv(A):
"""
Inverts N MxM matrices, A.shape = (M, M, N), inv(A).shape = (M, M, N).
"""
return np.array(map(np.linalg.inv, A.transpose(2, 0, 1))).transpose(1, 2, 0)
def eig(A):
"""
Calculate the eigenvalues and eigenvectors of N MxM matrices,
A.shape = (M, M, N), eig(A)[0].shape = (M, N), eig(A)[1].shape = (M,
M, N)
"""
tmp = zip(*map(np.linalg.eig, A.transpose(2, 0, 1)))
return (np.array(tmp[0]).swapaxes(0,1), np.array(tmp[1]).transpose(1,2,0))
The above uses "map" to fake a vector solution, but this is heinously
slow. Are there any better ways to do this without resorting to cython
or weave (would it even be faster (or possible) to use "np.linalg.eig"
and "np.linalg.inv" within cython)? I could write specialized versions
when M=2 or M=3, which could be fully vectorized, but I'd rather have
a general solution. Are there better algorithms than the ones used in
"np.linalg.inv" and "np.linalg.eig" when M < 9, say, that I could hand
code using numpy in a fully vectorized way?
My end goal is to implement a Riemann flux calculation as part of a
finite volume solver. This requires calculating "Abar = R E inv(R)"
where R are the right eigenvectors of a matrix A and E is the matrix
with the absolute value of the eigenvalues of A along the diagonal
(both R and E are sorted in ascending order of their corresponding
eigenvalues). As well as using "eig" and "inv" defined above, matrix
multiplication and sorted eigenvalues and eigenvectors are also
required. Fortunately, these can be vectorized as follows,
def sum(a, axis=0):
"""
Faster than using np.sum.
"""
return np.tensordot(np.ones(a.shape[axis], 'l'), a, (0, axis))
def mul(A, B):
"""
Matrix multiply N MxM matrices, A.shape = B.shape = (M, M, N),
mul(A, B).shape = (M, M, N)
"""
return sum(A.swapaxes(0,1)[:, :, np.newaxis] * B[:, np.newaxis], 0)
def sortedeig(A):
"""
Calculates the sorted eigenvalues and eigenvectors of N MxM
matrices, A.shape = (M, M, N), sortedeig(A)[0].shape = (M, N),
sortedeig(A)[1].shape = (M, M, N).
"""
N = A.shape[-1]
eigenvalues, R = eig(A)
order = eigenvalues.argsort(0).swapaxes(0, 1)
Nlist = [[i] for i in xrange(N)]
return (eigenvalues[order, Nlist].swapaxes(0, 1),
R[:, order, Nlist].swapaxes(1, 2))
The above two functions take very little time compared with "eig" and
"inv". Using "sum" helps too. Given A, calculating "Abar = R E inv(R)"
can then simply be written
eigenvalues, R = sortedeig(A)
E = abs(eigenvalues) * numerix.identity(eigenvalues.shape[0])[..., np.newaxis]
Abar = mul(mul(R, E), inv(R))
Maybe there is a better way to calculate "Abar" rather than explicitly
calculating the eigenvalues and eigenvectors for every matrix. Any
help is much appreciated.
--
Daniel Wheeler
_______________________________________________
NumPy-Discussion mailing list
[email protected]
http://mail.scipy.org/mailman/listinfo/numpy-discussion
