Computing expectation value for a given observable after DMRG simulation

Question

I am currently interested in computing the quantity

Tr(ρA)

where A is some pre-determined matrix with the same dimensions as the density matrix rho. Rho is given as the output from a DMRG simulation. However, I am unsure how to take the output psi (i.e., the updated value of randomMPS(sites) after DMRG) and transform it into a matrix to perform the appropriate multiplication with A. I assume there is a built-in feature for doing this, but as I am new to iTensor I am unsure how to do this. Also, please note that I am only using the Julia version.

EDIT: Also, I know it's not advantageous to work in array format, but this is just for benchmarking purposes for small systems.

Answer 1

Thanks for the question. There isn't exactly a built-in feature for what you are asking, but it's straightforward to do. I made a sample code below that hopefully show you the steps you need.

To transform psi into a "vector" (tensor with N indices in the full Hilbert space), you can contract all of the MPS tensors into a single tensor using either a for loop or the prod function (prod does something similar to the capital-Pi symbol in mathematics).

From what I understand of your question, there should be no need to really form the entire density matrix from psi. This is because it's a pure state so just an outer product of psi with itself and therefore much more efficient to deal with each copy of psi separately.

In the code below, after transforming psi into a single tensor psiT, I make another (arbitrary) tensor AT which has two copies of the "site" indices: one set with prime level 0 and the other with prime level 1. Performing the tensor contraction using the ITensor * operation contracts over the unprimed indices, and then the noprime function removes the primes from the remaining indices.

In the end, the resulting density matrix is formally defined as this new state Apsi outer producted with your original state psi (being sloppy here about conjugation, but that's easy to include if dealing with complex numbers or QN-conserving tensors).

Please reply below if you have any follow-up questions.

using ITensors

let
  N = 10
  chi = 2
  s = siteinds("S=1/2",N)

  # Make a random MPS with sites s
  # and bond dimension chi:
  psi = randomMPS(s,chi)
  # psi could come from a DMRG
  # calculation instead

  # Contract MPS tensors
  # together to make full
  # wavefunction tensor:
  psiT = ITensor(1.0)
  for j=1:N
    psiT = psiT*psi[j]
  end
  #         s1  s2  s3  ...
  #         |   |   |    |
  # psiT = ============== ...
  #
  # Note:
  # A shorter way to do the above
  # is just `psiT = prod(psi)`
  #

  # Make a random tensor with
  # s and s' indices:
  A = randn(2^N,2^N)
  AT = ITensor(A,prime.(s)...,s...)
  #
  #      s1' s2' s3' ...
  #      |   |   |    |
  # AT = ============== ...
  #      |   |   |    |
  #      s1  s2  s3  ...
  #

  ApsiT = noprime(AT*psiT)

  # Note: afterward can turn ApsiT
  # back into an MPS by calling
  # `Apsi = MPS(ApsiT,s...;cutoff=1E-10)`

  return
end

Comments

Thank you; so then if I want ApsiT to be in the form of a 2^N x 2^N Julia array, would I use array(ApsiT)? It seems as if the dimensions of this new array is (2, 2, 2, 2, 2, 2, 2, 2, 2, 2); would reshaping it give me the desired format?

---

Good question: yes array(ApsiT) will convert it back to a Julia array but as you said it is really a Julia tensor with N "indices" or dimensions. If you want it to be a one-dimensional array of length 2^N you can call reshape(array(ApsiT),2^N).

Finally, note that ApsiT would not be 2^N by 2^N in the way the code is set up above. It's only a 2^N "vector" (or 2^N x 1 "matrix) because I did not include the other copy of psi that makes up the density matrix. You could include it, but as you know it doesn't really "participate" in the calculation since it's only tensor producted on, so for computational efficiency it's best to only include it at a later step or in a formal way.

---

Thank you; so if I represent this density matrix in this manner, is there an easy way to perform multiplication of rho with some other matrix, call it B, and then taking the trace? Or could this just be by first taking the matrix form of rho?

---

So if I understand correctly, you would like Tr[|psi><psi| B] where B is an operator acting on the entire Hilbert space? If so, then the best way to get this would be to rewrite it as <psi|B|psi>. Then the code above is doing |eta> = B|psi> and so at the end you would compute <psi|eta> to get the result. That set of steps would be much more efficient than forming |psi><psi| as a dense matrix and working that way.

If those steps sound good to you, then to get the inner product <psi|eta> (where eta = ApsiT in the code above) you would just contract psiT*ApsiT and call the function `scalar` on the result to get a regular scalar instead of a zero-index ITensor.

---

Yes, that is what I was thinking of, and I think I see now. So if I have a specific 2^N by 2^N matrix in mind, I would use this instead of the randn(2^N,2^N), construct the new ApsiT, perform B = prod(psiT*ApsiT), then scalar(B) to get the regular scalar value I desire
?

---

Yes, you're correct and would use your own 2^N x 2^N matrix in place of "A" in the code above.

For obtaining the result, you would not use `Prod(psiT*ApsiT)` but actually just:

    answer = scalar(psiT*ApsiT)

(The reason is that the function `prod` takes the product over an array of tensors, whereas here psiT*ApsiT is just explicitly the contraction of two ITensors with each other, and there's no array of tensors to loop over or anything.)

---

Thank you; I believe everything works now.

---

Great! Feel free to post a follow up or new question if you run into more issues.