Optimization of nested "for" loops

Hi.

I would like to optimize my code. I examined it and found a serious bottle neck. Below is the chunk of code with nested "for" loops which in my opinion are very inefficient.

PHI    = zeros(length(X),length(X));
PHI_X  = zeros(length(X),length(X));
PHI_XX = zeros(length(X),length(X));
PHI_Y  = zeros(length(X),length(X));
PHI_YY = zeros(length(X),length(X));

for i=1:length(X)
    for j=1:length(Y)
        x = X(i)-X(j);
        y = Y(i)-Y(j);
        PHI    (i,j) =                               (x^2 + y^2 + e^2) ^  p;
        PHI_X  (i,j) = 2 * p * x                   * (x^2 + y^2 + e^2) ^ (p-1);
        PHI_Y  (i,j) = 2 * p * y                   * (x^2 + y^2 + e^2) ^ (p-1);
        PHI_XX (i,j) = 2 * p * ((2*p-1)*x^2 + e^2) * (x^2 + y^2 + e^2) ^ (p-2);
        PHI_YY (i,j) = 2 * p * ((2*p-1)*y^2 + e^2) * (x^2 + y^2 + e^2) ^ (p-2);
    end
end

X and Y are vectors while e and p are real-valued coefficients. This part of the code is very important because it is executed even several hundred times (for different values of e and p) during one code execution. I would like to optimize it somehow but I don't know how because I'm quite a new user of MATLAB.

Kind regards

Answer 1

[DGM]

Should be able to do something like this without the loop (don't need to preallocate either):

x = X'-X;
y = Y'-Y;
PHI    =                                 (x.^2 + y.^2 + e^2) .^  p;
PHI_X  = 2 * p * x                    .* (x.^2 + y.^2 + e^2) .^ (p-1);
PHI_Y  = 2 * p * y                    .* (x.^2 + y.^2 + e^2) .^ (p-1);
PHI_XX = 2 * p * ((2*p-1)*x.^2 + e^2) .* (x.^2 + y.^2 + e^2) .^ (p-2);
PHI_YY = 2 * p * ((2*p-1)*y.^2 + e^2) .* (x.^2 + y.^2 + e^2) .^ (p-2);

Those first couple lines will only work in R2016b or later, otherwise, you'd need to do:

x = bsxfun(@minus,X,X');

etc

[Olaf Popczyk]

Thank you!

Answer 2

[Matt J]

Replace the loop with,

x=X(:)-X(:).';
y=Y(:)-Y(:).';

xsq=x.^2;
ysq=y.^2;
xysq= xsq+ysq;

tmp1=xysq + e^2;

        PHI    =  tmp1 .^  p;
        
tmp2=tmp1 .^ (p-1);

        PHI_X  = 2 * p * x .*tmp2;
        PHI_Y  = 2 * p * y .*tmp2;
        
        clear tmp2
      
tmp3=tmp1 .^ (p-2);       

        PHI_XX = (2*p*(2*p-1)*xsq + 2*p*e^2) .* tmp3;
        PHI_YY = (2*p*(2*p-1)*ysq + 2*p*e^2) .* tmp3;

[DGM]

That's a big improvement too.