I want to calculate the entanglement entropy of a pure state. I find the julia version code
using ITensors
function entropy_von_neumann(psi::MPS, b::Int)
s = siteinds(psi)
orthogonalize!(psi, b)
_,S = svd(psi[b], (linkind(psi, b-1), s[b]))
SvN = 0.0
for n in 1:dim(S, 1)
p = S[n,n]^2
SvN -= p * log(p)
end
return SvN
end
N = 10
s = siteinds("S=1/2", N)
psi = randomMPS(s, 4)
SvN = entropy_von_neumann(psi, b)
But I have two naive questions about the sentence
_,S = svd(psi[b], (linkind(psi, b-1), s[b]))
Why we need to add the site index s[b] in the arguments of svd? Does svd(psi[b],linkind(psi,b-1)) not work?
Also, I think the svd returns three objects U,S,V. Why do we need just _, S here?
the reason for the two indices (linkind(psi, b-1), s[b]) is that here we want the entanglement across bond b (dividing sites 1,2,…b from sites b+1,b+2,…N). The index linkind(psi, b-1) represents the Hilbert space of all of the sites s[1], s[2], … , s[b-1] compressed together and then we also want to include s[b] in the set of sites that we are factorizing to the left (onto the “U” of the SVD). If we did not include s[b] in the svd function call, then it would go onto the “V” tensor afterward and the entanglement we would be computing would be for bond b-1 which isn’t what we want.
the reason only S is needed to compute entanglement goes back to the definition of entanglement entropy. For a bipartition of the system into subsystems A and B, entanglement is defined by first taking a state |\Psi\rangle and factorizing it as
|\Psi\rangle = \sum_n s_n |u_n\rangle |v_n\rangle
where the s_n are real, non-negative numbers such that \sum_n s^2_n=1 and the |u_n\rangle and |v_n\rangle are orthonormal bases for the A and B subsystems. This form is known as the Schmidt decomposition and it always exists. Then the entanglement is given by S_{vN} = -\sum_n s^2_n \log{s^2_n}. So if you look at the above Schmidt decomposition, it is nothing more than an SVD of the wavefunction (viewed as a huge matrix) and the von Neumann entanglement entropy S_{vN} only involves the singular values s_n of the SVD and not the U or V matrices which define the |u_n\rangle and |v_n\rangle states.
As to why this correspondence between singular values and Schmidt values (which define the entanglement) continues to hold for a state represented as an MPS when you just SVD a single MPS tensor, that has to do with the left and right orthogonality conditions imposed on the other MPS tensors which is done in the code above by calling orthogonalize!(psi,b). It’s not an obvious thing and I’m not explaining it fully here, but just pointing out that the MPS orthogonality is the reason you can just SVD a single tensor to get the Schmidt values (singular values) and compute the entanglement.
orthogonalize!(psi, b)
_,S = svd(psi[b], (linkind(psi, b-1), s[b]))
Is it correct or not? Suppose I still want to calculate the entanglement entropy between sites b and b+1.
I think the svd function returns three matrices U, S, V. If I use _,S = svd(...) then _ is my matrix U , right? Then where is matrix V ? (I guess this question is due to my ignorance about julia syntax).
Yes I believe those two pieces of code will give the same result. The only difference is if you also save the “U” tensor in both lines of code, they will be different (one of those code’s “U” will be the other code’s “V”)
I believe if you just write _,S = svd(...) or U,S = svd(...) then the V just gets dropped or ignored, similar to if you don’t catch the return value of a function
I’m very happy to be associated with you all, this is one of the best forums with the perfect solutions with an expert (@miles) team for any complex query.