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,

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!