Simplify nested loops with condition-CodePudding

I have a matlab program with 5 nested

for

loops and a

if

condition like this:

for x0=1:N
    for y0=1:N
        for k=1:N
            for x1=1:N
                for y1=1:N
                    if ~((y1-x1>N/2)||(x1-y1>N/2)) && ~((y0-x0>N/2)||(x0-y0>N/2))
                        A(x0,y0)=A(x0,y0) 2^(k*((x0-y0) (x1-y1)))*B(x1,y1)
                    end
                end
            end
        end
    end
end

where A and B are two matrices. How can I make this program run faster?

I've tried to use meshgrid but it seems doesn't work because there's a

if

condition.

CodePudding user response：

Lets be smart about loops and conditions first, as you are using the loop indices as condition variables.

We start with ~(y1-x1>N/2)||(x1-y1>N/2), or way clearer, abs(y1-x1)<N/2.

Instead of having an if condition, why not enforce y1 to be in range, always?

The last loop can be written as y1=max(x1-N/2,1):min(x1 N/2,N), and thus the entirety of the first part of the if condition is not needed. We can do the same for the other variables, of course:

for x0=1:N
    for y0=max(x0-N/2,1):min(x0 N/2,N)
        for k=1:N
            for x1=1:N
                for y1=max(x1-N/2,1):min(x1 N/2,N)
                     A(x0,y0)=A(x0,y0) 2^(k*((x0-y0) (x1-y1)))*B(x1,y1)
                end
            end
        end
    end
end

Now, for clarity, lets reshuffle and vectorize that k. There is no need for it to be the middle loop, in fact, its only feature as the middle loop is to confuse the person reading the code. But aside from that, there is no need for it to be a loop either.

k=1:N;
for x0=1:N
    for y0=max(x0-N/2,1):min(x0 N/2,N)
         for x1=1:N
             for y1=max(x1-N/2,1):min(x1 N/2,N)
                  A(x0,y0)=A(x0,y0) sum(2.^(k*((x0-y0) (x1-y1))))*B(x1,y1)
             end
         end
    end
end

Now, is this faster?

No. MATLAB is really good at optimizing your code, so it is not faster. But at least its way way clearer, so I guess you got that going for you. But you need it faster! Well.... I am not sure you can. You have a 5 nested loops, that is just super slow. I don't think you can meshgrid this, even without the conditions, because you intermingle all loops. meshgrid is good when well, you do operations on a mesh grid, but in your case you use all x1,y1 for every x0,y0 and thus its not a mesh operation.

CodePudding user response：

Update this solution may not be faster under Matlab, because the execution engine can optimise the loops in the original code. (It does provide a speedup under Octave.)

One trick to deal with if statements within loops is to turn the if statement (or part of it) into a logical matrix. You can then multiply the logical matrix elementwise by the matrix of values you are adding in each step. A false value will evaluate to zero and will not change the result.

This only works if each element can be computed independently of the others. It will generally make the actual calculation line slower, but in Matlab this is often outweighed by the huge speed improvement from the the removal of the for loops.

In your example, we can use this idea with meshgrid to remove the loops over x0 and y0. The calculation line needs to become an elementwise matrix caluclation, so elementwise operators .*, .^ and | replace *, ^ and |.

% Warning: Y0 and X0 are swapped in this example
[Y0, X0] = meshgrid(1:N,1:N);

% Create a logical matrix which represents part of the if statement
C = ~((Y0-X0>N/2) | (X0-Y0>N/2));

for k=1:N
    for x1=1:N
        for y1=1:N
           if ~((y1-x1>N/2)||(x1-y1>N/2))
               % Multiply by C elementwise
               A = A   C.*2.^(k*((X0-Y0) (x1-y1)))*B(x1,y1);
           end
        end
    end
end

You could even take this concept further and remove more loops using multidemensional ndgrids, but it becomes more complex (you have to start summing over dimensions) and the multidimensional arrays become unwieldy if N is large.

Note: be careful with index order. meshgrid defines y as rows and x as columns, so matrix indexing is A(y,x) but you are using A(x,y). So to make your example work I've switched x and y in the output of meshgrid.