MPS with 4 indices per site

I am trying to create an MPS that has 4 distinct indices per site. I need to contract individually so simply vectorizing these 4 indices into 1 (larger) index is not a solution. Is there already something built in that does this? I tried the following, but not sure whether this is the best thing:

# creates an MPS with 4 indices per site
using ITensors

function superMPO(sites::Vector{<:Index}, linkdims::Integer=1)
    T = ComplexF64
    N = length(sites)
    v = Vector{ITensor}(undef, N)
    if N == 1
        s = sites[1]
        v[1] = ITensor(T, s, s', s'', s''')
        return MPS(v)

    space = if hasqns(sites)
        [QN() => linkdims]

    l = [Index(space, "Link,l=$ii") for ii in 1:(N - 1)]
    for ii in eachindex(sites)
        s = sites[ii]
        if ii == 1
            v[ii] = ITensor(T, l[ii], s, s', s'', s''')
        elseif ii == N
            v[ii] = ITensor(T, dag(l[ii - 1]), s, s', s'', s''')
            v[ii] = ITensor(T, dag(l[ii - 1]), s, s', s'', s''', l[ii])
    return MPS(v)

function main()
    d = 2
    N = 5
    sites = siteinds(d, N)
    A1 = rand(d*ones(Int,N)...)
    A2 = rand(d*ones(Int,4*N)...)

    M1 = MPS(A1, sites)
    M2 = superMPO(sites)

    @show M1
    @show M2



I took inspiration from the implementation of the MPS constructor (compare ITensors.jl/mps.jl at efacea7b2e5bc1894b58ce75874352d825eb8954 · ITensor/ITensors.jl · GitHub).

In addition, what function exactly does MPS(A1, sites) call? I was not able to find it in the source code.

Good question, but no, there is not a built-in function that makes a generalized MPO with four indices per site. Usually with tensor networks one tries to keep the number of indices on each tensor low, since the number of entries of a tensor scales exponentially in the number of parameters, and your tensors here will have six indices.

So for example, instead of having four external indices, you could have an MPS with a “unit cell” of four sites that repeat like s1, s1’, s1’‘, s1’‘’, s2, s2’, s2’‘, s2’‘’. Just an idea.

As for where a function is in the ITensor source, Julia offers a really helpful macro called @which that you can use like @which MPS(A1,sites). On my machine, it gives the following result:

(::Type{MPST})(A::ITensor, sites; leftinds, orthocenter, kwargs...) where MPST<:ITensors.AbstractMPS in ITensors at /Users/mstoudenmire/.julia/dev/ITensors/src/mps/abstractmps.jl:1860

So you can see that the function called is on line 1860 of the file src/mps/abstractmps.jl of the ITensors.jl source code.

1 Like

Thanks for your reply. So you suggest instead of using a ‘superMPO’ of length N with 4 indices per site, I could also use an MPS with length 4*N? Here’s a picture for N=2.

I am wondering, does this not lead to a higher dimensional object since we also have the bond indices between, say, i_1 and i_1'?

And cool feature, didn’t know about @which.

Good question. I would say it is not a higher dimensional object, since the number of external or physical indices is the same. It is a longer object in the sense that there are more MPS tensors to store, but that is only a small linear-scaling cost.

The main thing to be wary of is tensors which have too many indices, since the number of elements of a tensor grows exponentially in the number of its indices. So having 4 physical indices plus 2 link indices for 6 total indies is a lot!

I think the MPS version you drew at the bottom will probably be more efficient. It might be slightly less convenient in some way to work with, though, depending on what you want to do. But say if you are applying two-site gates, then our apply function automatically handles further-neighbor gates for you so you wouldn’t have to worry about the technicalities of doing that.