How to extract block matrices along the diagonal entries without looping?
43 views (last 30 days)
Show older comments
I have a square matrix, A with dimension (N*k x N*k), and I want to extract the sum of all the diagonal entries of the block matrices (each block, Aii, with dimension NxN).
A =
[A11, A12, ..., A1k
A21, A22, ..., A2k
...
Ak1, Ak2, ..., Akk]
Then I want to obtain the sum of the main diagonal (A11 + A22 + ...), and (A21 + A32 + ...) and so on...
A is symmetric in the blocks, i.e. A12 = A21, so lower/upper entries suffice.
If you have a way to extract the sum just for the main diagonal it will be excellent as well, since deleting block rows and block columns will get me the desired diagonal.
Multiple loops is unfortunately not an option for me due to the size of the problem.
9 Comments
Star Strider
on 27 Jul 2014
So you have a matrix of 1400 (3x3) submatrices, and you want the sums of the (main) diagonals of the submatrices, probably returned as a matrix, or do I have it backwards?
Accepted Answer
per isakson
on 28 Jul 2014
Edited: per isakson
on 29 Jul 2014
Hopefully, I understood the task better this time.
- This function takes a bit less than three seconds to run on my old desktop (R2013a,64bit,Win7).
- Replacing the nested function by a subfunction and passing all arguments needed provides the same performance
- A double loop in the main function with the help variable ss - same performance again.
- However, a double loop without ss., i.e. working directly with S((RR-1)*N+(1:N),(1:N)), triples the execution time.
function S = cssm_v2
N = 3; % size of each submatrix
K = 1400; % number of submatrices in each direction
S = zeros( N*K, N );
% Create data, which gives me a chance to see if the result
% is as expected
X = zeros( K );
xx = 1 : K;
for RR = 1 : K
X( RR, : ) = RR*10 + xx;
end
A = kron( X, ones(N) );
% Solution based on loops
for RR = 1 : K
S( (RR-1)*N+(1:N) , (1:N) ) = sum_( );
end
function ss = sum_( )
ss = zeros( N );
for jj = RR : K
rr = (jj-1)*N+(1:N);
cc = rr - (RR-1)*N;
ss = ss + A( rr, cc );
end
end
end
I won't try to vectorize this.
2 Comments
Pierre Benoit
on 31 Jul 2014
Hi,
I have a slightly better solution if you want : (1,2sec vs 2sec for Per isakson's solution on my computer) Thanks to Per idea to make a subfunction for the sum, it gave better results on my computer.
I store all the block matrices only of the upper matrice in a 3-dimensionnal array and then proceed to do the sums needed.
function S = sum_diag
N = 3; % size of each submatrix
k = 1400; % number of submatrices in each direction
S = zeros( N*k,N );
A = ones( N*k, N*k );
% Create a 3-dimensionnal array with only blocks from the upper
% matrice, row wise
B = zeros(N,N,k*k);
for ii=1:k
temp = (ii-1)*k;
B(:,:,temp+1:temp + (k-ii+1) ) = reshape( A(N*(ii-1)+1 : N*ii , (ii-1)*N+1 : end) , [N N (k-ii+1)]);
end
for ii = 1:k
S(N*(ii-1)+1:N*ii,:) = sum_diag_sub();
end
function ss = sum_diag_sub()
ss = zeros(N,N);
for jj = 1:(k-ii+1)
ss = ss + B(:,:,ii+(jj-1)*k);
end
end
end
More Answers (5)
per isakson
on 28 Jul 2014
Edited: per isakson
on 28 Jul 2014
This function takes 3 seconds on my machine and I guess less than half of that on the machine of Image.
function S = cssm()
N = 3;
k = 1400;
A = rand( N*k, N*k );
S = zeros( N, k, k );
cc_blk = 0;
for cc = 1 : N : N*k
cc_blk = cc_blk + 1;
rr_blk = 0;
for rr = 1 : N : cc
rr_blk = rr_blk + 1;
S(:,rr_blk,cc_blk) = [ A(rr ,cc)+A(rr+1,cc+1)+A(rr+2,cc+2)
A(rr+1,cc)+A(rr+2,cc+1)
A(rr+2,cc) ];
end
end
end
1 Comment
Image Analyst
on 28 Jul 2014
Edited: Image Analyst
on 28 Jul 2014
On my Dell M6800:
Elapsed time is 1.102767 seconds.
Does this get the sum of the diagonals all the way down to the very bottom of the matrix like my toeplitz solution does?
Image Analyst
on 27 Jul 2014
Why not simply do
% Create sample data
A = rand(N*k, N*k);
N = 3;
k = 1400;
% Initialize mask
mask = zeros(N*k, N*k);
tic; % Start timer
for row = 1 : N : N*k
r1 = row;
r2 = r1 + N - 1;
mask(r1:r2, r1:r2) = 1;
end
% Extract from A
A_masked = A .* mask;
toc; % Takes 0.05 seconds on my computer.
It's simple and fast - only 0.05 seconds on my computer even with the sizes you gave.
2 Comments
Image Analyst
on 28 Jul 2014
After seeing per's answer and reading yours again, I think this might not be what you wanted. I had thought that k was the smaller dimension and you had a bunch of blocks k wide placed along the main diagonal.
Image Analyst
on 28 Jul 2014
Edited: Image Analyst
on 28 Jul 2014
An alternative way:
tic
N = 3;
k = 1400;
A = rand( N*k, N*k );
c = 1:k;
r = 1:k;
% Get the sums going down from the top row.
t = toeplitz(c,r);
upper_t = triu(t);
for c = 1 : k
diagVector = A(upper_t == c);
theDiagSumsTop(c) = sum(diagVector);
end
toc
Takes 4.6 seconds though.
0 Comments
Anders
on 28 Jul 2014
Edited: Anders
on 28 Jul 2014
4 Comments
Image Analyst
on 29 Jul 2014
It's not clear. All I asked for was a simple example array, like
m=[
1 2 3 0 0 0 0 0 0;...
2 1 5 0 0 0 0 0 0;...
3 5 1 0 0 0 0 0 0;...
0 0 0 9 8 3 0 0 0;...
0 0 0 8 5 2 0 0 0;...
0 0 0 3 2 7 0 0 0;...
0 0 0 0 0 0 7 3 1;...
0 0 0 0 0 0 3 4 5;...
0 0 0 0 0 0 1 5 1]
and then you tell me if the second sum is actually 3 sums: 2+5 (just over the first block), then 8+2, then 3+5, OR if it's just a single sum of 2+5+0+8+2+0+3+5 (diagonal crosses all 3 blocks). But if you don't want to, that's fine. Good luck. By the way, the toeplitz solution I gave covers all the blocks, not just one of them.
Swarup Rj
on 19 Jan 2015
Is it same as " FINDING THE SUM OF DIAGONAL BLOCKS OF A MATRIX."
I am solving an optimization problem on Clustering.
If you have the code .. do help.. :-)
0 Comments
See Also
Categories
Find more on Operating on Diagonal Matrices in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!