linsolve unstable under repeated application

[keyi-ray-liu]

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.