时间序列聚合效率
Time series aggregation efficiency
我通常需要使用给定的聚合函数（即总和、平均值等）来汇总时间不规则的时间序列。但是，我目前的解决方案似乎效率低下且速度慢。
取聚合函数：
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 s)。
将此与周期无关的情况进行比较：
tic; cumsum(a); toc
Elapsed time is 0.000930 seconds.
是否有更有效的方法来汇总这些数据？
计时结果
获取每个响应并将其放入单独的函数中，这是我使用 timeit 和 Matlab 2015b 在 Windows 7 上使用 Intel i7 获得的计时结果：
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的澄清
二维 groupIndex 的一个例子是为一组涵盖 1980-2015 的每日数据指定年份和周数：
a2 = rand(36*52*5, 10);
b2 = [sort(repmat(1980:2015, [1 52*5]))' repmat(1:52, [1 36*5])'];
因此，"year-week" 句点由一行 groupIndex 唯一标识。这可以通过调用 unique(groupIndex, 'rows') 并获取第三个输出来有效处理，因此请随意忽略问题的这一部分。
取消内部循环，即
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 倍）并且避免了错误，例如求和或均值产生，当他们遇到没有维度参数的单行数据，然后跨列而不是标签折叠时。
如果你只有一个组标签向量，即所有数据列的组标签相同，你可以进一步加快速度：
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 测试用例（此处还提供了组索引的 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 倍的加速。
方法#1
您可以创建对应于 grIdx 的掩码
groups 与 bsxfun(@eq,..). Now, for collapseFn as @sum, you can bring in 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) 到 O(n)。一些快速测试表明，这可以显着加快原始循环代码的速度。这是实现 -
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 似乎可以保存整数值，可以 滥用 通过考虑每一行来更有效地 M 创建groupIndex 作为索引元组，从而将 groupIndex 的每一行转换为标量，最终得到 groupIndex 的一维数组版本。这必须更有效，因为数据大小现在为 0(n)。 M 可以提供给 post 中列出的所有方法。所以，我们会有 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;
Mex 函数 1
锤击时间：Mex function to crush it：
使用问题原始代码的基本案例测试在我的机器上花费了 1.334139 秒。恕我直言， 是：
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) returns 0 和 norm(aggArray2 - aggArray3) returns 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' 的 and/or 开销快。它的算法复杂度与另一个版本基本相同。
[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 倍的加速。
好吧，我有一个几乎与 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
聚会有点晚了，但是使用 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 倍的加速！
我通常需要使用给定的聚合函数（即总和、平均值等）来汇总时间不规则的时间序列。但是，我目前的解决方案似乎效率低下且速度慢。
取聚合函数：
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 s)。
将此与周期无关的情况进行比较：
tic; cumsum(a); toc
Elapsed time is 0.000930 seconds.
是否有更有效的方法来汇总这些数据？
计时结果
获取每个响应并将其放入单独的函数中，这是我使用 timeit 和 Matlab 2015b 在 Windows 7 上使用 Intel i7 获得的计时结果：
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的澄清
二维 groupIndex 的一个例子是为一组涵盖 1980-2015 的每日数据指定年份和周数：
a2 = rand(36*52*5, 10);
b2 = [sort(repmat(1980:2015, [1 52*5]))' repmat(1:52, [1 36*5])'];
因此，"year-week" 句点由一行 groupIndex 唯一标识。这可以通过调用 unique(groupIndex, 'rows') 并获取第三个输出来有效处理，因此请随意忽略问题的这一部分。
取消内部循环，即
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 倍）并且避免了错误，例如求和或均值产生，当他们遇到没有维度参数的单行数据，然后跨列而不是标签折叠时。
如果你只有一个组标签向量，即所有数据列的组标签相同，你可以进一步加快速度：
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 测试用例（此处还提供了组索引的 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 倍的加速。
方法#1
您可以创建对应于 grIdx 的掩码
groups 与 bsxfun(@eq,..). Now, for collapseFn as @sum, you can bring in 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) 到 O(n)。一些快速测试表明，这可以显着加快原始循环代码的速度。这是实现 -
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 似乎可以保存整数值，可以 滥用 通过考虑每一行来更有效地 M 创建groupIndex 作为索引元组，从而将 groupIndex 的每一行转换为标量，最终得到 groupIndex 的一维数组版本。这必须更有效，因为数据大小现在为 0(n)。 M 可以提供给 post 中列出的所有方法。所以，我们会有 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;
Mex 函数 1
锤击时间：Mex function to crush it：
使用问题原始代码的基本案例测试在我的机器上花费了 1.334139 秒。恕我直言，
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) returns 0 和 norm(aggArray2 - aggArray3) returns 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' 的 and/or 开销快。它的算法复杂度与另一个版本基本相同。
[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 倍的加速。
好吧，我有一个几乎与 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
聚会有点晚了，但是使用 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 倍的加速！