logo标签列表

提高MATLAB的效率

performancematlabparfor
3
我有一个Matlab代码，它非常低效，我需要多次运行它。 该代码基本上是一个大的parfor循环，我认为几乎不可能绕过它。
该代码首先加载多个参数和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

这是运行代码所需矩阵的链接：http://www.filedropper.com/file_406。
有没有更好的方法可以显著减少计算时间？我猜没有......理想情况下，应该有一种向量化k = 1:nsim循环的方法。
2
你可能想要在这里提问：http://codereview.stackexchange.com/ - EBH
3个回答
6
为了测试你的代码，我去掉了parfor循环并用for循环代替，然后使用MATLAB分析器。 我在测试中使用了nsims = 500。
通过使用分析器，我发现你的代码中有两个关键瓶颈。 第一个是最嵌套的for循环（n-loop）中的find()函数。 第二个是三次调用interpn()函数。 这4行代码使用了计算时间的88％。 enter image description here 在这种情况下，由于函数调用的开销（尤其是考虑到嵌套循环中它接收的调用数量），以及内置的错误检查和管理，find函数比理想情况下要慢。将find函数替换为硬编码的二分搜索（如下所示）可大大提高性能，这仅替换了n循环中的find。对于nsims = 500，使用find的运行时间为29.8秒。使用二分搜索，运行时间为12.1秒。注意：这仅适用于数据已排序的情况，此代码无法替换每个实例中的find。编辑：使用@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函数（仍然使用未更改的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的调用，可以将总运行时间缩短到3.8秒，几乎比初始运行时间快了8倍。

还有三个需要注意的点：

1）我建议遵循先前答案中概述的建议，尽可能避免使用结构，并将您可以移出循环的任何内容移出（特别是ndgrid，因为这肯定只需要运行一次）。

2）我注意到，在使用nsims=5000时，此脚本使用的总内存接近2.5G。如果这接近您系统的总可用内存，则可能会导致显着的减速。在这种情况下，我建议执行较小的计算批次，保存结果，然后进行进一步的计算。

3)最后，在我的测试中，使用parfor比使用标准的for循环更慢。我认为这是因为在这种情况下，每个单独的循环都是一个相对较短的操作，而指定并行作业和管理集群工作人员所需的开销超过了并行处理的收益。如果您在一台大型计算机集群上，情况可能不同，但在我的单个（4核）机器上，parfor只会导致减速。如果您愿意，我建议在进行上述推荐更改后，使用您的实际工作代码对每种情况进行测试。

我所做的全部代码更改如下所示，其中不包括使用结构体或其他小优化所建议的任何更改。

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，我将三个插值器合并为一个，并通过连接K、V和P插值器的网格和值，成功地获得了额外的15%速度提升。在代码顶部（最好在循环外），我用以下代码替换了原始网格创建：

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

将3个对griddedInterpolant的调用替换为：

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);
你用interpn提出了一个很好的观点。看看我的答案，还有另一种更快的方法来替换find，而不需要假设数据已排序。+1 - EBH
@EBH 这是一个创造性的查找操作向量化，再加上1分给你！ - Tar
4
  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

使用：

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');

如果需要，在循环结束后，您可以将所有三个变量放入一个结构体中。

  1. 从循环中移除所有随机生成器。

在循环之前写：

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)
  1. 将最内层的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，其中10秒用于该循环。矢量化后，模拟仅花费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，而parfor在没有向量化的情况下可以显著提高速度。 - phdstudent
1
@volcompt 请查看我编辑后的答案，这可能会更进一步地向量化可能的内容。 - EBH
@volcompt 1) 对于计时，最好使用 timeit。2) 将模拟移动到另一个函数时，您还要将输出从结构更改为矩阵，这应该会带来一些改善。3) 我认为您应该先尝试向量化模拟函数，并将 nsim 维度留给 parfor。 - EBH
1
太好了！非常感谢！我会再把赏金留几天，以防其他人想要尝试，然后履行承诺并授予两个赏金。 - phdstudent
1
当然，这是我之前编写的one_sim函数遗留下的...我会编辑它。 - EBH
显示剩余8条评论
0

除了更高级别的优化（代码似乎可向量化，但我不知道没有运行它），您在循环内部进行了许多不必要的计算/内存重新分配。

您真的需要将xid、zid、sid、simx、zsim和simsx清零吗？如果它们具有垃圾值或零，因为您每次只访问当前索引，所以这并不会对其造成太大影响。更不用说，您可能会通过不清零来节省大量内存分配操作。

此外，没有必要重新计算ndgrid。

将所有这些东西移到循环外面。

转念一想，由于访问的并行性，零可能是必要的...但 ndgrid 绝对可以放在外面。 - Tasos Papastylianou
感谢您的评论。我肯定可以将ndgrid移到外面。但是这并没有显著提高效率。我试图理解如何使它可向量化，但是无法想出... - phdstudent
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接