在MATLAB中对2D矩阵中的块求和

Question

我正在使用Matlab，并想知道如何在一个大矩阵中添加项。具体来说，我有一个4914x4914的矩阵，并希望创建一个189x189的矩阵，其中每个项都等于每个26x26子集中项的总和。

举个例子，假设我有以下魔术4x4矩阵：

[16 2  3  13; 5  11 10 8; 9  7  6  12; 4  14 15 1]

我希望创建一个与原始魔法4x4中每个2x2矩阵之和相等的2x2矩阵，即：

[(16+2+5+11)   (3+13+10+8);

(9+7+4+14)  (6+12+15+1)]

非常感谢您的帮助和建议！

谢谢，

jake

- jake_matlab_novice

---

如果您拥有图像处理工具箱，您可以使用 blockproc 来进行此操作。 - Dan

Answer 1

假设输入 4914x4914 矩阵为 A，以下是一个高效（运行时间）的方法 -

sublen = 26; %// subset length
squeeze(sum(reshape(sum(reshape(A,sublen,[])),size(A,1)/sublen,sublen,[]),2))

对于通用的块大小，我们可以编写一个函数 -

function out = sum_blocks(A,block_nrows, block_ncols)
out = squeeze(sum(reshape(sum(reshape(A,block_nrows,[])),...
                    size(A,1)/block_nrows,block_ncols,[]),2));
return

样例运行 -

>> A = randi(9,4,6);
>> A
A =
     8     2     4     9     4     5
     3     3     8     3     6     8
     9     6     6     7     1     9
     4     5     5     7     1     2
>> sum_blocks(A,2,3)
ans =
    28    35
    35    27
>> sum(sum(A(1:2,1:3)))
ans =
    28
>> sum(sum(A(1:2,4:6)))
ans =
    35
>> sum(sum(A(3:4,1:3)))
ans =
    35
>> sum(sum(A(3:4,4:6)))
ans =
    27

如果您想避免挤压 -

sum(permute(reshape(sum(reshape(A,sublen,[])),size(A,1)/sublen,sublen,[]),[1 3 2]),3)

基准测试

希望您关心性能，在此为发布的所有解决方案提供基准测试结果。以下是我使用的基准测试代码 -

num_runs = 100; %// Number of iterations to run benchmarks
A = rand(4914);
for k = 1:50000
    tic(); elapsed = toc(); %// Warm up tic/toc
end

disp('---------------------- With squeeze + reshape + sum')
tic
for iter = 1:num_runs
    sublen = 26; %// subset length
    out1 = squeeze(sum(reshape(sum(reshape(A,sublen,[])),...
                                   size(A,1)/sublen,sublen,[]),2));
end
time1 = toc;
disp(['Avg. elapsed time = ' num2str(time1/num_runs) ' sec(s)']), clear out1 sublen

disp('---------------------- With kron + matrix multiplication')
tic
for iter = 1:num_runs
    n = 189; k = 26;
    B = kron(speye(k), ones(1,n));
    result = B*A*B';
end
time2 = toc;
disp(['Avg. elapsed time = ' num2str(time2/num_runs) ' sec(s)']),clear result n k B

disp('---------------------- With accumarray')
tic
for iter = 1:num_runs
    s = 26; n = size(A,1)/s;
    subs = kron(reshape(1:(n^2), n, n),ones(s));
    out2 = reshape(accumarray(subs(:), A(:)), n, n);
end
time2 = toc;
disp(['Avg. elapsed time = ' num2str(time2/num_runs) ' sec(s)']),clear s n subs out2

我在我的系统上得到的基准测试结果 -

---------------------- With squeeze + reshape + sum
Avg. elapsed time = 0.050729 sec(s)
---------------------- With kron + matrix multiplication
Avg. elapsed time = 0.068293 sec(s)
---------------------- With accumarray
Avg. elapsed time = 0.64745 sec(s)

- Divakar

Answer 2

如果您没有图像处理工具箱，则可以使用accumarray来完成此操作：

s = 26;
n = size(A,1)/s;
subs = kron(reshape(1:(n^2), n, n),ones(s)); 
reshape(accumarray(subs(:), A(:)), n, n) 

如果您决定聚合一些方式而不是简单的求和，那么这将是可重用的，例如中位数：

reshape(accumarray(subs(:), A(:), [], @median), n, n) 

- Dan

---

+1 我也差不多要写类似的东西了！但我想知道是否有 kron 的替代方案，因为这个函数总是那么慢... - Robert Seifert

Answer 3

当然你可以使用矩阵乘法：

n = 26;
k = 189;
B = kron(speye(k), ones(1,n));
result = B*A*B';

- Mohsen Nosratinia

Answer 4

另一种方法是将整个矩阵重塑为4D矩阵，并对第一维和第三维的元素求和：

result = squeeze(sum(sum(reshape(A,26,189,26,189),1),3));

- Mohsen Nosratinia

---

重塑为4维..太棒了 :) - Santhan Salai