# Memory usage for large matrices

I have a quick and probably naive question about how MMA (11.2) handles memory (I am very novice with numerical calculation using MMA..).

I have to work with large matrices (of sizes over $30000^2$). I obtain two symmetric matrices, say $A$ and $B$, from two lower triangular ones ($At$ and $Bt$), and I compute the eigenvalues of $A.B$, i.e.

A=At+Transpose[At]-DiagonalMatrix[Flatten[MapIndexed[#1[[#2]]&,At]]];
B=Bt+Transpose[Bt]-DiagonalMatrix[Flatten[MapIndexed[#1[[#2]]&,Bt]]];
Clear[At,Bt];

P=A.B;
Ev=Eigenvalues[P];
Clear[P];


I would like to know if stocking the symmetric matrices in "variables" $A$ and $B$ and their product as $P$ uses up more memory than the code below:

Ev=Eigenvalues[(At+Transpose[At]-DiagonalMatrix[Flatten[MapIndexed[#1[[#2]]&,At]]]).(Bt+Transpose[Bt]-DiagonalMatrix[Flatten[MapIndexed[#1[[#2]]&,Bt]]])];


I would say that the first code is more memory demanding than the second one. (BTW, is it OK to compute the eigenvalues of a sum of matrices directly? error-wise and time-wise I mean)

A subsidiary question: is the way I construct my symmetric matrices the fastest one?

• Do these matrices happen to be sparse? In that case, SparseArray would be more memory efficient. – Henrik Schumacher Jul 17 '18 at 6:20

Let me try to answer to the best of my understanding. First of all the function MemoryInUse reflects the amount of allocated kernel memory quite sensibly (I swapped identifiers by mistake in the test below, my Ad and A are your respective A and At), and the MaxMemoryUsed reports the maximum amount of memory used during the evaluation of its argument:

In[1]:= 30000^2*8

Out[1]= 7,200,000,000

In[2]:= A = RandomReal[{0, 1}, {30000, 30000}]; // MaxMemoryUsed // AbsoluteTiming

Out[2]= {7.3374, 7,200,000,496}

In[3]:= ByteCount[A]

Out[3]= 7,200,000,152

In[4]:= MemoryInUse[]

Out[4]= 7,252,233,296

In[5]:= B = RandomReal[{0, 1}, {30000, 30000}]; // MaxMemoryUsed // AbsoluteTiming
MemoryInUse[]

Out[5]= {7.50039, 7200000288}

Out[6]= 14,452,239,016


So the memory grows by the expected byte size of the two allocated matrices.

Note that unlike Timing or AbsoluteTiming, MaxMemoryUsed apparently has no provision for returning the result of the evaluation, so the evaluation has to have a side effect, practically.

Next of all, you can speed the construction of the diagonal matrix a little by using MMA's own Diagonal instead of your iterate-then-flatten approach:

In[7]:= Ad = A + Transpose[A] - DiagonalMatrix[Diagonal[A]]; // MaxMemoryUsed // AbsoluteTiming
Bd = B + Transpose[B] - DiagonalMatrix[Flatten[MapIndexed[#1[[#2]] &, B]]]; // MaxMemoryUsed // AbsoluteTiming

Out[7]= {21.4354, 21,600,001,760}

Out[8]= {26.9613, 21,600,000,888}


Note a missing possible optimization here. It seems that during each evaluation there was an allocation for two extra temporary matrices besides the one returned. Only one allocation instead of three would be required in e. g. a C++ program (Ad allocated and set to $-diagonal$, then A added in place to it twice: iterated once row first, and the second time column first for transposing). I do not really know if this is possible in MMA.

Clearing variables releases the memory immediately. This may not be necessarily true for smaller objects, but clearing the large matrices has an immediate effect:

In[9]:= Clear[Ad, Bd] // MaxMemoryUsed
MemoryInUse[]

Out[9]= 0

Out[10]= 14,452,252,200


Note that MaxMemoryUsed does not return a negative number when releasing memory.

In general, I think that the way you do it is likely the fastest. Modern CPUs are very good at vectorized operations, and MMA uses them quite extensively. This is why Diagonal is so much faster than your iterating function. A separate timing run shows this drastic speed difference in computing the argument to the last term using the two approaches:

In[11]:= Bdiag = Flatten[MapIndexed[#1[[#2]] &, B]]; // AbsoluteTiming

Out[11]= {4.72654, Null}

Out[12]= {0.0108354, Null}


Clearing variables releases memory, but MaxMemoryUsed without arguments reports the maximum memory usage since the start of the kernel session

In[13]:= Clear[A, B, Adiag, Bdiag];
MemoryInUse[]
MaxMemoryUsed[]

Out[14]= 52,262,896

Out[15]= {0, 0, 0, 0, 0, 0}

Out[16]= 43,252,250,256


Previous version of this answer mentioned that not all memory is always freed. I do not know why and when this happens. This may be a bug or a feature, but while running the updated script I do not see the memory left behind allocated. Out[14] above is only slightly higher than the starting kernel memory use.

To answer your last question, you should find it out by wrapping the evaluations in MaxMemoryUsed. I just did not dare to try solving the eigenproblem on a random 30K² matrix.

I will not be surprised if MMA uses some memory optimization, avoiding full allocation for partial expressions. From my experience with large matrices and Intel MKL (large in the same sense yours are medium-sized), eigenproblems and decompositions usually require the full matrix present in memory, but the dot product can be computed band-by-band, so I won't be surprised if MMA does not fully store the intermediate arguments to the Dot function. On the other hand, the result from evaluating In[7] above shows that obvious optimizations may be missed. In any case, wrapping your evaluation in MaxMemoryUsed will give the definitive answer. Often, splitting and/or re-arranging calculations in an expression results in different timing (and probably also different memory use, but in my experience I mostly cared about run time, so cannot tell much here).

Thanks to @b3m2a1 for the suggestions on making the answer more concise and @Mr.Wizard for improving the formatting.

I am pasting below only the input expressions from all of the above, so that you can repeat the tests quickly on your machine. I am using version 11.3 on Windows, for comparison.

30000^2*8

A = RandomReal[{0, 1}, {30000, 30000}]; // MaxMemoryUsed // AbsoluteTiming

ByteCount[A]

MemoryInUse[]

B = RandomReal[{0, 1}, {30000, 30000}]; // MaxMemoryUsed // AbsoluteTiming
MemoryInUse[]

Ad = A + Transpose[A] - DiagonalMatrix[Diagonal[A]]; // MaxMemoryUsed // AbsoluteTiming
Bd = B + Transpose[B] - DiagonalMatrix[Flatten[MapIndexed[#1[[#2]] &, B]]]; // MaxMemoryUsed // AbsoluteTiming

MemoryInUse[]

Bdiag = Flatten[MapIndexed[#1[[#2]] &, B]]; // AbsoluteTiming

• MaxMemoryUsed can also take an expression as an argument and tell you how much memory was used in the calculation, e.g. MaxMemoryUsed[Bd = B + ...]. I find that can sometimes be more helpful. – b3m2a1 Jul 17 '18 at 5:01