Unitary operators on evolved state

Question

I am trying to act some state with 2-layered unitaries. For instance, consider the following code:

function ITensors.op(::OpName"evoH", ::SiteType"S=1/2", s1::Index, s2::Index)
     mat = [1. 0. 0. 0.
             0. 1. 0. 0.
             0. 0. 0. 1.0*im
             0. 0. 1.0*im 0.]

     return itensor(mat, s2', s1', s2, s1)
end

function evolution_single_step(sites, psi)
    N = length(psi)
    os = [("evoH", 1, 2)]
    for site_os=3:2:N-1
        push!(os, ("evoH", site_os, site_os+1))
    end
    for site_es=2:2:N-1
        push!(os, ("evoH", site_es, site_es+1))
    end
    psi = apply(ops(os, sites), psi; cutoff = 1e-15)
    return psi
end

sites0 = siteinds("S=1/2", 10) 
state = ["Dn" for n in 1:10]
psi_init = productMPS(ComplexF64, sites0, state)

psi_init_dp = deepcopy(psi_init)
evolution_single_step(sites0, psi_init_dp)
println(inner(psi_init_dp, psi_init_dp))

As expected the inner product is 1.

But when I try to recreate the above using the following code:

function two_qubit_gate(psi, i, j)
      sind_i = siteind(psi, I)
      sind_j = siteind(psi, j)

      # construct the unitary
      tq_gate = ITensor(sind_i', sind_i, sind_j', sind_j)
      tq_gate = complex(tq_gate)

      tq_gate.store.data[1] = 1.
      tq_gate.store.data[4] = 1.
      tq_gate.store.data[14] = 1.0*im
      tq_gate.store.data[15] = 1.0*im

      return tq_gate
end


function evolution_single_step(psi)
      N = length(psi)
      for site_os=1:2:N-1
             orthogonalize!(psi, site_os)
             op_psi1 = psi[site_os]*psi[site_os+1]*two_qubit_gate(psi, site_os, site_os+1)
             noprime!(op_psi1)
             U,S,V = svd(op_psi1, siteind(psi, site_os))
             psi[site_os] = U
             psi[site_os+1] = S*V
             println(inner(psi, psi))
      end

      for site_es=2:2:N-1
             orthogonalize!(psi, site_es)
             op_psi2 = psi[site_es]*psi[site_es+1]*two_qubit_gate(psi, site_es, site_es+1)
             noprime!(op_psi2)
             U,S,V = svd(op_psi2, siteind(psi, site_es))
             psi[site_es] = U
             psi[site_es+1] = S*V
             println(inner(psi, psi))
      end
      return psi
end

sites0 = siteinds("S=1/2", 10)
state = ["Dn" for n in 1:10]
psi_init = productMPS(ComplexF64, sites0, state)

evolution_single_step(psi_init)
println(inner(psi_init, psi_init))

For the first layer (for the first for loop) the norm is conserved, it goes to 1 while for the action on the evolved state, the norm remains unconserved (grows on as an exponential, resulting in 2, 4, 8, 16). I would have assumed that the machinery of apply is prototyped as in the second function. It would be helpful to know if my comparison is right and if so why the action on the evolved state is resulting in an unconserved norm. Thanks!

Answer 1

Hi Amit,
It's hard to immediately guess what it causing the error you're having, but I agree the result shouldn't be as you say.

Here are two possible issues I spotted right away (with #2 being more likely the culprit, though #1 is also important for improving your code):

1. we would not recommend that you set elements of ITensors by accessing the storage directly, as you did in your second panel of code inside the function twoqubitgate. I'd recommend either using the pattern you used in your "evoH" gate where you initialize an ITensor with a Julia matrix or array, or else setting the elements using the indices of the ITensor such through code like:

   tqgate[sindi'=>1,sindi=>1,sindj'=>1,sind_j=>1] = 1.0

and similar. Accessing the storage is not part of the interface of an ITensor and is error prone and subject to arbitrary future breaking changes.

1. I noticed that in your SVD'd after applying a gate, you only put one of the site indices as the inded to go onto U. But once the MPS develops non-trivial link indices as the result of entanglement generated between sites, you would need to also include the appropriate link or bond index in the list of indices that should go onto U. Otherwise your algorithm will be incorrect and will not be equivalent to the standard "TEBD" algorithm for applying gates to an MPS.

Happy to discuss more -

Miles

Comments

Hi Miles,

Thanks for the quick response and the recommended usage. Regarding point 2, I was following the guide on applying two-qubit operators to a MPS as in the following reference: https://itensor.org/docs.cgi?vers=julia&page=formulas/mps_twosite_op. I started to look into the contents of the MPS obtained after the first `for` loop and I noticed that there were additional indices that were not contracted, though they share the same index (with the string: "Link,u", which has been mentioned in the docs of the SVD). I was of the thinking that these additional bonds between the indices would be automatically contracted as they share the same index. My apologies if I missed something trivial in the above, but any further pointers on how to improve/make the second snippet work would be helpful. Thanks!

---

Hi Amit,
Yes I can see why that example was a bit misleading, because there I got the indices to put into the SVD in sort of a "tricky" way. I've just improved the example to get the appropriate indices in a more generally-applicable way and given more explanation about it:

http://itensor.org/docs.cgi?vers=julia&page=formulas/mps_twosite_op

So I think if you use the same pattern of calling `uniqueinds` to get the indices to put into the SVD, then it should now work the way you want. It is crucial to have both the left site index and left-most link index end up on the U tensor so that U and (S*V) or else (U*S) and V are appropriate MPS tensors afterward.

Best,
Miles

---

Hi Miles,

Thanks for the solution. It seems to preserve the norm now :) but I have to test it further. Many thanks!

---

Great! And in addition to just checking that the norm is preserved, I would encourage you to print out some of the tensors or the whole MPS to confirm that the correct index structure of the tensors is being preserved too.

Also, I wanted to mention that if you do need to access the elements of an ITensor using “linear indexing” i.e. treating it as a flat array, the proper interface to use is like this:

    A = ITensor(i,j)
    A[1] = 0.1
    A[2] = -0.2
    # etc.

But for setting a small number of elements I would still suggest the index=>value approach I mentioned above. Or I would suggest making an array first and constructing the ITensor from the array.

Miles