Vectorize eigenvalue calculation in Numpy

Question

I would like a numpy-sh way of vectorizing the calculation of eigenvalues, such that I can feed it a matrix of matrices and it would return a matrix of the respective eigenvalues.

For example, in the code below, B is the block 6x6 matrix composed of 4 copies of the 3x3 matrix A. C is what I would like to see as output, i.e. an array of dimension (2,2,3) (because A has 3 eigenvalues).

This is of course a very simplified example, in the general case the matrices A can have any size (although they are still square), and the matrix B is not necessarily formed of copies of A, but different A1, A2, etc (all of same size but containing different elements).

import numpy as np
A = np.array([[0, 1, 0],
              [0, 2, 0],
              [0, 0, 3]])
B = np.bmat([[A, A], [A,A]])
C = np.array([[np.linalg.eigvals(B[0:3,0:3]),np.linalg.eigvals(B[0:3,3:6])],
              [np.linalg.eigvals(B[3:6,0:3]),np.linalg.eigvals(B[3:6,3:6])]])

Answer 1

Edit: if you're using a version of numpy >= 1.8.0, then np.linalg.eigvals operates over the last two dimensions of whatever array you hand it, so if you reshape your input to an (n_subarrays, nrows, ncols) array you'll only have to call eigvals once:

import numpy as np

A = np.array([[0, 1, 0],
              [0, 2, 0],
              [0, 0, 3]])

# the input needs to be an array, since matrices can only be 2D.
B = np.repeat(A[np.newaxis,...], 4, 0)

# for arbitrary input arrays you could do something like:
# B = np.vstack(a[np.newaxis,...] for a in input_arrays)
# but for this to work it will be necessary for each element in 
# 'input_arrays' to have the same shape

# eigvals will operate over the last two dimensions of the array and return
# a (4, 3) array of eigenvalues
C = np.linalg.eigvals(B)

# reshape this output so that it matches your original example
C.shape = (2, 2, 3)

If your input arrays don't all have the same dimensions, e.g. input_arrays[0].shape == (2, 2), input_arrays[1].shape == (3, 3) etc. then you could only vectorize this calculation across subsets with matching dimensions.

If you're using an older version of numpy then unfortunately I don't think there's any way to vectorize the calculation of the eigenvalues over multiple input arrays - you'll just have to loop over your inputs in Python instead.

Answer 2

You could just do something like this

C = np.array([[np.linalg.eigvals(B[i:i+3, j:j+3])
              for i in xrange(0, B.shape[0], 3)]
                for j in xrange(0, B.shape[1], 3)])

Perhaps a nicer approach is to use the block_view function from https://stackoverflow.com/a/5078155/1352250:

B_blocks = block_view(B)
C = np.array([[np.linalg.eigvals(m) for m in v] for v in B_blocks])

Update

As ali_m points out, this method is a form of syntactic sugar that will not reduce the overhead incurred from calling eigvals a large number of times. While this overhead should be small if each matrix it is applied to is large-ish, for the 6x6 matrices that the OP is interested in, it is not trivial (see the comments below; according to ali_m, there might be a factor of three difference between the version I give above, and the version he posted that uses Numpy >= 1.8.0).