[BUG] Some SVD truncations don't work

[mtfishman]

For example:

using BlockArrays: Block, blockedrange
using BlockSparseArrays: BlockSparseArray
using MatrixAlgebraKit: svd_trunc
r = blockedrange([2, 2])
a = BlockSparseArray{Float64}(undef, r, r)
a[Block(2, 2)] = randn(2, 2)
u, s, v = svd_trunc(a; trunc=(; atol=1e-3))
s

outputs something like:

2×2-blocked 4×4 BlockSparseMatrix{Float64, Matrix{Float64}, …, …}:
  ⋅    ⋅   │   ⋅        ⋅      
  ⋅    ⋅   │   ⋅        ⋅      
 ──────────┼───────────────────
  ⋅    ⋅   │  1.64506  0.0     
  ⋅    ⋅   │  0.0      0.359026

If Block(1, 1) is set instead, the truncation appears to work (i.e. it drops the unallocated block). Also, maxrank appears to work.

[mtfishman]

It looks like the root cause of this is that the MatrixAlgebraKit.TruncationKeepAbove truncation strategy assumes the spectrum is sorted, which isn't the case above:

julia> using MatrixAlgebraKit

julia> values = [0, 0, 0.5, 0.3]
4-element Vector{Float64}:
 0.0
 0.0
 0.5
 0.3

julia> strategy = MatrixAlgebraKit.TruncationKeepAbove(1e-3, 0)
MatrixAlgebraKit.TruncationKeepAbove{Float64}(0.001, 0.0)

julia> MatrixAlgebraKit.findtruncated(values, strategy)
1:4

In contrast, setting maxrank works because it doesn't assume the spectrum is sorted:

julia> strategy = MatrixAlgebraKit.TruncationKeepSorted(2, abs, true)
TruncationKeepSorted{typeof(abs)}(2, abs, true)

julia> MatrixAlgebraKit.findtruncated(values, strategy)
2-element Vector{Int64}:
 3
 4

[mtfishman]

Defining truncation as the following:

function MatrixAlgebraKit.findtruncated(values::AbstractVector, strategy::TruncationKeepAbove)
  atol = max(strategy.atol, strategy.rtol * first(values))
  if issorted(values; by=abs, rev=true)
    howmany = @something findlast(≥(atol), values) 0
    return 1:howmany
  else
    perm = sortperm(values; by=abs, rev=true)
    # `collect` fixes a bug when `values isa BlockSparseVector` with a graded axis,
    # we'll have to investigate that issue.
    howmany = @something findlast(≥(atol), collect(values)[perm]) 0
    return perm[1:howmany]
  end
end

fixes this issue. @lkdvos does something like that seem reasonable for MatrixAlgebraKit.jl?

[mtfishman]

I made a PR along those lines here: QuantumKitHub/MatrixAlgebraKit.jl#26.

[lkdvos]

I don't think this is great for generic vectors since it's type unstable always, but it seems reasonable to specialize this for types where the output isn't sorted, and then don't even do the check?

What I'm also thinking is that it might be reasonable to sort the values here before handing them off to matrixalgebrakit, since it is really fair to assume that singular values are sorted

[mtfishman]

Sorting before inputting into MatrixAlgebraKit.findtruncated could be an option, I was hoping we could use as much generic code from MatrixAlgebraKit.jl as possible but it seems like that would require overloading things at a higher level, say even overloading MatrixAlgebraKit.svd_trunc!.

Alternatively, if MatrixAlgebraKit.findtruncated accepted more arguments about the factorization, such as the matrix being truncated, the algorithm, etc. I could see that might give us enough context to specialize it on the fact that we are doing a blockwise block sparse SVD.

I would find it finicky to overload MatrixAlgebraKit.findtruncated(values::BlockSparseVector, strategy::TruncationKeepAbove) since we wouldn't know where that spectrum came from (it doesn't even need to be an SVD, it could have been some other factorization).

[mtfishman]

Fixed by #26.