# Problems with Contracting MPO with cutoff **URL:** \ **Category:** ITensor Julia Questions\ **Created:** [October 30, 2024, 2:35pm UTC](https://itensor.discourse.group/t/problems-with-contracting-mpo-with-cutoff/2096 "2024-10-30T14:35:20Z")\ **Posts on this page:** 9\ **Page:** 1
**Author:** ![sourav](https://avatars.discourse-cdn.com/v4/letter/s/a4c791/32.png) [@sourav](https://itensor.discourse.group/u/sourav)\ **Post date:** [October 30, 2024, 2:35pm UTC](https://itensor.discourse.group/t/problems-with-contracting-mpo-with-cutoff/2096/1 "2024-10-30T14:35:20Z")
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 ```auto 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 ```auto 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 ```auto contract(rho[2]', rho[3]) ``` and I get ```auto 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 ```auto contract(rho[2]', rho[3]; cutofff=1e-8) ``` I immediately get the output ```auto 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. ---
**Author:** ![mtfishman](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/mtfishman/32/11_2.png) [@mtfishman](https://itensor.discourse.group/u/mtfishman)\ **Post date:** [October 30, 2024, 2:43pm UTC](https://itensor.discourse.group/t/problems-with-contracting-mpo-with-cutoff/2096/2 "2024-10-30T14:43:03Z")
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? ---
**Author:** ![sourav](https://avatars.discourse-cdn.com/v4/letter/s/a4c791/32.png) [@sourav](https://itensor.discourse.group/u/sourav)\ **Post date:** [October 30, 2024, 2:48pm UTC](https://itensor.discourse.group/t/problems-with-contracting-mpo-with-cutoff/2096/3 "2024-10-30T14:48:53Z")
Thanks for the quick answer. Well, the original goal was to compute some quantity as shown in this screenshot: ![Screenshot 2024-10-30 at 15.45.29](https://global.discourse-cdn.com/free1/uploads/itensor/original/1X/0fc40657fb53978211026de98ba6deae316a52bd.png) 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 ---
**Author:** ![mtfishman](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/mtfishman/32/11_2.png) [@mtfishman](https://itensor.discourse.group/u/mtfishman)\ **Post date:** [October 30, 2024, 2:54pm UTC](https://itensor.discourse.group/t/problems-with-contracting-mpo-with-cutoff/2096/4 "2024-10-30T14:54:21Z")
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). ---
**Author:** ![sourav](https://avatars.discourse-cdn.com/v4/letter/s/a4c791/32.png) [@sourav](https://itensor.discourse.group/u/sourav)\ **Post date:** [October 30, 2024, 3:04pm UTC](https://itensor.discourse.group/t/problems-with-contracting-mpo-with-cutoff/2096/5 "2024-10-30T15:04:20Z")
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 ---
**Author:** ![mtfishman](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/mtfishman/32/11_2.png) [@mtfishman](https://itensor.discourse.group/u/mtfishman)\ **Post date:** [October 30, 2024, 3:07pm UTC](https://itensor.discourse.group/t/problems-with-contracting-mpo-with-cutoff/2096/6 "2024-10-30T15:07:54Z")
You could look here: [ITensorMPS.jl/src/mpo.jl at v0.3.1 · ITensor/ITensorMPS.jl · GitHub](https://github.com/ITensor/ITensorMPS.jl/blob/v0.3.1/src/mpo.jl#L823) ---
**Author:** ![sourav](https://avatars.discourse-cdn.com/v4/letter/s/a4c791/32.png) [@sourav](https://itensor.discourse.group/u/sourav)\ **Post date:** [October 31, 2024, 5:35pm UTC](https://itensor.discourse.group/t/problems-with-contracting-mpo-with-cutoff/2096/7 "2024-10-31T17:35:11Z")
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): ```auto u, s, v = svd(rho[2]'*rho[3], uniqueinds(rho[2]', rho[3]); cutoff = 1e-8) contracted_form = u*s*v ``` I checked with ```auto @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 ---
**Author:** ![mtfishman](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/mtfishman/32/11_2.png) [@mtfishman](https://itensor.discourse.group/u/mtfishman)\ **Post date:** [October 31, 2024, 5:38pm UTC](https://itensor.discourse.group/t/problems-with-contracting-mpo-with-cutoff/2096/8 "2024-10-31T17:38:29Z")
That’s the correct usage of `svd`, I’m not really sure what your question is. ---
**Author:** ![sourav](https://avatars.discourse-cdn.com/v4/letter/s/a4c791/32.png) [@sourav](https://itensor.discourse.group/u/sourav)\ **Post date:** [October 31, 2024, 5:40pm UTC](https://itensor.discourse.group/t/problems-with-contracting-mpo-with-cutoff/2096/9 "2024-10-31T17:40:54Z")
Yes, I wanted to ask if I am using svd correctly–sorry for using rather confusing sentences–thanks for the answer.