Block diagonal and sparse matrices

Specific Domains Numerics
1

I have a code that does something like this:

Z = inv(ABCDE)

Here, A is a diagonal matrix. B is block diagonal, its diagonal blocks are 2x2. C is block diagonal, its diagonal blocks are 4x4. Subsequent matrices have larger blocks, but the whole thing is “compatible” so that ABCDE is indeed block diagonal, and the inverse Z is therefore also block diagonal.

In MATLAB, I’d use sparse matrices throughout. In my experience, if a sparse matrix is block diagonal, MATLAB seems to do the right thing.

In Julia, the BlockDiagonals package doesn’t cope too well with these increasing block sizes. The sparse matrix package does, but I have yet to find a way to convince it to invert the resulting block diagonal matrix.

My favorite solution would be to somehow get Julia to invert these sparse matrices that happen to be block diagonal, but if I cannot do it, I’ll have to implement my own BlockDiag class that better understands blocks. Does anyone know how to efficiently invert a block diagonal matrix encoded as a sparse matrix in Julia?

Thanks,

S

2

The first question is: do you really need to compute the full inverse, or do you just need to solve linear systems with it?

2 Likes
3

I’m doing something with H matrices. You typically have to compute inv(ABCDE)FGH or, if you prefer, (ABCDE)\FGH, or whatever decomposition. However, I’ve already tried that. If I do Z = factorize(ABCDE), I obtain a sparse LU decomposition, but then L\F is dense even if F is sparse.

So, on the one hand, you can make do with a factorization, but on the other hand, you need way more than a factorization that lets you solve with dense vectors. In the end, you need to compute inv(ABCDE)F for some sparse matrix F, whether by factorization or by matrix inverse is immaterial.

I wasn’t able to do this in any way in Julia. I mean, if you can do sparse\sparse then you can compute a sparse inverse by doing sparse\spdiagm(0=>ones(n))

4

Have you tried factorizing A, B, C, D, E independently, and then computing E \ (D \ (C \ (B \ (A\ FGH))))?

1 Like
5

