向量化代码 - 如何减少 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 之间）。我的想法是：创建一个网格向量，其中包含在任何模拟中发生某些事情的所有时刻。然后，遍历模拟，遍历发生某事的时间步长，并递增与该特定时间片对应的所有计数器索引。

但是第二部分非常繁重（我的意思是在标准四核处理器上处理数天 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 秒，因此整个过程增加了一个多小时。不是很快...

您也可以尝试下面的一两个调整来压缩运行时间：

正如您一样，accumarray 似乎在小数组中比 histcouns 工作得更快。所以可能要切换到它。
另外，直接计算线性指数比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 测试 [=30 中的 100 或 1K 个元素=]), 所以你最好自己测试一下。然而，这仍然很长。

您可能要考虑的一件事是尝试离散化您的时间步长，这样您可以循环的独特元素就会少得多。在您的数据中，大约 8% 的时间间隔小于 1。如果您可以假设这足够短以被视为一个时间步长，那么您可以四舍五入 T_ 并仅获得 ~12.5K 独特元素，这需要大约一分钟的时间来循环。您可以对 0.1 个间隔（小于时间间隔的 1%）执行相同的操作，并让 122K 个元素循环，大约需要 8 个小时...

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

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

Vectorizing code - How to reduce MATLAB computational time

cpu

performance

time

matlab