向量化代码 - 如何减少MATLAB计算时间

Question

我有这样一段代码

N=10^4; 
for i = 1:N 
    [E,X,T] = fffun(); % Stochastic simulation. Returns every time three different vectors (whose length is 10^3). 
    X_(i,:)=X; 
    T_(i,:)=T; 
    GRID=[GRID T]; 
end 
GRID=unique(GRID); 
% Second part 
for i=1:N 
for j=1:(kmax) 
    f=find(GRID==T_(i,j) | GRID==T_(i,j+1)); 
    s=f(1); 
    e=f(2)-1; 

counter(X_(i,j), s:e)=counter(X_(i,j), s:e)+1; 
end 
end 

代码执行一个随机过程（它由10^3点的事件，在离散的时刻发生的N个不同的模拟（T （第二部分）作为时间的函数，我想知道在一个特定的状态下有多少个模拟（X取值在1到10之间）。我有这样的想法：创建一个具有所有时刻的网格矢量，然后，在模拟上循环，循环发生某些事情的时间步并递增与此p相对应的所有计数器余数时间的关键片段。

然而，这第二部分是非常重的（我的意思是在一个标准的四核CPU上处理几天）。它不应该。 是否有任何想法（可能是以更高效的方式比较向量）来缩短CPU时间？

这是一个独立 'SECOND_PART'

N=5000; 
counter=zeros(11,length(GRID)); 

for i=1:N 
    disp(['Counting sim #' num2str(i)]); 
    for j=1:(kmax) 
     f=find(GRID==T_(i,j) | GRID==T_(i,j+1),2); 
     s=f(1); 
     e=f(2)-1; 

     counter(X_(i,j), s:e)=counter(X_(i,j), s:e)+1; 

    end 
end 

counter=counter/N; 
stop=find(GRID==Tmin); 
stop=stop-1; 
plot(counter(:,(stop-500):stop)') 

具有相关联的伪数据（filedropper.com/data_38）。在真实的情况下，矩阵有2x行和10x列。

Answer 1

这里是我的理解：

T_是从N个模拟时间步的矩阵。
X_是在那些模拟中的模拟状态矩阵T_。

所以，如果你这样做：

[ut,~,ic]= unique(T_(:)); 

你得到ic这是指数在T_所有独特元素的矢量。那么你可以写：

counter = accumarray([ic X_(:)],1); 

并得到counter没有。的行作为你独特的时间步调，没有。作为X_（它们都是，必须是整数）中的唯一状态列。现在你可以说，对于每个时间步ut(k)模拟处于状态m的时间数是counter(k,m)。

在您的数据中，m和k的值大于1的唯一组合是(1,1)。

---

编辑：

从下面的评论，我知道你记录所有的状态变化，以及何时发生的时间步骤。然后每次模拟改变一个状态时，你想从所有模拟中收集所有状态，并计算每种类型有多少个状态。

这里的主要问题是，你的时间是连续的，所以基本上在T_每个元素都是独特，你有超过一百万的时间步长遍历。完全矢量化这样的过程将需要大约80GB的内存，这可能会卡住你的电脑。

所以我寻找矢量化和循环的时间步骤的组合。首先，我们要找出所有独特的间隔和预分配counter：

ut = unique(T_(:)); 
stt = 11; % no. of states 
counter = zeros(stt,numel(ut));r = 1:size(T_,1); 
r = 1:size(T_,1); % we will need that also later 

然后我们循环遍历所有ut元素，每一次找相关的时间步长T_在所有模拟的量化方法。最后，我们使用histcounts计算所有的状态：

for k = 1:numel(ut) 
    temp = T_<=ut(k); % mark all time steps before ut(k) 
    s = cumsum(temp,2); % count the columns 
    col_ind = s(:,end); % fins the column index for each simulation 
    % convert the coulmns to linear indices: 
    linind = sub2ind(size(T_),r,col_ind.'); 
    % count the states: 
    counter(:,k) = histcounts(X_(linind),1:stt+1); 
end 

这在我的电脑模拟1000大约需要4秒钟，因此将其添加到略多于一个小时的整个过程。不是很快......

您也可以尝试一个或两个以下调整的挤运行时间多一点点：

1. 正如你can read here，accumarray似乎在小数组更快的工作，然后histcouns。所以可能想切换到它。

2. 此外，直接计算线性索引是比sub2ind更快的方法，因此您可能需要尝试。

实施上述循环这些建议，我们得到：

R = size(T_,1); 
r = (1:R).'; 
for k = 1:K 
    temp = T_<=ut(k); % mark all time steps before ut(k) 
    s = cumsum(temp,2); % count the columns 
    col_ind = s(:,end); % fins the column index for each simulation 
    % convert the coulmns to linear indices: 
    linind = R*(col_ind-1)+r; 
    % count the states: 
    counter(:,k) = accumarray(X_(linind),1,[stt 1]); 
end 

在我的电脑切换到accumarray和或删除sub2ind获得略有改善，但它不是一致的（使用timeit，以测试在ut中有100个或1K个元素），所以你最好自己测试一下。但是，这仍然很长。你可能要考虑

---

有一件事是想你的离散时间步长，所以你必须遍历更独特的元素。在你的数据中，约8％的时间间隔小于1.如果你可以认为这个时间足够短，那么你可以围绕你的T_，只得到约12.5K的独特元素，一分钟循环。你可以在0.1的时间间隔（小于1％的时间间隔）内完成相同的操作，并获得122K个元素循环，大约需要8个小时...

当然，以上所有时间都是粗略估计使用相同的算法。如果你选择做选择循环时间可能有更好的方法来解决这个问题。