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.

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.

StackExchange Code Review Q#21234, answer score: 2