Matrix element of (tensor products of) Pauli operators in Julia

Question

Hi,

Is there an easy way of calculating 2^N matrix elements like:

A(i1,i2,...,iN)= < phi| X^(i1)⊗X^(i2),...,⊗X^(iN))|psi> ?

where X denotes the spin 1/2 Pauli X matrix and I get the 1D MPS phi and psi from a DMRG optimization routine. Also note i is a N-component binary vector.

Thanks in advance,
Arnab

Answer 1

Hi Arnab, yes this is a task that MPS (and ITensor) is perfectly suited for.

The simplest way to do this is to:
1. choose a target string i1, i2, i3, that you want to compute (i.e. fix the i’s)
2. apply the X^i_n operators to the MPS psi. Since these are single-site operators they can be applied very quickly and efficiently (even in parallel if you want to go there) by just multiplying each MPS by X on its physical index. In ITensor, you obtain the X operator for that site, use the * operator to contract it with the MPS, then call noprime to set the site’s prime level back to 0. Of course for sites where you act with the identity, you can just skip over such sites.
3. finally take the resulting MPS, calling it Xpsi, say, then call inner(phi,Xpsi) to obtain the overlap of the modified psi with phi.

Then you can repeat this for other i-strings.

Finally, if you are planning to obtain a large number, or all amplitudes this way, there should be a way to organize them where say you reuse substrings, i.e. overlap by hand the modified psi tensors with phi tensors and then attach the next ones in both combinations and so on. It would take a bit of thought to see how advantage there would be in doing this (because ultimately you’d still have to make exponentially many such “substring partial overlap” tensors anyway).

Let me know if any of those steps aren’t clear.

Oh also we have a new framework precisely for making it easier to apply various operators to MPS, but it’s not totally well documented yet. Here is a sample code file showing it off:

https://github.com/ITensor/ITensors.jl/blob/master/examples/gate_evolution/quantum_simulator.jl

But it’s overkill just for your case of single-site operators which you can definitely do with the usual ITensor interface.

Best,
Miles

Comments

Hi Miles,

Thank you for your answer. Actually my initial code does what you say upstairs. But the issue is I am indeed extracting all the amplitudes (all the 2^N of them) this way and thus it would be nice to not write a line of code for each fixed value of i.

Best,
Arnab

---

Hi Arnab,
Yes let's discuss this some more. So in the version I am thinking, you would write one line of code (or short block of code) not for *every* i, but which would work generically for any fixed, given i. Then you would put a "for" loop around this which loops over all possible values of i, and store the results in an array.

For example, inside this block of code (within the for loop), as you reach each site, the code would ask "if i_n == 1" then if it is true, it would apply an X on site n. Otherwise it would go onto the next site. So there would be an inner loop from i=1:N with such an if statement inside.

Is that what you already had in mind or was I misunderstanding the question?

---

In terms of how to loop over all bit strings in Julia, there are a few different ways if you don't already know them (just thought I would share what I know). You'll have to double check that these work properly:

(1) This bit of code:

  bits = [1:2 for n in 1:N]
  its = NTuple{N,Int}[]
  for it in Iterators.product(bits...)
  # ... code goes here ...
  end

(2) or the solutions in this forum post:
https://discourse.julialang.org/t/is-there-an-idiomatic-way-to-iterate-over-all-binary-strings-of-a-given-length/35487/11

---

Hi Miles,

Thank you for the additional comments. I am indeed not very familiar with Julia as a language. I will try the looping over bit string idea to see if that resolves the problem.

---

Hi Miles,

A quick follow-up to this: from what you outline in your answer, what exactly does the "noprime" function play here and how do you apply it? Is it possible to refer me to some documentation regarding this or just a sample line of code that applies X⊗I⊗X⊗X⊗I⊗I to a 6-qubit MPS state |psi> in Julia?

---

Hi Arnab, so ITensor indices (Index objects) have a "prime level" which can be used to distinguish them from other indices that would otherwise compare equal to them.

The noprime function, when called on an ITensor, returns a copy of that ITensor with all of the primelevels of the indices set to zero.

There is more information about all this in the following page of the ITensor book:
http://itensor.org/docs.cgi?vers=julia&page=book/index

Also here is the documentation for the `noprime` function (at least the version taking an Index, and there are similar versions taking an ITensor, an MPS, etc.):
https://itensor.github.io/ITensors.jl/stable/IndexType.html#ITensors.noprime-Tuple{Index}

To access this kind of documentation, you can also do: ?noprime
inside of the interactive mode of Julia.

Best,
Miles

---

Here is what doing ?noprime inside of Julia gives:

help?> noprime
search: noprime noprime!

  noprime(i::Index)

  Return a copy of Index i with its prime level set to zero.

  ───────────────────────────────────────────────────────────────────────────────────────────────────

  noprime[!](A::ITensor; <keyword arguments>) -> ITensor
  
  noprime(is::IndexSet; <keyword arguments>) -> IndexSet

  Set the prime level of the indices of an ITensor or IndexSet to zero.

  Optionally, only modify the indices with the specified keyword arguments.

...

---

Hi Miles,

Sorry for these lingering questions. I tried another method for terms involving pauli operators that have support at more than one site:

1. First define a MPO using the AutoMPO function :
N=8
sites = siteinds("S=1/2",N)
sample = AutoMPO()
sample+= 4,"Sx",4,"Sx",5;
Operator=MPO(ampo,sites);

2. Use Inner(phi,Operator,psi) to get the matrix element for X_4⊗X_5.

However this outputs nonsensical results like -7.968 for normalized states psi/phi. Can you point me to why this method doesn't work?

---

Hi Arnab,
I'm not sure why you are getting that result. The only thing that looks 'off' about your code to me is that the line: Operator = MPO(ampo,sites); seems to be taking "ampo" as input rather than "sample". So did you define another AutoMPO object earlier in the code called "ampo"?

---

When I run this code in Julia I get the following result:

julia> let
       N=8
       sites = siteinds("S=1/2",N)
       sample = AutoMPO()
       sample+= 4,"Sx",4,"Sx",5;
       Operator=MPO(sample,sites)
       A = randomMPS(sites,4); B = randomMPS(sites,4)
       inner(A,Operator,B)
       end

Output:
-0.06913588751318667

---

Oops. I am so stupid. Sorry about that. Thanks for pointing me to that.

---

Not stupid - these kind of bugs happen!