Problem of projector calculation for a complex MPS

Hello everyone!

I am interested in calculating from a \ket{\psi} state the projector associated with that state i.e. \ket{\psi}\bra{\psi}. However, I get a solution that seems incorrect when I perform this operation for an MPS \ket{\psi} containing complex coefficients. When I calculate the projector associated with the state I get an MPO. When I contract that MPO to see what is the rank of the matrix I obtain that the solution is not 1. This behavior only happens when the MPS \ket{\psi} that I define has complex coefficients, if it has only real coefficients then it works. Attached is the code to reproduce the results:

# library

using ITensors # version v0.7.13
using ITensorMPS # version v0.3.6
using LinearAlgebra 

# number of sites

N = 2 
sites = siteinds("S=1/2", 2)  

# create a normalize random real MPS state 

psi = randomMPS(sites) 

psi /= sqrt(inner(psi,psi))

# calculate the associated projector |psi><psi|

P = projector(psi; normalize = false)

# contract the MPO obtained and shape it into a matrix.

rank(reshape(array(contract(P)), (4,4)))

When the rank of the matrix is calculated, it is observed that its value is 1, which is the value that I expect. When we perform the same procedure for complex MPSs the result is different.

# library

using ITensors # version v0.7.13
using ITensorMPS # version v0.3.6
using LinearAlgebra 

# number of sites

N = 2 
sites = siteinds("S=1/2", 2)  

# create a normalize random complex MPS state 

coefficient = 2im

psi = randomMPS(sites)  + coefficient * randomMPS(sites)

psi /= sqrt(inner(psi,psi))

# calculate the associated projector |psi><psi|

P = projector(psi; normalize = false)

# contract the MPO obtained and shape it into a matrix.

rank(reshape(array(contract(P)), (4,4)))

In this case the rank of the matrix is 4 and not 1, so I am not getting the projector in the right way.

I would like to know what is going on, is there a possible bug in the functions I am using or am I doing something wrong?

1 Like

I think you can’t assume the MPO indices naturally form the matrix you’re after. If you call array with the correct order I get rank 1

julia> rank(reshape(array(contract(P), sites...,sites'...),(4,4)))
1

You could use a combiner to make this clearer without a reshape

julia> Callinds = combiner(sites)
ITensor ord=3 (dim=4|id=804|"CMB,Link") (dim=2|id=814|"S=1/2,Site,n=1") (dim=2|id=967|"S=1/2,Site,n=2")
NDTensors.Combiner

julia> combined_P = contract(P)*Callinds*Callinds'
ITensor ord=2 (dim=4|id=804|"CMB,Link") (dim=4|id=804|"CMB,Link")'
NDTensors.Dense{ComplexF64, Vector{ComplexF64}}

julia> c = combinedind(Callinds)
(dim=4|id=804|"CMB,Link")

julia> rank(array(combined_P,c,c'))
1

Thank you very much, that fixed my mistake, I was not considering the correct index contraction.

1 Like