0

Is there any obvious way in numpy to replace something like:

for x in X:
    xi, xj = meshgrid(x, x, indexing='ij')

with a single (and possibly more efficient) operation like:

Xi, Xj = multi_meshgrid(X, X, indexing='ij')

The example of X is the following:

X = np.array([[0,1,2,3,4,5], [5,6,7,8,9,10], [11,12,13,14,15], ...])

The main problem is that i can have tens and hundreds of thousands of entries in X and possibly operation is often repeated.

The problem arises from assembling global stiffness matrix K in finite element method. For each entry in X of length n i have a matrix "n x n" which i have to inscribe into this global matrix. This matrix is in scipy.sparse coordinate format.

Regards, Marek

mrkwjc
  • 1,080
  • 8
  • 17
  • 1
    Can you give some context about what you are trying to achieve? The answer is often "this can be done with broadcasting". An input example for `X` would also be helpful – YXD Jun 03 '14 at 16:48
  • Also a related but not identical question: http://stackoverflow.com/questions/1827489/numpy-meshgrid-in-3d – YXD Jun 04 '14 at 10:39

1 Answers1

2

I think this answers the question, though I'm not sure if this is the best for constructing the sparse matrix in the end.. Anyway the following code creates a "view" into X, so it's very efficient both computationally and memory-wise. Try it :)

from numpy.lib.stride_tricks import as_strided

m = 3
n = 4
X = np.arange(m*n).reshape((m,n))

sz = X.itemsize
Xi = as_strided(X, shape=(m,n,n), strides=(n*sz, sz, 0))
Xj = as_strided(X, shape=(m,n,n), strides=(n*sz, 0, sz))

This does, however, not work when X is not a regular matrix. E.g. in your example the third row has 5 elements whereas the others have 6.

  • This works great, though it's not so obvious. :) I'm about 50x faster at this step. Thanks! – mrkwjc Jun 04 '14 at 11:52