LBFGS Settings

src/api/cpd.jl

@@ -134,7 +134,7 @@ function _component_energy_summary(energies::AbstractVector{<:Real})
         top1 = ordered[1],
         top2 = sum(@view ordered[1:min(2, r)]),
         top3 = sum(@view ordered[1:min(3, r)]),
-        effective = hhi > 0 ? inv(hhi) : NaN,
+        effective = hhi > 0 ? one(hhi) / hhi : NaN,
         argmax_component = argmax(shares),
     )
 end
@@ -1000,7 +1000,8 @@ If `r` is omitted, uses the smallest tensor mode as a heuristic rank.
     - `warm_init = TuckerInit()`: Before finding the warm start initial point, this sets the good starting point for ALS.
     - `warm_steps = 500`: Once finding the best initial point from warm_init, it runs this many ALS iterations to refine the initial point.
 * `maxiter = 500`: Maximum number of Riemannian gradient descent iterations.
-* `stepsize = 1.0`: Initial step size for line search in Riemannian gradient descent.
+* `stepsize`: Initial step size for line search in Riemannian gradient descent. Defaults to `1.0` for ordinary CPD and `0.01` for the nonnegative route.
+* `armijo_alpha_min = 1e-8`: Minimum Armijo line-search step size for `solver = :rgd`.
 * `tol = 1e-6`: Convergence tolerance.
 * `gradient_mode = :riemannian`: Gradient rule for manifold solvers. 
     - If the model has a direct rgrad, it uses that.
@@ -1022,7 +1023,7 @@ If `r` is omitted, uses the smallest tensor mode as a heuristic rank.
 * `:squaring_metric` and `:softplus_metric` require `nonnegative = true`.
 * When `nonnegative = true`, `cpd(...)` routes to `nncpd(...)`. In that route:
     - if `solver != :als` and `geometry` is left at `:canonical`, the effective geometry becomes `:softplus_metric`
-    - if `stepsize` is left at `1.0`, the effective default becomes `0.01`
+    - if `stepsize` is left as `nothing`, the effective default becomes `0.01`
     - if `init = :tucker`, the effective initializer becomes `:alswarm`
 
 ## Example 
