Vectorising orthogonal-triangular decomposition for 3D matrix

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


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.

StackExchange Code Review Q#121075, answer score: 3