Memory/speed tradeoff for many small matrix inverses

Question

Problem

In the case of a finite element code, I have many small (order of 30x30) matrix inverses (or LU factorizations), one per finite element. These matrix inverses never change and must be applied to vectors repeatedly in the context of a matrix-free solver.

Storing all these inverses is the memory bottleneck of the code, because they have to be applied for each finite element in the mesh. On the other hand, computing the inverses 'on the fly' (using LAPACK or Numpy) is not practical, because these inverses must be applied at each iteration of a sparse iterative solve, so each time-step would require

(number of CG iterations) * (number of elements) * (time of inverse)

which is too expensive, even for problems far smaller than those I am interested in.

(possible) solution

I am considering performing the matrix inversions or LU factorizations once, and then storing them to disk, rather than keeping them in memory. Then the idea is to read the inverses into memory either one at a time or by chunks, and apply them, avoiding the repeated inversions. But I know that reading from disk can be expensive.

Are there good or accepted ways of doing this? Most out-of-core papers I've looked at are considering direct solves or much more complicated problems, so I would appreciate any pointers to relevant references as well.

EDIT: The code is a mix of Python and C, and I've found some posts indicating that perhaps the HDF5 library is a good alternative to rolling my own cache, as suggested above.

Answer 1

There is no faster way to storing things than to put them into memory. If your problem is such that you can't store all of the element matrix inverses in memory, then you have two options: (i) buy more memory for the machines you use, (ii) use more machines. There is also (iii) solve smaller problems. Trying to externalize the matrices to disk will make things vastly slower because disks are at least an order of magnitude slower than memory. Disks are also vastly slower than communication between machines.

But it's worth stepping back. If your solver scheme requires storing lots of local matrix inverses, then maybe that's not a great solver. I'm making the assumption here that in your scheme, storing the inverses is more expensive than storing the local matrices to begin with -- for example, because the local matrices are sparse but the inverse are not. Generally, my experience from 20 years of finite element computations is that if you are running out of memory, then you've already run out of CPU time as well. So I'm a bit surprised that runtimes don't also show up on the list of your worries.

Answer 2

One idea would be to store the low-rank approximation of the inverse of the matrix $A^{-1}\approx UV$ using SVD. If $A\in \mathbb{R}^{30\times 30}$, then $U\in \mathbb{R}^{30\times k}\ \& \ V\in \mathcal{R}^{k\times 30} $, where $k\approx 4,5$ is the rank. Use the $UV$ approximation as a preconditioner for solving the system $Ax=b$ using GMRES. You can use vary $k$ as a tuning parameter that trades off memory with the convergence rate.