Reivive this package?

[mipals]

Hi there,

I tried to use this package for some computations. While it works it does look to be pretty inefficient at times. E.g. multiplication allocates lots of smaller vectors, linear time indexing ( #121 ), and probably more.

The recent pull requests also seem stuck. As such I am in doubt whether to contribute here or simply create my own code. Any input would be greatly appreciated.

Cheers, Mikkel

[dkarrasch]

Another option you may wish to consider is to use LinearMaps.jl, which has limited "AbstractMatrix" support, though.

[mipals]

@dkarrasch I've used LinearMaps.jl quite a lot in the past with success (great package! 🥳). Do LinearMaps.jl have the option of providing an exact inverse map? This is quite a big reason for me to use the BlockDiagonal structure.

In any case I still think that this package could use some care 🙂

[dkarrasch]

Inverses and solves are not within the scope of LinearMaps.jl, but we do have an InverseMap type, whose multiplication corresponds to solves. Feel free to open an issue there to continue the discussion.

[mipals]

I ended up creating a small package instead. Two reasons: Nobody seemed to respond here (other than you @dkarrasch - Thanks a lot for your input), and secondly that I do not agree with the choice made in this package of not storing the block indicies (this makes every computation serial by design and indexing of large matrices slow as you always have to start at index (1,1) and then loop through the blocks).

[jishnub]

I think it'll be good to revive this package. I'll try to see if I can contribute

[mipals]

Okay, that kind of put me in a weird spot. As I said earlier I would be happy to also help out, but now I've created my own package. Anyhow, maybe we should discuss how to combine the efforts?

[nathanaelbosch]

I would also love to use this package (and I also ended up re-implementing the functionality I need recently here). There are some clear things I could contribute, but PRs do not seem to get merged right now. So I have the same question as OP: Should we revive this package?

[jishnub]

Yes please. I can try to push PRs, if someone is willing to contribute

[mipals]

We could revive this package, but before doing so it might be a good idea to briefly discuss what our issues with the current state of the package is. If these align, or at the very least are not opposing eachother, it would make sense to continue.

Personally I am interested in very large matrices with relatively small blocks. Such a matrix could be represented fairly efficiently using a sparse matrix. However if each small block has structure it would be best to utilize said structure in the computation - As such I needed a BlockDiagonal package.

What annoyed me with this package is that it does not include global indices. As such indexing is extremely slow for large matrices (note the plot here only show for small matrices - The problem keeps growing). At the same it result in most linear algebra being sequential by design and you therefore loose out on otherwise embarrassingly parallel computation. For me the bare minimum would be to extend the structure to behave something like this (but not necessarily exactly like this)

struct BlockDiagonal{T, V <: AbstractMatrix{T}} <: AbstractBlockDiagonal{T}
    blocks::Vector{V}
    block_row_indices::Vector{Int}
    block_col_indices::Vector{Int}
    is_block_square::Bool
end

where the rows and cols are separated in order include non-square blocks - and then update everything to utilize this additional information in the computations to fix the run-time of indexing, size and linear algebra routines.

[nathanaelbosch]

Personally I just want operations to be broadcast to operate on the blocks. I only need indexing for printing in the REPL and nothing else (I would even be fine with indexing throwing an error just to make sure that the matrix is never actually built; but that's not a serious proposal for this package of course). So for my purposes,

struct BlockDiagonal{T, V <: AbstractMatrix{T}} <: AbstractBlockDiagonal{T}
    blocks::Vector{V}
end

is fully sufficient.

I think the most important thing might be to actually define an AbstractBlockDiagonal type with a basic interface, and then define most of the current operations for this abstract type. Then if different people have different particular requirements, they can implement their own subtypes to do the one thing that's missing, but without re-implementing every Base function and every basic linear algebra overload. Would that be a reasonable direction forward?

[cvsvensson]

One thing that distinguishes this package from block diagonals in LinearMaps.jl is that it subtypes AbstractMatrix. As such, an efficient getindex would be very useful. This issue has also been discussed for LinearMaps in https://github.com/JuliaLinearAlgebra/LinearMaps.jl/issues/38 and https://github.com/JuliaLinearAlgebra/LinearMaps.jl/issues/180. The julia manual also sets up this expectation for AbstractArrays. Quoting from the manual:

It is recommended that these operations have nearly constant time complexity, as otherwise some array functions may be unexpectedly slow.

I wasn't aware of #121 and based on the linked discussions I didn't expect it. I like @mipals suggestion, but I don't know how development should be organized.

[mipals]

I now just realize that LinearMaps.jl actually have block-diagonal maps. I thought the earlier suggestion was to use the regular LinearMap 😅 Anyhow, as stated by @cvsvensson there are some difference in that this package subtypes an AbstractMatrix (and the blocks do as well!).

I already have code that fixes the indexing (It runs log(n) but I guess that is as close you we can get to constant in this case). Having the block information can also be used to make other operations efficient (i.e. size can be done in constant time and not linear w.r.t. the number of blocks). I believe that including the code will be non-breaking as long as the blocks of the new structure can be accessed as .block ?