start integration of GCP into quasi_Newton

Author: mateuszbaran

docs/src/references.bib

@@ -260,6 +260,20 @@ @book{Boumal:2023
     ISBN      = {978-1-00-916616-4}
 }
 
+@article{ByrdNocedalSchnabel:1994,
+	title = {Representations of quasi-{Newton} matrices and their use in limited memory methods},
+	volume = {63},
+	issn = {1436-4646},
+	doi = {10.1007/BF01582063},
+	number = {1},
+	urldate = {2025-09-06},
+	journal = {Mathematical Programming},
+	author = {Byrd, Richard H. and Nocedal, Jorge and Schnabel, Robert B.},
+	month = jan,
+	year = {1994},
+	pages = {129--156},
+}
+
 @article{ByrdLuNocedalZhu:1995,
 	title = {A {Limited} {Memory} {Algorithm} for {Bound} {Constrained} {Optimization}},
 	volume = {16},

src/helpers/LineSearchesTypes.jl

@@ -14,7 +14,8 @@ Wrapper for line searches available in the `LineSearches.jl` library.
 
 Wrap `linesearch` (for example [`HagerZhang`](https://julianlsolvers.github.io/LineSearches.jl/latest/reference/linesearch.html#LineSearches.HagerZhang)
 or [`MoreThuente`](https://julianlsolvers.github.io/LineSearches.jl/latest/reference/linesearch.html#LineSearches.MoreThuente)).
-The initial step selection from Linesearches.jl is not yet supported and the value 1.0 is used.
+The initial step selection from Linesearches.jl is not yet supported and `initial_guess` is
+always used (by default [`ConstantInitialGuess`](@ref)).
 
 # Keyword Arguments
 

src/plans/box_plan.jl

@@ -7,9 +7,24 @@ with a standard manifold. Otherwise return `false`.
 requires_gcp(::AbstractManifold) = false
 requires_gcp(M::ProductManifold) = requires_gcp(M.manifolds[1])
 
+@doc raw"""
+    mutable struct LimitedMemoryHessianApproximation end
+
+An approximation of Hessian of a scalar function of the form ``B_0 = θ I``,
+``B_{k+1} = B_k - W_k M_k W_k^{\mathrm{T}}``,
+where ``\theta > 0`` is an initial scaling guess.
+Matrix ``M_k = \begin{psmallmatrix}M_{11} & M_{21}^{\mathrm{T}}\\ M_{21} & M_{22}\end{psmallmatrix}``
+is stored using its blocks.
+Blocks ``W_k`` are (implicitly) composed from `memory_y` and `memory_s`.
+
+Initial scale ``\theta`` is stored in the field `initial_scale` but if the memory isn't empty,
+the current scale is set to squared norm of $s_k$ divided by inner product of ``s_k`` and ``y_k``
+where ``k`` is the oldest index for which the denominator is not equal to 0.
+
+See [ByrdNocedalSchnabel:1994](@cite) for details.
+"""
 mutable struct QuasiNewtonLimitedMemoryBoxDirectionUpdate{
         TDU <: QuasiNewtonLimitedMemoryDirectionUpdate,
-        TM <: AbstractManifold,
         F <: Real,
         T_HM <: AbstractMatrix,
         V <: AbstractVector,
@@ -22,18 +37,60 @@ mutable struct QuasiNewtonLimitedMemoryBoxDirectionUpdate{
     M_11::T_HM
     M_21::T_HM
     M_22::T_HM
-    # buffer for calculating stuff
+    # buffer for calculating W_k blocks
     coords_Sk_X::V
     coords_Sk_Y::V
     coords_Yk_X::V
     coords_Yk_Y::V
 end
 
+function QuasiNewtonLimitedMemoryBoxDirectionUpdate(
+        qn_du::QuasiNewtonLimitedMemoryDirectionUpdate{<:AbstractQuasiNewtonUpdateRule, T, F}
+    ) where {T, F <: Real}
+    memory_size = capacity(qn_du.memory_s)
+    M_11 = zeros(F, memory_size, memory_size)
+    M_21 = zeros(F, memory_size, memory_size)
+    M_22 = zeros(F, memory_size, memory_size)
+    coords_Sk_X = zeros(F, memory_size)
+    coords_Sk_Y = zeros(F, memory_size)
+    coords_Yk_X = zeros(F, memory_size)
+    coords_Yk_Y = zeros(F, memory_size)
+    return QuasiNewtonLimitedMemoryBoxDirectionUpdate{
+        typeof(qn_du), F, typeof(M_11), typeof(coords_Sk_X),
+    }(
+        qn_du,
+        qn_du.initial_scale,
+        M_11,
+        M_21,
+        M_22,
+        coords_Sk_X,
+        coords_Sk_Y,
+        coords_Yk_X,
+        coords_Yk_Y,
+    )
+end
+
 function initialize_update!(ha::QuasiNewtonLimitedMemoryBoxDirectionUpdate)
     initialize_update!(ha.qn_du)
     return ha
 end
 
+function (d::QuasiNewtonLimitedMemoryBoxDirectionUpdate)(
+        mp::AbstractManoptProblem, st
+    )
+    r = zero_vector(get_manifold(mp), get_iterate(st))
+    return d(r, mp, st)
+end
+function (d::QuasiNewtonLimitedMemoryBoxDirectionUpdate)(
+        r, mp::AbstractManoptProblem, st
+    )
+    d.qn_du(r, mp, st)
+    # TODO find_gcp_direction!
+    return r
+end
+
+get_update_vector_transport(u::QuasiNewtonLimitedMemoryBoxDirectionUpdate) = get_update_vector_transport(u.qn_du)
+
 function get_at_bound_index(M::ProductManifold, X, b)
     return get_at_bound_index(M.manifolds[1], submanifold_component(M, X, Val(1)), b)
 end
@@ -306,8 +363,8 @@ function (::GenericFPFPPUpdater)(M::AbstractManifold, p, old_f_prime, old_f_doub
 end
 
 function init_updater!(M::AbstractManifold, fpfpp_upd::LimitedMemoryFPFPPUpdater, d, ha::QuasiNewtonLimitedMemoryBoxDirectionUpdate)
-    fpfpp_upd.c_s .= 0
-    fpfpp_upd.c_y .= 0
+    fill!(fpfpp_upd.c_s, 0)
+    fill!(fpfpp_upd.c_y, 0)
     ii = 1
     for i in eachindex(ha.ρ)
         if iszero(ha.ρ[i])
@@ -420,7 +477,7 @@ function GCPFinder(
 end
 
 """
-    find_gcp!(gcp::GCPFinder, d_out, p, d, X)
+    find_gcp_direction!(gcp::GCPFinder, d_out, p, d, X)
 
 Find generalized Cauchy point looking from point `p` in direction `d` and save the tangent
 vector pointing at it to `d_out`. Gradient of the objective at `p` is `X`.

src/plans/stepsize/stepsize.jl

@@ -114,13 +114,19 @@ function (a::ArmijoLinesearchStepsize)(
     )
     p = get_iterate(s)
     grad = isnothing(gradient) ? get_gradient(mp, get_iterate(s)) : gradient
-    return a(mp, p, grad, η; initial_guess = a.initial_guess(mp, s, k, a.last_stepsize, η))
+    return a(mp, p, grad, η; initial_guess = a.initial_guess(mp, s, k, a.last_stepsize, η), kwargs...)
 end
 function (a::ArmijoLinesearchStepsize)(
-        mp::AbstractManoptProblem, p, X, η; initial_guess = 1.0, kwargs...
+        mp::AbstractManoptProblem, p, X, η; initial_guess::Real = 1.0, kwargs...
     )
     reset_messages!(a.messages)
     l = norm(get_manifold(mp), p, η)
+    local swse
+    if :stop_when_stepsize_exceeds in keys(kwargs)
+        swse = kwargs.stop_when_stepsize_exceeds
+    else
+        swse = (a.stop_when_stepsize_exceeds / l)
+    end
     a.last_stepsize = linesearch_backtrack!(
         get_manifold(mp),
         a.candidate_point,
@@ -133,7 +139,7 @@ function (a::ArmijoLinesearchStepsize)(
         gradient = X,
         retraction_method = a.retraction_method,
         stop_when_stepsize_less = (a.stop_when_stepsize_less / l),
-        stop_when_stepsize_exceeds = (a.stop_when_stepsize_exceeds / l),
+        stop_when_stepsize_exceeds = swse,
         stop_increasing_at_step = a.stop_increasing_at_step,
         stop_decreasing_at_step = a.stop_decreasing_at_step,
         additional_decrease_condition = a.additional_decrease_condition,
@@ -1251,7 +1257,7 @@ mutable struct NonmonotoneLinesearchStepsize{
             retraction_method::TRM = default_retraction_method(M),
             stepsize_reduction::R = 0.5,
             stop_when_stepsize_less::R = 0.0,
-            stop_when_stepsize_exceeds = real(max_stepsize(M)),
+            stop_when_stepsize_exceeds::R = real(max_stepsize(M)),
             stop_increasing_at_step::I = 100,
             stop_decreasing_at_step::I = 1000,
             storage::Union{Nothing, StoreStateAction} = StoreStateAction(

src/solvers/quasi_Newton.jl

@@ -339,6 +339,7 @@ function quasi_Newton!(
         O <: Union{AbstractManifoldFirstOrderObjective{E}, AbstractDecoratedManifoldObjective{E}},
     }
     keywords_accepted(quasi_Newton!; kwargs...)
+    local local_dir_upd
     if memory_size >= 0
         local_dir_upd = QuasiNewtonLimitedMemoryDirectionUpdate(
             M,
@@ -351,6 +352,9 @@ function quasi_Newton!(
             nonpositive_curvature_behavior = nonpositive_curvature_behavior,
             sy_tol = sy_tol,
         )
+        if requires_gcp(M)
+            local_dir_upd = QuasiNewtonLimitedMemoryBoxDirectionUpdate(local_dir_upd)
+        end
     else
         local_dir_upd = QuasiNewtonMatrixDirectionUpdate(
             M,
@@ -401,6 +405,7 @@ function step_solver!(mp::AbstractManoptProblem, qns::QuasiNewtonState, k)
     M = get_manifold(mp)
     get_gradient!(mp, qns.X, qns.p)
     qns.direction_update(qns.η, mp, qns)
+    # current_max_length = get_parameter(qns.direction_update, Val(:max_length))
     if !(qns.nondescent_direction_behavior === :ignore)
         qns.nondescent_direction_value = real(inner(M, qns.p, qns.η, qns.X))
         if qns.nondescent_direction_value > 0

test/solvers/test_quasi_Newton_box.jl

@@ -116,7 +116,7 @@ using LinearAlgebra: I, eigvecs, tr, Diagonal, dot
         f(M, p) = sum(p .^ 2)
         grad_f(M, p) = 2 .* p
         p0 = [0.0, 4.0, 10.0]
-        p_opt = quasi_Newton(M, f, grad_f, p0)
+        # p_opt = quasi_Newton(M, f, grad_f, p0)
     end
 
     @testset "requires_gcp" begin