Asked by KSMAT21
on 23 Jun 2019

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.

Answer by David Goodmanson
on 24 Jun 2019

Edited by 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

David Goodmanson
on 24 Jun 2019

Hi ks

If you only have a few thousand blocks it makes sense to just use a for loop for each one rather than something more complicated.

I modified the example code by using m0 for the row index and n0 for the column index, and by using ndgrid instead of meshgrid. Ndgrid is the more natural choice here. Compared to ndgrid, the two outputs of meshgrid are swapped which is good for plotting purposes but not as good for index calculations.

It appears that at lot of your blocks have the same shape. Clearly ndgrid only has to be done once for all of those, and you can then do a for loop on the m0 and n0 indices. (m0 and n0 could be columns in an Nx2 array, addressed by the array row index). The sample code below takes m0 and n0 out of the meshgrid (now ndgrid) calculation.

Then in the for loop you have to address the different block matrices themselves, which will of course be much easier if a lot of them are identical.

A = magic(4);

A(:,4) = [];

[Arow Acol] = size(A)

[irow icol] = ndgrid((0:Arow-1),(0:Acol-1)) % block indices

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

S = sparse(m0+irow(:),n0+icol(:),A(:))

% end

full(S)

A

KSMAT21
on 24 Jun 2019

Thanks a lot, David.

I am correct to guess that the sparse matrix must be added to the previous iterations in the loop to create multiple blocks. Editing from your above code:

n = 1e4;

S = spalloc(n,n,3*n)

%for loop here to index m0,n0, and A

S = S + sparse(m0+irow(:),n0+icol(:),A(:),n,n)

% end

Is that correct?

Because, if not, each iteration in the for loop will clear the previous entries of the sparse matrix.

Regards,

KS

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.

Sign in to comment.

Answer by John D'Errico
on 24 Jun 2019

Edited by 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)

John D'Errico
on 25 Jun 2019

This really is far more trivial than you may think, BUT you ABSOLUTELY need to learn to use sparse.

For example, a 3(5)x3(5) block matrix, with 3 blocks.

blocksize = 5;

nblocks = 3;

% The elements contained in the various blocks.

% If each block is the same, it is easy to build.

% I've just used random elements.

Z = rand(blocksize,blocksize,nblocks);

[subr,subc] = meshgrid(1:blocksize);

% The pattern matrix:

blockrc = [1 2;2 1;3 3];

rowind = subr + reshape((blockrc(:,1)-1)*blocksize,[1 1 nb]);

colind = subc + reshape((blockrc(:,2)-1)*blocksize,[1 1 nb]);

% create A in one call to sparse

A = sparse(rowind(:),colind(:),Z(:),blocksize*max(blockrc(:,1)),blocksize*max(blockrc(:,2)));

whos A

Name Size Bytes Class Attributes

A 15x15 1328 double sparse

spy(A)

So A is 15x15 sparse, with 3 blocks of size 5x5 each. Change the code to match your problem, in that all you need to do is set up your own pattern matrix, and assign the matrix Z with the elemnts that it contains. But see how trivial it really is to do, in ONE call to sparse. Again, learn to use sparse, and your code will be easy to write.

KSMAT21
on 26 Jun 2019

Thanks, this works great. Only thing, I used ndgrid instead of meshgrid as explained in one of David's comment.

The sparse function took 5.97 seconds created a sparse matrix of 60 millions non-zero entries. Is the time taken resonable for this size of matrix?

John D'Errico
on 26 Jun 2019

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 4 Comments

## KSMAT21 (view profile)

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/468511-create-block-sparse-matrix#comment_717871

## James Tursa (view profile)

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/468511-create-block-sparse-matrix#comment_717878

## KSMAT21 (view profile)

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/468511-create-block-sparse-matrix#comment_717882

## James Tursa (view profile)

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/468511-create-block-sparse-matrix#comment_718386

Sign in to comment.