提高效率matlab
Improve efficiency matlab
我有一个非常低效的 matlab 代码，我需要 运行 多次。
代码基本上是一个很大的 parfor loop ，我想几乎不可能绕过它。
代码首先加载几个参数和 4-D 矩阵，然后需要进行一些插值。所有都需要完成 5000 次（因此 parfor 循环）。
代码如下。我尽可能地简化了代码，没有去掉关键成分。
load file
nsim = 5000
T = 12;
N = 1000;
cumQx = cumsum(Qx);
cumQz = cumsum(Qz);
cumQs = cumsum(Qs);
for k=1:nsim
st(k).ksim = kstar*ones(N, T);
st(k).Vsim = zeros(N,T);
st(k).Psim = zeros(N,T);
end
parfor k = 1:nsim
sysrand = rand(T, 1);
idiorand = rand(N, T);
sigmarand = rand(T,1);
xid = zeros(T, 1);
zid = zeros(N, T);
sid = zeros(T,1);
xid(1) = 8;
zid(:, 1) = 5;
sid(1) = 1;
% Initializing the simulation
simx = zeros(T,1);
zsim = ones(N,T)*zbar;
simsx = zeros(T,1);
% Construct 3-D grid using 'ndgrid'
[ks,zs] = ndgrid(kgrid,z);
for j = 2:T
sid(j) = find(cumQs(:, sid(j-1)) >= sigmarand(j), 1);
simsx(j-1) = sigmax(sid(j));
xid(j) = find(cumQx(:, xid(j-1)) >= sysrand(j), 1);
simx(j-1) = x(xid(j));
for n = 1:N
zid(n, j) = find(cumQz(:, zid(n, j-1)) >= idiorand(n, j), 1);
zsim(n,j-1) = z(zid(n, j));
end
st(k).ksim(:,j) = interpn(ks, zs , squeeze(kprime(:,xid(j),:,sid(j))), st(k).ksim(:,j-1),zsim(:,j-1),'linear'); % K
st(k).Vsim(:,j) = interpn(ks, zs , squeeze(V(:,xid(j),:,sid(j))), st(k).ksim(:,j-1),zsim(:,j-1),'linear'); % V
st(k).Psim(:,j) = interpn(ks, zs , squeeze(P(:,xid(j),:,sid(j))), st(k).ksim(:,j-1),zsim(:,j-1),'linear'); % P
end;
end
这里是 link 矩阵 运行 代码：http://www.filedropper.com/file_406
有没有更好的方法可以显着减少计算时间？我的猜测是不...
理想情况下，将有一种向量化 k = 1:nsim 循环的方法。
撇开更高级别的优化（代码似乎是可向量化的，但我自己不 运行 就不知道），你在循环中做了很多不必要的计算/内存重新分配。
你真的需要清零吗xid, zid, sid, simx, zsim , 和 simsx?如果它们有垃圾值或零并没有什么坏处，因为你每次只访问当前索引。更不用说，您可能会因为不清零而节省大量内存分配操作。
此外，没有必要重新计算 ndgrid。
将所有这些东西移到循环之外。
- 尽量不要使用结构，这也会使您的代码更容易向量化（如果可能的话）。
在代码的开头，替换：
for k=1:nsim
st(k).ksim = kstar*ones(N, T);
st(k).Vsim = zeros(N,T);
st(k).Psim = zeros(N,T);
end
与：
ksim = kstar*ones(N,T,nsim);
Vsim = zeros(N,T,nsim);
Psim = zeros(N,T,nsim);
然后在最后而不是调用 st，而是调用：
ksim(:,j,k) = interpn(ks,zs,squeeze(kprime(:,xid(j),:,sid(j))),ksim(:,j-1,k),zsim(:,j-1),'linear');
Vsim(:,j,k) = interpn(ks,zs,squeeze(V(:,xid(j),:,sid(j))), Vsim(:,j-1,k),zsim(:,j-1),'linear');
Psim(:,j,k) = interpn(ks,zs,squeeze(P(:,xid(j),:,sid(j))), Psim(:,j-1,k),zsim(:,j-1),'linear');
如果需要，您可以在循环结束后将所有三个变量放入一个结构中。
- 从循环中取出所有随机生成器。
循环前写入：
sysrand = rand(T,nsim);
idiorand = rand(N, T,nsim);
sigmarand = rand(T,nsim);
然后在循环内替换：
sigmarand(j) --> sigmarand(j,k)
sysrand(j) --> sysrand(j,k)
idiorand(n,j) --> idiorand(n,j,k)
- 向量化最内层的 for 循环：
替换：
for n = 1:N
zid(n,j) = find(cumQz(:, zid(n, j-1)) >= idiorand(n,j), 1);
zsim(n,j-1) = z(zid(n, j));
end
与：
cumQz_k = cumQz(:, zid(:, j-1)).';
pos = bsxfun(@ge,cumQz_k,idiorand(:,j));
zid(:,j) = (sum(~cumsum(pos'))+1).';
zsim(:,j-1) = z(zid(:, j));
这将为您带来高达 90% 的显着提升，整体提升 60%。所有其他条件都相同（但不使用 parfor），nsim = 200 for 循环的整个模拟花费了 16s，其中 10s 专用于此循环。矢量化后，仿真只用了 6s，对于替代矢量化方法，他们大约需要一秒钟。
使用上面的所有备注，包括使用中建议的griddedInterpolant，将在不到一分钟的时间内将您带到5000，所以最好不要使用parfor ，这可能会因为开销通信（即使用大量内存）而减慢速度。
这是最终代码（请注意，我将您的 j 替换为 t）：
nsim = 5000;
T = 12;
N = 1000;
cumQx = cumsum(Qx);
cumQz = cumsum(Qz);
cumQs = cumsum(Qs);
ksim = kstar*ones(N,T,nsim);
Vsim = zeros(N,T,nsim);
Psim = zeros(N,T,nsim);
sysrand = rand(T,nsim);
idiorand = rand(N,T,nsim);
sigmarand = rand(T,nsim);
% Construct 3-D grid using 'ndgrid'
[ks,zs] = ndgrid(kgrid,z);
for k = 1:nsim
% Initializing the simulation:
xid = [8; zeros(T-1, 1)];
zid = [ones(N,1)*5 zeros(N,T-1)];
alter = [ones(N,1)*5 zeros(N,T-1)];
sid = [1; zeros(T-1,1)];
simx = zeros(T,1);
zsim = ones(N,T)*zbar;
simsx = zeros(T,1);
for t = 2:T
sid(t) = find(cumQs(:, sid(t-1)) >= sigmarand(t,k), 1);
simsx(t-1) = sigmax(sid(t));
xid(t) = find(cumQx(:, xid(t-1)) >= sysrand(t,k), 1);
simx(t-1) = x(xid(t));
cumQz_k = cumQz(:, zid(:, t-1)).';
pos = bsxfun(@ge,cumQz_k,idiorand(:,t,k));
zid(:,t) = (sum(~cumsum(pos'))+1).';
zsim(:,t-1) = z(zid(:, t));
F = griddedInterpolant(ks,zs,squeeze(kprime(:,xid(t),:,sid(t))), 'linear','none');
ksim(:,t,k) = F(ksim(:,t-1,k),zsim(:,t-1));
F = griddedInterpolant(ks,zs,squeeze(V(:,xid(t),:,sid(t))), 'linear','none');
Vsim(:,t,k) = F(ksim(:,t-1,k),zsim(:,t-1));
F = griddedInterpolant(ks,zs,squeeze(P(:,xid(t),:,sid(t))), 'linear','none');
Psim(:,t,k) = F(ksim(:,t-1,k),zsim(:,t-1));
end
end
非常欢迎您在评论中询问代码中任何不清楚的部分。
为了计时和测试您的代码，我删除了 parfor 循环并将其替换为 for 循环，然后使用了 MATLAB 探查器。我使用 nsims = 500 进行测试。
使用探查器，我在您的代码中发现了两个关键瓶颈。第一个是嵌套最多的 for 循环 (n-loop) 中的 find() 函数。第二个是 interpn() 函数的三个调用。 这 4 行使用了 +88% 的计算时间。
由于函数调用的开销（特别是考虑到它在嵌套循环中接收的调用次数）以及 built-in 错误检查和管理。将 find 函数替换为 hard-coded 二进制搜索（如下）大大提高了性能，而这只是替换了 n-loop 中的 find。对 nsims = 500 使用 find 导致 运行 时间为 29.8 秒。使用二进制搜索，运行 时间是 12.1 seconds.Note：这只有效，因为您的数据已排序，此代码不能在每个实例中替换查找。 编辑： 在@EBH 的另一个答案中使用替代方法是一种更简洁的方法。
%perform binary search (same as your find function)
searchVal=idiorand(n, j);
il = 1;
iu = sizeCumQZ; % should be defined outside the loop as size(cumQz,1)
while (il+1<iu)
lw=floor((il+iu)/2); % split the upper index
if cumQz(lw,zid(n, j-1)) >= searchVal
iu=lw; % decrease lower_index_b (whose x value remains \geq to lower bound)
else
il=lw; % increase lower_index_a (whose x value remains less than lower bound)
end
end
if cumQz(il,zid(n, j-1))>=searchVal
zid(n,j) = il;
else
zid(n,j) = iu;
end
interpn 函数由于方法检查、错误管理和数据格式管理而显着变慢。标准 interpn 中使用的大约 100 行代码可以减少到每次调用 2 行，并且通过知道我们只需要一种类型的插值并且我们的数据符合特定格式来显着提高性能。为此，我们直接使用 griddedInterpolant 函数（见下文）。再次使用nsims = 500，使用interpn函数（仍然使用un-changed find代码）有一个运行时间为29.8秒。使用下面的改进方法，将 运行 时间减少到 20.4 秒。
改进后的方法替换了对 interp 的调用，此处显示
st(k).ksim(:,j) = interpn(ks, zs , squeeze(kprime(:,xid(j),:,sid(j))), st(k).ksim(:,j-1),zsim(:,j-1),'linear'); % K
st(k).Vsim(:,j) = interpn(ks, zs , squeeze(V(:,xid(j),:,sid(j))), st(k).ksim(:,j-1),zsim(:,j-1),'linear'); % V
st(k).Psim(:,j) = interpn(ks, zs , squeeze(P(:,xid(j),:,sid(j))), st(k).ksim(:,j-1),zsim(:,j-1),'linear'); % P
更多直接调用 griddedInterpolant，如下所示：
F = griddedInterpolant(ks,zs,squeeze(kprime(:,xid(j),:,sid(j))), 'linear','none');
st(k).ksim(:,j) = F(st(k).ksim(:,j-1),zsim(:,j-1));
F = griddedInterpolant(ks,zs,squeeze(V(:,xid(j),:,sid(j))), 'linear','none');
st(k).Vsim(:,j) = F(st(k).ksim(:,j-1),zsim(:,j-1));
F = griddedInterpolant(ks,zs,squeeze(P(:,xid(j),:,sid(j))), 'linear','none');
st(k).Psim(:,j) = F(st(k).ksim(:,j-1),zsim(:,j-1));
将二进制搜索而不是 find 与对 griddedInterpolant 的调用而不是 interpn 相结合将总 运行 时间减少到 3.8 秒，比最初的 运行 时间提高了近 8 倍。
另外三点注意事项：
1) 我建议遵循之前答案中概述的建议，特别是尽可能远离结构，并尽可能将任何东西移出循环（特别是 ndgrid，因为这肯定只需要 运行一次）。
2) 我注意到，当使用 nsims=5000 时，此脚本使用的总内存接近 2.5gigs。如果这接近系统的总可用内存，则可能会导致速度显着下降。在这种情况下，我建议执行较小的批次计算，保存结果，然后继续进行进一步的计算。
3) 最后，在我的测试中，使用 parfor 实际上比使用标准的 for 循环要慢。我相信这是因为在这种情况下，每个单独的循环都是一个相对较短的操作，与指定并行作业和管理集群工作人员相关的开销超过了并行处理的收益。如果您在大型计算机集群上，情况可能并非如此，但在我的单台（4 核）机器上，parfor 仅导致 slow-down。如果您愿意，我会建议您在完成上述建议的更改后，使用您的实际工作代码自己测试每个案例。
我所做的完整代码更改如下所示，这些更改不包括结构使用方面的任何更改或其他仍然推荐的小优化。
load file
tic;
nsim = 500
T = 12;
N = 1000;
searchVal=1;
il = 1;
iu = 1;
cumQx = cumsum(Qx);
cumQz = cumsum(Qz);
cumQs = cumsum(Qs);
sizeCumQZ = size(cumQz,1);
for k=1:nsim
st(k).ksim = kstar*ones(N, T);
st(k).Vsim = zeros(N,T);
st(k).Psim = zeros(N,T);
end
%was parfor
for k = 1:nsim
sysrand = rand(T, 1);
idiorand = rand(N, T);
sigmarand = rand(T,1);
xid = zeros(T, 1);
zid = zeros(N, T);
sid = zeros(T,1);
xid(1) = 8;
zid(:, 1) = 5;
sid(1) = 1;
% Initializing the simulation
simx = zeros(T,1);
zsim = ones(N,T)*zbar;
simsx = zeros(T,1);
% Construct 3-D grid using 'ndgrid'
[ks,zs] = ndgrid(kgrid,z);
for j = 2:T
sid(j) = find(cumQs(:, sid(j-1)) >= sigmarand(j), 1);
simsx(j-1) = sigmax(sid(j));
xid(j) = find(cumQx(:, xid(j-1)) >= sysrand(j), 1);
simx(j-1) = x(xid(j));
for n = 1:N
%perform binary search (same as your find function)
searchVal=idiorand(n, j);
il = 1;
iu = sizeCumQZ;
while (il+1<iu)
lw=floor((il+iu)/2); % split the upper index
if cumQz(lw,zid(n, j-1)) >= searchVal
iu=lw; % decrease lower_index_b (whose x value remains \geq to lower bound)
else
il=lw; % increase lower_index_a (whose x value remains less than lower bound)
end
end
if cumQz(il,zid(n, j-1))>=searchVal
zid(n,j) = il;
else
zid(n,j) = iu;
end
zsim(n,j-1) = z(zid(n,j));
end
F = griddedInterpolant(ks,zs,squeeze(kprime(:,xid(j),:,sid(j))), 'linear','none');
st(k).ksim(:,j) = F(st(k).ksim(:,j-1),zsim(:,j-1));
F = griddedInterpolant(ks,zs,squeeze(V(:,xid(j),:,sid(j))), 'linear','none');
st(k).Vsim(:,j) = F(st(k).ksim(:,j-1),zsim(:,j-1));
F = griddedInterpolant(ks,zs,squeeze(P(:,xid(j),:,sid(j))), 'linear','none');
st(k).Psim(:,j) = F(st(k).ksim(:,j-1),zsim(:,j-1));
end;
end
toc;
编辑： 稍微弄乱 griddedInterpolant，通过将三个插值组合成一个，我能够额外增加 15% 的速度K、V 和 P 插值的网格和值的串联。在代码的顶部，最好在循环之外完成，我将原来的网格创建替换为：
zRange = max(z(:))-min(z(:))+1;
zzzGrid = [z;z+1*zRange;z+2*zRange];% for K, V, and P
[ksBig,zsBig] = ndgrid(kgrid,zzzGrid);
nZ = numel(z); %used below
valGrid = zeros(size(ksBig)); %used below
并将对 griddedInterpolant 的 3 次调用替换为：
valGrid(:,1:nZ) = squeeze(kprime(:,xid(j),:,sid(j)));%K
valGrid(:,nZ+1:2*nZ) = squeeze(V(:,xid(j),:,sid(j)));%V
valGrid(:,2*nZ+1:3*nZ) = squeeze(P(:,xid(j),:,sid(j)));%P
F = griddedInterpolant(ksBig,zsBig,valGrid, 'linear','none');
st(k).ksim(:,j) = F(st(k).ksim(:,j-1),zsim(:,j-1));
st(k).Vsim(:,j) = F(st(k).ksim(:,j-1),zsim(:,j-1)+zRange);
st(k).Psim(:,j) = F(st(k).ksim(:,j-1),zsim(:,j-1)+2*zRange);
或者如果我们可以用更复杂的设置换取 19% 的增加而不是 15% 的增加，我们可以将 griddedInterpolant 完全移出 j 循环。在代码的开头，设置网格如下：
zRange = max(z(:))-min(z(:))+1;
zzzGrid = [z;z+1*zRange;z+2*zRange];
zzzRange = max(zzzGrid(:))-min(zzzGrid(:))+1;
zzzTGrid = [];
for j = 2:T
zzzTGrid(end+1:end+numel(zzzGrid)) = zzzGrid+(j-2)*zzzRange;
end
[ksBig,zsBig] = ndgrid(kgrid,zzzTGrid);
nZ = numel(z); %used below
valGrid = zeros(size(ksBig)); %used below
并替换之前的griddedInterpolant如下图：
for j = 2:T
%%%%%
%...
%Other code in the j loop
%...
%%%%%
valGrid(:,(1:nZ)+3*nZ*(j-2)) = squeeze(kprime(:,xid(j),:,sid(j)));%K
valGrid(:,(nZ+1:2*nZ)+3*nZ*(j-2)) = squeeze(V(:,xid(j),:,sid(j)));%V
valGrid(:,(2*nZ+1:3*nZ)+3*nZ*(j-2)) = squeeze(P(:,xid(j),:,sid(j)));%P
end;
F = griddedInterpolant(ksBig,zsBig,valGrid, 'linear','none');
for j = 2:T
st(k).ksim(:,j) = F(stTest(k).ksim(:,j-1),zsim(:,j-1)+3*zRange*(j-2));
st(k).Vsim(:,j) = F(stTest(k).ksim(:,j-1),zsim(:,j-1)+zRange+3*zRange*(j-2));
st(k).Psim(:,j) = F(stTest(k).ksim(:,j-1),zsim(:,j-1)+2*zRange+3*zRange*(j-2));
end
如果我们想要更挑剔（速度提高 2%），我们可以将所有对 squeeze 的调用替换为不进行错误检查的版本 mySqueeze：
function b = mySqueeze(a)
%Trimmed down version of squeeze, a built-in MATLAB function, has no error-managment or case optimization
siz = size(a);
siz(siz==1) = []; % Remove singleton dimensions.
siz = [siz ones(1,2-length(siz))]; % Make sure siz is at least 2-D
b = reshape(a,siz);
我有一个非常低效的 matlab 代码，我需要 运行 多次。
代码基本上是一个很大的 parfor loop ，我想几乎不可能绕过它。
代码首先加载几个参数和 4-D 矩阵，然后需要进行一些插值。所有都需要完成 5000 次（因此 parfor 循环）。
代码如下。我尽可能地简化了代码，没有去掉关键成分。
load file
nsim = 5000
T = 12;
N = 1000;
cumQx = cumsum(Qx);
cumQz = cumsum(Qz);
cumQs = cumsum(Qs);
for k=1:nsim
st(k).ksim = kstar*ones(N, T);
st(k).Vsim = zeros(N,T);
st(k).Psim = zeros(N,T);
end
parfor k = 1:nsim
sysrand = rand(T, 1);
idiorand = rand(N, T);
sigmarand = rand(T,1);
xid = zeros(T, 1);
zid = zeros(N, T);
sid = zeros(T,1);
xid(1) = 8;
zid(:, 1) = 5;
sid(1) = 1;
% Initializing the simulation
simx = zeros(T,1);
zsim = ones(N,T)*zbar;
simsx = zeros(T,1);
% Construct 3-D grid using 'ndgrid'
[ks,zs] = ndgrid(kgrid,z);
for j = 2:T
sid(j) = find(cumQs(:, sid(j-1)) >= sigmarand(j), 1);
simsx(j-1) = sigmax(sid(j));
xid(j) = find(cumQx(:, xid(j-1)) >= sysrand(j), 1);
simx(j-1) = x(xid(j));
for n = 1:N
zid(n, j) = find(cumQz(:, zid(n, j-1)) >= idiorand(n, j), 1);
zsim(n,j-1) = z(zid(n, j));
end
st(k).ksim(:,j) = interpn(ks, zs , squeeze(kprime(:,xid(j),:,sid(j))), st(k).ksim(:,j-1),zsim(:,j-1),'linear'); % K
st(k).Vsim(:,j) = interpn(ks, zs , squeeze(V(:,xid(j),:,sid(j))), st(k).ksim(:,j-1),zsim(:,j-1),'linear'); % V
st(k).Psim(:,j) = interpn(ks, zs , squeeze(P(:,xid(j),:,sid(j))), st(k).ksim(:,j-1),zsim(:,j-1),'linear'); % P
end;
end
这里是 link 矩阵 运行 代码：http://www.filedropper.com/file_406
有没有更好的方法可以显着减少计算时间？我的猜测是不... 理想情况下，将有一种向量化 k = 1:nsim 循环的方法。
撇开更高级别的优化（代码似乎是可向量化的，但我自己不 运行 就不知道），你在循环中做了很多不必要的计算/内存重新分配。
你真的需要清零吗xid, zid, sid, simx, zsim , 和 simsx?如果它们有垃圾值或零并没有什么坏处，因为你每次只访问当前索引。更不用说，您可能会因为不清零而节省大量内存分配操作。
此外，没有必要重新计算 ndgrid。
将所有这些东西移到循环之外。
- 尽量不要使用结构，这也会使您的代码更容易向量化（如果可能的话）。
在代码的开头，替换：
for k=1:nsim
st(k).ksim = kstar*ones(N, T);
st(k).Vsim = zeros(N,T);
st(k).Psim = zeros(N,T);
end
与：
ksim = kstar*ones(N,T,nsim);
Vsim = zeros(N,T,nsim);
Psim = zeros(N,T,nsim);
然后在最后而不是调用 st，而是调用：
ksim(:,j,k) = interpn(ks,zs,squeeze(kprime(:,xid(j),:,sid(j))),ksim(:,j-1,k),zsim(:,j-1),'linear');
Vsim(:,j,k) = interpn(ks,zs,squeeze(V(:,xid(j),:,sid(j))), Vsim(:,j-1,k),zsim(:,j-1),'linear');
Psim(:,j,k) = interpn(ks,zs,squeeze(P(:,xid(j),:,sid(j))), Psim(:,j-1,k),zsim(:,j-1),'linear');
如果需要，您可以在循环结束后将所有三个变量放入一个结构中。
- 从循环中取出所有随机生成器。
循环前写入：
sysrand = rand(T,nsim);
idiorand = rand(N, T,nsim);
sigmarand = rand(T,nsim);
然后在循环内替换：
sigmarand(j) --> sigmarand(j,k)
sysrand(j) --> sysrand(j,k)
idiorand(n,j) --> idiorand(n,j,k)
- 向量化最内层的 for 循环：
替换：
for n = 1:N
zid(n,j) = find(cumQz(:, zid(n, j-1)) >= idiorand(n,j), 1);
zsim(n,j-1) = z(zid(n, j));
end
与：
cumQz_k = cumQz(:, zid(:, j-1)).';
pos = bsxfun(@ge,cumQz_k,idiorand(:,j));
zid(:,j) = (sum(~cumsum(pos'))+1).';
zsim(:,j-1) = z(zid(:, j));
这将为您带来高达 90% 的显着提升，整体提升 60%。所有其他条件都相同（但不使用 parfor），nsim = 200 for 循环的整个模拟花费了 16s，其中 10s 专用于此循环。矢量化后，仿真只用了 6s，对于替代矢量化方法，他们大约需要一秒钟。
使用上面的所有备注，包括使用griddedInterpolant，将在不到一分钟的时间内将您带到5000，所以最好不要使用parfor ，这可能会因为开销通信（即使用大量内存）而减慢速度。
这是最终代码（请注意，我将您的 j 替换为 t）：
nsim = 5000;
T = 12;
N = 1000;
cumQx = cumsum(Qx);
cumQz = cumsum(Qz);
cumQs = cumsum(Qs);
ksim = kstar*ones(N,T,nsim);
Vsim = zeros(N,T,nsim);
Psim = zeros(N,T,nsim);
sysrand = rand(T,nsim);
idiorand = rand(N,T,nsim);
sigmarand = rand(T,nsim);
% Construct 3-D grid using 'ndgrid'
[ks,zs] = ndgrid(kgrid,z);
for k = 1:nsim
% Initializing the simulation:
xid = [8; zeros(T-1, 1)];
zid = [ones(N,1)*5 zeros(N,T-1)];
alter = [ones(N,1)*5 zeros(N,T-1)];
sid = [1; zeros(T-1,1)];
simx = zeros(T,1);
zsim = ones(N,T)*zbar;
simsx = zeros(T,1);
for t = 2:T
sid(t) = find(cumQs(:, sid(t-1)) >= sigmarand(t,k), 1);
simsx(t-1) = sigmax(sid(t));
xid(t) = find(cumQx(:, xid(t-1)) >= sysrand(t,k), 1);
simx(t-1) = x(xid(t));
cumQz_k = cumQz(:, zid(:, t-1)).';
pos = bsxfun(@ge,cumQz_k,idiorand(:,t,k));
zid(:,t) = (sum(~cumsum(pos'))+1).';
zsim(:,t-1) = z(zid(:, t));
F = griddedInterpolant(ks,zs,squeeze(kprime(:,xid(t),:,sid(t))), 'linear','none');
ksim(:,t,k) = F(ksim(:,t-1,k),zsim(:,t-1));
F = griddedInterpolant(ks,zs,squeeze(V(:,xid(t),:,sid(t))), 'linear','none');
Vsim(:,t,k) = F(ksim(:,t-1,k),zsim(:,t-1));
F = griddedInterpolant(ks,zs,squeeze(P(:,xid(t),:,sid(t))), 'linear','none');
Psim(:,t,k) = F(ksim(:,t-1,k),zsim(:,t-1));
end
end
非常欢迎您在评论中询问代码中任何不清楚的部分。
为了计时和测试您的代码，我删除了 parfor 循环并将其替换为 for 循环，然后使用了 MATLAB 探查器。我使用 nsims = 500 进行测试。
使用探查器，我在您的代码中发现了两个关键瓶颈。第一个是嵌套最多的 for 循环 (n-loop) 中的 find() 函数。第二个是 interpn() 函数的三个调用。 这 4 行使用了 +88% 的计算时间。
由于函数调用的开销（特别是考虑到它在嵌套循环中接收的调用次数）以及 built-in 错误检查和管理。将 find 函数替换为 hard-coded 二进制搜索（如下）大大提高了性能，而这只是替换了 n-loop 中的 find。对 nsims = 500 使用 find 导致 运行 时间为 29.8 秒。使用二进制搜索，运行 时间是 12.1 seconds.Note：这只有效，因为您的数据已排序，此代码不能在每个实例中替换查找。 编辑： 在@EBH 的另一个答案中使用替代方法是一种更简洁的方法。
%perform binary search (same as your find function)
searchVal=idiorand(n, j);
il = 1;
iu = sizeCumQZ; % should be defined outside the loop as size(cumQz,1)
while (il+1<iu)
lw=floor((il+iu)/2); % split the upper index
if cumQz(lw,zid(n, j-1)) >= searchVal
iu=lw; % decrease lower_index_b (whose x value remains \geq to lower bound)
else
il=lw; % increase lower_index_a (whose x value remains less than lower bound)
end
end
if cumQz(il,zid(n, j-1))>=searchVal
zid(n,j) = il;
else
zid(n,j) = iu;
end
interpn 函数由于方法检查、错误管理和数据格式管理而显着变慢。标准 interpn 中使用的大约 100 行代码可以减少到每次调用 2 行，并且通过知道我们只需要一种类型的插值并且我们的数据符合特定格式来显着提高性能。为此，我们直接使用 griddedInterpolant 函数（见下文）。再次使用nsims = 500，使用interpn函数（仍然使用un-changed find代码）有一个运行时间为29.8秒。使用下面的改进方法，将 运行 时间减少到 20.4 秒。
改进后的方法替换了对 interp 的调用，此处显示
st(k).ksim(:,j) = interpn(ks, zs , squeeze(kprime(:,xid(j),:,sid(j))), st(k).ksim(:,j-1),zsim(:,j-1),'linear'); % K
st(k).Vsim(:,j) = interpn(ks, zs , squeeze(V(:,xid(j),:,sid(j))), st(k).ksim(:,j-1),zsim(:,j-1),'linear'); % V
st(k).Psim(:,j) = interpn(ks, zs , squeeze(P(:,xid(j),:,sid(j))), st(k).ksim(:,j-1),zsim(:,j-1),'linear'); % P
更多直接调用 griddedInterpolant，如下所示：
F = griddedInterpolant(ks,zs,squeeze(kprime(:,xid(j),:,sid(j))), 'linear','none');
st(k).ksim(:,j) = F(st(k).ksim(:,j-1),zsim(:,j-1));
F = griddedInterpolant(ks,zs,squeeze(V(:,xid(j),:,sid(j))), 'linear','none');
st(k).Vsim(:,j) = F(st(k).ksim(:,j-1),zsim(:,j-1));
F = griddedInterpolant(ks,zs,squeeze(P(:,xid(j),:,sid(j))), 'linear','none');
st(k).Psim(:,j) = F(st(k).ksim(:,j-1),zsim(:,j-1));
将二进制搜索而不是 find 与对 griddedInterpolant 的调用而不是 interpn 相结合将总 运行 时间减少到 3.8 秒，比最初的 运行 时间提高了近 8 倍。
另外三点注意事项：
1) 我建议遵循之前答案中概述的建议，特别是尽可能远离结构，并尽可能将任何东西移出循环（特别是 ndgrid，因为这肯定只需要 运行一次）。
2) 我注意到，当使用 nsims=5000 时，此脚本使用的总内存接近 2.5gigs。如果这接近系统的总可用内存，则可能会导致速度显着下降。在这种情况下，我建议执行较小的批次计算，保存结果，然后继续进行进一步的计算。
3) 最后，在我的测试中，使用 parfor 实际上比使用标准的 for 循环要慢。我相信这是因为在这种情况下，每个单独的循环都是一个相对较短的操作，与指定并行作业和管理集群工作人员相关的开销超过了并行处理的收益。如果您在大型计算机集群上，情况可能并非如此，但在我的单台（4 核）机器上，parfor 仅导致 slow-down。如果您愿意，我会建议您在完成上述建议的更改后，使用您的实际工作代码自己测试每个案例。
我所做的完整代码更改如下所示，这些更改不包括结构使用方面的任何更改或其他仍然推荐的小优化。
load file
tic;
nsim = 500
T = 12;
N = 1000;
searchVal=1;
il = 1;
iu = 1;
cumQx = cumsum(Qx);
cumQz = cumsum(Qz);
cumQs = cumsum(Qs);
sizeCumQZ = size(cumQz,1);
for k=1:nsim
st(k).ksim = kstar*ones(N, T);
st(k).Vsim = zeros(N,T);
st(k).Psim = zeros(N,T);
end
%was parfor
for k = 1:nsim
sysrand = rand(T, 1);
idiorand = rand(N, T);
sigmarand = rand(T,1);
xid = zeros(T, 1);
zid = zeros(N, T);
sid = zeros(T,1);
xid(1) = 8;
zid(:, 1) = 5;
sid(1) = 1;
% Initializing the simulation
simx = zeros(T,1);
zsim = ones(N,T)*zbar;
simsx = zeros(T,1);
% Construct 3-D grid using 'ndgrid'
[ks,zs] = ndgrid(kgrid,z);
for j = 2:T
sid(j) = find(cumQs(:, sid(j-1)) >= sigmarand(j), 1);
simsx(j-1) = sigmax(sid(j));
xid(j) = find(cumQx(:, xid(j-1)) >= sysrand(j), 1);
simx(j-1) = x(xid(j));
for n = 1:N
%perform binary search (same as your find function)
searchVal=idiorand(n, j);
il = 1;
iu = sizeCumQZ;
while (il+1<iu)
lw=floor((il+iu)/2); % split the upper index
if cumQz(lw,zid(n, j-1)) >= searchVal
iu=lw; % decrease lower_index_b (whose x value remains \geq to lower bound)
else
il=lw; % increase lower_index_a (whose x value remains less than lower bound)
end
end
if cumQz(il,zid(n, j-1))>=searchVal
zid(n,j) = il;
else
zid(n,j) = iu;
end
zsim(n,j-1) = z(zid(n,j));
end
F = griddedInterpolant(ks,zs,squeeze(kprime(:,xid(j),:,sid(j))), 'linear','none');
st(k).ksim(:,j) = F(st(k).ksim(:,j-1),zsim(:,j-1));
F = griddedInterpolant(ks,zs,squeeze(V(:,xid(j),:,sid(j))), 'linear','none');
st(k).Vsim(:,j) = F(st(k).ksim(:,j-1),zsim(:,j-1));
F = griddedInterpolant(ks,zs,squeeze(P(:,xid(j),:,sid(j))), 'linear','none');
st(k).Psim(:,j) = F(st(k).ksim(:,j-1),zsim(:,j-1));
end;
end
toc;
编辑： 稍微弄乱 griddedInterpolant，通过将三个插值组合成一个，我能够额外增加 15% 的速度K、V 和 P 插值的网格和值的串联。在代码的顶部，最好在循环之外完成，我将原来的网格创建替换为：
zRange = max(z(:))-min(z(:))+1;
zzzGrid = [z;z+1*zRange;z+2*zRange];% for K, V, and P
[ksBig,zsBig] = ndgrid(kgrid,zzzGrid);
nZ = numel(z); %used below
valGrid = zeros(size(ksBig)); %used below
并将对 griddedInterpolant 的 3 次调用替换为：
valGrid(:,1:nZ) = squeeze(kprime(:,xid(j),:,sid(j)));%K
valGrid(:,nZ+1:2*nZ) = squeeze(V(:,xid(j),:,sid(j)));%V
valGrid(:,2*nZ+1:3*nZ) = squeeze(P(:,xid(j),:,sid(j)));%P
F = griddedInterpolant(ksBig,zsBig,valGrid, 'linear','none');
st(k).ksim(:,j) = F(st(k).ksim(:,j-1),zsim(:,j-1));
st(k).Vsim(:,j) = F(st(k).ksim(:,j-1),zsim(:,j-1)+zRange);
st(k).Psim(:,j) = F(st(k).ksim(:,j-1),zsim(:,j-1)+2*zRange);
或者如果我们可以用更复杂的设置换取 19% 的增加而不是 15% 的增加，我们可以将 griddedInterpolant 完全移出 j 循环。在代码的开头，设置网格如下：
zRange = max(z(:))-min(z(:))+1;
zzzGrid = [z;z+1*zRange;z+2*zRange];
zzzRange = max(zzzGrid(:))-min(zzzGrid(:))+1;
zzzTGrid = [];
for j = 2:T
zzzTGrid(end+1:end+numel(zzzGrid)) = zzzGrid+(j-2)*zzzRange;
end
[ksBig,zsBig] = ndgrid(kgrid,zzzTGrid);
nZ = numel(z); %used below
valGrid = zeros(size(ksBig)); %used below
并替换之前的griddedInterpolant如下图：
for j = 2:T
%%%%%
%...
%Other code in the j loop
%...
%%%%%
valGrid(:,(1:nZ)+3*nZ*(j-2)) = squeeze(kprime(:,xid(j),:,sid(j)));%K
valGrid(:,(nZ+1:2*nZ)+3*nZ*(j-2)) = squeeze(V(:,xid(j),:,sid(j)));%V
valGrid(:,(2*nZ+1:3*nZ)+3*nZ*(j-2)) = squeeze(P(:,xid(j),:,sid(j)));%P
end;
F = griddedInterpolant(ksBig,zsBig,valGrid, 'linear','none');
for j = 2:T
st(k).ksim(:,j) = F(stTest(k).ksim(:,j-1),zsim(:,j-1)+3*zRange*(j-2));
st(k).Vsim(:,j) = F(stTest(k).ksim(:,j-1),zsim(:,j-1)+zRange+3*zRange*(j-2));
st(k).Psim(:,j) = F(stTest(k).ksim(:,j-1),zsim(:,j-1)+2*zRange+3*zRange*(j-2));
end
如果我们想要更挑剔（速度提高 2%），我们可以将所有对 squeeze 的调用替换为不进行错误检查的版本 mySqueeze：
function b = mySqueeze(a)
%Trimmed down version of squeeze, a built-in MATLAB function, has no error-managment or case optimization
siz = size(a);
siz(siz==1) = []; % Remove singleton dimensions.
siz = [siz ones(1,2-length(siz))]; % Make sure siz is at least 2-D
b = reshape(a,siz);