[Enhancement] Add singular value decomposition

[lkdvos]

This PR is a migration of ITensor/ITensors.jl#1572, which aims to add support for SVDs.

[codecov]

Codecov Report

Attention: Patch coverage is 22.41379% with 90 lines in your changes missing coverage. Please review.

Project coverage is 27.25%. Comparing base (7570c2c) to head (af4a595).
Report is 1 commits behind head on main.

Files with missing lines	Patch %	Lines
...tractblocksparsearray/abstractblocksparsematrix.jl	0.00%	53 Missing ⚠️
src/factorizations/svd.jl	60.46%	17 Missing ⚠️
src/BlockArraysExtensions/BlockArraysExtensions.jl	0.00%	10 Missing ⚠️
src/blocksparsearray/blockdiagonalarray.jl	0.00%	10 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #16      +/-   ##
==========================================
- Coverage   27.48%   27.25%   -0.23%     
==========================================
  Files          26       29       +3     
  Lines         866      987     +121     
==========================================
+ Hits          238      269      +31     
- Misses        628      718      +90     

Flag	Coverage Δ
docs	26.80% <22.41%> (-0.59%)	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[mtfishman]

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?

[mtfishman]

Related to my question about how the code is handling rectangular block structures, I also noticed it makes use of the full keyword argument, but I can't tell if that is handling rectangular block structures properly, i.e. as we discussed, if the block structure is rectangular and full=true it should probably construct some new random unitary blocks on U or V depending on the block structure.

That also reminds me that we will need to think about how to preserve sector information when the axes are graded but we don't necessarily have to do that in this PR.

[lkdvos]

So, as an intermediate update:

• I overhauled the svd implementation quite a bit, now no longer taking the BlockDiagonal path. I was trying to make that work, but the combination of rectangular blocks, as well as empty rows and or columns actually made me realize that I was trying too hard to make that work, and the code simplified quite a bit by letting that go. In summary, the current implementation is now a lot closer again to the original implementation that @mtfishman made.
• I'm pretty sure that the svd implementation is now correct, and that it should handle rectangularity as expected. There is an ambiguity in what kind of axes S has whenever there are multiple empty rows/columns, since we allow for a permutation that could arbitrarily make this "blockdiagonal", which is now not uniquely determinable. In this case, I simply keep the current order of the blocks as much as possible. (note that this is different from graded axes, where there is a more natural way of "blockdiagonalizing" based on the flux.
• The tests are completely failing since the resulting objects are not behaving. In particular, it is not clear how to recombine: with A = U * Diagonal(S) * Vt, S is an annoying Diagonal-type to work with (that we do not own), and other options break the interface of LinearAlgebra slightly. There are also some issues with missing Derive implementations

Given this, I would suggest to add a diagonal(::AbstractVector)::AbstractMatrix function in DiagonalArrays.jl to act similar to adjoint or transpose, in the sense that it should create a diagonal array but is not opinionated about the output type. This is opposed to Diagonal, which implicitly expect the return type to be <:Diagonal.
The default implementation there probably should just return Diagonal, but for diagonal(::BlockSparseVector) we can then decide on what type is most fitting.

The easiest here is probably BlockSparseMatrix with blocktype DiagonalArray of a sparse vector of DiagonalArray, since we own all these types and can make them function as desired.

Before I move on with this, I'd like to ask some feedback to see if I'm not overlooking anything, as well as seeing if anyone has better ideas, which I would gladly listen to 😃

[mtfishman]

@lkdvos this is getting a bit hard to review but I think it is good to go from my end assuming open issues we discussed are addressed in follow-up PRs.

[lkdvos]

Yeah, I am having the same problem keeping track of everything. Let me double check things are behaving as they should, add some comments about todos and open some issues, I'll update the draft status here once I have that figured out.

[lkdvos]

This depends on having ITensor/DerivableInterfaces.jl#11 merged, but afterwards I think I managed to resolve all issues.

To summarize:

In order to attempt to keep an overview, I chose to (temporarily) not attempt to make Diagonal(::Diagonal) work, since this is really a separate topic. The svd implementation here is working, and yields the correct results, and the discussion on the diagonal implementation is mostly a matter of how to recombine the tensors.

A final change I made is a result of testing more thoroughly the case where a row/column is fully empty. Then, it is necessary to add identity blocks in U and V to make them isometries again.

[mtfishman]

@lkdvos besides the comment I left about reorganizing the code logic around handling missing blocks, this looks good to me.