[SciPy-user] optimal way to solve large amount of sets of linear equations

Lubos Vrbka lists@vrbka....
Wed Jun 11 10:15:45 CDT 2008

hi guys,

i have to solve a set of equations, that can be in matrix notation 
written as
   {E-C}H = C
where E is identity matrix. actually as a result, i want to have 
function G, defined as H - C, i.e.
   G = H - C = {E-C}^(-1)C - C
the problem is, that each of the matrices H, C, G represent set of 
discretized functions. then, i have to solve this problem for every step
   G(1) = {E-C(1)}^(-1)C(1) - C(1)
   G(2) = {E-C(2)}^(-1)C(2) - C(2)
the matrices themselves are usually small (2x2, 3x3), but the number 
these equations is quite large (1024, 2048, 4096, or 8192; usually not 

i tried to use the following pieces of code, but don't know (since i'm a 
newbie to scipy), whether it cannot be done in a better or (less stupid) 

i'd be really very grateful for any comments.

G_f_ij is scipy.array((n,n,npoints))
E_ij is scipy.eye(npoints)
c_f_ij is scipy.array((n,n,npoints))

the slow way:
for i in range(npoints):
   # solve the matrix problem for all dr
   G_f_ij[:,:,i] = scipy.mat((E_ij - c_f_ij[:,:,i])).I * 
scipy.mat(c_f_ij[:,:,i]) - scipy.mat(c_f_ij[:,:,i])

the faster way:
for i in range(npoints):
   G_f_ij[:,:,dr] = numpy.linalg.solve(scipy.mat((E_ij - 
c_f_ij[:,:,dr])), scipy.mat(c_f_ij[:,:,dr])) - scipy.mat(c_f_ij[:,:,dr])

the fastest way:
forget the scipy.mat as it's not needed for linalg
for i in range(npoints):
   G_f_ij[:,:,dr] = numpy.linalg.solve(E_ij - c_f_ij[:,:,dr], 
c_f_ij[:,:,dr]) - c_f_ij[:,:,dr]

the last snippet of code is roughly 2x faster than the first one. for 
2x2x1024 points it takes 0.2, 1.3, 0.7 seconds on my computer. alone, 
that's not that bad - but when i need to repeat that 100-1000x during 
the execution of the program, then it starts to take a lot of time...

best regards,

Lubos _@_"

More information about the SciPy-user mailing list