The fastest method to perform sum of function

Suppose that I have a function fun(x,y,z) whose output is a 2 by 2 matrix (one example is given below). I want to calculate the sum of this function over a 3D grid of x, y, and z. Of course I can use for() loops, but that will be very slow. What is the fastest way to perform this task? Vectorizing, but how?

clear; clc;

dx = 0.1;
dy = 0.15;
dz = 0.05;


xmin = -2;
xmax = 4;

ymin = -2*pi;
ymax = 2*pi;

zmin = -1;
zmax = 2;


xs = xmin:dx:xmax;
ys = ymin:dy:ymax;
zs = zmin:dz:zmax;


answer = sum(fun(xs,ys,zs),'All')

Arrays have incompatible sizes for this operation.

Error in solution>fun (line 33)
hx = [0,1; 1,0]*(-J*S*(cos(y/2 - (3^(1/2)*x)/2) + cos(y/2 + (3^(1/2)*x)/2) + cos(y)));

% answer should be a 2-by-2 matrix

function out = fun(x,y,z)
J = 1;
S = 1;
D = 0.75;
Sigma0_R = -1i*0.1*eye(2);

h0 = eye(2)*(3*J*S);
hx = [0,1; 1,0]*(-J*S*(cos(y/2 - (3^(1/2)*x)/2) + cos(y/2 + (3^(1/2)*x)/2) + cos(y)));
hy = [0, -1i; 1i, 0]*(-J*S*(sin(y/2 - (3^(1/2)*x)/2) + sin(y/2 + (3^(1/2)*x)/2) - sin(y)));
hz = [1, 0; 0, -1]*(-2*D*S*(sin(3^(1/2)*x) + sin((3*y)/2 - (3^(1/2)*x)/2) - sin((3*y)/2 + (3^(1/2)*x)/2)));

out = inv( z*eye(2) - (h0 + hx + hy + hz) - Sigma0_R );
end

[Stephen23]

"but that will be very slow"

Why?

[Luqman Saleem]

@Stephen23 because I will be taken thousands of points between x,y,z ranges. that will make for loops very slow

[Walter Roberson]

I want to calculate the sum of this function for all values between some ranges of x, y, and z.

Then you want to create a 3D grid of x y z values to calculate over, rather than having three vectors (of different lengths)

[Luqman Saleem]

@Walter Roberson yes yes you are right. sorry, English is not my first language.

[Matt J]

Suppose that I have a function fun(x,y,z) whose output is a 2 by 2 matrix... I want to calculate the sum of this function over a 3D grid of x, y, and z.

Does that mean the final result will also be a 2x2 matrix, or will it be a scalar representing the sum over all 2*2*Nx*Ny*Nz values?

Answer 1

[Walter Roberson]

The fastest way of adding the values depends upon the fine details of your hardware, and on the details of how the data is laid out in memory, and depends upon what else is happening on your system to determine how calculations are scheduled by the CPU.

These are not suitable topics for discussion here.

[Luqman Saleem]

Ah, I believe you are correct. Let me revise my question slightly: What is a more efficient alternative to using a for() loop for performing this sum?

[Walter Roberson]

vectorize your function so that it creates a 2 x 2 x length(x) x length(y) x length(z) result, and sum() 'all'

[Luqman Saleem]

When I do sum( , 'All'), it take sum of all matrix elements. I want to get a 2 x 2 matrix at the end.

[Aaron]

In that case, you'd do sum(A, [3,4,5]) to just sum over the dimensions you're storing the dimensional dependence in.

[Walter Roberson]

sum(ARRAY, [3 4 5])

Answer 2

[Torsten]

I tried precomputing the output of your function symbolically and then evaluate it for numerical input. But it seems there is no advantage over a simple loop with the numerical function called.

I don't know of a fast way to rewrite your function so that it can deal with array inputs for x,y and z.

dx = 0.1;
dy = 0.15;
dz = 0.05;
xmin = -2;
xmax = 4;
ymin = -2*pi;
ymax = 2*pi;
zmin = -1;
zmax = 2;
xs = xmin:dx:xmax;
ys = ymin:dy:ymax;
zs = zmin:dz:zmax;
tic
% Compute result in a simple loop
s=0;
for i=1:numel(xs)
    for j=1:numel(ys)
        for k=1:numel(zs)
           s = s+fun(xs(i),ys(j),zs(k));
        end
    end
end
toc

Elapsed time is 2.185238 seconds.

s

s = 
   1.0e+05 *

  -0.4528 - 1.4036i   0.0005 - 0.0010i
   0.0007 - 0.0010i  -0.3424 - 1.6394i

