[BUG] SVD of rectangular block structure

[mtfishman]

using BlockArrays
using BlockSparseArrays
using BlockSparseArrays: blocksparse
using MatrixAlgebraKit

ax = (blockedrange([2, 2]), blockedrange([2, 2, 2]))
bs = Dict(
  Block(1, 1) => randn(2, 2),
  Block(2, 3) => randn(2, 2),
)
a = blocksparse(bs, ax)
u, s, v = svd_compact(a)

returns:

ERROR: LoadError: BoundsError: attempt to access 2-element Vector{Base.OneTo{Int64}} at index [3]
Stacktrace:
  [1] throw_boundserror(A::Vector{Base.OneTo{Int64}}, I::Tuple{Int64})
    @ Base ./essentials.jl:14
  [2] setindex!
    @ ./array.jl:986 [inlined]
  [3] initialize_output(::typeof(svd_compact!), A::BlockSparseMatrix{…}, alg::BlockSparseArrays.BlockPermutedDiagonalAlgorithm{…})
    @ BlockSparseArrays ~/.julia/dev/BlockSparseArrays/src/factorizations/svd.jl:65
  [4] svd_compact!
    @ ~/.julia/packages/MatrixAlgebraKit/Yrs13/src/algorithms.jl:189 [inlined]
  [5] #svd_compact!#67
    @ ~/.julia/packages/MatrixAlgebraKit/Yrs13/src/algorithms.jl:183 [inlined]
  [6] svd_compact!
    @ ~/.julia/packages/MatrixAlgebraKit/Yrs13/src/algorithms.jl:182 [inlined]
  [7] svd_compact(A::BlockSparseMatrix{Float64, Matrix{…}, SparseArraysBase.SparseMatrixDOK{…}, Tuple{…}})
    @ MatrixAlgebraKit ~/.julia/packages/MatrixAlgebraKit/Yrs13/src/algorithms.jl:178
  [8] top-level scope
    @ ~/.julia/dev/BlockSparseArrays/dev/rectangular_svd_bug.jl:12
  [9] include(fname::String)
    @ Main ./sysimg.jl:38
 [10] top-level scope
    @ REPL[8]:1
 [11] top-level scope
    @ none:1
in expression starting at /Users/mfishman/.julia/dev/BlockSparseArrays/dev/rectangular_svd_bug.jl:12
Some type information was truncated. Use `show(err)` to see complete types.

The error is thrown at this line:

BlockSparseArrays.jl/src/factorizations/svd.jl

Line 65 in 83ed436

u_axes[col] = infimum(brows[row], bcols[col])

. We have to think more carefully about which blocks should get set when the block structure is rectangular.

[mtfishman]

@lkdvos in general I've found the current block sparse SVD code to be a bit tough to reason about when trying to handle various corner cases like this one and also handle symmetries, maybe we should reconsider rewriting the logic in terms of first transforming into a block diagonal matrix with a block permutation and then taking the SVD of that matrix.

[lkdvos]

I'll have a look, but perhaps that means that it is now really time to just sit down and actually design a proper type for the BlockSparseDiagonal (or any other name) and go from there.

[mtfishman]

I'll have a look, but perhaps that means that it is now really time to just sit down and actually design a proper type for the BlockSparseDiagonal (or any other name) and go from there.

I think that would make sense. Though I would hope a lot of code would be agnostic about the block storage anyway.