How to get rid of this nested for cycles?

Here is my code:

% simulating similar inputs:
K = 3;
H = 250;
N = 50;
B = 10;
Q = cell(1,K);
rng('default')
Q = cellfun(@(x) rand(H, N, B) - rand(H, N, B),Q, 'UniformOutput', false);

% here starts my problem:
cumret  = cellfun(@cumsum,Q,'UniformOutput',false);

rev = cell(1,3);

% i would like to get rid of the following triple for clycles:
for k = 1:K
 for i = 1:N
   for b = 1:B
       rev{1,k}(:,b,i) = cumret{1,k}(:,i,b);
   end
 end
end

% and to get this very same result (whic is correct but slow to compute):
rev = cell2mat(rev);

% in other terms the following two plot must be the same:
figure(1)
for k = 1:K
    for b = 1:B
        plot(cumret{1,k}(:,1,b))
        hold on
    end
end

figure(2)
plot(rev(:,:,1))

Thank you in advance for your kind help.

Answer 1

[Bruno Luong]

There is no reason using cell, rather 4-dimensional array.

Admitely the code is not clear if you don't keep track which dimension goes where, but much shorter. Each variable is generate by a single line.

K = 3;
H = 250;
N = 50;
B = 10;

rng('default')

Q = rand(H, N, B, K) - rand(H, N, B, K);
% here your problem is solved
cumret  = cumsum(Q); % H x N x B x K
% rev is reorganization of cumret in this way
rev = permute(cumret(:,:,:),[1 3 2]);  % H x B*K x N

figure(1)
hold on
for k = 1:K
    for b = 1:B
        plot(cumret(:,1,b,k))
    end
end

figure(2)
plot(rev(:,:,1))

[Barbab]

Thank you for providing the alternative code.

I opted for cells because each of the matrices stored in the cell array Q represents distinct and different data. Using cells allowed me to keep them separate and retain their individuality, which was essential for my reasoning and analysis.

Additionally, I found the cell-based implementation to be more intuitive for my specific use case, as it closely mirrors the problem's structure and provides clarity in keeping track of dimensions.

Do you think that 4D matrices are more efficient from a computational POV?

[Bruno Luong]

" stored in the cell array Q represents distinct and different data."

This is purely psychologic. You purposingly create a barrier for your own good. There is no different to do

Qc{l}(i,j,k)

than

Qa(i,j,k,l)

"Do you think that 4D matrices are more efficient from a computational POV?"

Surely. If your data can be put on nd-array, then by allmean do it.

Just look at your cell code: it looks much more complicated and not efficient as mine code with 4D array.

Wth cell you will make efforts to convert back and forth to array with cell2mat, mat2cell; celfun, for-loop ... and you don't take advantage of MATLAB vectorized library for array manipulation.

[Barbab]

Ok I understand. Let me show you where the matrix Q comes from:

% this code is for the wild bootstrap technique (block sampling)

% simulating similar inputs:
K = 3;
H = 75;     % in the actual code H = 250
T = 500;    % in the actual code T = 2274
N = 50;     % in the actual code N = 58
B = 10;     % in the actual code B = 10000

rng('default')

% R is observed
R = rand(T, N, K) - rand(T, N, K);

% define the variables to be bootstrapped
Q = cell(1,K);

% variables pre allocation
for k = 1:K
    Q{k} = nan(H,N,B);
end

% this nested loop cycles takes ages with the "actual" H,T,N and B
parfor k = 1:K
    for n = 1:N
        for b = 1:B
            % randomly pick a number between 1 and T - H + 1
            idx = unidrnd(T - H + 1);
            % wild bootstrap: extract R's rows from idx to idx + H - 1
            Q{k}(:,n,b) = R(idx:(idx + H - 1),n,k);
        end
    end
end

% then my question code...
cumret  = cellfun(@cumsum,Q,'UniformOutput',false);

rev = cell(1,3);

% i would like to get rid of the following triple for clycles:
for k = 1:K
 for i = 1:N
   for b = 1:B
       rev{1,k}(:,b,i) = cumret{1,k}(:,i,b);
   end
 end
end

% and to get this very same result (whic is correct but slow to compute):
rev = cell2mat(rev);

