Separate Sparse Matrix Blocks into Matrices

General Usage
1

Is there an efficient way to convert a block diagonalizable sparse matrix into the separate dense/near dense sub-blocks?

As a small example, I have the following matrix:

julia> mat
5×5 SparseMatrixCSC{Float64, Int64} with 13 stored entries:
 14.5   ⋅         ⋅         0.3    ⋅ 
   ⋅   8.5       0.173205    ⋅    0.3
   ⋅   0.173205  7.5         ⋅    0.866025
  0.3   ⋅         ⋅        10.0    ⋅ 
   ⋅   0.3       0.866025    ⋅   14.5

And would like to automatically generate the permutations to produce the following sub-matrices

julia> sub1 = mat[[1;4],[1;4]]
2×2 SparseMatrixCSC{Float64, Int64} with 4 stored entries:
 14.5   0.3
  0.3  10.0

julia> sub2 = mat[[2;3;5],[2;3;5]]
3×3 SparseMatrixCSC{Float64, Int64} with 9 stored entries:
 8.5       0.173205   0.3
 0.173205  7.5        0.866025
 0.3       0.866025  14.5

Thank you in advance and apologies if this is similar to another post

2

Maybe related to the bandwidth minimization problem (Graph bandwidth - Wikipedia, Cuthill–McKee algorithm - Wikipedia, GitHub - PetrKryslUCSD/SymRCM.jl: Reverse Cuthill-McKee node-renumbering algorithm.). Though you may need additional permutations to identify larger blocks.

I’d be curious to hear other ideas, we are interested in this problem as well.

3

It looks like SymRCM.jl does the job, at least in this case:

julia> using SparseArrays

julia> using SymRCM

julia> Is = [(1, 1), (2, 2), (3, 3), (4, 4), (5, 5)]
5-element Vector{Tuple{Int64, Int64}}:
 (1, 1)
 (2, 2)
 (3, 3)
 (4, 4)
 (5, 5)

julia> Js = [(1, 4), (2, 3), (2, 5), (3, 5)]
4-element Vector{Tuple{Int64, Int64}}:
 (1, 4)
 (2, 3)
 (2, 5)
 (3, 5)

julia> a = spzeros(5, 5)
5×5 SparseMatrixCSC{Float64, Int64} with 0 stored entries:
  ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅ 

julia> foreach(I -> (a[I...] = 10I[1] + I[2]), Is)

julia> foreach(I -> (a[I...] = 10I[1] + I[2]), [Js; reverse.(Js)])

julia> a
5×5 SparseMatrixCSC{Float64, Int64} with 13 stored entries:
 11.0    ⋅     ⋅   14.0    ⋅ 
   ⋅   22.0  23.0    ⋅   25.0
   ⋅   32.0  33.0    ⋅   35.0
 41.0    ⋅     ⋅   44.0    ⋅ 
   ⋅   52.0  53.0    ⋅   55.0

julia> p = symrcm(a)
5-element Vector{Int64}:
 5
 3
 2
 4
 1

julia> a[p, p]
5×5 SparseMatrixCSC{Float64, Int64} with 13 stored entries:
 55.0  53.0  52.0    ⋅     ⋅ 
 35.0  33.0  32.0    ⋅     ⋅ 
 25.0  23.0  22.0    ⋅     ⋅ 
   ⋅     ⋅     ⋅   44.0  41.0
   ⋅     ⋅     ⋅   14.0  11.0
4

This looks great! Thank you!

Do you happen to know if there is a good way to split the permutation along the blocks of the matrix? My ultimate goal is to extra the blocks to diagonalize separately rather than diagonalizing the larger matrix

5

One way to do that could be by identifying the connected components of the graph that is derived by interpreting the matrix as an adjacency matrix:

julia> using Graphs

julia> using LinearAlgebra

julia> connected_components(Graph(Symmetric(a[p, p])))
2-element Vector{Vector{Int64}}:
 [1, 2, 3]
 [4, 5]

Note that appears to be another way to find the blocks, without needing bandwidth minimization to begin with:

julia> connected_components(Graph(Symmetric(a)))
2-element Vector{Vector{Int64}}:
 [1, 4]
 [2, 3, 5]

Again, I’d be curious to hear if anyone has other suggestions, there may be a more efficient/robust way to do this.