VUMPS for real-valued iMPS

[stecrotti]

Hi!

VUMPS truncation is not working for real-valued infinite MPSs.
MWE:

using TensorKit, MPSKit

m = 10  # original bond dimension
d = 5   # target bond dimension
q = 3   # physical dimension

A = rand(m, q, m)  # matrix A for iMPS to be truncated
ψ₀ = InfiniteMPS([TensorMap(A, (ℝ^size(A, 1) ⊗ ℝ^q), ℝ^size(A, 1))])

init = rand(d, q, d)   # initalization for output iMPS
ψ = InfiniteMPS([TensorMap(init, (ℝ^d ⊗ ℝ^size(A,2)), ℝ^d)])

# Pick the identity as MPO so that <MPS|MPO|MPS> reduces to the overlap <MPS|MPS> in the optimization
Id = DenseMPO([MPSKit.add_util_leg(id(storagetype(MPSKit.site_type(ψ₀)), physicalspace(ψ₀, i)))
    for i in 1:length(ψ₀)])

# truncate using VUMPS
ψ, = approximate(ψ, (Id, ψ₀), VUMPS(; maxiter=100, tol=1e-12))

raises an error whenever a certain dot product λ in the mixed_fixpoints subroutine happens to be negative, because its square root is taken. This is not a problem in the general case where λ is complex.

ERROR: LoadError: TaskFailedException

    nested task error: DomainError with -0.1618786374563001:
    sqrt was called with a negative real argument but will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
    Stacktrace:
     [1] throw_complex_domainerror(f::Symbol, x::Float64)
       @ Base.Math ./math.jl:33
     [2] sqrt
       @ ./math.jl:608 [inlined]
     [3] macro expansion
       @ ~/.julia/packages/MPSKit/diLPU/src/environments/permpoinfenv.jl:183 [inlined]

...

(julia v1.11.1, MPSKit v0.11.4, TensorKit v0.12.7)

Happy to contribute with a PR if you think it's worth doing something about this.

Background

I'm using MPSKit for VUMPS-based truncations of infinite MPSs that are used to parametrize the probability of a classical particle evolving in time. It's been very useful, so thanks, especially to @leburgel who helped previously! Should anyone be interested this is the article.

The main point is that my iMPSs are real-valued, because they directly represent probability distributions, instead of their "square roots" as in a quantum setting. We need to truncate these iMPSs to smaller bond dimension, while keeping them real-valued. I'm sure this is doable, since all the subroutines in VUMPS (computation of leading eigenvalues, SVDs, etc.) return real numbers if initialized within the reals. For now i'm working with a personal fork of MPSKit that does the easy fix sqrt(λ) -> sqrt(abs(λ)), but it would be nice to reference the official package as we make our code available after publication.

[lkdvos]

Thanks for raising the issue and explaining the background, that helps a lot!

The main point of these lines of code is to get some form of normalization convention for the environments, such that the overlap of the left and right part is 1. Currently, we do this by symmetrically distributing the factors over left and right (hence the sqrt), but in principle this is not really necessary.
One option would be to simply always pick either the left or the right side to absorb that factor, thus removing the need for the sqrt altogether.
The other is to split the magnitude evenly over the left and right side, but absorb the phase only left or right.
I would slightly favor the first option, if it's not symmetric anyways we might as well not go through the additional steps of trying to make them the same magnitude.
I could of course be missing something, maybe @maartenvd could have a quick look that my reasoning is sound?

[lkdvos]

I pushed a changed version that should have fixed this, which is now also tagged. Could you verify that this is now fixed for you, and let me know if the problem remains, or anything else shows up?

[stecrotti]

It works, thanks!