# Memory/speed tradeoff for many small matrix inverses

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.

• When you say it's the memory bottleneck, do you mean bandwidth/speed or capacity/size? If it's bandwidth, moving it to disk is going to make things worse. If it's a size thing.. are you sure this is really an issue? Storing a single matrix per element, of fixed size, is not really asymptotically worse than storing the state vector/iterate vector itself (though certainly worse by a large constant factor). If storing them really is using too much memory, I think you would be better served by looking at distributed memory parallelism than out-of-core techniques. – rchilton1980 Sep 18 '18 at 14:33
• @rchilton1980 I mean capacity/size. As for the matrix vs. vector, yes, the large constant factor is a problem, and only gets worse with higher polynomial order. I have written a distributed version of the solver that splits the problem between nodes. But on each node, the storage of these inverses or their computation is still the limiting factor. Of course, the distributed approach also introduces a communication overhead between nodes. – user3482876 Sep 18 '18 at 14:48
• It might end up happening that doing the same inverses on the fly is faster than storing them. One keyword you might be looking for regarding the out-of-core solver is io latency hiding. Its efficiency will depend on how your other operations during the iterative solve allow you to incorporate the disk IO without significant slowdown. – Anton Menshov Sep 18 '18 at 15:39
• Any repetition in the grid/elements that could be exploited? Can you maybe use a mostly-structured grid? – rchilton1980 Sep 18 '18 at 17:04
• @rchilton1980 in general, no; the large problems in which I am interested involve unstructured grids. It's a good suggestion, though. – user3482876 Sep 18 '18 at 18:38

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.