Say I have a very large vector or matrix M and I want to do M*M’ to obtain say a time by time matrix. This doesn’t work for huge matrices. I guess the issue is memory? Are there any common ways to solve this?
The answer depends on how you need it,
- if its to build an operator better build a SciMLOperator.jl matrix free operator (M M^T x) and solve what you need with Kyrlov
- if you need the matrix you will have to store it in disk and play one chunk at a time ie Mmap.jl.
- if you need some elements of the matrix you can go with LazyArrays.jl
- If you want something more “high-lever” i guess Dagger.jl + DiskArrays.jl / DiskArrayEngine.jl could do the trick but its very young and I cant comment on it.
I think you need to step back for a moment and think of the bigger picture.
First, if your matrix M fits in memory but your matrix MM^T does not, I’m guessing that M is a “tall” matrix (rows ≫ columns). Even though the textbook formula for your answer may involve MM^T, one typically does not explicitly form MM^T in practice — not only is the latter matrix huge, but it also squares the condition number (which may double the number of digits you lose to roundoff error). Typically one uses either QR or SVD factorizations of M instead.
(You can also form MM^T implicitly as a linear operator as @yolhan_mannes suggested above, then use an iterative solver, but this doesn’t solve the conditioning problem … and ill-conditioning may cause such an iterative solver to converge slowly.)
More generally, if even M itself is too big to fit in memory, then you usually want to look for structure in where the matrix comes from to let you implicitly construct it in terms of fast matrix–vector products. If M truly comes from unstructured data and doesn’t fit in memory, then you may be stuck with expensive “out-of-core” algorithms that involve reading M by chunks, but even then one often tries to shrink the data required by e.g. random sampling / “sketching”.
What problem are you trying to solve with MM^T? That is, what are you trying to do with this matrix? And how big is M and where does it come from?
Thanks for the thorough comments! Yeah, it’s a “tall” matrix as you mention. I want to try out different approaches that may be similar to what recurrence quantification analysis does (not necessarily the embedded delay addition part though), including approaches that probably should not be called like that because they don’t dichotomize the resulting matrices. And try different ways to construct the time by time matrix (self-, cross- similarity/distance). The problem is (lag-)synchrony estimation. M does fit in memory.
The point is that any algorithm that relies on constructing the whole Gram matrix G = MM^T will not scale, and you want to look for ways to skip over this to compute your final desired result implicitly without G. If you’re ultimately going to threshold G or something, for example, you might look into tree-based nearest-neighbor searches.