Acceleration of code by replacing two nested for-loops

Hello,

I am trying to accelerate my code, where I calculate the matrix containing the output image 'Out' by using two nested for loops as part of a simulation for imaging with a spatially varying point spread function. Basically, 'Out' is the sum of each pixel (ix,iy) of the input image 'A' multiplied with an image 'B(:,:,ix,iy)' corresponding to the point spread function of that pixel. The code to calculate that looks something like this (I defined A and B as zeros here just to show their size):

x_max = 256;
y_max = 256;
A = zeros(x_max,y_max);
B = zeros(x_max,y_max,x_max,y_max);
Out = zeros(x_max,y_max);
for ix=1:x_max
    for iy=1:y_max
        Out = Out + A(ix,iy)*B(:,:,ix,iy);
    end
end

My question is, can this code be accerated by replacing these nested for-loops by using matrix multiplication?

Many thanks,

Tom

[Aquatris]

Something does not feel right with your code. In here you are adding 1x1 matrix with a 256x256x1x1 matrix and rewrite the whole Out variable. Is that what you actually are trying to do?

Out = Out + A(ix,iy)*B(:,:,ix,iy);

[Torsten]

@Aquatris

It's an addition of matrices that have the form a_ij*B_ij.

[Tom Glosemeyer]

Yes, Torsten is right. The 2D output image is the sum of all the images obtained by multiplying each pixel (ix,iy) of A with the corresponding 2D-image B(:,:,ix,iy).

[Aquatris]

@Torsten Ohh I see. Thanks for pointing it out Torsten.

Answer 1

[Star Strider]

This could be difficult with your code because of the multidimensional nature of ‘B’, however the usual approach is to create matrices for the different arguments and then use those in the function, rather than creating nested loops.

To illustrate —

x = linspace(-5, 5, 20);
y = linspace(0,  10, 30);
z = @(x,y) exp(-(x.^2/10 + (y-5).^2/5));

[X, Y] = ndgrid(x,y);
Z = z(X,Y)*5;                                       % No (Explicit) Loops!

figure
surf(X, Y, Z)
colormap(turbo)
colorbar
xlabel('X')
ylabel('Y')
axis('equal')

.