# Create MPO using matrix description of observable?

Hello, I am a beginner to ITensor so this may be a basic question, but I did not find the answer to my question in the documentation or this forum.

I have a state psi, and I would like to compute the overlap <psi|rho|psi> where rho is a density matrix. psi can be represented as an MPS, and rho will have low entanglement, so I would like to convert it into an MPO to efficiently compute this overlap, however I am not sure how to do this conversion.

The form of my density matrix rho has the exact form of (S44) in https://arxiv.org/pdf/2002.08953.pdf (pg 25) i.e. a tensor product of local operators. How can I create a MPO to represent this? Can I take the matrix representation and convert it into a MPO? Or is there an alternate method more natural for ITensor.

Thanks for the question. We do have a system for making this kind of MPO, which is called the OpSum system. Most of the code examples in the documentation show it being used to make MPOs for local Hamiltonians, not density matrices, but if you density matrix is a sum of products of local operators, then this is the very thing OpSum is for.

Please look at this example where OpSum is being used:
https://itensor.github.io/ITensors.jl/stable/examples/DMRG.html

Hi Miles,

Thank you for your answer. It seems like with OpSum I could do this however I would have to first expand the state into products of local operators. Is there a way I can avoid this by exploiting the tensor product structure of the operator? It seems like this should be simple to represent as a bond dimension 1 MPO due to this

I see. I had thought in your question that “tensor product structure” or the operators and expansion in terms of local operators would be the same thing. Is there a reason they are different here? Or what is mean by tensor product structure? (Could you use Latex to write the math for the structure you wish to capture?)

Thanks -

My target MPO has the form \otimes_j (3 U_{j}^{\dagger} |b_j><b_j| U_j - \mathbb{I}_j), where |b_j><b_j| are projectors onto the |0> or |1> state, U_j are rotations into the X, Y or Z basis ,and \mathbb{I}_j is the identity on the jth qubit. Seeing as the state can be locally described on each qubit, is there a way of separately constructing each 3 U_{j}^{\dagger} |b_j><b_j| U_j - \mathbb{I} into a 1 qubit MPO, then tensoring them together into the final MPO?

Thanks for the extra details.

So actually you’re in luck and this is the simplest kind of MPO to construct. You could still do it using OpSum but it would be overkill in this case. All you have to do in this case is the following. Make an MPO for your N qubits and a collection of site indices you are going to use

N = 10
s = siteinds("Qubit",N)
M = MPO(s)


then just set the MPO tensors one by one to the values you want, following the ITensor convention that the j’th MPO tensor has the indices s[j], s[j]'.

So like this:

for j=1:N
T = ITensor(s[j],s[j]')
T[1,1] = ...
T[2,1] = ...
T[1,2] = ...
T[2,2] = ...
M[j] = T
end


You could either work out the matrix elements of T ahead of time and then set them one by one like I am suggesting in the code snippet above. Or you could use some built in ITensor facilities for making each MPO tensor. What I mean is, say you want to make the operator

3|x+\rangle\langle x\!+\!| - I

You could do this with the code

T = 3*state("X+",s[j]')*state("X+",s[j]) - op("I",s[j])


Here I haven’t included the U’s explicity, but based on what you said they are Pauli matrices so you could work out what U^\dagger_j |b_j\rangle is in advance without numerically acting U^\dagger_j onto a state.

You can also define your own custom overloads of the state and op functions if you wish: