Sum over blocks in a 2D matrix - MATLAB

Question

I'm working on Matlab and was wondering how I add terms within a large matrix. Specifically, I have a 4914x4914 matrix and would like to create a 189x189 matrix, where each term is equal to the sum of the terms in each 26x26 subset.

To illustrate, say I had the magic 4x4 matrix as follows:

[16 2  3  13;

 5  11 10 8;

 9  7  6  12;

 4  14 15 1]

and I wanted to create a 2x2 matrix equal to the sum of each 2x2 matrix within the original magic 4x4, i.e.:

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

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

Grateful for any advice! Thanks jake

Answer 1

Assuming A to be the input 4914x4914 matrix, this could be an efficient (in terms of runtime) approach -

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

For a generic block size, let's have a function -

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

Sample run -

>> 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

If you would like to avoid squeeze -

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

---

Benchmarking

Hoping you would care about performance, here are the benchmark results for all the solutions posted here. The benchmarking code that I have used -

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

The benchmarks results I got on my system -

---------------------- 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)

Answer 2

If you don't have the image processing toolbox then you can do this using accumarray:

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

this is reusable should you decide to aggregate some way other than a simple sum e.g. a median:

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

Answer 3

You can use matrix multiplication, of course:

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

Answer 4

An alternative way is to reshape the whole matrix into a 4D matrix and sum the elements over first and third dimension:

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