将功能应用于滚动窗口

Question

说我有一个长长的清单A值(比如长度为1000),为此我要计算std成对的100,也就是我想计算std(A(1:100)),std(A(2:101)),std(A(3:102)),..., std(A(901:1000)).

在Excel/VBA中,可以通过=STDEV(A1:A100)在一个单元格中写入然后一次填充来轻松实现此目的.现在我的问题是,如何在Matlab中有效地实现这一点,而无需使用任何昂贵的for循环.

编辑:是否也可以为时间序列列表执行此操作,例如,当A的尺寸为1000 x 4(即4个时间序列长度为1000)时？输出矩阵的尺寸应为901 x 4.

Answer 1

注意:要获得最快的解决方案,请参阅Luis Mendo的回答

所以首先使用for循环(特别是如果那些是你的实际尺寸)真的不会很贵.除非你使用的是旧版本的Matlab,否则JIT编译器(当然还有预分配)会使循环变得便宜.

其次 - 你试过循环吗？因为在开始过早优化之前,你应该先尝试一下天真的实现.

第三 - arrayfun可以使它成为一个单行,但它基本上只是一个带有额外开销的for循环,并且很可能比for循环慢,如果速度确实是你的关注.

最后一些代码:

n = 1000;
A = rand(n,1);
l = 100;

for loop(几乎不笨重,可能效率很高):

S = zeros(n-l+1,1);  %//Pre-allocation of memory like this is essential for efficiency!
for t = 1:(n-l+1)
    S(t) = std(A(t:(t+l-1)));
end

矢量化(内存效率不高!)解决方案:

[X,Y] = meshgrid(1:l)
S = std(A(X+Y-1))

一个可能更好的矢量化解决方案(和一个单行)但仍然没有内存效率:

