Efficiently store an N-dimensional array of mostly zeros in Matlab

Question

I implemented a finite differences algorithm to solve a PDE.

The grid is a structured 2D domain of size [Nx, Nz], solved Nt times.

I pre-allocate the object containing all solutions:

sol = zeros(Nx, Nz, Nt, 'single') ;

This becomes very easily too large and I get a 'out of memory' error. Unfortunately sparse doesn't work for N-dimensional arrays.

For the sake of the question it's not important to know the values, it goes without saying that the RAM usage grows exponentially with decreasing the grid spacing and increasing the simulation time.

I am aware that I do not need to store each time instant for the purpose of the advancement of the solution. It would be sufficient to just store the previous two time steps. However, for post-processing reasons I need to access the solution at all time-steps (or at least at a submultiple of the total number).It might help to specify that, even after the solution, the grid remains predominantly populated by zeros.

Am I fighting a lost battle or is there a more efficient way to proceed (other type of objects, vectorization...)?

Thank you.

Answer 1

You could store the array in sparse, linear form; that is, a column vector with length equal to the product of dimensions:

sol = sparse([], [], [], Nx*Nz*Nt, 1); % sparse column vector containing zeros

Then, instead of indexing normally,

sol(x, z, t),

you need to translate the indices x, z, t into the corresponding linear index:

• For scalar indices you use

  sol(x + Nx*(z-1) + Nx*Nz*(t-1))
  

  You can define a helper function for convenience:

  ind = @(sol, x, y, t) sol(x + Nx*(z-1) + Nx*Nz*(t-1))
  

  so the indexing becomes more readable:

  ind(sol, x, z, t)
  

• For general (array) indices you need to reshape the indices along different dimensions so that implicit expansion produces the appropriate linear index:

  sol(reshape(x,[],1,1) + Nx*(reshape(z,1,[],1)-1) + Nx*Nz*(reshape(t,1,1,[])-1))
  

  which of course could also be encapsulated into a function.

Check that the conversion to linear indexing works (general case, using non-sparse array to compare with normal indexing):

Nx = 15; Nz = 18; Nt = 11;
sol = randi(9, Nx, Nz, Nt);
x = [5 6; 7 8]; z = 7; t = [4 9 1];
isequal(sol(x, z, t), ...
    sol(reshape(x,[],1,1) + Nx*(reshape(z,1,[],1)-1) + Nx*Nz*(reshape(t,1,1,[])-1)))

gives

ans =
  logical
   1

Answer 2

You can create a a cell array of sparse matrices to store the results. However computations can be performed on full matrices if working with a full matrix is faster than sparse matrix and convert the full matrix to sparse matrix and place it in the cell.