Calculating the posterior distribution in multivariate Gaussian processes

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.

-

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

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

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