Example of the structure. In reality, might consist of a few thousands of blocks.
Create block sparse matrix
45 views (last 30 days)
Show older comments
How can I create a block sparse matrix in MATLAB using the 'sparse' funcion, such as S = sparse(I,J,Z)?
Instead of Z being a vector of non-zero entries, Z contains dense matrix blocks (e.g. 100*100). The size of the full matrix S can be very large, e.g 1 million * 1 million.
Thank you.
4 Comments
James Tursa
on 26 Jun 2019
How is the subdomain information stored? I think this could be done in a mex routine without sorting or copying the data more than once as long as the subdomain information is stored in a reasonable manner. The problem I see with the methods already posted is that they will require copying the data more than once and sorting the data as well, which will take extra time.
Answers (2)
David Goodmanson
on 24 Jun 2019
Edited: David Goodmanson
on 26 Jun 2019
HI ks,
Here is one way.
MODIFIED (see comment thread below)
showing two possible methods to fill the matrix by means of a for loop appoach. The second way is significantly faster, but both methods change the sparse matrix S as many times as you have blocks. For N blocks the overall time for the process appears to scale like N^2, so it's all right for a few hundred blocks. But when the numbers get large you will need to gin up some large index vectors and call sparse at most just a few times, preferably once as John has indicated.
The following may be useful.
From the [...] = ndgrid line in the code below, irow(:) is the row indices of the block itself (zero based). Then if m0vec is a column vector of N row indices for the upper left corner of the N blocks, and if you have a new enough version of Matlab
irowfull = irow(:) + m0vec';
irowfull = irowfull(:);
has the properties you want, similarly for column indices.
% create a matrix with known elements, not square
A = magic(4);
A(:,4) = [];
[Arow Acol] = size(A)
[irow icol] = ndgrid((0:Arow-1),(0:Acol-1)) % block indices
s_row = 10; s_col = 10; n_s = 100; % small example
S = spalloc(s_row, s_col,n_s)
m0 = 3; n0 = 5; % these are the row and column indices
% of the upper left corner of the block
% for loop here to index m0,n0, and A
% first way
S = S + sparse(m0+irow(:),n0+icol(:),A(:),s_row,s_col)
% second way, creating new sparse matrix not required
S(m0:m0+Arow-1,n0:n0+Acol-1) = A
% end
full(S)
A
4 Comments
David Goodmanson
on 24 Jun 2019
Edited: David Goodmanson
on 24 Jun 2019
Hi KS,
That works, but gratifyingly enough, after (defining S with spalloc) so does
S(m0:m0+Arow-1,n0:n0+Acol-1) = A;
which obviates having to use ndgrid and also obviates having to create a new (albeit small) sparse matrix each time through the loop. I modified the original answer but I don't know yet how this would do speedwise.
John D'Errico
on 24 Jun 2019
Edited: John D'Errico
on 24 Jun 2019
Are you asking how to create a block diagonal sparse matrix? Easy peasy.
tic
N = 100;
M = 10000;
Z = sparse(rand(N,N*M));
Zc = mat2cell(Z,N,repmat(N,1,M));
A = blkdiag(Zc{:});
toc
Elapsed time is 6.487005 seconds.
whos A
Name Size Bytes Class Attributes
A 1000000x1000000 1608000008 double sparse
So A is a block diagonal sparse matrix, of size 1e6x1e6, with 100x100 blocks on the diagonal, 10,000 such blocks. 6 seconds seems reasonable to build it, since almost 50% of that time was just in creating the original random matrix Z.
tic,Z = sparse(rand(N,N*M));toc
Elapsed time is 2.936146 seconds.
spy(A)
5 Comments
John D'Errico
on 26 Jun 2019
That time taken is very reasonable, to create a sparse matrix with 60 million entries. Without doing a test right now, the sparse block diagonal matrix I created in my answer had 100 million entries, taking also about 6 seconds to create. So I would not be at all surprised at your result. So 5 seconds or so seems reasonable.
See Also
Categories
Find more on Creating and Concatenating 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!