build a BlockTensorMap with input sub-Tensors with product space

[ZongYongyue]

Hi Lukas,

I would like to build a series of block tensors like [1 B], [1 B; 0 1], [B, 1], where B=e*e, and e is the electron creation operator:

pspace = Vect[I]((0,0,0)=>1, (1,1,1//2)=>1, (0,2,0)=>1)
vspace = Vect[I]((1,1,1//2)=>1)
ee = add_single_util_leg(e)⊗add_single_util_leg(e)
ee = permute(ee, ((1,3,2,4),(5,7,6,8)))

i1  = isomorphism(storagetype(ee), oneunit(vspace)*oneunit(vspace)*pspace*pspace, pspace*pspace*oneunit(vspace)*oneunit(vspace))
i2  = isomorphism(storagetype(ee), vspace*vspace*pspace*pspace, pspace*pspace*vspace*vspace)

W1 = [i1, ee]
W2 = [i1 ee; 0 i2]

BlockTensorMap([i1, ee], SumSpace(oneunit(vspace)*oneunit(vspace))*pspace*pspace,  pspace*pspace*SumSpace(oneunit(vspace)*oneunit(vspace), vspace*vspace))

BlockTensorMap([i1 ee; 0 i2], SumSpace(oneunit(vspace)*oneunit(vspace), vspace*vspace)*pspace*pspace,  pspace*pspace*SumSpace(oneunit(vspace)*oneunit(vspace), vspace*vspace))

and did I build them in a right way?

[lkdvos]

Hi Zong,

I'm not sure if I'm fully following what you are trying to construct here (I cant really deduce what e is, or how the add_single_util_leg operates -- mostly by way of where all the legs go).

Let me first add some general comments:

It's not really possible to make SumSpace(ProductSpace(...)) work in the current way I'm handling things.
The entire formalism really depends on the underlying array of N1,N2 tensormaps having N1+N2 dimensions, as it would become quite unwieldy to handle mixed cases, and I am not sure if in that case I can always uniquely line up array dimensions with tensor legs.

When constructing the input tensors of N1,N2 legs, therefore BlockTensorMap is expecting an N1 + N2-dimensional array of tensors, with sizes corresponding to each of these legs.
One way to verify this is to explicitly construct the sum space V, and then the array would have size(V).

Then, am I correct in assuming this is effectively a local tensor of a representation of a kind of "MPO" on a ladder?
If so, I think what you want is something more like either of the following solutions:

The first option really separates out the sumspaces for each leg:

V1 = oneunit(vspace) * oneunit(vspace) * pspace * pspace \leftarrow  pspace * pspace * SumSpace(oneunit(vspace), vspace) * SumSpace(oneunit(vspace), vspace))
# in finite state language: filling level (1, (2, 1)) and (1, (1, 2))
W1 = reshape([i1 e; e 0], size(V1)) # (1, 1, 1, 1, 1, 1, 2, 2)

BlockTensorMap(W1, SumSpace(oneunit(vspace)*oneunit(vspace))*pspace*pspace,  pspace*pspace*SumSpace(oneunit(vspace)*oneunit(vspace), vspace*vspace))

...

The second would be to explicitly fuse the legs together:

V1 = SumSpace(fuse(oneunit(vspace)*oneunit(vspace))) * pspace * pspace \leftarrow pspace * pspace * SumSpace(fuse(oneunit(vspace)*oneunit(vspace)), fuse(vspace*vspace)))
...

I might also really be misunderstanding what you would like to construct, but hopefully this gives you a little bit of extra information and insight in what I wanted to support with this package!

[ZongYongyue]

I’m very sorry that I didn't clearly explain the MPO I wanted to construct at the beginning, as I hadn't expected it to be somewhat tricky.

The issue is as follows:

I would like to construct the infinite-temperature state of a canonical ensemble. I am following Sec VI of the paper :Matrix product purifications for canonical ensembles and quantum number distributions, which is referred and implemented in TeNPy. In this paper, one needs to apply an MPO B_tot , built from the tensor product of the creation operators of the physical Hilbert space and the auxiliary Hilbert space, to the vacuum state |vac〉 ⊗ |vac〉N times for N particles.

The vacuum state can be easy to construct by FiniteMPS(f, elt, [pspase*pspase for _ in 1:L], vspases) with set the right and left charge to zero, but the B_tot reads as:

B_tot = [1 b1⊗b1][1 b2⊗b2; 0 1][1 b3⊗b3; 0 1][b4⊗b4; 1]

which is not a MPOHamiltonian and b1⊗b1 is not a MPOTensor, so I can not use MPOHamiltonian directly.

Therefore, I try to write a new type for B_tot and write a method for B_tot* |vac〉 ⊗ |vac〉. So, I construct ee = e⊗e first, where rank-3 tensor e is the electron creation operator and add_single_util_leg is just add a trivial leg for e to be a rank-4 tensor. Hence ,the left virtual space of ee is oneunit(vspace)*oneunit(vspace), the right virtual space of ee is vspace*vspace, the up- or down- physical-auxiliary-mixed space is pspace*pspace.

And then, I try to write every matrix of B_tot as the role of Ws in MPOHamiltonian, where Ws[1]=[1 b1⊗b1] , Ws[2]=[1 b2⊗b2; 0 1] ..., which need to build BlockTensors and then return to the beginning of the story...

Thank you for taking the time to go through my line of thought, can you give me some advice or help?

[ZongYongyue]

This is actually a follow-up to the question I raised last time about restricting tr(ρA) to a given particle-number sector. At that time, I think it should use a projection operator there. But later, I looked up some articles on this topic and found that the common approach—also the one adopted in TeNPy—is to construct the infinite-temperature state of the canonical ensemble to take place the role of the direct product state of identity operators in the grand-canonical system.

[lkdvos]

Thanks for the clarification!

I went through the paper and I think I now understand what you want to achieve. If I got everything correctly, I think the resulting block tensor would most naturally be a 6-leg object, with one virtual leg on the left and right, and two physical legs on the top and bottom.
Keeping a bit with the MPOTensor convention, literally translating all of that would boil down to:

vspace = # auxiliary space of e \otimes e, fused (so presumably something with U1Irrep(2) for particle number
Vleft = oneunit(vspace)
V1 = SumSpace(oneunit(vspace), vspace)
...
VN = SumSpace(oneunit(vspace), fuse(vspace^N))
...
Vright = fuse(vspace^L)

W1 = BlockTensorMap(
    reshape([i1, ee], 1, 1, 1, 1, 1, 2),
    Vleft * pspace * pspace \leftarrow pspace * pspace * V1
)
W2 = BlockTensorMap(
    reshape([i1, ee, 0, i2], 2, 1, 1, 1, 1, 2)
    V1 * pspace * pspace \leftarrow pspace * pspace * V2
)
...

Importantly, here i1, i2, ... are all different: the first one passes along a trivial leg, the second one a singly charged leg, the third one a triply charged one, ...

---

However, after thinking this through a bit more, I got to thinking if it would be possible to not have to deal with the auxiliary spaces at all. I'm not sure if the following makes any sense, but somehow the total result has to be neutrally charged, i.e. the projector, or the density matrix itself transforms trivially under the symmetry, since exp(H)_Q should just be the restriction of an invariant operator to a given charged subspace, which still yields an invariant operator.
Therefore, I think it might actually be the case that you want to start from vac \otimes dual(val), and act with creation operators on vac and annihilations on dual(vac).
In other words, one of the spaces is {0, 1, 2, ...}, while the other one is flip({0, 1, 2, ...}) = {0, -1, ...}.
This should be equivalent, but then the combined auxiliary space of e \otimes e is trivial, ie the combination is charge neutral, and the auxiliary space Vright = oneunit(vspace) again.

Again, this is my intuition, I have not worked this through, but that seems to be sensible in my head, and at least the charge balance makes sense then.
Additionally, this coincides with this being proportional to N_projector = creation^N * vac_projector * annihilation^N (note how here we also get a mix of creation and annihilation operators)
I also think this then coincides with the Dicke state construction, just slightly more efficiently immediately on the operator level instead of the state level.
(Disclaimer that this really might all be wrong though!)

[ZongYongyue]

Thank you for your great help!
I follow your summary, build the B^+_tot and define * for B^+_tot and FiniteMPS as following:

function physical_auxiliary_plus(physical_plus::AbstractTensorMap, auxiliary_plus::AbstractTensorMap=physical_plus)
    if (length(domain(physical_plus))==1)&&(length(codomain(physical_plus))==1)
        physical_plus = add_util_leg(physical_plus)
        auxiliary_plus = add_util_leg(auxiliary_plus)
    elseif (length(domain(physical_plus))==2)&&(length(codomain(physical_plus))==1)
        physical_plus = add_single_util_leg(physical_plus)
        auxiliary_plus = add_single_util_leg(auxiliary_plus)
    else
        throw(ArgumentError("An invalid creation operator!"))
    end
    V1, V2 = codomain(physical_plus, 1), domain(physical_plus, 2)
    A1, A2 = codomain(auxiliary_plus, 1), domain(auxiliary_plus, 2)
    iso1 = isomorphism(storagetype(physical_plus), fuse(V1, A1), V1*A1)
    iso2 = isomorphism(storagetype(auxiliary_plus), fuse(V2, A2), V2*A2)
    @plansor pap[-1 -2 -3; -4 -5 -6] := iso1[-1; 1 2] * physical_plus[1 3; -4 5] * τ[2 -2; 3 4] * 
                                τ[5 6; -5 7] * auxiliary_plus[4 -3; 6 8] * conj(iso2[-6; 7 8])
    return pap
end


function physical_auxiliary_identities(pap::AbstractTensorMap{T, S, 3, 3}) where {T, S}
    id1 = isomorphism(storagetype(pap), codomain(pap)[1]*codomain(pap)[2]*codomain(pap)[3], domain(pap)[1]*domain(pap)[2]*codomain(pap)[1])
    id2 = isomorphism(storagetype(pap), domain(pap)[3]*codomain(pap)[2]*codomain(pap)[3], domain(pap)[1]*domain(pap)[2]*domain(pap)[3])
    return id1, id2
end

ep = e_plus(elt, SU2Irrep, U1Irrep; filling=filling)
pap = physical_auxiliary_plus(ep)
id1, id2 = physical_auxiliary_identities(pap)
vspace = domain(pap, 3)
pspace = domain(ep, 1)
z0 = TensorMap(zeros, vspace*pspace*pspace, pspace*pspace*oneunit(vspace))
V1 = oneunit(vspace)
V2 = SumSpace(oneunit(vspace), vspace)
V3 = vspace

W1 = BlockTensorMap(reshape([id1, pap], 1, 1, 1, 1, 1, 2), V1*pspace*pspace ← pspace*pspace*V2)
W2 = BlockTensorMap(reshape([id1 pap; z0 id2], 2, 1, 1, 1, 1, 2), V2*pspace*pspace ← pspace*pspace*V2)
W3 = BlockTensorMap(reshape([pap id2], 2, 1, 1, 1, 1, 1), V2*pspace*pspace ← pspace*pspace*V3)


const DoubleMPOTensor{S} = AbstractTensorMap{T, S, 3, 3} where {T}
import MPSKit: right_virtualspace, left_virtualspace
right_virtualspace(O::DoubleMPOTensor) = space(O, 6)'
left_virtualspace(O::DoubleMPOTensor) = space(O, 1)

struct DoubleMPOHamiltonian{TO <: DoubleMPOTensor, V <: AbstractVector{TO}} 
    W::V
end

Base.length(dmpo::DoubleMPOHamiltonian) = length(dmpo.W)
Base.getindex(dmpo::DoubleMPOHamiltonian,  i::Int) = getindex(dmpo.W, i)

function DoubleMPOHamiltonian(Ws::AbstractVector{O}) where {O <: DoubleMPOTensor}
    for i in eachindex(Ws)[1:(end - 1)]
        right_virtualspace(Ws[i]) == left_virtualspace(Ws[i + 1]) ||
            throw(ArgumentError("The virtual spaces of the MPO tensors at site $i do not match."))
    end
    return DoubleMPOHamiltonian{O, AbstractVector{O}}(Ws)
end
import TensorKit: scalartype
scalartype(::Type{DoubleMPOHamiltonian{O, AbstractVector{O}}}) where {O} = scalartype(O)
dmpo = DoubleMPOHamiltonian([W1,W2,W2,W3])

using MPSKit: check_length

function Base.:*(H::DoubleMPOHamiltonian, mps::FiniteMPS)
    N = check_length(H, mps)
    @assert N > 2 "MPS should have at least three sites, to be implemented otherwise"
    A = convert.(BlockTensorMap, [mps.AC[1]; mps.AR[2:end]])
    A′ = similar(
        A,
        tensormaptype(
            spacetype(mps), numout(eltype(mps)), numin(eltype(mps)),
            promote_type(scalartype(H), scalartype(mps))
        )
    )
    # left to middle
    U = ones(scalartype(H), left_virtualspace(H[1]))
    @plansor a[-1 -2 -3; -4 -5] := A[1][-1 2 3; -4] * H[1][1 -2 -3; 2 3 -5] * conj(U[1])
    Q, R = leftorth!(a; alg = QR())
    A′[1] = convert(TensorMap, Q)

    for i in 2:(N ÷ 2)
        @plansor a[-1 -2 -3; -4 -5] := R[-1; 1 2] * A[i][1 3 4; -4] * H[i][2 -2 -3; 3 4 -5]
        Q, R = leftorth!(a; alg = QR())
        A′[i] = convert(TensorMap, Q)
    end

    # right to middle
    F = isomorphism(storagetype(H[N]), right_virtualspace(A[end])*right_virtualspace(H[N]), fuse(right_virtualspace(A[end])*right_virtualspace(H[N])))
    @plansor a[-1 -2; -5 -4 -3] := A[end][-1 1 2; 3] * H[N][-2 -3 -4; 1 2 4] * F[3 4; -5]
    L, Q = rightorth!(a; alg = LQ())
    A′[end] = transpose(convert(TensorMap, Q), ((1, 4, 3), (2,)))

    for i in (N - 1):-1:(N ÷ 2 + 2)
        @plansor a[-1 -2; -5 -4 -3] := A[i][-1 1 2; 3] * H[i][-2 -3 -4; 1 2 4] * L[3 4; -5]
        L, Q = rightorth!(a; alg = LQ())
        A′[i] = transpose(convert(TensorMap, Q), ((1, 4, 3), (2,)))
    end

    # connect pieces
    @plansor a[-1 -2 -3; -4] := R[-1; 1 2] * A[N ÷ 2 + 1][1 3 4; 5] * H[N ÷ 2 + 1][2 -2 -3; 3 4 6] *
        L[5 6; -4]
    A′[N ÷ 2 + 1] = convert(TensorMap, a)

    return FiniteMPS(A′)
end

I don't know this is right or not, but at least this code can run with no error throwing, so I would like to share this here.

However, I still have some questions about it:

Firstly, the $b^\dagger\otimes b^\dagger$ is a 6-leg object or a 8-leg object ? In that paper, the author said that the $|\rho_Q\rangle \in H_Q \otimes H_Q$. It means the physical and auxiliary space conserved Q respectively. But if we fuse the two virtual leg of the two $b^\dagger$ (from 8-leg to 6-leg), it seems like the total physical and auxiliary space have 2Q conserved, which means the charge can transfer from physical to auxiliary space that seems to be a grand-canonical ensemble not a canonical ensemble.

Secondly, if I start from the vac $\otimes$ dual(vac) and act creation operators on vac and annihilations on dual(vac), should I connect the two virtual leg of the creation and annihilation operators or fuse them? From the point of neutrally charged, it seems that one should connect them that return a trivial right leg (0,0) rather than fuse them got a ((0,0), (0,1)) right leg. But whatever, it also means the charge can transfer from a space to another space, and I suspect that, in this way, will we return back to the identity operator as them as the initial state of the grand-canonical ensemble eventually.

Thirdly, should I truncate the $|\rho_Q\rangle$ beacuse the bond dimension seems to grow exponentially with N ? But the author said the growth of the bond dimension is linear in N, I'm not sure about it.

Of course, I know these questions are beyond the topic of this issue, but I would be very happy and grateful if you were willing to continue to discuss them with me.

[ZongYongyue]

Oh I remember that you said it is impossible to make SumSpace(ProductSpace(...)), so it seems like one must fuse the two virtual legs of the two creation operators, but how to deal with the conserved Q in physical and auxiliary spaces respectively, or it is actually not conserved Q respectively ?