Playing with MPS/ ITensors Indexing

Question

Hi,

I have a couple of questions which relate to indexing in ITensor.

Firstly, I have an MPS in a spin-1/2 1D chain denoted by T. I can get the local 2-D tensors (at each) by doing setting the site-index to 1 or 2. I am pretty sure about this but can you confirm that this would give me the components in the computational (Sz) basis ? Now, I would like to change modify the tensor so that I can get this in the (Sx) basis. This basically amounts to adding and subtracting the local 2-D tensors. However, I am not sure how to put them back in the original tensor without messing up the indexing. Here is some minimal code:

using ITensors
N=20
sites = siteinds("S=1/2",N)
bond_dim=24
T = randomMPS(sites,bond_dim);
for j=1:length(T)
T0=T[j]*setelt(sites[j]=>1)
T1=T[j]*setelt(sites[j]=>2)
Tx0=(T0+T1)/sqrt(2)
Tx1=(T0-T1)/sqrt(2)
#I would like to assign Tx0 to something like T[j]*setelt(sites[j]=>1) but how?
end

Finally, I have been recently doing a lot of calculations involving contraction of ITensor objects and I was wondering if there is some in-built function that returns the different index of any ITensor object ( somewhat like how "linkind" works for MPS objects ).

In general, is the handling of ITensor indexing in Julia well documented somewhere?

Comments

Hi Arnab,
Regarding ITensor indexing, it is documented here:
https://itensor.github.io/ITensors.jl/stable/ITensorType.html#Getting-and-setting-elements-1

Please let us know if that part of the documentation leaves out some information you were looking for.

Answer 1

Hi Arnab,
Thanks for the question. I'll give some brief answers and we can discuss more in comments below.

(1) "I can get the local 2-D tensors (at each) by doing setting the site-index to 1 or 2." This is not correct unless the MPS is a product (unentangled) state. I think what you are trying to do is to obtain the probability of finding some site j to be in the up state versus the down state. So you are looking for these two numbers, call them pup and pdn. These numbers are the diagonal elements of the reduced density matrix for site j. This reduced density matrix can be found by (a) first calling orthogonalize!(T,j) to make site j the orthogonality center of your MPS then (b) by creating the density matrix rhoj = T[j] * prime(dag(T[j]),"Site"). This density matrix will be a 2x2 matrix and its diagonal entries will be pup and p_dn. Please check that they add to 1.0.

(2) to obtain the probabilities of up and down states in a different basis, such as the x basis, you can perform a change of basis by acting with a unitary operator on site j. First, make sure your MPS is orthogonalized to site j as in (1a) above. Then obtain the j'th MPS tensor: Aj = T[j]. Now act with the Hadamard gate which will change the basis to the X basis: Ajx = noprime(op("H",sites,j)*Aj). Finally, make the density matrix in the X basis: rhox_j = Ajx * prime(dag(Ajx),"Site"). The diagonal entries of rhox_j will be the probabilities for up and down in the X basis.

Best,
Miles

Comments

Hi Miles,

Thank you for your detailed answer. I do agree with you that in general if I make a measurement at a site then I need to compute the reduced density matrix to get the corresponding probabilities. However, I am basically running a quantum wire protocol in measurement based model where one can just get away with dealing the tensor components in a specified basis (upto post processing).

(2) About applying the Hadamard: yes, I agree that this can be done to convert from the Z to X basis. However, for converting to an arbitrary basis (which I need for my setting) I would need to define an arbitrary local operator. Is this possible in ITensor?

Best,
Arnab

---

Hi Arnab, regarding the second part of your question, yes you can definitely create an arbitrary tensor in ITensor. To do so, construct an ITensor with the indices you want it to have (so for an operator it will have indices i and I’ where i=sites[n] is one of the physical indices of your MPS), then set the elements to the values you want them to have. There are many examples of how to set tensor elements on the website and in the docs, but the basic syntax is:

T[i=>a,i’=>b] = x

Where here one is setting the a,b element to the value x.

---

If you want to get slightly fancier, and make other operators from strings similar to how you can make a Hadamard from the string “H”, we have a tutorial about this:
http://itensor.org/docs.cgi?vers=julia&page=formulas/sitetype_extending

---

Hi Miles,

Thank you for your detailed answer. I tried some things out and most of them worked. However, I am still a little struck on the following:

I would like to contract a part of a MPS "T" (for example say from site 1 to site k)  after having projected the physical leg onto the computational basis . The most tedious way of doing this is to do: T[1]*setelt(sites[1]=>(arr[1]+1))*T[2]*setelt(sites[2]=>(arr[2]+1))*...*T[k]*setelt(sites[k]=>(arr[k]+1)) where I write down (almost) the same term k times ( "arr" contains the measurement outcomes) . However, in my case k is around 150 and I was wondering if there is an option in ITensor to just automate this.

Edit: Also, are there subtleties involved when a local operator has complex entries? My script throws up "InexactError" whenever I apply such a newly  defined operator to a MPS but everything works perfectly when the imaginary part is set to zero.

---

Hi Arnab, so the solution to your issue is to use a for loop. Something roughly like this:

R = ITensor(1)
for n=1:k
R = R*(T[n]*setelt(sites[n]=>(arr[n]+1)))
end

@show inds(R)

Regarding the InexactError, I might need to see the back trace and the line of code causing it to help, but it could be related to trying to convert a complex number to a real number in some way that would result in losing the non-zero imaginary part, versus actually just asking for the real part explicitly.

---

Hi Miles, thank you so much for the loop solution! Here is the local operator I am trying to define:

using ITensors
beta=1
function ITensors.op!(Op::ITensor,
                      ::OpName"Rz",
                      ::SiteType"S=1/2",
                      s::Index)
  Op[s'=>1,s=>1] = cos(beta)-sin(beta)*im
  Op[s'=>2,s=>2] = cos(beta)+sin(beta)*im
  end
But then when I try to apply it:

using ITensors
N=6
bond_dim=24
sites = siteinds("S=1/2",N)
T = randomMPS(sites,bond_dim)
@show T
# apply Rz on site 3
s = siteind(T,3)
newpsi= op(s,"Rz")*T[3]
noprime!(newpsi)
T[3]= newpsi
@show T

It throws up an error--  InexactError: Float64(0.5403023058681398 -0.8414709848078965im)

---

Hi Arnab, good question. So in the current version of ITensor you aren't allowed to set an element of an ITensor to a complex value if that ITensor was constructed having real storage internally. (We are thinking of lifting this restriction in the future.)

For now, you can fix this two ways:

1) in the first line of your op! function, you can put the line complex!(Op) which will switch the storage data in Op to be complex-valued, then you can successfully set your elements

2) you can use the other style of op function, which is op without the exclamation "!" mark. This version lets you construct the Op tensor any way you want and then you return it from that function. So like this:

function ITensors.op(Op::ITensor,
                      ::OpName"Rz",
                      ::SiteType"S=1/2",
                      s::Index)
  Op = ITensor(ComplexF64,s',dag(s))
  Op[s'=>1,s=>1] = cos(beta)-sin(beta)*im
  Op[s'=>2,s=>2] = cos(beta)+sin(beta)*im
  return Op
  end

As you can see in that version, when I constructed Op I made its element type ComplexF64 from the beginning. So that version (2) is a bit faster than (1) which makes the elements to be real zeros then has to remake them to be complex zeros.

---

Hello,
I am having a similar issue, so I thought I would just add to this thread! I have tried to use the complex!() function when I define my complex op:

ITensors.op(::OpName"YRot", ::SiteType"S=1/2", s::Index) = complex!((sqrt(2))*(op("Sy", s) + op("Sz", s)))

But I also get an InExactError thrown when I try and define and MPO with the op:

YRotMPO = MPO(sites, "YRot")

InexactError: Float64(0.0 - 0.7071067811865475im)

I have also tried defining a function, following the example above:

function ITensors.op(Op::ITensor, ::OpName"YRotfunc", ::SiteType"S=1/2", s::Index)
    Op = ITensor(ComplexF64,s',dag(s))
    Op[s'=>1,s=>1] = 1/sqrt(2)
    Op[s'=>2,s=>1] = -im/sqrt(2)
    Op[s'=>1,s=>2] = im/sqrt(2)
    Op[s'=>2,s=>2] = 1/sqrt(2)
    return Op
end

But this also throws the same error when I try and construct an MPO from it:

YRotMPOfunc = MPO(sites, "YRotfunc")

InexactError: Float64(0.0 - 0.7071067811865475im)

Any advice would be much appreciated! I am running v0.1 of ITensor -- I have been using Julia Pro, which for some reason only installs the v0.1 version by default.

---

Hi, so this issue is fixed as of the latest version of ITensor, version 0.2.3. Could you please update to that version?

Also, in your longer version of `op`, there should not be an ITensor `Op` as the first argument. This version ought to work:

function ITensors.op( ::OpName"YRotfunc", ::SiteType"S=1/2", s::Index)
    Op = ITensor(ComplexF64,s',dag(s))
    Op[s'=>1,s=>1] = 1/sqrt(2)
    Op[s'=>2,s=>1] = -im/sqrt(2)
    Op[s'=>1,s=>2] = im/sqrt(2)
    Op[s'=>2,s=>2] = 1/sqrt(2)
    return Op
end

---

I apologize if this is a naive question -- I'm new to Julia and ITensors -- but do you know if version 0.2.3 is compatible with Julia Pro? When I use the normal command:

pkg> update ITensors

no update occurs; the package remains at version 0.1.21.

---

Based on the article linked below, it sounds like JuliaPro used to block certain package updates but that the latest versions of JuliPro (as of August 2020) now let you use any package versions you want.

However, there can be other reasons a package doesn't update, mainly that you have other packages which themselves aren't updated so will block the update of a certain package. So the thing to try is just doing "update" which should update all packages to their latest versions, if I remember correctly. Please be sure you want this to happen, though, and that it might take some minutes to complete.

https://juliacomputing.com/blog/2020/08/juliapro-v15/

---

Thank you so much for all your help! I had to go from Julia Pro v. 1.5.1 to v. 1.5.4, and now installing ITensors automatically gives v. 2.3. Hopefully everything will fall into place now!

---

Great! Glad that was all it took