Sparse matrix storage is great, but the process for building a sparse matrices can consume a lot of memory. Here's how I dealt with the problem.

Update: I've written a small mex function for Matlab, that builds a sparse matrix directly from data stored on disk. See the article on direct sparse matrix construction for more information.

I'm currently running simulations of large, highly recurrent neural networks in matlab. Large means attempting to simulate all neurons within a 3×3 mm region of the superficial layers of cortex — in the order of 180,000 neurons. Since the simulation is highly recurrent, potentially any neuron can connect to any other. This implies a weight matrix with around 3×10¹⁰ entries, requiring around 240 GB of storage. It's clearly difficult to manipulate that amount of data.

Luckily each neuron makes only around 5,000 synapses, implying that we only need 9×10⁸ non-zero entries in a very sparse matrix (reducing our memory consuption to only 6 GB). Still big, but much more manageable.

Triple the storage, triple the problems!

The problem comes when you have to generate that sparse matrix. Matlab has good sparse support (but only for double-precision or logical values, which is a shame). The way to construct a sparse matrix in matlab is by using the function sparse:

   S = sparse(i,j,s,m,n,nzmax);

To use this function, you have to generate a set of row vectors i, j and s that contain the row and column indices and values that should be assigned to each non-zero entry in the sparse matrix. I was using code to loop over each source neuron, and composing the vectors into cell arrays, something like this:

   cRowInd = {};
   cColInd = {};
   cValues = {};

   for (nNeuron = 1:nNumNeurons)
      [cRowInd{end+1}, ...
       cColInd{end+1}, ...
       cValues{end+1}] = ...
         generate_synapses();
   end

   mfWeights = sparse(vertcat(cRowInd{:}), ...
                  vertcat(cColInd{:}), ...
                  vertcat(cValues{:}), ...
                  nNumNeurons, nNumNeurons);

The loop had to be there since I couldn't generate the synapse probability distributions simultaneously for all neurons, again due to memory limitations. I could have pre-allocated the cell arrays, but that wouldn't make much difference to the execution speed. The main problem is in the final call to sparse: the amount of memory required to build the sparse matrix was more than triple that required by the sparse matrix itself. That's because each of the vectors i, j and s were stored once in the cell arrays, once again when using vertcat to build the vectors for the call to sparse, and finally within the sparse matrix itself.

Being efficient

Luckily I knew the maximum number of non-zero elements I should expect in the sparse matrix (nNumNeurons×nNumSynapsesPerNeurons). I could use this to pre-allocate the i, j and s vectors directly:

   nNonZeroEntries = nNumNeurons * nNumSynapsesPerNeuron;
   vnRowInd = zeros(nNonZeroEntries, 1);
   vnColInd = zeros(nNonZeroEntries, 1);
   vfValues = zeros(nNonZeroEntries, 1);
   nLastInd = 0;

   for (nNeuron = 1:nNumNeurons)
      [vnTheseRowInd, vnTheseColInd, vfTheseValues] = ...
         generate_synapses();

      nNumTheseSynapses = numel(vfTheseValues);
      vnThisWindow = (nLastInd+1):(nLastInd+nNumTheseSynapses);
      vnRowInd(vnThisWindow) = vnTheseRowInd;
      vnColInd(vnThisWindow) = vnTheseColInd;
      vfValues(vnThisWindow) = vfTheseValues;

      nLastInd = nLastInd + nNumTheseSynapses;
   end

   mfWeights = sparse(vnRowInd(1:nLastInd), ...
                     vnColInd(1:nLastInd), ...
                     vfValues(1:nLastInd), ...
                     nNumNeurons, nNumNeurons);

This is slightly wasteful, since I may generate fewer synapses than I expected and so the vectors will be over-allocated. But now we're only using twice as much memory as we need to! The final finesse I made was to make the vectors MappedTensors and so push the data to disk. This involves a performance hit, but it's better than thrashing virtual memory. All this would be much improved if Matlab supported uint classes for providing the indices, which require much less storage than double.

One final option would be to de- and re-compose the sparse matrix on every loop iteration, using something like

   [i,j,s] = find(S);
   S = sparse([i; i2], [j j2]; [s; s2]);

but that is very time consuming, and re-allocates a lot of memory on each iteration.

If you have another suggestion of how to generate sparse matrices without wasting so much space, I'd love to hear it!