Make compatible with BlockBandedMatrices.jl

[dlfivefifty]

Note https://github.com/JuliaMatrices/BlockBandedMatrices.jl exists. We should make sure these packages are compatible in terms of the "block-banded" interface: e.g., blockbandwidths(::BlockDIagonal) should return (0,0).

Also, there is a BlockArray with Diagonal backend which should also be supported:

https://github.com/JuliaArrays/BlockArrays.jl/blob/78fcde4b63a607ae6e63d76e85f4685e33637ea9/test/test_blockarrayinterface.jl#L79

[nickrobinson251]

Thanks for opening this! I was thinking that it'd be nice for this package to work nicely within the BlockArrays/Matrices ecosystem, but it's not fully decided, and a little extra guidance would be appreciated!

compatible in terms of the "block-banded" interface

I see that blockbandwidths(::BlockDiagonal) currently gives the wrong answer... but I'm not sure a dependency on BlockBandedMatrices.jl is worth it for extending a single method. What else would be gained and are there other methods that'd need implementing? (Sorry, i've not used BlockBandedMatrices myself)

there is a BlockArray with Diagonal backend which should also be supported

What did you have in mind specifically?

I wonder if instead of having BlockDiaognal as part of the BlockArray ecosystem, we should instead just improve LinearAlgebra's Diagonal{<:AbstractMatrix} (related JuliaLang/LinearAlgebra.jl#610). BlockDiagonals.jl would be a useful package until LinerAlgebra is improved, then it'd become redundant, and the BlockArrays ecosystem would still have BlockMatrix{T, Diagonal{<:AbstractMatrix{T})} as it does now.

Would welcome feedback/comments/suggestions and of course PRs!

pinging @ararslan who has contributed here and may or may not have opinions :)

[dlfivefifty]

I see that blockbandwidths(::BlockDiagonal) currently gives the wrong answer... but I'm not sure a dependency on BlockBandedMatrices.jl is worth it for extending a single method.

Definitely not. We can move blockbandwidths to BlockArrays.jl

What did you have in mind specifically?

I haven't thought it through much, but I suppose making

const BlockDiagonal{T,MT<:AbstractMatrix{T}} = BlockMatrix{T,Diagonal{MT},DefaultBlockSizes{2}}

might be a good way.

Note DefaultBlockSizes should be replaced at some point with a special block sizes type that captures that the row and column blocksizes are the same.

[nickrobinson251]

Hmm, interesting, there's quite a significant conceptual difference between BlockDiagonal{T, ::AbstractMatrix{T}} and BlockMatrix{T, ::Diagonal{<:AbstractMatrix{T}}} right now:

The BlockDiagonal assumes the underlying structure is a AbstractMatrix{<:Number} and "blocks" are useful ways to refer to some of those Numbers. e.g. diag returns a Vector of Numbers

The "BlockMatrix{Diagonal}" assumes the underlying structure is AbstractMatrix{<:AbstractMatrix} and "blocks" are useful ways to refer to those "inner" Matrices e.g. diag returns a Vector of Matrices

[edited for clarity]

[dlfivefifty]

That is not correct;

julia> A = BlockArray(randn(6,6), 1:3,1:3)
3×3-blocked 6×6 BlockArray{Float64,2}:
 -0.84975   │   2.20661   -0.920463  │   1.9042      -2.02121    -0.0278585 
 ───────────┼────────────────────────┼──────────────────────────────────────
 -1.38429   │  -0.410338  -0.311795  │   1.47326     -0.36143     1.64884   
  0.831705  │  -0.658984  -1.27877   │   1.48502      0.180917   -0.270828  
 ───────────┼────────────────────────┼──────────────────────────────────────
 -1.05782   │  -0.468753  -1.10017   │  -0.00996191  -0.0705051  -0.146721  
 -1.9989    │   0.856402   1.71011   │   0.143385     1.43957    -0.00442993
  0.266811  │  -0.174526   0.142927  │  -0.970129    -0.9923      1.1756    

julia> diag(A)
6-element Array{Float64,1}:
 -0.849750393735063   
 -0.41033810902489665 
 -1.2787711809201605  
 -0.009961911521250116
  1.4395683893644051  
  1.1756042915582254  

[nickrobinson251]

ah, sorry, i didn't mean that it was true of BlockArrays in general, only of the "Diagonal BlockArray" BlockMatrix{T, Diagonal{Matrix{T}}}

julia> blocks = [rand(2, 2), rand(3, 3)];

julia> BD = BlockDiagonal(blocks);

julia> diag(BD)
5-element Array{Float64,1}:
 0.4462094596368906
 0.6789518568486066
 0.8103392123794502
 0.3730887997881511
 0.25490783919204985

julia> BA = BlockArray(Diagonal(blocks));

julia> diag(BA)
2-element Array{Array{Float64,2},1}:
 [0.446209 0.173019; 0.720358 0.678952]
 [0.810339 0.895358 0.722163; 0.439333 0.373089 0.271466; 0.794175 0.673562 0.254908]

[dlfivefifty]

That is a bug.

[dlfivefifty]

Actually, it’s not a bug: you want mortar(Diagonal(blocks)). BlockArray(A::AbstractArray) == A should hold true.

[nickrobinson251]

Thanks a lot for your help! I'm going to play around with this a bit more and try figure out what a better integration with BlockArrays might be :D

[ararslan]

It seems to me that this package's functionality could be fully folded into BlockArrays, though I recall some hesitance from @eperim about that. After all, a block diagonal matrix is a straightforward special case of a block matrix.

[dlfivefifty]

I would be fine with that. In fact, I think all Block* where * exists in LinearAlgebra.jl (BlockDiagonal, BlockTridiagonal, BlockSymTridiagonal, BlockUpperTriangular, etc.) can live in BlockArrays.jl

[eperim]

My hesitation was mainly in adding extra dependencies to this package, but that has already happened, so, as long as we don't lose any of the functionalities here, I have nothing against it.

[nickrobinson251]

I've not had time to work on this for a while, but at some point i would like to both...

1. Contribute this into BlockArrays
2. Make this alternative package that is separate for BlockArrays (does not depend on it / implement the interface) in case anyone ever wanted a light weight package for block diagonals

of course, i'd be delighted to see someone else do this first! I hope to have time later this summer :)

[nickrobinson251]

I think I will close this and open an issue on BlockArrays.jl for "Add more "Block Diagonal" functionality to BlockArrays.jl" But would be happy to re-open if anyone thinks other things should change in this repo. I found this discussion really help :) thanks!