how to vectorize a for-loop or use parfor to speed it up?

I have this for loop right now:

     for x = 2:Nx_max-1
          for y = 2:Ny_max-1
              lp(x,y,t1) = (1/(2*h^2))*dxy(x,y)*p(x+1,y+1,t1) + ...
              (dxx(x,y)/(h^2) - vx(x,y)/(2*h) - (x-Nx-1)/(2*gamma))* p(x+1,y,t1) - ...
              (1/(2*h^2))*dxy(x,y)*p(x+1,y-1,t1) + (dyy(x,y)/(h^2) - ...
              vy(x,y)/(2*h) - (y-Ny-1)/(2*gamma))*p(x,y+1,t1) - ...
              (2*dxx(x,y) + 2*dyy(x,y))*(1/h^2)*p(x,y,t1) + ...
              (dyy(x,y)/(h^2) + vy(x,y)/(2*h) +(y-Ny-1)/(2*gamma))*p(x,y-1,t1) - ...
              (dxy(x,y)/(2*h^2))*p(x-1,y+1,t1) +(dxx(x,y)/(h^2) +vx(x,y)/(2*h) +...
              (x-Nx-1)/(2*gamma))*p(x-1,y,t1) + (dxy(x,y)/(2*h^2))*p(x-1,y-1,t1);

          end
      end

I asked a question previously where I didn't actually write out the full function here, and I was recommended to vectorize it for speed by replacing x by 2:(Nx_max-1) and y by 2:(Ny_max-1). This is causing problems because it looks like the multiplications would cause it to be matrix multiplication instead of multiplying the specific elements of the matrix that I'm interested in mylitplying together. And I also have a section where there ends up being a vector (2:(Nx_max-1)) multiplying a matrix. I don't want the end result to be a vector (as it would in regular matlab notation), I want each element of the matrix in that case to be multiplied by its x value.

Is there a way to vectorize this or somehow speed up the process?

I had thought about parallel computation, but it looks like I can't nest parfor loops.

Any suggestions would be appreciated

Answer 1

[Jan]

An elementwise multiplication is performed by the .* operator, while * is a matrix multiplication. This should solve your problem already. Please post your trial to vectorize the code instead of letting us create it from scratch.

An optimization of code depends on the input values. It matters if Nx_max is 10 or 1e7. It matters if dxy and gamma are arrays or functions. So please post a version of your code, which can be started by copy&paste. Otherwise suggestions contains a lot of guessing.

But if guessing is required, consider the standard tips fopr loop optimization at first:

• avoid repeated calculations
• pre-allocate the output

E.g. (1/(2*h^2)) is calculated nearly Ny_max * Nx_max times. So better create a variable before the loop to carry the result. So befpre creating the vectorized version, measure the speed with a cleaner version like this:

hs =  h^2;
h2  = 2 * h;
hs2 = 2 * hs;
c1  = 1/hs2;
g2  = 2 * gamma;

lp = zeros(Nx_max-1, Ny_max-1, t1);
for x = 2:Nx_max-1
  c2 = (x-Nx-1)/g2;
  for y = 2:Ny_max-1
    lp(x,y,t1) = c1 * dxy(x,y) * p(x+1,y+1,t1) + ...
       (dxx(x,y) / hs - vx(x,y) / h2 - c2) * p(x+1,y,t1) - ...
       c1 * dxy(x,y) * p(x+1,y-1,t1) + ...
       (dyy(x,y)/hs - vy(x,y)/h2 - (y-Ny-1)/g2) * p(x,y+1,t1) - ...
       (2 * dxx(x,y) + 2 * dyy(x,y)) * p(x,y,t1) / hs + ...
       (dyy(x,y)/hs + vy(x,y)/h2 +(y-Ny-1)/g2) * p(x,y-1,t1) - ...
       dxy(x,y) / hs2 * p(x-1,y+1,t1) + ...
       (dxx(x,y)/hs + vx(x,y) / h2 + c2) * p(x-1,y,t1) + ...
       dxy(x,y) / hs2 * p(x-1,y-1,t1);
  end
end

Then, in a next step, replace the loop over x by replacing all "x" by "(2:Nx_max-1)" and the operators "*" and "/" by their elementwise versions ".*" and "./". Test, if the result is still not changed except for rounding errors. Fix "(2:Nx_max-1)+1" to the faster "3:Nx_max", because in the 2nd version the index vector is not calculated exlicitly.

Finally you can try to replace the inner loop also, but perform a speed test afterwards. A complete vectorization is not necessarily the fastest version.

[sasha]

