How to avoid the for loops?

I want to make the following code faster by avoiding the two for loops (if possible). I tried to make a vector multiplication .* operator which is not working here.

Since, the arguments / variables (Order, X and Y) are small, the Y_grid cell array can be calculated so fast here in the given script. However, I would like to make it as a function for later use. (So that, I can input different and higher arguments). Then the two for loops would make the programme so slow.

Inside the loops, there is a function too - (matFunction.m)*, is a small function here which needs varying input arguments. I don't want to make any change in this function.

Please have a look and suggest me some ideas! Many thanks in advance******

Order = 5; 
X  = 0:3;
Y  = 0:5;
Y_grid {length(Y), length(X)}=zeros; % Initialization of cell array by zeros.
tic
  for i=1:length(X)
      for j=1:length(Y)
         Y_grid{j,i}= matFunction( Order,Y(j),X(i) );
      end
  end
toc

(The matFunction.m* is given below. ):

function [ Mat ] = matFunction( Order , Y , X )  
Mat=[0:1:Order]*X*Y;
Mat=repmat(Mat,Order,Order);    
end

Answer 1

[Stephen23]

Here is a solution using matrix multiplication and one bsxfun call:

Order = 5; 
X  = 0:3;
Y  = 0:5;
A = ones(Order,1)*(0:Order);
B = bsxfun(@times,A,reshape(Y(:)*X,1,1,[],numel(X)));
C = squeeze(num2cell(repmat(B,1,Order),[1,2]));

Under some circumstances this code is faster than the nested loops and separate function, but this depends on factors such as the matrix sizes, etc. With the values you give it is around 1/3 to 1/2 faster:

Elapsed time is 0.002530 seconds. % nested loops + function
Elapsed time is 0.001254 seconds. % matrix + bsxfun

To really make the code faster I would avoid cell arrays entirely, and just stick to using the dimensions of the numeric array. This makes writing vectorized code much more intuitive and easier.

[Jacky Jo]

Thanks for your answer Stephen Cobeldick.

But the situation is different. I actually working with a nested functions for spherical harmonics. My question (which I asked here) is only the main function of this nesting (Inside that I have two more functions). For simplification, I coded the above. Maybe I can post the real functions here now. (It might be complicated).

I would like to optimise the entire coding by avoiding the for loops wherever its possible. Here I think, only in the first function for loops can be avoided. In the function - 3 is recursion which could not be avoided.

More importantly, I would like to keep these three functions separate itself.

Could you have a rough look (or can do further, if you are intersted.) and give me an idea.

Many thanks!...

Function - 1 : Ylm_grid

function [ Y_grid ] = Ylm_grid( L_or_M , Co_latitude_grid , Longitude_grid )
tic
Y_grid {length(Co_latitude_grid), length(Longitude_grid)}=zeros; % Initialization of grid by zeros.
  for i=1:length(Longitude_grid)
      for j=1:length(Co_latitude_grid)
         Y_grid{j,i}=Ylm( L_or_M,Co_latitude_grid(j),Longitude_grid(i));
      end
  end
toc
end

Function - 2 : Ylm

function [ RSlm ] = Ylm( Degree_or_Order , Co_latitude , Longitude )                     
  P_mm_theta  = Plt(Degree_or_Order,Co_latitude);        
  Rlm = P_mm_theta.*repmat( cosd( [0:1:Degree_or_Order]*Longitude ), [Degree_or_Order+1, 1] );  
  Slm = rot90(P_mm_theta.*repmat( sind( [0:1:Degree_or_Order]*Longitude ), [Degree_or_Order+1, 1] ),2); 
  Slm_new = [zeros(Degree_or_Order+1, 1),  Slm(:,1:Degree_or_Order)];    
  RSlm = Rlm  + Slm_new ;
end

Function - 3 : Plt

function [ P_bar ] = Plt( Degree_or_Order, Theta_In_Degree )
if (Degree_or_Order > 0 || Degree_or_Order ==0)  % Program starts only when Degree is greater than Order
  grid_length=length(Theta_In_Degree);
  Cos(1,1,1:grid_length)=cosd(Theta_In_Degree);     
  Sin(1,1,1:grid_length)=sind(Theta_In_Degree); 
  P_bar(1,1,1:grid_length)=1/sqrt(4*pi) ; % P_bar(0,0) initialization
  r=0;   
 for m=0:Degree_or_Order         % The two for loops fill the matrix from M=0 to M
          format long e;         %                                    -column-wise-
          if m==0;  del=1;         
          else     del=0;                             
          end                  
          for L=r:Degree_or_Order        
              if and(L==m,m<Degree_or_Order) % Check the diagonal elements and assign values to array.       
                  P_bar(m+2,m+2,1:grid_length)=Sin* sqrt( (((2*m)+3)*(1+del)) /...
                                       ((2*m)+2) ).*P_bar(m+1,m+1,1:grid_length);
              elseif L==(m+1) % Check the lower diagonal elements and assign values to array. 
                  P_bar(m+2,m+1,1:grid_length)=Cos*sqrt((2*m)+3).*P_bar(m+1,m+1,1:grid_length);
              elseif and(L>m,L~=(m+1)) % Fills rest of the  values to array.
                  P_bar(L+1,m+1,1:grid_length)=( Cos* sqrt( (((2*(L-1))+3)*((2*(L-1))+1)) / (((L-1)-m+1)*((L-1)+m+1)) ).*P_bar(L,m+1,1:grid_length) )...
                                       -( sqrt( (((2*(L-1))+3)*(((L-1)^2)-(m*m))) / (((2*(L-1))-1)*((((L-1)+1)^2)-(m*m))) ) * P_bar(L-1,m+1,1:grid_length) );
              end
          end
          r=r+1;
 end
else    
  fprintf('\t\t\tError!\n\tPlease enter +ve number!\n'); 
end
end

[Stephen23]

Interesting task... and not an easy question for us to answer. The first question: is it worth trying to optimize the code? If it takes a few hours of re-coding in order to save ten minutes of computation, then is this a good use of time? How much faster do you expect the code to be able to run?

One way that you might be able to speed it up is to consider removing any operations from inside the loops that are always the same or that are applied to the entire output of the loops, and placing these before/after the loops. In other words: pre-compute everything possible before the loops. This would be quite a fundamental change in how the code works, and might make it less intuitive to follow.