Hi! Hope this finds you well. At the high level, I am trying to do the following:
Starting with an MPO of length L, choose two sites, A and B. These sites can be located anywhere in the MPO; as such, they can be separated by an arbitrary number of sites
Construct three reduced density matrices: \rho_A, \rho_B, and \rho_{AB}
Using these reduced density matrices, compute the mutual information I(A:B)
Unfortunately, my mutual information calculation is giving me all sorts of nonsensical output. Sometimes S_A will be negative, for example, or the mutual information of two non-entangled sites will be nonzero. I feel relatively confident in the computation of the mutual information itself (I am just using example von Neumann entropy code from the ITensors documentation). This leads me to believe that my error might be in calculating my reduced density matrices. Below, I have attached my code for computing \rho_A, \rho_B, and \rho_{AB} given an input density matrix \rho and two integer index locations, A and B.
function twosite_density_matrix(ρ::MPO, A::Int, B::Int)
ρ = copy(ρ)
sites = siteinds(ρ)
@assert A != B
# the code is meant to work for A < B
if B < A
A, B = B, A
end
# we start with the tensors at A and B
ρ_A = ρ[A]
ρ_B = ρ[B]
# trace out everything on the left and multiply into ρ_A
Lside = ITensor(1.0)
for i in 1:(A-1)
Lside = Lside * ρ[i] * delta(sites[i][1], sites[i][2])
end
ρ_A *= Lside
# trace out everything between A and B and multiply into ρ_A
mid = ITensor(1.0)
for i in (A+1):(B-1)
mid = mid * ρ[i] * delta(sites[i][1], sites[i][2])
end
ρ_A *= mid
# trace out everything to the right of B and multiply into ρ_B
Rside = ITensor(1.0)
for i in (B+1):length(ρ)
Rside = Rside * ρ[i] * delta(sites[i][1], sites[i][2])
end
ρ_B *= Rside
# the composite density matrix
ρ_AB = ρ_A * ρ_B
end
Am I making any obvious mistakes? I am not able to figure out where I am going wrong.
Thank you very much for all the support! I really appreciate it
These things can be tricky but your basic approach and code patterns look good.
I think the key thing to do when writing code like this (which I always do) is to print out each ITensor (actually just its indices, by doing @show inds(T) or similar) at each step of the for loop, then drawing diagrams of each of these intermediate tensors to confirm that they have the indices you expect.
Just to say that the approach posted in Is this a correct way to calculate one-site and two-site entanglement spectrum and entropy? for calculating the n-site reduced density matrix works correctly. I used it to calculate the mutual information for an example study involving the transverse field Ising model and it successfully benchmarked with a paper that I was following for the results.
I’d found this thread from a search related to my query about correctly extracting the reduced density matrices. Hence thought of posting about what worked for me.
@rominafsh - pls check the link I posted (above your post) two months ago post that directs to a useful, easy approach to construct one- or two- sites (generalizable also to few-sites) reduced density matrix (RDM) proposed by a different user here… Once one has that, calculating most other things from it (mutual info, entanglement entropy, etc) is straightforward.
Only caution - that approach creates RDM as an itensor object and not an MPO directly. If one must work with an MPO, one can directly convert it to an MPS first and then an MPO by calling said functions in-built in iTensors.
Thanks for your reply. I checked that post and tried implementing a function based on that. But I am confused about effectively generalizing the indices. Testing the code on a random MPS, sometimes I get NaN + NaN*im results. Otherwise the resuls are mostly zero. Here’s my code:
using ITensors
using LinearAlgebra
function rdm(psi::MPS, a::Int, b::Int)
N = length(psi)
psi1=psi @assert 1 ≤ a < b ≤ length(psi) #reduced density matrix of the site a
orthogonalize!(psi1,a)
rhoa=prime(psi1[a],sites[a])*dag(psi1[a])
#reduced density matrix of the site b
orthogonalize!(psi1,b)
rhob=prime(psi1[b],sites[b])*dag(psi1[b])
#two-site reduced density matrix from a site to b site
psi1dag=prime(dag(psi1),linkinds(psi1))
rhoab = prime(psi1[a],linkinds(psi1,a-1))*prime(psi1dag[a],sites[a])
for i=a+1:N
rhoab *= psi1[i]
rhoab *= psi1dag[i]
end
rhoab *= prime(psi1[b],linkinds(psi1,b+1))*prime(psi1dag[b],sites[b])
return rhoa, rhob, rhoab
end
function entropy(ρ)
D, _ = eigen(ρ)
S = 0.0
for n=1:dim(D, 1)
p = D[n,n]^2
S -= p * log2(p)
end
return S
end
In your example above, provide some starting linkdims to the randomMPS, because by default randomMPS creates MPS with linkdims=1 and so there’s no entanglement to begin with… always a good idea to give the initial randomMPS some number of linkdims with for example randomMPS(sites, linkdims=5)