Why using weighted singular values for the truncation?

[Phy-David-Zhang]

In the file truncation.jl, I saw the following function:

# auxiliary function
function _findnexttruncvalue(Σdata, truncdim::SectorDict{I,Int}, p::Real) where {I<:Sector}
    # early return
    (isempty(Σdata) || all(iszero, values(truncdim))) && return nothing

    # find some suitable starting candidate
    cmin = findfirst(>(0), truncdim)
    weightedσmin = dim(cmin)^inv(p) * Σdata[cmin][truncdim[cmin]]

    # find the actual minimum singular value
    for (c, σs) in Σdata
        if truncdim[c] > 0
            weightedσ = dim(c)^inv(p) * σs[truncdim[c]]
            if weightedσ < weightedσmin
                cmin, weightedσmin = c, weightedσ
            end
        end
    end
    return cmin
end

It seems to me that the truncation is based on the weighted singluar values, as

weightedσmin = dim(cmin)^inv(p) * Σdata[cmin][truncdim[cmin]]

and

weightedσ = dim(c)^inv(p) * σs[truncdim[c]]

In the non-Abelian case, singular values are weighted by the square root of the degeneracy of the corresponding sector. I am wondering why the singular values should be weighted in that way? Usually, one only truncates based on the singular values without weights. I ask because in some cases, the spin-0 sectors have the largest singular values, but smallest degeneracy. So sometimes the truncation truncates out the most important spin-0 sector...

[Jutho]

The motivation for this comes from minimising the truncation error. The total truncation error using the $p$-norm is

$$\left(\sum_{\text{sectors $c$}} d_c \sum_{\text{truncated $\sigma$ in sector $c$}} \sigma^p\right)^{1/p} = \left(\sum_{\text{sectors $c$}} \sum_{\text{truncated $\sigma$ in sector $c$}} (d_c^{1/p} \sigma)^p \right)^{1/p}$$

For example, for $p=2$, if you have a singular value 10e-4 in the spin=0 sector, and a singular value 8e-4 in the spin 1/2 sector, truncating the latter will actually have a larger effect on the total truncation error then truncating the former.

Note that there is also TruncationCutoff with the truncbelow constructor to just truncate singular values based their actual value, without taking degeneracies into account.
Alternatively, you can also use p=Inf to kill the effect of taking the degeneracy into account.

[lkdvos]

To be fair, I do agree that it is somewhat of an unconventional choice. I think that a somewhat contrived counterexample would be to have a tensor with singular values SectorDict(SU2Irrep(0) => fill(0.5, 20), SU2Irrep(3//2) => fill(0.49, 1)). In the case where you want to now truncate 4 values, the algorithm would find that it is beneficial to throw away 4 values in the trivial sector, while clearly discarding to SU2Irrep(3//2) irrep would lead to a lower truncation error.
The point being that since our approach is not a global minimizer anyways, the choice might be a bit odd, as truncating tensors with different subgroups imposed would happen differently.

[Jutho]

Ok I had not considered that, so the problem is that selecting for removal one-by-one is not leading the optimal solution. Not a surprise in hindsight, but something I surely overlooked in the design/development.

[lkdvos]

I think it is reasonable to just always select the smallest ones, and just overcut a little bit. At the very least, this would give behavior that’s somewhat consistent, and doesn’t change the way the values are selected between symmetric and non-symmetric simulations

[Gertian]

I second that. In my (humble) opinion, the default should behave the same for symmetric and non-symmetric code.
It is also my understanding the manual is written with that philosophy in mind ? If not this should at least be made more clear.

[Jutho]

Well it's never going to be possible to be exactly the same, since in the symmetric case you can cut away complete multiplets. But I am not bound to this behaviour and will thus change it (and maybe give it a few moments of thought if I can come up with a better strategy to determine what to cut altogether).