Stack multiple MPOs on an MPS

Dear Julia team,

I want to speed up repeated MPO-MPS multiplications by first stacking a bunch of MPOs (maybe lazily) and only then contracting them with the MPS, instead of contracting them one at a time.

I am not an expert, but it feels like this method will better suit a GPU processor, also.

My problem is to get the indices right… I get confused with the prime levels and my code has bugs that I can’t seem to fix on my own.

The solution mustn’t assume that the quantum systems are solely spin or bosonic, since in my system each quantum system contains both.

Hereby I attach a minimal code that I wrote to show what I am trying to achieve (and how badly my code is written, also). I emphasize that this code does not work, however, It shows what I want to compute- a series of MPOs that should be contracted simultaneously:

#This code takes an MPS of spins in the all-down state and flips them all one at a time 

using ITensors

#State preparation
function surface_tensor(dimension)
    tens = zeros(1, dimension)
    tens[1, 1] = 1 
    return tens

function bulk_tensor(dimension)
    tens = zeros(1, 1, dimension)
    tens[1, 1, 1] = 1 
    return tens

function gs_product_state(s, N)

    # Initial state's MPS representation
    α = [Index(1, "α$n") for n in 1:(N-1)] # bond indices

    ψgs = MPS(N)
    for n in 1:N
        if n == 1 # the first ion
            ψgs[n] = ITensor(surface_tensor(dim(s[n])), α[n], s[n])
        elseif n == N # the last ion
            ψgs[n] = ITensor(surface_tensor(dim(s[n])), α[n - 1], s[n])
        else # all the other sites
            ψgs[n] = ITensor(bulk_tensor(dim(s[n])), α[n - 1], α[n], s[n])
    return ψgs

#MPOs generation
σx = [0 1;1 0]
I2 = [1 0;0 1]

function σx_onsite(N, s, site_ind)
    # Apply σx to the spin on site numbered by site_ind.
    gate = MPO(s, [if n==site_ind σx else I2 end for n in 1:N])
    return gate

# concatenate N gates with sequential σxs:
N = 2
s = [Index(2, "s$n") for n in 1:N] # site indices
ψ0 = gs_product_state(s, N) # Initialize the MPS

gates = ITensor[]
for n in 1:N
    gate = σx_onsite(N, s(n), n)
    append!(gates, gate)

ψ = apply(gates, ψ0)


Hi Yotam,
Is your question about how to use ITensor to multiply an MPO times an MPS? Because if you know how to do that correctly, you can apply multiple MPO’s to an MPS by repeating this operation.

Thank you for the question.

MPO-MPS multiplication works well for me, but it works sequentially.
I wondered if I could gain a speed-up from pre-generating a tensor network built up from multiple MPOs acting on each other and finally on the MPS?
Would that be more efficient than contracting each MPO at a time?

For instance, I might contract simultaneously each of two adjacent MPOs and then contract the rest.
Is there a pre-built algorithm designed for such problems?
Can the GPU package figure this out in a more performant manner?

Thanks again,