+1 vote
asked by (240 points)
edited by

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.

1 Answer

+1 vote
answered by (70.1k points)
selected by
 
Best answer

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
commented by (240 points)
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?
commented by (70.1k points)
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.
commented by (240 points)
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?
commented by (70.1k points)
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.
commented by (240 points)
edited by
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
?
commented by (70.1k points)
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.)
commented by (240 points)
Thank you; I believe everything works now.
commented by (70.1k points)
Great! Feel free to post a follow up or new question if you run into more issues.
Welcome to ITensor Support Q&A, where you can ask questions and receive answers from other members of the community.

Formatting Tips:
  • To format code, indent by four spaces
  • To format inline LaTeX, surround it by @@ on both sides
  • To format LaTeX on its own line, surround it by $$ above and below
  • For LaTeX, it may be necessary to backslash-escape underscore characters to obtain proper formatting. So for example writing \sum\_i to represent a sum over i.
If you cannot register due to firewall issues (e.g. you cannot see the capcha box) please email Miles Stoudenmire to ask for an account.

To report ITensor bugs, please use the issue tracker.

Categories

...