Open quantum systems

I wanted to work with MPS as density matrices and apply gates on it. So I follow a similar post (ITensor for Open quantum systems - #5 by Fra_perciavalle) and vectorize a density matrix and would then proceed with a (enlarged) MPS, see code below.

using ITensors

N = 3  # number of sites
sites = siteinds("S=1/2", N)
states = ["Up" for n=1:N]
psi = MPS(sites, states)

# https://itensor.discourse.group/t/itensor-for-open-quantum-systems/271/5
rho = outer(psi', psi)
rho_vec = convert(MPS, rho)

Questions:

  1. This code takes about 37 seconds to execute on my laptop. Is it to be expected or is there any improvement?
  2. Probably a very basic question (sorry): How do I get the indices of rho_vec? I now that some of the indices are sites, but how do I get the others?
  3. Is there a simple way to see print the elements of rho_vec? I have read about MPS and MPO Examples · ITensors.jl, but I was wondering whether there exists a simpler approach to print all elements.

Hi,

  1. This code takes about 37 seconds to execute on my laptop. Is it to be expected or is there any improvement?

Does it take as long to run it if you repeat these exact lines of code once again? The (rough) way Julia works is that it compiles the code just ahead of time, which is to say that it compiles all the functions that you’ve called the first time you call them, so that every subsequent use is very fast. On my laptop, the first run took about as much time as yours, but running the exact same code again takes 0.04 seconds. This is to be expected.

One way to speed this up would be to use a system image using ITensors.compile() and following the instructions given here

  1. Probably a very basic question (sorry): How do I get the indices of rho_vec? I now that some of the indices are sites, but how do I get the others?

I use the siteinds function to get all the site indices. If you need all the inds (link indices as well), I would just iterate through each of the ITensors rho_vec[1], rho_vec[2], etc. and use inds(). Alternatively, if you want the link index between two sites j and j+1 I would use linkind.

  1. Is there a simple way to see print the elements of rho_vec? I have read about MPS and MPO Examples · ITensors.jl, but I was wondering whether there exists a simpler approach to print all elements.

Again with a case like this, I use a loop like

N = length(rho_vec)
for i in 1:length(rho_vec)
    print(rho_vec[i])
end

This might be a naive or too “low-level” a solution, but if there are better ways to do this, I’m sure somebody will add to it. Hope this helps!

1 Like

Thanks for your quick answer.

  1. It even took about 37 seconds when I repeated the same line of codes. But ITensors.compile() massively speeds it up :slight_smile:

  2. Cool, that was what I needed.

  3. Okay, satisfies my needs.

Regarding the time of the code, were you running it like

$ julia mycode.jl

on the command line? Or were you loading it into the Julia console (or “REPL”) like this:

$ julia
julia> include("mycode.jl")

?

It’s only in the second case where it could be faster when you run it twice. In the first case (running from the command line) that Julia has to recompile everything each time. But glad the ITensors.compile() command helped a lot, and we expect to be able to compile more and more parts of ITensor in the future as the Julia ecosystem continues to develop.

Of course! If you want a slightly more succinct answer to 3, I just realised that you could also do the following:

for Itens in rho_vec
    print(Itens)
end

This sort of holds (semi-) generally, in that sites, Indexsets (like inds(psi[3])) or even ITensors themselves like psi[1] can be directly iterated over. So if you wanted to see what was inside any of these objects, you would just do something like

for a in ObjOfInterest
    print(a)
end

where ObjOfInterest would be a variable of any of the above types.

I ran it like julia mycode.jl.