performance - 提高效率Matlab

Question

我有一个matlab代码，效率很低，我需要运行几次。
代码基本上是一个大的parfor loop，我想这几乎是不可能绕过的。
代码首先加载几个参数和四维矩阵，然后需要进行一些插值所有操作都需要执行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循环的方法。

Answer 1

为了计时和测试代码，我删除了parfor循环，并用for循环替换它，然后使用MATLAB profiler我用nsims = 500来做测试。
使用profiler，我发现了代码中的两个关键瓶颈第一个是嵌套最多的for循环（n-loop）中的find()函数第二个是对interpn()函数的三次调用这4行使用了+88%的计算时间
在这种情况下，由于函数调用的开销（特别是考虑到它在嵌套循环中接收的调用数）以及内置的错误检查和管理，find函数的速度比期望的要慢将find函数替换为硬编码二进制搜索（如下所示）可以极大地提高性能，而这只是在n循环中替换find对find使用nsims = 500可获得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而不是 interpn相结合，将总运行时间减少到3.8秒，比初始运行时间提高了近8倍。
还有三点需要注意：
1）我建议遵循前面答案中概述的建议，特别是尽可能远离结构，并将您能移动的东西移出循环（特别是ndgrid，因为这肯定只需要运行一次）。
2）我注意到，当使用 nsims=5000时，此脚本期间使用的总内存接近2.5gig如果这接近系统的总可用内存，则可能导致显著的速度减慢在这种情况下，我建议执行较小的计算批，保存结果，然后继续进行进一步的计算。
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

并将对 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的调用，取而代之的是：

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