Using parfor to evaluate integrations

Probably a stupid question.

I want to integrate a function f(x,y). Can I use parfor() as shown below to compute the sum? I am confused because I consider that each loop over one "ix" value in the following parfor() is run independently from other "ix" values, which means that I must obtain a different "sum_f" for each "ix" value. Right?

clear; clc;

xmin = -2; xmax = 1;
ymin = 1; ymax = 3;

dx = 0.001;
dy = dx;

xs = xmin:dx:xmax; Nx = length(xs);
ys = ymin:dy:ymax; Ny = length(ys);


f = @(x,y) x^2+y^2;

sum_f = 0;
parfor ix = 1:Nx
    x = xs(ix);
    for iy = 1:Ny
        y = ys(iy);
        sum_f = sum_f + f(x,y)*dx*dy;
    end
end

I = sum(sum_f(:))

[Dyuman Joshi]

Any particular reason why you are using a double for loop, instead of vectorizing the function handle and the sum?

[Luqman Saleem]

@Dyuman Joshi umm.. I don't get you. I am not familiar with "vectorizing" the functions. Do you mean rewriting f(x,y) as

f = @(x,y) arrayfun(@(x,y) f(x,y), x, y)

Answer 1

[Dyuman Joshi]

See - Array vs Matrix Operations

xmin = -2; xmax = 1;
ymin = 1; ymax = 3;

dx = 0.001;
dy = dx;

xs = xmin:dx:xmax; Nx = length(xs);
ys = ymin:dy:ymax; Ny = length(ys);

%vectorizing the function handle
f = @(x,y) x.^2+y.^2;


tic
sum_f = 0;
for ix = 1:Nx
    x = xs(ix);
    for iy = 1:Ny
        y = ys(iy);
        sum_f = sum_f + f(x,y)*dx*dy;
    end
end

I1 = sum(sum_f(:))

I1 = 32.0317

toc

Elapsed time is 0.364607 seconds.

tic
I2 = sum(f(xs,ys.').*dx.*dy, 'all')

I2 = 32.0317

toc

Elapsed time is 0.033681 seconds.

%Using tolerance to compare floating point numbers
abs(I1-I2)<1e-10

ans = logical
   1

[Luqman Saleem]

@Dyuman Joshi :-O areee bhai.. that's much more faster. wow. thank you so much

[Mike Croucher]

Ahhh loops vs vectors....my old friend. We meet again.

I'm going to let you in on a secret......sometimes loops are faster!

• Original loop with a function call: 0.218 seconds
• Optimised loop with function call removed: 0.025 seconds
• vectorised loop: 0.0275 seconds

So here, my loop version is slightly faster than the vectorised version. What you see might be dependent on machine, problem size and MATLAB version

Observe:

xmin = -2; xmax = 1;
ymin = 1; ymax = 3;
dx = 0.001;
dy = dx;
xs = xmin:dx:xmax; Nx = length(xs);
ys = ymin:dy:ymax; Ny = length(ys);
%vectorizing the function handle
f = @(x,y) x.^2+y.^2;
tic
sum_f = 0;
for ix = 1:Nx
    x = xs(ix);
    for iy = 1:Ny
        y = ys(iy);
        sum_f = sum_f + f(x,y)*dx*dy;
    end
end
I1 = sum(sum_f(:))

I1 = 32.0317

toc

Elapsed time is 0.218968 seconds.

disp('Inline the function call in the loop')

Inline the function call in the loop

xmin = -2; xmax = 1;
ymin = 1; ymax = 3;
dx = 0.001;
dy = dx;
xs = xmin:dx:xmax; Nx = length(xs);
ys = ymin:dy:ymax; Ny = length(ys);

tic
sum_f = 0;
for ix = 1:Nx
    x = xs(ix);
    for iy = 1:Ny
        y = ys(iy);
        sum_f = sum_f + (x.^2+y.^2)*dx*dy;
    end
end
I1 = sum(sum_f(:))

I1 = 32.0317

toc

Elapsed time is 0.025205 seconds.

disp('vectorised version')

vectorised version

tic
I2 = sum(f(xs,ys.').*dx.*dy, 'all')

I2 = 32.0317

toc

Elapsed time is 0.027567 seconds.