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.

• Adriaan: How about an N-by-4 array with columns Nx, Nz, Nt, Sol? In case that also becomes too large for your RAM, you might have to write it to disk, in which case you could e.g. store a sparse sol(Nx, Nz) on each time step.

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

• nekomatic: I probably wouldn't override the builtin get though.
• Wolfie: would sub2ind help avoid having to define custom indexing functions? Docs: Convert subscripts to linear indices. Looks to handle N-dimensional indices
• Luis Mendo: @nekomatic Whoops, of course! Thanks
• Luis Mendo: @Wolfie I'm so used to doing it manually that I hadn't thought of sub2ind. Anyway, that requires indices of the same size. It doesn't do the "Cartesian product" of indices (akin to implicit expansion) like normal indexing does

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.

• Cris Luengo: This is what I wanted to suggest as well. It's likely not necessary to keep everything in the same matrix, and computations are likely more efficient if each matrix represents one time step.