Sum of two MPOs

EmilioRui · March 16, 2023, 3:18pm

The summation of two MPOs is theoretically very easy, given two MPOs:
H_1= \sum A_{\alpha_{i-1},\alpha_i}^{s_i s_i'}\ket{s_i}\bra{s_i'}
H_2= \sum B_{\beta_{i-1},\beta_i}^{s_i s_i'}\ket{s_i}\bra{s_i'}
H_1 + H_2 = \sum C_i^{s_i s_i'}\ket{s_i}\bra{s_i'}
where
C_i = A_i \oplus B_i

Is there a built in way to sum two MPOs?

phx · March 16, 2023, 3:53pm

I believe you can just sum them with a +, like, x + y given two MPOs x and y, otherwise use add(x, y; kwargs...) if you want to specify some additional arguments such as a cutoff for the compression of the result.
I can’t find this function in the documentation, but it is implemented here (at line 1375) in the package code.

EmilioRui · March 16, 2023, 4:23pm

Thank you!
I kept on receiving errors, but it was due to a bad definition of the MPOs.
I was using the same link index for the two MPOs (same id) and so it was contracting them in unwanted manner I guess

EmilioRui · March 16, 2023, 4:51pm

Actually, I cannot find where I am making mistakes

sites = siteinds("Boson",3,dim=7)

If I create random MPO I can sum them

MPO_1 = randomMPO(sites)
MPO_2 = randomMPO(sites)

