Computing Renyi entropy of a sub-system in the middle of a 1D chain

[tamoghnaRay]

Hi, I have a spin- 1/2 XXZ chain of length L. I have a sub-system of length L_S, in the middle of the chain. Say L = 10 and the system goes from site 4 to 7 (L_S = 4). I want to compute S_n = \frac{1}{1-n}\log[\text{Tr} \rho_S^n]. I am using the following code to compute this -

Function to compute partial trace -

function PartialTrace(rho::MPO, # MPO of the total density matrix
                        start::Int, # start index of the sub-system
                        stop::Int # end index of the sub-system
                        )
    L = length(rho) # total system size
    M = stop - start + 1 # sub-system size
    L1 = ITensor(1.)
    for i = 1:start-1
        L1 *= rho[i] # contracting all the link indices of the matrices on the left
        L1 = tr(L1) # taking trace over the left part
    end
    # L1 = tr(L1) # taking trace over the left part
    rho[start] *= L1 # absorbing the left contraction to the first site of the sub-system
    R1 = ITensor(1.)
    for i = stop+1:L
        R1 *= rho[i] # contracting all the link indices of the matrices on the right
        R1 = tr(R1) # taking trace over the right part
    end
    # R1 = tr(R1) # taking trace over the right part
    rho[stop] *= R1 # absorbing the right contraction to the last site of the sub-system
    rho_reduced = MPO(M) # initializing the reduced density matrix MPO
    for i=1:M
        rho_reduced[i] = rho[start-1+i]
    end 
    return rho_reduced 
end

Computing reduced density matrix-

rho = outer(psi', psi) # psi is the MPS of the full system
rho_sys = PartialTrace(rho, start, stop) # reduced density matrix of the sub-system

There are two ways that I have used to compute S_2. Method 2 seems faster, at least when L_S is small.

# method - 1 
rho_sys_sq = apply(rho_sys, rho_sys)
tr_rho_sys_sq = tr(rho_sys_sq)
S2 = -log(real(tr_rho_sys_sq))

# method - 2
rho_sys_mat = prod(rho_sys)
D,U = eigen(rho_sys_mat, ishermitian=true)  
u_idx = commonind(U, D)
eigenvalues = [D[u_idx => k, prime(u_idx) => k] for k in 1:dim(u_idx)]
S2 = -log(sum(eigenvalues.^2))

The step that takes the most time is when I am computing rho = outer(psi', psi) . This time is quite significant compared to the rest of the steps, especially when the bond dimension is high. I have also used rho = projector(psi). Both take a similar time to run.

1. Is there a more efficient way to compute the MPO rho, or is there a way to compute S_2 directly from the MPS psi?
2. Another step that takes a lot of time is rho_sys_mat = prod(rho_sys), where I compute the dense tensor from the MPO. Is there a better way to do this?

[miles]

Yes, I’m not surprised that outer(psi', psi) is a bottleneck. The short answer is that making a density matrix that way is not an efficient way to do it and there are much more efficient ways.

One thing that might immediately help you is the ITensorEntropyTools.jl package:
https://itensor.discourse.group/t/introducing-itensorentropytools-jl/2192
One of the main features of this package is that it automatically determines the best strategy for computing density matrices of the type you are computing (there are roughly two good ways to compute these and it figures out which one).

If that package doesn’t cover your case, please reply below and we could give you additional pointers about how to do this calculation.

[tamoghnaRay]

Hi @miles, thanks for the suggestion. This package works fine for my case.

[miles]

Glad to hear it!