Obtaining eigenstates of reduced density matrix

Question

Hi Miles,

In literature people often obtain eigenstates of reduced density matrix from a DMRG calculation (e.g. Cincio & Vidal, PRL 110, 067208 (2013)). While there is a nice code formula to calculate the entanglement spectrum, I couldn't find one for this seemingly different task. Is there a way of implemeting this calculation using ITensor? Thanks!

Best,
Chengshu

Comments

From linear algegra, for an hermitian matrix H, if we can find a unitary matrix U so that U^{\dagger} H U = diag{h_1, h_2, ..., h_n}, then basically the eigenstate corresponds to h_i is the i'th column of U.

Naïvely I expect the same thing happens for the reduced density matrix. Here after we "gauge" the MPS to, say, site i, then we automatically have something like \rho = U^{\dagger} D U, where D is diagonal and "right orthogonality" is used. So I guess all we need to do then is to extract the "column vector" from the U. Is that the correct way of doing this? Thanks!

Answer 1

Hi Chengshu,
This is kind of a tricky question to briefly answer on a message board and I'd really need to refer you to an article such as Schollwock's review on DMRG and MPS to give you the full answer. But the basic answer is that if you properly gauge an MPS, then the MPS tensors themselves already give you the eigenstates of the reduced density matrix of a block of sites starting from either the left or right edge of the system, up to the site that is the gauge "center". (You can see why this is hard to explain in writing - it's really something that's easier to draw with diagrams.)

On the other hand if the reduced density matrix you are asking about (you didn't specify) involves tracing out a block of sites that doesn't start from the left or right edge of the system (as defined by the MPS path) then it's highly non-trivial to get the reduced density matrix eigenstates. An example of such a hard case would be a 1d system where the r.d.m. you want is the one obtained by tracing out every even numbered site. In that case you can technically get information about the reduced density matrix by using a sampling procedure but it would be very hard to get the eigenstates of such an r.d.m.

Hope that helps you to understand the issue/question better.

Miles

Comments

Hi Miles,

Thanks for your reply! Yes I am just calculating the r.d.m from the edge so it's basically the MPS. And I think I understand why the gauged MPS is just the eigenstate from the diagram. Now I want to extract the eigenstate corresponding to some particular eigenvalues, which I think I need to multiply the MPS by something like "setElt". How should one define the setElt to do this? Thanks!

---

Hi Chengshu,
Great I think from the sound of it that you understand the part of my answer that was hard to express in a brief reply on a forum (i.e. without diagrams).

So for the part of setting an index to a particular value, then yes setElt is what you need. The index sticking out of the MPS tensor that labels the rdm eigenvalues is the thing you want to fix to a particular value (rather than it being a whole index that could take a range of values). So call setElt(I(n)) where I is the index and n is the value you want to set. All setElt does is make a tensor with the single index I and all zero elements except for a "1" in the n'th position (i.e. a standard basis vector e_n). Contracting with the tensor setElt returns will set the index sticking out to the value that corresponds to the entry with the "1" of course (the nth entry).

The only possible complication is that if you are working with IQTensors and IQIndices you have to keep in mind the arrow direction of the IQIndex, so you may need to use
dag(setElt(I(n))
(or you can use dag(I) in place of I). You'll know about the arrows if you just draw the diagrams with the arrows.

Miles

---

Hi Miles,

Thanks for your explancations. To define the setElt I think the index should somehow match that of the existing MPS. How can I declare the index?

---

Hi, so you don't want to "declare" the index because declare means making a new one. What you want to do is obtain the existing one. A good way to do this is to use the commonIndex function. The function commonIndex works by taking in two ITensors (or IQTensors) and then returning the Index (or IQIndex) they share in common (assuming they share exactly one index, as should be the case for two neighboring MPS tensors). Let me know if that doesn't give you what you need and perhaps I can suggest something else.

Miles

---

After some effort the code is now totally working. Thanks very much!

---

Glad to help and to hear that you're making progress!