linsolve unstable under repeated application

Previously posted on the wrong Issues page

Hi,

I've been trying to implement a power iteration algorithm using the linsolve function. Here's a minimum example:
That is, I'm trying to repeatedly apply: solving (H - λ) ψ = ϕ, ψ → ϕ, i.e. solving (H - λ)^ (-N) ϕ

using ITensors
using ITensorTDVP: linsolve

function test()
    L = 12
    N = 6
    t = 1.0

    λ = - 10.0
    sites = siteinds("Fermion", L, conserve_qns=true)
    ampo = OpSum()

    for i=1:L-1
        ampo += -t , "C",i,"Cdag",i+1
        ampo += -t , "C",i+1,"Cdag",i
    end 

    H = MPO(ampo,sites)
    state = append!([ "Occ" for n=1:N] , ["Emp" for n=1:L -N])
    ϕ = randomMPS(sites,state)

    # outer iteration, solving excited state eigen energy immediately above λ
    for ex=1:10

        # power iteration, repeatedly solving (H - λ) ψ = ϕ, ψ → ϕ
        for cnt = 1:100

            # generate an initial guess ψ0
            #state = append!([ "Occ" for n=1:N] , ["Emp" for n=1:L -N])
            #ψ0 = randomMPS(sites,state)
            ψ0 = ϕ
            ψ = linsolve(H, ϕ, ψ0, -λ, 1.0)
            ϕ = ψ

            energy = inner(ϕ', H, ϕ) / inner( ϕ', ϕ)
            println("energy = ", energy)

        end 

        λ = energy
        println("λ =  ", λ)
    end 

end 

test()

In practice, this power iteration algorithm is extremely unstable, as I frequently encounter the following errors:

