Playing with MPS/ ITensors Indexing

+1 vote

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?

commented by (53.8k points)
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.

selected by

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

commented by (630 points)
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
commented by (53.8k points)
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.
commented by (53.8k points)
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
commented by (630 points)
edited by
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.
commented by (53.8k points)
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.
commented by (630 points)
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)
commented by (53.8k points)
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.