% in other terms the following two plot must be the same:
figure(1)
for k = 1:K
    for b = 1:B
        plot(cumret{1,k}(:,1,b))
        hold on
    end
end
figure(2)
plot(rev(:,:,1))

So, my main questions are:

1. How can I efficiently obtain the matrix Q of size H x B x N x K directly during the wild bootstrap process? As you can see in the code, the current nested loop implementation takes a considerable amount of time due to the large values of H, T, N, and B.
2. Additionally, I am looking for suggestions on how to optimize and speed up my code as much as possible. Currently, the computation takes hours, if not days, to complete on my laptop, which is not brand new anymore.
3. Finally, how do you store such large matrices? Do you suggest to save it as .mat files? Currently I have more than one file.mat larger than 13.5gb and they make MATALB crash every time I load them...

I would greatly appreciate any insights or alternative approaches that could help me improve the performance of this wild bootstrap technique. Thank you in advance for your assistance!

[Bruno Luong]

You can do the same with array

for k = 1:K
    for n = 1:N
        for b = 1:B
            % randomly pick a number between 1 and T - H + 1
            idx = unidrnd(T - H + 1);
            % wild bootstrap: extract R's rows from idx to idx + H - 1
            Q(:,n,b,k) = R(idx:(idx + H - 1),n,k);
        end
    end
end

And yes I can accelerate such loop without using parfor.

But if you are more comfortable wih cell then stay with it.

[Barbab]

May I know how can I accelerate the loop?

And what about the storage problem? I have six matrices similar to Q (250 x 10000 x 58 x 3) I am currently using the "matfile" function without actually loading the variables but is very very slow.

Thank you so much for your help.

[Bruno Luong]

@Barbab "May I know how can I accelerate the loop?"

Here the vectorization version. On TMW Online serve it's 14 times faster than for-loop and 22 ties faster than cell with parfor (your code) tested with the toy dimension you have provided.

It is also important regarding to computational POV that you pay attention in designing the right order of the dimensions to your data so that to avoid permute and accessing to data in the contiguous way.

K = 3;
H = 75;     % in the actual code H = 250
T = 500;    % in the actual code T = 2274
N = 50;     % in the actual code N = 58
B = 10;     % in the actual code B = 10000

R = rand(T, N, K) - rand(T, N, K);

%t1=timeit(@() BootstrapCellParLoop(R, H, B), 1)
t2=timeit(@() BootstrapLoop(R, H, B), 1)

t2 = 0.0026

t3=timeit(@() BootstrapVec(R, H, B), 1)

t3 = 1.8133e-04

try
    r13 = t1/t3;
    fprintf('BootstrapVec is %g faster than BootstrapCellParLoop\n', r13)
end

r23 = t2/t3;
fprintf('BootstrapVec is %g faster than BootstrapLoop\n', r23)

BootstrapVec is 14.4863 faster than BootstrapLoop

function Q = BootstrapCellParLoop(R, H, B)
[T, N, K] = size(R);
Q = cell(1,K);
% variables pre allocation
for k = 1:K
    Q{k} = nan(H,N,B);
end
parfor k = 1:K
    for n = 1:N
        for b = 1:B
            % randomly pick a number between 1 and T - H + 1
            idx = randi(T - H + 1);
            % wild bootstrap: extract R's rows from idx to idx + H - 1
            Q{k}(:,n,b) = R(idx:(idx + H - 1),n,k);
        end
    end
end
end

function Q = BootstrapLoop(R, H, B)
[T, N, K] = size(R);
Q = zeros(H,N,B,K);
for k = 1:K
    for n = 1:N
        for b = 1:B
            % randomly pick a number between 1 and T - H + 1
            idx = randi([1 T-H+1]);
            % wild bootstrap: extract R's rows T,Nrom idx to idx + H - 1
            Q(:,n,b,k) = R(idx:(idx + H - 1),n,k);
        end
    end
end
end


function Q = BootstrapVec(R, H, B)
[T, N, K] = size(R);
IDX = randi([1 T-H+1], [1 B]);
IDX = IDX + (0:H-1).';
Q = R(IDX,:);
Q = reshape(Q,[H,B,N,K]);
Q= permute(Q,[1 3 2 4]);
end