syms x y z
J = 1;
S = 1;
D = 0.75;
Sigma0_R = -1i*0.1*eye(2);
h0 = eye(2)*(3*J*S);
hx = [0,1; 1,0]*(-J*S*(cos(y/2 - (3^(1/2)*x)/2) + cos(y/2 + (3^(1/2)*x)/2) + cos(y)));
hy = [0, -1i; 1i, 0]*(-J*S*(sin(y/2 - (3^(1/2)*x)/2) + sin(y/2 + (3^(1/2)*x)/2) - sin(y)));
hz = [1, 0; 0, -1]*(-2*D*S*(sin(3^(1/2)*x) + sin((3*y)/2 - (3^(1/2)*x)/2) - sin((3*y)/2 + (3^(1/2)*x)/2)));
out = inv( z*eye(2) - (h0 + hx + hy + hz) - Sigma0_R )

out = 

out = matlabFunction(out)

out = function_handle with value:
    @(x,y,z)reshape([((z.*-1.0e+1+sin(sqrt(3.0).*x).*1.5e+1+sin(y.*(3.0./2.0)-(sqrt(3.0).*x)./2.0).*1.5e+1-sin(y.*(3.0./2.0)+(sqrt(3.0).*x)./2.0).*1.5e+1+3.0e+1-1i).*1.0e+1)./(z.*(6.0e+2-2.0e+1i)+cos(y./2.0-(sqrt(3.0).*x)./2.0).*cos(y./2.0+(sqrt(3.0).*x)./2.0).*2.0e+2+sin(sqrt(3.0).*x).^2.*2.25e+2+cos(y./2.0-(sqrt(3.0).*x)./2.0).*cos(y).*2.0e+2+cos(y./2.0+(sqrt(3.0).*x)./2.0).*cos(y).*2.0e+2+sin(y./2.0-(sqrt(3.0).*x)./2.0).*sin(y./2.0+(sqrt(3.0).*x)./2.0).*2.0e+2-sin(y.*(3.0./2.0)-(sqrt(3.0).*x)./2.0).*sin(y.*(3.0./2.0)+(sqrt(3.0).*x)./2.0).*4.5e+2-sin(y./2.0-(sqrt(3.0).*x)./2.0).*sin(y).*2.0e+2-sin(y./2.0+(sqrt(3.0).*x)./2.0).*sin(y).*2.0e+2+cos(y./2.0-(sqrt(3.0).*x)./2.0).^2.*1.0e+2+cos(y./2.0+(sqrt(3.0).*x)./2.0).^2.*1.0e+2+cos(y).^2.*1.0e+2+sin(y./2.0-(sqrt(3.0).*x)./2.0).^2.*1.0e+2+sin(y./2.0+(sqrt(3.0).*x)./2.0).^2.*1.0e+2+sin(y.*(3.0./2.0)-(sqrt(3.0).*x)./2.0).^2.*2.25e+2+sin(y.*(3.0./2.0)+(sqrt(3.0).*x)./2.0).^2.*2.25e+2+sin(y).^2.*1.0e+2-z.^2.*1.0e+2+sin(sqrt(3.0).*x).*sin(y.*(3.0./2.0)-(sqrt(3.0).*x)./2.0).*4.5e+2-sin(sqrt(3.0).*x).*sin(y.*(3.0./2.0)+(sqrt(3.0).*x)./2.0).*4.5e+2+-8.99e+2+6.0e+1i),((cos(y./2.0-(sqrt(3.0).*x)./2.0)+cos(y./2.0+(sqrt(3.0).*x)./2.0)+cos(y)+sin(y./2.0-(sqrt(3.0).*x)./2.0).*1i+sin(y./2.0+(sqrt(3.0).*x)./2.0).*1i-sin(y).*1i).*1.0e+2)./(z.*(6.0e+2-2.0e+1i)+cos(y./2.0-(sqrt(3.0).*x)./2.0).*cos(y./2.0+(sqrt(3.0).*x)./2.0).*2.0e+2+sin(sqrt(3.0).*x).^2.*2.25e+2+cos(y./2.0-(sqrt(3.0).*x)./2.0).*cos(y).*2.0e+2+cos(y./2.0+(sqrt(3.0).*x)./2.0).*cos(y).*2.0e+2+sin(y./2.0-(sqrt(3.0).*x)./2.0).*sin(y./2.0+(sqrt(3.0).*x)./2.0).*2.0e+2-sin(y.*(3.0./2.0)-(sqrt(3.0).*x)./2.0).*sin(y.*(3.0./2.0)+(sqrt(3.0).*x)./2.0).*4.5e+2-sin(y./2.0-(sqrt(3.0).*x)./2.0).*sin(y).*2.0e+2-sin(y./2.0+(sqrt(3.0).*x)./2.0).*sin(y).*2.0e+2+cos(y./2.0-(sqrt(3.0).*x)./2.0).^2.*1.0e+2+cos(y./2.0+(sqrt(3.0).*x)./2.0).^2.*1.0e+2+cos(y).^2.*1.0e+2+sin(y./2.0-(sqrt(3.0).*x)./2.0).^2.*1.0e+2+sin(y./2.0+(sqrt(3.0).*x)./2.0).^2.*1.0e+2+sin(y.*(3.0./2.0)-(sqrt(3.0).*x)./2.0).^2.*2.25e+2+sin(y.*(3.0./2.0)+(sqrt(3.0).*x)./2.0).^2.*2.25e+2+sin(y).^2.*1.0e+2-z.^2.*1.0e+2+sin(sqrt(3.0).*x).*sin(y.*(3.0./2.0)-(sqrt(3.0).*x)./2.0).*4.5e+2-sin(sqrt(3.0).*x).*sin(y.*(3.0./2.0)+(sqrt(3.0).*x)./2.0).*4.5e+2+-8.99e+2+6.0e+1i),((cos(y./2.0-(sqrt(3.0).*x)./2.0)+cos(y./2.0+(sqrt(3.0).*x)./2.0)+cos(y)-sin(y./2.0-(sqrt(3.0).*x)./2.0).*1i-sin(y./2.0+(sqrt(3.0).*x)./2.0).*1i+sin(y).*1i).*1.0e+2)./(z.*(6.0e+2-2.0e+1i)+cos(y./2.0-(sqrt(3.0).*x)./2.0).*cos(y./2.0+(sqrt(3.0).*x)./2.0).*2.0e+2+sin(sqrt(3.0).*x).^2.*2.25e+2+cos(y./2.0-(sqrt(3.0).*x)./2.0).*cos(y).*2.0e+2+cos(y./2.0+(sqrt(3.0).*x)./2.0).*cos(y).*2.0e+2+sin(y./2.0-(sqrt(3.0).*x)./2.0).*sin(y./2.0+(sqrt(3.0).*x)./2.0).*2.0e+2-sin(y.*(3.0./2.0)-(sqrt(3.0).*x)./2.0).*sin(y.*(3.0./2.0)+(sqrt(3.0).*x)./2.0).*4.5e+2-sin(y./2.0-(sqrt(3.0).*x)./2.0).*sin(y).*2.0e+2-sin(y./2.0+(sqrt(3.0).*x)./2.0).*sin(y).*2.0e+2+cos(y./2.0-(sqrt(3.0).*x)./2.0).^2.*1.0e+2+cos(y./2.0+(sqrt(3.0).*x)./2.0).^2.*1.0e+2+cos(y).^2.*1.0e+2+sin(y./2.0-(sqrt(3.0).*x)./2.0).^2.*1.0e+2+sin(y./2.0+(sqrt(3.0).*x)./2.0).^2.*1.0e+2+sin(y.*(3.0./2.0)-(sqrt(3.0).*x)./2.0).^2.*2.25e+2+sin(y.*(3.0./2.0)+(sqrt(3.0).*x)./2.0).^2.*2.25e+2+sin(y).^2.*1.0e+2-z.^2.*1.0e+2+sin(sqrt(3.0).*x).*sin(y.*(3.0./2.0)-(sqrt(3.0).*x)./2.0).*4.5e+2-sin(sqrt(3.0).*x).*sin(y.*(3.0./2.0)+(sqrt(3.0).*x)./2.0).*4.5e+2+-8.99e+2+6.0e+1i),((z.*1.0e+1+sin(sqrt(3.0).*x).*1.5e+1+sin(y.*(3.0./2.0)-(sqrt(3.0).*x)./2.0).*1.5e+1-sin(y.*(3.0./2.0)+(sqrt(3.0).*x)./2.0).*1.5e+1+-3.0e+1+1i).*-1.0e+1)./(z.*(6.0e+2-2.0e+1i)+cos(y./2.0-(sqrt(3.0).*x)./2.0).*cos(y./2.0+(sqrt(3.0).*x)./2.0).*2.0e+2+sin(sqrt(3.0).*x).^2.*2.25e+2+cos(y./2.0-(sqrt(3.0).*x)./2.0).*cos(y).*2.0e+2+cos(y./2.0+(sqrt(3.0).*x)./2.0).*cos(y).*2.0e+2+sin(y./2.0-(sqrt(3.0).*x)./2.0).*sin(y./2.0+(sqrt(3.0).*x)./2.0).*2.0e+2-sin(y.*(3.0./2.0)-(sqrt(3.0).*x)./2.0).*sin(y.*(3.0./2.0)+(sqrt(3.0).*x)./2.0).*4.5e+2-sin(y./2.0-(sqrt(3.0).*x)./2.0).*sin(y).*2.0e+2-sin(y./2.0+(sqrt(3.0).*x)./2.0).*sin(y).*2.0e+2+cos(y./2.0-(sqrt(3.0).*x)./2.0).^2.*1.0e+2+cos(y./2.0+(sqrt(3.0).*x)./2.0).^2.*1.0e+2+cos(y).^2.*1.0e+2+sin(y./2.0-(sqrt(3.0).*x)./2.0).^2.*1.0e+2+sin(y./2.0+(sqrt(3.0).*x)./2.0).^2.*1.0e+2+sin(y.*(3.0./2.0)-(sqrt(3.0).*x)./2.0).^2.*2.25e+2+sin(y.*(3.0./2.0)+(sqrt(3.0).*x)./2.0).^2.*2.25e+2+sin(y).^2.*1.0e+2-z.^2.*1.0e+2+sin(sqrt(3.0).*x).*sin(y.*(3.0./2.0)-(sqrt(3.0).*x)./2.0).*4.5e+2-sin(sqrt(3.0).*x).*sin(y.*(3.0./2.0)+(sqrt(3.0).*x)./2.0).*4.5e+2+-8.99e+2+6.0e+1i)],[2,2])

