Running out of memory selecting parts of a large sparse matrix

I have a sparse matrix of around 400k x 400k rows/columns containing similarities (symmetric, sparse, double values between 0 and 1, diagonal is always 1). I'm hoping to eventually be able to work with 1000k x 1000 k matrices.

I need to select all entries above a given threshold. My problem is that I keep running out of memory (I have 64 GB available).

I found that I get the best results so far if I do not set values below the threshold to zero, but rather if I select all values at or above the threshold and build a new sparse matrix with those. However with matrix sizes increasing, I'm hitting the memory limit.

Can you point me to a way to reduce memory consumption for this process?

With the help of this forum, I tried several ways:

% Get relevant values from matrix
rel_values = matrix(matrix>threshold);

% Get index to relevant values
[i,j] = find(matrix>threshold);

% Form sparse matrix
matrix = sparse(i,j,rel_values,columns_mat,columns_mat);

I have run into limits both with matrix(matrix>threshold) and with find(matrix>threshold).

To avoid matrix(matrix>threshold) I tried:

% get relevant positions
rel_values_pos = matrix>threshold;

% Get index to relevant values
[i,j] = find(rel_values_pos);

% Get relevant values (i.e. values over threshold)
idx = sub2ind(size(matrix),i,j);
rel_values = matrix(idx);

% Form sparse matrix
matrix = sparse(i,j,rel_values,columns_mat,columns_mat);

This doesn't work either, both first lines have run into memory shortages.

My hope was to use the symmetry and only consider the following:

matrix = triu(matrix,1); % I know the diagonal to be one, so I can neglect those entries and append them later

However, triu(matrix,1) runs out of memory itself.

Using the symmetry seems to me the most potent approach, as it reduces the amount of entries by more than half (considering that I drop the diagonal), however I'm unaware of a smooth way to select the upper right triangle without triu.

If there is an aspect I'm missing entirely, I'm thankful for those hints as well of course.

[James Tursa]

Do you want the result to be the same size as the original, but only contain the values of interest?

[Benjamin Lender]

Yes that's the goal. The original problem was "Set all values below threshold to zero", but manipulating the matrix was more memory intense than selecting the values of interest and build a new one.

[David Goodmanson]

Hi Benjamin,

how many nonzero values are there typically?

[Benjamin Lender]

800 entries: 98,44 % zeros

8000 entries: 99,85 % zeros

80000 entries: 99,94 % zeros

[Bjorn Gustavsson]

I think you might get a problem that might be avoidable by "set all non-zero values below the threshold to zero" instead of "set all the values below threshold to zero".

[James Tursa]

A mex routine can generate this without needing any large intermediate variables. I will try to find some time to write this tonight. That being said, even a very sparse 1000K x 1000K matrix is still going to need an internal 1000K integer array to hold all of the column indexing data.

What version of MATLAB are you using?

[Benjamin Lender]

Please excuse my absence, I'm running 2020b

[James Tursa]

I have some prototype code written but it is untested. Specically I don't have R2020b at the moment to test with, but I will see what I can do ...

[Benjamin Lender]

@James Tursa I don't know what a mex routine is, but not needing any large intermediate variables sounds like something that would help a lot. I'm trying to avoid those at every step, I hope I can do better yet.

[James Tursa]

A mex routine is compiled C/C++ or Fortran code that, once compiled, can be called like a normal MATLAB function for inputs/outputs. If I can find a way to test it, I would be supplying you with C code. You would need to install a supported C compiler and compile the code at your end. E.g., see this link:

https://www.mathworks.com/help/matlab/matlab_external/install-mingw-support-package.html?searchHighlight=mingw&s_tid=srchtitle

[Benjamin Lender]

Hello world,

is there any new developmet, is there additional information I could provide that would help? @James Tursa a routine that would not necessitate large intermediate variables sounds really really good, nevermind that I hadn't heard of mex routines at all previously.

Answer 1

[David Goodmanson]

Hi Benjamin,

I believe this works as required. This example is 4e5 x 4e5, with 99.99% zeros. On my pc it takes 3 sec to make the initial sparse matrix and 1 sec to make the new one.

tic
n = 4e5;
fracnz = 1e-4;
m = round(n^2*fracnz);
r = randi(n,m,1);
c = randi(n,m,1);
p = rand(m,1);
s = sparse(r,c,p,n,n);
nnz(s)/n^2   % check
toc

tic
F = find(s~=0);
S = full(s(F));
ind = S<1/2;      % or whatever the required condition is; 
                  % these elements are eliminated
F(ind) = [];
S(ind) = [];
[r1 c1] = ind2sub([n n],F);
snew = sparse(r1,c1,S,n,n);
toc

[Benjamin Lender]

Hi David,

thank you! I have tried your code. At the point where I can make the calculation task manager shows 14 GB available memory. I run out of memory at:

F = find(s~=0)

@Bjorn Gustavsson with this in mind my problem is getting to the subset above zero.