On Mon, Jul 11, 2011 at 05:01:07PM -0400, Daniel Wheeler wrote:
> 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.
If you really care about speed, for matrices of this size you shouldn't
call a linear algebra pack, but simply precompute the close-form
solution. Here is sympy code that I used a while ago to generate fast
code to do inverse of SPD matrices.
G
"""
Generate the code for fast_inv
"""
from itertools import chain
from sympy import Symbol, cse, Matrix
from sympy.printing.str import StrPrinter, precedence, S
# First create a special to print squares and multiplications, rather
# than pows, as cython is not able to optimize them
class MyPrinter(StrPrinter):
def _print_Pow(self, expr):
PREC = precedence(expr)
if expr.exp is S.NegativeOne:
return '1/%s'%(self.parenthesize(expr.base, PREC))
elif expr.exp == 2:
return '%s*%s'%(self.parenthesize(expr.base, PREC),
self.parenthesize(expr.base, PREC))
else:
return '%s**%s'%(self.parenthesize(expr.base, PREC),
self.parenthesize(expr.exp, PREC))
my_printer = MyPrinter()
def code(n):
""" Generate the code for inversing a C matrix of rank n.
"""
l = list()
for i in range(n):
ll = list()
for j in range(n):
if j < i:
ll.append(l[j][i])
else:
ll.append(Symbol('a[%i, %i]' % (i, j)))
l.append(ll)
m = Matrix(l)
a = m.inverse_LU()
vars, expr = cse([poly#.radsimp()
for poly in chain(*a.tolist())])
var_str = '\n '.join(['cdef double %s = %s' %
(name, my_printer.doprint(v)) for name, v in vars])
out_str = """
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.profile(False)
@cython.cdivision(True)
@cython.cdivision_warnings(False)
cdef _inv_%i(np.ndarray[np.float_t, ndim=2, mode="c"] a):
# Intermediate variable definition
%s
# Final matrix
return np.array(%s).reshape((%i, %i))
""" % (n, var_str, my_printer.doprint(expr), n, n)
return out_str
out_file = file('fast_inv.pyx', 'w')
out_file.write("""
'''
Code to do a fast inverse on small symmetric matrices
'''
cimport cython
from scipy import linalg
import numpy as np
cimport numpy as np
from fast_inv import inv
""")
N = 7
for n in range(2, N+1):
out_file.write(code(n))
out_file.flush()
print 'Done order %i' % n
out_file.write('''
def sym_inv(a):
""" Return the inverse of the symetric matrix 'a'.
This code is optimized for small symmetric matrices.
"""
cdef int n = a.shape[0]
if n == 1:
return 1./a
%s
else:
return inv(a)
''' % ''.join(["""
elif n == %i:
return _inv_%i(np.ascontiguousarray(a))"""
% (i, i) for i in range(2, N+1)])
)
out_file.close()
_______________________________________________
NumPy-Discussion mailing list
[email protected]
http://mail.scipy.org/mailman/listinfo/numpy-discussion