[Xs,Ys,Zs]=meshgrid(xs,ys,zs);
% Compute result from the precomputed function
tic
answer = arrayfun(@(x,y,z)out(x,y,z),Xs(:),Ys(:),Zs(:),'UniformOutput',0);
s=sum(cat(3,answer{:}),3)

s = 
   1.0e+05 *

  -0.4528 - 1.4036i   0.0005 - 0.0010i
   0.0007 - 0.0010i  -0.3424 - 1.6394i

toc

Elapsed time is 2.888949 seconds.

function out = fun(x,y,z)
J = 1;
S = 1;
D = 0.75;
Sigma0_R = -1i*0.1*eye(2);
h0 = eye(2)*(3*J*S);
hx = [0,1; 1,0]*(-J*S*(cos(y/2 - (3^(1/2)*x)/2) + cos(y/2 + (3^(1/2)*x)/2) + cos(y)));
hy = [0, -1i; 1i, 0]*(-J*S*(sin(y/2 - (3^(1/2)*x)/2) + sin(y/2 + (3^(1/2)*x)/2) - sin(y)));
hz = [1, 0; 0, -1]*(-2*D*S*(sin(3^(1/2)*x) + sin((3*y)/2 - (3^(1/2)*x)/2) - sin((3*y)/2 + (3^(1/2)*x)/2)));
out = inv( z*eye(2) - (h0 + hx + hy + hz) - Sigma0_R );
end

