如何在不循环的情况下累积子矩阵（子阵列平滑）？

Question

在 Matlab 中，我需要累积一个大矩阵的重叠对角块。示例代码如下。

由于这段代码需要多次运行，所以会消耗大量资源。该过程用于阵列信号处理，用于所谓的子阵列平滑或空间平滑。有什么方法可以更快地做到这一点？

% some values for parameters
M = 1000; % size of array
m = 400; % size of subarray
n = M-m+1; % number of subarrays 
R = randn(M)+1i*rand(M);

% main code
S = R(1:m,1:m);
for i = 2:n
    S = S + R(i:m+i-1,i:m+i-1);
end

尝试：

1) 我尝试了以下替代矢量化版本，但不幸的是它变得慢了很多！

[X,Y] = meshgrid(1:m);
inds1 = sub2ind([M,M],Y(:),X(:));
steps = (0:n-1)*(M+1);
inds = repmat(inds1,1,n) + repmat(steps,m^2,1);
RR = sum(R(inds),2);
S = reshape(RR,m,m); 

2) 我使用 Matlab 编码器创建了一个 MEX 文件，它变得慢了很多！

Answer 1

我个人最近不得不整理我的代码的某些部分。完全不是专家，我建议尝试以下方法：

1) 矢量化：

摆脱for循环

S = R(1:m,1:m);
for i = 2:n
    S = S + R(i:m+i-1,i:m+i-1)
end

并将其替换为基于cumsum 的替代方案应该是这里的方法。

注意：将在未来的编辑中尝试并采用这种方法

2) 生成 MEX 文件：

在某些情况下，您可以简单地启动 Matlab Coder app（假设您在当前的 Matlab 版本中有它）。 这应该会为您生成一个.mex 文件，您可以调用它，因为它是您要替换的函数。

---

无论您选择 (1) 还是 2))，您都应该使用 tic; my_function(); toc; 分析您当前的实现，以获得相当数量的函数调用，并将其与您当前的实现进行比较：

my_time = zeros(1,10000);
for count = 1:10000
    tic;
    my_function();
    my_time(count) = toc;
end
mean(my_time)