10

我想对以下 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 <= r <= 20 ）。convnfft.m 声明 r 应该 > 20 才能体现任何好处。

4

4 回答 4

11

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

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

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

如果 C 和 r 很大，您可能需要查看 Matlab File Exchange 中的CONVNFFT以加快计算速度。

于 2010-05-19T11:58:06.940 回答
3

基于JS的想法，正如Jonas在评论中指出的那样，这可以通过IM2COL一些数组操作在两行中完成：

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

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

将此与 Jonas 的基于卷积的解决方案进行比较，这不执行任何矩阵乘法，仅索引...

于 2010-05-21T20:56:52.653 回答
2

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

编辑

我没有设法找到一种方法来干净地对其进行矢量化，但我确实找到了 function 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]);

create 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和其他方法来提高效率；为了清楚起见，我选择在这里进行布局。不过，我还没有对此进行基准测试，以查看它与双循环方法的比较；我得先回家吃饭！

于 2010-05-19T09:08:05.717 回答
2

这是你要找的吗？这个函数添加对角线并将它们放入一个向量中，类似于函数“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

干杯

于 2013-10-17T05:43:58.767 回答