Extract sparse block diagonal without loop

Hello, I'm interested in pulling the block diagonal out of a sparse matrix and putting each of those diagonal matrices into the depth index of a 3D matrix. I can do this in a loop easily, but would like to avoid the loop if possible.

The following is an example of what I'm trying to do completed in a loop:

A = [1 3 4 5 0 0 0 0;
     6 7 8 9 0 0 0 0;
     3 4 1 2 0 0 0 0;
     0 0 0 0 9 9 8 7;
     0 0 0 0 5 4 3 2;
     0 0 0 0 1 9 7 7];

n = 2;  % number of block matrices in diagonal

[R, C] = size(A);

r = R/n;  % number of rows in each block matrix
c = C/n;  % number of columns in each block matrix

B(r, c, n) = 0;

i = 1; j = 1; k = 1;
while i < R && j < C
   B(:, :, k) = A(i:i+r-1, j:j+c-1);
   i = i + r;
   j = j + c;
   k = k + 1;
end

Is there a way to vectorize this assignment and make it quicker/more efficient?

Thanks!

Answer 1

[Stephen]

A slightly more simple (and convoluted) loop for doing the same thing:

ind = [1, 1];  % [i, j]
for k = 1:n
   B(:, :, k) = A( ind(1):ind(1) + r - 1, ind(2):ind(2) + c - 1 );
   ind = [ind(1) + r, ind(2) + c];
end

Help me dump this loop!

Answer 2

[Azzi Abdelmalek]

ii=A(A~=0);
B=reshape(ii,3,3,size(A,1)/3)

[John D'Errico]

No. Sorry, but this answer presumes that the blocks are strictly non-zero. This will throw an error if any of those blocks have any zeros in them.

[Stephen]

John, that's what I was thinking as well. Otherwise I believe it would work. Thanks for your input Azzi.

[Brian Wessel]

use kron to make the map of indices, then pull from your block diagonal matrix using the code above:

A  = [1 3 4 5 0 0 0 0;
      6 7 8 9 0 0 0 0;
      3 4 1 2 0 0 0 0;
      0 0 0 0 9 9 8 7;
      0 0 0 0 5 4 3 2;
      0 0 0 0 1 9 7 7];

blk_nrow = 3;
blk_ncol = 4;  
nblks    = size(A,1) / blk_nrow;  
      
kmap     = kron(eye(nblks), ones(blk_nrow,blk_ncol))
blk_list = A(kmap ~= 0);
B        = reshape(blk_list,blk_nrow,blk_ncol,size(A,1)/blk_nrow)

Answer 3

[Azzi Abdelmalek]

r=4;
c=3;
n=100;
[nr,mc]=size(A);
jdx=repmat(1:mc,r,1);
idx=permute(reshape(repmat(1:nr,c,1),c,r,[]),[2 1 3]);
ii=sub2ind(size(A),idx(:),jdx(:));
B=reshape(A(ii),r,c,[]);

Or

[n1,m1]=size(A);
h=bsxfun(@times,ones(1,r*c),(0:n-1)')'*r;
g=bsxfun(@plus,(0:m1-1)*n1,(1:r)');
V=reshape(A(g(:)+h(:)),r,c,n);

[Azzi Abdelmalek]

I think your for loop is faster