-1

I'm currently working on many small 6x6 matrices: shape A = (N, N, N, 6, 6) with N is about 500. I store these matrices in a HDF5 file by Pytables (http://www.pytables.org).

I want to do some calculations on these matrices, say inverting, transposing, multiplication, etc... It's quite easy while N is not very big, by example numpy.linalg.inv(A) should do the trick without loop. But in my case, it works very slow and sometimes I have a memory's problem.

Could you suggest me an approach to do this more efficiently?

Chris Seymour
  • 83,387
  • 30
  • 160
  • 202
user2863620
  • 635
  • 3
  • 10
  • 17
  • 1
    This is not an answer, but if you are up on your mathematics, ["Don't invert that matrix"](http://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/) (including comments) might help. Also I assume you have seen it, but [this question](http://stackoverflow.com/questions/1053928/python-numpy-very-large-matrices) and particularly [this answer](http://stackoverflow.com/a/1054113/377366) seem to be in line with your problem. (Summary - use 64 bit os/python/numpy; be smarter with calculations and manipulation) – KobeJohn Mar 02 '14 at 15:30
  • Do a for-loop around the first dimension, so that you don't need to load the NxNxNx6x6 array to memory at once, but work e.g. on slices of 1xNxNx6x6. – pv. Mar 02 '14 at 18:57

1 Answers1

0

A 6x6 matrix has 36 8 byte double values, or 288 bytes; we'll say 0.5KB for ease of calculation and overhead. If you accept that, then 500 of them would only represent 250KB - not much memory.

If you keep all those inverses in memory it's still not a lot - just 500KB.

Can you calculate the amount of RAM your program is consuming to confirm this?

What are you doing - finite element analysis? Are these stiffness, mass, and damping matricies for 500 elements?

If yes, you should not be inverting element matricies. You have to assemble those into global matricies, which will consume more memory, and then solve that. Inverse still isn't calculated - LU decomposition in place is the usual way to do it.

I wouldn't consider a 500 element mesh to be a large problem. I saw meshes with thousands and ten of thousands of elements when I stopped doing that kind of work in 1995. I'm sure that hardware makes even larger problems possible today.

You're doing something else wrong or there are details that you aren't providing.

duffymo
  • 305,152
  • 44
  • 369
  • 561
  • I think you made a mistake: I have 500x500x500 of 6x6 matrices which are about 36 GB. – user2863620 Mar 02 '14 at 15:51
  • I was afraid of that. No, I didn't want to assume you were that silly. Do you have 36GB of RAM available to you? I'd question why you need that many matricies. You still haven't given more details about your problem and what you're doing. – duffymo Mar 02 '14 at 16:14
  • I don’t do FE but something similar. The inverting and transposing 6x6 matrix for each node is required. – user2863620 Mar 02 '14 at 16:30
  • So you'll agree that efficiency has nothing to do with it. You can't find X of anything into a Y container when X > Y. I'd start by writing each matrix out to disk as you create it. You should only bring them back into memory as you need them; no sense to keep them all. Otherwise you'd need > 36GB of RAM on a 64 bit machine, because the OS and VM take their share, too. – duffymo Mar 02 '14 at 16:37