TEBD Truncation Error

[Joacop16]

Dear ITensor Team!

I have implemented TEBD for a research project, and now I am wondering if there is a way to track the accumulated error.

For Illustration purposes, consider a very simple code

using ITensors
using ITensorMPS

N = 10
sites = siteinds("Fermion", N)
ψ = random_mps(sites, linkdims = 1)
TEBD_Gates = [exp(-1im*(op("Cdag", sites[i])*op("C", sites[i+1]) + op("C", sites[i+1])*op("Cdag", sites[i]))) for i = 1:N-1];

ψ_t = ψ
for i = 1:20
    ψ_t = apply(TEBD_Gates, ψ_t; maxdim = 32, cutoff = 0.0)
end

For N = 10, the maximum bond dimension would be 2 ^5 = 32. There is no truncation error (although some error may arise from the Taylor expansion when comparing with the exact time evolution operator). However, depending on the values of cutoff and maxdim, this can change.

I know that cutoff is the truncation error considered by svd() insideapply(). I was testing the SVD function from ITensors, and I realized that given some cutoff, svd() keeps as many singular values as it is required to have \cfrac{\sum_{k=\chi + 1}^{32} \lambda_k^2 }{\sum_{k=1}^{32} \lambda_k^2} = cutoff.

I assume that if svd() works in this way, then apply() should work similarly. So, if I apply all TEBD Gates, it tries to take at most singular values to guarantee a given cutoff. Now in my real code, I have a system of electrons with N=48 (it is impossible to keep all information), I am not sure what maxdim would be enough to represent the state, and I need to perform a very long time evolution, so I usually do something like apply(TEBD_Gates, ψ_t; maxdim = 800, cutoff = 1e-10) because it has a decent running time. After some sweeps, the bond dimension of the MPS becomes 800 in almost all bonds, so at that point, the code is actually truncating by maxdim instead of cutoff. So, my truncation error is not really cutoff (If I do not include maxdim inside apply, the bond dimension keeps increasing without limit).

Is there a way to know what is the truncation error from the truncation done by apply with maxdim? I do not have access to all singular values, so I can not calculate it by myself. Am I missing something? Is there a better measurement of accumulated error for TEBD? I just want some notion of truncation error to know if I can trust on my TEBD results or if I should increase the bond dimension.

[mtfishman]

It’s a good question. Ideally the apply function would add up all of the truncation errors of all of the gate applications and output that, but we haven’t added that functionality yet (I think someone started working on it at some point but it wasn’t completed).

I believe the only thing that would need to be changed would be to output the truncation error from this line: ITensorMPS.jl/src/abstractmps.jl at v0.3.25 · ITensor/ITensorMPS.jl · GitHub and then add up all of the truncation errors.

[mtfishman]

It turns out it was pretty easy to implement, see: [WIP] Access truncation errors in apply by mtfishman · Pull Request #168 · ITensor/ITensorMPS.jl · GitHub .

You would use it like this:

using ITensorMPS: apply, maxlinkdim, op, random_mps, siteinds
using StableRNGs: StableRNG

n = 30
nlayers = 20
cutoff = 1e-3
s = siteinds("S=1/2", n)
rng = StableRNG(123)
x = random_mps(rng, s)
U1 = [op("CX", s[i], s[i + 1]) for i in 1:2:(n - 1)]
U2 = [op("CX", s[i], s[i + 1]) for i in 2:2:(n - 2)]
U = [U1; U2]
U = reduce(vcat, repeat(U, nlayers))

# Define a callback function that will be called after the `factorize` function
# where the truncation is performed.
truncerr = Ref(0.0)
post_factorize_callback(; truncation_error) = (truncerr[] += truncation_error)

# Call `apply` with the callback function to accumulate the truncation errors.
Ux = apply(U, x; cutoff, post_factorize_callback)

@show maxlinkdim(Ux)
@show truncerr[]
@show length(U)
@show cutoff * length(U)
@show truncerr[] < cutoff * length(U)

which outputs:

maxlinkdim(Ux) = 153
truncerr[] = 0.23065869540430445
length(U) = 580
cutoff * length(U) = 0.58
truncerr[] < cutoff * length(U) = true

[Joacop16]

Thank you so much for your rapid response.

I tested your code, and it works. However, when I tried to implement it in my research project, I got something weird.

If I understood correctly cutoff*length(TEBD_GATES) would be the maximum error that we can have and still trust on the results. I was expecting to see that truncerr[] becomes bigger than this value at some point during the time evolution, probably when the bond dimension stops growing. However, it is happening even at the first step of the time evolution when the bond dimension is still smaller than the maxdim I have in the parameters.

So, I am running:

truncerr = Ref(0.0)
post_factorize_callback(; truncation_error) = (truncerr[] += truncation_error)        
Ψ_t = apply(TEBD_GATES, Ψ_t; cutoff=1e-10, maxdim=800, post_factorize_callback); 
@show truncerr[], maxlinkdim(Ψ_t)

In this case length(TEBD_GATES) = 18, so cutoff*length(TEBD_GATES) = 1.8e-9. In contrast to this, my output says:

(truncerr[], maxlinkdim(Ψ_t)) = (2.3708737771223453e-9, 59)

i.e. truncerr[] > cutoff*length(TEBD_GATES), but maxlinkdim(Ψ_t)<800, so in principle the computer could have taken more singular values to guarantee a better truncation error, why only 59?

I was thinking and I believe it may be related with the fact that each TEBD Gate is acting on 4 electron sites (i.e each gate is a 256x256 matrix) instead of 2 like in the common TEBD algorithm. Usually contracting a 2-site gate with an MPS requires 1 SVD, I supposed that a 4-site gate maybe requires more than 1 SVD, so probably the maximum error would not be simply cutoff*length(TEBD_GATES). Does that makes sense? How does it works apply() for 4-site gates?

[mtfishman]

For 4-site gates, it would contract the gate with 4 sites of the MPS and then perform 3 SVDs (one for each internal bond) to split it back up into 4 MPS tensors to insert back into the MPS.

You could get more information by storing all of the truncation errors:

truncerrs = Float64[]
post_factorize_callback(; truncation_error) = push!(truncerrs, truncation_error)
Ux = apply(U, x; cutoff, post_factorize_callback)