Calculating observables from density matrices

Question

Hey everyone, I am rather new to MPS so please forgive the naive question. I want to replicate the results from the paper, https://arxiv.org/pdf/1702.04210.pdf, to learn MPS methods. The primary quantity of interests is the spin density,
$$ s(x, t) = \text{tr} (\rho (t) S^{z}_{x})$$, where x is the spin-site. In this work, they prepare the system in a "soft" domain wall density matrix.

My approach: The way I went about this is by thinking of the density matrix as an MPO and then evolve this in time using trotterisation methods. The first problem, I have is computing the $$\text{tr} (\rho (t) S^{z}_{x})$$. You can't just use the "*" method to contract the MPO's and then trace.

My first question is how to perform this operation correctly?

My solution attempt: I just went ahead and calculated real(inner(rho(t), S^z)) for each site which matches small scale exact diagonalization methods, but I am not sure if this is still correct. Here I just used a density matrix made from a Neel state

My second question is how to go about making this initial density matrix as an MPO?

My solution attempt: I first wrote a function in Julia that can create the matrix for this soft domain wall state. I then define the necessary site indices and convert to an MPO but it doesn't have the correct structure and also this won't work for large N because it explicitly makes the matrix.

Any feedback or ides is much appreciated

Comments

Hi, so I think your question has two parts: 1. computing local observables from a density matrix MPO and 2. preparing an MPO to be some specific initial density matrix, correct?

Regarding 1, your method of using inner(rho(t),S^z) is a correct procedure. The inner function for MPO's computes the quantity inner(A,B) = Tr[A^\dagger B] which is known as the Frobenius inner product. (You can also view the MPO's as having their physical indices all folded into the same direction to "vectorize" them and then the inner product is computed by treating them as MPS with two indices per site.)

If you are just measuring S^z_j on individual sites j, there is a more efficient way however than representing S^z_j as an MPO. Our documentation could benefit from a code example about this, but what you would do is trace all of the other sites of the MPO and contract the MPO tensor at site j with your operator S^z_j represented as a single ITensor with two indices.

Finally, regarding how to prepare an initial density matrix, the answer depends on which one you are trying to prepare. Could you provide more detail about that?

---

Note that I've started an issue related to this question here: https://github.com/ITensor/ITensors.jl/issues/730

I've made a simple example code in that issue for computing the expectation value of an MPO.

---

Hey Matt,

Thank you for posting that little code. I think in your section where you say 'used like:' you have an error. I think it should be, expect(M, "Sx", 3) instead of "op(M, "Sx", 3)". Perhaps I am not understanding correctly.

---

Yes that's right, thanks for pointing that out, I've fixed it.

---

Thanks, you guys for the insightful comments let me address Miles question. The initial density matrix I want to prepare is the following,
$$\rho \sim (1+\mu S_{z})^{\otimes N/2} \otimes (1-\mu S_{z})^{\otimes N/2}$$, where $$\mu \in [-1,1]$$. You can see that for $$\mu \sim \pi/1800$$ this is effectively the infinite temperature state on each site.

Also, how exactly can I control the maximum bond dimension with this kind of evolution? The reason I ask is that if I do apply(gates, rho; cutoff=cutoff, maxdim=maxdim, apply_dag=true) for evolution I get the same curve if I change maxdim. I would expect that for larger max bond dim I get a better approximation to the dynamics.

Chris

---

