ogauthe
commented
2 days ago Looks nice, I'm eager to try.
Considering truncation, indeed the simplest is to compute dense block SVD and truncate afterwards. In the case of non-abelian symmetries, there are choices to be made on how to define the ranks of each block.
When non-abelian symmetries are involved, one does not control $D$ directly but $D^*$, the number of kept multiplets. One can either directly set $D^*$ as a target or else set a $D$ target and keep adding mutliplets until $D$ is reached. This choice matters because when simulating a given model, depending on the parameters such as the ratio of 2 coupling constants $r$, the multiplet decomposition of a singular spectrum may dramatically change.
If one sets a $D^*$ as a target, any truncation is straightforward. The tradeoff is that $D$ is not controlled and may reach very different values when tuning $r$. Moreover, it makes it hard to benchmark simulations done with different symmetries (for instance making use of an extended symmetry at an exceptional point) as the map from $D^*$ to $D$ is symmetry dependent.
If one sets $D$ as a target, it may not be possible to reach it exactly. One then has to decide how to set $D^*$: possible choices are largest value before $D$ is reached, smallest value above $D$, or closest value to $D$.
A key missing feature in
BlockSparseArrays, as mentioned in #1336, is blockwise matrix factorizations. I'm starting an issue here to sketch out an implementation plan for that.Here is an initial implementation of a blockwise SVD, using the new submodules of
NDTensors.jl:Here's a demonstration that it works:
There are a number of improvements to make to the code above:
Sas aBlockSparseMatrixwhere the blocks areDiagonalorDiagonalMatrixinstead ofMatrix. That would be relatively easy to support, I just haven't tested that and there are probably a few issues to work out inBlockSparseArraysto support that.Matrixwhich is on CPU.BlockSparseMatrixrepresentations of $U$, $S$, and $V^{\dagger}$. In the code above I first store the results of SVDing each block in dictionaries that store the blocks of $U$, $S$, and $V^{\dagger}$, and then build theBlockSparseMatrixrepresentations of $U$, $S$, and $V^{\dagger}$ from those blocks and also axes determined from the input matrix and the results of the SVD. I think that logic is pretty good and simple but maybe it can be simplified a little bit. For example, we could store the blocks in aSparseMatrixDOKand then convert it to aBlockSparseMatrix, which could automatically determine the block sizes.BlockSparseMatrixinputs that invariant under group symmetries. Really the only thing that needs to be changed above is adding on symmetry labels/irrep labels to the blocks of the axes that get attached to the internal dimensions (those connecting $S$ to $U$ and $V^{\dagger}$, calledr_sin the code above). The logic for that would be pretty simple, and be uniquely determined using the convention that $U$ and $V^{\dagger}$ are in trivial symmetry sectors (have zero flux).For some references on other implementations of blockwise matrix factorizations, see https://github.com/ITensor/ITensors.jl/blob/v0.6.16/NDTensors/src/lib/BlockSparseArrays/src/backup/qr.jl and https://github.com/JuliaArrays/BlockDiagonals.jl/blob/master/src/linalg.jl.
@ogauthe @emstoudenmire