S = std(A(bsxfun(@plus, 0:l-1, (1:l)')))

请注意,使用所有这些方法,您可以替换std为任何函数,只要它适用于矩阵的列(这是Matlab中的标准)

去2D:

要去2D,我们需要去3D

n = 1000;
k = 4;
A = rand(n,k);
l = 100;

ind = bsxfun(@plus, permute(o:n:(k-1)*n, [3,1,2]), bsxfun(@plus, 0:l-1, (1:l)'));    %'
S = squeeze(std(A(ind)));
M = squeeze(mean(A(ind)));
%// etc...

要么

[X,Y,Z] = meshgrid(1:l, 1:l, o:n:(k-1)*n);
ind = X+Y+Z-1;
S = squeeze(std(A(ind)))
M = squeeze(mean(A(ind)))
%// etc...

要么

ind = bsxfun(@plus, 0:l-1, (1:l)');                                                  %'
for t = 1:k
    S = std(A(ind));
    M = mean(A(ind));
    %// etc...
end

或者(取自Luis Mendo的回答 - 在他的回答中注意到他显示了这个简单循环的更快的替代方案)

S = zeros(n-l+1,k);
M = zeros(n-l+1,k);
for t = 1:(n-l+1)
    S(t,:) = std(A(k:(k+l-1),:));
    M(t,:) = mean(A(k:(k+l-1),:));
    %// etc...
end

Answer 2

你正在做的基本上是一个过滤操作.

如果您有权访问图像处理工具箱,

stdfilt(A,ones(101,1)) %# assumes that data series are in columns

会做的伎俩(无论是哪个维度A).请注意,如果您还可以访问并行计算工具箱,则可以让这些过滤器操作在GPU上运行,尽管您的问题可能太小而无法产生明显的加速.

Answer 3

为了最大限度地减少操作次数,您可以利用标准偏差可以计算为涉及第二和第一时刻的差异的事实,

通过累积总和(使用cumsum)有效地获得滚动窗口上的时刻:

A = randn(1000,4); %// random data
N = 100; %// window size
c = size(A,2);
A1 = [zeros(1,c); cumsum(A)];
A2 = [zeros(1,c); cumsum(A.^2)];
S = sqrt( (A2(1+N:end,:)-A2(1:end-N,:) ...
    - (A1(1+N:end,:)-A1(1:end-N,:)).^2/N) / (N-1) ); %// result

标杆

这是基于循环的解决方案的比较,使用timeit.循环方法与Dan的解决方案一样,但适用于2D情况,std以矢量化方式展示沿每列工作的事实.

%// File loop_approach.m
function S = loop_approach(A,N);
[n, p] = size(A);
S = zeros(n-N+1,p);
for k = 1:(n-N+1)
    S(k,:) = std(A(k:(k+N-1),:));
end

%// File bsxfun_approach.m
function S = bsxfun_approach(A,N);
[n, p] = size(A);
ind = bsxfun(@plus, permute(0:n:(p-1)*n, [3,1,2]), bsxfun(@plus, 0:n-N, (1:N).')); %'
S = squeeze(std(A(ind)));

%// File cumsum_approach.m
function S = cumsum_approach(A,N);
c = size(A,2);
A1 = [zeros(1,c); cumsum(A)];
A2 = [zeros(1,c); cumsum(A.^2)];
S = sqrt( (A2(1+N:end,:)-A2(1:end-N,:) ...
    - (A1(1+N:end,:)-A1(1:end-N,:)).^2/N) / (N-1) );

%// Benchmarking code
clear all
A = randn(1000,4); %// Or A = randn(1000,1);
N = 100;
t_loop   = timeit(@() loop_approach(A,N));
t_bsxfun = timeit(@() bsxfun_approach(A,N));
t_cumsum = timeit(@() cumsum_approach(A,N));
disp(' ')
disp(['loop approach:   ' num2str(t_loop)])
disp(['bsxfun approach: ' num2str(t_bsxfun)])
disp(['cumsum approach: ' num2str(t_cumsum)])
disp(' ')
disp(['bsxfun/loop gain factor: ' num2str(t_loop/t_bsxfun)])
disp(['cumsum/loop gain factor: ' num2str(t_loop/t_cumsum)])

结果

我正在使用Matlab R2014b,Windows 7 64位,双核处理器,4 GB RAM:

• 4柱案例:

  loop approach:   0.092035
  bsxfun approach: 0.023535
  cumsum approach: 0.0002338
  
  bsxfun/loop gain factor: 3.9106
  cumsum/loop gain factor: 393.6526
  

• 单柱情况:

  loop approach:   0.085618
  bsxfun approach: 0.0040495
  cumsum approach: 8.3642e-05
  
  bsxfun/loop gain factor: 21.1431
  cumsum/loop gain factor: 1023.6236
  

因此,cumsum基于-Based的方法似乎是最快的:比4列情况下的循环快约400倍,而单列情况下快1000倍.

Answer 4

Matlab 中有几个函数可以有效地完成这项工作。

一方面，您可以使用colfilt或 等函数nlfilter，它对滑块执行计算。colfilt比 更有效nlfilter，但仅当块内元素的顺序无关紧要时才可以使用。以下是如何在您的数据上使用它：

S = colfilt(A, [100,1], 'sliding', @std);

或者

S = nlfilter(A, [100,1], @std);

在您的示例中，您可以清楚地看到性能差异。但有一个技巧：两个函数都会填充输入数组，以便输出向量与输入数组具有相同的大小。要仅获取输出向量的相关部分，您需要跳过第一个floor((100-1)/2) = 49元素并获取1000-100+1值。

S(50:end-50)

但还有另一种解决方案，接近colfilt，效率更高。colfilt调用col2im将输入向量重塑为矩阵，并在每个不同的列上应用给定的函数。这会将大小为 [1000,1] 的输入向量转换为大小为 [100,901] 的矩阵。但是colfilt用 0 或 1 填充输入数组，您就不需要它了。因此，您可以colfilt在没有填充步骤的情况下运行，然后应用于std每一列，这很容易，因为std应用于矩阵会返回列标准差的行向量。最后，如果需要，可以将其转置以获得列向量。简而言之：

S = std(im2col(X,[100 1],'sliding')).';

备注：如果您想应用更复杂的功能，请参见colfilt第 144 和 147 行的代码（适用于 v2013b）。