How to Ensure Partial Trace of an MPO Returns an MPO?

[ktrdcp]

I’m working with MPOs in Julia using the ITensor library and need some help with a function for performing a partial trace. My goal is to trace out a specified number of sites from both the left and right sides of the MPO and obtain the result as an MPO. However, I am encountering issues where the output of the partial trace function is an ITensor rather than an MPO.

Background

I have an MPO W. I need to compute the partial trace over a certain number of sites from the left and right, and then return the resulting reduced density matrix as an MPO.

Problem

The partial trace function I’ve written does not return the result as an MPO. Instead, it returns an ITensor, which is not what I need. I want to ensure that the final result is an MPO that correctly represents the reduced density matrix.

Current Code

Here’s the current implementation of my partial trace function:

function PartialTrace(W::MPO, l::Int, r::Int)
    # Get the total number of sites
    N = length(W)
    
    # Validate input
    if l + r > N
        error("The sum l + r cannot be greater than the total number of sites N.")
    end

    # Determine the indices for the sites to keep
    keep_sites = l+1:N-r
    
    # Initialize the reduced density matrix (RDM)
    rdm = W[keep_sites[1]]
    L = R = ITensor(1.)
    
    # Contract the MPO tensors for the middle sites
    for j in 2:length(keep_sites)
        rdm *= W[keep_sites[j]]
    end
    
    # Trace out the leftmost l sites
    for site in 1:l
        L *= tr(W[site])
    end
    
    # Trace out the rightmost r sites
    for site in N-r+1:N
        R *= tr(W[site])
    end

    rdm = L * rdm * R

    return rdm
end

Question

Is there a standard or more efficient approach in ITensor for performing this operation and maintaining the original MPO structure in the middle?

I tried to make \texttt{siteinds} inherited from the original one, but I haven’t completed it yet.

I would greatly appreciate any guidance or suggestions from the community on how to resolve this issue.

Thank you in advance for your help!

[ktrdcp]

I think I’ve found a solution. Here’s renewed code:

function PartialTrace(W::MPO, l::Int, r::Int)
    # Get the total number of sites
    N = length(W)
    
    # Validate input
    if l + r > N
        error("The sum l + r cannot be greater than the total number of sites N.")
    end

    # Determine the indices for the sites to keep
    keep_sites = l+1:N-r
    
    # Initialize the reduced density matrix (RDM) as an MPO
    rdm = W[keep_sites[1]]
    L = R = ITensor(1.)
    
    # Contract the MPO tensors for the middle sites
    for j in 2:length(keep_sites)
        rdm *= W[keep_sites[j]]
    end
    
    # Trace out the leftmost n sites
    for site in 1:l
        L *= tr(W[site])
    end
    
    # Trace out the rightmost m sites
    for site in N-r+1:N
        R *= tr(W[site])
    end

    rdm = L * rdm * R

    M = MPO(rdm, sites[l+1:N-r])

    return M
end

After organizing my thoughts to make the post, I gained a deeper understanding of the problem. In any case, thank you very much!

[miles]

Yes, that’s the solution I would have suggested too.

One small twist on it is that instead of making the “rdm” part in the middle to be a single tensor, which is good when the region is small but exponentially expensive when the region is large, is to leave the middle as an MPO (with separate tensors for each site) and just absorb your L tensor into the first one, and R into the last one.