我有一个3D网格，X，Y，Z。我想创建一个新的3D数组，它是X，Y和amp；Z.该函数包括位于不同点的几个3D高斯的和。目前，我有一个for循环，它在我有高斯的不同点上运行，并且我有一组中心位置r0(n高斯，1:3)
[X,Y,Z]=meshgrid(-10:.1:10);
Psi=0*X;
for index = 1:nGauss
Psi = Psi + Gauss3D(X,Y,Z,[r0(index,1),r0(index,2),r0(index,3)]);
end
其中我的3D高斯函数是
function output=Gauss3D(X,Y,Z,r0)
output=exp(-(X-r0(1)).^2 + (Y-r0(2)).^2 + (Z-r0(3)).^2);
end
我很高兴重新设计这个函数，这是我代码中最慢的部分，必须多次执行，但我不知道如何将其矢量化，使其运行得更快。如有任何建议，将不胜感激
*****注意，原始函数中有一个平方根，并进行了修改，使其成为实际的高斯***
注意我修改了你的代码，创建了一个高斯，它是：
output=exp(-sqrt((X-r0(1)).^2 + (Y-r0(2)).^2 + (Z-r0(3)).^2));
这不构成高斯。我把它改成：
output = exp(-((X-r0(1)).^2 + (Y-r0(2)).^2 + (Z-r0(3)).^2));
(注意编号sqrt)。这是一个sigma=sqrt(1/2)的高斯。
如果这不是你想要的，那么这个答案可能对你不是很有用，因为你的函数没有高斯函数那么快变为0，因此更难截断，而且它是不可分离的。
矢量化这个代码是毫无意义的，正如其他答案所证明的那样。MATLAB的JIT完全能够以最快的速度运行它。但你可以通过注意高斯很快变为零，并且是可分离的来显著减少计算量：
您在这里进行的大多数
exp评估产生的数字都很小。你不需要计算这些，只需填写0即可。exp(-x.^2-y.^2)与exp(-x.^2).*exp(-y.^2)相同，计算成本要低得多。
让我们测试一下这两件事。这是测试代码：
function gaussian_test
N = 100;
r0 = rand(N,3)*20 - 10;
% Original
tic
[X,Y,Z] = meshgrid(-10:.1:10);
Psi1 = zeros(size(X));
for index = 1:N
Psi1 = Psi1 + Gauss3D(X,Y,Z,r0(index,:));
end
t = toc;
fprintf('original, time = %fn',t)
% Fast, large truncation
tic
[X,Y,Z] = deal(-10:.1:10);
Psi2 = zeros(numel(X),numel(Y),numel(Z));
for index = 1:N
Psi2 = Gauss3D_fast(Psi2,X,Y,Z,r0(index,:),5);
end
t = toc;
fprintf('tuncation = 5, time = %fn',t)
fprintf('mean abs error = %fn',mean(reshape(abs(Psi2-Psi1),[],1)))
fprintf('mean square error = %fn',mean(reshape((Psi2-Psi1).^2,[],1)))
fprintf('max abs error = %fn',max(reshape(abs(Psi2-Psi1),[],1)))
% Fast, smaller truncation
tic
[X,Y,Z] = deal(-10:.1:10);
Psi3 = zeros(numel(X),numel(Y),numel(Z));
for index = 1:N
Psi3 = Gauss3D_fast(Psi3,X,Y,Z,r0(index,:),3);
end
t = toc;
fprintf('tuncation = 3, time = %fn',t)
fprintf('mean abs error = %fn',mean(reshape(abs(Psi3-Psi1),[],1)))
fprintf('mean square error = %fn',mean(reshape((Psi3-Psi1).^2,[],1)))
fprintf('max abs error = %fn',max(reshape(abs(Psi3-Psi1),[],1)))
% DIPimage, same smaller truncation
tic
Psi4 = newim(201,201,201);
coords = (r0+10) * 10;
Psi4 = gaussianblob(Psi4,coords,10*sqrt(1/2),(pi*100).^(3/2));
t = toc;
fprintf('DIPimage, time = %fn',t)
fprintf('mean abs error = %fn',mean(reshape(abs(Psi4-Psi1),[],1)))
fprintf('mean square error = %fn',mean(reshape((Psi4-Psi1).^2,[],1)))
fprintf('max abs error = %fn',max(reshape(abs(Psi4-Psi1),[],1)))
end % of function gaussian_test
function output = Gauss3D(X,Y,Z,r0)
output = exp(-((X-r0(1)).^2 + (Y-r0(2)).^2 + (Z-r0(3)).^2));
end
function Psi = Gauss3D_fast(Psi,X,Y,Z,r0,trunc)
% sigma = sqrt(1/2)
x = X-r0(1);
y = Y-r0(2);
z = Z-r0(3);
mx = abs(x) < trunc*sqrt(1/2);
my = abs(y) < trunc*sqrt(1/2);
mz = abs(z) < trunc*sqrt(1/2);
Psi(my,mx,mz) = Psi(my,mx,mz) + exp(-x(mx).^2) .* reshape(exp(-y(my).^2),[],1) .* reshape(exp(-z(mz).^2),1,1,[]);
% Note! the line above uses implicit singleton expansion. For older MATLABs use bsxfun
end
这是我机器上的输出，为了可读性而重新排序(我仍在使用MATLAB R2017a)：
| time(s) | mean abs | mean sq. | max abs
--------------+----------+----------+----------+----------
original | 5.035762 | | |
tuncation = 5 | 0.169807 | 0.000000 | 0.000000 | 0.000005
tuncation = 3 | 0.054737 | 0.000452 | 0.000002 | 0.024378
DIPimage | 0.044099 | 0.000452 | 0.000002 | 0.024378
正如你所看到的，使用高斯的这两个特性，我们可以将时间从5.0秒减少到0.17秒，速度提高了30倍，几乎没有明显的差异(在5*sigma处截断)。通过允许一个小误差，可以进一步获得3倍的加速。截断值越小，速度就越快，但误差就越大。
我添加了最后一个方法，DIPimage中的gaussianblob函数(我是一名作者)，只是为了在需要从代码中挤出额外时间的情况下显示该选项。该函数是用C++实现的。我使用的这个版本你需要自己编译。我们目前的官方版本仍然在M文件代码中实现了这个功能，而且速度没有那么快。
如果坐标的小数部分总是相同(相对于像素网格)，则有进一步的改进机会。在这种情况下，可以绘制一次高斯，然后将其移动到每个质心。
另一种替代方案涉及在稍大的尺度上计算一次高斯，并对其进行插值，以生成生成输出所需的每个1D高斯。我没有实现这一点，我不知道它是否会更快，或者时差是否会很大。在过去，exp是昂贵的，我不确定现在是否仍然如此。
所以，我正在构建我上面的答案@Durkee。我喜欢这类问题，所以我想了一点如何使每个展开式都是隐式的，下面是一行函数。使用这个功能，我把通话时间缩短了.11秒，这完全可以忽略不计。看起来你的相当不错。我的唯一优势可能是代码如何在更精细的网格上缩放。
xLin = [-10:.1:10]';
tic
psi2 = sum(exp(-sqrt((permute(xLin-r0(:,1)',[3 1 4 2])).^2 ...
+ (permute(xLin-r0(:,2)',[1 3 4 2])).^2 ...
+ (permute(xLin-r0(:,3)',[3 4 1 2])).^2)),4);
toc
我电脑上的相对运行时间是(所有东西都保持不变)：
Original - 1.234085
Other - 2.445375
Mine - 1.120701
所以这是一个有点不寻常的问题，在我的计算机上，未矢量化的代码实际上比矢量化的编码工作得更好，这是我的脚本
clear
[X,Y,Z]=meshgrid(-10:.1:10);
Psi=0*X;
nGauss = 20; %Sample nGauss as you didn't specify
r0 = rand(nGauss,3); % Just make this up as it doesn't really matter in this case
% Your original code
tic
for index = 1:nGauss
Psi = Psi + Gauss3D(X,Y,Z,[r0(index,1),r0(index,2),r0(index,3)]);
end
toc
% Vectorize these functions so we can use implicit broadcasting
X1 = X(:);
Y1 = Y(:);
Z1 = Z(:);
tic
val = [X1 Y1 Z1];
% Change the dimensions so that r0 operates on the right elements
r0_temp = permute(r0,[3 2 1]);
% Perform the gaussian combination
out = sum(exp(-sqrt(sum((val-r0_temp).^2,2))),3);
toc
% Check to make sure both functions match
sum(abs(vec(Psi)-vec(out)))
function output=Gauss3D(X,Y,Z,r0)
output=exp(-sqrt((X-r0(1)).^2 + (Y-r0(2)).^2 + (Z-r0(3)).^2));
end
function out = vec(in)
out = in(:);
end
正如您所看到的，这可能是尽可能矢量化的。整个功能是使用广播和矢量化操作来完成的，这通常会将性能提高一百倍。然而，在这种情况下，这不是我们所看到的
Elapsed time is 1.876460 seconds.
Elapsed time is 2.909152 seconds.
这实际上表明未经系数化的版本速度更快。
这可能有几个原因，而我绝非专家。
- MATLAB现在使用JIT编译器，这意味着for循环不再低效
- 您的代码已经合理地向量化了，您一次操作800万个元素
- 除非nGauss是1000左右，否则你不会循环那么多，在这一点上，矢量化意味着你将耗尽内存
- 我可能达到了内存阈值，因为我使用了太多内存，这使我的代码效率低下。我注意到，当我降低网格上的分辨率时，矢量化版本效果更好
顺便说一句，我在GTX 1060 GPU上测试了单精度(单精度比大多数GPU上的双精度快10倍)
Elapsed time is 0.087405 seconds.
Elapsed time is 0.241456 seconds.
再一次，非系数化版本更快，很抱歉我不能帮你，但除非你降低网格上的公差，否则你的代码似乎和你想要的一样好。