3
$\begingroup$

Given a large sparse (square) matrix in MATLAB, how can I extract the banded or the block-diagonal parts (of fixed size) of it efficiently?

These are useful operations when prototyping and testing preconditioners, or if you want a quick and dirty (but still relatively efficient) preconditioner. Sometimes these approximations are competitive against more sophisticated preconditioners, like ILU.

However, there are hardly any implementations available and the answers I could find require converting the sparse matrix into a full matrix first (which is not viable for large matrices) or employ for loops (it is slow - one of those cases that the JIT compiler cannot optimize yet).

$\endgroup$
3
$\begingroup$

(This answer is valid for both MATLAB and Octave, even though I mainly refer to MATLAB)

There are two beasts to slay; but let's first understand the underlying data structure. MATLAB and Octave store sparse matrices in the coordinate (COO) format, i.e. a sparse matrix S is a collection of three arrays of the length equal to the number of nonzero entries of S: row contains the row coordinates of the nonzero entries, col contains the column coordinates of the nonzero entries and finally val contains the values of the nonzero entries. For example, if row(i)=1, col(i)=10 and val(i)=3.14 then we know that the first row - tenth column entry of the sparse matrix S is 3.14 (which also is the i-th nonzero).

We can get this information from MATLAB using the find command:

[i,j,v] = find(S);

Using these arrays, it is trivial to find any banded part of S. We simply look for the entries which are at the distance of n or less from the main diagonal:

idx = find(abs(i-j)<n);

Then the banded part B of S can be constructed:

B = sparse(i(idx),j(idx),v(idx));

So it is just three steps, no conversion to a full matrix is required. These three steps can be replaced by two calls to spdiags: first, use spdiags to extract the diagonal and sub/superdiagonals, then use spdiags together with the extracted data to reconstruct the banded matrix. The spdiags approach works significantly faster in most cases, but the former approach will let us understand the answer to the second part better.

Notice that the n-block diagonal part of a matrix is contained in the n-banded part of the matrix. So a similar strategy may work to efficiently extract the n-block diagonal part of the matrix. The main observation is that depending on the row, there are some column entries we don't want. If it is the first row of the block, we don't want any of the entries to the left of the diagonal entry and if it is the last row of the block, we don't want any of the entries to the right of the diagonal entry. These changes linearly with respect to the row of the block. This should already remind us modular arithmetic, if mod(i(k),n)==0 we are at the last row of a block, so we should not accept any nonzero entry k s.t. j(k)>i(k). We can derive similar conditions for each row, here is a short (somehow cryptic) solution I came up with to summarize these rules:

d = j - i;
idx=find( (d< n-(mod(i-1,n))) .* (d>=-(mod(i-1,n))));

Now, the array idx contains the indices of all the nonzero entries which fall into any diagonal block. We can construct the block diagonal part of S exactly the same way as above:

B = sparse(i(idx),j(idx),v(idx));

I should mention that MATLAB also provides the command blkdiag. This command is more flexible in some sense, because it allows one to create block diagonal matrices with variable block sizes (steps above assume a fixed block size n). However, there is no built-in way to extract diagonal blocks of a sparse matrix. So one would have to extract the block themselves, and then use blkdiag. Even then it is not trivial to get a sparse matrix back; if the input is a cell-structure of dense matrices then MATLAB returns a dense matrix. You should make sure that at least one of them is sparse to guarantee that MATLAB gives you a sparse matrix back. Hence, to my best knowledge, the steps I described above are the best approach to do these operations. Someone might have it implemented in MEX which could run faster or be more memory-efficient, but I couldn't find it on the small section of the internet I looked at.

I would love to hear if anyone has ideas to improve this solution.

$\endgroup$
3
  • $\begingroup$ Did you just give a comprehensive answer to your own question? :-) $\endgroup$ Jun 10 at 2:02
  • $\begingroup$ @WolfgangBangerth, Yes :) But mainly because I couldn't find a community wiki for scicomp.SE to submit this piece of knowledge. Then I looked if it was OK to answer your own questions. meta.SE said it was encouraged (given that you obey common courtesy rules and not accept your own answer to increase your points). I think this may be useful to other people, as it is useful to me occasionally, and I don't want it to be forgotten and lost. Hopefully, my behaviour is okay on scicomp.SE. Otherwise, if there is a wiki, please let me know. The worst case, I can delete the question/answer. $\endgroup$ Jun 10 at 2:18
  • $\begingroup$ I don't find this objectionable :-) Interesting idea that one can accept one's own answer! $\endgroup$ Jun 10 at 4:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.