Site ordering for MPS

ITensor Julia Questions
, ,
1

Hi there,

Julia and ITensor version info:

julia> versioninfo()
Julia Version 1.10.2
Commit bd47eca2c8a (2024-03-01 10:14 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin22.4.0)
  CPU: 8 × Apple M1
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, apple-m1)
Threads: 1 default, 0 interactive, 1 GC (on 4 virtual cores)

julia> using Pkg; Pkg.status("ITensors")
Status `~/.julia/environments/v1.10/Project.toml`
  [9136182c] ITensors v0.6.21

Apologies for asking a possibly naive question. I read the MPS ITensor examples in the docs and a few forum questions/responses to put together the following code for calculating the von Neumann entropy of an MPS partitioned between site b and b+1.

# Import packages
using ITensors
using ITensorMPS

# Define function to calculate the entanglement entropy of MPS psi partitioned betweem site b and b+1
function MPS_entanglement_entropy(psi::MPS, b::Int)
    # Extract the site indices of psi
    s = siteinds(psi)
    # Set the orthogonal centre of the MPS to site b
    psi =  orthogonalize!(psi, b)
    # Extract the bond index between sites b-1 and b
    l = linkind(psi, b-1)
    # Define ls as the set of indices to appear as the rows in the SVD (if loop for consistency with b=1)
    ls = isnothing(l) ? (s[b],) : (l, s[b])
    # Extract singular values w.r.t chosen grouping of legs from MPS tensor at site b
    _,S = svd(psi[b], ls)
    # Calculate and return von Neumann entropy
    SvN = 0.0
    for n=1:dim(S, 1)
        p = S[n,n]^2
        if p > 0
            SvN -= p * log(2,p)    
        end
    end    
    return SvN
end

To test the code I tried printing the output for some standard pure states with all possible values of b:

# Calculate entanglement entropies for |0000>+|1100> (Bell-state otimes |00>)
let
    A = [ 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 ]
    A = A/norm(A) # Normalise state
    sites = siteinds("Qubit",4)
    cutoff = 1E-8
    maxdim = 10
    M = MPS(A,sites;cutoff=cutoff,maxdim=maxdim)
    println("Results for |0000>+|1100>:")
    println("SvN for b=1: ",round(MPS_entanglement_entropy(M,1),digits=8))
    println("SvN for b=2: ",round(MPS_entanglement_entropy(M,2),digits=8))
    println("SvN for b=3: ",round(MPS_entanglement_entropy(M,3),digits=8))
    println()
end

# Calculate entanglement entropies for |0000>+|0011> (|00> otimes Bell-state)
let
    A = [ 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 ]
    A = A/norm(A) # Normalise state
    sites = siteinds("Qubit",4)
    cutoff = 1E-8
    maxdim = 10
    M = MPS(A,sites;cutoff=cutoff,maxdim=maxdim)
    println("Results for |0000>+|0011>:")
    println("SvN for b=1: ",round(MPS_entanglement_entropy(M,1),digits=8))
    println("SvN for b=2: ",round(MPS_entanglement_entropy(M,2),digits=8))
    println("SvN for b=3: ",round(MPS_entanglement_entropy(M,3),digits=8))
    println()
end

This gives me the output:

Results for |0000>+|1100>:
SvN for b=1: 0.0
SvN for b=2: 0.0
SvN for b=3: 1.0

Results for |0000>+|0011>:
SvN for b=1: 1.0
SvN for b=2: -0.0
SvN for b=3: -0.0

I would have expected the results for the first state to be 1,0,0, and for the second to be 0,0,1. This is because I assumed that b=1 refers to the left-most MPS site, and that the standard basis vectors are ordered like {|0000>,|0001>,|0010>,|0011>…} when specifying A. Obviously the above output means I have misunderstood something, are you able to help me make sense of this? E.g. maybe b=1 refers to the rightmost MPS site, or maybe the basis vectors are ordered like {|0000>,|1000>,|0100>,|1100>,…}.

Thanks very much!

2

This is an interesting “code experiment”.

I would guess most likely from the fact that the results are reversed from what you’d expect that the issue is with the input arrays called “A”. Perhaps you are inputting the amplitudes in a kind of ‘row-major’ order whereas our system expects a ‘column-major’ order. That is, please test on a two-qubit example that your understanding of how the entries of the array map to the entries of the MPS is correct.

3

Thanks for the advice Miles! Indeed I misinterpreted how the entries of A were being placed into the MPS. In case others have the same issue, I think the following code with 3 MPS sites makes it clear that the basis is ordered with the least significant bit on the left {|000>,|100>,|010>,|110>,…} when defining A (I suppose this is a result of ‘column-major’ order):

let
    N = 3
    A = [ 1 2 3 4 5 6 7 8 ]
    sites = siteinds("Qubit",N)
    cutoff = 1E-8
    maxdim = 10
    M = MPS(A,sites;cutoff=cutoff,maxdim=maxdim)
    println("Results for A = [1,2,3,4,5,6,7,8]:")
    for i3=1:2
        for i2=1:2
            for i1=1:2
                index = [i1,i2,i3]
                V = ITensor(1.)
                for j=1:N
                  V *= (M[j]*state(sites[j],index[j]))
                end
                v = scalar(V)
                println("MPS element [",i1-1,",",i2-1,",",i3-1,"] is: ",round(v,digits=1))
            end
        end
    end
end

which gives output

Results for A = [1,2,3,4,5,6,7,8]:
MPS element [0,0,0] is: 1.0
MPS element [1,0,0] is: 2.0
MPS element [0,1,0] is: 3.0
MPS element [1,1,0] is: 4.0
MPS element [0,0,1] is: 5.0
MPS element [1,0,1] is: 6.0
MPS element [0,1,1] is: 7.0
MPS element [1,1,1] is: 8.0

Thanks again for the help :slight_smile:

4

Glad that solved it!