Problems with Contracting MPO with cutoff

sourav · October 30, 2024, 2:35pm

Hi,
I am having a strange problem with contracting ITensors using the contract function. I had pointed this out in a previous post in a different context. Since this problem seems to go beyond the context in which I had asked it previously, I thought to open a new topic which might be useful. Let me state my problem with an example.

I give a simple example where I am evolving a density matrix by MPO evolution in TEBD. A simple code is

N, cf, tau, ttotal= 6, 1e-10, 0.1, 4.0
s = siteinds("S=1/2", N)
chi = 60
δt=0.1
ttotal=5.0

gates = ITensor[]
for j in 1:(N - 1)
    s1 = s[j]
    s2 = s[j + 1]
    hj =
    op("Sz", s1) * op("Sz", s2) +
    1 / 2 * op("S+", s1) * op("S-", s2) +
    1 / 2 * op("S-", s1) * op("S+", s2)
    Gj = exp(-im * tau / 2 * hj)
    push!(gates, Gj)
end

Nsteps = Int(ttotal / δt)

append!(gates, reverse(gates))
psi = MPS(s, n -> isodd(n) ? "Up" : "Dn")
c = div(N, 2) 
rho = MPO(psi)
t = 0.0
for step in 1:Nsteps
    rho = apply(gates, rho; cutoff=cf, apply_dag=true)
    t += δt
end

With this I have resulting MPO with structure

