Original code:
for m = 1 : M
cholA = chol(exchsig(N1(m)-1, N2(m)-1, th));
epsim12 = zeros(N1(m), N2(m), R);
epsimR12 = reshape(mtimesx(repmat(cholA', [1 1 R]), reshape(epsimR(psa(m)+1: psa(m)+(N1(m)-1)*(N2(m)-1), :, :), [(N1(m)-1)*(N2(m)-1) 1 R])), [N1(m)-1 N2(m)-1 R]);
epsim12(2:end, 2:end, :) = epsimR12;
epsim_cell{m} = reshape(epsim12, [N1(m)*N2(m) R]); % Added by Wei
f(psa(m)+1 : psa(m+1), :) = f(psa(m)+1 : psa(m+1), :) - reshape(epsim12, [N1(m)*N2(m) R])
end