How does MATLAB allocate memory while solving large (10m-ish) linear systems?

I am trying to solve a 'sparse' linear system of equations (Ax=b) with number of variables in the range of 1-10 million. Its a basic finite element simulation for solid mechanics problems, with a sparse and symmetric matrix (A), and sparse right side (b). The matrix is such that 0.01% of the total entries are non zero. I am using the '\' operator.

MATLAB was able to solve a 2.3m variable problem on my personal laptop (R2021b, Windows, 16GB RAM, i7-7700 2.8 GHz) in about 6 hours.

I have access to a (Linux) cluster that has a good number of processors and I can request any amount of memory (using PBS). However, when I try to run the same code, with 5m variables (a finer mesh), 300 GB of memory allocated, the job was quit right at the linear solution step, with memory error.

I need to know what I am doing wrong, and how can I fix this problem. 300 GB is outrageous, 16 GB is nothing. Is there a memory setting that MATLAB uses that is different for Linux and windows? I think MATLAB uses my hard drive to temporarily store data during the solution process, and it might not be the case with Linux. If that is the case, how do I switch that setting on inside my code.

I also welcome suggestions on different solvers or approaches I could use to solve problems such as these. Let me know if you guys need other information from me, and sorry if this has already been answered, just direct me there.

Thank you for your time!

Answer 1

[Matt J]

With 5m variables and a density of 0.01%, your matrix is going to consume about 56GB even before mldivide starts to do any work.

(5e6)^2*0.01/100*8*3/2^30

ans = 55.8794

You should check with the memory command how much memory is available to Matlab on your system (it wouldn't be the total 300Gb available to the system as a whole). However, it doesn't seem unreasonable that mldivide() could consume another 100 GB at least, and I can imagine that exceeding the amount of memory that Matlab has.

For a problem this size, you should probably be using an iterative solver, like pcg, rather than mldivide().

[Stephen23]

Links to some relevant docs:

https://www.mathworks.com/help/matlab/sparse-matrices.html (see section "Iterative Methods")

https://www.mathworks.com/help/matlab/math/iterative-methods-for-linear-systems.html

[Yash Agrawal]

As far as I know, there is no memory command on Linux, but I am pretty sure I am allocating 300 GB memory on the cluster to my job. How can I check the memory allocated to MATLAB in Linux?

Beside that, I was interested in knowing how is the 2.3m variable problem solved on my much inferior laptop (with 16 GB memory), and can I replicate those settings on the Linux cluster.

[Matt J]

Beside that, I was interested in knowing how is the 2.3m variable problem solved on my much inferior laptop (with 16 GB memory)

The RAM consumption of a 2.3m variable problem is much less. Repeating the above calculation:

(2.3e6)^2*0.01/100*8*3/2^30

ans = 11.8241

[Stephen23]

"With the memory command."

The memory documentation states "The memory function is available only on Microsoft® Windows® platforms."

[Matt J]

but I am pretty sure I am allocating 300 GB memory on the cluster to my job.

But the node you are working on does not have 300 GB, I suspect. I think you would need the MATLAB Parallel Server to make mldivide split its work across a cluster.

[John D'Errico]

I recall that the sparse solver in MLDIVIDE is not parallelized? If that is the case, then a parallel server would be of little help.

[Walter Roberson]

https://linuxhint.com/check_memory_usage_process_linux/

https://stackoverflow.com/questions/131303/how-can-i-measure-the-actual-memory-usage-of-an-application-or-process

[Yash Agrawal]

Thanks for the links and your help, guys! I will go with PCG.

[Yash Agrawal]

Hey guys, wanted to ask one more thing along the same line. How do we define that the matrix is symmetric in MATLAB, storing only one half of it. I have come across answers that might be outdated, saying that this functionality is not available. Is it available now?

[Matt J]

No, it is not.

[Walter Roberson]

I see to recall that someone posted a Symmetric Matrix class to File Exchange. I don't think it was sparse though.

In the case of dense (non-sparse) matrices, then squareform() converts back and forth between the two representations. However, there are hardly any other functions that work with that representation.

[Yash Agrawal]

I see, thanks guys! Have a wonderful time.

[Bruno Luong]

@Yash Agrawal "Hey guys, wanted to ask one more thing along the same line. How do we define that the matrix is symmetric in MATLAB, storing only one half of it. I have come across answers that might be outdated, saying that this functionality is not available. Is it available now?"

It seems CHOL use only upper-half of the input matrix.

A=rand(5)

A = 5×5

0.6082 0.1098 0.4983 0.2333 0.6290 0.1717 0.0780 0.9934 0.1774 0.1036 0.6496 0.0495 0.1380 0.4528 0.6525 0.9646 0.5738 0.0109 0.0844 0.9123 0.0663 0.7908 0.7017 0.5123 0.6362

A=triu(A'*A)

A = 5×5

1.7563 0.7182 0.6203 0.5818 1.7464 0 0.9752 0.7001 0.5154 1.1361 0 0 1.7465 0.7153 0.9627 0 0 0 0.5605 0.8634 0 0 0 0 2.0692

L=chol(A)

L = 5×5

1.3252 0.5420 0.4681 0.4390 1.3178 0 0.8255 0.5408 0.3361 0.5110 0 0 1.1113 0.2952 0.0625 0 0 0 0.4094 0.2312 0 0 0 0 0.1186

L'*L

ans = 5×5

1.7563 0.7182 0.6203 0.5818 1.7464 0.7182 0.9752 0.7001 0.5154 1.1361 0.6203 0.7001 1.7465 0.7153 0.9627 0.5818 0.5154 0.7153 0.5605 0.8634 1.7464 1.1361 0.9627 0.8634 2.0692

So you might use CHOL or DECOMPOSITION with 'upper' argument. However tat doesn't save at all memory for dense matrix. If your matruix is sparse then that might be a tip to save memory.

But as Matt's suggestion, for sparse matrix you should use iterative method. You can program your own matrix x vector function using only upper-storage.

[Yash Agrawal]

Thanks, Bruno!

I see, this is a nice way to save memory for sparse matrices when you want to go for Cholesky. What about sparse and symmetric matrices for the pcg algorithm (iterative solvers). Do you know any way to save memory over there.

[Bruno Luong]

As I just wrote above you use iterative solver by passing the function handle, for input x vector that provide

Ax := A*x

by this calculation:

Ax = (x'*U)'+U*x-diag(U).*x

where U is triu of A. So you need to store U = triu of A, and not A.

[Yash Agrawal]

I get that, thanks for helping me out!