Extract sparse block diagonal without loop

4 次查看（过去 30 天）
Stephen
Stephen 2014-4-15
编辑： Brian Wessel 2020-1-9
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!

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

回答（3 个）

Stephen
Stephen 2014-4-15
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!
Azzi Abdelmalek
Azzi Abdelmalek 2014-4-15
ii=A(A~=0);
B=reshape(ii,3,3,size(A,1)/3)
  3 个评论
显示 1更早的评论
Stephen
Stephen 2014-4-15
John, that's what I was thinking as well. Otherwise I believe it would work. Thanks for your input Azzi.
Brian Wessel
Brian Wessel 2020-1-9
编辑：Brian Wessel 2020-1-9
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)

请先登录，再进行评论。

Azzi Abdelmalek
Azzi Abdelmalek 2014-4-15
编辑：Azzi Abdelmalek 2014-4-15
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);

类别

在 Help Center 和 File Exchange 中查找有关 Creating and Concatenating Matrices 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by