From 26f6c38103b05843229de1927fcf9c8d78c6a453 Mon Sep 17 00:00:00 2001 From: Jutho Haegeman Date: Wed, 25 Mar 2026 22:11:24 +0100 Subject: [PATCH 1/2] add fixedpoint method with 2 algorithms --- src/OptimKit.jl | 68 +++++++++-- src/anderson.jl | 258 +++++++++++++++++++++++++++++++++++++++++ src/simpleiteration.jl | 92 +++++++++++++++ 3 files changed, 411 insertions(+), 7 deletions(-) create mode 100644 src/anderson.jl create mode 100644 src/simpleiteration.jl diff --git a/src/OptimKit.jl b/src/OptimKit.jl index 6415200..f6ab04f 100644 --- a/src/OptimKit.jl +++ b/src/OptimKit.jl @@ -11,12 +11,13 @@ const LS_MAXITER = ScopedValue(10) const LS_MAXFG = ScopedValue(20) const LS_VERBOSITY = ScopedValue(1) -const GRADTOL = ScopedValue(1e-8) +const GRADTOL = ScopedValue(1.0e-8) const MAXITER = ScopedValue(1_000_000) const VERBOSITY = ScopedValue(1) # Default values for the manifold structure _retract(x, d, α) = (add(x, d, α), d) +_invretract(x, y) = add(y, x, -1) _inner(x, v1, v2) = v1 === v2 ? norm(v1)^2 : real(inner(v1, v2)) _transport!(v, xold, d, α, xnew) = v _scale!(v, α) = scale!!(v, α) @@ -26,7 +27,7 @@ _precondition(x, g) = deepcopy(g) _finalize!(x, f, g, numiter) = x, f, g # Default structs for new convergence and termination keywords -@kwdef struct DefaultHasConverged{T<:Real} +@kwdef struct DefaultHasConverged{T <: Real} gradtol::T end @@ -44,6 +45,7 @@ end # Optimization abstract type OptimizationAlgorithm end +abstract type FixedPointAlgorithm end const _xlast = Ref{Any}() const _glast = Ref{Any}() @@ -104,11 +106,61 @@ Also see [`GradientDescent`](@ref), [`ConjugateGradient`](@ref), [`LBFGS`](@ref) """ function optimize end +""" + fixedpoint(fp, x, alg; + finalize! = _finalize!, + shouldstop = DefaultShouldStop(alg.maxiter), + hasconverged = DefaultHasConverged(alg.gradtol), + retract = _retract, invretract = _invretract, + inner = _inner, (transport!) = _transport!, + (scale!) = _scale!, (add!) = _add!, + isometrictransport = (transport! == _transport! && inner == _inner)) + -> x, g, numfp, history + +Find a fixed point of the function `fp` starting from an initial point `x₀` and using the fixed point algorithm `alg`, +which is an instance of `SimpleIteration` or `AndersonMixing`. + +Returns the final point `x`, the residual `g`, the total number of calls to `fp`, +and the history of the norm of the residual across the different iterations. + +The algorithm is run until either `hasconverged(x, false, g, normg)` returns `true` or +`shouldstop(x, false, g, numfp, numiter, time)` returns `true`. The latter case happening before +the former is considered to be a failure to converge, and a warning is issued. + +The keyword arguments are: +- `finalize!::Function`: A function that takes the final point `x`, a dummy variable `f = false`, + the residual `g`, and the iteration number, and returns a possibly modified values for + `x`, `f` and `g`. By default, the identity is used. + It is the user's responsibility to ensure that the modified values do not lead to + inconsistencies within the optimization algorithm. +- `hasconverged::Function`: A function that takes the current point `x`, a dummy variable `f = false`, + the residual `r`, and the norm of the residual, and returns a boolean indicating whether + the fixed point iteration has converged. By default, the norm of the residual is compared to the + tolerance `gradtol` as encoded in the algorithm instance. +- `shouldstop::Function`: A function that takes the current point `x`, a dummy variable `f = false`, + the residuals `g`, the number of calls to `fg`, the iteration number, and the time spent + so far, and returns a boolean indicating whether the optimization should stop. By default, + the number of iterations is compared to the maximum number of iterations as encoded in the + algorithm instance. + +Check the README of this package for further details on creating an algorithm instance, +as well as for the meaning of the remaining keyword arguments and their default values. + +!!! Warning + + The default values of `hasconverged` and `shouldstop` are provided to ensure continuity + with the previous versions of this package. However, this behaviour might change in the + future. + +Also see [`SimpleIteration`](@ref) and [`AndersonMixing`](@ref). +""" +function fixedpoint end + function format_time(t::Float64) - if t < 1e-3 - return @sprintf("%5.1f μs", 1e6*t) + if t < 1.0e-3 + return @sprintf("%5.1f μs", 1.0e6 * t) elseif t < 1 - return @sprintf("%5.1f ms", 1e3*t) + return @sprintf("%5.1f ms", 1.0e3 * t) elseif t < 60 return @sprintf("%5.2f s", t) elseif t < 3600 @@ -122,7 +174,8 @@ include("linesearches.jl") include("gd.jl") include("cg.jl") include("lbfgs.jl") - +include("simpleiteration.jl") +include("anderson.jl") const gd = GradientDescent() const cg = ConjugateGradient() const lbfgs = LBFGS() @@ -131,6 +184,7 @@ export optimize, gd, cg, lbfgs, optimtest export GradientDescent, ConjugateGradient, LBFGS export FletcherReeves, HestenesStiefel, PolakRibiere, HagerZhang, DaiYuan export HagerZhangLineSearch +export fixedpoint, SimpleIteration, AndersonMixing """ optimtest(fg, x, [d]; alpha = -0.1:0.001:0.1, retract = _retract, inner = _inner) @@ -142,7 +196,7 @@ Test the compatibility between the computation of the gradient, the retraction a It is up to the user to check that the values in `dfs1` and `dfs2` match up to expected precision, by inspecting the numerical values or plotting them. If these values don't match, the linesearch in `optimize` cannot be expected to work. """ -function optimtest(fg, x, d=fg(x)[2]; alpha=-0.1:0.001:0.1, retract=_retract, inner=_inner) +function optimtest(fg, x, d = fg(x)[2]; alpha = -0.1:0.001:0.1, retract = _retract, inner = _inner) # evaluate function at given edge points fs_edges = map(alpha) do a f, = fg(retract(x, d, a)[1]) diff --git a/src/anderson.jl b/src/anderson.jl new file mode 100644 index 0000000..1d9327f --- /dev/null +++ b/src/anderson.jl @@ -0,0 +1,258 @@ +""" + struct AndersonMixing{T <: Real} <: FixedPointAlgorithm + AndersonMixing(m::Int; + damping::Real = 1, + maxiter::Int=MAXITER[], # 1_000_000 + gradtol::Real=GRADTOL[], # 1e-8 + verbosity::Int=VERBOSITY[]) # 1 + +Anderson mixing, also known as Anderson acceleration, for fixed point problems. + +## Parameters +- `m::Int`: The number of previous iterates to use for Anderson extrapolation. +- `damping::Real`: The damping parameter for Anderson extrapolation; a value of 1 corresponds to no damping, while a value between 0 and 1 corresponds to under-relaxation. +- `maxiter::Int`: The maximum number of iterations. +- `gradtol::T`: The tolerance for the norm of the residual. +- `verbosity::Int`: The verbosity level of the optimization algorithm. + +The verbosity level use the following scheme: +- 0: no output +- 1: only warnings upon non-convergence +- 2: convergence information at the end of the algorithm +- 3: progress information after each iteration +""" +struct AndersonMixing{T <: Real} <: FixedPointAlgorithm + m::Int + damping::T + maxiter::Int + gradtol::T + verbosity::Int +end +function AndersonMixing(m::Int=8; + damping::Real = 1, + maxiter::Int=MAXITER[], + gradtol::Real=GRADTOL[], + verbosity::Int=VERBOSITY[]) + damping′, gradtol′ = promote(damping, gradtol) + return AndersonMixing(m, damping′, maxiter, gradtol′, verbosity) +end + +using LinearAlgebra +function fixedpoint( + fp::F, x₀, alg::AndersonMixing; + (finalize!) = _finalize!, + shouldstop = DefaultShouldStop(alg.maxiter), + hasconverged = DefaultHasConverged(alg.gradtol), + retract = _retract, invretract = _invretract, + inner = _inner, (transport!) = _transport!, + (scale!) = _scale!, (add!) = _add!, + isometrictransport = (transport! == _transport! && inner == _inner) + ) where {F} + + t₀ = time() + verbosity = alg.verbosity + x = x₀ + fx = fp(x) + numfp = 1 + numiter = 0 + g = invretract(x, fx) # residual, i.e. direction from x to fx, similar to F(x) - x in the vector case + x, _, g = finalize!(x, false, g, numiter) + innergg = inner(x, g, g) + normg = sqrt(innergg) + normghistory = [normg] + verbosity >= 2 && + @info @sprintf("Anderson: initializing with ‖x - F(x)‖ = %.4e", normg) + t = time() - t₀ + _hasconverged = hasconverged(x, false, g, normg) + _shouldstop = shouldstop(x, false, g, numfp, numiter, t) + + if !(_hasconverged || _shouldstop) # first iteration is special + xprev = x + gprev = g + Δxprev = deepcopy(g) + told = t + x, Δx = retract(x, g, 1) + fx = fp(x) + numfp += 1 + numiter += 1 + g = invretract(x, fx) # residual, i.e. direction from x to fx, similar to F(x) - x in the vector case + x, _, g = finalize!(x, false, g, numiter) + innergg = inner(x, g, g) + normg = sqrt(innergg) + push!(normghistory, normg) + t = time() - t₀ + Δt = t - told + _hasconverged = hasconverged(x, false, g, normg) + _shouldstop = shouldstop(x, false, g, numfp, numiter, t) + end + TangentType = typeof(g) + H = AndersonHistory(alg.m, TangentType[], TangentType[]) + overlap_full = zeros(typeof(innergg), alg.m, alg.m) + rhs_full = zeros(typeof(innergg), alg.m) + while !(_hasconverged || _shouldstop) + verbosity >= 3 && + @info @sprintf("Anderson acceleration: iter %4d, Δt %s: ‖f(x) - x‖ = %.4e", numiter, format_time(Δt), normg) + + told = t + # compute new update direction using Anderson extrapolation + gprev = transport!(gprev, xprev, Δxprev, 1, x) + Δg = add!(deepcopy(g), gprev, -1) + # Δx = transport!(deepcopy(Δx), xprev, Δx, 1, x) # transport previous Anderson extrapolation step to current point + + mold = length(H) + push!(H, (Δg, Δx)) + for k in 1:(length(H) - 1) + (Δgₖ, Δxₖ) = H[k] + Δgₖ = transport!(Δgₖ, xprev, Δxprev, 1, x) # transport stored residuals to current point + Δxₖ = transport!(Δxₖ, xprev, Δxprev, 1, x) + H[k] = (Δgₖ, Δxₖ) # update stored residuals and steps to current point + end + m = length(H) + + overlap = view(overlap_full, 1:m, 1:m) + rhs = view(rhs_full, 1:m) + if isometrictransport + if mold == m + @inbounds for j in 1:(m - 1) + for i in 1:(m - 1) + overlap[i, j] = overlap[i + 1, j + 1] + end + end + end + @inbounds for i in 1:(m - 1) + (Δgᵢ, Δxᵢ) = H[i] + overlap[m, i] = overlap[i, m] = inner(x, Δgᵢ, Δg) + end + overlap[m, m] = inner(x, Δg, Δg) + else + @inbounds for j in 1:m + (Δgⱼ, Δxⱼ) = H[j] + for i in 1:(j - 1) + (Δgᵢ, Δxᵢ) = H[i] + overlap[j, i] = overlap[i, j] = inner(x, Δgᵢ, Δgⱼ) + end + overlap[j, j] = inner(x, Δgⱼ, Δgⱼ) + end + end + @inbounds for i in 1:m + (Δgᵢ, Δxᵢ) = H[i] + rhs[i] = inner(x, Δgᵢ, g) + end + # overlapX = [inner(x, H[i][2], H[j][2]) for i in 1:m, j in 1:m] + # overlap += sqrt(eps(one(eltype(overlap)))) * overlapX + overlap += LinearAlgebra.tr(overlap) / size(overlap, 1) * sqrt(eps(one(eltype(overlap)))) * I + Γ = overlap \ rhs # solve least squares problem to get Anderson coefficients + ḡ = deepcopy(g) + Δx = scale!(deepcopy(Δx), 0) + @inbounds for k in 1:m + (Δgₖ, Δxₖ) = H[k] + ḡ = add!(ḡ, Δgₖ, -Γ[k]) + Δx = add!(Δx, Δxₖ, -Γ[k]) # compute Anderson extrapolation direction + end + @show sqrt(inner(x, ḡ, ḡ)) + Δx = add!(Δx, ḡ, alg.damping) # add damping to Anderson extrapolation direction + + # store current quantities as previous quantities + xprev = x + Δxprev = Δx + gprev = g + _xlast[] = x # store result in global variables to debug failures + _glast[] = g + + # take step and evaluate new point + x, Δx = retract(x, Δx, 1) + fx = fp(x) + numfp += 1 + numiter += 1 + g = invretract(x, fx) # Direction from x to fx, i.e. similar to F(x) - x in the vector case + x, _, g = finalize!(x, false, g, numiter) + innergg = inner(x, g, g) + normg = sqrt(innergg) + push!(normghistory, normg) + t = time() - t₀ + Δt = t - told + _hasconverged = hasconverged(x, false, g, normg) + _shouldstop = shouldstop(x, false, g, numfp, numiter, t) + + # check stopping criteria and print info + if _hasconverged || _shouldstop + break + end + end + if _hasconverged + verbosity >= 2 && + @info @sprintf("Anderson acceleration: converged after %d iterations and time %s: ‖f(x) - x‖ = %.4e", numiter, format_time(t), normg) + else + verbosity >= 1 && + @warn @sprintf("Anderson acceleration: not converged to requested tol after %d iterations and time %s: ‖f(x) - x‖ = %.4e", numiter, format_time(t), normg) + end + return x, g, numfp, normghistory +end + +mutable struct AndersonHistory{TangentType} + maxlength::Int + length::Int + first::Int + Δgesiduals::Vector{TangentType} + Δpositions::Vector{TangentType} + function AndersonHistory{T}(maxlength::Int, Δgesiduals::Vector{T}, Δpositions::Vector{T}) where {T} + l = length(Δgesiduals) + @assert l == length(Δpositions) "AndersonHistory: Δgesiduals and Δpositions must have the same length" + @assert l <= maxlength "AndersonHistory: initial history length cannot exceed maxlength" + Δgesiduals = resize!(copy(Δgesiduals), maxlength) + Δpositions = resize!(copy(Δpositions), maxlength) + return new{T}(maxlength, l, 1, Δgesiduals, Δpositions) + end +end +function AndersonHistory(maxlength::Int, Δgesiduals::Vector{T}, Δpositions::Vector{T}) where {T} + return AndersonHistory{T}(maxlength, Δgesiduals, Δpositions) +end + +Base.length(H::AndersonHistory) = H.length + +@inline function Base.getindex(H::AndersonHistory, i::Int) + @boundscheck if i < 1 || i > H.length + throw(BoundsError(H, i)) + end + n = H.maxlength + idx = H.first + i - 1 + idx = ifelse(idx > n, idx - n, idx) + return (getindex(H.Δgesiduals, idx), getindex(H.Δpositions, idx)) +end + +@inline function Base.setindex!(H::AndersonHistory, (Δg, Δx), i) + @boundscheck if i < 1 || i > H.length + throw(BoundsError(H, i)) + end + idx = mod1(H.first + i - 1, H.maxlength) + setindex!(H.Δgesiduals, Δg, idx) + setindex!(H.Δpositions, Δx, idx) + return (Δg, Δx) +end + +@inline function Base.push!(H::AndersonHistory, value) + if H.length < H.maxlength + H.length += 1 + else + H.first = mod1(H.first + 1, H.maxlength) + end + @inbounds setindex!(H, value, H.length) + return H +end +@inline function Base.pop!(H::AndersonHistory) + @inbounds v = H[H.length] + H.length -= 1 + return v +end +@inline function Base.popfirst!(H::AndersonHistory) + @inbounds v = H[1] + H.first = mod1(H.first + 1, H.maxlength) + H.length -= 1 + return v +end + +@inline function Base.empty!(H::AndersonHistory) + H.length = 0 + H.first = 1 + return H +end diff --git a/src/simpleiteration.jl b/src/simpleiteration.jl new file mode 100644 index 0000000..8c4eccd --- /dev/null +++ b/src/simpleiteration.jl @@ -0,0 +1,92 @@ +""" + struct SimpleIteration{T <: Real} <: FixedPointAlgorithm + SimpleIteration(; + maxiter::Int=MAXITER[], # 1_000_000 + gradtol::Real=GRADTOL[], # 1e-8 + verbosity::Int=VERBOSITY[]) # 1 + +Simple iteration for fixed point problems. + +## Parameters +- `maxiter::Int`: The maximum number of iterations. +- `gradtol::Real`: The tolerance for the norm of the residual. +- `verbosity::Int`: The verbosity level of the optimization algorithm. + +The verbosity level use the following scheme: +- 0: no output +- 1: only warnings upon non-convergence +- 2: convergence information at the end of the algorithm +- 3: progress information after each iteration +""" +struct SimpleIteration{T <: Real} <: FixedPointAlgorithm + maxiter::Int + gradtol::T + verbosity::Int +end +function SimpleIteration(; + maxiter::Int=MAXITER[], + gradtol::Real=GRADTOL[], + verbosity::Int=VERBOSITY[]) + return SimpleIteration(maxiter, gradtol, verbosity) +end + + +function fixedpoint( + fp, x₀, alg::SimpleIteration; + (finalize!) = _finalize!, + shouldstop = DefaultShouldStop(alg.maxiter), + hasconverged = DefaultHasConverged(alg.gradtol), + retract = _retract, invretract = _invretract, + inner = _inner, (transport!) = _transport!, + (scale!) = _scale!, (add!) = _add!, + isometrictransport = (transport! == _transport! && inner == _inner) + ) + + t₀ = time() + verbosity = alg.verbosity + x = x₀ + fx = fp(x) + numfp = 1 + numiter = 0 + g = invretract(x, fx) # residual, i.e. direction from x to fx, similar to F(x) - x in the vector case + x, _, g = finalize!(x, false, g, numiter) + innergg = inner(x, g, g) + normg = sqrt(innergg) + normghistory = [normg] + verbosity >= 2 && + @info @sprintf("SimpleIteration: initializing with ‖x - F(x)‖ = %.4e", normg) + t = time() - t₀ + _hasconverged = hasconverged(x, false, g, normg) + _shouldstop = shouldstop(x, false, g, numfp, numiter, t) + + while !(_hasconverged || _shouldstop) + told = t + x = fx + fx = fp(x) + numfp += 1 + numiter += 1 + g = invretract(x, fx) # residual, i.e. direction from x to fx, similar to F(x) - x in the vector case + x, _, g = finalize!(x, false, g, numiter) + innergg = inner(x, g, g) + normg = sqrt(innergg) + push!(normghistory, normg) + t = time() - t₀ + Δt = t - told + _hasconverged = hasconverged(x, false, g, normg) + _shouldstop = shouldstop(x, false, g, numfp, numiter, t) + # check stopping criteria and print info + if _hasconverged || _shouldstop + break + end + verbosity >= 3 && + @info @sprintf("SimpleIteration: iter %4d, Δt %s: ‖f(x) - x‖ = %.4e", numiter, format_time(Δt), normg) + end + if _hasconverged + verbosity >= 2 && + @info @sprintf("SimpleIteration: converged after %d iterations and time %s: ‖f(x) - x‖ = %.4e", numiter, format_time(t), normg) + else + verbosity >= 1 && + @warn @sprintf("SimpleIteration: not converged to requested tol after %d iterations and time %s: ‖f(x) - x‖ = %.4e", numiter, format_time(t), normg) + end + return x, g, numfp, normghistory +end From 70df4a6e423f1da65572dd35d33a6f6bdbe45448 Mon Sep 17 00:00:00 2001 From: Jutho Haegeman Date: Fri, 3 Jul 2026 00:46:23 +0200 Subject: [PATCH 2/2] add correctness test and formatting --- src/anderson.jl | 39 +++++----- src/cg.jl | 90 +++++++++++++---------- src/gd.jl | 68 ++++++++++-------- src/lbfgs.jl | 97 ++++++++++++++----------- src/linesearches.jl | 159 ++++++++++++++++++++++++++--------------- src/linesearches2.jl | 40 +++++++++++ src/simpleiteration.jl | 7 +- test/runtests.jl | 73 ++++++++++++++----- test/sphere.jl | 42 +++++++++++ 9 files changed, 412 insertions(+), 203 deletions(-) create mode 100644 src/linesearches2.jl create mode 100644 test/sphere.jl diff --git a/src/anderson.jl b/src/anderson.jl index 1d9327f..03503c5 100644 --- a/src/anderson.jl +++ b/src/anderson.jl @@ -8,6 +8,8 @@ Anderson mixing, also known as Anderson acceleration, for fixed point problems. +WARNING: Experimental implementation – subject to change. + ## Parameters - `m::Int`: The number of previous iterates to use for Anderson extrapolation. - `damping::Real`: The damping parameter for Anderson extrapolation; a value of 1 corresponds to no damping, while a value between 0 and 1 corresponds to under-relaxation. @@ -28,11 +30,13 @@ struct AndersonMixing{T <: Real} <: FixedPointAlgorithm gradtol::T verbosity::Int end -function AndersonMixing(m::Int=8; - damping::Real = 1, - maxiter::Int=MAXITER[], - gradtol::Real=GRADTOL[], - verbosity::Int=VERBOSITY[]) +function AndersonMixing( + m::Int = 8; + damping::Real = 1, + maxiter::Int = MAXITER[], + gradtol::Real = GRADTOL[], + verbosity::Int = VERBOSITY[] + ) damping′, gradtol′ = promote(damping, gradtol) return AndersonMixing(m, damping′, maxiter, gradtol′, verbosity) end @@ -100,7 +104,7 @@ function fixedpoint( # Δx = transport!(deepcopy(Δx), xprev, Δx, 1, x) # transport previous Anderson extrapolation step to current point mold = length(H) - push!(H, (Δg, Δx)) + push!(H, (Δg, Δx)) for k in 1:(length(H) - 1) (Δgₖ, Δxₖ) = H[k] Δgₖ = transport!(Δgₖ, xprev, Δxprev, 1, x) # transport stored residuals to current point @@ -141,7 +145,7 @@ function fixedpoint( # overlapX = [inner(x, H[i][2], H[j][2]) for i in 1:m, j in 1:m] # overlap += sqrt(eps(one(eltype(overlap)))) * overlapX overlap += LinearAlgebra.tr(overlap) / size(overlap, 1) * sqrt(eps(one(eltype(overlap)))) * I - Γ = overlap \ rhs # solve least squares problem to get Anderson coefficients + Γ = LinearAlgebra.cholesky(overlap) \ rhs # solve least squares problem to get Anderson coefficients ḡ = deepcopy(g) Δx = scale!(deepcopy(Δx), 0) @inbounds for k in 1:m @@ -149,7 +153,6 @@ function fixedpoint( ḡ = add!(ḡ, Δgₖ, -Γ[k]) Δx = add!(Δx, Δxₖ, -Γ[k]) # compute Anderson extrapolation direction end - @show sqrt(inner(x, ḡ, ḡ)) Δx = add!(Δx, ḡ, alg.damping) # add damping to Anderson extrapolation direction # store current quantities as previous quantities @@ -193,19 +196,19 @@ mutable struct AndersonHistory{TangentType} maxlength::Int length::Int first::Int - Δgesiduals::Vector{TangentType} + Δresiduals::Vector{TangentType} Δpositions::Vector{TangentType} - function AndersonHistory{T}(maxlength::Int, Δgesiduals::Vector{T}, Δpositions::Vector{T}) where {T} - l = length(Δgesiduals) - @assert l == length(Δpositions) "AndersonHistory: Δgesiduals and Δpositions must have the same length" + function AndersonHistory{T}(maxlength::Int, Δresiduals::Vector{T}, Δpositions::Vector{T}) where {T} + l = length(Δresiduals) + @assert l == length(Δpositions) "AndersonHistory: Δresiduals and Δpositions must have the same length" @assert l <= maxlength "AndersonHistory: initial history length cannot exceed maxlength" - Δgesiduals = resize!(copy(Δgesiduals), maxlength) + Δresiduals = resize!(copy(Δresiduals), maxlength) Δpositions = resize!(copy(Δpositions), maxlength) - return new{T}(maxlength, l, 1, Δgesiduals, Δpositions) + return new{T}(maxlength, l, 1, Δresiduals, Δpositions) end end -function AndersonHistory(maxlength::Int, Δgesiduals::Vector{T}, Δpositions::Vector{T}) where {T} - return AndersonHistory{T}(maxlength, Δgesiduals, Δpositions) +function AndersonHistory(maxlength::Int, Δresiduals::Vector{T}, Δpositions::Vector{T}) where {T} + return AndersonHistory{T}(maxlength, Δresiduals, Δpositions) end Base.length(H::AndersonHistory) = H.length @@ -217,7 +220,7 @@ Base.length(H::AndersonHistory) = H.length n = H.maxlength idx = H.first + i - 1 idx = ifelse(idx > n, idx - n, idx) - return (getindex(H.Δgesiduals, idx), getindex(H.Δpositions, idx)) + return (getindex(H.Δresiduals, idx), getindex(H.Δpositions, idx)) end @inline function Base.setindex!(H::AndersonHistory, (Δg, Δx), i) @@ -225,7 +228,7 @@ end throw(BoundsError(H, i)) end idx = mod1(H.first + i - 1, H.maxlength) - setindex!(H.Δgesiduals, Δg, idx) + setindex!(H.Δresiduals, Δg, idx) setindex!(H.Δpositions, Δx, idx) return (Δg, Δx) end diff --git a/src/cg.jl b/src/cg.jl index 97426ea..3bab1b8 100644 --- a/src/cg.jl +++ b/src/cg.jl @@ -41,7 +41,7 @@ The `flavor` parameter can take the values - `PolakRibiere(; pos = true)`: Polak-Ribiere formula for β - `DaiYuan()`: Dai-Yuan formula for β """ -struct ConjugateGradient{F<:CGFlavor,T<:Real,L<:AbstractLineSearch} <: OptimizationAlgorithm +struct ConjugateGradient{F <: CGFlavor, T <: Real, L <: AbstractLineSearch} <: OptimizationAlgorithm flavor::F restart::Int maxiter::Int @@ -50,29 +50,33 @@ struct ConjugateGradient{F<:CGFlavor,T<:Real,L<:AbstractLineSearch} <: Optimizat linesearch::L end function ConjugateGradient(; - flavor::CGFlavor=HagerZhang(), - restart::Int=typemax(Int), - maxiter::Int=MAXITER[], - gradtol::Real=GRADTOL[], - verbosity::Int=VERBOSITY[], - ls_maxiter::Int=LS_MAXITER[], - ls_maxfg::Int=LS_MAXFG[], - ls_verbosity::Int=LS_VERBOSITY[], - linesearch::AbstractLineSearch=HagerZhangLineSearch(; - maxiter=ls_maxiter, - maxfg=ls_maxfg, - verbosity=ls_verbosity)) + flavor::CGFlavor = HagerZhang(), + restart::Int = typemax(Int), + maxiter::Int = MAXITER[], + gradtol::Real = GRADTOL[], + verbosity::Int = VERBOSITY[], + ls_maxiter::Int = LS_MAXITER[], + ls_maxfg::Int = LS_MAXFG[], + ls_verbosity::Int = LS_VERBOSITY[], + linesearch::AbstractLineSearch = HagerZhangLineSearch(; + maxiter = ls_maxiter, + maxfg = ls_maxfg, + verbosity = ls_verbosity + ) + ) return ConjugateGradient(flavor, restart, maxiter, gradtol, verbosity, linesearch) end -function optimize(fg, x, alg::ConjugateGradient; - precondition=_precondition, - (finalize!)=_finalize!, - shouldstop=DefaultShouldStop(alg.maxiter), - hasconverged=DefaultHasConverged(alg.gradtol), - retract=_retract, inner=_inner, (transport!)=_transport!, - (scale!)=_scale!, (add!)=_add!, - isometrictransport=(transport! == _transport! && inner == _inner)) +function optimize( + fg, x, alg::ConjugateGradient; + precondition = _precondition, + (finalize!) = _finalize!, + shouldstop = DefaultShouldStop(alg.maxiter), + hasconverged = DefaultHasConverged(alg.gradtol), + retract = _retract, inner = _inner, (transport!) = _transport!, + (scale!) = _scale!, (add!) = _add!, + isometrictransport = (transport! == _transport! && inner == _inner) + ) t₀ = time() verbosity = alg.verbosity f, g = fg(x) @@ -112,11 +116,15 @@ function optimize(fg, x, alg::ConjugateGradient; if mod(numiter, alg.restart) == 0 β = zero(α) else - β = oftype(α, - let x = x - alg.flavor(g, gprev, Pg, Pgprev, ηprev, - (η₁, η₂) -> inner(x, η₁, η₂)) - end) + β = oftype( + α, + let x = x + alg.flavor( + g, gprev, Pg, Pgprev, ηprev, + (η₁, η₂) -> inner(x, η₁, η₂) + ) + end + ) η = add!(η, ηprev, β) end @@ -130,9 +138,11 @@ function optimize(fg, x, alg::ConjugateGradient; _xlast[] = x # store result in global variables to debug linesearch failures _glast[] = g _dlast[] = η - x, f, g, ξ, α, nfg = alg.linesearch(fg, x, η, (f, g); - initialguess=α, - retract=retract, inner=inner) + x, f, g, ξ, α, nfg = alg.linesearch( + fg, x, η, (f, g); + initialguess = α, + retract = retract, inner = inner + ) numfg += nfg numiter += 1 x, f, g = finalize!(x, f, g, numiter) @@ -150,8 +160,10 @@ function optimize(fg, x, alg::ConjugateGradient; break end verbosity >= 3 && - @info @sprintf("CG: iter %4d, Δt %s: f = %.12e, ‖∇f‖ = %.4e, α = %.2e, β = %.2e, nfg = %d", - numiter, format_time(Δt), f, normgrad, α, β, nfg) + @info @sprintf( + "CG: iter %4d, Δt %s: f = %.12e, ‖∇f‖ = %.4e, α = %.2e, β = %.2e, nfg = %d", + numiter, format_time(Δt), f, normgrad, α, β, nfg + ) # transport gprev, ηprev and vectors in Hessian approximation to x gprev = transport!(gprev, xprev, ηprev, α, x) @@ -167,22 +179,26 @@ function optimize(fg, x, alg::ConjugateGradient; end if _hasconverged verbosity >= 2 && - @info @sprintf("CG: converged after %d iterations and time %s: f = %.12e, ‖∇f‖ = %.4e", - numiter, format_time(t), f, normgrad) + @info @sprintf( + "CG: converged after %d iterations and time %s: f = %.12e, ‖∇f‖ = %.4e", + numiter, format_time(t), f, normgrad + ) else verbosity >= 1 && - @warn @sprintf("CG: not converged to requested tol after %d iterations and time %s: f = %.12e, ‖∇f‖ = %.4e", - numiter, format_time(t), f, normgrad) + @warn @sprintf( + "CG: not converged to requested tol after %d iterations and time %s: f = %.12e, ‖∇f‖ = %.4e", + numiter, format_time(t), f, normgrad + ) end history = [fhistory normgradhistory] return x, f, g, numfg, history end -struct HagerZhang{T<:Real} <: CGFlavor +struct HagerZhang{T <: Real} <: CGFlavor η::T θ::T end -HagerZhang(; η::Real=4 // 10, θ::Real=1 // 1) = HagerZhang(promote(η, θ)...) +HagerZhang(; η::Real = 4 // 10, θ::Real = 1 // 1) = HagerZhang(promote(η, θ)...) function (HZ::HagerZhang)(g, gprev, Pg, Pgprev, dprev, inner) dd = inner(dprev, dprev) diff --git a/src/gd.jl b/src/gd.jl index 6b9c3ac..1447d8d 100644 --- a/src/gd.jl +++ b/src/gd.jl @@ -28,34 +28,38 @@ Both `verbosity` and `ls_verbosity` use the following scheme: - 3: progress information after each iteration - 4: more detailed information (only for the linesearch) """ -struct GradientDescent{T<:Real,L<:AbstractLineSearch} <: OptimizationAlgorithm +struct GradientDescent{T <: Real, L <: AbstractLineSearch} <: OptimizationAlgorithm maxiter::Int gradtol::T verbosity::Int linesearch::L end function GradientDescent(; - maxiter::Int=MAXITER[], - gradtol::Real=GRADTOL[], - verbosity::Int=VERBOSITY[], - ls_maxiter::Int=LS_MAXITER[], - ls_maxfg::Int=LS_MAXFG[], - ls_verbosity::Int=LS_VERBOSITY[], - linesearch::AbstractLineSearch=HagerZhangLineSearch(; - maxiter=ls_maxiter, - maxfg=ls_maxfg, - verbosity=ls_verbosity)) + maxiter::Int = MAXITER[], + gradtol::Real = GRADTOL[], + verbosity::Int = VERBOSITY[], + ls_maxiter::Int = LS_MAXITER[], + ls_maxfg::Int = LS_MAXFG[], + ls_verbosity::Int = LS_VERBOSITY[], + linesearch::AbstractLineSearch = HagerZhangLineSearch(; + maxiter = ls_maxiter, + maxfg = ls_maxfg, + verbosity = ls_verbosity + ) + ) return GradientDescent(maxiter, gradtol, verbosity, linesearch) end -function optimize(fg, x, alg::GradientDescent; - precondition=_precondition, - (finalize!)=_finalize!, - shouldstop=DefaultShouldStop(alg.maxiter), - hasconverged=DefaultHasConverged(alg.gradtol), - retract=_retract, inner=_inner, (transport!)=_transport!, - (scale!)=_scale!, (add!)=_add!, - isometrictransport=(transport! == _transport! && inner == _inner)) +function optimize( + fg, x, alg::GradientDescent; + precondition = _precondition, + (finalize!) = _finalize!, + shouldstop = DefaultShouldStop(alg.maxiter), + hasconverged = DefaultHasConverged(alg.gradtol), + retract = _retract, inner = _inner, (transport!) = _transport!, + (scale!) = _scale!, (add!) = _add!, + isometrictransport = (transport! == _transport! && inner == _inner) + ) t₀ = time() verbosity = alg.verbosity f, g = fg(x) @@ -87,9 +91,11 @@ function optimize(fg, x, alg::GradientDescent; _xlast[] = x # store result in global variables to debug linesearch failures _glast[] = g _dlast[] = η - x, f, g, ξ, α, nfg = alg.linesearch(fg, x, η, (f, g); - initialguess=α, - retract=retract, inner=inner) + x, f, g, ξ, α, nfg = alg.linesearch( + fg, x, η, (f, g); + initialguess = α, + retract = retract, inner = inner + ) numfg += nfg numiter += 1 x, f, g = finalize!(x, f, g, numiter) @@ -107,20 +113,26 @@ function optimize(fg, x, alg::GradientDescent; break end verbosity >= 3 && - @info @sprintf("GD: iter %4d, Δt %s: f = %.12e, ‖∇f‖ = %.4e, α = %.2e, nfg = %d", - numiter, format_time(Δt), f, normgrad, α, nfg) + @info @sprintf( + "GD: iter %4d, Δt %s: f = %.12e, ‖∇f‖ = %.4e, α = %.2e, nfg = %d", + numiter, format_time(Δt), f, normgrad, α, nfg + ) # increase α for next step α = 2 * α end if _hasconverged verbosity >= 2 && - @info @sprintf("GD: converged after %d iterations and time %s: f = %.12e, ‖∇f‖ = %.4e", - numiter, format_time(t), f, normgrad) + @info @sprintf( + "GD: converged after %d iterations and time %s: f = %.12e, ‖∇f‖ = %.4e", + numiter, format_time(t), f, normgrad + ) else verbosity >= 1 && - @warn @sprintf("GD: not converged to requested tol after %d iterations and time %s: f = %.12e, ‖∇f‖ = %.4e", - numiter, format_time(t), f, normgrad) + @warn @sprintf( + "GD: not converged to requested tol after %d iterations and time %s: f = %.12e, ‖∇f‖ = %.4e", + numiter, format_time(t), f, normgrad + ) end history = [fhistory normgradhistory] return x, f, g, numfg, history diff --git a/src/lbfgs.jl b/src/lbfgs.jl index 68bc51d..9bbfc81 100644 --- a/src/lbfgs.jl +++ b/src/lbfgs.jl @@ -30,7 +30,7 @@ Both `verbosity` and `ls_verbosity` use the following scheme: - 3: progress information after each iteration - 4: more detailed information (only for the linesearch) """ -struct LBFGS{T<:Real,L<:AbstractLineSearch} <: OptimizationAlgorithm +struct LBFGS{T <: Real, L <: AbstractLineSearch} <: OptimizationAlgorithm m::Int maxiter::Int gradtol::T @@ -38,29 +38,34 @@ struct LBFGS{T<:Real,L<:AbstractLineSearch} <: OptimizationAlgorithm verbosity::Int linesearch::L end -function LBFGS(m::Int=8; - acceptfirst::Bool=true, - maxiter::Int=MAXITER[], - gradtol::Real=GRADTOL[], - verbosity::Int=VERBOSITY[], - ls_maxiter::Int=LS_MAXITER[], - ls_maxfg::Int=LS_MAXFG[], - ls_verbosity::Int=LS_VERBOSITY[], - linesearch::AbstractLineSearch=HagerZhangLineSearch(; - maxiter=ls_maxiter, - maxfg=ls_maxfg, - verbosity=ls_verbosity)) +function LBFGS( + m::Int = 8; + acceptfirst::Bool = true, + maxiter::Int = MAXITER[], + gradtol::Real = GRADTOL[], + verbosity::Int = VERBOSITY[], + ls_maxiter::Int = LS_MAXITER[], + ls_maxfg::Int = LS_MAXFG[], + ls_verbosity::Int = LS_VERBOSITY[], + linesearch::AbstractLineSearch = HagerZhangLineSearch(; + maxiter = ls_maxiter, + maxfg = ls_maxfg, + verbosity = ls_verbosity + ) + ) return LBFGS(m, maxiter, gradtol, acceptfirst, verbosity, linesearch) end -function optimize(fg, x, alg::LBFGS; - precondition=_precondition, - (finalize!)=_finalize!, - shouldstop=DefaultShouldStop(alg.maxiter), - hasconverged=DefaultHasConverged(alg.gradtol), - retract=_retract, inner=_inner, (transport!)=_transport!, - (scale!)=_scale!, (add!)=_add!, - isometrictransport=(transport! == _transport! && inner == _inner)) +function optimize( + fg, x, alg::LBFGS; + precondition = _precondition, + (finalize!) = _finalize!, + shouldstop = DefaultShouldStop(alg.maxiter), + hasconverged = DefaultHasConverged(alg.gradtol), + retract = _retract, inner = _inner, (transport!) = _transport!, + (scale!) = _scale!, (add!) = _add!, + isometrictransport = (transport! == _transport! && inner == _inner) + ) t₀ = time() verbosity = alg.verbosity f, g = fg(x) @@ -105,11 +110,13 @@ function optimize(fg, x, alg::LBFGS; _xlast[] = x # store result in global variables to debug linesearch failures _glast[] = g _dlast[] = η - x, f, g, ξ, α, nfg = alg.linesearch(fg, x, η, (f, g); - initialguess=one(f), - acceptfirst=alg.acceptfirst, - # for some reason, line search seems to converge to solution alpha = 2 in most cases if acceptfirst = false. If acceptfirst = true, the initial value of alpha can immediately be accepted. This typically leads to a more erratic convergence of normgrad, but to less function evaluations in the end. - retract=retract, inner=inner) + x, f, g, ξ, α, nfg = alg.linesearch( + fg, x, η, (f, g); + initialguess = one(f), + acceptfirst = alg.acceptfirst, + # for some reason, line search seems to converge to solution alpha = 2 in most cases if acceptfirst = false. If acceptfirst = true, the initial value of alpha can immediately be accepted. This typically leads to a more erratic convergence of normgrad, but to less function evaluations in the end. + retract = retract, inner = inner + ) numfg += nfg numiter += 1 x, f, g = finalize!(x, f, g, numiter) @@ -127,8 +134,10 @@ function optimize(fg, x, alg::LBFGS; break end verbosity >= 3 && - @info @sprintf("LBFGS: iter %4d, Δt %s: f = %.12e, ‖∇f‖ = %.4e, α = %.2e, m = %d, nfg = %d", - numiter, format_time(Δt), f, normgrad, α, length(H), nfg) + @info @sprintf( + "LBFGS: iter %4d, Δt %s: f = %.12e, ‖∇f‖ = %.4e, α = %.2e, m = %d, nfg = %d", + numiter, format_time(Δt), f, normgrad, α, length(H), nfg + ) # transport gprev, ηprev and vectors in Hessian approximation to x gprev = transport!(gprev, xprev, ηprev, α, x) @@ -192,18 +201,22 @@ function optimize(fg, x, alg::LBFGS; end if _hasconverged verbosity >= 2 && - @info @sprintf("LBFGS: converged after %d iterations and time %s: f = %.12e, ‖∇f‖ = %.4e", - numiter, format_time(t), f, normgrad) + @info @sprintf( + "LBFGS: converged after %d iterations and time %s: f = %.12e, ‖∇f‖ = %.4e", + numiter, format_time(t), f, normgrad + ) else verbosity >= 1 && - @warn @sprintf("LBFGS: not converged to requested tol after %d iterations and time %s: f = %.12e, ‖∇f‖ = %.4e", - numiter, format_time(t), f, normgrad) + @warn @sprintf( + "LBFGS: not converged to requested tol after %d iterations and time %s: f = %.12e, ‖∇f‖ = %.4e", + numiter, format_time(t), f, normgrad + ) end history = [fhistory normgradhistory] return x, f, g, numfg, history end -mutable struct LBFGSInverseHessian{TangentType,ScalarType} +mutable struct LBFGSInverseHessian{TangentType, ScalarType} maxlength::Int length::Int first::Int @@ -211,20 +224,24 @@ mutable struct LBFGSInverseHessian{TangentType,ScalarType} Y::Vector{TangentType} ρ::Vector{ScalarType} α::Vector{ScalarType} # work space - function LBFGSInverseHessian{T1,T2}(maxlength::Int, S::Vector{T1}, Y::Vector{T1}, - ρ::Vector{T2}) where {T1,T2} + function LBFGSInverseHessian{T1, T2}( + maxlength::Int, S::Vector{T1}, Y::Vector{T1}, + ρ::Vector{T2} + ) where {T1, T2} @assert length(S) == length(Y) == length(ρ) l = length(S) S = resize!(copy(S), maxlength) Y = resize!(copy(Y), maxlength) ρ = resize!(copy(ρ), maxlength) α = similar(ρ) - return new{T1,T2}(maxlength, l, 1, S, Y, ρ, α) + return new{T1, T2}(maxlength, l, 1, S, Y, ρ, α) end end -function LBFGSInverseHessian(maxlength::Int, S::Vector{T1}, Y::Vector{T1}, - ρ::Vector{T2}) where {T1,T2} - return LBFGSInverseHessian{T1,T2}(maxlength, S, Y, ρ) +function LBFGSInverseHessian( + maxlength::Int, S::Vector{T1}, Y::Vector{T1}, + ρ::Vector{T2} + ) where {T1, T2} + return LBFGSInverseHessian{T1, T2}(maxlength, S, Y, ρ) end Base.length(H::LBFGSInverseHessian) = H.length @@ -276,7 +293,7 @@ end return H end -function (H::LBFGSInverseHessian)(g, precondition, inner, add!, scale!; α=H.α) +function (H::LBFGSInverseHessian)(g, precondition, inner, add!, scale!; α = H.α) q = deepcopy(g) for k in length(H):-1:1 s, y, ρ = H[k] diff --git a/src/linesearches.jl b/src/linesearches.jl index ba7fac9..56a9fb8 100644 --- a/src/linesearches.jl +++ b/src/linesearches.jl @@ -2,7 +2,7 @@ # (Hager & Zhang, ACM Transactions on Mathematical Software, Vol 32 (2006)) abstract type AbstractLineSearch end -struct LineSearchPoint{T<:Real,X,G} +struct LineSearchPoint{T <: Real, X, G} α::T # step length ϕ::T # local function value dϕ::T # local directional derivative of cost function: dϕ/dα @@ -23,7 +23,7 @@ end # Algorithm 851: CG_DESCENT # (Hager & Zhang, ACM Transactions on Mathematical Software, Vol 32 (2006)) -struct HagerZhangLineSearch{T<:Real} <: AbstractLineSearch +struct HagerZhangLineSearch{T <: Real} <: AbstractLineSearch c₁::T # parameter for the (approximate) first Wolfe condition (Armijo rule), controlling sufficient decay of function value: c₁ < 1/2 < c₂ c₂::T # parameter for the second Wolfe condition (curvature contition), controlling sufficient decay of slope: c₁ < 1/2 < c₂ ϵ::T # parameter for expected accuracy of objective function, controlling maximal allowed increase in function value @@ -58,15 +58,17 @@ This method returns a `HagerZhangLineSearch` object `ls`, that can then be can b to perform a line search for function `fg` that computes the objective function and its gradient at a given point, starting from `x₀` in direction `η₀`. """ -function HagerZhangLineSearch(; c₁::Real=1 // 10, - c₂::Real=9 // 10, - ϵ::Real=1 // 10^6, - θ::Real=1 // 2, - γ::Real=2 // 3, - ρ::Real=5 // 1, - maxiter::Int=LS_MAXITER[], - maxfg::Int=LS_MAXFG[], - verbosity::Int=LS_VERBOSITY[]) +function HagerZhangLineSearch(; + c₁::Real = 1 // 10, + c₂::Real = 9 // 10, + ϵ::Real = 1 // 10^6, + θ::Real = 1 // 2, + γ::Real = 2 // 3, + ρ::Real = 5 // 1, + maxiter::Int = LS_MAXITER[], + maxfg::Int = LS_MAXFG[], + verbosity::Int = LS_VERBOSITY[] + ) return HagerZhangLineSearch(promote(c₁, c₂, ϵ, θ, γ, ρ)..., maxiter, maxfg, verbosity) end @@ -111,13 +113,15 @@ Perform a Hager-Zhang line search to find a step length that satisfies the (appr - `α`: Step length that satisfies the (approximate) Wolfe conditions. - `numfg`: Number of function evaluations performed. """ -function (ls::HagerZhangLineSearch)(fg, x₀, η₀, fg₀=fg(x₀); - retract=_retract, inner=_inner, - initialguess::Real=one(fg₀[1]), - acceptfirst::Bool=false, - maxiter::Int=ls.maxiter, - maxfg::Int=ls.maxfg, - verbosity::Int=ls.verbosity) +function (ls::HagerZhangLineSearch)( + fg, x₀, η₀, fg₀ = fg(x₀); + retract = _retract, inner = _inner, + initialguess::Real = one(fg₀[1]), + acceptfirst::Bool = false, + maxiter::Int = ls.maxiter, + maxfg::Int = ls.maxfg, + verbosity::Int = ls.verbosity + ) (f₀, g₀) = fg₀ ϕ₀ = f₀ dϕ₀ = inner(x₀, g₀, η₀) @@ -127,8 +131,10 @@ function (ls::HagerZhangLineSearch)(fg, x₀, η₀, fg₀=fg(x₀); end p₀ = LineSearchPoint(zero(ϕ₀), ϕ₀, dϕ₀, x₀, f₀, g₀, η₀) - iter = HagerZhangLineSearchIterator(fg, retract, inner, p₀, η₀, initialguess, - acceptfirst, verbosity, ls) + iter = HagerZhangLineSearchIterator( + fg, retract, inner, p₀, η₀, initialguess, + acceptfirst, verbosity, ls + ) verbosity >= 2 && @info @sprintf("Linesearch start: dϕ₀ = %.2e, ϕ₀ = %.2e", dϕ₀, ϕ₀) next = iterate(iter) @@ -139,17 +145,23 @@ function (ls::HagerZhangLineSearch)(fg, x₀, η₀, fg₀=fg(x₀); (x, f, g, ξ, α, dϕ), state = next a, b, numfg, done = state verbosity >= 3 && - @info @sprintf("Linesearch iteration step %d, function evaluation count %d:\n[a,b] = [%.2e, %.2e], dϕᵃ = %.2e, dϕᵇ = %.2e, ϕᵃ - ϕ₀ = %.2e, ϕᵇ - ϕ₀ = %.2e", - k, numfg, a.α, b.α, a.dϕ, b.dϕ, a.ϕ - ϕ₀, b.ϕ - ϕ₀) + @info @sprintf( + "Linesearch iteration step %d, function evaluation count %d:\n[a,b] = [%.2e, %.2e], dϕᵃ = %.2e, dϕᵇ = %.2e, ϕᵃ - ϕ₀ = %.2e, ϕᵇ - ϕ₀ = %.2e", + k, numfg, a.α, b.α, a.dϕ, b.dϕ, a.ϕ - ϕ₀, b.ϕ - ϕ₀ + ) if done verbosity >= 2 && - @info @sprintf("Linesearch converged after %d iterations and %d function evaluations:\nα = %.2e, dϕ = %.2e, ϕ - ϕ₀ = %.2e", - k, numfg, α, dϕ, f - ϕ₀) + @info @sprintf( + "Linesearch converged after %d iterations and %d function evaluations:\nα = %.2e, dϕ = %.2e, ϕ - ϕ₀ = %.2e", + k, numfg, α, dϕ, f - ϕ₀ + ) return x, f, g, ξ, α, numfg elseif k >= maxiter || numfg >= maxfg verbosity >= 1 && - @warn @sprintf("Linesearch not converged after %d iterations and %d function evaluations:\nα = %.2e, dϕ = %.2e, ϕ - ϕ₀ = %.2e", - k, numfg, α, dϕ, f - ϕ₀) + @warn @sprintf( + "Linesearch not converged after %d iterations and %d function evaluations:\nα = %.2e, dϕ = %.2e, ϕ - ϕ₀ = %.2e", + k, numfg, α, dϕ, f - ϕ₀ + ) return x, f, g, ξ, α, numfg else next = iterate(iter, state) @@ -157,14 +169,15 @@ function (ls::HagerZhangLineSearch)(fg, x₀, η₀, fg₀=fg(x₀); k += 1 end end + return end # Hager-Zhang Line Search Algorithm implemented as iterator -struct HagerZhangLineSearchIterator{T₁<:Real,F₁,F₂,F₃,X,G,T₂<:Real} +struct HagerZhangLineSearchIterator{T₁ <: Real, F₁, F₂, F₃, X, G, T₂ <: Real} fdf::F₁ # computes function value and gradient for a given x, i.e. f, g = f(x) retract::F₂ # function used to step in direction η₀ with step size α, i.e. x, ξ = retract(x₀, η₀, α) where x = Rₓ₀(α*η₀) is the new position and ξ = D Rₓ₀(α*η₀)[η₀] is the derivative or tangent of x to α at the position x inner::F₃ # function used to compute inner product between gradient and direction, i.e. dϕ = inner(x, g, d); can depend on x (i.e. metric on a manifold) - p₀::LineSearchPoint{T₁,X,G} # initial position, containing x₀, f₀, g₀ + p₀::LineSearchPoint{T₁, X, G} # initial position, containing x₀, f₀, g₀ η₀::G # search direction α₀::T₁ # initial guess for step size acceptfirst::Bool # whether or not the initial guess can be accepted (e.g. LBFGS) @@ -186,38 +199,48 @@ function Base.iterate(iter::HagerZhangLineSearchIterator) ewolfe = checkexactwolfe(c, p₀, c₁, c₂) awolfe = checkapproxwolfe(c, p₀, c₁, c₂, ϵ) verbosity >= 4 && - @info @sprintf(" Linesearch initial step: c = %.2e, dϕᶜ = %.2e, ϕᶜ - ϕ₀ = %.2e, exact wolfe = %d, approx wolfe = %d", - c.α, c.dϕ, c.ϕ - p₀.ϕ, ewolfe, awolfe) + @info @sprintf( + " Linesearch initial step: c = %.2e, dϕᶜ = %.2e, ϕᶜ - ϕ₀ = %.2e, exact wolfe = %d, approx wolfe = %d", + c.α, c.dϕ, c.ϕ - p₀.ϕ, ewolfe, awolfe + ) if ewolfe || awolfe return (c.x, c.f, c.∇f, c.ξ, c.α, c.dϕ), (c, c, numfg, true) end else verbosity >= 4 && - @info @sprintf(" Linesearch initial step (cannot be accepted): c = %.2e, dϕᶜ = %.2e, ϕᶜ - ϕ₀ = %.2e", - c.α, c.dϕ, c.ϕ - p₀.ϕ) + @info @sprintf( + " Linesearch initial step (cannot be accepted): c = %.2e, dϕᶜ = %.2e, ϕᶜ - ϕ₀ = %.2e", + c.α, c.dϕ, c.ϕ - p₀.ϕ + ) end # L0 in the Line Search Algorithm: find initial bracketing interval a, b, nfg = bracket(iter, c) verbosity >= 4 && - @info @sprintf(" Linesearch initial bracket: [a,b] = [%.2e, %.2e], dϕᵃ = %.2e, dϕᵇ = %.2e, ϕᵃ - ϕ₀ = %.2e, ϕᵇ - ϕ₀ = %.2e", - a.α, b.α, a.dϕ, b.dϕ, a.ϕ - p₀.ϕ, b.ϕ - p₀.ϕ) + @info @sprintf( + " Linesearch initial bracket: [a,b] = [%.2e, %.2e], dϕᵃ = %.2e, dϕᵇ = %.2e, ϕᵃ - ϕ₀ = %.2e, ϕᵇ - ϕ₀ = %.2e", + a.α, b.α, a.dϕ, b.dϕ, a.ϕ - p₀.ϕ, b.ϕ - p₀.ϕ + ) numfg += nfg if a.α == b.α return (a.x, a.f, a.∇f, a.ξ, a.α, a.dϕ), (a, b, numfg, true) elseif (b.α - a.α) < eps(one(a.α)) verbosity >= 1 && - @warn @sprintf(" Linesearch bracket converged to a point without satisfying Wolfe conditions: [a,b] = [%.2e, %.2e], dϕᵃ = %.2e, dϕᵇ = %.2e, ϕᵃ - ϕ₀ = %.2e, ϕᵇ - ϕ₀ = %.2e", - a.α, b.α, a.dϕ, b.dϕ, a.ϕ - p₀.ϕ, b.ϕ - p₀.ϕ) + @warn @sprintf( + " Linesearch bracket converged to a point without satisfying Wolfe conditions: [a,b] = [%.2e, %.2e], dϕᵃ = %.2e, dϕᵇ = %.2e, ϕᵃ - ϕ₀ = %.2e, ϕᵇ - ϕ₀ = %.2e", + a.α, b.α, a.dϕ, b.dϕ, a.ϕ - p₀.ϕ, b.ϕ - p₀.ϕ + ) return (a.x, a.f, a.∇f, a.ξ, a.α, a.dϕ), (a, b, numfg, true) else return (a.x, a.f, a.∇f, a.ξ, a.α, a.dϕ), (a, b, numfg, false) end end -function Base.iterate(iter::HagerZhangLineSearchIterator, - state::Tuple{LineSearchPoint,LineSearchPoint,Int,Bool}) +function Base.iterate( + iter::HagerZhangLineSearchIterator, + state::Tuple{LineSearchPoint, LineSearchPoint, Int, Bool} + ) c₁ = iter.parameters.c₁ c₂ = iter.parameters.c₂ ϵ = iter.parameters.ϵ @@ -249,8 +272,10 @@ function Base.iterate(iter::HagerZhangLineSearchIterator, a, b = A, B end verbosity >= 4 && - @info @sprintf(" Linesearch updated bracket (secant2): [a,b] = [%.2e, %.2e], dϕᵃ = %.2e, dϕᵇ = %.2e, ϕᵃ - ϕ₀ = %.2e, ϕᵇ - ϕ₀ = %.2e", - a.α, b.α, a.dϕ, b.dϕ, a.ϕ - p₀.ϕ, b.ϕ - p₀.ϕ) + @info @sprintf( + " Linesearch updated bracket (secant2): [a,b] = [%.2e, %.2e], dϕᵃ = %.2e, dϕᵇ = %.2e, ϕᵃ - ϕ₀ = %.2e, ϕᵇ - ϕ₀ = %.2e", + a.α, b.α, a.dϕ, b.dϕ, a.ϕ - p₀.ϕ, b.ϕ - p₀.ϕ + ) if a.α == b.α return (a.x, a.f, a.∇f, a.ξ, a.α, a.dϕ), (a, b, numfg, true) end @@ -260,15 +285,19 @@ function Base.iterate(iter::HagerZhangLineSearchIterator, a, b, nfg = update(iter, a, b, (a.α + b.α) / 2) numfg += nfg verbosity >= 4 && - @info @sprintf(" Linesearch updated bracket (bisection): [a,b] = [%.2e, %.2e], dϕᵃ = %.2e, dϕᵇ = %.2e, ϕᵃ - ϕ₀ = %.2e, ϕᵇ - ϕ₀ = %.2e", - a.α, b.α, a.dϕ, b.dϕ, a.ϕ - p₀.ϕ, b.ϕ - p₀.ϕ) + @info @sprintf( + " Linesearch updated bracket (bisection): [a,b] = [%.2e, %.2e], dϕᵃ = %.2e, dϕᵇ = %.2e, ϕᵃ - ϕ₀ = %.2e, ϕᵇ - ϕ₀ = %.2e", + a.α, b.α, a.dϕ, b.dϕ, a.ϕ - p₀.ϕ, b.ϕ - p₀.ϕ + ) end if a.α == b.α return (a.x, a.f, a.∇f, a.ξ, a.α, a.dϕ), (a, b, numfg, true) elseif (b.α - a.α) < eps(one(a.α)) verbosity >= 1 && - @warn @sprintf(" Linesearch bracket converged to a point without satisfying Wolfe conditions: [a,b] = [%.2e, %.2e], dϕᵃ = %.2e, dϕᵇ = %.2e, ϕᵃ - ϕ₀ = %.2e, ϕᵇ - ϕ₀ = %.2e", - a.α, b.α, a.dϕ, b.dϕ, a.ϕ - p₀.ϕ, b.ϕ - p₀.ϕ) + @warn @sprintf( + " Linesearch bracket converged to a point without satisfying Wolfe conditions: [a,b] = [%.2e, %.2e], dϕᵃ = %.2e, dϕᵇ = %.2e, ϕᵃ - ϕ₀ = %.2e, ϕᵇ - ϕ₀ = %.2e", + a.α, b.α, a.dϕ, b.dϕ, a.ϕ - p₀.ϕ, b.ϕ - p₀.ϕ + ) return (a.x, a.f, a.∇f, a.ξ, a.α, a.dϕ), (a, b, numfg, true) else return (a.x, a.f, a.∇f, a.ξ, a.α, a.dϕ), (a, b, numfg, false) @@ -291,8 +320,10 @@ secant(a, b, fa, fb) = (a * fb - b * fa) / (fb - fa) # update bracketing interval: interval has (a.dϕ < 0, a.ϕ <= f₀+ϵ), (b.dϕ >= 0) # the bisection step is separtely implemented in the bisect function -function update(iter::HagerZhangLineSearchIterator, a::LineSearchPoint, b::LineSearchPoint, - αc) +function update( + iter::HagerZhangLineSearchIterator, a::LineSearchPoint, b::LineSearchPoint, + αc + ) p₀ = iter.p₀ c₁ = iter.parameters.c₁ c₂ = iter.parameters.c₂ @@ -311,8 +342,10 @@ function update(iter::HagerZhangLineSearchIterator, a::LineSearchPoint, b::LineS ewolfe = checkexactwolfe(c, p₀, c₁, c₂) awolfe = checkapproxwolfe(c, p₀, c₁, c₂, ϵ) verbosity >= 4 && - @info @sprintf(" Linesearch update step (U1-U2): c = %.2e, dϕᶜ = %.2e, ϕᶜ - ϕ₀ = %.2e, exact wolfe = %d, approx wolfe = %d", - c.α, c.dϕ, c.ϕ - ϕ₀, ewolfe, awolfe) + @info @sprintf( + " Linesearch update step (U1-U2): c = %.2e, dϕᶜ = %.2e, ϕᶜ - ϕ₀ = %.2e, exact wolfe = %d, approx wolfe = %d", + c.α, c.dϕ, c.ϕ - ϕ₀, ewolfe, awolfe + ) if ewolfe || awolfe return c, c, 1 end @@ -338,8 +371,10 @@ function bisect(iter::HagerZhangLineSearchIterator, a::LineSearchPoint, b::LineS while true if (b.α - a.α) <= eps(one(a.α))^(3 // 4) || numfg >= iter.parameters.maxfg if verbosity >= 1 - @warn @sprintf(" Linesearch bisection failure: [a, b] = [%.2e, %.2e], b-a = %.2e, dϕᵃ = %.2e, dϕᵇ = %.2e, (ϕᵇ - ϕᵃ)/(b-a) = %.2e", - a.α, b.α, b.α - a.α, a.dϕ, b.dϕ, (b.ϕ - a.ϕ) / (b.α - a.α)) + @warn @sprintf( + " Linesearch bisection failure: [a, b] = [%.2e, %.2e], b-a = %.2e, dϕᵃ = %.2e, dϕᵇ = %.2e, (ϕᵇ - ϕᵃ)/(b-a) = %.2e", + a.α, b.α, b.α - a.α, a.dϕ, b.dϕ, (b.ϕ - a.ϕ) / (b.α - a.α) + ) end return a, b, numfg end @@ -349,8 +384,10 @@ function bisect(iter::HagerZhangLineSearchIterator, a::LineSearchPoint, b::LineS ewolfe = checkexactwolfe(d, p₀, c₁, c₂) awolfe = checkapproxwolfe(d, p₀, c₁, c₂, ϵ) verbosity >= 4 && - @info @sprintf(" Linesearch bisection update (U3): d = %.2e, dϕᵈ = %.2e, ϕᵈ - ϕ₀ = %.2e, exact wolfe = %d, approx wolfe = %d", - d.α, d.dϕ, d.ϕ - p₀.ϕ, ewolfe, awolfe) + @info @sprintf( + " Linesearch bisection update (U3): d = %.2e, dϕᵈ = %.2e, ϕᵈ - ϕ₀ = %.2e, exact wolfe = %d, approx wolfe = %d", + d.α, d.dϕ, d.ϕ - p₀.ϕ, ewolfe, awolfe + ) if ewolfe || awolfe return d, d, numfg end @@ -362,6 +399,7 @@ function bisect(iter::HagerZhangLineSearchIterator, a::LineSearchPoint, b::LineS b = d end end + return end # bracket function for finding an initial bracketing interval from a given initial step @@ -382,25 +420,30 @@ function bracket(iter::HagerZhangLineSearchIterator{T}, c::LineSearchPoint) wher c = takestep(iter, α) numfg += 1 verbosity >= 4 && - @info @sprintf(" Linesearch bracket step: c = %.2e, dϕᶜ = %.2e, ϕᶜ - ϕ₀ = %.2e", - c.α, c.dϕ, c.ϕ - p₀.ϕ) + @info @sprintf( + " Linesearch bracket step: c = %.2e, dϕᶜ = %.2e, ϕᶜ - ϕ₀ = %.2e", + c.α, c.dϕ, c.ϕ - p₀.ϕ + ) end - c.dϕ >= 0 && return a, c, numfg# B1 + c.dϕ >= 0 && return a, c, numfg # B1 # from here: c.dϕ < 0 if c.ϕ > fmax # B2 a, b, nfg = bisect(iter, iter.p₀, c) return a, b, numfg + nfg - else# B3 + else # B3 a = c α *= iter.parameters.ρ c = takestep(iter, α) numfg += 1 verbosity >= 4 && - @info @sprintf(" Linesearch bracket step: c = %.2e, dϕᶜ = %.2e, ϕᶜ - ϕ₀ = %.2e", - c.α, c.dϕ, c.ϕ - p₀.ϕ) + @info @sprintf( + " Linesearch bracket step: c = %.2e, dϕᶜ = %.2e, ϕᶜ - ϕ₀ = %.2e", + c.α, c.dϕ, c.ϕ - p₀.ϕ + ) # if checkexactwolfe(c, p₀, c₁, c₂) || checkapproxwolfe(c, p₀, c₁, c₂, ϵ) # return c, c, numfg # end end end + return end diff --git a/src/linesearches2.jl b/src/linesearches2.jl new file mode 100644 index 0000000..8473ddd --- /dev/null +++ b/src/linesearches2.jl @@ -0,0 +1,40 @@ +using LineSearches: LineSearches + +struct LineSearchWrapper{A <: LineSearches.AbstractLineSearch} <: OptimKit.AbstractLineSearch + alg::A +end + +function make_ϕdϕ(fg, retract, inner, x₀, η₀) + function ϕ(α) + x, ξ = retract(x₀, η₀, α) + f, ∇f = fg(x) + return f + end + function dϕ(α) + x, ξ = retract(x₀, η₀, α) + f, ∇f = fg(x) + return inner(x, ∇f, ξ) + end + function ϕdϕ(α) + x, ξ = retract(x₀, η₀, α) + f, ∇f = fg(x) + return f, inner(x, ∇f, ξ) + end + return ϕ, dϕ, ϕdϕ +end + +function (ls::LineSearchWrapper)( + fg, x₀, η₀, fg₀ = fg(x₀); + retract = OptimKit._retract, inner = OptimKit._inner, + initialguess::Real = one(fg₀[1]), + acceptfirst::Bool = false + ) + ϕ, dϕ, ϕdϕ = make_ϕdϕ(fg, retract, inner, x₀, η₀) + ϕ₀ = fg₀[1] + dϕ₀ = inner(x₀, fg₀[2], η₀) + α, _ = ls.alg(ϕ, dϕ, ϕdϕ, initialguess, ϕ₀, dϕ₀) + x, ξ = retract(x₀, η₀, α) + f, g = fg(x) + numfg = 0 + return x, f, g, ξ, α, numfg +end diff --git a/src/simpleiteration.jl b/src/simpleiteration.jl index 8c4eccd..c2908e5 100644 --- a/src/simpleiteration.jl +++ b/src/simpleiteration.jl @@ -24,9 +24,10 @@ struct SimpleIteration{T <: Real} <: FixedPointAlgorithm verbosity::Int end function SimpleIteration(; - maxiter::Int=MAXITER[], - gradtol::Real=GRADTOL[], - verbosity::Int=VERBOSITY[]) + maxiter::Int = MAXITER[], + gradtol::Real = GRADTOL[], + verbosity::Int = VERBOSITY[] + ) return SimpleIteration(maxiter, gradtol, verbosity) end diff --git a/test/runtests.jl b/test/runtests.jl index 4904773..837b128 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,35 +5,41 @@ using Random: Random Random.seed!(1234) # Test linesearches -@testset "Linesearch" for (fg, x₀) in [(x -> (sin(x) + x^4, cos(x) + 4 * x^3), 0.0), - (x -> (x^2, 2 * x), 1.0), - (x -> (exp(x) - x^3, exp(x) - 3 * x^2), 3.9), # should trigger bisection - (x -> (exp(x) - x^3, exp(x) - 3 * x^2), 2), # should trigger infinities - (x -> (sum(x .^ 2), 2 * x), [1.0, 2.0]), - (x -> (2 * x[1]^2 + x[2]^4 - x[3]^2 + x[3]^4, - [4 * x[1], 4 * x[2]^3, -2 * x[3] + 4 * x[3]^3]), - [1.0, 2.0, 3.0])] +@testset "Linesearch" for (fg, x₀) in [ + (x -> (sin(x) + x^4, cos(x) + 4 * x^3), 0.0), + (x -> (x^2, 2 * x), 1.0), + (x -> (exp(x) - x^3, exp(x) - 3 * x^2), 3.9), # should trigger bisection + (x -> (exp(x) - x^3, exp(x) - 3 * x^2), 2), # should trigger infinities + (x -> (sum(x .^ 2), 2 * x), [1.0, 2.0]), + ( + x -> ( + 2 * x[1]^2 + x[2]^4 - x[3]^2 + x[3]^4, + [4 * x[1], 4 * x[2]^3, -2 * x[3] + 4 * x[3]^3], + ), + [1.0, 2.0, 3.0], + ), + ] f₀, g₀ = fg(x₀) for i in 1:100 c₁ = 0.5 * rand() c₂ = 0.5 + 0.5 * rand() - ls = HagerZhangLineSearch(; c₁=c₁, c₂=c₂, ϵ=0, ρ=1.5, maxfg=100, maxiter=100) - x, f, g, ξ, α, numfg = ls(fg, x₀, -g₀; verbosity=4) + ls = HagerZhangLineSearch(; c₁ = c₁, c₂ = c₂, ϵ = 0, ρ = 1.5, maxfg = 100, maxiter = 100) + x, f, g, ξ, α, numfg = ls(fg, x₀, -g₀; verbosity = 4) @test f ≈ fg(x)[1] @test g ≈ fg(x)[2] @test ξ == -g₀ @test dot(ξ, g) >= c₂ * dot(ξ, g₀) @test f <= f₀ + α * c₁ * dot(ξ, g₀) || (2 * c₁ - 1) * dot(ξ, g₀) > dot(ξ, g) - x, f, g, ξ, α, numfg = ls(fg, x₀, -g₀; initialguess=1e-4, verbosity=2) # test extrapolation phase + x, f, g, ξ, α, numfg = ls(fg, x₀, -g₀; initialguess = 1.0e-4, verbosity = 2) # test extrapolation phase @test f ≈ fg(x)[1] @test g ≈ fg(x)[2] @test ξ == -g₀ @test dot(ξ, g) >= c₂ * dot(ξ, g₀) @test f <= f₀ + α * c₁ * dot(ξ, g₀) || (2 * c₁ - 1) * dot(ξ, g₀) > dot(ξ, g) - x, f, g, ξ, α, numfg = ls(fg, x₀, -g₀; initialguess=1e4) # test infinities + x, f, g, ξ, α, numfg = ls(fg, x₀, -g₀; initialguess = 1.0e4) # test infinities @test f ≈ fg(x)[1] @test g ≈ fg(x)[2] @test ξ == -g₀ @@ -72,10 +78,10 @@ algorithms = (GradientDescent, ConjugateGradient, LBFGS) A = A' * A fg = quadraticproblem(A, y) x₀ = randn(n) - alg = algtype(; verbosity=2, gradtol=1e-12, maxiter=10_000_000) + alg = algtype(; verbosity = 2, gradtol = 1.0e-12, maxiter = 10_000_000) x, f, g, numfg, normgradhistory = optimize(fg, x₀, alg) - @test x ≈ y rtol = cond(A) * 1e-12 - @test f < 1e-12 + @test x ≈ y rtol = cond(A) * 1.0e-12 + @test f < 1.0e-12 n = 1000 y = randn(n) @@ -85,11 +91,40 @@ algorithms = (GradientDescent, ConjugateGradient, LBFGS) # well conditioned, all eigenvalues between 1 and 2 fg = quadratictupleproblem(A' * A, (y, 1.0)) x₀ = (randn(n), 2.0) - alg = algtype(; verbosity=3, gradtol=1e-8) + alg = algtype(; verbosity = 3, gradtol = 1.0e-8) x, f, g, numfg, normgradhistory = optimize(fg, x₀, alg) - @test x[1] ≈ y rtol = 1e-7 - @test x[2] ≈ 1 rtol = 1e-7 - @test f < 1e-12 + @test x[1] ≈ y rtol = 1.0e-7 + @test x[2] ≈ 1 rtol = 1.0e-7 + @test f < 1.0e-12 +end + +include("sphere.jl") + +@testset "Manifold correctness" begin + D = Diagonal(collect(0.1:0.1:1.0)) + x₀ = randn(10) + x₀ = x₀ / norm(x₀) + function fg(x) + f = -dot(x, D * x) / 2 + g = -D * x - 2 * f * x + return f, g + end + alg = GradientDescent(; verbosity = 0, gradtol = 1.0e-6, maxiter = 10) + x, f, g, numfg, normgradhistory = optimize(fg, x₀, alg; retract = Sphere.retract, transport! = Sphere.transport) + alg = ConjugateGradient(; verbosity = 0, gradtol = 1.0e-6, maxiter = 10) + x, f, g, numfg, normgradhistory = optimize(fg, x₀, alg; retract = Sphere.retract, transport! = Sphere.transport) + alg = LBFGS(5; verbosity = 0, gradtol = 1.0e-6, maxiter = 10) + x, f, g, numfg, normgradhistory = optimize(fg, x₀, alg; retract = Sphere.retract, transport! = Sphere.transport) + + function fp(x) + Dx = D * x + return Dx / norm(Dx) + end + + alg = SimpleIteration(; verbosity = 0, gradtol = 1.0e-6, maxiter = 10) + x, g, numfp, normgradhistory = fixedpoint(fp, x₀, alg; retract = Sphere.retract, transport! = Sphere.transport, invretract = Sphere.invretract) + alg = AndersonMixing(5; verbosity = 0, gradtol = 1.0e-6, maxiter = 10) + x, g, numfp, normgradhistory = fixedpoint(fp, x₀, alg; retract = Sphere.retract, transport! = Sphere.transport, invretract = Sphere.invretract) end @testset "Aqua" verbose = true begin diff --git a/test/sphere.jl b/test/sphere.jl new file mode 100644 index 0000000..bccfb5f --- /dev/null +++ b/test/sphere.jl @@ -0,0 +1,42 @@ +module Sphere +using Test +using VectorInterface +function retract(x, g, α) + @test norm(x) ≈ 1 + @test abs(inner(x, g)) < sqrt(eps(eltype(x))) + gnorm = norm(g) + ug = g / gnorm + θ = α * gnorm + xnew = cos(θ) * x + sin(θ) * ug + gnew = -sin(θ) * gnorm * x + cos(θ) * g + gnew = add!(gnew, xnew, -inner(xnew, gnew)) + @test norm(xnew) ≈ 1 + @test abs(inner(xnew, gnew)) < sqrt(eps(eltype(x))) + return xnew, gnew +end +function invretract(x₀, x₁) + @test norm(x₀) ≈ 1 + @test norm(x₁) ≈ 1 + cosθ = inner(x₀, x₁) + θ = acos(clamp(cosθ, -1, 1)) + g = (x₁ - cosθ * x₀) / sinc(θ / pi) + return g +end +function transport(v, x₀, g₀, α, x₁) + @test norm(x₀) ≈ 1 + @test norm(x₁) ≈ 1 + @test abs(inner(x₀, v)) < sqrt(eps(eltype(x₀))) + @test abs(inner(x₀, g₀)) < sqrt(eps(eltype(x₀))) + @test retract(x₀, g₀, α)[1] ≈ x₁ + gnorm = norm(g₀) + θ = α * gnorm + ug = g₀ / gnorm + gv = inner(ug, v) + v₁ = v - gv * ug + vnew = gv * (-sin(θ) * x₀ + cos(θ) * ug) + v₁ + @test abs(inner(x₁, vnew)) < sqrt(eps(eltype(x₀))) + @test norm(vnew) ≈ norm(v) + vnew = add!(vnew, x₁, -inner(x₁, vnew)) + return vnew +end +end