`MPO` of a matrix doubling my matrix dimensions?

ITensor Julia Questions
1

I posted a week ago about comparing ED vs DMRG, which I’m still doing.
What I want to do now is study the magnetization of the spin chain, and then compare it with the already given MPO of the “Sz” operator. The MWE for this bug

using ITensors, ITensorMPS
function op_matrix(L; op=[ 1 0 ; 0 -1 ])
    Sz = fill(ID, L)
    Sz[1] = op

    M = zeros(Int, 2^L, 2^L)
    for i=1:L
        M += foldl(kron, Sz)
        Sz = circshift(Sz, 1)
    end
    M/L
end   

L = 3
sites = siteinds("S=1/2", L)
Mz_matrix = op_matrix(L) # Returns an 2³ x 2³ matrix
Mz_MPO = MPO(Mz_matrix, sites)

Yields ERROR: DimensionMismatch: In ITensor(::AbstractArray, inds), length of AbstractArray (64) must match total dimension of IndexSet (8). For some reason MPO is doubling the input matrix.

Does anyone have any idea why that could be?

2

(Unfortunately I can’t run your code because ID is not defined. But I think it’s just a stand-in for arrays that get overwritten so I changed it to zeros(2,2).)

I see where the error is coming from: if you look at line 2006 of abstractmps.jl in the ITensorMPS source, what that constructor does is forward the array and the array of indices (in your case, sites) to the ITensor constructor that makes an ITensor from an array and a collection of indices. That leads to a somewhat awkward interface here, and anyway that code path seems to be broken by another function later on.

Here is a workaround you can try:

  1. make an array of length 2L consisting of sites and primed sites, something like
mpo_sites = [ [s,s'] for s in sites]
mpo_sites = collect(Iterators.flatten(mpo_sites))
  1. turn your array into an ITensor (this is what that MPO constructor was trying to do anyway)
  A = ITensor(Mz_matrix,mpo_sites)
  1. construct the MPO from the ITensor and original sites array
  Mz_MPO = MPO(A, sites)

You may have to play around with the ordering of the indices in mpo_sites to get the mapping from the array to the ITensor to be correct.