Vectorising orthogonal-triangular decomposition for 3D matrix

Question

I am trying to optimise the computation time of my code (for the second time after a first optimisation that gives very good results). Currently, this code is very time-consuming depending on the size of the data (i.e. sometimes very big matrices). So if anyone has some ideas to help me to optimize this code, I will be very grateful!

matrice1=0+(20-0).*rand(500,5);
matrice2=datasample(matrice1,100);

[n1,~]=size(matrice1);
[n2,p2]=size(matrice2); 
% Pre-calculate replicated B and the indices to be modified at each iteration
B_rep=repmat(matrice2,[1 1 n2]);
B_idx=bsxfun(@plus,((0:p2-1)*n2+1)',(0:n2-1)*(n2*p2+1));
B=mean(B_rep,1);
out=zeros(n1,n2); % initialize output array
for i=1:n1
    B_rep(B_idx)=repmat(matrice1(i,:)',1,n2);   
    B_meanB=bsxfun(@minus,B_rep,B); % B minus mean values of B
    A_B_meanB=matrice2'-reshape(B,p2,[]); % A minus B_meanB    

    for j=1:n2
        [~,R]=qr(B_meanB(:,:,j),0);   
        out(i,j)=sum((R'\A_B_meanB(:,j)).^2)*(n2-1);    
    end
end

Answer 1

Your code looks fairly good, but there are a few things I would do differently:

• It's very confusing that you call B_rep for B, and call B the mean of B_rep. The comment and code here looks very strange. You should call B_mean = mean(B_rep, 1) to stick with the general naming convention in your code.

  B_meanB=bsxfun(@minus,B_rep,B); % B minus mean values of B
  

• bsxfun performs better than repmat, so instead of B_rep=repmat(matrice2,[1 1 n2]); you can do:

  B_rep = bsxfun(@times, matrice2, ones(1,1,n2)); 
  

• Instead of ', you should use .' when transposing an array. The first one is the complex conjugated transpose.

• i and j are bad variable names in Matlab.

• You should try to use more spaces. It makes the code much easier to read.

---

I don't have the statistical toolbox, so I can't test your code myself. I'm not sure how it can be vectorized, so I can't help with much when it comes performance gain I'm afraid.