Linesearches for Optimizers

In GeometricOptimizers we typically build the search direction by multiplying the gradient with a Hessian. When starting at $x_k$ we take:

\[ p_k = H_{x_k}^{-1}(\nabla_{x_k}f),\]

where $[H_{x_k}]_{ij} = \partial^2{}f\partial{}x_i\partial{}x_j|_{x_k}$ is the Hessian. Note that we often use approximations of this Hessian in practice (such as the HessianBFGS).

The linesearch objective is then built as

\[f^\mathrm{ls}(\alpha) = f(x_k + \alpha{}p_k).\]

For manifolds [8] defining a Hessian, equivalently to defining a gradient, requires a Riemannian metric and the associated Levi-Civita connection $\nabla$:

\[\mathrm{Hess}(f) := \nabla\nabla{}f = \nabla{}df \in \Gamma(T^*\mathcal{M}\otimes{}T^*\mathcal{M}).\]

For specific vector fields $\xi, \eta \in \Gamma(T\mathcal{M})$ we can write this as:

\[\langle \mathrm{Hess}(f)[\xi], \eta \rangle = \xi(\eta{}f) - (\nabla_\xi\eta)f.\]

Example

We look at the following example:

f(x::Union{T, Vector{T}}) where {T<:Number} = exp.(x) .* (x .^ 3 .- 5x .+ 2x) .+ 2one(T)
f!(y::AbstractVector{T}, x::AbstractVector{T}) where {T} = y .= f.(x)
F!(y::AbstractVector{T}, x::AbstractVector{T}, params) where {T} = f!(y, x)

F! (generic function with 1 method)

We hence use linesearch_problem not for a SimpleSolvers.NewtonSolver, but for an Optimizer:

x₀ = [0., .1, .2]
x = copy(x₀)
obj = OptimizerProblem(sum∘f, x₀)
grad = GradientAutodiff{Float64}(obj.F, length(x₀))
_cache = NewtonOptimizerCache(x₀)
state = NewtonOptimizerState(x₀)
hess = HessianAutodiff(obj, x₀)
update!(state, grad, x₀)
update!(_cache, state, grad, hess, x₀)
params = (x = state.x, state = state)
# the retraction is how a trial step is taken; on an `AbstractVector` like this one it is
# never consulted, but `linesearch_problem` needs it for the manifold case
ls_obj = linesearch_problem(obj, grad, _cache, Cayley())

fˡˢ(alpha) = ls_obj.F(alpha, params)
∂fˡˢ∂α(alpha) = ls_obj.D(alpha, params)

We now again want to find the minimum with quadratic line search and repeat the procedure above:

p₀ = fˡˢ(0.)

4.9464834626645615

p₁ = ∂fˡˢ∂α(0.)

6.452258296067367

params = (x = state.x, state = state)
α₀ = bracket_minimum_with_fixed_point(ls_obj, params, 0.)[1]
y = fˡˢ(α₀)
p₂ = (y - p₀ - p₁*α₀) / α₀^2
p(α) = p₀ + p₁ * α + p₂ * α^2
α₁ = -p₁ / (2p₂)

-1.2922280958762398

We now again move the original $x$ in the Newton direction with step length $\alpha_1$:

sum∘f(x)

sum ∘ [2.0, 1.6695538954953815, 1.2769295671691796]

compute_new_iterate!(x, α₁, direction(_cache))

3-element Vector{Float64}:
 0.6461140479381199
 0.8491210578860486
 1.0708493689600747

sum∘f(x)

sum ∘ [-1.183900159425661, -2.5235695736019492, -3.7907361578835586]