MPSKit.jl icon indicating copy to clipboard operation
MPSKit.jl copied to clipboard

Published 20 hours ago • verified publisher icon maartenvd

negative Grassmann inner product

Open lkdvos opened this issue 2 years ago • 2 comments

Due to numerical errors, it seems like the inner product of a grassman gradient with itself can sometimes take on negative values that are of the order -eps, which then fails when OptimKit wants to take the square root.

MBA:

using MPSKit, MPSKitModels
psi = InfiniteMPS(2, 12)
H = transverse_field_ising(; hx=0.0) # exact product state solution

julia> psi, envs, delta = find_groundstate(psi, H)
[ Info: vumps @iteration 1 galerkin = 0.5339115698990319
[ Info: vumps @iteration 2 galerkin = 0.0012648372627602012
[ Info: vumps @iteration 3 galerkin = 1.1562975109758041e-11
ERROR: DomainError with -1.413494684364987e-21:
sqrt will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
Stacktrace:
 [1] throw_complex_domainerror(f::Symbol, x::Float64)
   @ Base.Math ./math.jl:33
 [2] sqrt
   @ ./math.jl:591 [inlined]
 [3] optimize(fg::typeof(MPSKit.GrassmannMPS.fg), x::MPSKit.GrassmannMPS.ManifoldPoint{InfiniteMPS{TensorKit.TensorMap{TensorKit.ComplexSpace, 2, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}, TensorKit.TensorMap{TensorKit.ComplexSpace, 1, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}}, MPSKit.MPOHamInfEnv{MPOHamiltonian{TensorKit.ComplexSpace, TensorKit.TensorMap{TensorKit.ComplexSpace, 2, 2, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}, ComplexF64}, TensorKit.TensorMap{TensorKit.ComplexSpace, 2, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}, InfiniteMPS{TensorKit.TensorMap{TensorKit.ComplexSpace, 2, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}, TensorKit.TensorMap{TensorKit.ComplexSpace, 1, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}}, KrylovKit.GMRES{KrylovKit.ModifiedGramSchmidt2, Float64}}, PeriodicArray{TensorKitManifolds.Grassmann.GrassmannTangent{TensorKit.TensorMap{TensorKit.ComplexSpace, 2, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}, TensorKit.TensorMap{TensorKit.ComplexSpace, 2, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}, TensorKit.TensorMap{TensorKit.ComplexSpace, 1, 1, TensorKit.Trivial, Matrix{Float64}, Nothing, Nothing}, TensorKit.TensorMap{TensorKit.ComplexSpace, 1, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}}, 1}, Vector{TensorKit.TensorMap{TensorKit.ComplexSpace, 1, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}}}, alg::OptimKit.ConjugateGradient{OptimKit.HagerZhang{Rational{Int64}}, Float64, OptimKit.HagerZhangLineSearch{Rational{Int64}}}; precondition::typeof(OptimKit._precondition), finalize!::typeof(OptimKit._finalize!), retract::Function, inner::typeof(MPSKit.GrassmannMPS.inner), transport!::typeof(MPSKit.GrassmannMPS.transport!), scale!::typeof(MPSKit.GrassmannMPS.scale!), add!::typeof(MPSKit.GrassmannMPS.add!), isometrictransport::Bool)
   @ OptimKit ~/.julia/packages/OptimKit/xpmbV/src/cg.jl:27
 [4] find_groundstate(Ψ::InfiniteMPS{TensorKit.TensorMap{TensorKit.ComplexSpace, 2, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}, TensorKit.TensorMap{TensorKit.ComplexSpace, 1, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}}, H::MPOHamiltonian{TensorKit.ComplexSpace, TensorKit.TensorMap{TensorKit.ComplexSpace, 2, 2, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}, ComplexF64}, alg::GradientGrassmann, envs::MPSKit.MPOHamInfEnv{MPOHamiltonian{TensorKit.ComplexSpace, TensorKit.TensorMap{TensorKit.ComplexSpace, 2, 2, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}, ComplexF64}, TensorKit.TensorMap{TensorKit.ComplexSpace, 2, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}, InfiniteMPS{TensorKit.TensorMap{TensorKit.ComplexSpace, 2, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}, TensorKit.TensorMap{TensorKit.ComplexSpace, 1, 1, TensorKit.Trivial, Matrix{ComplexF64}, Nothing, Nothing}}, KrylovKit.GMRES{KrylovKit.ModifiedGramSchmidt2, Float64}})
   @ MPSKit ~/Projects/JuliaProjects/MPSKit.jl/src/algorithms/groundstate/gradient_grassmann.jl:53
 [5] find_groundstate
   @ ~/Projects/JuliaProjects/MPSKit.jl/src/algorithms/unionalg.jl:25 [inlined]
 [6] #find_groundstate#381
   @ ~/Projects/JuliaProjects/MPSKit.jl/src/algorithms/groundstate/find_groundstate.jl:41 [inlined]
 [7] find_groundstate (repeats 2 times)
   @ ~/Projects/JuliaProjects/MPSKit.jl/src/algorithms/groundstate/find_groundstate.jl:19 [inlined]
 [8] top-level scope
   @ REPL[12]:1

lkdvos avatar May 06 '23 16:05 lkdvos

So I've been trying to fix this bug by implementing what we discussed together, namely:

function inner(g1::PrecGrad, g2::PrecGrad, rho=one(g1.rho))
   
    if g1 === g2
        max(real(unsafe_inner(g1,g2,rho)),eps())
    else   
        unsafe_inner(g1,g2,rho)
    end
end
function unsafe_inner(g1::PrecGrad, g2::PrecGrad, rho=one(g1.rho))
    Grassmann.base(g1) == Grassmann.base(g2) || throw(ArgumentError("incompatible base"))
    if g1.rho == rho
        Grassmann.dot(g1.g.Z, g2.Pg.Z)
    elseif g2.rho == rho
        Grassmann.dot(g1.Pg.Z, g2.g.Z)
    else
        Grassmann.dot(g1.Pg.Z, g2.Pg.Z * rho)
    end

end

This fixes the bug for the sqrt but introduces another later down the line

ERROR: LoadError: linesearch was not given a descent direction!
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] (::OptimKit.HagerZhangLineSearch{Rational{Int64}})(fg::Function, x₀::MPSKit.GrassmannMPS.ManifoldPoint{InfiniteMPS{TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 2, 1, Matrix{ComplexF64}}, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 1, 1, Matrix{ComplexF64}}}, MPSKit.MPOHamInfEnv{MPOHamiltonian{TensorKit.ComplexSpace, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 2, 2, Matrix{ComplexF64}}, ComplexF64}, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 2, 1, Matrix{ComplexF64}}, InfiniteMPS{TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 2, 1, Matrix{ComplexF64}}, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 1, 1, Matrix{ComplexF64}}}, KrylovKit.GMRES{KrylovKit.ModifiedGramSchmidt2, Float64}}, PeriodicArray{TensorKitManifolds.Grassmann.GrassmannTangent{TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 2, 1, Matrix{ComplexF64}}, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 2, 1, Matrix{ComplexF64}}, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 1, 1, Matrix{Float64}}, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 1, 1, Matrix{ComplexF64}}}, 1}, Vector{TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 1, 1, Matrix{ComplexF64}}}}, η₀::Vector{MPSKit.GrassmannMPS.PrecGrad{TensorKitManifolds.Grassmann.GrassmannTangent{TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 2, 1, Matrix{ComplexF64}}, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 2, 1, Matrix{ComplexF64}}, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 1, 1, Matrix{Float64}}, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 1, 1, Matrix{ComplexF64}}}, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 1, 1, Matrix{ComplexF64}}}}, fg₀::Tuple{Float64, Vector{MPSKit.GrassmannMPS.PrecGrad{TensorKitManifolds.Grassmann.GrassmannTangent{TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 2, 1, Matrix{ComplexF64}}, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 2, 1, Matrix{ComplexF64}}, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 1, 1, Matrix{Float64}}, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 1, 1, Matrix{ComplexF64}}}, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 1, 1, Matrix{ComplexF64}}}}}; retract::Function, inner::typeof(MPSKit.GrassmannMPS.inner), initialguess::Float64, acceptfirst::Bool)
   @ OptimKit ~/.julia/packages/OptimKit/xpmbV/src/linesearches.jl:255
 [3] optimize(fg::typeof(MPSKit.GrassmannMPS.fg), x::MPSKit.GrassmannMPS.ManifoldPoint{InfiniteMPS{TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 2, 1, Matrix{ComplexF64}}, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 1, 1, Matrix{ComplexF64}}}, MPSKit.MPOHamInfEnv{MPOHamiltonian{TensorKit.ComplexSpace, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 2, 2, Matrix{ComplexF64}}, ComplexF64}, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 2, 1, Matrix{ComplexF64}}, InfiniteMPS{TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 2, 1, Matrix{ComplexF64}}, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 1, 1, Matrix{ComplexF64}}}, KrylovKit.GMRES{KrylovKit.ModifiedGramSchmidt2, Float64}}, PeriodicArray{TensorKitManifolds.Grassmann.GrassmannTangent{TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 2, 1, Matrix{ComplexF64}}, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 2, 1, Matrix{ComplexF64}}, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 1, 1, Matrix{Float64}}, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 1, 1, Matrix{ComplexF64}}}, 1}, Vector{TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 1, 1, Matrix{ComplexF64}}}}, alg::OptimKit.ConjugateGradient{OptimKit.HagerZhang{Rational{Int64}}, Float64, OptimKit.HagerZhangLineSearch{Rational{Int64}}}; precondition::typeof(OptimKit._precondition), finalize!::typeof(OptimKit._finalize!), retract::Function, inner::typeof(MPSKit.GrassmannMPS.inner), transport!::typeof(MPSKit.GrassmannMPS.transport!), scale!::typeof(MPSKit.GrassmannMPS.scale!), add!::typeof(MPSKit.GrassmannMPS.add!), isometrictransport::Bool)
   @ OptimKit ~/.julia/packages/OptimKit/xpmbV/src/cg.jl:68
 [4] find_groundstate(ψ::InfiniteMPS{TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 2, 1, Matrix{ComplexF64}}, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 1, 1, Matrix{ComplexF64}}}, H::MPOHamiltonian{TensorKit.ComplexSpace, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 2, 2, Matrix{ComplexF64}}, ComplexF64}, alg::GradientGrassmann, envs::MPSKit.MPOHamInfEnv{MPOHamiltonian{TensorKit.ComplexSpace, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 2, 2, Matrix{ComplexF64}}, ComplexF64}, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 2, 1, Matrix{ComplexF64}}, InfiniteMPS{TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 2, 1, Matrix{ComplexF64}}, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 1, 1, Matrix{ComplexF64}}}, KrylovKit.GMRES{KrylovKit.ModifiedGramSchmidt2, Float64}})
   @ MPSKit ~/.julia/packages/MPSKit/atykv/src/algorithms/groundstate/gradient_grassmann.jl:54
 [5] find_groundstate
   @ ~/.julia/packages/MPSKit/atykv/src/algorithms/unionalg.jl:25 [inlined]
 [6] #find_groundstate#380
   @ ~/.julia/packages/MPSKit/atykv/src/algorithms/groundstate/find_groundstate.jl:40 [inlined]
 [7] find_groundstate
   @ ~/.julia/packages/MPSKit/atykv/src/algorithms/groundstate/find_groundstate.jl:19 [inlined]
 [8] find_groundstate(ψ::InfiniteMPS{TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 2, 1, Matrix{ComplexF64}}, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 1, 1, Matrix{ComplexF64}}}, H::MPOHamiltonian{TensorKit.ComplexSpace, TensorKit.TrivialTensorMap{TensorKit.ComplexSpace, 2, 2, Matrix{ComplexF64}}, ComplexF64})
   @ MPSKit ~/.julia/packages/MPSKit/atykv/src/algorithms/groundstate/find_groundstate.jl:19
 [9] top-level scope
   @ ~/Desktop/GitHub/MyGradientFix.jl/test/runtests.jl:12
in expression starting at /Users/daan/Desktop/GitHub/MyGradientFix.jl/test/runtests.jl:12

I made a small package called MyGradientFix.jl that overloads the inner function in GrasmannMPS, hence this stacktrace. The line in question inside OptimKit that throws the error message is

df₀ = inner(x₀, g₀, η₀)
    if df₀ >= zero(df₀)
        error("linesearch was not given a descent direction!")
    end

I do not understand OptimKit well enough to know why this check is required and why it fails.

DaanMaertens avatar Dec 06 '23 14:12 DaanMaertens

optimkit needs to be given a gradient and a descent direction (in gradient descent you would make the descent direction be - gradient). A descent direction has to satsify that inner(x,gradient,direction) is smaller than zero, which is what optimkit tests. Honestly I wouldn't mind if optimkit turns it into a warning ,and just flips the sign.

Gradientgrassman does something wrong in the case that C is numerically degenerate. A real fix would probably entail going through the derivation again, and making sure that everything is defined when C is no longer invertible.

maartenvd avatar Dec 07 '23 09:12 maartenvd