From 613c7baa487497cb4952598c981976e15e6e4e43 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 1 Oct 2026 16:50:11 -0400 Subject: [PATCH 1/3] Guard optimizers and linesearch against non-descent directions and zero steps MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A non-descent search direction (for example from a preconditioner that is indefinite in finite precision) made the linesearch return a zero step, which none of the optimizers handled: - ConjugateGradient: the next Hager-Zhang β is 0/0 = NaN; the resulting NaN direction passed the descent check of the linesearch, whose bracketing phase then loops forever. A zero initial guess (2 * 0) also loops forever there. - LBFGS and GradientDescent: the same zero step is retried until `maxiter`. The optimizers now check for descent before the linesearch. ConjugateGradient restarts from the preconditioned gradient when β is not finite or the conjugate direction is not a descent direction, and LBFGS resets its inverse Hessian approximation. If the preconditioned gradient itself is not a descent direction, or the linesearch makes no progress along it, they stop with a warning. The linesearch rejects invalid initial guesses, treats a NaN slope as non-descent, and bounds its bracket expansion by `maxfg`. Co-Authored-By: Claude Opus 5.5 --- src/cg.jl | 26 ++++++++++++++++++++++++++ src/gd.jl | 14 ++++++++++++++ src/lbfgs.jl | 28 +++++++++++++++++++++++++++- src/linesearches.jl | 15 ++++++++++++++- test/runtests.jl | 40 ++++++++++++++++++++++++++++++++++++++++ 5 files changed, 121 insertions(+), 2 deletions(-) diff --git a/src/cg.jl b/src/cg.jl index 3bab1b8..4118210 100644 --- a/src/cg.jl +++ b/src/cg.jl @@ -125,8 +125,28 @@ function optimize( ) end ) + # e.g. 0/0 after a zero step, where g == gprev + isfinite(β) || (β = zero(α)) η = add!(η, ηprev, β) end + dϕ = inner(x, g, η) + if !(dϕ < 0) && !iszero(β) + verbosity >= 2 && + @info "CG: not a descent direction, restarting with the preconditioned gradient" + β = zero(α) + η = scale!(deepcopy(Pg), -1) + dϕ = inner(x, g, η) + end + if !(dϕ < 0) + verbosity >= 1 && + @warn @sprintf( + "CG: preconditioned gradient is not a descent direction (dϕ = %.2e), stopping", + dϕ + ) + break + end + # after a zero step, doubling would leave the initial guess at zero + iszero(α) && (α = 1 / sqrt(-dϕ)) # store current quantities as previous quantities xprev = x @@ -174,6 +194,12 @@ function optimize( end ηprev = transport!(deepcopy(ηprev), xprev, ηprev, α, x) + if iszero(α) && iszero(β) + verbosity >= 1 && + @warn "CG: linesearch made no progress along the preconditioned gradient, stopping" + break + end + # increase α for next step α = 2 * α end diff --git a/src/gd.jl b/src/gd.jl index 1447d8d..aed0a8d 100644 --- a/src/gd.jl +++ b/src/gd.jl @@ -86,6 +86,15 @@ function optimize( # compute new search direction Pg = precondition(x, deepcopy(g)) η = scale!(Pg, -1) # we don't need g or Pg anymore, so we can overwrite it + dϕ = inner(x, g, η) + if !(dϕ < 0) + verbosity >= 1 && + @warn @sprintf( + "GD: preconditioned gradient is not a descent direction (dϕ = %.2e), stopping", + dϕ + ) + break + end # perform line search _xlast[] = x # store result in global variables to debug linesearch failures @@ -118,6 +127,11 @@ function optimize( numiter, format_time(Δt), f, normgrad, α, nfg ) + if iszero(α) + verbosity >= 1 && @warn "GD: linesearch made no progress, stopping" + break + end + # increase α for next step α = 2 * α end diff --git a/src/lbfgs.jl b/src/lbfgs.jl index 9bbfc81..28d3a57 100644 --- a/src/lbfgs.jl +++ b/src/lbfgs.jl @@ -95,11 +95,26 @@ function optimize( H(g, ξ -> precondition(x, ξ), (ξ1, ξ2) -> inner(x, ξ1, ξ2), add!, scale!) end η = scale!(Hg, -1) - else + if !(inner(x, g, η) < 0) + verbosity >= 2 && + @info "LBFGS: not a descent direction, resetting the inverse Hessian approximation" + empty!(H) + end + end + if length(H) == 0 Pg = precondition(x, deepcopy(g)) normPg = sqrt(inner(x, Pg, Pg)) η = scale!(Pg, -0.01 / normPg) # initial guess: scale invariant end + dϕ = inner(x, g, η) + if !(dϕ < 0) + verbosity >= 1 && + @warn @sprintf( + "LBFGS: preconditioned gradient is not a descent direction (dϕ = %.2e), stopping", + dϕ + ) + break + end # store current quantities as previous quantities xprev = x @@ -139,6 +154,17 @@ function optimize( numiter, format_time(Δt), f, normgrad, α, length(H), nfg ) + if iszero(α) + if length(H) == 0 + verbosity >= 1 && + @warn "LBFGS: linesearch made no progress along the preconditioned gradient, stopping" + break + end + # the same direction would be proposed again + empty!(H) + continue + end + # transport gprev, ηprev and vectors in Hessian approximation to x gprev = transport!(gprev, xprev, ηprev, α, x) for k in 1:length(H) diff --git a/src/linesearches.jl b/src/linesearches.jl index 56a9fb8..eddc052 100644 --- a/src/linesearches.jl +++ b/src/linesearches.jl @@ -124,8 +124,11 @@ function (ls::HagerZhangLineSearch)( ) (f₀, g₀) = fg₀ ϕ₀ = f₀ + if !(isfinite(initialguess) && initialguess > 0) + throw(ArgumentError("initial guess for the step length should be positive and finite, got $initialguess")) + end dϕ₀ = inner(x₀, g₀, η₀) - if dϕ₀ >= zero(dϕ₀) + if !(dϕ₀ < zero(dϕ₀)) @warn "Linesearch was not given a descent direction: returning zero step length" return x₀, f₀, g₀, η₀, zero(one(f₀)), 0 end @@ -425,6 +428,11 @@ function bracket(iter::HagerZhangLineSearchIterator{T}, c::LineSearchPoint) wher c.α, c.dϕ, c.ϕ - p₀.ϕ ) end + if !(isfinite(c.ϕ) && isfinite(c.dϕ)) + verbosity >= 1 && + @warn " Linesearch bracket: no finite function value or slope within the allowed function evaluations" + return a, a, numfg + end c.dϕ >= 0 && return a, c, numfg # B1 # from here: c.dϕ < 0 if c.ϕ > fmax # B2 @@ -432,6 +440,11 @@ function bracket(iter::HagerZhangLineSearchIterator{T}, c::LineSearchPoint) wher return a, b, numfg + nfg else # B3 a = c + if numfg >= iter.parameters.maxfg + verbosity >= 1 && + @warn " Linesearch bracket: slope still negative after the allowed function evaluations" + return a, a, numfg + end α *= iter.parameters.ρ c = takestep(iter, α) numfg += 1 diff --git a/test/runtests.jl b/test/runtests.jl index 837b128..0673628 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -98,6 +98,46 @@ algorithms = (GradientDescent, ConjugateGradient, LBFGS) @test f < 1.0e-12 end +# an indefinite preconditioner eventually produces a non-descent direction +@testset "Non-descent direction $algtype" for algtype in algorithms + fg = quadraticproblem(Matrix(1.0I, 2, 2), zeros(2)) + precondition(x, g) = [g[1], -g[2] / 2] + alg = algtype(; verbosity = 1, gradtol = 1.0e-12, maxiter = 100) + x, f, g, numfg, history = @test_logs (:warn, r"not a descent direction") (:warn, r"not converged") optimize( + fg, [1.0, 0.1], alg; precondition + ) + @test all(isfinite, x) + @test size(history, 1) - 1 < 100 +end + +struct ConstantFlavor <: OptimKit.CGFlavor + β::Float64 +end +(flavor::ConstantFlavor)(args...) = flavor.β + +# a bad β should restart from the preconditioned gradient instead of breaking the optimization +@testset "ConjugateGradient restarts for β = $β" for β in (NaN, -1.0e3) + n = 10 + y = randn(n) + A = randn(n, n) + A = A' * A + I + fg = quadraticproblem(A, y) + alg = ConjugateGradient(; flavor = ConstantFlavor(β), verbosity = 0, gradtol = 1.0e-8, maxiter = 10_000) + x, f, g, numfg, history = optimize(fg, randn(n), alg) + @test all(isfinite, x) + @test x ≈ y rtol = cond(A) * 1.0e-8 +end + +@testset "Linesearch rejects invalid initial guess $α" for α in (0.0, -1.0, NaN, Inf) + fg = x -> (x^2, 2 * x) + @test_throws ArgumentError HagerZhangLineSearch()(fg, 1.0, -2.0; initialguess = α) +end + +@testset "Linesearch with non-finite slope" begin + fg = x -> (x^2, 2 * x) + @test_logs (:warn, r"not given a descent direction") match_mode = :any HagerZhangLineSearch()(fg, 1.0, NaN) +end + include("sphere.jl") @testset "Manifold correctness" begin From de383a1838fc3d00475ade29ed3b723950aa3243 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 1 Oct 2026 17:52:25 -0400 Subject: [PATCH 2/3] Update src/cg.jl Co-authored-by: Jutho --- src/cg.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/cg.jl b/src/cg.jl index 4118210..db0b063 100644 --- a/src/cg.jl +++ b/src/cg.jl @@ -126,8 +126,11 @@ function optimize( end ) # e.g. 0/0 after a zero step, where g == gprev - isfinite(β) || (β = zero(α)) - η = add!(η, ηprev, β) + if isfinite(β) + η = add!(η, ηprev, β) + else + β = zero(α) + end end dϕ = inner(x, g, η) if !(dϕ < 0) && !iszero(β) From fd8870dd6e8cc2b9cf7003ffb2069a489ff3ec81 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Fri, 2 Oct 2026 09:52:43 -0400 Subject: [PATCH 3/3] explain NaN and positive/negative --- src/OptimKit.jl | 8 ++++++++ src/cg.jl | 6 ++++-- src/gd.jl | 3 ++- src/lbfgs.jl | 6 ++++-- src/linesearches.jl | 5 +++-- 5 files changed, 21 insertions(+), 7 deletions(-) diff --git a/src/OptimKit.jl b/src/OptimKit.jl index 426c08a..d2db986 100644 --- a/src/OptimKit.jl +++ b/src/OptimKit.jl @@ -15,6 +15,14 @@ const GRADTOL = ScopedValue(1.0e-8) const MAXITER = ScopedValue(1_000_000) const VERBOSITY = ScopedValue(1) +# `!isnegative(x)` is not equivalent to `x >= 0`: it is also `true` for NaN +@static if !isdefined(Base, :isnegative) # added to Base in Julia 1.13 + isnegative(x::Real) = x < 0 +end +@static if !isdefined(Base, :ispositive) # added to Base in Julia 1.13 + ispositive(x::Real) = x > 0 +end + # Default values for the manifold structure _retract(x, d, α) = (add(x, d, α), d) _invretract(x, y) = add(y, x, -1) diff --git a/src/cg.jl b/src/cg.jl index db0b063..24218a5 100644 --- a/src/cg.jl +++ b/src/cg.jl @@ -133,14 +133,16 @@ function optimize( end end dϕ = inner(x, g, η) - if !(dϕ < 0) && !iszero(β) + # not equivalent to `dϕ >= 0`: a NaN slope must also count as non-descent + if !isnegative(dϕ) && !iszero(β) verbosity >= 2 && @info "CG: not a descent direction, restarting with the preconditioned gradient" β = zero(α) η = scale!(deepcopy(Pg), -1) dϕ = inner(x, g, η) end - if !(dϕ < 0) + # not equivalent to `dϕ >= 0`: a NaN slope must also count as non-descent + if !isnegative(dϕ) verbosity >= 1 && @warn @sprintf( "CG: preconditioned gradient is not a descent direction (dϕ = %.2e), stopping", diff --git a/src/gd.jl b/src/gd.jl index aed0a8d..45bf212 100644 --- a/src/gd.jl +++ b/src/gd.jl @@ -87,7 +87,8 @@ function optimize( Pg = precondition(x, deepcopy(g)) η = scale!(Pg, -1) # we don't need g or Pg anymore, so we can overwrite it dϕ = inner(x, g, η) - if !(dϕ < 0) + # not equivalent to `dϕ >= 0`: a NaN slope must also count as non-descent + if !isnegative(dϕ) verbosity >= 1 && @warn @sprintf( "GD: preconditioned gradient is not a descent direction (dϕ = %.2e), stopping", diff --git a/src/lbfgs.jl b/src/lbfgs.jl index 28d3a57..7a69974 100644 --- a/src/lbfgs.jl +++ b/src/lbfgs.jl @@ -95,7 +95,8 @@ function optimize( H(g, ξ -> precondition(x, ξ), (ξ1, ξ2) -> inner(x, ξ1, ξ2), add!, scale!) end η = scale!(Hg, -1) - if !(inner(x, g, η) < 0) + # not equivalent to `inner(x, g, η) >= 0`: a NaN slope must also count as non-descent + if !isnegative(inner(x, g, η)) verbosity >= 2 && @info "LBFGS: not a descent direction, resetting the inverse Hessian approximation" empty!(H) @@ -107,7 +108,8 @@ function optimize( η = scale!(Pg, -0.01 / normPg) # initial guess: scale invariant end dϕ = inner(x, g, η) - if !(dϕ < 0) + # not equivalent to `dϕ >= 0`: a NaN slope must also count as non-descent + if !isnegative(dϕ) verbosity >= 1 && @warn @sprintf( "LBFGS: preconditioned gradient is not a descent direction (dϕ = %.2e), stopping", diff --git a/src/linesearches.jl b/src/linesearches.jl index eddc052..ee53fc5 100644 --- a/src/linesearches.jl +++ b/src/linesearches.jl @@ -124,11 +124,12 @@ function (ls::HagerZhangLineSearch)( ) (f₀, g₀) = fg₀ ϕ₀ = f₀ - if !(isfinite(initialguess) && initialguess > 0) + if !(isfinite(initialguess) && ispositive(initialguess)) throw(ArgumentError("initial guess for the step length should be positive and finite, got $initialguess")) end dϕ₀ = inner(x₀, g₀, η₀) - if !(dϕ₀ < zero(dϕ₀)) + # not equivalent to `dϕ₀ >= 0`: a NaN slope must also count as non-descent + if !isnegative(dϕ₀) @warn "Linesearch was not given a descent direction: returning zero step length" return x₀, f₀, g₀, η₀, zero(one(f₀)), 0 end