6

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])]])
6502
  • 112,025
  • 15
  • 165
  • 265
Andrei Berceanu
  • 351
  • 2
  • 4
  • 10

2 Answers2

5

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.

ali_m
  • 71,714
  • 23
  • 223
  • 298
  • 1
    Which edition of numpy are you using? `C = np.linalg.eigvals(B)` gives `LinAlgError: 3-dimensional array given. Array must be two-dimensional `. `eigval` gives values for one square matrix. – hpaulj Oct 19 '13 at 19:32
  • You're right - it must be a relatively new feature. I'm using `1.9.0.dev-f665c61` which works, but `1.7.1` fails with that `LinAlgError` exception. – ali_m Oct 19 '13 at 19:37
  • 1
    That feature was added in `1.8.0rc2`. Here's the commit: https://github.com/numpy/numpy/commit/15a9c3b25c0aff2799c927ea9e602fc3c134d3fc – ali_m Oct 19 '13 at 19:40
  • so, any solution for numpy 1.7? – Andrei Berceanu Oct 19 '13 at 19:49
  • @AndreiBerceanu Unfortunately I don't think there is any way to vectorize the eigenvalue decomposition in 1.7. I think you'll have to either upgrade to a newer version of numpy or else just live with looping over your input arrays and calling `eigvals` on each one. How many input arrays do you have, and how big are they? – ali_m Oct 19 '13 at 19:55
  • If all of the sub matrices are 3x3 (or 2x2) you could vectorize the algebraic solution. – hpaulj Oct 19 '13 at 20:04
  • my matrix "B" has dimensions (6,6,512,512), i.e. is formed by 512**2 6x6 matrices that i want to diagonalize. – Andrei Berceanu Oct 19 '13 at 20:07
  • as for upgrading numpy, i would gladly do it, but was unable to find any ubuntu binaries for v1.8 or realiable instructions on how to do it from source.. – Andrei Berceanu Oct 19 '13 at 20:08
  • 3
    The overhead of 36 iterations, compared to the actual processing time of the 512x512 matrices is going to be negligible: use a couple of nested for loops, or use `np.vectorize`. If you had 512x512x6x6, then you would want to switch to 1.8. – Jaime Oct 19 '13 at 20:10
  • 1
    but this is what i have, in fact. A is 6x6 and B is 512x512xA, so i need to call np.eigenvalues 512**2 times. – Andrei Berceanu Oct 19 '13 at 20:28
  • The main difficulty with installing numpy on ubuntu is making sure you have all of the dependencies installed. I already installed 1.7, so upgrading to 1.9 was straight forward. – hpaulj Oct 20 '13 at 00:42
1

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).

Community
  • 1
  • 1
Legendre17
  • 621
  • 5
  • 13
  • that's syntactically a bit nicer, but you're not really vectorizing the eigenvalue decomposition - you still call `np.linalg.eigvals` separately on each submatrix – ali_m Oct 19 '13 at 20:19
  • I'm not sure what you mean. I don't think there's any other way to get the eigenvalues of each block of a matrix than to actually calculate them for each block -- whether this is done under the hood by some built-in function or not. I don't think such a Numpy built-in function exists, so the above is a way to implement it. Note that there isn't really any performance issue, since slicing and `block_view`'s stride tricks don't deep copy the submatrices. – Legendre17 Oct 19 '13 at 20:54
  • 1
    `np.linalg.eigvals` is a wrapper that calls the `dgeev` LAPACK function internally, which indeed operates on `(N, N)` 2D slices of the input array. However, there is a substantial performance benefit to looping over `M` slices in C, which is what happens when you call on an `(M ,N, N)` array (at least using numpy >= 1.8), versus using Python loops like in your example. 'Vectorization' usually only refers to cases where the loops are pushed down to a low-level language such as C or Fortran rather than in a high-level language such as Python where they carry much more overhead. – ali_m Oct 19 '13 at 21:16
  • Is this true? How big is the overhead? I don't have a new-enough version of Numpy to test this. – Legendre17 Oct 19 '13 at 21:24
  • 1
    The significance of that overhead will depend on the ratio between how many matrices you're looping over and how large they are. If you have a small number of large matrices then a greater fraction of your time will be spent computing eigenvalues rather than looping, so the benefits of vectorization will be negligible. However, if you have a large number of small matrices then there may be a more significant advantage - for a `(500, 6, 6)` random array it takes roughly twice as long looping over each `(6, 6)` submatrix and calling `eigvals` on it versus calling `eigvals` on the whole array. – ali_m Oct 19 '13 at 21:40
  • Yeah, that's what I meant -- I was wondering how much the overhead matters for the 6x6 matrices that the OP is interested in. Of course it depends on the size. A factor of two is certainly significant, but if updating to a newer Numpy isn't possible, I don't think there's an alternative (short of writing your own C module). – Legendre17 Oct 19 '13 at 22:02
  • It's actually more like a factor of 3 for a `(512**2, 6, 6)` array like the OP has. – ali_m Oct 19 '13 at 22:17