From 91faa8ae55bc313e45debb382a09ea280dc4ed99 Mon Sep 17 00:00:00 2001 From: PBrdng Date: Thu, 18 Jun 2026 10:39:56 +0200 Subject: [PATCH 1/6] bug fix --- src/solvers/lbfgs.jl | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/src/solvers/lbfgs.jl b/src/solvers/lbfgs.jl index 1d4c17d..ffbbafe 100644 --- a/src/solvers/lbfgs.jl +++ b/src/solvers/lbfgs.jl @@ -29,8 +29,11 @@ function LBFGSSolver(; 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, @@ -47,9 +50,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 From 08c4a80c295e9907878fd35b3ca7c2b47dff3d72 Mon Sep 17 00:00:00 2001 From: PBrdng Date: Thu, 18 Jun 2026 13:25:39 +0200 Subject: [PATCH 2/6] make stepsize in rgd consistent --- src/solvers/rgd.jl | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/solvers/rgd.jl b/src/solvers/rgd.jl index e7517b5..64a0743 100644 --- a/src/solvers/rgd.jl +++ b/src/solvers/rgd.jl @@ -35,8 +35,7 @@ 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(1e-8) tol_g = setup.dual_grad_tol dual_stop = StopWhenCostRelChangeAndGradientLess(T(tol), tol_g) stopping = _manopt_stopping( @@ -55,9 +54,9 @@ function solve_rgd( p0_local, setup.solver_grad, retraction_method, - stepsize_eff_base; + T(stepsize); alpha_min = armijo_alpha_min, - ) : stepsize_eff_base + ) : 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 = @@ -180,7 +179,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], From 93218fb40533ab2db7d2e0fae17e879e36fc50a4 Mon Sep 17 00:00:00 2001 From: PBrdng Date: Thu, 18 Jun 2026 13:29:25 +0200 Subject: [PATCH 3/6] make armijo_alpha_min a keyword --- src/solvers/rgd.jl | 31 +++++++++++++++++++++++-------- src/solvers/solve_dispatch.jl | 3 ++- 2 files changed, 25 insertions(+), 9 deletions(-) diff --git a/src/solvers/rgd.jl b/src/solvers/rgd.jl index 64a0743..6b5c199 100644 --- a/src/solvers/rgd.jl +++ b/src/solvers/rgd.jl @@ -20,7 +20,11 @@ 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,14 +39,14 @@ function solve_rgd( p0_local = setup.p0 T = setup.T retraction_method = _solver_retraction_method(M, p0_local) - armijo_alpha_min = T(1e-8) + 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,12 +59,12 @@ function solve_rgd( setup.solver_grad, retraction_method, T(stepsize); - alpha_min = armijo_alpha_min, + 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 = @@ -73,7 +77,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, @@ -121,7 +125,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 @@ -208,15 +215,22 @@ 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) = @@ -254,6 +268,7 @@ function run_first_order_solver( iteration_callbacks, grad_tol, normalized_objective, + armijo_alpha_min = solver.armijo_alpha_min, ) end diff --git a/src/solvers/solve_dispatch.jl b/src/solvers/solve_dispatch.jl index 9bdc801..728092f 100644 --- a/src/solvers/solve_dispatch.jl +++ b/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...) From 313d405f9b00036e30ab4ca07298a8735e298d55 Mon Sep 17 00:00:00 2001 From: PBrdng Date: Thu, 18 Jun 2026 13:45:43 +0200 Subject: [PATCH 4/6] fix more inconsistencies --- src/api/cpd.jl | 15 ++++++++------- src/solvers/lbfgs.jl | 9 --------- src/solvers/solve_dispatch.jl | 5 ----- 3 files changed, 8 insertions(+), 21 deletions(-) diff --git a/src/api/cpd.jl b/src/api/cpd.jl index e03ffab..b7b4ecc 100644 --- a/src/api/cpd.jl +++ b/src/api/cpd.jl @@ -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, diff --git a/src/solvers/lbfgs.jl b/src/solvers/lbfgs.jl index ffbbafe..297b8db 100644 --- a/src/solvers/lbfgs.jl +++ b/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,7 +19,6 @@ 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, ) @@ -39,7 +35,6 @@ function LBFGSSolver(; memory_size, cautious_update, Float64(initial_scale), - nonpositive_curvature_behavior, linesearch, preconditioner, ) @@ -85,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, @@ -169,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 @@ -210,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, diff --git a/src/solvers/solve_dispatch.jl b/src/solvers/solve_dispatch.jl index 728092f..4b9bd48 100644 --- a/src/solvers/solve_dispatch.jl +++ b/src/solvers/solve_dispatch.jl @@ -35,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), ) From 161dc7077bb8d235b4c84f89cb2e74ead452f0cd Mon Sep 17 00:00:00 2001 From: Se Eun Choi Date: Thu, 18 Jun 2026 18:27:15 +0200 Subject: [PATCH 5/6] minor fixes: replace inv with one/t for type stability. --- src/api/cpd.jl | 2 +- src/core/types.jl | 2 +- src/cpd/core/cp_cost.jl | 2 +- src/manifolds/softplus_metric.jl | 5 ++--- src/manifolds/squaring_metric.jl | 3 +-- src/solvers/manopt_helpers.jl | 2 +- 6 files changed, 7 insertions(+), 9 deletions(-) diff --git a/src/api/cpd.jl b/src/api/cpd.jl index b7b4ecc..3ad4058 100644 --- a/src/api/cpd.jl +++ b/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 diff --git a/src/core/types.jl b/src/core/types.jl index d39ee6a..9549dcc 100644 --- a/src/core/types.jl +++ b/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 """ diff --git a/src/cpd/core/cp_cost.jl b/src/cpd/core/cp_cost.jl index c99547b..a10ef67 100644 --- a/src/cpd/core/cp_cost.jl +++ b/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 diff --git a/src/manifolds/softplus_metric.jl b/src/manifolds/softplus_metric.jl index 6edfb39..e6c4010 100644 --- a/src/manifolds/softplus_metric.jl +++ b/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) = diff --git a/src/manifolds/squaring_metric.jl b/src/manifolds/squaring_metric.jl index 696d1c3..3a16497 100644 --- a/src/manifolds/squaring_metric.jl +++ b/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) diff --git a/src/solvers/manopt_helpers.jl b/src/solvers/manopt_helpers.jl index 67b8ac8..e6bc167 100644 --- a/src/solvers/manopt_helpers.jl +++ b/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), From 2def375fff65cae1f6ca49c10409165b7ccae853 Mon Sep 17 00:00:00 2001 From: Se Eun Choi Date: Thu, 18 Jun 2026 17:37:37 +0000 Subject: [PATCH 6/6] apply JuliaFormatter --- src/solvers/rgd.jl | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/solvers/rgd.jl b/src/solvers/rgd.jl index 6b5c199..cd45cbd 100644 --- a/src/solvers/rgd.jl +++ b/src/solvers/rgd.jl @@ -22,9 +22,8 @@ function solve_rgd( 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"), - ) + 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, @@ -226,9 +225,8 @@ struct RGDSolver <: AbstractFirstOrderSolver 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"), - ) + 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