Is there a way to project a dense MPS into a symmetry sector?

DMRG and Numerical Methods
1

Hi, ITensor Teams！

First, a little bit of context, for a study in open quantum systems, I need to vectorize the identity:

I = \sum_{n_{1}, \cdots, n_{N} = { 0,1}} |n_{1} \cdots n_{n} \rangle \langle n_{1} \cdots n_{n} | \rightarrow |I \rangle =\sum_{n_{1}, \cdots, n_{N} = { 0,1}} |n_{1} n_{1}\cdots n_{n} n_{n} \rangle

Basically, I am adding ancilla sites at the right of each physical site, hence doubling the total number of physical sites. Implementing |I \rangle as dense MPS is straightforward for spinless fermions because it is a product state:

| I \rangle = \bigotimes _{i =1}^{M} \left( |0 0 \rangle + | 1 1 \rangle \right) 
using ITensors
using ITensorMPS

function Build_left_vacuum(sites::Vector{Index{Int64}})
    #The left vacuum |I> will be used to calculate the expectation values and the norms.

    M = Int(length(sites)/2)
    left_vacuum  = MPS(sites,"0") #we can start with a simple state like this, and apply an operator to create the superposition that characterizes |I>.
    #Each site will be converted into a superposition upup + dndn in order to create the left_vacuum.
    Create_Superposition_Matrix = [1 0 0 0; 0 0 0 0; 0 0 0 0; 1 0 0 0]  #Two-site operator that converts |0,0> into (|0,0> + |1,1>). Note that it is unnormalized.
    
    for i = 1:M 
        k = site_SF(i) 
        Create_Superposition_OP = op(Create_Superposition_Matrix, [sites[k], sites[k+1]]) #it affects the site and its ancilla
        left_vacuum = apply(Create_Superposition_OP, left_vacuum)
    end

    return left_vacuum
end

M = 10
sites = siteinds("S=1/2",2*M)
I_vec = Build_left_vacuum(sites);

This build |I \rangle as a dense MPS and works for most of my calculations. However, in some calculations, my Liouvillian conserves the total number of fermions, so I can perform TEBD using symmetries; clearly, an MPS |\psi \rangle from that process would be a block-sparse MPS.

In general I need to project |\psi \rangle against \langle I | for calculating expectation values, so far what I am doing for that is taking the dense version of |\psi \rangle by using NDTensors.dense(), that works but it is not very efficient and sometimes has numerical errors from considering very small singular values. Note that I can not define | I \rangle directly in half-filling because its product state form has contributions from all particle number symmetry sectors (is based on the completeness of the Fock basis).

I want to know if there is a way to project a dense MPS into a symmetry sector (say, half-filling). In ED, this would be done by defining the projector \hat{P} with the eigenstates of \hat{N} that have N/2 particles as columns. However, I am not sure how to do it in the language of tensor networks. I am not sure if it is possible to build an MPO with dense input indices and block-sparse output indices. Any idea that can help me?