[SciPy-User] Vectorize matrix inversion

[Jonathan Stickel]

On 5/8/11 00:11 , scipy-user-request at scipy.org wrote:
> Date: Sat, 7 May 2011 16:15:54 -0700 (PDT)
> From:"davide"
> Subject: [SciPy-User] Vectorize matrix inversion
>
> Hey all,
>
> I have a NxNxP array. Let's call it foo. foo[:,:,p] contains a
> matrix
> I want to invert. Is there a pure python/scipy way to compute an array
> bar without loops such that it
> would be equivalent to the following?
>
> import scipy.linalg as la
> import numpy as np
>
> bar = np.zeros((N, N, P)) + 0j
>
> for i in range(0,P):
>          bar[:,:,i] = la.inv(foo[:,:,i])
>
> I realize I could write some fortran code to do this or even use
> cython, but it would be nice if I could do this without needing to
> compile some extra code. As a summary, does anyone know how to
> compute
> the above (either example, but preferably both) without using loops?
>
> Cheers,
>
> Davide


Recently I wanted to invert a block diagonal matrix where each block was 
a vandermonde matrix.  I found it reasonably fast to put each block in a 
list and then use list comprehension to invert each block:

# form the block matrices
B = [np.vander(x[j:j+nn],nn) for j in xrange(N-n)]
# compute the inverse for each block matrix
Binv = [scipy.linalg.inv(B[j]) for j in xrange(N-n)]

Note exactly "vectorized", but faster than an explicit for loop.

HTH,
Jonathan