How can I use vectorization to speed up the following code

I have the following code. x, y and t are variables. M is a big number and I discovered that the double loop will take a longer time. I need help in making it faster, all other example have seen are partially helpful. Thanks a lot.

function f=exact(x,y,t,M,j)

f=0;

    for i=1:M
         for j=1:M
             f=f+(16*(sin(j*pi)*sin(i*pi)*l*i^2*pi^3+4*sin(j*pi)*cos(i*pi)*j*i*pi^2+2*sin(i*pi)*cos(j*pi)*i^2*pi^2+2*sin(j*pi)*j*i*pi^2-2*i^2*pi^2*sin(i*pi)...
                 -6*sin(j*pi)*sin(i*pi)*j*pi+8*cos(i*pi)*cos(j*pi)*i*pi-8*i*pi*cos(i*pi)+4*cos(j*pi)*i*pi-12*sin(i*pi)*cos(j*pi)-4*i*pi+12*sin(j*pi)))/(i^4*pi^7*j^3)...
                   *cos(sqrt(i^2+j^2)*pi*t).*sin(i*pi*x).*sin(j*pi*y);

         end
    end

[Jan]

The expected size of M matters: Completely different approaches are required for M=4 and M=4*1e7 . "Faster" ist not clear enough: Do you have a fixed time limit, such that it is worth to spend an hour to create a C-Mex function? Or doy you want to save some minutes only, such that spending an hour would mean a loss of total time?

The input "j" is not used, because it is overwritten by the FOR loop counter. So perhaps this is a typo: The variable "l" is undefinde. Then accelerating buggy code is fruitless.

[Jan]

The formula looks very strange: You evaluate a lot of "sin(i*pi)" and "cos(j*pi)", but "i" and "j" are intergers. Then rounding errors dominate the result: Mathematically "sin(i*pi)" is exactly zeros, but the limited precision leads to a value of 1.2246e-016. Nevertheless, you can omit all terms, which contain a "sin(i*pi)" or "sin(j*pi)", while the cos-expressions can be replaced by a multiplication by 1 - which can be omitted.

I have the strong impression, that you did not post the real code. But optimizing pseudo-code is meaningless. Please provide a valid input also.

Answer 1

[Jan]

Although I'm in doubt that the shown code is correct, I've cleaned it up a little bit moving repeated calculations outside the inner loop:

function f=exact(x,y,t,M,l)  % "l" instead of "j"?!

pi3 = pi^3;
pi2 = pi^2;
pi7 = pi^7;

f = 0;
for i=1:M
   i2 = i^2;
   sin_i_pi = sin(i * pi);
   cos_i_pi = cos(i * pi);
   sin_i_pi_x = sin(i * pi * x);

     C1 = sin_i_pi * i2;
     C2 = 2 * i2 * pi2 * sin_i_pi;
     C3 = 8 * i * pi * cos_i_pi;
     C4 = 4 * i * pi;
     C5 = i^4 * pi7;
     C6 = 4 * cos_i_pi * i * pi2;
     C7 = 8 * cos_i_pi * i * pi;
     C8 = 6 * sin_i_pi * pi;

     for j=1:M
        sin_j_pi = sin(j * pi);
        cos_j_pi = cos(j * pi);

        f = f + (sin_j_pi * C1 * l * pi3 + ...
           sin_j_pi * j * C6 + ...
           2 * C1 * cos_j_pi * pi2 + ...
           2 * sin_j_pi * j * i * pi2 - ...
           C2 - ...
           sin_j_pi * j * C8 + ...
           cos_j_pi * C7 - ...
           C3 + ...
           4 * cos_j_pi * i * pi - ...
           12* sin_i_pi * cos_j_pi - ...
           C4 + 12 * sin_j_pi) / ...
           (C5 * j^3) * ...
           cos(sqrt(i2 + j^2) * pi * t) .* sin_i_pi_x .* sin(j*pi*y);
     end
  end

f = f * 16;

This takes 1.10 sec instead of 3.15 sec on my old Matlab R2011b/64/Win7/Core2Duo.

But as written in my comment already: All sin terms can be omitted, because they are 0, while the cos expression is 1, such that the formula can be simplified substantially.

[Oluwaseun Lijoka]

Dear Jan, Thanks for spotting the error in the function. The correct function declaration is function f=exact(x,y,t,M). Thanks a lot.

