fix a few issues with quasi-Newton

Author: mateuszbaran

src/plans/box_plan.jl

@@ -85,6 +85,7 @@ end
 
 function initialize_update!(ha::QuasiNewtonLimitedMemoryBoxDirectionUpdate)
     initialize_update!(ha.qn_du)
+    ha.last_gcp_result = :not_searched
     return ha
 end
 
@@ -652,8 +653,7 @@ function find_generalized_cauchy_direction!(gcp::GeneralizedCauchyDirectionFinde
     dt_min = max(dt_min, 0.0)
     t_old = t_old + dt_min
     d_out .*= t_old
-    set_stepsize_bound!(M, d_out, p, ts, t_current)
-
+    set_stepsize_bound!(M, d_out, p, ts, t_old)
     if has_finite_limit
         return :found_limited
     else

src/solvers/quasi_Newton.jl

@@ -389,6 +389,15 @@ function quasi_Newton!(
 end
 calls_with_kwargs(::typeof(quasi_Newton!)) = (decorate_objective!, decorate_state!)
 
+function _get_max_stepsize(M::AbstractManifold, qns::QuasiNewtonState)
+    current_max_stepsize = get_parameter(qns.direction_update, Val(:max_stepsize))
+    if !isnothing(current_max_stepsize) && !isfinite(current_max_stepsize)
+        return max_stepsize(M, qns.p) / norm(qns.η)
+    else
+        return current_max_stepsize
+    end
+end
+
 function initialize_solver!(amp::AbstractManoptProblem, qns::QuasiNewtonState)
     M = get_manifold(amp)
     get_gradient!(amp, qns.X, qns.p)
@@ -401,10 +410,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_stepsize = get_parameter(qns.direction_update, Val(:max_stepsize))
-    if !isnothing(current_max_stepsize) && !isfinite(current_max_stepsize)
-        current_max_stepsize = max_stepsize(M, qns.p) / norm(qns.η)
-    end
+    current_max_stepsize = _get_max_stepsize(M, qns)
     if !(qns.nondescent_direction_behavior === :ignore)
         qns.nondescent_direction_value = real(inner(M, qns.p, qns.η, qns.X))
         if qns.nondescent_direction_value >= 0
@@ -416,6 +422,12 @@ function step_solver!(mp::AbstractManoptProblem, qns::QuasiNewtonState, k)
             if qns.nondescent_direction_behavior === :reinitialize_direction_update
                 initialize_update!(qns.direction_update)
             end
+            # update direction after reinitialization to get a valid one
+            if qns.nondescent_direction_behavior === :step_towards_negative_gradient ||
+                    qns.nondescent_direction_behavior === :reinitialize_direction_update
+                qns.direction_update(qns.η, mp, qns)
+                current_max_stepsize = _get_max_stepsize(M, qns)
+            end
         end
     end
     local α # COV_EXCL_LINE
@@ -802,6 +814,17 @@ function update_hessian!(
     end
 
     if reforming_required
+        # we need to move first vectors in memory too because they most likely won't be
+        # overwritten by new pairs
+        if start == 2
+            vector_transport_to!(
+                M, d.memory_s[1], p_old, d.memory_s[1], p, d.vector_transport_method
+            )
+            vector_transport_to!(
+                M, d.memory_y[1], p_old, d.memory_y[1], p, d.vector_transport_method
+            )
+            fill_rho_i!(M, p, d, 1)
+        end
         _drop_zero_rho_vectors!(d)
     end