Hi! Hope you are doing well. I am running into a strange issue when I try to calculate the trace of my reduced density matrices. Here are some of the cases and associated outcomes:
Case 1: I construct the full density matrix of the system and take the trace. In this case, I get the expected result \text{Tr}[\rho] = 1. This result does not depend on whether the system has an odd or even number of sites.
Case 2: I construct a reduced density matrix (see code below) with an odd number of sites and then take the trace. Again, I get the expected behaviour with \text{Tr}[\rho] = 1.
Case 3: Finally, if I try to build a reduced density matrix with an even number of sites, I run into this strange behaviour where my density matrix has a trace of -1.
My impression: I believe that this might have something to do with the hanging link indices of dimension 1 that are left over on the left- and right-hand sides of the MPO after I calculate the reduced density matrix. For example, in the correct full density matrix formulation, I have only one link index on the left-most MPO element, and again one link index on the right-most MPO element. But in my construction of the reduced density matrix, I end up with two link indices on my edge sites, one of which is just hanging out.
Here is my code. I’ve included everything to reproduce my results save for the code to get the ground state of my system. The toy example I am working with is the 100-site TFIM with an applied h field of 0.5
function make_density_matrix(ψ::MPS, indslist::Vector{Int})
sort!(indslist)
sites = physical_indices(ψ, indslist)
ψdag = dag(prime(ψ, sites))
prime!(ψdag, "Link")
ρ = outer(ψ, ψdag)
# get rid of hanging indices
# put density matrix into MPO form
orthogonalize!(ρ,indslist[1])
ρ.llim = 0
ρ.rlim = 2
ρ.data = ρ[indslist]
return ρ
end
## PARAMS FOR TRANSVERSE-FIELD ISING ##
N = 100 # number of sites
h = 0.5 # field
# Do DMRG
ψ0 = find_groundstate(N,h)
# Select the indices I want to RETAIN in the reduced density matrix
# (NOT the indices I trace over; note the slightly weird convention)
inds_to_keep = collect(1:length(ψ0))
ρ = make_density_matrix(ψ0, inds_to_keep)
@show ρ
@show length(inds_to_keep)
@show tr(ρ)
inds_to_keep = collect(2:3)
ρA = make_density_matrix(ψ0, inds_to_keep)
@show ρA
@show length(inds_to_keep)
@show tr(ρA)
inds_to_keep = [2]
ρB = make_density_matrix(ψ0, inds_to_keep)
@show ρB
@show length(inds_to_keep)
@show tr(ρB)
Hi Nicole,
Thanks for the question. If I am correctly understanding the goal of your code, you are intending to take two copies of an MPS and sum over (trace out) a subset of the site indices, resulting in a reduced density matrix.
However, it looks like you’re trying to do this using the function outer which doesn’t do what you need. Specifically, outer does not contract any tensor indices or sum out / trace any of the site indices of the MPS’s given to it.
I think your best bet is to compute the r.d.m. by taking individual ITensors out of the MPS and then contracting them together in a certain pattern so that you end up with the r.d.m. as an ITensor. Or if you want to represent the r.d.m. as an MPO for some purpose please let me know and we can discuss how to do that.
My goal is to actually create a bipartite reduced density matrix in the form of an MPO. Suppose I have an MPS that is 100 sites long. Then I would want to create a reduced density matrix \rho_{AB} where region A might correspond to sites 5-10 and region B might be sites 90-92, for example. The main point here is that there exists a significant buffer between these two regions that needs tracing out, which I think makes it incorrect to just orthogonalize appropriately and then extract tensors out of the MPO(?)
Thanks for clarifying. One thing I don’t understand about the example you just gave though is that region A and region B don’t add back up to cover the whole system. So would a better example be that region A is sites 5-10 and B is all of the other sites 1-4 and 6-100 ?
Apologies for the lack of clarity! What I mean to say is that the original system might be composed of three parts–say, A, B, and C. Subsystems A and B are separated by C. My goal is to trace out C, thus leaving behind only A and B (in this example, I am also taking C to be some sites at the beginning and end of the chain; basically, everything else that is not in A or B). Here is some quick code I wrote up that I believe(?) may compute the partial trace, but I am not entirely sure. At least the trace of the reduced density matrix is 1, as expected:
findnearest(A,x) = argmin(abs.(A .- x))
function density_matrix(ψ::MPS)
sites = siteinds(ψ)
ψdag = dag(prime(ψ, sites))
prime!(ψdag, "Link")
return outer(ψ, ψdag)
end
function partial_trace(ρ::MPO, indslist::Vector{Int})
sort!(indslist)
ρ̃ = copy(ρ)
sites = physical_indices(ρ)[indslist]
# Go through all the traced-out sites and tie them together
for (i,s) in zip(indslist, sites)
ρ̃[i] *= delta(s, prime(s))
end
# Retrieve all the sites we are keeping
to_keep = setdiff(collect(1:length(ρ)), indslist)
sites_to_keep = physical_indices(ρ)[to_keep]
# For each site being traced out, find the nearest site that is being kept
blocks = [to_keep[findnearest(to_keep, k)] for k in indslist]
# Orthogonalize around the first site that is being kept
orthogonalize!(ρ̃, minimum(to_keep))
for i in 1:length(indslist)
# Multiply the link indices into the nearest site that we keep
nearest_site = blocks[i]
ρ̃[nearest_site] *= ρ̃[indslist[i]]
end
ρ_new = randomMPO(sites_to_keep)
orthogonalize!(ρ_new, 1)
ρ_new.data = ρ̃[to_keep]
return ρ_new
end
## PARAMS FOR TRANSVERSE-FIELD ISING ##
N = 100 # number of sites
h = 0.5 # field
# Do DMRG
ψ0 = find_groundstate(N,h)
A = vcat(5:13)
B = vcat(89:91)
ρ = density_matrix(ψ0)
ρA = partial_trace(ρ, setdiff(collect(1:length(ψ)), [A...,B...]))
Thanks again for all the help! I really appreciate it
I see - that makes more sense now. So before discussing any code, as far as I am aware, computing a reduced density matrix (rdm) of a region other than 1:L or L+1:N (for some L < N) is not something that can be done efficiently with an MPS.
I say “efficiently” because you can write a code to do it that will work in theory, such as by using outer to make an MPO from an MPS and then partial tracing it. However, many of the steps will scale as the fourth power of the bond dimension of the MPS, whereas DMRG (for example) scales as the third power. If the bond dimension is 100, say, that means your rdm computation would have a cost like 100 DMRG calculations.
So just mentioning that you might want to find a different quantity to compute unless you are planning to use MPS with rather small bond dimensions.