[Jan]

The "l" is undefined in the first term:

sin(j*pi) * sin(i*pi) * l * i^2 * pi^3

What about the vanishing value of "sin(i*pi)"?

[Oluwaseun Lijoka]

Dear Jan, The code is an exact solution to a solved wave equation in which when M is tending to infinity, that is why I stated M to be a big number. I want to compare this exact solution with a numerical approximation in terms of error and convergence. Looping i and j over M each time could be waste of time. Here is the whole code:

function f=exact(x,y,t,M,j)

if (nargin<5) j=0; end f=0; % tic if (j==0) % exact soln for m=1:M for l=1:M f=f+((16*(sin(l*pi)*sin(m*pi)*l*m^2*pi^3+4*sin(l*pi)*cos(m*pi)*l*m*pi^2+2*sin(m*pi)*cos(l*pi)*m^2*pi^2+2*sin(l*pi)*l*m*pi^2-2*m^2*pi^2*sin(m*pi)... -6*sin(l*pi)*sin(m*pi)*l*pi+8*cos(m*pi)*cos(l*pi)*m*pi-8*m*pi*cos(m*pi)+4*cos(l*pi)*m*pi-12*sin(m*pi)*cos(l*pi)-4*m*pi+12*sin(m*pi)))/(m^4*pi^7*l^3))... .*cos(sqrt(m^2+l^2)*pi*t).*sin(m*pi*x).*sin(l*pi*y);

         end
    end
elseif (j==1)  % derivative w.r.t time

    for m=1:M
        for l=1:M
            f=f+((16*(sin(l*pi)*sin(m*pi)*l*m^2*pi^3+4*sin(l*pi)*cos(m*pi)*l*m*pi^2+2*sin(m*pi)*cos(l*pi)*m^2*pi^2+2*sin(l*pi)*l*m*pi^2-2*m^2*pi^2*sin(m*pi)...
                 -6*sin(l*pi)*sin(m*pi)*l*pi+8*cos(m*pi)*cos(l*pi)*m*pi-8*m*pi*cos(m*pi)+4*cos(l*pi)*m*pi-12*sin(m*pi)*cos(l*pi)-4*m*pi+12*sin(m*pi)))/(m^4*pi^7*l^3))...
                   *(-sqrt(m^2+l^2)*pi).*sin(sqrt(m^2+l^2)*pi*t).*sin(m*pi*x).*sin(l*pi*y);
        end
    end

elseif (j==2) % derivative w.r.t variable x for m=1:M for l=1:M f=f+((16*(sin(l*pi)*sin(m*pi)*l*m^2*pi^3+4*sin(l*pi)*cos(m*pi)*l*m*pi^2+2*sin(m*pi)*cos(l*pi)*m^2*pi^2+2*sin(l*pi)*l*m*pi^2-2*m^2*pi^2*sin(m*pi)... -6*sin(l*pi)*sin(m*pi)*l*pi+8*cos(m*pi)*cos(l*pi)*m*pi-8*m*pi*cos(m*pi)+4*cos(l*pi)*m*pi-12*sin(m*pi)*cos(l*pi)-4*m*pi+12*sin(m*pi)))/(m^4*pi^7*l^3))... *(m*pi).*cos(sqrt(m^2+l^2)*pi*t).*cos(m*pi*x).*sin(l*pi*y); end

    end

elseif (j==3) % derivative w r t variable y

    for m=1:M
         for l=1:M
             f=f+((16*(sin(l*pi)*sin(m*pi)*l*m^2*pi^3+4*sin(l*pi)*cos(m*pi)*l*m*pi^2+2*sin(m*pi)*cos(l*pi)*m^2*pi^2+2*sin(l*pi)*l*m*pi^2-2*m^2*pi^2*sin(m*pi)...
                 -6*sin(l*pi)*sin(m*pi)*l*pi+8*cos(m*pi)*cos(l*pi)*m*pi-8*m*pi*cos(m*pi)+4*cos(l*pi)*m*pi-12*sin(m*pi)*cos(l*pi)-4*m*pi+12*sin(m*pi)))/(m^4*pi^7*l^3))...
                   *(l*pi).*cos(sqrt(m^2+l^2)*pi*t).*sin(m*pi*x).*cos(l*pi*y);

         end
    end

end Thanks a lot