# Matrix Vector multiplication for very big matrices (>>RAM)

Problem

Suppose a matrix-vector product of the form $M\vec{v}$ should be calculated where the amount of storage needed to store $M$ is noticeably larger than the available RAM on a machine. What is the fastest way to perform such an operation?

Slow solution approach

Since Mathematica takes care of such problems by the provision of InputStream object I thought of using them and implemented a version where the matrix $M$ is given as (very long) stream. I Dot it row-wise with the vector $\vec{v}$ and collect the results as final result vector

Table[Dot[NextStreamEntries[matrixStr, VectorLength], v], {ii,VectorLength}].


Here matrixStr is the matrix as an InputStream and VectorLength is the Lengthof the vector $\vec{v}$. By the way, I use a Table to do exactly the same thing VectorLength times. This seems weird but I did not find another fast (!) solution.

The performance of this is a catastrophy. Using the RuntimeToolsProfile I found that the performance is worse because of this method

NextStreamEntries[stream_, count_] := BinaryReadList[stream, "Real64", count];


with which I travel through the rows of the matrix.

On request: What does the matrix store and where does it come from?

The matrix comes from outside MMA since performance matters. A C++ routine writes it to disk as a binary file and I bit by bit read in in from MMA with BinaryReadList. The matrix elements are quantum mechanical expectation values (needed in context of the iterated equations of motion approach to describe non-equilibrium physics) with very few elements equal to zero.

• ParallelMap[#.v &, m] Aug 19, 2017 at 12:40
• @Edmund Thanks for your answer. Can you elaborate on it in detail? The problem I see is that the stream is not guaranteed to be read in the correct order which will lead to a wrong result. Maybe I did not get your point.
– pbx
Aug 19, 2017 at 12:43
• can you tell me how many elements has he matrix? Aug 19, 2017 at 12:51
• @Alucard It depends on many different factors. A typical example is a quadratic matrix of $70k$ entries of floating numbers that leads to $64*(70000^2) \mathrm{bit} \approx 39.2 \mathrm{GB}$. And this doesn't fit into (at least) my memory.
– pbx
Aug 19, 2017 at 13:01
• @pbx Okay. That's quite likely too many nonzeros to be exploited. Still, I am curious: Where does this matrix come from? How is it produced? Do you need many matrix-vector multiplications or only a few? How expensive is the creation of the matrix rows compared to saving them on HD and even piping them through string expressions? A binary format for storage should be much faster since the bottleneck here should be the bandwidth of your mass storage... Aug 19, 2017 at 17:17

You may use ParallelTable and other functions from the Parallel Computation Setup & Configuration guide.

With

m = LowerTriangularize@ConstantArray[1, {10, 10}];
v = ConstantArray[1, 10];
MatrixForm /@ {m, v} This m and v is selected to demonstrate that order is preserved.

Then create the stream.

stream = StringToStream@ExportString[m, "Table"];


and setup for parallel evaluation

LaunchKernels[];
SetSharedVariable[stream];
sig = ConstantArray[Real, Length@v];
UnsetShared[stream];
MatrixForm@b Clean-up.

Close[stream]


As you can see the order is preserved as the stream is shared over the kernels. You can use Number instead of Real if you have integer data.

I am getting a warning that the stream is not open but it is and the code does execute in parallel. I will look into that. You can silence the warning with Quiet@Read.

Hope this helps.

• This is highly interesting. Up to now I did not even know about SetSharedVariable which does the magic here undoubtedly. +1 and of course an answer acceptance. Thanks. :)
– pbx
Aug 19, 2017 at 13:14
• @pbx Just added an update so that the Read signature is not created on each read. Should be ever so slightly faster. May result in a big difference in your case. Aug 19, 2017 at 13:17
• @pbx Since your matrix is so huge you can look into modifying this answer (19542) for ParallelTable` so that you can track the progress as it evaluates. Aug 19, 2017 at 13:33