[Bruno Luong]

And what about the storage problem

Inside the RAM the total amount of cell or ND array is almost the same (ND has less overhead) but it requires contiguous block of RAM to store it.

In "matfile" it is probably doesn't matter but I have little experience with it to advise you.

Answer 2

[Voss]

I assume that you cannot change Q, and you need to calculate rev from Q.

% simulating similar inputs:
K = 3;
H = 250;
N = 50;
B = 10;
Q = cell(1,K);
rng('default')
Q = cellfun(@(x) rand(H, N, B) - rand(H, N, B),Q, 'UniformOutput', false);

Here is your original code for calculating rev from Q:

% here starts my problem:
cumret  = cellfun(@cumsum,Q,'UniformOutput',false);

rev = cell(1,3);

% i would like to get rid of the following double for clycles:
for k = 1:K
    for i = 1:N
        for b = 1:B
            rev{1,k}(:,b,i) = cumret{1,k}(:,i,b);
        end
    end
end

% and to get this very same result (whic is correct but slow to compute):
rev = cell2mat(rev);

Here is another way to calculate rev from Q. See the comments in the code for an explanation of each step and the links for a detailed description of each function used:

rev_test = permute( reshape( cumsum( cat(4, Q{:}), 1), [H N B*K]), [1 3 2]);
%                                    ^^^^^^^^^^^^                             concatenate the contents of the cells of Q, making a 250x50x10x3 array
%                            ^^^^^^^^^^^^^^^^^^^^^^^^                         then cumsum that 4D array along its 1st dimension
%                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^             reshape that 4D cumsum'ed array to be 250x50x30
%          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^  permute to 250x30x50

% and it gives the same result:
isequal(rev_test,rev)

ans = logical
   1

https://www.mathworks.com/help/matlab/ref/double.cat.html

https://www.mathworks.com/help/matlab/matlab_prog/comma-separated-lists.html

https://www.mathworks.com/help/matlab/ref/cumsum.html#d124e293537

https://www.mathworks.com/help/matlab/ref/reshape.html

https://www.mathworks.com/help/matlab/ref/permute.html

[Barbab]

Thank for your response.

Yes, exactly, I cannot easily change Q. I can see that you efficiently used a combination of array concatenation, cumulative sum, reshape, and permutation functions to achieve the desired result. While I use MATLAB regularly, I sometimes find it challenging to fully grasp the functionalities of reshape and permute functions as well.

One thing: when you cite my code, the triple for cycle is different from mine. Is there a specific reason for that?

[Bruno Luong]

More to the point, you see in the Voss's code the following,

... cat(4, Q{:}) ...

It simply first converts cell array to 4D array.

If you insist on using cell fomat you would have to do it many many times.

[Voss]

"when you cite my code, the triple for cycle is different from mine. Is there a specific reason for that?"

No reason, except that I made a mistake. I began by editing your code to remove some loops, and I mistakenly didn't put it back the way you had it. I've edited my answer so that your original code should be accurate now.

Answer 3

[Bruno Luong]

If you want to work with cell

% simulating similar inputs:
K = 3;
H = 250;
N = 50;
B = 10;
Q = cell(1,K);
rng('default')
Q = cellfun(@(x) rand(H, N, B) - rand(H, N, B),Q, 'UniformOutput', false);

% here starts my problem:
cumret  = cellfun(@cumsum,Q,'UniformOutput',false);

% rev = cell(1,3);
% for k = 1:K
%  for i = 1:N
%    for b = 1:B
%        rev{1,k}(:,b,i) = cumret{1,k}(:,i,b);
%    end
%  end
% end
rev = cellfun(@(a) permute(a, [1 3 2]), cumret, 'UniformOutput',false);

% and to get this very same result (whic is correct but slow to compute):
% rev2 = cell2mat(rev);
rev = cat(2, rev{:}); % faster

% in other terms the following two plot must be the same:
figure(1)
for k = 1:K
    for b = 1:B
        plot(cumret{1,k}(:,1,b))
        hold on
    end
end
figure(2)
plot(rev(:,:,1))