How to extract block matrices along the diagonal entries without looping?

Question

Anders 2014-7-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/143535-how-to-extract-block-matrices-along-the-diagonal-entries-without-looping

回答： Swarup Rj 2015-1-19

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 个评论
显示 7更早的评论隐藏 7更早的评论

Star Strider 2014-7-27

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?

Anders 2014-7-28

No, that's exactly my question. Preferably returned as an (Nk x N) matrix (vector of blocks - where each block is (NxN) submatrix, which is the sum of the diagonal blocks).

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

per isakson 2014-7-28

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/143535-how-to-extract-block-matrices-along-the-diagonal-entries-without-looping#answer_146697

编辑：per isakson 2014-7-29

在 MATLAB Online 中打开

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 个评论
显示无隐藏无

Anders 2014-7-31

Thank you for the solution. :-) I modified it to my code, and as you point out, the ss-subfunction slows the execution. I'm gonna go with the double-loop solution in the main.

Pierre Benoit 2014-7-31

在 MATLAB Online 中打开

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

请先登录，再进行评论。

Answer 2

per isakson 2014-7-28

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/143535-how-to-extract-block-matrices-along-the-diagonal-entries-without-looping#answer_146584

编辑：per isakson 2014-7-28

在 MATLAB Online 中打开

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 个评论
显示 -1更早的评论隐藏 -1更早的评论

Image Analyst 2014-7-28

编辑：Image Analyst 2014-7-28

在 MATLAB Online 中打开

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?

请先登录，再进行评论。

Answer 3

Image Analyst 2014-7-27

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/143535-how-to-extract-block-matrices-along-the-diagonal-entries-without-looping#answer_146577

在 MATLAB Online 中打开

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 个评论
显示无隐藏无

Image Analyst 2014-7-28

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.

Anders 2014-7-28

编辑：Anders 2014-7-28

It's part of the way though. I still need to extract the matrix sum of the block diagonal of A_masked. And again, I need to do this for each i =1,...,k if I want to extract all (lower or upper) diagonal blocks from A. So perhaps this is not the most efficient approach?

请先登录，再进行评论。

Answer 4

Image Analyst 2014-7-28

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/143535-how-to-extract-block-matrices-along-the-diagonal-entries-without-looping#answer_146588

编辑：Image Analyst 2014-7-28

在 MATLAB Online 中打开

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 5

Anders 2014-7-28

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/143535-how-to-extract-block-matrices-along-the-diagonal-entries-without-looping#answer_146614

编辑：Anders 2014-7-28

在 MATLAB Online 中打开

Thank you both for the effort. But I think I haven't expressed myself clear enough.

@per isakson: This code of yours sums the diagonals inside each submatrix, but I need to have a sum of the submatrices instead. I have included a small sample to illustrate which sums I am interested in:

    N = 3; % size of each submatrix 
    i = 3; % number of submatrices in each direction
    phi = rand(3,3);
    sigma = rand(3,3);
    O = zeros(N*i, N);
      for s = 1:i
        O(s*N-(N-1):N*s,1:N) = (eye(N) - phi)^(i-1-(s-1)) * sigma;
        % (I - phi) decreases in power at every step starting from i-1,   
      end
       A = O * O'; % Multiplied together to get all the cross-term
                  % S is an (Ni x Ni) 
     % I want to extract these diagonal sums in the lower diagonal (NxN- matrices) 
     % preferably in a Ni X N matrix:
     S_0 = A(1:3, 1:3) + A(4:6, 4:6) + A(7:9, 7:9); % main diagonal     
     S_1 = A(4:6, 1:3) + A(7:9, 4:6); % -1 diagonal block
     S_2 = A(7:9, 1:3); % -2 diagonal block
     S = [S_0; S_1; S_2]; % Collected in an [Ni x N]-matrix

4 个评论
显示 2更早的评论隐藏 2更早的评论

Anders 2014-7-28

@Image Analyst

Yes, my query is about submatrices/blocks. Perhaps this was not all clear from the start. The output should be sums of (NxN)-submatrices from the matrix A above. And when I say along the diagonal, I meant along the diagonal-blocks (treating blocks as elements) and not along the elements within each submatrix. It should be clear from the small example above, where the S_i's (i = 0,1,2) are the sums in question.

@ per isakson Yes, that's a typo. I meant "A is an (Ni x Ni) - matrix".

The matrix, S, could just be one possible way to report the sums of the diagonal submatrices - if you have a better way, i don't mind...

Image Analyst 2014-7-29

在 MATLAB Online 中打开

It's not clear. All I asked for was a simple example array, like

m=[
2 3 0 0 0 0 0 0;...
1 5 0 0 0 0 0 0;...
5 1 0 0 0 0 0 0;...
0 0 9 8 3 0 0 0;...
0 0 8 5 2 0 0 0;... 
0 0 3 2 7 0 0 0;...
0 0 0 0 0 7 3 1;...
0 0 0 0 0 3 4 5;...
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.

请先登录，再进行评论。