all things are arrays except for the term that has just a value of x or y. For the vectorization I just replaced x and y by a vector 2:Nx_max-1. The elementwise multiplication didn't work because most of my terms are 2D arrays, but there is one term (the one where it's just x or y) that is a column vector. I can't do elementwise multiplication on that.

I have already pre-allocated all of the arrays.

Nx_max and Ny_max are both 201. It is this large size that slows it down a lot.

I will try replacing the h^2 and 2*h expressions for values. I did not think that this would be the cause of the large time lag. Thanks

[Jan]

I still do not understand, why the elementwise mutliplications "did not work". Please post your code and the corresponding error message. Actually replacing the outer loop by inserting the index vector should be trivial. And I'd actually did this for you, if you post a running example, such that I could debug my suggestion.

[sasha]

element-wise did not work because x - Ny-1 would be written as a vector of size Nx_max-2 by 1. While all other terms are Nx_max-2 by Ny_max-2. This gives an error stating that matrix dimensions do not match.

There are 4 of these problematic terms, the first of which is in the second line of the definition of my new matrix (the last term of that line).

Here are the two things I tried for that specific term (ignoring the 1/(2*gamma) because that is just a constant and is not causing any problems):

((2:(Nx_max-1)) -Nx-1) * p((2:(Nx_max-1))+1, (2:(Ny_max-1)),t1)

The problem with this isn't that it gave me an error message. What it did was give me a column vector (not what I want) my result should be a matrix of the same dimensions as p. What I want to happen is that for all the x=2 terms of p to be multiplied by (2-Nx-1), for all the x=3 terms of p to be multiplied by (3-Nx-1), for all the x=4 terms of p to be multiplied by (4-Nx-1), etc. This should be a matrix with the same dimensions as p.

Also, because it ends with a vector it can not be added to any of the other terms that end up as matrices (because they need to be matrices).

Then I tried:

x1 = 2:(Nx_max-1);
y1 = 2:(Ny_max-1);

x1 .* p(:,:,t1) - (Nx+1) * p(:,:,t1)

The second part of this term is fine. But the first part gives an error message saying that matrix dimensions need to match. As you can see, they don't because x1 is a vector and p is a matrix. Thus element wise multiplication can't happen because the number of elements do not match.

Does this make more sense now?

Note: the actual expressions for vx, vy, p, dxx, dxy, and dyy don't really matter. What is important is that they are all Nx_max by Ny_max matrices (and then I only want to deal with the bulk elements, ie none of the edges of the matrices) so I work with Nx_max-2 by Ny_max-2 values.

I know that the for loop gives me the correct answer, but it is slow. I have been testing the new attempts against the values of the for loop. Basically running both attempts at the same time (giving them different names) and then seeing if they give the same result by subtracting them from one another.

[Jan]

Here a guessed version of a vectorized outer loop:

hs =  h^2;
h2  = 2 * h;
hs2 = 2 * hs;
g2  = 2 * gamma;
Nx1 = Nx_max - 1;
Nx2 = Nx_max - 2;
c2  = ((2 - Nx-1):(Nx_max - 2 - Nx)) / g2;
lp  = zeros(Nx1, Ny_max-1, t1);
for y = 2:Ny_max-1
  c3 = dxy(2:Nx1,y) ./ hs2;
  c4 = dxy(2:Nx1,y) ./ hs2;
  c5 = vy(2:Nx1,y)  ./ h2;
  c6 = vx(2:Nx1,y)  ./ h2;
  c7 = dyy(2:Nx1,y) ./ hs;
  lp(2:Nx1, y, t1) = ...
    c3 .* p(3:Nx_max,y+1,t1) + ...
    (dxx(2:Nx1,y) ./ hs - c6 - c2) .* p(3:Nx_max,y,t1) - ...
    c3 .* p(3:Nx_max,y-1,t1) + ...
    (c7 - c5 - (y-Ny-1) ./ g2)    .* p(2:Nx1,y+1,t1) - ...
    (dxx(2:Nx1,y) + dyy(2:Nx1,y)) .* p(2:Nx1,y,t1) .* (2 ./ hs) + ...
    (c7 + c5 + (y-Ny-1) ./ g2)    .* p(2:Nx1,y-1,t1) - ...
    c4 .* p(1:Nx2,y+1,t1) + ...
    (dxx(2:Nx1,y) ./ hs + c6 + c2) .* p(1:Nx2,y,t1) + ...
    c4 .* p(1:Nx2,y-1,t1);
end

Unfortunately I'm in doubt that this produces still the same results, because this kind of editing is prone to typos. But you see the strategy and can replace the repeated calculations step by step by your own and check, if teh results are not changed.