向量化代码 - 如何减少 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

代码执行随机过程的 N 个不同模拟（由 10^3 个事件组成，发生在离散时刻（T 向量），具体取决于具体模拟。 现在（第二部分）我想知道，作为时间的函数，有多少模拟处于特定状态（X 假设值在 1 到 10 之间）。我的想法是：创建一个网格向量，其中包含在任何模拟中发生的所有时刻。然后，遍历模拟，遍历发生某些事情的时间步长，并递增与该特定时间片相对应的所有计数器 indeces。

但是，第二部分非常繁重（我的意思是在标准四核 CPU 上处理的天数）。它不应该。 是否有任何想法（也许是关于以更有效的方式比较向量）来减少 CPU 时间？

这是一个独立的“第二部分”

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

---

编辑：

从下面的 cmets 中，我了解到您记录了所有状态变化，以及它们发生的时间步长。然后，每次模拟更改状态时，您都希望从所有模拟中收集所有状态并计算每种类型有多少状态。

这里的主要问题是你的时间是连续的，所以基本上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 小时...

当然，以上所有时间都是使用相同算法的粗略估计。如果您确实选择舍弃时间，可能会有更好的方法来解决这个问题。