@@ -1046,7 +1047,7 @@ function cpd(
     solver = :rgd,
     geometry = :canonical,
     maxiter = 500,
-    stepsize = 1.0,
+    stepsize = nothing,
     tol = 1e-6,
     gradient_mode = :riemannian,
     normalization = :auto,
@@ -1060,13 +1061,12 @@ function cpd(
     kwargs...,
 ) where {T<:AbstractFloat,N}
     if nonnegative
-        solver_obj = _solver_object(solver, stepsize; kwargs...)
         # Align effective defaults with nncpd() on the nonnegative route.
-        # Explicitly passed non-default values are preserved.
+        stepsize_nn = isnothing(stepsize) ? 0.01 : stepsize
+        solver_obj = _solver_object(solver, stepsize_nn; kwargs...)
         init_nn = _cpd_nonnegative_init(init)
         warm_steps_nn = warm_steps
         geometry_nn = _cpd_nonnegative_geometry(solver_obj, geometry)
-        stepsize_nn = stepsize == 1.0 ? 0.01 : stepsize
         return nncpd(
             A,
             r;
@@ -1090,6 +1090,7 @@ function cpd(
             kwargs...,
         )
     end
+    stepsize_eff = isnothing(stepsize) ? 1.0 : stepsize
     return _cpd_impl(
         A,
         r;
@@ -1100,7 +1101,7 @@ function cpd(
         solver = solver,
         geometry = geometry,
         maxiter = maxiter,
-        stepsize = stepsize,
+        stepsize = stepsize_eff,
         tol = tol,
         gradient_mode = gradient_mode,
         normalization = normalization,

src/core/types.jl

@@ -446,7 +446,7 @@ function LinearAlgebra.cond(R::CPDResult{T}) where {T<:AbstractFloat}
     U = hcat(Us...)
     s = svdvals(U)
     smin = minimum(s)
-    return iszero(smin) ? T(Inf) : inv(smin)
+    return iszero(smin) ? T(Inf) : one(T) / smin
 end
 
 """

src/cpd/core/cp_cost.jl

@@ -69,7 +69,7 @@ end
 end
 
 @inline function _softplus_derivative(x::Real)
-    x >= 0 ? inv(one(x) + exp(-x)) : begin
+    x >= 0 ? one(x) / (one(x) + exp(-x)) : begin
         ex = exp(x)
         ex / (one(x) + ex)
     end

src/manifolds/softplus_metric.jl

@@ -7,7 +7,7 @@
 
 export SoftplusEuclidean, softplus_metric_inverse
 
-@inline _sp_sigmoid(x::Real) = x >= 0 ? inv(one(x) + exp(-x)) : begin
+@inline _sp_sigmoid(x::Real) = x >= 0 ? one(x) / (one(x) + exp(-x)) : begin
     ex = exp(x)
     ex / (one(x) + ex)
 end
@@ -57,8 +57,7 @@ function softplus_metric_diag(M::SoftplusEuclidean, p::AbstractVector)
 end
 
 function softplus_metric_inverse(M::SoftplusEuclidean, p::AbstractVector, X::AbstractVector)
-    g_inv = inv.(softplus_metric_diag(M, p))
-    return g_inv .* X
+    return X ./ softplus_metric_diag(M, p)
 end
 
 pullback_metric_inverse(M::SoftplusEuclidean, p::AbstractVector, X::AbstractVector) =

src/manifolds/squaring_metric.jl

@@ -72,8 +72,7 @@ function pullback_metric_diag(M::SqEuclidean, p::AbstractVector)
 end
 
 function pullback_metric_inverse(M::SqEuclidean, p::AbstractVector, X::AbstractVector)
-    g_inv = 1.0 ./ pullback_metric_diag(M, p)
-    return g_inv .* X
+    return X ./ pullback_metric_diag(M, p)
 end
 
 function ManifoldsBase.inner(M::SqEuclidean, p, X::AbstractVector, Y::AbstractVector)

src/solvers/lbfgs.jl

@@ -3,8 +3,6 @@ export LBFGSSolver
 
 """
     LBFGSSolver(; memory_size=1, cautious_update=true, initial_scale=1.0,
-        nonpositive_curvature_behavior=:ignore,
-        nondescent_direction_behavior=:reinitialize_direction_update,
         linesearch=:wolfe, preconditioner=nothing)
 
 Limited-memory Riemannian BFGS wrapper built on `Manopt.quasi_Newton`.
@@ -13,7 +11,6 @@ struct LBFGSSolver <: AbstractSecondOrderROSolver
     memory_size::Int
     cautious_update::Bool
     initial_scale::Float64
-    nonpositive_curvature_behavior::Symbol
     linesearch::Symbol
     preconditioner::Any
 end
@@ -22,21 +19,22 @@ function LBFGSSolver(;
     memory_size::Int = 1,
     cautious_update::Bool = true,
     initial_scale::Real = 1.0,
-    nonpositive_curvature_behavior::Symbol = :ignore,
     linesearch::Symbol = :wolfe,
     preconditioner = nothing,
 )
     memory_size >= 1 || throw(ArgumentError("memory_size must be >= 1, got $memory_size"))
     initial_scale > 0 ||
         throw(ArgumentError("initial_scale must be > 0, got $initial_scale"))
-    linesearch in (:wolfe, :hagerzhang) || throw(
-        ArgumentError("Unsupported linesearch=$linesearch. Use :wolfe or :hagerzhang."),
+    supported = _lbfgs_supported_linesearches()
+    linesearch in supported || throw(
+        ArgumentError(
+            "Unsupported linesearch=$linesearch. Use one of " * join(supported, ", ") * ".",
+        ),
     )
     return LBFGSSolver(
         memory_size,
         cautious_update,
         Float64(initial_scale),
-        nonpositive_curvature_behavior,
         linesearch,
         preconditioner,
     )
@@ -47,9 +45,20 @@ solver_symbol(::LBFGSSolver) = :lbfgs
 second_order_diagnostics_recorder(::LBFGSSolver) =
     _SolverDiagnosticsRecorder(line_search_enabled = true)
 
+function _lbfgs_supported_linesearches()
+    base = (:wolfe,)
+    return isdefined(Manopt, :HagerZhangLinesearch) ? (base..., :hagerzhang) : base
+end
+
 @inline function _lbfgs_linesearch(kind::Symbol)
-    kind === :hagerzhang && return Manopt.HagerZhangLinesearch()
-    kind === :wolfe && return Manopt.WolfePowellLinesearch()
+    kind === :wolfe && return Manopt.WolfePowellLinesearch(
+        sufficient_curvature = 0.9,
+        stop_when_stepsize_less = 1e-8,
+        stop_decreasing_at_step = 100,
+    )
+    if kind === :hagerzhang && isdefined(Manopt, :HagerZhangLinesearch)
+        return getproperty(Manopt, :HagerZhangLinesearch)()
+    end
     throw(ArgumentError("Unsupported linesearch kind $kind."))
 end
 
@@ -71,7 +80,6 @@ function solve_lbfgs(
     memory_size::Int = 1,
     cautious_update::Bool = true,
     initial_scale::Real = 1.0,
-    nonpositive_curvature_behavior::Symbol = :ignore,
     linesearch::Symbol = :wolfe,
     preconditioner = nothing,
     grad_tol = nothing,
@@ -155,10 +163,8 @@ function solve_lbfgs(
             memory_size = memory_size,
             cautious_update = cautious_update,
             initial_scale = initial_scale,
-            nonpositive_curvature_behavior = nonpositive_curvature_behavior,
             linesearch = linesearch,
             has_preconditioner = !isnothing(preconditioner),
-            uses_nonpositive_curvature_behavior = false,
         ),
     )
 end
@@ -196,7 +202,6 @@ function run_second_order_solver(
         memory_size = solver.memory_size,
         cautious_update = solver.cautious_update,
         initial_scale = solver.initial_scale,
-        nonpositive_curvature_behavior = solver.nonpositive_curvature_behavior,
         linesearch = solver.linesearch,
         preconditioner = solver.preconditioner,
         normalized_objective,

src/solvers/manopt_helpers.jl

@@ -300,7 +300,7 @@ end
 function _relative_solver_functions(model_cost, model_grad, scale::Real)
     scale > 0 || return model_cost, model_grad, false
     scale == one(scale) && return model_cost, model_grad, false
-    inv_scale = inv(scale)
+    inv_scale = one(scale) / scale
     return (
         (M, p) -> model_cost(M, p) * inv_scale,
         (M, p) -> _scale_solver_tangent(model_grad(M, p), inv_scale),

src/solvers/rgd.jl

@@ -20,7 +20,10 @@ function solve_rgd(
     iteration_callbacks = (),
     grad_tol = nothing,
     normalized_objective::Bool = true,
+    armijo_alpha_min::Real = 1e-8,
 )
+    armijo_alpha_min > 0 ||
+        throw(ArgumentError("armijo_alpha_min must be > 0, got $armijo_alpha_min"))
     setup = _prepare_manopt_solver_functions(
         model_cost,
         model_egrad,
@@ -35,15 +38,14 @@ function solve_rgd(
     p0_local = setup.p0
     T = setup.T
     retraction_method = _solver_retraction_method(M, p0_local)
-    stepsize_eff_base = T(stepsize) * setup.objective_scale
-    armijo_alpha_min = T(1e-8) * setup.objective_scale
+    armijo_alpha_min_T = T(armijo_alpha_min)
     tol_g = setup.dual_grad_tol
     dual_stop = StopWhenCostRelChangeAndGradientLess(T(tol), tol_g)
     stopping = _manopt_stopping(
         maxiter,
         setup.grad_stop_tol,
         dual_stop;
-        extra = (StopWhenStepsizeLess(armijo_alpha_min),),
+        extra = (StopWhenStepsizeLess(armijo_alpha_min_T),),
     )
 
     use_squaring_armijo = _contains_sqeuclidean_manifold(M)
@@ -55,13 +57,13 @@ function solve_rgd(
             p0_local,
             setup.solver_grad,
             retraction_method,
-            stepsize_eff_base;
-            alpha_min = armijo_alpha_min,
-        ) : stepsize_eff_base
+            T(stepsize);
+            alpha_min = armijo_alpha_min_T,
+        ) : T(stepsize)
     armijo_contraction = use_squaring_armijo ? T(0.5) : T(0.85)
     armijo_sufficient_decrease = use_squaring_armijo ? T(1e-4) : T(1e-3)
     armijo_stop_decreasing =
-        _armijo_max_decreases(initial_stepsize_eff, armijo_contraction, armijo_alpha_min)
+        _armijo_max_decreases(initial_stepsize_eff, armijo_contraction, armijo_alpha_min_T)
     armijo_max_step = use_strict_sqeuclidean ? initial_stepsize_eff : Inf
     armijo_stop_increasing = use_strict_sqeuclidean ? 0 : 100
     armijo_additional_decrease =
@@ -74,7 +76,7 @@ function solve_rgd(
         initial_stepsize = initial_stepsize_eff,
         contraction_factor = armijo_contraction,
         sufficient_decrease = armijo_sufficient_decrease,
-        stop_when_stepsize_less = armijo_alpha_min,
+        stop_when_stepsize_less = armijo_alpha_min_T,
         stop_when_stepsize_exceeds = armijo_max_step,
         stop_increasing_at_step = armijo_stop_increasing,
         stop_decreasing_at_step = armijo_stop_decreasing,
@@ -122,7 +124,10 @@ function solve_rgd(
         return_stats,
         verbose,
         normalized_objective = setup.uses_relative_objective,
-        solver_info_extra = (initial_stepsize_eff = Float64(initial_stepsize_eff),),
+        solver_info_extra = (
+            initial_stepsize_eff = Float64(initial_stepsize_eff),
+            armijo_alpha_min = Float64(armijo_alpha_min_T),
+        ),
     )
 end
 
@@ -180,7 +185,7 @@ function solve_rgd_fixed(
         setup.solver_grad,
         p0_local;
         retraction_method = retraction_method,
-        stepsize = Manopt.ConstantStepsize(M, T(stepsize) * setup.objective_scale),
+        stepsize = Manopt.ConstantStepsize(M, T(stepsize)),
         stopping_criterion = stopping,
         debug = callbacks.debug_actions,
         count = [:Cost, :Gradient],
@@ -209,15 +214,21 @@ end
 # ========== RGDSolver (AbstractFirstOrderSolver) ==========
 
 """
-    RGDSolver(stepsize=1.0)
+    RGDSolver(stepsize=1.0; armijo_alpha_min=1e-8)
 
 Riemannian gradient descent with Armijo backtracking line search. Call via
 `solve(RGDSolver(...), model; init=:random, gradient_mode=:riemannian)`.
 """
 struct RGDSolver <: AbstractFirstOrderSolver
     stepsize::Float64
+    armijo_alpha_min::Float64
+end
+function RGDSolver(stepsize::Real = 1.0; armijo_alpha_min::Real = 1e-8)
+    stepsize > 0 || throw(ArgumentError("stepsize must be > 0, got $stepsize"))
+    armijo_alpha_min > 0 ||
+        throw(ArgumentError("armijo_alpha_min must be > 0, got $armijo_alpha_min"))
+    return RGDSolver(Float64(stepsize), Float64(armijo_alpha_min))
 end
-RGDSolver() = RGDSolver(1.0)
 
 solver_symbol(::RGDSolver) = :rgd
 first_order_diagnostics_recorder(::RGDSolver) =
@@ -255,6 +266,7 @@ function run_first_order_solver(
         iteration_callbacks,
         grad_tol,
         normalized_objective,
+        armijo_alpha_min = solver.armijo_alpha_min,
     )
 end
 

src/solvers/solve_dispatch.jl

@@ -15,7 +15,8 @@ end
 _solver_object(solver::AbstractSolver, ::Real; kwargs...) = solver
 
 _solver_object(::Val{:als}, ::Real; kwargs...) = ALSSolver()
-_solver_object(::Val{:rgd}, stepsize::Real; kwargs...) = RGDSolver(stepsize)
+_solver_object(::Val{:rgd}, stepsize::Real; kwargs...) =
+    RGDSolver(stepsize; armijo_alpha_min = get(kwargs, :armijo_alpha_min, 1e-8))
 _solver_object(::Val{:rgd_fixed}, stepsize::Real; kwargs...) = RGDFixedSolver(stepsize)
 
 function _solver_object(::Val{:rcg}, ::Real; kwargs...)
@@ -34,11 +35,6 @@ function _solver_object(::Val{:lbfgs}, ::Real; kwargs...)
         memory_size = get(kwargs, :memory_size, 1),
         cautious_update = get(kwargs, :cautious_update, true),
         initial_scale = get(kwargs, :initial_scale, 1.0),
-        nonpositive_curvature_behavior = get(
-            kwargs,
-            :nonpositive_curvature_behavior,
-            :ignore,
-        ),
         linesearch = get(kwargs, :linesearch, :wolfe),
         preconditioner = get(kwargs, :preconditioner, nothing),
     )