ERROR: BoundsError: attempt to access 0-element Vector{Pair{QN, Int64}} at index [1]
Stacktrace:
  [1] getindex
    @ ./array.jl:861 [inlined]
  [2] combineblocks(qns::Vector{Pair{QN, Int64}})
    @ ITensors ~/.julia/packages/ITensors/OjQuG/src/qn/qnindex.jl:407
  [3] combineblocks
    @ ~/.julia/packages/ITensors/OjQuG/src/qn/qnindex.jl:469 [inlined]
  [4] combiner(inds::Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ITensors ~/.julia/packages/ITensors/OjQuG/src/qn/qnitensor.jl:364
  [5] combiner
    @ ~/.julia/packages/ITensors/OjQuG/src/qn/qnitensor.jl:360 [inlined]
  [6] #combiner#216
    @ ~/.julia/packages/ITensors/OjQuG/src/itensor.jl:1575 [inlined]
  [7] combiner(::Index{Vector{Pair{QN, Int64}}}, ::Index{Vector{Pair{QN, Int64}}})
    @ ITensors ~/.julia/packages/ITensors/OjQuG/src/itensor.jl:1575
  [8] svd(A::ITensor, Linds::Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}; kwargs::Base.Pairs{Symbol, Any, NTuple{10, Symbol}, NamedTuple{(:which_decomp, :tags, :maxdim, :mindim, :cutoff, :eigen_perturbation, :ortho, :normalize, :svd_alg, :alg), Tuple{Nothing, TagSet, Int64, Int64, Float64, Nothing, String, Bool, String, String}}})
    @ ITensors ~/.julia/packages/ITensors/OjQuG/src/decomp.jl:112
  [9] factorize_svd(A::ITensor, Linds::Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}; kwargs::Base.Pairs{Symbol, Any, NTuple{9, Symbol}, NamedTuple{(:which_decomp, :tags, :maxdim, :mindim, :cutoff, :eigen_perturbation, :ortho, :normalize, :svd_alg), Tuple{Nothing, TagSet, Int64, Int64, Float64, Nothing, String, Bool, String}}})
    @ ITensors ~/.julia/packages/ITensors/OjQuG/src/decomp.jl:404
 [10] factorize(A::ITensor, Linds::Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}; kwargs::Base.Pairs{Symbol, Any, NTuple{9, Symbol}, NamedTuple{(:which_decomp, :tags, :maxdim, :mindim, :cutoff, :eigen_perturbation, :ortho, :normalize, :svd_alg), Tuple{Nothing, TagSet, Int64, Int64, Float64, Nothing, String, Bool, String}}})
    @ ITensors ~/.julia/packages/ITensors/OjQuG/src/decomp.jl:524
 [11] replacebond!(M::MPS, b::Int64, phi::ITensor; kwargs::Base.Pairs{Symbol, Any, NTuple{8, Symbol}, NamedTuple{(:maxdim, :mindim, :cutoff, :eigen_perturbation, :ortho, :normalize, :which_decomp, :svd_alg), Tuple{Int64, Int64, Float64, Nothing, String, Bool, Nothing, String}}})
    @ ITensors ~/.julia/packages/ITensors/OjQuG/src/mps/mps.jl:457
 [12] tdvp_site_update!(nsite_val::Val{2}, reverse_step_val::Val{false}, solver::ITensorTDVP.var"#linsolve_solver#68"{ITensorTDVP.var"#linsolve_solver#67#69"{Float64, Float64}}, PH::ITensorTDVP.ProjMPO_MPS2, psi::MPS, b::Int64; current_time::Float64, outputlevel::Int64, time_step::Float64, normalize::Bool, direction::Base.Order.ForwardOrdering, noise::Float64, which_decomp::Nothing, svd_alg::String, cutoff::Float64, maxdim::Int64, mindim::Int64, maxtruncerr::Float64)
    @ ITensorTDVP ~/.julia/packages/ITensorTDVP/p6CK4/src/tdvp_step.jl:311
 [13] tdvp_site_update!(solver::ITensorTDVP.var"#linsolve_solver#68"{ITensorTDVP.var"#linsolve_solver#67#69"{Float64, Float64}}, PH::ITensorTDVP.ProjMPO_MPS2, psi::MPS, b::Int64; nsite::Int64, reverse_step::Bool, current_time::Float64, outputlevel::Int64, time_step::Float64, normalize::Bool, direction::Base.Order.ForwardOrdering, noise::Float64, which_decomp::Nothing, svd_alg::String, cutoff::Float64, maxdim::Int64, mindim::Int64, maxtruncerr::Float64)
    @ ITensorTDVP ~/.julia/packages/ITensorTDVP/p6CK4/src/tdvp_step.jl:155
 [14] tdvp_sweep(direction::Base.Order.ForwardOrdering, solver::Function, PH::ITensorTDVP.ProjMPO_MPS2, time_step::Float64, psi::MPS; kwargs::Base.Pairs{Symbol, Real, NTuple{7, Symbol}, NamedTuple{(:current_time, :reverse_step, :sweep, :maxdim, :mindim, :cutoff, :noise), Tuple{Float64, Bool, Int64, Int64, Int64, Float64, Float64}}})
    @ ITensorTDVP ~/.julia/packages/ITensorTDVP/p6CK4/src/tdvp_step.jl:79
 [15] tdvp_step(order::ITensorTDVP.TDVPOrder{2, Base.Order.ForwardOrdering()}, solver::Function, PH::ITensorTDVP.ProjMPO_MPS2, time_step::Float64, psi::MPS; current_time::Float64, kwargs::Base.Pairs{Symbol, Real, NTuple{6, Symbol}, NamedTuple{(:reverse_step, :sweep, :maxdim, :mindim, :cutoff, :noise), Tuple{Bool, Int64, Int64, Int64, Float64, Float64}}})
    @ ITensorTDVP ~/.julia/packages/ITensorTDVP/p6CK4/src/tdvp_step.jl:9
 [16] macro expansion
    @ ~/.julia/packages/ITensorTDVP/p6CK4/src/tdvp_generic.jl:84 [inlined]
 [17] macro expansion
    @ ./timing.jl:299 [inlined]
 [18] tdvp(solver::Function, PH::ITensorTDVP.ProjMPO_MPS2, t::Float64, psi0::MPS; kwargs::Base.Pairs{Symbol, Bool, Tuple{Symbol}, NamedTuple{(:reverse_step,), Tuple{Bool}}})
    @ ITensorTDVP ~/.julia/packages/ITensorTDVP/p6CK4/src/tdvp_generic.jl:83
 [19] linsolve(A::MPO, b::MPS, x₀::MPS, a₀::Float64, a₁::Float64; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ITensorTDVP ~/.julia/packages/ITensorTDVP/p6CK4/src/linsolve.jl:47
 [20] linsolve(A::MPO, b::MPS, x₀::MPS, a₀::Float64, a₁::Float64)
    @ ITensorTDVP ~/.julia/packages/ITensorTDVP/p6CK4/src/linsolve.jl:22
 [21] test()
    @ Main path-to-script/problem.jl:33
 [22] top-level scope
    @ path-to-script/problem.jl:47

Of course the severity of this problem varies with the initial guesses, however I would almost always encounter this issue somewhere down the iterations.

I'm assuming this error occurs because the variational linear solver locally encountered a matrix with no solution.

Start by reproducing the example and tracing linsolve in linsolve.jl:47 into tdvp_step.jl:311, where the failure reaches qnindex.jl:407 through SVD and bond replacement. Determine why repeated applications produce an empty QN block, then make the solver handle the case without the reported BoundsError. Done means the example no longer crashes under repeated application and the relevant behavior is covered by a regression test.