`cutoff` and `maxdim` are the correct ways to control the truncation/maximum bond dimension. Have you tried adjusting both the `cutoff` and `maxdim`, or just `maxdim`? We would need more details about your calculation to understand what trouble you might be having with the truncation/convergence with bond dimension (for example, I'm not sure what curve you are trying to compute).

Note that in general with real time evolution of tensor networks, the entanglement grows linearly with time, so generally you need exponentially large bond dimensions to reach longer and longer times, though the system you are looking at may be different.

---

Sorry I am rather new to tensor networks so I may not understand everything correctly yet. But what I want to calculate is the following:
$$C(x,t)=\langle \rho (t) S^{z}_{x} \rangle$$, where x indicated the site of the chain and $$\rho(t)$$ is the evolved MPO discussed above. Specifically, I want to fix the site position and plot this curve for different bond dimensions to see the convergence.

---

I see. So when you plot @@C(x, t)@@ for different bond dimensions, you are seeing that the curve does not change with bond dimension? Is that true for every time @@t@@? At later times, the state likely becomes more entangled, so you may not be running to long enough time scales to see the entanglement saturate the bond dimension you set. Also, it is helpful to look at the entanglement of @@\rho(t)@@ as well as the smallest singular value of @@\rho(t)@@ as a function of time to get an idea for how complex your state is to represent.

---

Hey Matt,

I based my TEBD code off the gate evolution example for an MPO, but when I try to use conserve_qns=true, I get a weird svd error which says, "In 'svd' the left or right indices are empty." Should I open a new issue for this?

thanks!

---

Hi, I’d say this is a good candidate for opening an issue. We may find that it’s not a bug, but rather something about the setup of your code or inputs, etc. but it would be good in this case for us to have a look. Thanks.

Please do see if you can find the simplest case or version of your code that reproduces this issue.

Miles

Answer 1

See the discussion above.

Answer 2

For the part of your question (in comments above) about preparing the initial state MPO, since it is a tensor product of local operators, you can make each operator separately as an ITensor and then "load" them into an MPO by assigning the tensors of the MPO one by one (in a for loop of course) to be the operators you want.

For example, just showing the first site:

N = 100
s = siteinds("S=1/2",N)

rho = MPO(N)

rho[1] = op("Id",s[1]) + mu*op("Sz",s[1])

and similar for sites 2,...N/2 in a for loop (use div(N,2) to get N/2 in Julia). Then another loop to go from N/2+1 to N with the other sign for the Sz term.

Some of our algorithms may complain if there are not virtual or link indices between the tensors of the MPO. If that happens, let us know and we can recommend a course of action. Thanks!

Miles

Comments

Great thank you for an alternative method. Currently I just create the local MPO using the "op" method and create a list of the ops 1...N/2 and then N/2+1...N. Ill see which of the two ways is more efficient.

A final question if I define a function, say, "op(::OpName"Sx", ::SiteType"S=1/2"). Whenever I do this everything after the "SiteType" turns red and always tries to change the indentation within the function body. Do you know why this occurs?

---

I see. So are you putting a “ character in front of op or was that just a typo in your comment above? Otherwise, it might just be an issue with the syntax highlighting in your text editor. Did you specifically install or enable Julia syntax highlighting for .jl files?

---

No, that was just to distinguish it from the sentence. I only get this in julia when I have ::SiteType"S=1/2" in the function. No other function gives my this behavior.

---

Yes, so I think it is probably your text editor settings. Actually, it must be of course because it is an issue about colorization and indentation of text, correct?

---

It seems whenever I define a function that has an equal sign within quotation, it turns everything after that red. I am guessing it doesn't hurt the code, just thought it was weird.

For example,

function physics(N, ::thing"S=bleh")

once you have ".. = .." it turns red, just weird.

---

Yes, it’s totally an independent issue from your code, and a matter of how your text editor renders code. Most likely it is trying to do code highlighting but for a different language than Julia. So the soliton would be to enable Julia-specific highlighting.

---

Yeah I see. Thank you so much for your help, I just have one final question.

Is there a to specify the bond dimension of an MPO. For example, take the initial state we discussed above. Can I specify the bond dimension to be a particular value? Generically, because the state is a product of local operators it will have a bond dimension 1 before evolving it. Also, If I want to check how the bond dimension grows is doing maxlinkdim(rho) at each time-step sufficient?

---

Could you please say a bit more about your goal of why you want a larger bond dimension? The time evolution code will increase it for you. For the initial state, a dimension of 1 is sufficient as you observed. So is it for a different purpose?

---

Yeah, I am trying to figure out how to optimize this code because I was profiling my script and for N = 4, dt = 0.5, and total = 5, cutoff=1e-10 the memory estimate is 118.36MiB and allocates around 373185. Which is a massive amount for such a small system and time scale. As a result, if I put my code on a cluster for N > 50 it runs out of memory quickly. I am having issues finding where all these allocations occur so I was playing with changing bond dimensions to see how these values change.

---

One aspect of Julia that can be confusing is that the language and Base library performs a huge amount of allocations, but most of these are small, fast, and highly optimized (I think the language runtime probably uses some custom allocators). So that means those 373185 type numbers which are very common are not too helpful or meaningful in a practical sense. I wouldn’t worry about a number like that.

On the other hand, if you are running out of memory for N>50 then of course that’s a real problem. Please check if all of the parts of your code either (1) are controlled by a cutoff and/or maxdim such as when you do calls to things like `contract` or `svd` etc. in ITensor and (2) that you don’t have a bug where certain tensors end up with more indices than expected. For the second one it requires just printing out a lot of tensors to verify correctness.

---

Hm I see this will come from more experience. I want to TEBD evolution, so I am just creating the two sites gates and applying across the chain and then doing the reverse. I then apply the dagger to get Heisenberg evolution of an MPO. But after each time-step is it necessary to orthogonalize the MPO?

---

Hi, so in TEBD is it very important to keep the MPS or MPO in an orthogonal form. The part about applying the gates will be fine and always numerically exact, but the part where you factorize and truncate the state afterward can incur uncontrolled errors if the other tensors to the left and right are not left- and right-orthogonal.

Note that when doing TEBD on an MPO, one is really treating it as an MPS with two physical indices on every MPS tensor, in terms of how the canonical forms are defined. There aren’t really so many useful canonical forms specific to MPOs besides this MPS-derived one.

---

I see, thank you for taking all this time to answer probably rather simple questions. Currently, after I evolve the MPO I take the expectation value with some spin operator on each site. What I did was to write a function that takes the (MPO, "Sz", sitesindex, site) where the MPO is the evolved state. I always first orthogonalize the MPO before computing the expectation value. It also seems the using number conserving for my initial density matrix does not work as well.

---

I see - sounds good. I had thought you were asking about keeping an orthogonal form during the time evolution itself, which is very important to do, as I think you saw already from my answer above.

---

Correct, for the evolution I just followed the example given on one of the pages for the quantum simulator on github.