JuliaArrays / BlockDiagonals.jl

Functionality for working efficiently with block diagonal matrices.
MIT License
49 stars 11 forks source link

`BlockDiagonal` name is misleading #94

Closed mzgubic closed 2 years ago

mzgubic commented 2 years ago

This package allows a BlockDiagonal matrix to have non-square blocks on the diagonal, which is in contrast to Block Diagonal Matrix, a term reserved for the subset of BlockDiagonals with square blocks.

Note that this package has nothing to do with Block Matrices, which are implemented in BlockArrays.jl.

This identity confusion has led to bugs, e.g. the diag implementation which leads to https://github.com/invenia/BlockDiagonals.jl/issues/67 assumes that the blocks are square, and similarly for svdvals

So we need a name for:

  1. the matrix where the blocks are square
  2. the matrix where the blocks can be rectangular
  3. an abstract supertype for the two

Currently, we have:

  1. no name
  2. BlockDiagonal
  3. no name

Math convention is:

  1. Block Diagonal
  2. no name
  3. no name

We could have instead:

  1. BlockDiagonal
  2. RectangularBlockDiagonal
  3. AbstractBlockDiagonal

This is a breaking change since it introduces an extra constraint on the existing type (BlockDiagonal).

nickrobinson251 commented 2 years ago

This identity confusion has led to bugs

That's not really how these bugs came about; we used to support only square blocks and when we dropped that restriction we just missed a few functions that needed updating.

Supporting non-square blocks was discussed at length in https://github.com/invenia/BlockDiagonals.jl/pull/22, and we conclude that there is a use-case for it and no objections to supporting it. As I say, we missed a few function that needed updating alongside that change (really we should have tested everything on non-square blocks as part of the change), so it's good to get those bugs fixed.

Block Diagonal Matrix, a term reserved for the subset of BlockDiagonals with square blocks.

This package's BlockDiagonal has the same behaviour as Scipy's block_diag, so there is precedent for software supporting non-square blocks (and under the "block diagonal" name), and they note in the docstring that "Block Diagonal Matrix" usually refeserve for square blocks -- do you think improved documentation would be sufficient for this package too?

has nothing to do with Block Matrices, which are implemented in BlockArrays.jl.

This is false (or i'm somehow misunderstanding this sentence).

This package concerns a subset of block matrices; those that have the form image where each A_i is a matrix (not necessarily square, as noted above), and each 0 is a zero-matrix.

This class of matrix can also be represented using BlockArrays.jl. I think the BlockArrays.jl spelling is mortar(Diagonal([A1, A2])) whereas for BlockDiaonals.jl it'd be BlockDiagonal([A1, A2]). Again, A1, A2 need not be square for BlockArrays.jl either. BlockDiaonals.jl and BlockArrays.jl both only stores the A_i matrices (not the zero matrices).

We could have instead: ... introduces an extra constraint on the existing type

Would these two concrete types actually store the data differently? If not, what would be the benefit of having the different types, if the only differences is an assertion in the constructor?

Would there be a benefit to being able to dispatch on "all the blocks are square", e.g. do we have lots of code like foo(b::BlockDiagonal) = all(is_square, blocks(b)) ? foo_square(b) : foo_nonsquare(b)?

Overall, i could be convinced of introducing some kind of distinction between the square-block case and the non-square-block case, but I think there'd have to be a compelling practical benefit

mzgubic commented 2 years ago

You've convinced me that changing the public API is unnecessary, as long as we change the docs and add tests.

As I say, we missed a few function that needed updating alongside that change (really we should have tested everything on non-square blocks as part of the change)

Ah, I see, have opened an issue: https://github.com/invenia/BlockDiagonals.jl/issues/97

This package's BlockDiagonal has the same behaviour as Scipy's block_diag, ... do you think improved documentation would be sufficient for this package too?

I imagine most people will still expect square blocks but if we have clear docstrings then I think that's sufficient given the scipy precedent.

has nothing to do with Block Matrices, which are implemented in BlockArrays.jl.

Ah, sorry, poor wording. I meant to say that this package does not want to support all Block Matrices, which in general have non zero off-diagonal entries. I did not realise though that BlockArrays.jl supports sparse zero offidiagonal blocks.

Overall, i could be convinced of introducing some kind of distinction between the square-block case and the non-square-block case, but I think there'd have to be a compelling practical benefit

I've thought about this over the weekend a bit. Indeed the data storage would be done in the same way, and one practical benefit is that we could simplify dispatch on more efficient implementations e.g. for svdvals. I don't think its too common but it does feel wrong even if it is for a small number of cases. Perhaps adding a parameter to the type could be the least invasive way of adding this distinction?

mzgubic commented 2 years ago

closing as cleared up, the type parameter is the only unresolved bit, and can be discussed separately in https://github.com/invenia/BlockDiagonals.jl/issues/107