how to reduce computation time for nested loops

[N,M] = size(fft_img_test(:,:,1));
B = zeros(N,M,N,M);

for k=1:No_frames
for vy = 1:N
    for vx = 1:M
        for uy = 1:N
            for ux = 1:M
                
                u = uy+vy-2; v = ux+vx-2;
                u = mod(u,N); v = mod(v,M);
                u = u + 1; v = v + 1;
              
                B(uy,ux,vy,vx) =  B(uy,ux,vy,vx) + fft_img_test(uy,ux,k) * fft_img_test(vy,vx,k)...
                    * conj(fft_img_test(u,v,k)); 
            end
        end
    end
end
end
phase_B = (angle(B/k));

[Adam Danz]

It's easier if we can actually run the code but there are variable values missing that is preventing this. If you attach a mat file containing all needed variables to run the code, it would be helpful.

[abdelelah alzahed]

I have attached the m-files that are related to my question.The nested loops are presented in the Bispec file.

[Adam Danz]

You've got more than 536.8 million iterations (M * N * M * N * No_frames = 536870912).

Vectorization sometimes speeds up computation time but it's not always faster than loops and the expansion that may be required to do these computations using vectorization may exceed memory limits.

Could you explain in words the logic of this section below? It's not making much sense to me. When vx and vy equal 1, u and v are the same values as ux and uy. But then it starts to change when vx and vy are no longer 1.

u = uy+vy-2;
u = mod(u,N);
v = ux+vx-2;
v = mod(v,M);
u = u + 1;
v = v + 1;

conj(fft_img_test(u,v,k))

[abdelelah alzahed]

ux and uy or vx and vy are two coordinates representing spatial frequencies of an image. Here, u and v are optimized coordinates for the calculation of the last term (the conjugate of an image) where v = ux+vx and u = uy+vy. They are modified as indicated above to maintain the sum within the maximum number of spatial frequency points N or M. Therefore, once u or v are greater than N and M, respectively, the command lines, as stated above, will restart the calculation of the conjugate term to the starting point (1,1). This calculation is also known as the image bispectrum.

[Adam Danz]

"Therefore, once u or v are greater than N and M, respectively, the command lines, as stated above, will restart the calculation of the conjugate term to the starting point (1,1). "

Then shouldn't u and v always equal uy and ux (which isn't the case)?

Answer 1

[Adam Danz]

Your loops sum to over 536.8 million iterations but due to the size of your variables (B has >268.4 million elements) the expansion that would be required to vectorize this would very likely exceed memory capacity.

One thing I changed below to slightly speed up the code is to move the complex conjugate calculations outside of the loop and index the results rather than performing the complex conjugate within the loop.

I've also added a multi-loop waitbar found on the file exchange (here) that will show you the progress of each loop so you don't have to wonder about the progress of the execution.

function [phase_B]=Bispec(fft_img_test, No_frames)

[N,M] = size(fft_img_test(:,:,1));
B = zeros(N,M,N,M);
multiWaitbar('k-loop',0);
multiWaitbar('vy-loop',0);
multiWaitbar('vx-loop',0);
multiWaitbar('uy-loop',0);
multiWaitbar('ux-loop',0);

fft_img_testConj = conj(fft_img_test); 
for k=1:No_frames
    multiWaitbar( 'k-loop', k/No_frames);
    for vy = 1:N
        multiWaitbar( 'vy-loop', vy/N);
        for vx = 1:M
            multiWaitbar( 'vx-loop', vx/M);
            for uy = 1:N
                multiWaitbar( 'uy-loop', uy/N);
                for ux = 1:M
                    multiWaitbar( 'ux-loop', ux/M);
                    u = uy+vy-2;
                    u = mod(u,N);
                    v = ux+vx-2;
                    v = mod(v,M);
                    u = u + 1;
                    v = v + 1;
                    
                    B(uy,ux,vy,vx) =  B(uy,ux,vy,vx) + fft_img_test(uy,ux,k) * fft_img_test(vy,vx,k) * fft_img_testConj(u,v,k);
                    % Show current indices
                    % fprintf('uy,ux,vy,vx,u,v = %d %d %d %d %d %d\n', uy,ux,vy,vx,u,v);
                end
            end
        end
    end
end
multiWaitbar('closeall')
phase_B = angle(B/No_frames);  % old version: (angle(B/k))

I played around with some vectorization which does speed up the code but since I'm not certain about how the last term using the u and v indices fit in, I stopped here.

fft_img_testConj = conj(fft_img_test);
for k=1:No_frames
    multiWaitbar( 'k-loop', k/No_frames);
    for vy = 1:N
        multiWaitbar( 'vy-loop', vy/N);
        for vx = 1:M
            multiWaitbar( 'vx-loop', vx/M);
            
            % Vectorization attempt
            u = (1:N)+vy-2;
            u = mod(u,N);
            v = (1:N)+vx-2;
            v = mod(v,M);
            u = u + 1;
            v = v + 1;
            B(:,:,vy,vx) =  B(:,:,vy,vx) + fft_img_test(:,:,k) .* fft_img_test(vy,vx,k) .* fft_img_testConj(u,v,k);
            
        end
    end
end

[abdelelah alzahed]

Thank you Adam Danz !