# Reduce memory cost with three matrices multiplication

My program is eating up all my RAM and the operation is related to 3 matrices multiplication, as the code shown below.

example code:

{T, W, H} = {4, 10, 10}; (*crashed with larger values *)
matT = RandomReal[10, {T + 1, 2 T + 1}];
data = RandomReal[10, {W, H}];
matX = Table[RandomReal[i]*data, {i, 0, T}];
matX // Dimensions
matY = Table[RandomComplex[n+n I]*data, {n, -T, T}];
matY // Dimensions
Total[matX matT.matY]


After research I have detected the line of code which is killing a great deal of memory:

Total[matX matT.matY]


Where the dimension of matX is {T+1,W,H} , matY: {2T+1,W,H} and matT: {T+1,2T+1}, T is any integer but preferred even numbers.

edit:

one of the reason the RAM is being eaten up is due to the fact matY has imaginary part which is obtained from Fourier angular functions. The byte count is double size of real counterpart matX which only has real numbers.

This code works fine when {T,W,H}are small. However,my goal is to calculate the case when {T,W,H}={1000,1000,1000}, but now I only half way to it.

I guess MMA is copying a lot of data in the RAM to do this Total calculation and I think it can be improved in terms of reducing RAM usage. I would really appreciate your helps to give me some ideas.

• Where do the matrices matX, matY, and matT come from? Maybe that structure can be exploited. And what are the typical valued of T, W, and H that you want to use in practice? – Henrik Schumacher Oct 8 '19 at 12:27
• thx for reply. I have added more clarifications, I have tried to minimize the structruce of those matrix such as using symmetry etc. But at least I known {W,H } cannot be any less while I can split T dimensions into many parts. It's only solution I have now. It's a bit of tedious but it works anyway. – cj9435042 Oct 8 '19 at 12:59
• Btw., I suspect you want ot use RandomComplex[n + n I] and not RandomComplex[n + n i] – Henrik Schumacher Oct 8 '19 at 13:07

This is the actual data, everything else is derived from it:

{T, W, H} = {1000, 100, 100};
matT = RandomReal[10, {T + 1, 2 T + 1}];
data = RandomReal[10, {W, H}];
u = RandomReal[{0, 1}, T + 1];
v = Table[RandomComplex[n + n I], {n, -T, T}];


That is how you compute the result:

matX = TensorProduct[u, data];
matY = TensorProduct[v, data];
result = Total[matX matT.matY]; // MaxMemoryUsed // AbsoluteTiming


{1.67506, 512368624}

Here is how I would compute it:

result2 = u.matT.v data^2; // MaxMemoryUsed // AbsoluteTiming
Norm[result - result2, "Frobenius"]/Norm[result, "Frobenius"]


{0.001187, 240432}

The relative error is small:

Norm[result - result2, "Frobenius"]/Norm[result, "Frobenius"]


7.58384*10^-15

Now it is up to you to think deeply and to figure out why this produces the correct result more than a thousand times faster and with the two thousandth part of memory. After that, you won't make that mistake again.

;)

• thanks for this solution. I love it! But my challenge is how do I acquire u and v by given matX and matY, and data aren't same for each matrix. Is it possible to use "inverse" TensorProduct to get the u and v? Because I assume main point of your idea is that do not explicit manipulate whole matrices in memory, so I guess obtaining u and v is the key to collapse the dimensions during computations. – cj9435042 Oct 9 '19 at 12:46
• Well, I deduced u and v from the way you built matX and matY. If you build them withRandomReal[{0, 10}, {T + 1, W, H}], there is probably nothing that one can exploit... So, basically, what we exploit here is that Flatten/@matX and Flatten/@matY have a low rank decomposition. Sometimes such decompositions can computed efficiently by exploiting how exactly these tensors are construct. For completely random tensors, I don't know of a better method than singular value decomposition, but that is much too expensive. – Henrik Schumacher Oct 9 '19 at 13:21