MPO_1
MPO
[1] ((dim=7|id=528|"Boson,Site,n=1")', (dim=7|id=528|"Boson,Site,n=1"), (dim=1|id=994|"Link,l=1"))
[2] ((dim=7|id=826|"Boson,Site,n=2")', (dim=7|id=826|"Boson,Site,n=2"), (dim=1|id=850|"Link,l=2"), (dim=1|id=994|"Link,l=1"))
[3] ((dim=7|id=755|"Boson,Site,n=3")', (dim=7|id=755|"Boson,Site,n=3"), (dim=1|id=850|"Link,l=2"))

MPO_2
MPO
[1] ((dim=7|id=528|"Boson,Site,n=1")', (dim=7|id=528|"Boson,Site,n=1"), (dim=1|id=312|"Link,l=1"))
[2] ((dim=7|id=826|"Boson,Site,n=2")', (dim=7|id=826|"Boson,Site,n=2"), (dim=1|id=455|"Link,l=2"), (dim=1|id=312|"Link,l=1"))
[3] ((dim=7|id=755|"Boson,Site,n=3")', (dim=7|id=755|"Boson,Site,n=3"), (dim=1|id=455|"Link,l=2"))

But if I create myself the MPO I cannot

a
MPO
[1] ((dim=7|id=528|"Boson,Site,n=1"), (dim=7|id=528|"Boson,Site,n=1"), (dim=1|id=797|"Link,l=1"))
[2] ((dim=7|id=826|"Boson,Site,n=2"), (dim=7|id=826|"Boson,Site,n=2"), (dim=1|id=348|"Link,l=2"), (dim=1|id=797|"Link,l=1"))
[3] ((dim=7|id=755|"Boson,Site,n=3"), (dim=7|id=755|"Boson,Site,n=3"), (dim=1|id=348|"Link,l=2"))



a.data 
3-element Vector{ITensor}:
 ITensor ord=3
Dim 1: (dim=7|id=272|"Boson,Site,n=1")
Dim 2: (dim=7|id=272|"Boson,Site,n=1")
Dim 3: (dim=1|id=118|"Link,l=1")
NDTensors.Dense{ComplexF64, Vector{ComplexF64}}...

a+a



Trying to contract a tensor with indices:

((dim=7|id=528|"Boson,Site,n=1"), (dim=7|id=528|"Boson,Site,n=1"), (dim=1|id=797|"Link,l=1"))

and labels:

(-2, -4, 5)

with a combiner tensor with indices:

((dim=49|id=961|"Boson,Site,n=1"), (dim=7|id=528|"Boson,Site,n=1"), (dim=7|id=528|"Boson,Site,n=1"))

and labels:

(6, -3, -4).

This is not a valid combiner contraction.

If you are combining, the combined index of the combiner should be the only one uncontracted.

If you are uncombining, the combined index of the combiner should be the only one contracted.

By convention, the combined index should be the index in position 1 of the combiner tensor.


Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] invalid_combiner_contraction_error(tensor::NDTensors.DenseTensor{ComplexF64, 3, Tuple{Index{Int64}, Index{Int64}, Index{Int64}}, NDTensors.Dense{ComplexF64, Vector{ComplexF64}}}, tensor_labels::Tuple{Int64, Int64, Int64}, combiner_tensor::NDTensors.Tensor{Number, 3, NDTensors.Combiner, Tuple{Index{Int64}, Index{Int64}, Index{Int64}}}, combiner_tensor_labels::Tuple{Int64, Int64, Int64})
    @ NDTensors ~/.julia/packages/NDTensors/tAAFJ/src/combiner/combiner.jl:237
  [3] check_valid_combiner_contraction
    @ ~/.julia/packages/NDTensors/tAAFJ/src/combiner/combiner.jl:216 [inlined]
  [4] is_combining
    @ ~/.julia/packages/NDTensors/tAAFJ/src/combiner/combiner.jl:194 [inlined]
  [5] contract!!(output_tensor::NDTensors.DenseTensor{ComplexF64, 2, Tuple{Index{Int64}, Index{Int64}}, NDTensors.Dense{ComplexF64, Vector{ComplexF64}}}, output_tensor_labels::Tuple{Int64, Int64}, combiner_tensor::NDTensors.Tensor{Number, 3, NDTensors.Combiner, Tuple{Index{Int64}, Index{Int64}, Index{Int64}}}, combiner_tensor_labels::Tuple{Int64, Int64, Int64}, tensor::NDTensors.DenseTensor{ComplexF64, 3, Tuple{Index{Int64}, Index{Int64}, Index{Int64}}, NDTensors.Dense{ComplexF64, Vector{ComplexF64}}}, tensor_labels::Tuple{Int64, Int64, Int64})
    @ NDTensors ~/.julia/packages/NDTensors/tAAFJ/src/combiner/combiner.jl:109
  [6] contract!!
    @ ~/.julia/packages/NDTensors/tAAFJ/src/combiner/combiner.jl:156 [inlined]
  [7] contract(tensor1::NDTensors.DenseTensor{ComplexF64, 3, Tuple{Index{Int64}, Index{Int64}, Index{Int64}}, NDTensors.Dense{ComplexF64, Vector{ComplexF64}}}, labelstensor1::Tuple{Int64, Int64, Int64}, tensor2::NDTensors.Tensor{Number, 3, NDTensors.Combiner, Tuple{Index{Int64}, Index{Int64}, Index{Int64}}}, labelstensor2::Tuple{Int64, Int64, Int64}, labelsoutput_tensor::Tuple{Int64, Int64})
    @ NDTensors ~/.julia/packages/NDTensors/tAAFJ/src/generic_tensor_operations.jl:98
  [8] contract(::Type{NDTensors.CanContract{NDTensors.DenseTensor{ComplexF64, 3, Tuple{Index{Int64}, Index{Int64}, Index{Int64}}, NDTensors.Dense{ComplexF64, Vector{ComplexF64}}}, NDTensors.Tensor{Number, 3, NDTensors.Combiner, Tuple{Index{Int64}, Index{Int64}, Index{Int64}}}}}, tensor1::NDTensors.DenseTensor{ComplexF64, 3, Tuple{Index{Int64}, Index{Int64}, Index{Int64}}, NDTensors.Dense{ComplexF64, Vector{ComplexF64}}}, labels_tensor1::Tuple{Int64, Int64, Int64}, tensor2::NDTensors.Tensor{Number, 3, NDTensors.Combiner, Tuple{Index{Int64}, Index{Int64}, Index{Int64}}}, labels_tensor2::Tuple{Int64, Int64, Int64})
    @ NDTensors ~/.julia/packages/NDTensors/tAAFJ/src/generic_tensor_operations.jl:76
  [9] contract
    @ ~/.julia/packages/SimpleTraits/l1ZsK/src/SimpleTraits.jl:331 [inlined]
 [10] _contract(A::NDTensors.DenseTensor{ComplexF64, 3, Tuple{Index{Int64}, Index{Int64}, Index{Int64}}, NDTensors.Dense{ComplexF64, Vector{ComplexF64}}}, B::NDTensors.Tensor{Number, 3, NDTensors.Combiner, Tuple{Index{Int64}, Index{Int64}, Index{Int64}}})
    @ ITensors ~/.julia/packages/ITensors/GbZMQ/src/tensor_operations/tensor_algebra.jl:3
 [11] _contract(A::ITensor, B::ITensor)
    @ ITensors ~/.julia/packages/ITensors/GbZMQ/src/tensor_operations/tensor_algebra.jl:9
 [12] contract(A::ITensor, B::ITensor)
    @ ITensors ~/.julia/packages/ITensors/GbZMQ/src/tensor_operations/tensor_algebra.jl:104
 [13] *
    @ ~/.julia/packages/ITensors/GbZMQ/src/tensor_operations/tensor_algebra.jl:91 [inlined]
 [14] _broadcast_getindex_evalf
    @ ./broadcast.jl:670 [inlined]
 [15] _broadcast_getindex
    @ ./broadcast.jl:643 [inlined]
 [16] getindex
    @ ./broadcast.jl:597 [inlined]
 [17] macro expansion
    @ ./broadcast.jl:961 [inlined]
 [18] macro expansion
    @ ./simdloop.jl:77 [inlined]
 [19] copyto!
    @ ./broadcast.jl:960 [inlined]
 [20] copyto!
    @ ./broadcast.jl:913 [inlined]
 [21] copyto!(ψ::MPO, b::Base.Broadcast.Broadcasted{Base.Broadcast.Style{MPO}, Tuple{Base.OneTo{Int64}}, typeof(*), Tuple{MPO, Vector{ITensor}}})
    @ ITensors ~/.julia/packages/ITensors/GbZMQ/src/mps/abstractmps.jl:2354
 [22] materialize!
    @ ./broadcast.jl:871 [inlined]
 [23] materialize!(dest::MPO, bc::Base.Broadcast.Broadcasted{Base.Broadcast.Style{MPO}, Nothing, typeof(*), Tuple{MPO, Vector{ITensor}}})
    @ Base.Broadcast ./broadcast.jl:868
 [24] +(::NDTensors.Algorithm{:densitymatrix}, ::MPO, ::Vararg{MPO}; cutoff::Float64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ITensors ~/.julia/packages/ITensors/GbZMQ/src/mps/abstractmps.jl:205
 [25] +(::NDTensors.Algorithm{:densitymatrix}, ::MPO, ::MPO)
    @ ITensors ~/.julia/packages/ITensors/GbZMQ/src/mps/abstractmps.jl:1430
 [26] +(::MPO, ::Vararg{MPO}; alg::NDTensors.Algorithm{:densitymatrix}, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ITensors ~/.julia/packages/ITensors/GbZMQ/src/mps/abstractmps.jl:1536
 [27] +(::MPO, ::MPO)
    @ ITensors ~/.julia/packages/ITensors/GbZMQ/src/mps/abstractmps.jl:1535
 [28] top-level scope
    @ In[61]:1

And I have no clue of why

That is how I define the link indeces

l_indeces = [Index(1, "Link,l=$ii") for ii in 1:3-1]

and here how I create the MPO

for m in 1:3
    if m==1
        H_MPO[m] = ITensor(mats[m],sites[m],sites[m],l_indeces[m])
    elseif m ==n_modes
        H_MPO[m] = ITensor(mats[m],sites[m],sites[m],l_indeces[m-1])
    else
        H_MPO[m] = ITensor(mats[m],sites[m],sites[m],l_indeces[m],l_indeces[m-1])
    end
end

where mats is a tuple with all the matrices

phx · March 16, 2023, 5:18pm

I think there’s something wrong in the indices of the MPO you defined. Shouldn’t some indices be primed? I would write

H_MPO[m] = ITensor(mats[m], sites[m], sites[m]', l_indeces[m])

and so on in your loop, otherwise you end up with two identical indices in the ITensor object.

EmilioRui · March 16, 2023, 5:30pm

I surely did not,
Could you guide me to some documentation to understand what ‘priming’ does?

Thank you so much!

EmilioRui · March 16, 2023, 5:53pm

By the way it works.
Thank you so so much

miles · March 16, 2023, 5:55pm

Thanks @phx for your helpful answers.

@EmilioRui your questions are a good reminder for us that we need to improve the documentation on adding MPS and MPO, like at least showing some more examples and things and also making sure the add and sum functions appear in the docs.

@EmilioRui also for the philosophy of priming, a good resource is the ITensor Paper
https://www.scipost.org/SciPostPhysCodeb.4
Priming is discussed on page 8 and in other places.

In a nutshell, the main feature of ITensor is that identical indices automatically contract (and you are therefore free to not think about the ordering of indices; the interface is mostly independent of this choice). As a result, you are not supposed to have two copies of the same Index on an ITensor because otherwise the automatic contraction feature wouldn’t know which one you want to contract.

But of course sometimes you need essentially the same index to appear twice, like e.g. the “bra” version of the index and the “ket” version. Or for other reasons. So priming is a way of “disambiguating” multiple copies of the same index in a lightweight way. You can raise and lower the “prime level” (number of primes) of an index. Two indices with otherwise identical data but different prime levels are treated as distinct by the ITensor system.

Finally, and importantly, our convention for operators throughout the ITensor ecosystem is that operators carry two copies of their physical indices: one copy has no primes and the other has one prime.

EmilioRui · March 24, 2023, 8:53am

Thank you so much!

Topic		Replies	Views
Stack multiple MPOs on an MPS ITensor Julia Questions julia , mpo	2	146	February 19, 2024
Adding MPO's, defined on different sites ITensor Julia Questions	4	286	September 18, 2023
Bug in multiplying MPOs and tracing MPOs ITensor Julia Questions julia , mpo	2	331	October 23, 2023
Implementing the exponential of a sum of local operators as an MPO ITensor Julia Questions	11	83	September 28, 2024
Manually constructing MPO from set of ITensors ITensor Julia Questions julia , dmrg , mpo	2	137	July 15, 2024

Sum of two MPOs

Related topics