向量化矩阵中不同对角线的和

Question

我想对以下 MATLAB 代码进行矢量化处理。我认为它一定很简单，但我发现它仍然令人困惑。

r = some constant less than m or n
[m,n] = size(C);
S = zeros(m-r,n-r);
for i=1:m-r+1
    for j=1:n-r+1
        S(i,j) = sum(diag(C(i:i+r-1,j:j+r-1)));
    end
end

代码从另一个分数表 C 为动态规划算法计算分数表 S。
对角线求和是为用于生成 C 的各个数据片段生成分数，用于所有可能的片段（大小为 r）。

提前感谢您的任何回答！对不起，如果这个应该很明显......

注意
内置的 conv2 比 convnfft 更快，因为我的 eye(r) 非常小（ 5 20 才能体现任何好处。

Answer 1

我认为您可能需要将 C 重新排列为 3D 矩阵，然后再将其沿其中一个维度求和。我会尽快发布答案。

编辑

我没有设法找到一种干净矢量化它的方法，但我确实找到了函数accumarray，这可能会有所帮助。等我回家再详细看看。

编辑#2

通过使用线性索引找到了一个更简单的解决方案，但这可能会占用大量内存。

在C(1,1)处，我们要求和的索引是1+[0, m+1, 2*m+2, 3*m+3, 4*m+4, ... ]，或 (0:r-1)+(0:m:(r-1)*m)

sum_ind = (0:r-1)+(0:m:(r-1)*m);

创建S_offset，一个 (mr) by (nr) by r 矩阵，使得 S_offset(:,:,1) = 0, S_offset(:,:,2) = m+1, S_offset(:, :,3) = 2*m+2，以此类推。

S_offset = permute(repmat( sum_ind, [m-r, 1, n-r] ), [1, 3, 2]);

创建S_base，一个基数组地址矩阵，将从中计算偏移量。

S_base = reshape(1:m*n,[m n]);
S_base = repmat(S_base(1:m-r,1:n-r), [1, 1, r]);

最后，使用S_base+S_offset 来处理C 的值。

S = sum(C(S_base+S_offset), 3);

你当然可以使用bsxfun等方法来提高效率；为了清楚起见，我选择在这里进行布局。不过，我还没有对此进行基准测试，以查看它与双循环方法的比较；我得先回家吃饭了！

Answer 2

如果我理解正确，您正在尝试计算 C 的每个子数组的对角线和，其中您已删除 C 的最后一行和列（如果您不应该删除行/列，则需要循环到m-r+1，你需要将整个数组 C 传递给我下面解决方案中的函数）。

您可以通过convolution 执行此操作，如下所示：

S = conv2(C(1:end-1,1:end-1),eye(r),'valid');

如果 C 和 r 很大，您可能需要查看来自 Matlab 文件交换的 CONVNFFT 以加快计算速度。

Answer 3

基于JS 的想法，正如Jonas 在cmets 中指出的那样，这可以使用IM2COL 和一些数组操作分两行完成：

B = im2col(C, [r r], 'sliding');
S = reshape( sum(B(1:r+1:end,:)), size(C)-r+1 );

基本上B包含矩阵C上所有大小为r×r的滑块的元素。然后我们取每个块 B(1:r+1:end,:) 对角线上的元素，计算它们的总和，并将结果重新整形为预期的大小。

---

与 Jonas 的基于卷积的解决方案相比，它不执行任何矩阵乘法，仅索引...

Answer 4

这就是你要找的吗？此函数将对角线相加并将它们放入一个向量中，类似于函数“sum”将矩阵中的所有列相加并将它们放入一个向量中。

function [diagSum] = diagSumCalc(squareMatrix, LLUR0_ULLR1)
% 
% Input: squareMatrix: A square matrix.
%        LLUR0_ULLR1:  LowerLeft to UpperRight addition = 0     
%                      UpperLeft to LowerRight addition = 1
% 
% Output: diagSum: A vector of the sum of the diagnols of the matrix.
% 
% Example: 
% 
% >> squareMatrix = [1 2 3; 
%                    4 5 6;
%                    7 8 9];
% 
% >> diagSum = diagSumCalc(squareMatrix, 0);
% 
% diagSum = 
% 
%       1 6 15 14 9
% 
% >> diagSum = diagSumCalc(squareMatrix, 1);
% 
% diagSum = 
% 
%       7 12 15 8 3
% 
% Written by M. Phillips
% Oct. 16th, 2013
% MIT Open Source Copywrite
% Contact mphillips@hmc.edu fmi.
% 

if (nargin < 2)
    disp('Error on input. Needs two inputs.');
    return;
end

if (LLUR0_ULLR1 ~= 0 && LLUR0_ULLR1~= 1)
    disp('Error on input. Only accepts 0 or 1 as input for second condition.');
    return;
end

[M, N] = size(squareMatrix);

if (M ~= N)
    disp('Error on input. Only accepts a square matrix as input.');
    return;
end

diagSum = zeros(1, M+N-1);

if LLUR0_ULLR1 == 1
    squareMatrix = rot90(squareMatrix, -1);
end

for i = 1:length(diagSum)
    if i <= M
        countUp = 1;
        countDown = i;
        while countDown ~= 0
            diagSum(i) = squareMatrix(countUp, countDown) + diagSum(i);
            countUp = countUp+1;
            countDown = countDown-1;
        end
    end
    if i > M
        countUp = i-M+1;
        countDown = M;
        while countUp ~= M+1
            diagSum(i) = squareMatrix(countUp, countDown) + diagSum(i);
            countUp = countUp+1;
            countDown = countDown-1;
        end
    end
end

干杯