How to get rid of this nested for cycles?
2 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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.
0 Commenti
Risposta accettata
Bruno Luong
il 31 Lug 2023
Modificato: Bruno Luong
il 31 Lug 2023
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))
7 Commenti
Bruno Luong
il 1 Ago 2023
Modificato: Bruno Luong
il 2 Ago 2023
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)
t3=timeit(@() BootstrapVec(R, H, B), 1)
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)
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
il 1 Ago 2023
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.
Più risposte (2)
Voss
il 31 Lug 2023
Modificato: Voss
il 1 Ago 2023
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)
3 Commenti
Bruno Luong
il 1 Ago 2023
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
il 1 Ago 2023
"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.
Bruno Luong
il 1 Ago 2023
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))
0 Commenti
Vedere anche
Categorie
Scopri di più su Matrix Indexing in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!