[Enhancement] Add singular value decomposition #16
Conversation
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #16 +/- ##
==========================================
- Coverage 27.30% 24.65% -2.65%
==========================================
Files 26 29 +3
Lines 857 949 +92
==========================================
Hits 234 234
- Misses 623 715 +92
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Sentry. |
|
|
||
| # SVD implementation | ||
| function eigencopy_oftype(A::BlockDiagonal, S) | ||
| diag = map(Base.Fix2(eigencopy_oftype, S), A.blocks.diag) |
There was a problem hiding this comment.
We should make a proper function/API to access A.blocks.diag.
Some ideas are:
eachstoredblockDiagonalArrays.diagview(blocks(a))(from https://github.com/ITensor/DiagonalArrays.jl)
There was a problem hiding this comment.
I think ideally we would have eachstoredblock special cased to output the same thing as diagview(blocks(a)) in the case of block diagonal arrays, though that may be a more involved thing beyond this PR. It may just require defining storedvalues(a::AbstractDiagonalArray/Diagonal) = diagview(a), and then eachstoredblock(a) is already defined as storedvalues(blocks(a)) so then eachstoredblock(a) would output diagview(blocks(a)).
| if isnothing(A′) | ||
| # not block-diagonal, fall back to dense case | ||
| Adense = eigencopy_oftype(A, T) | ||
| return svd!(Adense; full, alg) |
There was a problem hiding this comment.
Does this output an SVD object?
There was a problem hiding this comment.
(I see from the definition below that it does.)
| # SVD is implemented by trying to | ||
| # 1. Attempt to find a block-diagonal implementation by permuting | ||
| # 2. Fallback to AbstractBlockArray implementation via BlockedArray | ||
| function svd( |
There was a problem hiding this comment.
So I see this is designed around a private BlockSparseArrays.svd function, is the plan to overload LinearAlgebra.svd? Probably we should define this as:
@interface ::AbstractBlockSparseArrayInterface function svd(...)(I would like to make the AbstractArrayInterface interface types take an ndims parameter so it could be AbstractBlockSparseMatrixInterface but that isn't supported yet) and then derive that function for LinearAlgebra.svd(::AnyAbstractBlockSparseMatrix) in the style of Derive.jl.
| A::AbstractBlockSparseMatrix; full::Bool=false, alg::Algorithm=default_svd_alg(A) | ||
| ) | ||
| T = LinearAlgebra.eigtype(eltype(A)) | ||
| A′ = try_to_blockdiagonal(A) |
There was a problem hiding this comment.
I found it confusing that this is call A′ but if it is a blocked permutation it is actually a composite object of the block diagonal matrix and the block permutations, we should make that clearer in the variable name, i.e. A′_and_blockperms.
| # compute block-by-block and permute back | ||
| A″, (I, J) = A′ | ||
| F = svd!(eigencopy_oftype(A″, T); full, alg) | ||
| return SVD(F.U[Block.(I), Block.(J)], F.S, F.Vt) |
There was a problem hiding this comment.
Can I and J already be designed as lists of Block?
There was a problem hiding this comment.
Actually it looks like that is already the case so it seems like Block.(...) may not be necessary.
There was a problem hiding this comment.
I'm a bit confused by this code, I would have though it would be something like:
SVD(F.U[I, :], F.S, F.Vt[:, J])i.e. I and J are block permutations of the outer spaces but maybe I have the wrong mental model for how this function is being designed.
|
Overall this code looks really nice and simple, thanks @lkdvos. I think this is a nice payoff for the BlockSparseArrays.jl design. I have more of a conceptual question, how does this code handle rectangular block structures? |
This PR is a migration of ITensor/ITensors.jl#1572, which aims to add support for SVDs.