Unexpected error with number-conserving MPO when result is zero

ITensor Julia Questions
, ,
1

The following code gives an error which I don’t understand.
I set up two different states with number conservation active:

using ITensors, ITensorMPS
N = 3
s = siteinds("Electron", N; conserve_nf=true)
x = MPS(s, ["Emp", "Emp", "Emp"])
y = MPS(s, ["Emp", "Emp", "Up"])
ops = OpSum()
ops += "c†↑", 1, "c↑", 3
A = MPO(ops, s)

Now, if I run

apply(A, x)

I get the following error:

ERROR: ArgumentError: collection must be non-empty
Stacktrace:
  [1] first
    @ ./abstractarray.jl:473 [inlined]
  [2] eigen(T::LinearAlgebra.Hermitian{…}; min_blockdim::Nothing, mindim::Int64, maxdim::Int64, cutoff::Float64, use_absolute_cutoff::Nothing, use_relative_cutoff::Nothing)
    @ NDTensors /opt/julia/packages/NDTensors/Q9Zpu/src/blocksparse/linearalgebra.jl:233
  [3] eigen(A::ITensor, Linds::Vector{…}, Rinds::Vector{…}; mindim::Int64, maxdim::Int64, cutoff::Float64, use_absolute_cutoff::Nothing, use_relative_cutoff::Nothing, ishermitian::Bool, tags::ITensors.TagSets.GenericTagSet{…}, lefttags::Nothing, righttags::Nothing, plev::Nothing, leftplev::Nothing, rightplev::Nothing)
    @ ITensors /opt/julia/packages/ITensors/AOPBm/src/tensor_operations/matrix_decomposition.jl:378
  [4] contract(::NDTensors.BackendSelection.Algorithm{…}, A::MPO, ψ::MPS; cutoff::Float64, maxdim::Int64, mindim::Int64, normalize::Bool, kwargs::@Kwargs{})
    @ ITensorMPS /opt/julia/packages/ITensorMPS/3OORN/src/mpo.jl:743
  [5] contract(::NDTensors.BackendSelection.Algorithm{:densitymatrix, @NamedTuple{}}, A::MPO, ψ::MPS)
    @ ITensorMPS /opt/julia/packages/ITensorMPS/3OORN/src/mpo.jl:692
  [6] product(alg::NDTensors.BackendSelection.Algorithm{:densitymatrix, @NamedTuple{}}, A::MPO, ψ::MPS; kwargs::@Kwargs{})
    @ ITensorMPS /opt/julia/packages/ITensorMPS/3OORN/src/mpo.jl:605
  [7] product
    @ /opt/julia/packages/ITensorMPS/3OORN/src/mpo.jl:604 [inlined]
  [8] #apply#453
    @ /opt/julia/packages/ITensorMPS/3OORN/src/mpo.jl:601 [inlined]
  [9] product(A::MPO, ψ::MPS)
    @ ITensorMPS /opt/julia/packages/ITensorMPS/3OORN/src/mpo.jl:600
 [10] top-level scope
    @ REPL[42]:1

while the following

apply(A, y)

works fine. I guess that the difference in the behaviour lies in the fact that Ay = 0 whereas Ax \neq 0, but why does it throw an error? The c_1^\dagger c_3 operator preserves the number of particles so applying it to x or y shouldn’t be an issue.

When QNs are not active, both apply(A,x) and apply(A,y) throw no errors, with \lVert Ay\rVert = 0 as expected.

2

On our side, we should improve this error message.

For your application, when your calculation results in a state with zero norm, further down in our system, a tensor with no blocks getting passed into eigen. In that situation, we don’t have a good way, even mathematically, to know or guess what the QN flux of the result should be. (We could also in the future add a system to ‘hint’ the QN of the result to work around cases like this but probably that’s too fancy.)

Do you need to be able to do apply(A,y) when the result will definitely end up being zero? I think at an outer level, if this is happening due to the QN sector of y, you could detect this by calling flux(y) and seeing if it’s even possible for Ay to be non-zero.

3

I see. It does make sense that the zero vector behaves as you describe, unfortunately.

My goal was to replicate ITensor’s correlation_function method for an MPS that does not represent a pure state but a mixed one, as in “Eigen currently only supports block diagonal matrices.” when applying an MPO to a MPS with QNs. For this kind of object I defined custom methods to compute expectation values, using the trace and such instead of the usual inner product. The correlation matrix would be defined using these new methods.

Since I’m dealing with fermions, I use the OpSum → MPO method to construct the observable so that I don’t have to put in Jordan-Wigner strings by hand. That’s why I’m using apply in this situation.