MPO
[1] ((dim=2|id=650|"S=1/2,Site,n=1")', (dim=2|id=650|"S=1/2,Site,n=1"), (dim=4|id=122|"Link,n=1"))
[2] ((dim=4|id=122|"Link,n=1"), (dim=2|id=731|"S=1/2,Site,n=2")', (dim=16|id=814|"Link,n=1"), (dim=2|id=731|"S=1/2,Site,n=2"))
[3] ((dim=2|id=850|"S=1/2,Site,n=3")', (dim=64|id=61|"Link,n=1"), (dim=2|id=850|"S=1/2,Site,n=3"), (dim=16|id=814|"Link,n=1"))
[4] ((dim=2|id=631|"S=1/2,Site,n=4")', (dim=16|id=251|"Link,n=1"), (dim=2|id=631|"S=1/2,Site,n=4"), (dim=64|id=61|"Link,n=1"))
[5] ((dim=2|id=720|"S=1/2,Site,n=5")', (dim=4|id=464|"Link,n=1"), (dim=2|id=720|"S=1/2,Site,n=5"), (dim=16|id=251|"Link,n=1"))
[6] ((dim=2|id=938|"S=1/2,Site,n=6")', (dim=2|id=938|"S=1/2,Site,n=6"), (dim=4|id=464|"Link,n=1"))

Now if I want to contract the local two tensors corresponding to site 2 and 3 of the resulting MPO, I simply do

contract(rho[2]', rho[3])

and I get

ITensor ord=8 (dim=4|id=122|"Link,n=1")' (dim=2|id=731|"S=1/2,Site,n=2")'' (dim=16|id=814|"Link,n=1")' (dim=2|id=731|"S=1/2,Site,n=2")' (dim=2|id=850|"S=1/2,Site,n=3")' (dim=64|id=61|"Link,n=1") (dim=2|id=850|"S=1/2,Site,n=3") (dim=16|id=814|"Link,n=1")
NDTensors.Dense{ComplexF64, Vector{ComplexF64}}

which is all fine. Problem starts if I put a cutoff for the contractions. So if I do

contract(rho[2]', rho[3]; cutofff=1e-8)

I immediately get the output

ERROR: StackOverflowError:
Stacktrace:
     [1] mapfoldl_impl(f::typeof(identity), op::ITensors.var"#213#215"{…}, nt::Base._InitialValue, itr::Tuple{…})
       @ Base ./reduce.jl:44
     [2] mapfoldl(f::Function, op::Function, itr::Tuple{ITensor, ITensor}; init::Base._InitialValue)
       @ Base ./reduce.jl:175
     [3] mapfoldl
       @ ./reduce.jl:175 [inlined]
     [4] foldl(op::Function, itr::Tuple{ITensor, ITensor}; kw::@Kwargs{})
       @ Base ./reduce.jl:198
     [5] contract(As::Tuple{ITensor, ITensor}; sequence::String, kwargs::@Kwargs{cutofff::Float64})
       @ ITensors ~/.julia/packages/ITensors/23eRw/src/tensor_operations/tensor_algebra.jl:158
     [6] contract(::ITensor, ::Vararg{ITensor}; kwargs::@Kwargs{cutofff::Float64})
       @ ITensors ~/.julia/packages/ITensors/23eRw/src/tensor_operations/tensor_algebra.jl:168
     [7] (::ITensors.var"#213#215"{@Kwargs{cutofff::Float64}})(A::ITensor, B::ITensor)
       @ ITensors ~/.julia/packages/ITensors/23eRw/src/tensor_operations/tensor_algebra.jl:158
     [8] (::Base.BottomRF{ITensors.var"#213#215"{@Kwargs{cutofff::Float64}}})(acc::ITensor, x::ITensor)
       @ Base ./reduce.jl:86
     [9] afoldl(::Base.BottomRF{ITensors.var"#213#215"{@Kwargs{cutofff::Float64}}}, ::Base._InitialValue, ::ITensor, ::ITensor)
       @ Base ./operators.jl:545
    [10] _foldl_impl(op::Base.BottomRF{ITensors.var"#213#215"{…}}, init::Base._InitialValue, itr::Tuple{ITensor, ITensor})
       @ Base ./reduce.jl:68
    [11] foldl_impl(op::Base.BottomRF{ITensors.var"#213#215"{@Kwargs{…}}}, nt::Base._InitialValue, itr::Tuple{ITensor, ITensor})
       @ Base ./reduce.jl:48
--- the last 11 lines are repeated 3083 more times ---
 [33925] mapfoldl_impl(f::typeof(identity), op::ITensors.var"#213#215"{…}, nt::Base._InitialValue, itr::Tuple{…})
       @ Base ./reduce.jl:44
 [33926] mapfoldl(f::Function, op::Function, itr::Tuple{ITensor, ITensor}; init::Base._InitialValue)
       @ Base ./reduce.jl:175
 [33927] mapfoldl
       @ ./reduce.jl:175 [inlined]
 [33928] foldl(op::Function, itr::Tuple{ITensor, ITensor}; kw::@Kwargs{})
       @ Base ./reduce.jl:198
 [33929] contract(As::Tuple{ITensor, ITensor}; sequence::String, kwargs::@Kwargs{cutofff::Float64})
       @ ITensors ~/.julia/packages/ITensors/23eRw/src/tensor_operations/tensor_algebra.jl:158
Some type information was truncated. Use `show(err)` to see complete types.

I have tried may ways to understand what is happening / what might be going wrong. But I could not resolve this. I would be grateful if you can help me out.

Thanks,
Sourav.

mtfishman · October 30, 2024, 2:43pm

When you call contract(rho[2]', rho[3]), you are contracting the individual ITensors of the MPO (the ones on sites 2 and 3) which each other, and we don’t support approximate contraction of ITensors. What’s your goal with contracting those two tensors together?

sourav · October 30, 2024, 2:48pm

Thanks for the quick answer. Well, the original goal was to compute some quantity as shown in this screenshot:

For this objective I need some contractions (of course here the contraction is between two ITensors sitting a site for a density matrix and its two copies–which can be accordingly implemented). But essentially needed the ‘contractions’ for this reason.

Best,
Sourav

mtfishman · October 30, 2024, 2:54pm

Thanks, that’s helpful. Have you tried performing the contractions without approximations? I guess you were hoping to perform approximate contraction to get better performance? If you want to perform the contractions approximately, probably the best way will be to call contract/* to do exact contractions, and then call low rank factorizations like svd yourself on the intermediate tensors (that’s essentially how approximate MPS/MPO contraction methods in ITensorMPS.jl are implemented).

sourav · October 30, 2024, 3:04pm

Without approximation everything works but I can’t do that for large systems because my system has local Hilbert space dimension 9 (it is originally 3 state qubit but it is open systems and followed by vectorization the local dim becomes 3^2=9) and hence tensors are pretty big–that’s why I wanted to put the cutoff. Can you give me one example of the method you pointed out?
Thanks in adv.

Best,
Sourav

mtfishman · October 30, 2024, 3:07pm

You could look here: ITensorMPS.jl/src/mpo.jl at v0.3.1 · ITensor/ITensorMPS.jl · GitHub

sourav · October 31, 2024, 5:35pm

thanks…
But I was thinking if I could use svd (instead of factorization as shown in the example you pointed out) like the following (say using my previous simple case example):


u, s, v = svd(rho[2]'*rho[3], uniqueinds(rho[2]', rho[3]); cutoff = 1e-8) 

contracted_form = u*s*v

I checked with

@show rho[2]'*rho[3] ≈ u*s*v

and this yields ‘true’–But still I thought to confirm that I may not be missing something important which might be easy to overlook?

Thanks,
Sourav

mtfishman · October 31, 2024, 5:38pm

That’s the correct usage of svd, I’m not really sure what your question is.

sourav · October 31, 2024, 5:40pm

Yes, I wanted to ask if I am using svd correctly–sorry for using rather confusing sentences–thanks for the answer.

Topic		Replies	Views
Bug in multiplying MPOs and tracing MPOs ITensor Julia Questions julia , mpo	2	340	October 23, 2023
TEBD on 2D fermionic system ITensor Julia Questions	3	131	September 12, 2024
Exponential of a three-site MPO ITensor Julia Questions julia , mpo	7	874	January 23, 2024
Constructing an MPO using apply ITensor Julia Questions	2	57	October 14, 2024
How to contract 2 MPS objects ITensor Julia Questions	2	543	January 30, 2023

Problems with Contracting MPO with cutoff

Related topics