memory problems with ndgrid

Question

In Matlab, I am using "ndgrid" to create a 6D matrix. Here is the code:

for i=1:3
    dd{i}=[0 0 0 0 0 0 0 1 1 1 1 1];
    ss{i}=[0 0 0 0 0 0 1 1 1 1 1 1];
end 
[D1 D2 D3 S1 S2 S3] = ndgrid(dd{1},dd{2},dd{3},ss{1},ss{2},ss{3});
out = D1.*S1.*D2.*S2.*D3.*S3;   

The problem that I am having is that although I have plenty of memory to store one or two 6-D matrices, I do not have enough memory to store all 6 matrices:

[D1 D2 D3 S1 S2 S3] 

As you can see, the matrix D1, D2... and "out" are sparse in nature, but the "sparse" functionality in Matlab does not work on multidimensional arrays. I have searched for other "grid" functionality in Matlab, but I can not find a "grid" function that helps me avoid the intermediary step of calculating D1, D2 etc.

In general, I want to allow dd{1} to be different than dd{2}. Also, this 6-D case I posted doesn't take up too much memory, but the 8-D case does, which is where I am running into an issue.

Any help making more efficient use of memory in this situation is much appreciated.

Answer 1

Here's a case where a loop will almost certainly be faster (and will fit!), since you're memory-bound. Just use 6 (or 8) nested for-loops, and drop in each element of out, one at a time. Use the nested loop indices to reference directly into your source vectors, and just be sure you match the correct index with the correct vector. This will limit you only to the maximum size of your output array.

Editing with a new idea:

Ooooh, I just thought of bsxfun, which would do exactly the right thing, except it only accepts two arguments. It essentially replicates as necessary across any singleton dimensions. Pairwise-nesting bsxfun calls would be pointless, as the intermediate results would be back to full size. But it is a great way to reduce your dimension by half, which for a power term is a big deal:

[D1 D2 D3] = ndgrid(dd{1},dd{2},dd{3});
[S1 S2 S3] = ndgrid(ss{1},ss{2},ss{3});
out = bsxfun(@times, D1.*D2.*D3, reshape(S1.*S2.*S3, [1 1 1 12 12 12]));

Obviously the 12s have to be replaced by dynamic sizes. Now, instead of D*N^D elements, you only need D*N^(D/2). You won't have a memory problem until you're back up to 14 or 16 dimensions.