MATLAB parfor is slower than for -- any help to speed up of my loop cacluations within KDL function?

II=9;
JJ=9;
M=11;
W=rand(II+1,JJ+1,3,M);
k_d=01e12;
syms x y
for i=1:II+1
    for j=1:JJ+1
        lgij(i,j) = legendreP(i-1, x)*legendreP(j-1, y);
    end
end 
wxy2 = sym('wxy2',[1 M]);
wxy3 = sym('wxy3',[1 M]);
wxy2(1:M) = sym('0');
wxy3(1:M) = sym('0');  
for r=1:M
    for i=1:II+1
        for j=1:JJ+1
            wxy2(r) = W(i, j, 2, r)*lgij(i,j) + wxy2(r); 
            wxy3(r) = W(i, j, 3, r)*lgij(i,j) + wxy3(r);
        end
    end
end 
wxxyy2=simplify(wxy2);
wxxyy3=simplify(wxy3);
for r=1:M
    for o=1:M
        Wijklmo(r,o)= expand(wxxyy2(r)*wxxyy2(o) + wxxyy3(r)*wxxyy3(o)- wxxyy2(r)*wxxyy3(o) - wxxyy2(o)*wxxyy3(r));
    end
end
dt = 3.8655e-08;
T=1e9*dt;
q = zeros(M, floor(taim/dt)+1);
RA=111;
s=1;
for tn=t__0:dt:T
    p_hat(1:M,s) = KDL(wxxyy2, wxxyy3, M, q, RA, k_d, Wijklmo, s)
    s=s+1
end

a

function kdl = KDL( wxxyy2, wxxyy3, M, q, RA, k_d, Wijklmo, s)
xx = sym('xx');
yy = sym('yy');
Wxy2 = sym('0');
Wxy3 = sym('0');
for r=1:M
    Wxy2 = wxxyy2(r)*q(r, s) + Wxy2;
    Wxy3 = wxxyy3(r)*q(r, s) + Wxy3;
end 
w23 = expand(Wxy3-Wxy2);
kdl=zeros(M,1);
xspan = [-1 1];
yspan = [-1 1];
Gy0 = 0;
H(xx,yy) =(0.5*(1+tanh(k_d*RA*w23)));
g = matlabFunction(H,'Vars',[xx yy]);
if w23==0
   kdl(:,1) =0;
else
   tic
   for i=1:M
   %parfor i=1:M    
       ff(i,1)= vpa(Wijklmo(i,:)*q(1:M, s));
       f = matlabFunction(ff(i,1),'Vars',[xx yy]);
       D = @(xx,yy)f(xx,yy).*(g(xx ,yy)>0);
       if w23==0
          kdl(i,1) =0;
       else
       [~,G] = ode45(@(y,Gy)fun(y,Gy,D,xspan),yspan,Gy0,odeset('RelTol',1e-10,'AbsTol',1e-10));
       kdl(i,1) = k_d*real(G(end));
   end
end
toc
end

a

function dGydy = fun(y,Gy,g,xspan)
  % Compute the x-integrals at y = y
  Gx0 = 0;
  [~,Gx] = ode45(@(x,~)g(x,y),xspan,Gx0,odeset('RelTol',1e-10,'AbsTol',1e-10));
  dGydy = Gx(end);
end

[Torsten]

It can be decided outside the for loop whether w23 == 0 and thus kdl(1:M) = 0. So you can take the if-statement out of the for-loop.

Same for the definition of g.

[Torsten]

I meant

function kdl = KDL( wxxyy2, wxxyy3, M, q, RA, k_d, Wijklmo, s)
  xx = sym('xx');
  yy = sym('yy');
  Wxy2 = sym('0');
  Wxy3 = sym('0');
  for r=1:M
    Wxy2 = wxxyy2(r)*q(r, s) + Wxy2;
    Wxy3 = wxxyy3(r)*q(r, s) + Wxy3;
  end 
  w23 = expand(Wxy3-Wxy2);
  kdl = zeros(M,1);
  if w23==0
    return
  end
  H(xx,yy) =(0.5*(1+tanh(k_d*RA*w23)));
  g = matlabFunction(H,'Vars',[xx yy]);
  parfor i=1:M    
       ff(i,1)= vpa(Wijklmo(i,:)*q(1:M, s));
       f = matlabFunction(ff(i,1),'Vars',[xx yy]);
       D = @(xx,yy)f(xx,yy).*(g(xx,yy)>0);
       [~,G] = ode45(@(y,Gy)fun(y,Gy,D,xspan),yspan,Gy0,odeset('RelTol',1e-10,'AbsTol',1e-10));
       kdl(i,1) = k_d*real(G(end));
  end
end

[Mehdi]

I did this, but did not speed-up. The problem is that parfor is slower than for.

[Torsten]

Ok, then I think @Matt J is correct: the gain in speed is overcompensated by the overhead of parallel computing.

Answer 1

[Matt J]

M=11 is a rather small number of iterations. If this is a realistic value, there may not be enough iterative work to make the overhead of parfor worthwhile. It might be better to move the parfor to the outer loop that calls KDL.

Tn=t__0:dt:T; %where are these values used???
parfor s=1:numel(Tn)
    
    p_hat(:,s) = KDL(wxxyy2, wxxyy3, M, q, RA, k_d, Wijklmo, s);
    
end

A few other remarks,

(1) Things like this loop

for r=1:M
    Wxy2 = wxxyy2(r)*q(r, s) + Wxy2;
    Wxy3 = wxxyy3(r)*q(r, s) + Wxy3;
end 

can be replaced with vectorized statements,

 Wxy2 = wxxyy2*q(:, s);
 Wxy3 = wxxyy3*q(:, s);

(2) Also, q(:, s) will probably evaluate faster than q(1:M,s).

[Mehdi]

faced error:

Error: Unable to classify the variable 'p_hat' in the body of the parfor-loop. For more information, see Parallel for Loops in MATLAB, "Solve Variable Classification Issues in parfor-Loops".

[Matt J]

You should pre-allocate p_hat

M=11;
N=5;
p_hat=zeros(M,N);

parfor s=1:N
    
    p_hat(:,s) = KDL(s);
    
end

Starting parallel pool (parpool) using the 'Processes' profile ...
Parallel pool using the 'Processes' profile is shutting down.

p_hat

p_hat = 11×5

1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5

function s=KDL(s)

end

[Mehdi]

this method does not work for me since in my original problem KDL(s+1)=f(KDL(s)).(because of complexities I have not brought this part of my code here.