# Matlab: Can this triple loop be simplified or eliminated

Matlab profiler throws this horrible bottleneck which is called 700000 times:

for k =1:model.nlf,
for r =1:model.nlf,
KuyDinvKyu = zeros(model.k,model.k);
for q =1:model.nout,
KuyDinvKyu = KuyDinvKyu + model.KuyDinv{k,q}*model.Kyu{q,r};
end
if (k == r)
model.A{k,r} = model.Kuu{k} + KuyDinvKyu;
else
model.A{k,r} = KuyDinvKyu;
end
end
end


Even if the math is correct, there must be a faster way. Thanks for any solution or direction of search.

Thanks to Rody for suggesting the modified loop excluding IF and with better preallocation.

KuyDinvKyu_0 = zeros(model.k);
for k =1:model.nlf,
for r =1:model.nlf,
KuyDinvKyu = KuyDinvKyu_0;
for q =1:model.nout,
KuyDinvKyu = KuyDinvKyu + model.KuyDinv{k,q}*model.Kyu{q,r};
end
model.A{k,r} = KuyDinvKyu;
end
end
for c =1:model.nlf,
model.A{c,c} = model.A{c,c} + model.Kuu{c};
end


Result: The Matlab profiler upgraded this code from red to pink which shows that the inner loop is really the tough nut.

Further up in the same routine which updates two matrices A and D inside a model array containing all the structural, logical, mathematical information about the model, here are two additional loops red-flagged by the profiler.

Dimensions: nlf 4; nout 4; KuuinvKuy 4x4 cell, each cell 50x650 double; Kuuinv 1x4 cell, each cell 50x50 double; Kyu 4x4 cell, each cell 650x50 double; Kyy 4x1 cell, each cell 650x1 double.

for r =1: model.nlf,
for k =1: model.nout,
model.KuuinvKuy{r,k} = model.Kuuinv{r}*model.Kyu{k,r}'; %'
end
end

for k =1: model.nout,
KyuKuuinvKuy = zeros(size(model.Kyy{k},1),size(model.Kyy{k},1));
for r =1: model.nlf,
KyuKuuinvKuy = KyuKuuinvKuy + model.Kyu{k,r}*model.KuuinvKuy{r,k};
end
switch model.approx
case 'dtc'
model.D{k} = sparseDiag(1/model.beta(k)*ones(size(model.X{k+model.nlf},1),1));
model.Dinv{k} = sparseDiag(model.beta(k)*ones(size(model.X{k+model.nlf},1),1));
model.logDetD{k} = -size(model.X{k+model.nlf},1)*log(model.beta(k));
model.Ktilde{k} = model.Kyy{k} - diag(KyuKuuinvKuy);
otherwise
error('Unknown approximation type')
end

end


All similar inner loops have same built-in inefficiencies. The real question is: do the usual tips and tricks of vectorisation apply on cells array.

Thank you again for your help

-

## migrated from stackoverflow.comFeb 2 '13 at 18:01

This question came from our site for professional and enthusiast programmers.

some dimensions and sizes can be useful here. What are size(model.KuyDinv{k,q}), size(model.Kyu{q,r}) and size(model.Kuu{k})? –  Shai Jan 27 '13 at 17:50
I also think that a bit of intuition into what this code is trying to compute (very high level motivation) can be of assistance. –  Shai Jan 27 '13 at 17:51
Is there a specific reason that all variables are cells? –  natan Jan 27 '13 at 18:45
As above; this questions is impossible to answer without knowing the dimensions of model.KuyDinv{k,q}. Are these vectors? Scalars? –  Pete Jan 27 '13 at 19:19
concerning dimensions model.nlf is number of latent forces: up to 20 lf; nout is number of outputs also up to 20. –  user2015897 Jan 27 '13 at 21:00

You will get a speed up by eliminating the if-statements in the inner loop. E.g. by splitting the inner loop via

    for r = setdiff(1:model.nlf,k)
%do the stuff to setup Model.A{k,r}
end

model.A{k,k} = ...

-
Thank you for the hint. It has been confirmed and implemented with the second answer below. The profiler has improved from red to pink, a step in the right direction. Many thanks for your help –  user2015897 Jan 29 '13 at 19:55

Assuming model.KuyDinv{k,q} and model.Kyu{q,r} are a matrices (as the name would imply), there is little you can do. You can move the initialization of KuyDinvKyu outside the loop, and eliminate the branch to compute the diagonals:

% define it once -- saves many calls to 'zeros'
KuyDinvKyu_0 = zeros(model.k);

for k = 1:model.nlf
for r = 1:model.nlf

KuyDinvKyu = KuyDinvKyu_0;

for q = 1:model.nout
KuyDinvKyu = KuyDinvKyu + model.KuyDinv{k,q}*model.Kyu{q,r};
end

% removed IF
model.A{k,r} = KuyDinvKyu;

end
end

% without IF
for c = 1:model.nlf
model.A{c,c} = model.A{c,c} + model.Kuu{c};
end


If model.KuyDinv{k,q} and/or model.Kyu{q,r} contain scalars, well then we can optimize this much further for sure. So I need to know the size and type of data that model.Kyu{q,r} and model.KuyDinv{k,q} will contain.

It might also be that your overall data design (e.g., the choice to use cells) is flawed and causes inefficiencies. So it could also be helpful to see some more surrounding code, so I can determine if there is some improvement to be made in that respect.

-
Thanks Rody for this very helpful reply. I implemented your solution and improved the profiler color code for that loop from red to pink. model.KuyDinv and all the others are matrices. I inherited this routine and I would agree that the overall data design using cells is the problem. –  user2015897 Jan 29 '13 at 20:04
@user2015897: I think you suggested an edit to my answer, that got rejected -- you should indeed post this as part of the question (rather than my answer), so that not only I but everyone here can take a look :) –  Rody Oldenhuis Jan 29 '13 at 20:36