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

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.

[Image Analyst]

What is the size exactly? Is Nxk more than 10 thousand or so?

[Anders]

Nk is around 4000

[Star Strider]

The obvious next question: what are N and k?

[Anders]

:-)

N is quite small with a value of 3 and k is around 1400...

[per isakson]

"Without looping" &nbsp - even if that is faster?

[Image Analyst]

That's a rather small array. The "size of the problem" is smaller than a run of the mill color digital image. Why are you afraid of a for loop in such a small array? I'd use a for loop just because it's so simple and intuitive. Maybe someone will come along and give you a one liner that uses some exotic matrix like toeplitz, hankle, circulant, or something. But a for loop to do a few million iterations should complete in less than a second.

[Anders]

The subroutine above is part of a function, which again is part of an optimization routine. So I was exactly hoping that someone came along with an exotic matrix to manipulate the blocks in order to save over-all run time. :-)

[Star Strider]

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]

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]

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.

[Anders]

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]

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]

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

[Image Analyst]

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]

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.

[Image Analyst]

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]

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]

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.

Answer 5

[Anders]

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

[Image Analyst]

I did sum the diagonals with my toeplitz method. Now you're saying something about blocks and submatrixes so I don't know what you're wanting to sum. Please give a small example with say a dozen rows and columns and 3 blocks and tell what you want summed.

[per isakson]

"% S is an (Ni x Ni)" &nbsp I assume that this comment is a mistake. It contradicts the comment, "% Collected in an [Ni x N]-matrix"

[Anders]

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

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.

Answer 6

[Swarup Rj]

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