When solving Ax=b with a sparse A, the output (along with the inverse A^{-1}) is generically dense even if b is sparse. (See also the thread Sparse solve with sparse rhs - #3 by stevengj) Julia’s generic sparse-matrix type assumes that this is the case and therefore always computes a dense output from A \ b or lu(A) \ b.

Of course, for the extremely special sparsity pattern of block-diagonal A this is not the case, and you indeed have a sparse (block-diagonal) inverse and sparse solutions for sparse right-hand-sides. But Julia’s SparseMatrixCSC will not exploit such “structured sparsity” — it is designed exclusively for unstructured sparsity.

Offhand, I can see at least three options for you:

  • Write your own block-diagonal arithmetic specialized for your case.
  • Submit a PR to BlockDiagonals.jl to support multiplying block-diagonal matrices with unequal — but compatible — block sizes, e.g. m \times m blocks with km \times km blocks. This is more work in some ways because the BlockDiagonals package probably supports more generic block sizes than you need, but may have the greatest long-term benefit.
  • Multiply the matrices as unstructured sparse matrices, but write a constructor that then converts the result to block-diagonal. (This is probably the least-efficient option.)
6

If L is the sparse lower triangular factor and F is sparse too, then L\F ends up being dense when I do it here, like I wrote above, so it doesn’t work.

7

Thanks Steven for your suggestions, I wrote this:

function sparseinv(A::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti<:Integer}
    a = [minimum(A.rowval[A.colptr[j]:A.colptr[j+1]-1]) for j=1:A.n]
    for j=A.n-1:-1:1
        a[j] = min(a[j],a[j+1])
    end
    b = [maximum(A.rowval[A.colptr[j]:A.colptr[j+1]-1]) for j=1:A.n]
    for j=2:A.n
        b[j] = max(b[j],b[j-1])
    end
    c = vcat(0,findall(b[1:end-1] .< a[2:end]),A.n)
    B = [sparse(inv(Matrix{Tv}(A[(c[j]+1):c[j+1],(c[j]+1):c[j+1]]))) for j=1:length(c)-1]
    blockdiag(B...)
end

However, I’m worried it won’t be enough. I suspect some of my matrices may end up being merely “morally” block diagonal, i.e. A[P,P] is block diagonal for some permutation P. I’m pretty sure MATLAB handles this fine too.

The problem is the algorithms implemented in the sparse libraries used by Julia (e.g. UMFPACK) don’t seem to have implemented the “natural” algorithm for sparse solves with a sparse right-hand-side, nor the “natural” algorithm for inverting a sparse matrix. You’re supposed to just do Gaussian elimination in the obvious manner and then be grateful that it stayed sparse. That’ll handle any encrypted block structure.

8

What is this “natural” algorithm? Generically, the inverse of a sparse matrix is dense.

Not quite — they spend a lot of effort permuting the rows/cols in an effort to minimize the fill-in of the L and U factors. (This is how nearly all unstructured sparse-direct libraries work, AFAIK.)

9

Perhaps this is not directly relevant, but I’ve tried to write a couple of BlockDiagonal helper functions:

using IterTools
using BlockDiagonals

isblockcompat(src, dst) = 
  accumulate((x,y)->x.+y, blocksizes(dst)) ⊆ accumulate((x,y)->x.+y, blocksizes(src))

function reblock(vals, shape)
    isblockcompat(vals, shape) || error("Block sizes incompatible")       
    return BlockDiagonal(collect(
      Iterators.map(((x,y),)->@view(vals[(:).(x.+1, y)...]), 
      IterTools.partition([(0,0);accumulate((x,y)->x.+y, blocksizes(shape))],2,1))
    ))
end

with these, BlockDiagonal matrices can be recast to a compatible BlockDiagonal shape:

julia> # prepare some matrices

julia> bm = BlockDiagonal([rand(1, 2), ones(2, 1), rand(2,2)])
5×5 BlockDiagonal{Float64, Matrix{Float64}}:
 0.658255  0.928167  0.0  0.0       0.0 
 0.0       0.0       1.0  0.0       0.0
 0.0       0.0       1.0  0.0       0.0
 0.0       0.0       0.0  0.904518  0.51962
 0.0       0.0       0.0  0.385899  0.994389

julia> cm = BlockDiagonal([rand(3,3), rand(2,2)])
5×5 BlockDiagonal{Float64, Matrix{Float64}}:
 0.439086  0.324064  0.614866  0.0       0.0
 0.855466  0.412557  0.82457   0.0       0.0
 0.770512  0.55656   0.362719  0.0       0.0
 0.0       0.0       0.0       0.968817  0.960156
 0.0       0.0       0.0       0.157498  0.390152

julia> # use functions

julia> isblockcompat(bm, cm)
true

julia> typeof(reblock(bm, cm))
BlockDiagonal{Float64, Matrix{Float64}}

julia> bm*cm
5×5 Matrix{Float64}:
 1.08305   0.596239  1.17008   0.0       0.0
 0.770512  0.55656   0.362719  0.0       0.0
 0.770512  0.55656   0.362719  0.0       0.0
 0.0       0.0       0.0       0.958152  1.07121
 0.0       0.0       0.0       0.530479  0.758486

julia> reblock(bm, cm)*cm
5×5 BlockDiagonal{Float64, Matrix{Float64}}:
 1.08305   0.596239  1.17008   0.0       0.0
 0.770512  0.55656   0.362719  0.0       0.0
 0.770512  0.55656   0.362719  0.0       0.0
 0.0       0.0       0.0       0.958152  1.07121
 0.0       0.0       0.0       0.530479  0.758486

The last expression is now a BlockDiagonal instead of a dense Matrix.

10

I guess I’m not entirely sure MATLAB has implemented this, but the algorithm I refer to works as follows. First, find the decomposition A = LU using your favorite algorithm. Then, it suffices to show how to compute C = L\B for a sparse matrix B. As everybody knows, in general, C will be dense, but in practical cases it may not be dense. We’ll thus do it sparsely.

Let b be some column of B, and c the corresponding column of C. If we denote by j \in J the entries b_j that are nonzero, and if we think of L as the adjacency matrix of some graph, then c_j \neq 0 is only possible if j is in the transitive closure (or you might say “connected component”) C(J) of J. So you have to run a little graph algorithm to compute C(J). Now all you have to do is compute those entries c_j with j \in C(J). If k = \#C(J) \ll n then you can gain a lot of performance, this is the case of my morally block diagonal matrices.

The easiest way to do this is to preallocate and pre-zero dense n-dimensional vectors x and y. You then copy sparse b into dense y vector, in the right places. You perform forward-substitution but only on the index set C(J) which requires something like k^2 FLOPS, this produces a dense representation x of the desired sparse vector c. You then copy the entries x_j into c_j for each j \in C(J).

You need to reset the vectors x and y to zero to proceed to the next columns of B and C but you can do this in O(k) time, because at this point the only entries y_j \neq 0 correspond to j \in J, and the only entries x_j \neq 0 correspond to j \in C(J).

I think this algorithm may be in the book of Tim Davis, and indeed it may be in SuiteSparse here. I think SuiteSparse used to be in Julia so perhaps Julia used to have this capability and lost it?

11

I don’t think Julia ever exported an interface using this function. There was a discussion in 2013 about possibly exporting it, but it never went anywhere: Inverse of SparseMatrixCSC · Issue #4439 · JuliaLang/julia · GitHub … it hasn’t come up very often, I think, because it’s so rare to use unstructured sparse matrix formats in the unusual cases where the inverse is sparse — structured formats like BlockDiagonal are usually more efficient in these cases. But you’re right that it could be convenient in cases like yours.

It should be straightforward to expose it, at least for the matrix types supported in UMFPACK, via careful use of ccall.

12

I think I’ll bite the bullet and either ccall it, as you write, or otherwise implement something. Looking at my current project, and from previous experiences elsewhere, when I write complicated codes to handle whatever special case I have, instead of relying on robust matrix primitives, I always regret it. It’s also bitten me when I expected the underlying primitive to behave a certain way. A while back, I was writing a preconditioner for PetSc and I was using the various matrix primitives built into PetSc. I think I multiplied some sparse matrices together. I knew that theoretically this could be done efficiently, and I could have written special-purpose code to make sure it was done efficiently, but I naively assumed PetSc’s primitive would be clever enough and unfortunately it didn’t work out. That preconditioner ended up being slow.

13

Thanks for this, you are right, this is also one of the possible solutions I’m currently considering.

Cheers,

S

14

Side note: If you have a preconditioner ABC that is a product of sparse matrices, normally I would express this as a linear operator x \mapsto A(B(Cx)) to avoid the risk of losing sparsity in the product. (I think in Petsc you can do this as a “composite” preconditioner, or by defining a custom “shell” preconditioner with your own matvec. In Julia, you can implement such “lazy” products easily with LinearMaps.jl or LinearOperators.jl.)