时间序列聚合效率

Question

我通常需要使用给定的聚合函数（即总和、平均值等）来总结具有不规则时序的时间序列。但是，我目前的解决方案似乎效率低下且速度慢。

取聚合函数：

function aggArray = aggregate(array, groupIndex, collapseFn)

groups = unique(groupIndex, 'rows');
aggArray = nan(size(groups, 1), size(array, 2));

for iGr = 1:size(groups,1)
    grIdx = all(groupIndex == repmat(groups(iGr,:), [size(groupIndex,1), 1]), 2);
    for iSer = 1:size(array, 2)
      aggArray(iGr,iSer) = collapseFn(array(grIdx,iSer));
    end
end

end

请注意，array 和 groupIndex 都可以是二维的。 array 中的每一列都是要聚合的独立系列，但应将 groupIndex 的列放在一起（作为一行）以指定一个周期。

那么当我们给它带上一个不规则的时间序列时（注意第二个周期比它长一个基周期），时序结果很差：

a = rand(20006,10);
b = transpose([ones(1,5) 2*ones(1,6) sort(repmat((3:4001), [1 5]))]);

tic; aggregate(a, b, @sum); toc
Elapsed time is 1.370001 seconds.

使用分析器，我们可以发现grpIdx 行大约需要执行时间的 1/4 (.28 s)，iSer 循环大约需要总时间的 3/4 (1.17 s) ( 1.48 秒）。

将此与周期无关的情况进行比较：

tic; cumsum(a); toc
Elapsed time is 0.000930 seconds.

有没有更有效的方法来聚合这些数据？

---

计时结果

获取每个响应并将其放入一个单独的函数中，这是我在 Windows 7 和 Intel i7 上使用 Matlab 2015b 的 timeit 得到的计时结果：

    original | 1.32451
      felix1 | 0.35446
      felix2 | 0.16432
    divakar1 | 0.41905
    divakar2 | 0.30509
    divakar3 | 0.16738
matthewGunn1 | 0.02678
matthewGunn2 | 0.01977

澄清groupIndex

2D groupIndex 的一个示例是为一组涵盖 1980-2015 年的每日数据指定年号和周号：

a2 = rand(36*52*5, 10);
b2 = [sort(repmat(1980:2015, [1 52*5]))' repmat(1:52, [1 36*5])'];

因此，“年-周”期间由一行groupIndex 唯一标识。这可以通过调用unique(groupIndex, 'rows') 并获取第三个输出来有效处理，因此请随意忽略这部分问题。

Answer 1

方法#1

您可以创建与grIdx 对应的掩码 groups 与 bsxfun(@eq,..) 合二为一。现在，对于 collapseFn 和 @sum，您可以引入 matrix-multiplication，从而获得完全矢量化的方法，就像这样 -

M = squeeze(all(bsxfun(@eq,groupIndex,permute(groups,[3 2 1])),2))
aggArray = M.'*array

对于collapseFn 和@mean，你需要做更多的工作，如下所示-

M = squeeze(all(bsxfun(@eq,groupIndex,permute(groups,[3 2 1])),2))
aggArray = bsxfun(@rdivide,M,sum(M,1)).'*array

---

方法#2

如果您使用的是通用collapseFn，您可以使用通过先前方法创建的二维掩码M 来索引array 的行，从而将复杂性从O(n^2) 更改为@987654337 @。一些快速测试表明，这可以明显加快原始循环代码的速度。这是实现 -

n = size(groups,1);
M = squeeze(all(bsxfun(@eq,groupIndex,permute(groups,[3 2 1])),2));
out = zeros(n,size(array,2));
for iGr = 1:n
    out(iGr,:) = collapseFn(array(M(:,iGr),:),1);
end

请注意，collapseFn(array(M(:,iGr),:),1) 中的1 表示将沿其应用collapseFn 的维度，因此1 在那里必不可少。

---

奖金

顾名思义，groupIndex 似乎会保存整数值，可能滥用通过将groupIndex 的每一行视为索引元组来创建更有效的M，因此将groupIndex 的每一行转换为标量，最终得到groupIndex 的一维数组版本。这必须更有效，因为数据大小现在是0(n)。这个M 可以用于本文中列出的所有方法。所以，我们会有M 像这样 -

dims = max(groupIndex,[],1);
agg_dims = cumprod([1 dims(end:-1:2)]);
[~,~,idx] = unique(groupIndex*agg_dims(end:-1:1).'); %//'

m = size(groupIndex,1);
M = false(m,max(idx));
M((idx-1)*m + [1:m]') = 1;

Answer 2

Mex 函数 1

锤击时间：Mex function to crush it： 使用问题中的原始代码进行的基本案例测试在我的机器上花费了 1.334139 秒。恕我直言，2nd fastest answer from @Divakar 是：

groups2 = unique(groupIndex); 
aggArray2 = squeeze(all(bsxfun(@eq,groupIndex,permute(groups,[3 2 1])),2)).'*array; 

经过的时间是 0.589330 秒。

然后是我的 MEX 函数：

[groups3, aggArray3] = mg_aggregate(array, groupIndex, @(x) sum(x, 1));

经过的时间是 0.079725 秒。

测试我们得到相同的答案：norm(groups2-groups3) 返回 0 和 norm(aggArray2 - aggArray3) 返回 2.3959e-15。结果也与原始代码一致。

生成测试条件的代码：

array = rand(20006,10);
groupIndex = transpose([ones(1,5) 2*ones(1,6) sort(repmat((3:4001), [1 5]))]);

为了纯粹的速度，去墨西哥。如果编译 c++ 代码/复杂性的想法太痛苦了，请选择 Divakar 的答案。另一个免责声明：我没有对我的功能进行稳健的测试。

墨西哥方法 2

让我有些惊讶的是，在某些情况下，这段代码甚至比完整的 Mex 版本还要快（例如，在这个测试中大约需要 0.05 秒）。它使用mex function mg_getRowsWithKey 来计算组的索引。我认为这可能是因为我在完整的 mex 函数中复制数组的速度不如预期的快和/或调用“feval”的开销。它与其他版本的算法复杂度基本相同。

[unique_groups, map] = mg_getRowsWithKey(groupIndex);

results = zeros(length(unique_groups), size(array,2));

for iGr = 1:length(unique_groups)
   array_subset             = array(map{iGr},:);

   %// do your collapse function on array_subset. eg.
   results(iGr,:)           = sum(array_subset, 1);
end

当您使用array(groups(1)==groupIndex,:) 提取与整个组关联的数组条目时，您正在搜索整个 groupIndex 长度。如果您有数百万行条目，这将完全糟糕透顶。 array(map{1},:) 效率更高。

---

在折叠函数上调用“feval”仍然存在不必要的内存复制和其他开销。如果您在 C++ 中高效地实现聚合器函数，避免复制内存，则可能可以再实现 2 倍的加速。

Answer 3

聚会晚了一点，但使用accumarray 的单个循环会产生巨大的不同：

function aggArray = aggregate_gnovice(array, groupIndex, collapseFn)

  [groups, ~, index] = unique(groupIndex, 'rows');
  numCols = size(array, 2);
  aggArray = nan(numel(groups), numCols);
  for col = 1:numCols
    aggArray(:, col) = accumarray(index, array(:, col), [], collapseFn);
  end

end

在 MATLAB R2016b 中使用 timeit 对问题中的样本数据进行计时，结果如下：

original | 1.127141
 gnovice | 0.002205

超过 500 倍的加速！

Answer 4

去掉内循环，即

function aggArray = aggregate(array, groupIndex, collapseFn)

groups = unique(groupIndex, 'rows');
aggArray = nan(size(groups, 1), size(array, 2));

for iGr = 1:size(groups,1)
    grIdx = all(groupIndex == repmat(groups(iGr,:), [size(groupIndex,1), 1]), 2);
   aggArray(iGr,:) = collapseFn(array(grIdx,:));
end

并使用尺寸参数调用折叠函数

res=aggregate(a, b, @(x)sum(x,1));

已经提供了一些加速（在我的机器上是 3 倍）并避免了错误，例如sum 或 mean 产生，当他们遇到没有维度参数的单行数据然后跨列而不是标签折叠时。

如果您只有一个组标签向量，即所有数据列的组标签相同，您可以进一步加快速度：

function aggArray = aggregate(array, groupIndex, collapseFn)

ng=max(groupIndex);
aggArray = nan(ng, size(array, 2));

for iGr = 1:ng
    aggArray(iGr,:) = collapseFn(array(groupIndex==iGr,:));
end

后一个函数为您的示例提供相同的结果，速度提高了 6 倍，但无法处理每个数据列的不同组标签。

假设组索引的 2D 测试用例（此处还提供了 groupIndex 的 10 个不同列：

a = rand(20006,10);
B=[]; % make random length periods for each of the 10 signals
for i=1:size(a,2)
      n0=randi(10);
      b=transpose([ones(1,n0) 2*ones(1,11-n0) sort(repmat((3:4001), [1 5]))]);
      B=[B b];
end
tic; erg0=aggregate(a, B, @sum); toc % original method 
tic; erg1=aggregate2(a, B, @(x)sum(x,1)); toc %just remove the inner loop
tic; erg2=aggregate3(a, B, @(x)sum(x,1)); toc %use function below

经过的时间是 2.646297 秒。 经过的时间是 1.214365 秒。 经过的时间是 0.039678 秒 (!!!!)。

function aggArray = aggregate3(array, groupIndex, collapseFn)

[groups,ix1,jx] = unique(groupIndex, 'rows','first');
[groups,ix2,jx] = unique(groupIndex, 'rows','last');

ng=size(groups,1);
aggArray = nan(ng, size(array, 2));

for iGr = 1:ng
    aggArray(iGr,:) = collapseFn(array(ix1(iGr):ix2(iGr),:));
end

我认为这与不使用 MEX 的速度一样快。感谢 Matthew Gunn 的建议！ 分析表明，'unique' 在这里真的很便宜，并且只取出 groupIndex 中重复行的第一个和最后一个索引会大大加快速度。通过这次聚合迭代，我获得了 88 倍的加速。

Answer 5

嗯，我有一个几乎和 mex 一样快但只使用 matlab 的解决方案。 逻辑与上述大部分相同，创建一个虚拟二维矩阵，但我没有使用@eq，而是从一开始就初始化了一个逻辑数组。

我的经过时间是 0.172975 秒。 Divakar 的经过时间 0.289122 秒。

function aggArray = aggregate(array, group, collapseFn)
    [m,~] = size(array);
    n = max(group);
    D = false(m,n); 
    row = (1:m)';
    idx = m*(group(:) - 1) + row;
    D(idx) = true;
    out = zeros(m,size(array,2));
    for ii = 1:n
        out(ii,:) = collapseFn(array(D(:,ii),:),1);
    end
end