Purifying MPO into MPS

ITensor Julia Questions
1

I’m trying to take an MPO representing a density operator and turn it into a doubled MPS via contracting Bell states to the unprimed indexes.

However, when I try to split the newly contracted tensors into two using SVD to get “skinny” indexes some are the right dimension (aka dim 2 for qubits) but others are dim 4.

Here is the code snippet I’m using to do this: I really appreciate the help!

function purify_MPO_to_MPS(M::MPO)::MPS
    N = length(M)
    
    s_unprimed_vec = [siteind(M, j, plev=0) for j in 1:N]
    s_primed_vec = [siteind(M, j, plev=1) for j in 1:N]

   
    doubled_sites = [addtags(s_unprimed_vec[j], "doubled") for j in 1:N]

    doubled_mps_tensors = ITensor[]
    
    for j in 1:N
        mpo_tensor = M[j]
        s_unprimed = s_unprimed_vec[j]
        s_primed = s_primed_vec[j]
        s_doubled_new = doubled_sites[j]

       
        bell_state = ITensor(s_doubled_new, s_primed)
        for k in 1:dim(s_doubled_new)
            bell_state[s_doubled_new(k), s_primed(k)] = 1.0
        end

        
        combined_tensor = mpo_tensor * bell_state
        
        
         U, S, V = svd(combined_tensor, s_unprimed, s_doubled_new)
         lefttensor=U
         righttensor=S*V
        
        push!(doubled_mps_tensors, lefttensor)
        push!(doubled_mps_tensors, righttensor)
    end

    
    purified_psi = MPS(doubled_mps_tensors)
    
    return normalize(purified_psi)
end
2

Nothing looks immediately wrong with your code, though I haven’t tested it myself. Some things to check:

  • what are the combined dimensions (products of dimensions) of the indices you are passing in? As you may know, if you SVD a tensor which is m \times n when reinterpreted as a matrix (which is how SVD of a tensor works underneath), the maximum possible rank will be \min(m,n). Does the rank you get match this bound?
  • if you are expecting the actual or approximate rank to be below the theoretical maximum, then you may have meant to do a truncated SVD, where singular values equal to zero or approximately zero get discarded, along with the corresponding slices of U and V. To do this, pass either the maxdim keyword or the cutoff keyword or both to the svd function.
3

Actually, I think it might be working correctly after all. However, I do want to run something by you @miles. The point of this code is to run a quantum circuit, trace out some ancillas and then compare the “purified” state to a state evolved via a Lindbladian operator in the Doubled space. I haven’t written the Doubled Linbladian evolution yet, but I was wondering if in context, this code for preparing the initial state to be a Doubled MPS would actually work for this.

using ITensors
using ITensorMPS
using Zygote




function purify_MPO_to_MPS(M::MPO)::MPS
    N = length(M)
    print(siteinds(M))
    s_unprimed_vec = [siteind(M, j, plev=0) for j in 1:N]
    s_primed_vec = [siteind(M, j, plev=1) for j in 1:N]

   
    doubled_sites = [addtags(s_unprimed_vec[j], "doubled") for j in 1:N]

    doubled_mps_tensors = ITensor[]
    
    for j in 1:N
        mpo_tensor = M[j]
        s_unprimed = s_unprimed_vec[j]
        s_primed = s_primed_vec[j]
        s_doubled_new = doubled_sites[j]

       
        bell_state = ITensor(s_doubled_new, s_primed)
        for k in 1:dim(s_doubled_new)
            bell_state[s_doubled_new(k), s_primed(k)] = 1.0
        end

        
        combined_tensor = mpo_tensor * bell_state
        
        
         U, S, V = svd(combined_tensor, s_unprimed, s_doubled_new)
         lefttensor=U
         righttensor=S*V
        
        push!(doubled_mps_tensors, lefttensor)
        push!(doubled_mps_tensors, righttensor)
    end

    
    purified_psi = MPS(doubled_mps_tensors)
    
    return normalize(purified_psi)
end




function partialTrace(rho::MPO,leftsite::Int,rightsite::Int)

  R = ITensor(1.0)
  for site in leftsite:rightsite
        R *= tr(rho[site])
  end
    rho[leftsite-1]*=R
  sites=siteinds("Qubit",length(rho)- (rightsite-leftsite+1))
  newrho=MPO(sites)
  for j in 1:leftsite-1
        newrho[j] = rho[j]
  end
  return newrho

  end


function forward(theta,phi)
    #initialize 6 qubit MPS GHZ state
    s=siteinds("Qubit",6)
    state = fill("0", 6)
    ghz = productMPS(s, state)
    gates::Vector{ITensor} = []
    #construct the circuit
    push!(gates, op("H", s[1]))
    push!(gates, op("CNOT", s[1], s[2]))
    push!(gates, op("CNOT", s[2], s[3]))
    for j=1:3
         push!(gates, op("Rz", s[j]; θ=theta))  
         push!(gates, op("Ry", s[j+3]; θ=phi))   
        #push!(gates, op("SWAP", s[j+1],s[j+3]))
         push!(gates, op("CZ", s[j], s[j+3]))
         #push!(gates, op("SWAP", s[j+1],s[j+3]))

    end
    #apply the circuit to the GHZ state
    ghz= apply(gates, ghz)
    #compute the reduced density matrix by partial trace -> turn into MPO
    rho = outer(ghz, dag(ghz))
   
    
    rho=partialTrace(rho,4,6)
    

   
    rho=purify_MPO_to_MPS(rho)
    return rho


end



ghz=forward(0.1,0.2)
println("Reduced density matrix after partial trace:")
println(siteinds(ghz))

I appreciate your help. My main concern, is that the dimensions of this output state would not be consistent with the output state of the Lindbladian evolution.

4

Actually something weird is going on in purify_MPO_to_MPS. If I print out the inds(lefttensor) it includes the doubled index. But the final tensor does not include these indeces.