[ITensorMPS] [BUG] Default algorithm for adding MPSs leading to large bond dimensions

[shinaoka]

Thank you for your efforts in the development.

Description of bug

Adding an MPS with an MPS representing 0 with the default algorithm leads to a larger bond dimension than that obtained with the directsum algorithm

Minimal code demonstrating the bug or unexpected behavior

Minimal runnable code

using ITensors
using Random

function zeromps(::Type{T}, sites) where {T<:Number}
    M = MPS(T, sites; linkdims=1)
    l = linkinds(M)
    for n in eachindex(M)
        if n == 1
            M[n] = ITensor(T, sites[n], l[n])
        elseif n == length(M)
            M[n] = ITensor(T, l[n - 1], sites[n])
        else
            M[n] = ITensor(T, l[n - 1], sites[n], l[n])
        end
        M[n] .= zero(T)
    end
    return M
end

L = 10
sites = [Index(2, "s=$n") for n in 1:L]

Random.seed!(1234)
A = zeromps(ComplexF64, sites)
B = randomMPS(ComplexF64, sites; linkdims=10)
AB = +(A, B; cutoff=1e-25, maxdim=typemax(Int))
@show maxlinkdim(AB)
truncate!(AB)
@show maxlinkdim(AB)

AB_ref = +(A, B; alg="directsum")
truncate!(AB_ref; cutoff=1e-25, maxdim=typemax(Int))
@show maxlinkdim(AB_ref)

Expected output or behavior

Since A=0, we should be able to restore the original bond dimension (=10).

Actual output or behavior

The bond dimension of the result is much larger than the expected one, even after an explicit truncation.

Output of minimal runnable code

maxlinkdim(AB) = 51
maxlinkdim(AB) = 25
maxlinkdim(AB_ref) = 10

Version information

• Output from versioninfo():

julia> versioninfo()
Julia Version 1.10.4
Commit 48d4fd48430 (2024-06-04 10:41 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin22.4.0)
  CPU: 8 × Apple M2
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, apple-m1)
Threads: 1 default, 0 interactive, 1 GC (on 4 virtual cores)
Environment:
  JULIA_NUM_THREADS = 1
  JULIA_PKG_USE_CLI_GIT = true
  JULIA_PROJECT = @.

• Output from using Pkg; Pkg.status("ITensors"):

julia> using Pkg; Pkg.status("ITensors")
  [9136182c] ITensors v0.6.21

[mtfishman]

Thanks for the report @shinaoka, what if you set the cutoff to something bigger than eps(Float64), say 1e-15?

[shinaoka]

Thank you for the quick reply. With 1e-16, I got the correct bond dimension 10.

BTW, I am sitting on the 9th floor of Flatiron Institute now.

[mtfishman]

Thank you for the quick reply. With 1e-16, I got the correct bond dimension 10.

Ok, thanks for confirming. As I suspected, I think this is due to floating point roundoff errors, but probably we can hard code an upper bound on the final bond dimensions based on the input bond dimensions, say as the sum of the input bond dimensions on each bond.

BTW, I am sitting on the 9th floor of Flatiron Institute now.

Sounds good, how long will you be around for? I'm on parental leave, I'll be back to work officially October 28.

[shinaoka]

probably we can hard code an upper bound on the final bond dimensions based on the input bond dimensions

Thank you for the suggestion. But, we would rather use the directsum algorithm because we observed this one is more robust against round errors in many other cases. Do you expect any drawbacks to using the directsum algorithm?

I will be at CCQ until the middle of December.

[mtfishman]

probably we can hard code an upper bound on the final bond dimensions based on the input bond dimensions

Thank you for the suggestion. But, we would rather use the directsum algorithm because we observed this one is more robust against round errors in many other cases. Do you expect any drawbacks to using the directsum algorithm?

EDIT (since I misread your question originally): The "densitymatrix" algorithm has better scaling than the "directsum" algorithm, which is why it is the default. Probably there is worse roundoff errors because it involves forming and then truncating the density matrix of the sum of states, which involves squaring the singular values of the sum of states and truncating according to those squared singular values, which loses precision beyond roughly sqrt(eps(scalartype)). Probably the best algorithm in terms of scaling and accuracy would be to first use the density matrix algorithm, and then use a variational fitting algorithm to improve the accuracy, though we don't have a fitting algorithm implemented right now for sums of states. I think it would be relatively easy to implement based on the generic alternating update solver code in ITensorTDVP.jl.

