Calculate reduced density matrix from QN MPS

ITensor Julia Questions
1

For an MPO representing the density matrix of an QN-preserving MPS, I would like to trace out some of the indices to be left with the reduced density matrix.
It was not obvious how to do this in the MPO’s native form, however I can get around this by first transforming the QN MPS to a dense MPS then proceeding. Below shows a simple example of this for a 3 qubit MPS where the first qubit is traced out.

However, doing this I obviously lose the computational benefit of working with QN MPS (commenting out the line ‘psi = dense(psi)’ gives an error about the indices requiring opposite directions to contract).
Is it possible to do this partial trace without first transforming to a dense MPS?

n = 3

sites = siteinds("Qubit", n; conserve_number=true)

state = [i<=2 ? "Dn" : "Up" for i=1:n]

psi = randomMPS(sites, state, linkdims=2)

psi = dense(psi)

sites_dense = siteinds(psi)

rho = outer(psi', psi)

rho[1] = rho[1]*delta(sites_dense[1],sites_dense[1]')
rho[2] *= rho[1]

rho_out = MPO(2)

for i in 1:2
    rho_out[i] = rho[i+1] 
end

T = ITensor(1.)
for j=1:length(rho_out)
  T *= rho_out[j]
end

A = Array(T, sites_dense[2]', sites_dense[3]', sites_dense[2], sites_dense[3])

A = reshape(A, 4, 4)

display(A)

;
2

Hi Joex, have a look at docs for dag function which can be used to change directions on indices. For these kinds of problems it is good practice the inspect the indices at each step to see what is going on. I have added to @show inds() lines to your example:

using ITensors
let
n = 3
sites = siteinds("Qubit", n; conserve_number=true)
@show sites
state = [i<=2 ? "Dn" : "Up" for i=1:n]
psi = randomMPS(sites, state, linkdims=2)
@show [inds(psi[i],tags="Site") for i in 1:n]
rho = outer(psi', psi)
@show [inds(rho[i],tags="Site") for i in 1:n]
delta1=delta(sites[1],dag(sites[1]'))
@show inds(delta1)
T = rho[1]*delta1 * rho[2] * rho[3]
A = Array(T, sites[2]', sites[3]', sites[2], sites[3])
A = reshape(A, 4, 4)
display(A)
end

You can see one use of dag() in the definition of delta1 is sufficient to resolve the error.
I hope this helps.
Regards
Jan

1 Like
3

This solved the problem, thanks for the help!

2 Likes