MPO construction from two site custom op

ITensor Julia Questions
1

Hi,
I want to build an MPO using a custom defined ITensor operator. The MPO can be constructed using
MPO(::sites, ::list of op) (reference to definition : https://itensor.github.io/ITensors.jl/dev/MPSandMPO.html).
However, the caveat is that the list of operators are two-body. For example, for spin half systems, this two body operator could be Z \otimes Z or the one shown below.

using ITensors
using LinearAlgebra

sites = siteinds("Qubit", 6)

function ITensors.op(::OpName"TwoSiteOp", ::SiteType"Qubit"; θ1, θ2)
    
    M = cos(θ1)*Diagonal(ones(4)) - im*sin(θ2)*[1 0 0 0; 0 -1 0 0; 0 0 -1 0; 0 0 0 1]

    return M
end

Now, the question is how can one construct an MPO for the defined ‘siteinds’, sites?
The op list will be made of various operators with different \theta 1 and \theta 2.

Any suggestion for this will be very helpful.

2

See Miles’ answer from here:

Essentially the OpSum system can’t do that at the moment, however if you can write it as a string (O_i \otimes O_j \dots ) or a sum of strings as in the post then that can be used. And for generic 2 site matrices you would have to decompose it by hand into that format

3

Hi Ryan,
Thank you very much for the reply. To clarify, the two-site operator I mentioned above, e^{-i \theta_1 \mathbb Z_{1} \otimes \mathbb Z_{2}}, can be decomposed as a product of single-site tensors as
\begin{pmatrix} \cos(\theta_1) \mathbb{I}_{1} &-i\sin(\theta_1) \mathbb Z_{1} \end{pmatrix} \otimes \begin{pmatrix} \mathbb{I}_{2} \\ \mathbb Z_{2} \end{pmatrix}.
Is there a way to construct an MPO by using such tensors as input?

4

If you have the operator O = \cos(\theta_1) I_1 \otimes I_2 -i\sin(\theta_1) Z_1 \otimes Z_2 then the corresponding OpSum call would be

os = OpSum()
os += cos(θ1), "I",1,"I",2      # the second I is not needed really
os += -1im*sin(θ2),"Z",1,"Z",2
O = MPO(os,sites)
5

Hi Ryan,
thank you so much! There might still be a miscommunication though. Just to be sure: I need to build an MPO for the operator
e^{-i\theta\mathbb Z_{1} \otimes \mathbb Z_{2}} \otimes e^{-i\theta\mathbb Z_{3} \otimes \mathbb Z_{4}} \otimes ....... \otimes e^{-i\theta\mathbb Z_{N-1} \otimes \mathbb Z_{N}} .
This is a product rather than sum of two-site operators (unlike the link with Miles’ answer).
Taking your suggestion to use OpSum, I tried the following code :

using ITensors

N = 8
sites = siteinds("Qubit", N)
U = MPO(sites)

θ = 0.5
os = OpSum()
os += cos(θ), "I", 1, "I", 2
os += -1im * sin(θ), "Z", 1, "Z", 2

for j in 1:2:N-1
    U[j:j+1] = MPO(os, sites[j:j+1])
end

However, this does not seem to produce the correct results. Could you help me correct it?

6

You’ll need to decompose your operator into a list of single site operators
O = \sum_n O_1^n \otimes O_2^n \otimes \dots, which for each n you would write the os+= line as above. You’d then call MPO just once to form the MPO.

Then one trick you can do, if you want say all "Z"s like A*Z_1 \otimes Z_2 \otimes \dots \otimes Z_N where A is a constant, then you can write

A = 0.1 # etc
term = [("Z",i) for i=1:N]
add!(os, A, (term...)...)
7

Thank you very much.

1 Like
8

If you expand everything out in terms of single site operators as @ryanlevy suggested then you’ll have 2^{N / 2} different products which won’t scale well. I think you, @Bhargava_Balaganchi, had the right idea in your code example but didn’t set the elements of U correctly.

using ITensors

N = 8
sites = siteinds("Qubit", N)
U = MPO(sites, "I")

θ = 0.5
os = OpSum()
os += cos(θ), "I", 1, "I", 1
os += -1im * sin(θ), "Z", 1, "Z", 2

for j in 1:2:N-1
  Ui = MPO_new(os, sites[j:(j + 1)])
  U[j] = Ui[1]
  U[j + 1] = Ui[2]
end

# The purpose of orthogonalize! is to insert the missing link dims of size 1. 
orthogonalize!(U, 1)
@show linkdims(U)

linkdims(U) = [2, 1, 2, 1, 2, 1, 2]

Alternatively, a more general approach is to construct MPOs for the individual terms and then multiply them together. For this you’ll want to use the ITensorMPOConstruction library since ITensors’ implementation adds extra links for the starting and ending identity operator which can really add up when you multiply a bunch of MPOs together.

using ITensors
using ITensorMPOConstruction

N = 8
θ = 0.5

sites = siteinds("Qubit", N)
U = MPO(sites, "I")

for j in 1:2:N-1
  os = OpSum()
  os += cos(θ), "I", j, "I", j + 1
  os += -1im * sin(θ), "Z", j, "Z", j + 1
  Ui = MPO_new(os, sites)
  U = apply(Ui, U; alg="naive", truncate=false)
end

@show linkdims(U)

linkdims(U) = [2, 1, 2, 1, 2, 1, 2]

If you were to use ITensors’ MPO construction algorithm you would get

linkdims(U) = [32, 16, 32, 16, 32, 16, 32]

and in general find the the link dimension grows exponentially with N.