I will be at CCQ until the middle of December.

Great, I'll see you at CCQ then.

[shinaoka]

a variational fitting algorithm to improve the accuracy, though we don't have a fitting algorithm implemented right now for sums of states.

Thank you for the explanation. Yes, I have an implementation of the variational algorithm for sums of states. I will give a try to densitymatrix + fitting.

[mtfishman]

a variational fitting algorithm to improve the accuracy, though we don't have a fitting algorithm implemented right now for sums of states.

Thank you for the explanation. Yes, I have an implementation of the variational algorithm for sums of states. I will give a try to densitymatrix + fitting.

Let me know how it goes, I'll be interested to hear if that works well.

[shinaoka]

Sorry, just one more question.

The "densitymatrix" algorithm has better scaling than the "directsum" algorithm, which is why it is the default.

What specifically does “better scaling” refer to? I naively assumed that the computational time of both algorithms scales as $O(\chi^3)$. The scaling of the directsum algorithm should be $O(d^3 \chi^3)$, where $d$ is the local dimension. Does the densitymatrix algorithm have smaller powers for $\chi$ or $d$?

[mtfishman]

Yeah, saying it has better scaling may not have been quite precise enough, I think the difference is in the prefactors. But the basic idea is that with the "directsum" method, for adding two MPS with bond dimension $\chi$, you first construct the state with bond dimension $2\chi$, and then truncate with scaling $O((2\chi)^3)$. While with "densitymatrix", the truncation is performed on the fly, so if the sum of the two states leads to a state with bond dimension $\chi$ you never have to work with a state of bond dimension $2\chi$. So maybe the improvement is roughly a factor of $8$ but I would have to work it out more carefully.

If you are summing many MPS, say $n$ MPS, you see even more of an improvement in the cost.

[shinaoka]

I observed that "densitymatrix" + "fit" restores the correct bond dimension.

using ITensors
using Random
import FastMPOContractions as FMPOC

L = 10
sites = [Index(2, "s=$n") for n in 1:L]

Random.seed!(1234)
A = 0.0 * randomMPS(ComplexF64, sites; linkdims=10)
B = randomMPS(ComplexF64, sites; linkdims=10)
AB = +(A, B; cutoff=1e-25, maxdim=typemax(Int))
@show maxlinkdim(AB)
truncate!(AB)
@show maxlinkdim(AB)

function _fitsum(
    input_states::AbstractVector{T},
    init::T;
    coeffs::AbstractVector{<:Number}=ones(Int, length(input_states)),
    kwargs...,
) where {T}
    if !(:nsweeps ∈ keys(kwargs))
        kwargs = Dict{Symbol,Any}(kwargs)
        kwargs[:nsweeps] = 1
    end
    Ψs = [MPS(collect(x)) for x in input_states]
    init_Ψ = MPS(collect(init))
    res = FMPOC.fit(Ψs, init_Ψ; coeffs=coeffs, kwargs...)
    return T(collect(res))
end

AB_ref = +(A, B; alg="directsum")
truncate!(AB_ref; cutoff=1e-25, maxdim=typemax(Int))
@show maxlinkdim(AB_ref)

AB_fit = _fitsum([A, B], AB; cutoff=1e-25, maxdim=typemax(Int))
@show maxlinkdim(AB_fit)

maxlinkdim(AB) = 52
maxlinkdim(AB) = 29
maxlinkdim(AB_ref) = 10
maxlinkdim(AB_fit) = 10

[mtfishman]

Thanks for confirming that. For now I think the simplest solution which would fix the original problem you brought up is to set a maximum dimension of the output link indices as the sum of the bond dimensions of the input link indices, analogous to what we do in the "densitymatrix" backend of MPO-MPS contraction: https://github.com/ITensor/ITensors.jl/blob/v0.6.22/src/lib/ITensorMPS/src/mpo.jl#L753.