The fastest method to perform sum of function

9 次查看（过去 30 天）
Luqman Saleem
Luqman Saleem 2024-2-2
评论： Matt J 2024-2-3
Ran in:
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
  5 个评论
显示 3更早的评论
Luqman Saleem
Luqman Saleem 2024-2-2
@Walter Roberson yes yes you are right. sorry, English is not my first language.
Matt J
Matt J 2024-2-3
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?

请先登录，再进行评论。

请先登录，再回答此问题。

回答（3 个）

Walter Roberson
Walter Roberson 2024-2-2
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.
  5 个评论
显示 3更早的评论
Aaron
Aaron 2024-2-2
In that case, you'd do sum(A, [3,4,5]) to just sum over the dimensions you're storing the dimensional dependence in.

请先登录，再进行评论。

Torsten
Torsten 2024-2-2
编辑：Torsten 2024-2-2
Ran in:
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
Matt J
Matt J 2024-2-2
编辑：Matt J 2024-2-3
Ran in:
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

类别

在 Help Center 和 File Exchange 中查找有关 Sparse Matrices 的更多信息

标签

产品


版本

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by