tall array introduces significant overhead if I call gather() at every iteration in a loop

I want to figure out the spread of the invariants of a symmetric 3x3 matrix using the code below

%input
nPts = 2;
ub = 2.0;
%create combinations of the SIX independent components
comb = combinations(linspace(0.1, ub, nPts), ...
                    linspace(0.1, ub, nPts), ...
                    linspace(0.1, ub, nPts), ...
                    linspace(0.1, ub, nPts), ...
                    linspace(0.1, ub, nPts), ...
                    linspace(0.1, ub, nPts)).Variables;
%results array
res = zeros(nPts^6, 2);
tic
for idx=1:nPts^6
    C = [comb(idx,1), comb(idx,4), comb(idx,5); 
        comb(idx,4), comb(idx,2), comb(idx,6); 
        comb(idx,5), comb(idx,6), comb(idx,3)];
    res(idx, :) =  [trace(C), trace(inv(C))];
end
toc
%creater scatter plot
scatter(res(:,1), res(:,2));

This works, however, I want to set nPts=60 for instance. So I will definitely end up with array sizes that do not fit into memory of my local machine anymore.

As a workaround, I thought storing comb as tall array

comb = tall( combinations(linspace(0.1, ub, nPts), ...
                        linspace(0.1, ub, nPts), ...
                        linspace(0.1, ub, nPts), ...
                        linspace(0.1, ub, nPts), ...
                        linspace(0.1, ub, nPts), ...
                        linspace(0.1, ub, nPts)).Variables );

and then

C = gather(C)
res(idx, :) =  [trace(C), trace(inv(C))];

in the loop.

However, this code takes 170 seconds on my screen for just 2^6 = 64 combinations.

So is tall array not appropriate in my application or do I just use it not correctly?

Answer 1

[Harald]

Hi,

tall arrays are more typically used when you have a single file or a set of files that is too large to be imported into memory at once.

One way I see to solve the challenge you are facing is to create only parts (e.g., leaving one of the variables constant for the time would "only" require 60^6 * 8 * 2 bytes = about 12 GB rather than 60^6 * 8 * 2 bytes = about 746 GB - sizes corrected compared to original post) of the data at a time and do the preprocessing you will need for the scatter plot. For example, you could use histcounts2 to bin the data, accumulate this over the parts, and then use heatmap for visualization.

That's just the idea - please let me know if you want to pursue this and need additional help with the implementation.

Another question will be computational efficiency. For this example, you can explicitly calculate the trace of the inverse, and it is a not overly complicated formula:

syms a b c d e f
M = [a d e; d b f; e f c];
simplify(trace(inv(M)))

That way, you can calculate the trace of a huge number of matrices and their inverses in a vectorized way.

This may also give you additional insights, such as: trace of inverse will be large, when denominator of that expression is close to 0.

Best wishes,

Harald

[SA-W]

Thank you for your answer!

One way I see to solve the challenge you are facing is to create only parts (e.g., leaving one of the variables constant for the time would "only" require < 1 GB rather than about 46 GB)

You mean something like...

a = linspace(0.1,2.0,50) 
for idx=1:numel(a)
    %allocate res with all combinations of the vectors b,c,d,e,f
    res_tmp = ...;
end

?

For example, you could use histcounts2 to bin the data, accumulate this over the parts, and then use heatmap for visualization.

This procedure is not clear to me yet.

So say I have filled a res_tmp array like shown above. How would I implement your strategy to finally come up with a heatmap?

[Harald]

Hi,

I must have edited my answer while you were writing. Please see that for some extra thoughts.

heatmap is not a good way of visualizing the result, scatter seems better. I have given it a shot. I must have made a mistake somewhere because results look different - perhaps a pair of fresh eyes can discover it quickly.

Best wishes,

Harald

nPts = 4;
a = linspace(0.1,ub,nPts);
N = 0;
xstep = 0.1;
ystep = 0.1;
for idx=1:numel(a)
    % all combinations of the vectors b,c,d,e,f
    comb = combinations(linspace(0.1, ub, nPts), ...
                    linspace(0.1, ub, nPts), ...
                    linspace(0.1, ub, nPts), ...
                    linspace(0.1, ub, nPts), ...
                    linspace(0.1, ub, nPts)).Variables;
   
    res_tmp = zeros(nPts^5, 2);
    for idx2=1:nPts^5
        C = [a(idx), comb(idx2,3), comb(idx2,4); 
        comb(idx2,3), comb(idx2,1), comb(idx2,5); 
        comb(idx2,4), comb(idx2,5), comb(idx2,2)];
        res_tmp(idx2, :) =  [trace(C), trace(inv(C))];
    end
    N = N + histcounts2(res_tmp(:,1), res_tmp(:,2), 0:xstep:6, -10:ystep:10);
end
x = repmat((0+xstep/2):xstep:(6-xstep/2), size(N,2), 1);
y = repmat(((-10+ystep/2):ystep:(10-ystep/2))', 1, size(N,1));
N(N == 0) = NaN;
scatter(x(:), y(:), [], N(:))
colorbar

[SA-W]

@Harald

Thanks for your idea.

I must have made a mistake somewhere because results look different

What did you compare your code with? The code from my question plots all points in a scatter, the solution in your last comment plots only a subset of them.

However, an eye-catcher for me is

x = repmat((0+xstep/2):xstep:(6-xstep/2), size(N,2), 1);
y = repmat(((-10+ystep/2):ystep:(10-ystep/2))', 1, size(N,1));

Why do you start and end at (0+xstep/2) and (6-xtsep/2), respectivly, and not 0 and 6, respectively?

[Harald]

I had compared this to your code, using nPts = 4; there as well. Also, I had constrained the y-axes to [-10, 10], since I need to specify a binning and I thought that something like 10^37 is probably just "close to inf" for practical purposes.

For histcounts2, you specify the edges of the bins. For the scatter plot, I wanted to use the centers of the bins. Whether you adjust it this way or another way can be argued, but you need to make some adjustment to overcome the dimension mismatch.

I hope this helps.

Best wishes,

Harald

[SA-W]

@Harald I think I found the mistake in your code:

I had to transpose the N matrix before passing it to scatter, i.e.,

N = N';
scatter(x(:), y(:), [], N(:))
colorbar

Does that make sense?

[Harald]

Good catch, that's it!

I hope that your question is then answered. If so, please consider marking the answer as "accepted".

Thanks and best wishes,

Harald