Answer 3

[Matt J]

dx = 0.1;
dy = 0.15;
dz = 0.05;
xmin = -2;
xmax = 4;
ymin = -2*pi;
ymax = 2*pi;
zmin = -1;
zmax = 2;
xs = xmin:dx:xmax;
ys = ymin:dy:ymax;
zs = zmin:dz:zmax;

tic;
answer = sum(fun(xs,ys,zs),3);
toc

Elapsed time is 0.049273 seconds.

function out = fun(x,y,z)
J = 1;
S = 1;
D = 0.75;
fnc=@(q) reshape(q,1,1,[]);
fn=@(q) shiftdim(q,-3);

x=x(:); y=y(:).'; z=fnc(z);


Sigma0_R = -1i*0.1*eye(2);


 h0 = eye(2)*(3*J*S);
 hx = [0,1; 1,0].*fn((-J.*S.*(cos(y./2 - (3.^(1./2).*x)./2) + cos(y./2 + (3.^(1./2).*x)./2) + cos(y))));
 hy = [0, -1i; 1i, 0].*fn((-J.*S.*(sin(y./2 - (3.^(1./2).*x)./2) + sin(y./2 + (3.^(1./2).*x)./2) - sin(y))));
 hz = [1, 0; 0, -1].*fn((-2.*D.*S.*(sin(3.^(1./2).*x) + sin((3.*y)./2 - (3.^(1./2).*x)./2) - sin((3.*y)./2 + (3.^(1./2).*x)./2))));

 tmp=fn(z);

 out = inv2( reshape( tmp.*eye(2) - (h0 + hx + hy + hz) - Sigma0_R ,2,2,[]));

end


function X=inv2(A) 
%Inverts 2x2xN matrix stack A. 
%Based on FEX contribution
%https://www.mathworks.com/matlabcentral/fileexchange/27762-small-size-linear-solver?s_tid=srchtitle
%by Bruno Luong

    A = reshape(A, 4, []).';
    X = [A(:,4) -A(:,2) -A(:,3) A(:,1)];
         
    % Determinant
    D = A(:,1).*A(:,4) - A(:,2).*A(:,3);
    X =  X./D;
    X = reshape(X.', 2, 2, []);

end