From 3dc0a97e7520b1b1577e708ee80cb6c4fc1f292f Mon Sep 17 00:00:00 2001
From: PBrdng <pbrdng@gmail.com>
Date: Wed, 17 Jun 2026 15:54:36 +0200
Subject: [PATCH 01/18] add lbfgs

---
 src/api/cpd.jl | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/src/api/cpd.jl b/src/api/cpd.jl
index 12f5094..d239a71 100644
--- a/src/api/cpd.jl
+++ b/src/api/cpd.jl
@@ -579,13 +579,14 @@ end
 function _validate_cpd_solver_supported(solver::AbstractSolver)
     throw(
         ArgumentError(
-            "Unsupported CPD solver $(typeof(solver)). Use :als, :rgd, :rgd_fixed, or :rcg.",
+            "Unsupported CPD solver $(typeof(solver)). Use :als, :rgd, :rgd_fixed, :rcg, or :lbfgs.",
         ),
     )
 end
 
-_validate_cpd_solver_supported(::Union{ALSSolver,RGDSolver,RGDFixedSolver,RCGSolver}) =
-    nothing
+_validate_cpd_solver_supported(
+    ::Union{ALSSolver,RGDSolver,RGDFixedSolver,RCGSolver,LBFGSSolver},
+) = nothing
 
 function _validate_cpd_solver_options(
     solver::AbstractSolver,
@@ -753,8 +754,9 @@ function _cpd_manifold_grad_tol(
     tol::Real,
 )
     inner = cpd_model(model)
-    inner.nonnegative || return nothing
-    return tol
+    inner.nonnegative && return tol
+    T = eltype(tensor(model))
+    return sqrt(T(tol)) * sqrt(sum(abs2, tensor(model)))
 end
 
 function _cpd_point_rel_error(model, p)
@@ -992,6 +994,7 @@ If `r` is omitted, uses the smallest tensor mode as a heuristic rank.
     - `rgd` (default): Riemannian gradient descent
     - `rgd_fixed`: Riemannian gradient descent with fixed step size
     - `rcg`: Riemannian conjugate gradient
+    - `lbfgs`: Limited-memory Riemannian quasi-Newton
     - `als`: Alternating Least Squares
 
 ## Extended Options

From 45e5138e90c9192e4c46c2787d6ee59a7e4ab6ba Mon Sep 17 00:00:00 2001
From: PBrdng <pbrdng@gmail.com>
Date: Wed, 17 Jun 2026 16:27:13 +0200
Subject: [PATCH 02/18] change the objective function to relative norm

---
 src/api/cpd.jl                |   6 +-
 src/solvers/abstract.jl       |   4 ++
 src/solvers/lbfgs.jl          |  26 ++++++---
 src/solvers/rcg.jl            |  25 ++++++--
 src/solvers/rgd.jl            | 104 +++++++++++++++++++++++++++-------
 src/solvers/solve_dispatch.jl |   2 +
 test/basic_tests.jl           |  26 +++++++++
 7 files changed, 154 insertions(+), 39 deletions(-)

diff --git a/src/api/cpd.jl b/src/api/cpd.jl
index d239a71..e2a8de6 100644
--- a/src/api/cpd.jl
+++ b/src/api/cpd.jl
@@ -753,10 +753,7 @@ function _cpd_manifold_grad_tol(
     solver::Union{RGDSolver,RGDFixedSolver,RCGSolver,LBFGSSolver},
     tol::Real,
 )
-    inner = cpd_model(model)
-    inner.nonnegative && return tol
-    T = eltype(tensor(model))
-    return sqrt(T(tol)) * sqrt(sum(abs2, tensor(model)))
+    return tol
 end
 
 function _cpd_point_rel_error(model, p)
@@ -865,6 +862,7 @@ function _run_cpd_solver(
         refinement_verbose = verbose,
         vector_transport_method,
         grad_tol = _cpd_manifold_grad_tol(model, solver, tol),
+        normalized_objective = solver isa Union{RGDSolver,RGDFixedSolver,RCGSolver,LBFGSSolver},
         iteration_callbacks,
         kwargs...,
     )
diff --git a/src/solvers/abstract.jl b/src/solvers/abstract.jl
index d7a2e20..dac468d 100644
--- a/src/solvers/abstract.jl
+++ b/src/solvers/abstract.jl
@@ -225,6 +225,7 @@ function solve(
     return_stats::Bool = false,
     vector_transport_method::Union{ManifoldsBase.AbstractVectorTransportMethod,Nothing} = nothing,
     grad_tol = nothing,
+    normalized_objective::Bool = false,
     iteration_callbacks = (),
 ) where {T<:AbstractFloat}
     setup = _prepare_solver_problem(model; init, p0, gradient_mode, verbose)
@@ -247,6 +248,7 @@ function solve(
         return_stats,
         vector_transport_method,
         grad_tol,
+        normalized_objective,
         post_step_callback,
         diagnostics_recorder,
         iteration_callbacks,
@@ -266,6 +268,7 @@ function solve(
     return_stats::Bool = false,
     vector_transport_method::Union{ManifoldsBase.AbstractVectorTransportMethod,Nothing} = nothing,
     grad_tol = nothing,
+    normalized_objective::Bool = false,
     iteration_callbacks = (),
 ) where {T<:AbstractFloat}
     setup = _prepare_solver_problem(model; init, p0, gradient_mode)
@@ -288,6 +291,7 @@ function solve(
         return_stats,
         vector_transport_method,
         grad_tol,
+        normalized_objective,
         post_step_callback,
         diagnostics_recorder,
         iteration_callbacks,
diff --git a/src/solvers/lbfgs.jl b/src/solvers/lbfgs.jl
index 992a759..fa0009c 100644
--- a/src/solvers/lbfgs.jl
+++ b/src/solvers/lbfgs.jl
@@ -75,21 +75,28 @@ function solve_lbfgs(
     linesearch::Symbol = :wolfe,
     preconditioner = nothing,
     grad_tol = nothing,
+    normalized_objective::Bool = false,
 )
     p0_local = _solver_point(M, p0)
     T = _scalar_eltype(p0_local)
     model_grad_raw = isnothing(model_grad) ? grad(model_egrad) : model_grad
     model_grad_local = _layout_adapt_gradient(model_grad_raw)
+    objective_scale =
+        normalized_objective && !isnothing(normA2) && normA2 > 0 ? T(normA2) : one(T)
+    solver_cost, solver_grad, uses_relative_objective =
+        _relative_solver_functions(model_cost, model_grad_local, objective_scale)
     retraction_method = _solver_retraction_method(M, p0_local)
     transport =
         isnothing(vector_transport_method) ?
         _default_vector_transport_method(M, p0_local, retraction_method) :
         vector_transport_method
+    grad_stop_tol = isnothing(grad_tol) ? T(tol) : T(grad_tol)
     tol_g = _dual_stop_grad_tol(T, tol, grad_tol)
+    tol_g_raw = uses_relative_objective ? tol_g * objective_scale : tol_g
     dual_stop = StopWhenCostRelChangeAndGradientLess(T(tol), tol_g)
     stopping = StopWhenAny(
         StopAfterIteration(maxiter),
-        StopWhenGradientNormLess(T(tol)),
+        StopWhenGradientNormLess(grad_stop_tol),
         dual_stop,
     )
     progress =
@@ -106,16 +113,17 @@ function solve_lbfgs(
         _solver_diagnostics_callback(diagnostics_recorder)
     progress_callback = _solver_progress_callback(
         progress,
-        model_cost,
-        model_grad_local,
+        solver_cost,
+        solver_grad,
         M;
+        normA2 = uses_relative_objective ? normA2 : nothing,
         diagnostics_recorder,
     )
 
     state = Manopt.quasi_Newton(
         M,
-        model_cost,
-        model_grad_local,
+        solver_cost,
+        solver_grad,
         p0_local;
         cautious_update = cautious_update,
         direction_update = Manopt.InverseBFGS(),
@@ -175,8 +183,9 @@ function solve_lbfgs(
         tol_T = T(tol),
         maxiter,
         solver = :lbfgs,
-        tiny_grad_tol = tol_g,
+        tiny_grad_tol = tol_g_raw,
         solver_info,
+        use_state_gradient = !uses_relative_objective,
     ) :
            _solver_stats(
         model_cost,
@@ -188,8 +197,9 @@ function solve_lbfgs(
         tol_T = T(tol),
         maxiter,
         solver = :lbfgs,
-        tiny_grad_tol = tol_g,
+        tiny_grad_tol = tol_g_raw,
         solver_info,
+        use_state_gradient = !uses_relative_objective,
     )
 end
 
@@ -205,6 +215,7 @@ function run_second_order_solver(
     diagnostics_recorder,
     iteration_callbacks,
     grad_tol = nothing,
+    normalized_objective::Bool = false,
 )
     return solve_lbfgs(
         setup.model_cost,
@@ -228,5 +239,6 @@ function run_second_order_solver(
         nonpositive_curvature_behavior = solver.nonpositive_curvature_behavior,
         linesearch = solver.linesearch,
         preconditioner = solver.preconditioner,
+        normalized_objective,
     )
 end
diff --git a/src/solvers/rcg.jl b/src/solvers/rcg.jl
index be3a705..a1ac2d9 100644
--- a/src/solvers/rcg.jl
+++ b/src/solvers/rcg.jl
@@ -106,23 +106,29 @@ function solve_rcg(
     diagnostics_recorder = nothing,
     iteration_callbacks = (),
     grad_tol = nothing,
+    normalized_objective::Bool = false,
 )
     p0_local = _solver_point(M, p0)
     T = _scalar_eltype(p0_local)
     model_grad_raw = isnothing(model_grad) ? grad(model_egrad) : model_grad
     model_grad_local = _layout_adapt_gradient(model_grad_raw)
+    objective_scale =
+        normalized_objective && !isnothing(normA2) && normA2 > 0 ? T(normA2) : one(T)
+    solver_cost, solver_grad, uses_relative_objective =
+        _relative_solver_functions(model_cost, model_grad_local, objective_scale)
     retraction_method = _solver_retraction_method(M, p0_local)
-    stopping = StopWhenAny(StopAfterIteration(maxiter), StopWhenGradientNormLess(T(tol)))
     transport =
         isnothing(vector_transport_method) ?
         _default_vector_transport_method(M, p0_local, retraction_method) :
         vector_transport_method
+    grad_stop_tol = isnothing(grad_tol) ? T(tol) : T(grad_tol)
     tol_g = _dual_stop_grad_tol(T, tol, grad_tol)
+    tol_g_raw = uses_relative_objective ? tol_g * objective_scale : tol_g
     dual_stop = StopWhenCostRelChangeAndGradientLess(T(tol), tol_g)
 
     stopping = StopWhenAny(
         StopAfterIteration(maxiter),
-        StopWhenGradientNormLess(T(tol)),
+        StopWhenGradientNormLess(grad_stop_tol),
         dual_stop,
     )
     progress =
@@ -134,16 +140,17 @@ function solve_rcg(
         _solver_diagnostics_callback(diagnostics_recorder)
     progress_callback = _solver_progress_callback(
         progress,
-        model_cost,
-        model_grad_local,
+        solver_cost,
+        solver_grad,
         M;
+        normA2 = uses_relative_objective ? normA2 : nothing,
         diagnostics_recorder,
     )
 
     state = conjugate_gradient_descent(
         M,
-        model_cost,
-        model_grad_local,
+        solver_cost,
+        solver_grad,
         p0_local;
         retraction_method = retraction_method,
         vector_transport_method = transport,
@@ -185,7 +192,9 @@ function solve_rcg(
         tol_T = T(tol),
         maxiter,
         solver = :rcg,
+        tiny_grad_tol = tol_g_raw,
         solver_info,
+        use_state_gradient = !uses_relative_objective,
     ) :
            _solver_stats(
         model_cost,
@@ -197,7 +206,9 @@ function solve_rcg(
         tol_T = T(tol),
         maxiter,
         solver = :rcg,
+        tiny_grad_tol = tol_g_raw,
         solver_info,
+        use_state_gradient = !uses_relative_objective,
     )
 end
 
@@ -227,6 +238,7 @@ function run_first_order_solver(
     diagnostics_recorder,
     iteration_callbacks,
     grad_tol = nothing,
+    normalized_objective::Bool = false,
 )
     return solve_rcg(
         setup.model_cost,
@@ -244,5 +256,6 @@ function run_first_order_solver(
         diagnostics_recorder,
         iteration_callbacks,
         grad_tol,
+        normalized_objective,
     )
 end
diff --git a/src/solvers/rgd.jl b/src/solvers/rgd.jl
index 5c755c2..93b1d73 100644
--- a/src/solvers/rgd.jl
+++ b/src/solvers/rgd.jl
@@ -259,6 +259,32 @@ end
     end
 end
 
+@inline _scale_solver_tangent(x::Number, scale::Real) = x * scale
+_scale_solver_tangent(x::AbstractArray, scale::Real) = x .* scale
+_scale_solver_tangent(x::ArrayPartition, scale::Real) =
+    ArrayPartition(map(part -> _scale_solver_tangent(part, scale), x.x)...)
+_scale_solver_tangent(x::Tuple, scale::Real) =
+    map(part -> _scale_solver_tangent(part, scale), x)
+
+function _scale_solver_tangent(x, scale::Real)
+    try
+        return x .* scale
+    catch
+        return scale * x
+    end
+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)
+    return (
+        (M, p) -> model_cost(M, p) * inv_scale,
+        (M, p) -> _scale_solver_tangent(model_grad(M, p), inv_scale),
+        true,
+    )
+end
+
 @inline _solver_has_converged(state) = Manopt.has_converged(state)
 
 
@@ -293,12 +319,13 @@ function _solver_stats(
     solver::Symbol,
     tiny_grad_tol = nothing,
     solver_info = (;),
+    use_state_gradient::Bool = true,
 )
     T = typeof(tol_T)
     final_cost = model_cost(M, p_opt)
     cost_for_error = max(T(0), T(2) * final_cost)
     rel_error = sqrt(cost_for_error)
-    grad_state = _solver_gradient(state)
+    grad_state = use_state_gradient ? _solver_gradient(state) : nothing
     grad_from_state = !isnothing(grad_state)
     grad_final =
         isnothing(grad_state) ? model_grad(M, p_opt) :
@@ -341,12 +368,13 @@ function _solver_stats(
     solver::Symbol,
     tiny_grad_tol = nothing,
     solver_info = (;),
+    use_state_gradient::Bool = true,
 )
     T = typeof(tol_T)
     final_cost = model_cost(M, p_opt)
     cost_for_error = max(T(0), T(2) * final_cost)
     rel_error = _relative_error_frob_sq(cost_for_error, T(normA2))
-    grad_state = _solver_gradient(state)
+    grad_state = use_state_gradient ? _solver_gradient(state) : nothing
     grad_from_state = !isnothing(grad_state)
     grad_final =
         isnothing(grad_state) ? model_grad(M, p_opt) :
@@ -414,16 +442,21 @@ function _solver_progress_callback(
     model_cost,
     model_grad,
     M;
+    normA2 = nothing,
     diagnostics_recorder = nothing,
 )
     progress isa NoMethodProgress && return nothing
+    has_relative_scale = !isnothing(normA2) && normA2 > 0
+    target_norm = has_relative_scale ? sqrt(normA2) : nothing
     return function (problem, state, k)
         k <= 0 && return nothing
         p = get_iterate(state)
         c = model_cost(M, p)
         g = model_grad(M, p)
         gnorm = norm(M, p, g)
-        showvalues = Any[("Iter", k), ("Cost", c), ("Grad norm", gnorm)]
+        c_display = has_relative_scale ? sqrt(max(2 * c, zero(c))) : c
+        gnorm_display = has_relative_scale ? gnorm * target_norm : gnorm
+        showvalues = Any[("Iter", k), ("Cost", c_display), ("Grad norm", gnorm_display)]
         if !isnothing(diagnostics_recorder)
             step = diagnostics_recorder.accepted_stepsize_history
             trials = diagnostics_recorder.line_search_trial_history
@@ -580,19 +613,27 @@ function solve_rgd(
     diagnostics_recorder = nothing,
     iteration_callbacks = (),
     grad_tol = nothing,
+    normalized_objective::Bool = false,
 )
     p0_local = _solver_point(M, p0)
     T = _scalar_eltype(p0_local)
     model_grad_raw = isnothing(model_grad) ? grad(model_egrad) : model_grad
     model_grad_local = _layout_adapt_gradient(model_grad_raw)
+    objective_scale =
+        normalized_objective && !isnothing(normA2) && normA2 > 0 ? T(normA2) : one(T)
+    solver_cost_base, solver_grad, uses_relative_objective =
+        _relative_solver_functions(model_cost, model_grad_local, objective_scale)
     retraction_method = _solver_retraction_method(M, p0_local)
-    armijo_alpha_min = T(1e-8)
+    stepsize_eff_base = T(stepsize) * objective_scale
+    armijo_alpha_min = T(1e-8) * objective_scale
+    grad_stop_tol = isnothing(grad_tol) ? T(tol) : T(grad_tol)
     tol_g = _dual_stop_grad_tol(T, tol, grad_tol)
+    tol_g_raw = uses_relative_objective ? tol_g * objective_scale : tol_g
     dual_stop = StopWhenCostRelChangeAndGradientLess(T(tol), tol_g)
 
     stopping = StopWhenAny(
         StopAfterIteration(maxiter),
-        StopWhenGradientNormLess(T(tol)),
+        StopWhenGradientNormLess(grad_stop_tol),
         StopWhenStepsizeLess(armijo_alpha_min),
         dual_stop,
     )
@@ -604,11 +645,11 @@ function solve_rgd(
         _adaptive_initial_stepsize(
             M,
             p0_local,
-            model_grad_local,
+            solver_grad,
             retraction_method,
-            T(stepsize);
+            stepsize_eff_base;
             alpha_min = armijo_alpha_min,
-        ) : T(stepsize)
+        ) : stepsize_eff_base
     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 =
@@ -617,7 +658,8 @@ function solve_rgd(
     armijo_stop_increasing = use_strict_sqeuclidean ? 0 : 100
     armijo_additional_decrease =
         use_strict_sqeuclidean ? ((M, q) -> _all_finite(q)) : ((M, q) -> true)
-    solver_cost = use_strict_sqeuclidean ? _safe_cost_function(model_cost) : model_cost
+    solver_cost =
+        use_strict_sqeuclidean ? _safe_cost_function(solver_cost_base) : solver_cost_base
     armijo = Manopt.ArmijoLinesearch(
         M;
         retraction_method = retraction_method,
@@ -641,15 +683,16 @@ function solve_rgd(
     progress_callback = _solver_progress_callback(
         progress,
         solver_cost,
-        model_grad_local,
+        solver_grad,
         M;
+        normA2 = uses_relative_objective ? normA2 : nothing,
         diagnostics_recorder,
     )
 
     state = gradient_descent(
         M,
         solver_cost,
-        model_grad_local,
+        solver_grad,
         p0_local;
         retraction_method = retraction_method,
         stepsize = armijo,
@@ -684,7 +727,7 @@ function solve_rgd(
     end
     return isnothing(normA2) ?
            _solver_stats(
-        solver_cost,
+        model_cost,
         model_grad_local,
         M,
         p_opt,
@@ -693,11 +736,12 @@ function solve_rgd(
         tol_T = T(tol),
         maxiter,
         solver = :rgd,
-        tiny_grad_tol = tol_g,
+        tiny_grad_tol = tol_g_raw,
         solver_info,
+        use_state_gradient = !uses_relative_objective,
     ) :
            _solver_stats(
-        solver_cost,
+        model_cost,
         model_grad_local,
         M,
         p_opt,
@@ -706,8 +750,9 @@ function solve_rgd(
         tol_T = T(tol),
         maxiter,
         solver = :rgd,
-        tiny_grad_tol = tol_g,
+        tiny_grad_tol = tol_g_raw,
         solver_info,
+        use_state_gradient = !uses_relative_objective,
     )
 end
 
@@ -727,14 +772,22 @@ function solve_rgd_fixed(
     diagnostics_recorder = nothing,
     iteration_callbacks = (),
     grad_tol = nothing,
+    normalized_objective::Bool = false,
 )
     p0_local = _solver_point(M, p0)
     T = _scalar_eltype(p0_local)
     model_grad_raw = isnothing(model_grad) ? grad(model_egrad) : model_grad
     model_grad_local = _layout_adapt_gradient(model_grad_raw)
+    objective_scale =
+        normalized_objective && !isnothing(normA2) && normA2 > 0 ? T(normA2) : one(T)
+    solver_cost, solver_grad, uses_relative_objective =
+        _relative_solver_functions(model_cost, model_grad_local, objective_scale)
     retraction_method = _solver_retraction_method(M, p0_local)
-    tiny_grad_tol = isnothing(grad_tol) ? T(1e-5) : T(grad_tol)
-    stopping = StopWhenAny(StopAfterIteration(maxiter), StopWhenGradientNormLess(T(tol)))
+    grad_stop_tol = isnothing(grad_tol) ? T(tol) : T(grad_tol)
+    tiny_grad_tol =
+        isnothing(grad_tol) ? T(1e-5) :
+        (uses_relative_objective ? T(grad_tol) * objective_scale : T(grad_tol))
+    stopping = StopWhenAny(StopAfterIteration(maxiter), StopWhenGradientNormLess(grad_stop_tol))
     progress =
         maxiter > 0 ?
         make_rgd_fixed_progress(maxiter; enabled = verbose, phase = :refinement, dt = 0.2) :
@@ -744,18 +797,19 @@ function solve_rgd_fixed(
         _solver_diagnostics_callback(diagnostics_recorder)
     progress_callback = _solver_progress_callback(
         progress,
-        model_cost,
-        model_grad_local,
+        solver_cost,
+        solver_grad,
         M;
+        normA2 = uses_relative_objective ? normA2 : nothing,
         diagnostics_recorder,
     )
     state = gradient_descent(
         M,
-        model_cost,
-        model_grad_local,
+        solver_cost,
+        solver_grad,
         p0_local;
         retraction_method = retraction_method,
-        stepsize = Manopt.ConstantStepsize(M, T(stepsize)),
+        stepsize = Manopt.ConstantStepsize(M, T(stepsize) * objective_scale),
         stopping_criterion = stopping,
         debug = _solver_debug_actions(
             verbose,
@@ -796,6 +850,7 @@ function solve_rgd_fixed(
         solver = :rgd_fixed,
         tiny_grad_tol = tiny_grad_tol,
         solver_info,
+        use_state_gradient = !uses_relative_objective,
     ) :
            _solver_stats(
         model_cost,
@@ -809,6 +864,7 @@ function solve_rgd_fixed(
         solver = :rgd_fixed,
         tiny_grad_tol = tiny_grad_tol,
         solver_info,
+        use_state_gradient = !uses_relative_objective,
     )
 end
 
@@ -841,6 +897,7 @@ function run_first_order_solver(
     diagnostics_recorder,
     iteration_callbacks,
     grad_tol = nothing,
+    normalized_objective::Bool = false,
 )
     return solve_rgd(
         setup.model_cost,
@@ -859,6 +916,7 @@ function run_first_order_solver(
         diagnostics_recorder,
         iteration_callbacks,
         grad_tol,
+        normalized_objective,
     )
 end
 
@@ -894,6 +952,7 @@ function run_first_order_solver(
     diagnostics_recorder,
     iteration_callbacks,
     grad_tol = nothing,
+    normalized_objective::Bool = false,
 )
     return solve_rgd_fixed(
         setup.model_cost,
@@ -911,5 +970,6 @@ function run_first_order_solver(
         diagnostics_recorder,
         iteration_callbacks,
         grad_tol,
+        normalized_objective,
     )
 end
diff --git a/src/solvers/solve_dispatch.jl b/src/solvers/solve_dispatch.jl
index fbb986a..fcce422 100644
--- a/src/solvers/solve_dispatch.jl
+++ b/src/solvers/solve_dispatch.jl
@@ -65,6 +65,7 @@ function _solve_with_solver(
     verbose::Bool,
     vector_transport_method::Union{ManifoldsBase.AbstractVectorTransportMethod,Nothing} = nothing,
     grad_tol = nothing,
+    normalized_objective::Bool = false,
     iteration_callbacks = (),
     kwargs...,
 )
@@ -81,6 +82,7 @@ function _solve_with_solver(
         return_stats = true,
         vector_transport_method,
         grad_tol,
+        normalized_objective,
         iteration_callbacks,
     )
 end
diff --git a/test/basic_tests.jl b/test/basic_tests.jl
index f225ccf..041a458 100644
--- a/test/basic_tests.jl
+++ b/test/basic_tests.jl
@@ -272,6 +272,32 @@ end
     @test cpd_rgd_object isa CPDResult
     @test cpd_rgd_object.solver == :rgd
 
+    cpd_lbfgs_symbol = cpd(
+        A,
+        2;
+        solver = :lbfgs,
+        init = :alswarm,
+        warm_steps = 2,
+        maxiter = 2,
+        verbose = false,
+    )
+    @test cpd_lbfgs_symbol isa CPDResult
+    @test cpd_lbfgs_symbol.solver == :lbfgs
+    @test cpd_lbfgs_symbol.solver_info.memory_size == 1
+
+    cpd_lbfgs_object = cpd(
+        A,
+        2;
+        solver = LBFGSSolver(memory_size = 3),
+        init = :alswarm,
+        warm_steps = 2,
+        maxiter = 2,
+        verbose = false,
+    )
+    @test cpd_lbfgs_object isa CPDResult
+    @test cpd_lbfgs_object.solver == :lbfgs
+    @test cpd_lbfgs_object.solver_info.memory_size == 3
+
     cpd_als_object =
         cpd(A, 2; solver = ALSSolver(), init = :tucker, maxiter = 1, verbose = false)
     @test cpd_als_object isa CPDResult

From 36622d0722abf7f04f4ed25e89793ae21d6b766a Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Wed, 17 Jun 2026 15:02:27 +0000
Subject: [PATCH 03/18] fix: apply JuliaFormatter output for cpd and rgd files

---
 src/api/cpd.jl     | 3 ++-
 src/solvers/rgd.jl | 3 ++-
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/api/cpd.jl b/src/api/cpd.jl
index e2a8de6..2d19de5 100644
--- a/src/api/cpd.jl
+++ b/src/api/cpd.jl
@@ -862,7 +862,8 @@ function _run_cpd_solver(
         refinement_verbose = verbose,
         vector_transport_method,
         grad_tol = _cpd_manifold_grad_tol(model, solver, tol),
-        normalized_objective = solver isa Union{RGDSolver,RGDFixedSolver,RCGSolver,LBFGSSolver},
+        normalized_objective = solver isa
+                               Union{RGDSolver,RGDFixedSolver,RCGSolver,LBFGSSolver},
         iteration_callbacks,
         kwargs...,
     )
diff --git a/src/solvers/rgd.jl b/src/solvers/rgd.jl
index 93b1d73..af180e9 100644
--- a/src/solvers/rgd.jl
+++ b/src/solvers/rgd.jl
@@ -787,7 +787,8 @@ function solve_rgd_fixed(
     tiny_grad_tol =
         isnothing(grad_tol) ? T(1e-5) :
         (uses_relative_objective ? T(grad_tol) * objective_scale : T(grad_tol))
-    stopping = StopWhenAny(StopAfterIteration(maxiter), StopWhenGradientNormLess(grad_stop_tol))
+    stopping =
+        StopWhenAny(StopAfterIteration(maxiter), StopWhenGradientNormLess(grad_stop_tol))
     progress =
         maxiter > 0 ?
         make_rgd_fixed_progress(maxiter; enabled = verbose, phase = :refinement, dt = 0.2) :

From 3e17731a685f2b3362521744fa2aec7d6650d9e9 Mon Sep 17 00:00:00 2001
From: PBrdng <pbrdng@gmail.com>
Date: Wed, 17 Jun 2026 20:36:40 +0200
Subject: [PATCH 04/18] make progress meter render correctly both phases

---
 src/core/progress.jl | 59 +++++++++++++++++++++++++++++++++++---------
 1 file changed, 48 insertions(+), 11 deletions(-)

diff --git a/src/core/progress.jl b/src/core/progress.jl
index 249deb4..30c514a 100644
--- a/src/core/progress.jl
+++ b/src/core/progress.jl
@@ -177,6 +177,32 @@ end
     return nothing
 end
 
+@inline function _force_visible_phase_finish(tracker::PhaseProgress, progress)
+    return tracker.phase == :refinement &&
+           _was_rendered(tracker.initialization) &&
+           !_was_rendered(progress)
+end
+
+function _render_unrendered_completion!(meter, showvalues)
+    # ProgressMeter does not render a meter that reaches 100% before its first
+    # visible update. Give it one display-only step before completion.
+    PM.update!(
+        meter,
+        meter.n;
+        showvalues,
+        force = true,
+        max_steps = meter.n + 1,
+    )
+    return nothing
+end
+
+function _force_visible_unrendered_progress!(progress, meter, showvalues)
+    _was_rendered(progress) && return nothing
+    _render_unrendered_completion!(meter, showvalues)
+    _mark_rendered!(progress)
+    return nothing
+end
+
 update_progress!(::NoMethodProgress, args...; kwargs...) = nothing
 
 function update_progress!(
@@ -192,14 +218,17 @@ function update_progress!(
         set_phase!(tracker, progress.phase)
     end
     t = time()
-    if force || current >= meter.n || t > meter.tlast + meter.dt
+    renders_by_time = t > meter.tlast + meter.dt
+    if force || current >= meter.n || renders_by_time
         showvalues_with_method = if isnothing(showvalues)
             Any[("Method", _method_name(progress))]
         else
             Any[("Method", _method_name(progress)); showvalues]
         end
-        PM.update!(meter, current; showvalues = showvalues_with_method)
-        _mark_rendered!(progress)
+        PM.update!(meter, current; showvalues = showvalues_with_method, force)
+        if force || current < meter.n || _was_rendered(progress)
+            _mark_rendered!(progress)
+        end
     end
     return nothing
 end
@@ -222,17 +251,14 @@ function finish_progress!(
         Any[("Method", _method_name(progress)); showvalues]
     end
 
-    force_refinement_finish =
-        tracker.phase == :refinement &&
-        _was_rendered(tracker.initialization) &&
-        !_was_rendered(progress)
-
-    if force_refinement_finish
-        PM.update!(meter, meter.n; showvalues = showvalues_with_method)
-        _mark_rendered!(progress)
+    if _force_visible_phase_finish(tracker, progress)
+        _force_visible_unrendered_progress!(progress, meter, showvalues_with_method)
+        PM.finish!(meter; showvalues = showvalues_with_method)
         return nothing
     end
 
+    _was_rendered(progress) || return nothing
+
     PM.finish!(meter; showvalues = showvalues_with_method)
     _mark_rendered!(progress)
     return nothing
@@ -251,12 +277,23 @@ function finish_progress!(
         if active_progress(tracker) === progress
             return finish_progress!(tracker; current, showvalues)
         end
+        if _force_visible_phase_finish(tracker, progress)
+            showvalues_with_method = if isnothing(showvalues)
+                Any[("Method", _method_name(progress))]
+            else
+                Any[("Method", _method_name(progress)); showvalues]
+            end
+            _force_visible_unrendered_progress!(progress, meter, showvalues_with_method)
+            PM.finish!(meter; showvalues = showvalues_with_method)
+            return nothing
+        end
     end
     showvalues_with_method = if isnothing(showvalues)
         Any[("Method", _method_name(progress))]
     else
         Any[("Method", _method_name(progress)); showvalues]
     end
+    _was_rendered(progress) || return nothing
     PM.finish!(meter; showvalues = showvalues_with_method)
     _mark_rendered!(progress)
     return nothing

From a1cffff4f08ceda6e1059f273ab77c3af70d0d1c Mon Sep 17 00:00:00 2001
From: PBrdng <pbrdng@gmail.com>
Date: Wed, 17 Jun 2026 20:39:44 +0200
Subject: [PATCH 05/18] export manopt helpers to a separate file

---
 src/backend.jl                |   1 +
 src/solvers/manopt_helpers.jl | 594 ++++++++++++++++++++++++++++++++++
 src/solvers/rgd.jl            | 594 ----------------------------------
 3 files changed, 595 insertions(+), 594 deletions(-)
 create mode 100644 src/solvers/manopt_helpers.jl

diff --git a/src/backend.jl b/src/backend.jl
index ab39e34..5a96538 100644
--- a/src/backend.jl
+++ b/src/backend.jl
@@ -3,6 +3,7 @@ include("solvers/nncp_updates.jl")
 include("solvers/cp_als.jl")
 include("solvers/btd_als.jl")
 include("solvers/rals.jl")
+include("solvers/manopt_helpers.jl")
 include("solvers/rgd.jl")
 include("solvers/btd_tsd.jl")
 include("solvers/rcg.jl")
diff --git a/src/solvers/manopt_helpers.jl b/src/solvers/manopt_helpers.jl
new file mode 100644
index 0000000..a107e27
--- /dev/null
+++ b/src/solvers/manopt_helpers.jl
@@ -0,0 +1,594 @@
+
+struct _SolverDebugSink <: IO end
+Base.isopen(::_SolverDebugSink) = true
+Base.write(::_SolverDebugSink, ::UInt8) = 1
+Base.write(::_SolverDebugSink, s::Union{String,SubString{String}}) = sizeof(s)
+Base.unsafe_write(::_SolverDebugSink, ::Ptr{UInt8}, n::UInt) = Int(n)
+
+const _SOLVER_DEBUG_SINK = _SolverDebugSink()
+
+mutable struct StopWhenCostRelChangeAndGradientLess{T<:Real} <: Manopt.StoppingCriterion
+    tol_cost::T
+    tol_grad::T
+    prev_cost::T
+    last_cost_rel_change::T
+    last_grad_norm::T
+    at_iteration::Int
+end
+
+function StopWhenCostRelChangeAndGradientLess(tol_cost::T, tol_grad::T) where {T<:Real}
+    return StopWhenCostRelChangeAndGradientLess{T}(
+        tol_cost,
+        tol_grad,
+        T(Inf),
+        T(Inf),
+        T(Inf),
+        -1,
+    )
+end
+
+function (c::StopWhenCostRelChangeAndGradientLess)(problem, state, i)
+    if i == 0
+        c.prev_cost = Manopt.get_cost(problem, Manopt.get_iterate(state))
+        c.last_cost_rel_change = oftype(c.tol_cost, Inf)
+        c.last_grad_norm = oftype(c.tol_grad, Inf)
+        c.at_iteration = -1
+        return false
+    end
+    M = Manopt.get_manifold(problem)
+    p = Manopt.get_iterate(state)
+    cost_val = Manopt.get_cost(problem, p)
+    grad_val = Manopt.get_gradient(problem, p)
+    grad_norm = norm(M, p, grad_val)
+    rel_change = abs(c.prev_cost - cost_val) / max(abs(c.prev_cost), one(cost_val))
+    c.prev_cost = cost_val
+    c.last_cost_rel_change = rel_change
+    c.last_grad_norm = grad_norm
+    if rel_change < c.tol_cost && grad_norm < c.tol_grad
+        c.at_iteration = i
+        return true
+    end
+    return false
+end
+
+function Manopt.get_reason(c::StopWhenCostRelChangeAndGradientLess)
+    if c.at_iteration >= 0
+        return "At iteration $(c.at_iteration) the relative cost change ($(c.last_cost_rel_change)) " *
+               "is below $(c.tol_cost) and the gradient norm ($(c.last_grad_norm)) " *
+               "is below $(c.tol_grad).\n"
+    end
+    return ""
+end
+
+function Manopt.status_summary(c::StopWhenCostRelChangeAndGradientLess)
+    has_stopped = c.at_iteration >= 0
+    status = has_stopped ? "reached" : "not reached"
+    return "cost rel change < $(c.tol_cost) and |grad f| < $(c.tol_grad): $status"
+end
+
+Manopt.indicates_convergence(::StopWhenCostRelChangeAndGradientLess) = true
+
+function Base.show(io::IO, c::StopWhenCostRelChangeAndGradientLess)
+    return print(
+        io,
+        "StopWhenCostRelChangeAndGradientLess($(c.tol_cost), $(c.tol_grad))\n    $(Manopt.status_summary(c))",
+    )
+end
+
+
+function _tk_get_solver_result(state)
+    try
+        return Manopt.get_solver_result(state)
+    catch
+    end
+    while hasproperty(state, :state)
+        state = state.state
+    end
+    for key in (:p, :x, :point)
+        hasproperty(state, key) && return getproperty(state, key)
+    end
+    throw(
+        ArgumentError(
+            "Cannot extract point from state $(typeof(state)). Properties: $(propertynames(state))",
+        ),
+    )
+end
+
+
+@inline _align_layout_like_point(p, x) =
+    hasproperty(p, :x) ?
+    (hasproperty(x, :x) ? x : (x isa Tuple ? ArrayPartition(x...) : x)) :
+    (hasproperty(x, :x) ? Tuple(getproperty(x, :x)) : x)
+
+
+function _to_array_partition(x)
+    if x isa ArrayPartition
+        return ArrayPartition(map(_to_array_partition, x.x)...)
+    elseif hasproperty(x, :x)
+        return ArrayPartition(map(_to_array_partition, getproperty(x, :x))...)
+    elseif x isa Tuple
+        return ArrayPartition(map(_to_array_partition, x)...)
+    end
+    return x
+end
+
+
+function _solver_point(M, p0)
+    M2 = _unwrap_solver_manifold(M)
+    return M2 isa ProductManifold ? _to_array_partition(p0) : p0
+end
+
+
+function _contains_sqeuclidean_manifold(M)
+    M2 = _unwrap_solver_manifold(M)
+    if M2 isa SqEuclidean || M2 isa SoftplusEuclidean
+        return true
+    elseif M2 isa ProductManifold
+        return any(_contains_sqeuclidean_manifold, M2.manifolds)
+    elseif hasproperty(M2, :native) && (
+        getproperty(M2, :native) isa SqEuclidean ||
+        getproperty(M2, :native) isa SoftplusEuclidean
+    )
+        return true
+    end
+    return false
+end
+
+function _contains_strict_sqeuclidean_manifold(M)
+    M2 = _unwrap_solver_manifold(M)
+    if M2 isa SqEuclidean
+        return true
+    elseif M2 isa ProductManifold
+        return any(_contains_strict_sqeuclidean_manifold, M2.manifolds)
+    elseif hasproperty(M2, :native) && (getproperty(M2, :native) isa SqEuclidean)
+        return true
+    end
+    return false
+end
+
+
+function _armijo_max_decreases(initial_stepsize::Real, contraction::Real, alpha_min::Real)
+    initial_stepsize <= alpha_min && return 0
+    (contraction <= 0 || contraction >= 1) && return 1000
+    n = floor(Int, log(alpha_min / initial_stepsize) / log(contraction))
+    return max(n, 0)
+end
+
+
+function _adaptive_initial_stepsize(
+    M,
+    p0,
+    model_grad,
+    retraction_method,
+    base_stepsize::T;
+    alpha_min::T,
+    scale_c::T = one(T),
+    clamp_low_factor::T = T(0.1),
+    clamp_high_factor::T = T(10),
+    delta_scale::T = T(1e-3),
+) where {T<:AbstractFloat}
+    g0 = model_grad(M, p0)
+    d = -copy(g0)
+    dnorm = norm(M, p0, d)
+    (!isfinite(dnorm) || dnorm <= sqrt(eps(T))) && return base_stepsize
+    δ = delta_scale / max(dnorm, one(T))
+    q = try
+        retract(M, p0, δ .* d, retraction_method)
+    catch
+        return base_stepsize
+    end
+    _all_finite(q) || return base_stepsize
+    gq = model_grad(M, q)
+    _all_finite(gq) || return base_stepsize
+    κ_num = inner(M, p0, gq .- g0, d)
+    κ_den = δ * dnorm^2
+    (!isfinite(κ_num) || !isfinite(κ_den) || κ_den <= eps(T)) && return base_stepsize
+    κ = max(κ_num / κ_den, eps(T))
+    α_raw = scale_c / κ
+    α_low = max(alpha_min, clamp_low_factor * base_stepsize)
+    α_high = clamp_high_factor * base_stepsize
+    α = clamp(α_raw, α_low, α_high)
+    if !isfinite(α)
+        return base_stepsize
+    end
+    return α
+end
+
+function _all_finite(x)
+    if x isa Number
+        return isfinite(x)
+    elseif hasproperty(x, :x)
+        return all(_all_finite, getproperty(x, :x))
+    elseif x isa AbstractArray
+        return all(isfinite, x)
+    elseif x isa Tuple
+        return all(_all_finite, x)
+    elseif x isa Manifolds.TuckerPoint
+        return _all_finite(x.hosvd.core) && all(_all_finite, x.hosvd.U)
+    elseif x isa Manifolds.TuckerTangentVector
+        return _all_finite(getproperty(x, :Ċ)) && all(_all_finite, getproperty(x, :U̇))
+    end
+    try
+        return all(isfinite, x)
+    catch
+        return false
+    end
+end
+
+
+function _safe_cost_function(model_cost)
+    return function (M, p)
+        _all_finite(p) || return Inf
+        c = model_cost(M, p)
+        return isfinite(c) ? c : Inf
+    end
+end
+
+
+function _layout_adapt_gradient(model_grad)
+    return function (M, p)
+        g = model_grad(M, p)
+        return _align_layout_like_point(p, g)
+    end
+end
+
+function _scalar_eltype(p)
+    if hasproperty(p, :x) || p isa AbstractVector || p isa Tuple
+        parts = point_parts(p)
+        isempty(parts) && return Float64
+        return _scalar_eltype(first(parts))
+    elseif p isa Manifolds.TuckerPoint
+        return eltype(p.hosvd.core)
+    elseif p isa Manifolds.TuckerTangentVector
+        return eltype(getproperty(p, :Ċ))
+    elseif p isa Real
+        return typeof(p)
+    else
+        return eltype(p)
+    end
+end
+
+@inline function _solver_gradient(state)
+    try
+        return Manopt.get_gradient(state)
+    catch
+        return nothing
+    end
+end
+
+@inline _scale_solver_tangent(x::Number, scale::Real) = x * scale
+_scale_solver_tangent(x::AbstractArray, scale::Real) = x .* scale
+_scale_solver_tangent(x::ArrayPartition, scale::Real) =
+    ArrayPartition(map(part -> _scale_solver_tangent(part, scale), x.x)...)
+_scale_solver_tangent(x::Tuple, scale::Real) =
+    map(part -> _scale_solver_tangent(part, scale), x)
+
+function _scale_solver_tangent(x, scale::Real)
+    try
+        return x .* scale
+    catch
+        return scale * x
+    end
+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)
+    return (
+        (M, p) -> model_cost(M, p) * inv_scale,
+        (M, p) -> _scale_solver_tangent(model_grad(M, p), inv_scale),
+        true,
+    )
+end
+
+@inline _solver_has_converged(state) = Manopt.has_converged(state)
+
+
+function _solver_iterations(state, maxiter::Int)
+    if isdefined(Manopt, :stopped_at)
+        try
+            k = Manopt.stopped_at(state)
+            return k > 0 ? Int(k) : maxiter
+        catch
+        end
+    end
+    while hasproperty(state, :state)
+        state = state.state
+    end
+    return hasproperty(state, :stop) && hasproperty(state.stop, :at_iteration) ?
+           state.stop.at_iteration : maxiter
+end
+
+@inline function _solver_iteration_source()
+    return isdefined(Manopt, :stopped_at) ? :stopped_at : :stop_at_iteration_fallback
+end
+
+function _solver_stats(
+    model_cost,
+    model_grad,
+    M,
+    p_opt,
+    state,
+    ::Nothing;
+    tol_T,
+    maxiter::Int,
+    solver::Symbol,
+    tiny_grad_tol = nothing,
+    solver_info = (;),
+    use_state_gradient::Bool = true,
+)
+    T = typeof(tol_T)
+    final_cost = model_cost(M, p_opt)
+    cost_for_error = max(T(0), T(2) * final_cost)
+    rel_error = sqrt(cost_for_error)
+    grad_state = use_state_gradient ? _solver_gradient(state) : nothing
+    grad_from_state = !isnothing(grad_state)
+    grad_final =
+        isnothing(grad_state) ? model_grad(M, p_opt) :
+        _align_layout_like_point(p_opt, grad_state)
+    grad_norm = norm(M, p_opt, grad_final)
+    iterations = _solver_iterations(state, maxiter)
+    converged_grad =
+        grad_norm < tol_T || (!isnothing(tiny_grad_tol) && grad_norm < tiny_grad_tol)
+    converged_state = _solver_has_converged(state)
+    solver_info = merge(
+        solver_info,
+        (
+            gradient_source = grad_from_state ? :state : :recomputed,
+            has_converged_state = converged_state,
+            converged_by_gradient_threshold = converged_grad,
+            iteration_source = _solver_iteration_source(),
+        ),
+    )
+    return (
+        point = p_opt,
+        cost = final_cost,
+        rel_error = rel_error,
+        grad_norm = grad_norm,
+        iterations = iterations,
+        converged = converged_state,
+        solver = solver,
+        solver_info = solver_info,
+    )
+end
+
+function _solver_stats(
+    model_cost,
+    model_grad,
+    M,
+    p_opt,
+    state,
+    normA2::Real;
+    tol_T,
+    maxiter::Int,
+    solver::Symbol,
+    tiny_grad_tol = nothing,
+    solver_info = (;),
+    use_state_gradient::Bool = true,
+)
+    T = typeof(tol_T)
+    final_cost = model_cost(M, p_opt)
+    cost_for_error = max(T(0), T(2) * final_cost)
+    rel_error = _relative_error_frob_sq(cost_for_error, T(normA2))
+    grad_state = use_state_gradient ? _solver_gradient(state) : nothing
+    grad_from_state = !isnothing(grad_state)
+    grad_final =
+        isnothing(grad_state) ? model_grad(M, p_opt) :
+        _align_layout_like_point(p_opt, grad_state)
+    grad_norm = norm(M, p_opt, grad_final)
+    iterations = _solver_iterations(state, maxiter)
+    converged_grad =
+        grad_norm < tol_T || (!isnothing(tiny_grad_tol) && grad_norm < tiny_grad_tol)
+    converged_state = _solver_has_converged(state)
+    solver_info = merge(
+        solver_info,
+        (
+            gradient_source = grad_from_state ? :state : :recomputed,
+            has_converged_state = converged_state,
+            converged_by_gradient_threshold = converged_grad,
+            iteration_source = _solver_iteration_source(),
+        ),
+    )
+    return (
+        point = p_opt,
+        cost = final_cost,
+        rel_error = rel_error,
+        grad_norm = grad_norm,
+        iterations = iterations,
+        converged = converged_state,
+        solver = solver,
+        solver_info = solver_info,
+    )
+end
+
+_solver_debug_callbacks(callbacks...) = Any[cb for cb in callbacks if !isnothing(cb)]
+
+# Allow callers to pass `nothing` (e.g., when verbose/debug is omitted)
+_solver_debug_actions(::Nothing, callbacks...) = _solver_debug_callbacks(callbacks...)
+
+function _solver_debug_actions(verbose::Bool, callbacks...)
+    callback_actions = _solver_debug_callbacks(callbacks...)
+    if verbose
+        io = _SOLVER_DEBUG_SINK
+        init_group = Manopt.DebugGroup([
+            Manopt.DebugDivider("Initial "; io, at_init = true),
+            Manopt.DebugCost(; io, format = "f(x): %.6e", at_init = true),
+            Manopt.DebugGradientNorm(; io, format = "|grad f(p)|:%.6e", at_init = true),
+            Manopt.DebugDivider("\n"; io, at_init = true),
+        ])
+        iter_group = Manopt.DebugEvery(
+            Manopt.DebugGroup([
+                Manopt.DebugIteration(; io, format = "# %-6d"),
+                Manopt.DebugDivider(" "; io, at_init = true),
+                Manopt.DebugCost(; io, format = "f(x): %.6e", at_init = true),
+                Manopt.DebugGradientNorm(; io, format = "|grad f(p)|:%.6e", at_init = true),
+                Manopt.DebugDivider("\n"; io, at_init = true),
+            ]),
+            100,
+        )
+        iteration_actions = Any[iter_group]
+        append!(iteration_actions, callback_actions)
+        return Any[:Start=>Any[init_group], :Iteration=>iteration_actions]
+    end
+    return callback_actions
+end
+
+function _solver_progress_callback(
+    progress,
+    model_cost,
+    model_grad,
+    M;
+    normA2 = nothing,
+    diagnostics_recorder = nothing,
+)
+    progress isa NoMethodProgress && return nothing
+    has_relative_scale = !isnothing(normA2) && normA2 > 0
+    target_norm = has_relative_scale ? sqrt(normA2) : nothing
+    return function (problem, state, k)
+        k <= 0 && return nothing
+        p = get_iterate(state)
+        c = model_cost(M, p)
+        g = model_grad(M, p)
+        gnorm = norm(M, p, g)
+        c_display = has_relative_scale ? sqrt(max(2 * c, zero(c))) : c
+        gnorm_display = has_relative_scale ? gnorm * target_norm : gnorm
+        showvalues = Any[("Iter", k), ("Cost", c_display), ("Grad norm", gnorm_display)]
+        if !isnothing(diagnostics_recorder)
+            step = diagnostics_recorder.accepted_stepsize_history
+            trials = diagnostics_recorder.line_search_trial_history
+            !isempty(step) && push!(showvalues, ("Accepted α", step[end]))
+            !isempty(trials) && push!(showvalues, ("Line-search trials", trials[end]))
+        end
+        update_progress!(progress, k; showvalues)
+        return nothing
+    end
+end
+
+mutable struct _SolverDiagnosticsRecorder
+    first_accepted_stepsize::Float64
+    min_accepted_stepsize::Float64
+    first_line_search_trials::Int
+    line_search_trial_count::Int
+    function_evaluations::Int
+    gradient_evaluations::Int
+    prev_function_evaluations::Int
+    prev_gradient_evaluations::Int
+    line_search_enabled::Bool
+    fallback_stepsize::Float64
+    accepted_stepsize_history::Vector{Float64}
+    line_search_trial_history::Vector{Int}
+end
+
+function _SolverDiagnosticsRecorder(;
+    line_search_enabled::Bool,
+    fallback_stepsize::Real = NaN,
+)
+    return _SolverDiagnosticsRecorder(
+        NaN,
+        Inf,
+        0,
+        0,
+        0,
+        0,
+        0,
+        0,
+        line_search_enabled,
+        Float64(fallback_stepsize),
+        Float64[],
+        Int[],
+    )
+end
+
+function _solver_eval_count(problem, sym::Symbol)
+    try
+        count = get_count(get_objective(problem), sym)
+        return count < 0 ? 0 : Int(count)
+    catch
+        return 0
+    end
+end
+
+function _solver_diagnostics_callback(recorder::_SolverDiagnosticsRecorder)
+    return function (problem, state, k)
+        fe = _solver_eval_count(problem, :Cost)
+        ge = _solver_eval_count(problem, :Gradient)
+        if k == 0
+            recorder.function_evaluations = fe
+            recorder.gradient_evaluations = ge
+            recorder.prev_function_evaluations = fe
+            recorder.prev_gradient_evaluations = ge
+            return nothing
+        end
+        step =
+            recorder.line_search_enabled ? get_last_stepsize(problem, state, k) :
+            recorder.fallback_stepsize
+        delta_fe = max(fe - recorder.prev_function_evaluations, 0)
+        step_f = Float64(step)
+        ls_trials = recorder.line_search_enabled ? max(delta_fe - 1, 0) : 0
+        if isnan(recorder.first_accepted_stepsize)
+            recorder.first_accepted_stepsize = step_f
+            recorder.first_line_search_trials = ls_trials
+        end
+        recorder.min_accepted_stepsize = min(recorder.min_accepted_stepsize, step_f)
+        recorder.line_search_trial_count += ls_trials
+        push!(recorder.accepted_stepsize_history, step_f)
+        push!(recorder.line_search_trial_history, ls_trials)
+        recorder.function_evaluations = fe
+        recorder.gradient_evaluations = ge
+        recorder.prev_function_evaluations = fe
+        recorder.prev_gradient_evaluations = ge
+        return nothing
+    end
+end
+
+function _solver_info(recorder::_SolverDiagnosticsRecorder, iterations::Int)
+    return (
+        total_iterations = iterations,
+        first_accepted_stepsize = recorder.first_accepted_stepsize,
+        min_accepted_stepsize = recorder.min_accepted_stepsize,
+        first_line_search_trials = recorder.first_line_search_trials,
+        line_search_trial_count = recorder.line_search_trial_count,
+        function_evaluations = recorder.function_evaluations,
+        gradient_evaluations = recorder.gradient_evaluations,
+        accepted_stepsize_history = recorder.accepted_stepsize_history,
+        line_search_trial_history = recorder.line_search_trial_history,
+    )
+end
+
+function _solver_post_step_callback(
+    model::AbstractDecompositionModel,
+    M,
+    normalization::AbstractNormalizationPolicy,
+    solver_sym::Symbol,
+)
+    normalization isa NoNormalization && return nothing
+    return function (problem, state, k)
+        p_old = get_iterate(state)
+        # Backend postprocessing (for example normalization) runs in canonical
+        # CP coordinates and then packs back into the solver's current layout.
+        p_new = post_step!(
+            model,
+            p_old;
+            normalization,
+            solver = solver_sym,
+            problem,
+            state,
+            iteration = k,
+        )
+        p_new = _align_layout_like_point(p_old, p_new)
+        p_new === p_old && return nothing
+        set_iterate!(state, M, p_new)
+        if solver_sym == :rcg && hasproperty(state, :X)
+            get_gradient!(problem, state.X, get_iterate(state))
+            if hasproperty(state, :δ)
+                state.δ = -copy(M, get_iterate(state), state.X)
+            end
+            hasproperty(state, :β) && (state.β = zero(typeof(state.β)))
+            if hasproperty(state, :coefficient) && hasproperty(state.coefficient, :storage)
+                update_storage!(state.coefficient.storage, problem, state)
+            end
+        end
+        return nothing
+    end
+end
\ No newline at end of file
diff --git a/src/solvers/rgd.jl b/src/solvers/rgd.jl
index af180e9..be850a6 100644
--- a/src/solvers/rgd.jl
+++ b/src/solvers/rgd.jl
@@ -2,600 +2,6 @@
 export RGDSolver, RGDFixedSolver
 using Manopt
 
-struct _SolverDebugSink <: IO end
-Base.isopen(::_SolverDebugSink) = true
-Base.write(::_SolverDebugSink, ::UInt8) = 1
-Base.write(::_SolverDebugSink, s::Union{String,SubString{String}}) = sizeof(s)
-Base.unsafe_write(::_SolverDebugSink, ::Ptr{UInt8}, n::UInt) = Int(n)
-
-const _SOLVER_DEBUG_SINK = _SolverDebugSink()
-
-mutable struct StopWhenCostRelChangeAndGradientLess{T<:Real} <: Manopt.StoppingCriterion
-    tol_cost::T
-    tol_grad::T
-    prev_cost::T
-    last_cost_rel_change::T
-    last_grad_norm::T
-    at_iteration::Int
-end
-
-function StopWhenCostRelChangeAndGradientLess(tol_cost::T, tol_grad::T) where {T<:Real}
-    return StopWhenCostRelChangeAndGradientLess{T}(
-        tol_cost,
-        tol_grad,
-        T(Inf),
-        T(Inf),
-        T(Inf),
-        -1,
-    )
-end
-
-function (c::StopWhenCostRelChangeAndGradientLess)(problem, state, i)
-    if i == 0
-        c.prev_cost = Manopt.get_cost(problem, Manopt.get_iterate(state))
-        c.last_cost_rel_change = oftype(c.tol_cost, Inf)
-        c.last_grad_norm = oftype(c.tol_grad, Inf)
-        c.at_iteration = -1
-        return false
-    end
-    M = Manopt.get_manifold(problem)
-    p = Manopt.get_iterate(state)
-    cost_val = Manopt.get_cost(problem, p)
-    grad_val = Manopt.get_gradient(problem, p)
-    grad_norm = norm(M, p, grad_val)
-    rel_change = abs(c.prev_cost - cost_val) / max(abs(c.prev_cost), one(cost_val))
-    c.prev_cost = cost_val
-    c.last_cost_rel_change = rel_change
-    c.last_grad_norm = grad_norm
-    if rel_change < c.tol_cost && grad_norm < c.tol_grad
-        c.at_iteration = i
-        return true
-    end
-    return false
-end
-
-function Manopt.get_reason(c::StopWhenCostRelChangeAndGradientLess)
-    if c.at_iteration >= 0
-        return "At iteration $(c.at_iteration) the relative cost change ($(c.last_cost_rel_change)) " *
-               "is below $(c.tol_cost) and the gradient norm ($(c.last_grad_norm)) " *
-               "is below $(c.tol_grad).\n"
-    end
-    return ""
-end
-
-function Manopt.status_summary(c::StopWhenCostRelChangeAndGradientLess)
-    has_stopped = c.at_iteration >= 0
-    status = has_stopped ? "reached" : "not reached"
-    return "cost rel change < $(c.tol_cost) and |grad f| < $(c.tol_grad): $status"
-end
-
-Manopt.indicates_convergence(::StopWhenCostRelChangeAndGradientLess) = true
-
-function Base.show(io::IO, c::StopWhenCostRelChangeAndGradientLess)
-    return print(
-        io,
-        "StopWhenCostRelChangeAndGradientLess($(c.tol_cost), $(c.tol_grad))\n    $(Manopt.status_summary(c))",
-    )
-end
-
-
-function _tk_get_solver_result(state)
-    try
-        return Manopt.get_solver_result(state)
-    catch
-    end
-    while hasproperty(state, :state)
-        state = state.state
-    end
-    for key in (:p, :x, :point)
-        hasproperty(state, key) && return getproperty(state, key)
-    end
-    throw(
-        ArgumentError(
-            "Cannot extract point from state $(typeof(state)). Properties: $(propertynames(state))",
-        ),
-    )
-end
-
-
-@inline _align_layout_like_point(p, x) =
-    hasproperty(p, :x) ?
-    (hasproperty(x, :x) ? x : (x isa Tuple ? ArrayPartition(x...) : x)) :
-    (hasproperty(x, :x) ? Tuple(getproperty(x, :x)) : x)
-
-
-function _to_array_partition(x)
-    if x isa ArrayPartition
-        return ArrayPartition(map(_to_array_partition, x.x)...)
-    elseif hasproperty(x, :x)
-        return ArrayPartition(map(_to_array_partition, getproperty(x, :x))...)
-    elseif x isa Tuple
-        return ArrayPartition(map(_to_array_partition, x)...)
-    end
-    return x
-end
-
-
-function _solver_point(M, p0)
-    M2 = _unwrap_solver_manifold(M)
-    return M2 isa ProductManifold ? _to_array_partition(p0) : p0
-end
-
-
-function _contains_sqeuclidean_manifold(M)
-    M2 = _unwrap_solver_manifold(M)
-    if M2 isa SqEuclidean || M2 isa SoftplusEuclidean
-        return true
-    elseif M2 isa ProductManifold
-        return any(_contains_sqeuclidean_manifold, M2.manifolds)
-    elseif hasproperty(M2, :native) && (
-        getproperty(M2, :native) isa SqEuclidean ||
-        getproperty(M2, :native) isa SoftplusEuclidean
-    )
-        return true
-    end
-    return false
-end
-
-function _contains_strict_sqeuclidean_manifold(M)
-    M2 = _unwrap_solver_manifold(M)
-    if M2 isa SqEuclidean
-        return true
-    elseif M2 isa ProductManifold
-        return any(_contains_strict_sqeuclidean_manifold, M2.manifolds)
-    elseif hasproperty(M2, :native) && (getproperty(M2, :native) isa SqEuclidean)
-        return true
-    end
-    return false
-end
-
-
-function _armijo_max_decreases(initial_stepsize::Real, contraction::Real, alpha_min::Real)
-    initial_stepsize <= alpha_min && return 0
-    (contraction <= 0 || contraction >= 1) && return 1000
-    n = floor(Int, log(alpha_min / initial_stepsize) / log(contraction))
-    return max(n, 0)
-end
-
-
-function _adaptive_initial_stepsize(
-    M,
-    p0,
-    model_grad,
-    retraction_method,
-    base_stepsize::T;
-    alpha_min::T,
-    scale_c::T = one(T),
-    clamp_low_factor::T = T(0.1),
-    clamp_high_factor::T = T(10),
-    delta_scale::T = T(1e-3),
-) where {T<:AbstractFloat}
-    g0 = model_grad(M, p0)
-    d = -copy(g0)
-    dnorm = norm(M, p0, d)
-    (!isfinite(dnorm) || dnorm <= sqrt(eps(T))) && return base_stepsize
-    δ = delta_scale / max(dnorm, one(T))
-    q = try
-        retract(M, p0, δ .* d, retraction_method)
-    catch
-        return base_stepsize
-    end
-    _all_finite(q) || return base_stepsize
-    gq = model_grad(M, q)
-    _all_finite(gq) || return base_stepsize
-    κ_num = inner(M, p0, gq .- g0, d)
-    κ_den = δ * dnorm^2
-    (!isfinite(κ_num) || !isfinite(κ_den) || κ_den <= eps(T)) && return base_stepsize
-    κ = max(κ_num / κ_den, eps(T))
-    α_raw = scale_c / κ
-    α_low = max(alpha_min, clamp_low_factor * base_stepsize)
-    α_high = clamp_high_factor * base_stepsize
-    α = clamp(α_raw, α_low, α_high)
-    if !isfinite(α)
-        return base_stepsize
-    end
-    return α
-end
-
-function _all_finite(x)
-    if x isa Number
-        return isfinite(x)
-    elseif hasproperty(x, :x)
-        return all(_all_finite, getproperty(x, :x))
-    elseif x isa AbstractArray
-        return all(isfinite, x)
-    elseif x isa Tuple
-        return all(_all_finite, x)
-    elseif x isa Manifolds.TuckerPoint
-        return _all_finite(x.hosvd.core) && all(_all_finite, x.hosvd.U)
-    elseif x isa Manifolds.TuckerTangentVector
-        return _all_finite(getproperty(x, :Ċ)) && all(_all_finite, getproperty(x, :U̇))
-    end
-    try
-        return all(isfinite, x)
-    catch
-        return false
-    end
-end
-
-
-function _safe_cost_function(model_cost)
-    return function (M, p)
-        _all_finite(p) || return Inf
-        c = model_cost(M, p)
-        return isfinite(c) ? c : Inf
-    end
-end
-
-
-function _layout_adapt_gradient(model_grad)
-    return function (M, p)
-        g = model_grad(M, p)
-        return _align_layout_like_point(p, g)
-    end
-end
-
-function _scalar_eltype(p)
-    if hasproperty(p, :x) || p isa AbstractVector || p isa Tuple
-        parts = point_parts(p)
-        isempty(parts) && return Float64
-        return _scalar_eltype(first(parts))
-    elseif p isa Manifolds.TuckerPoint
-        return eltype(p.hosvd.core)
-    elseif p isa Manifolds.TuckerTangentVector
-        return eltype(getproperty(p, :Ċ))
-    elseif p isa Real
-        return typeof(p)
-    else
-        return eltype(p)
-    end
-end
-
-@inline function _solver_gradient(state)
-    try
-        return Manopt.get_gradient(state)
-    catch
-        return nothing
-    end
-end
-
-@inline _scale_solver_tangent(x::Number, scale::Real) = x * scale
-_scale_solver_tangent(x::AbstractArray, scale::Real) = x .* scale
-_scale_solver_tangent(x::ArrayPartition, scale::Real) =
-    ArrayPartition(map(part -> _scale_solver_tangent(part, scale), x.x)...)
-_scale_solver_tangent(x::Tuple, scale::Real) =
-    map(part -> _scale_solver_tangent(part, scale), x)
-
-function _scale_solver_tangent(x, scale::Real)
-    try
-        return x .* scale
-    catch
-        return scale * x
-    end
-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)
-    return (
-        (M, p) -> model_cost(M, p) * inv_scale,
-        (M, p) -> _scale_solver_tangent(model_grad(M, p), inv_scale),
-        true,
-    )
-end
-
-@inline _solver_has_converged(state) = Manopt.has_converged(state)
-
-
-function _solver_iterations(state, maxiter::Int)
-    if isdefined(Manopt, :stopped_at)
-        try
-            k = Manopt.stopped_at(state)
-            return k > 0 ? Int(k) : maxiter
-        catch
-        end
-    end
-    while hasproperty(state, :state)
-        state = state.state
-    end
-    return hasproperty(state, :stop) && hasproperty(state.stop, :at_iteration) ?
-           state.stop.at_iteration : maxiter
-end
-
-@inline function _solver_iteration_source()
-    return isdefined(Manopt, :stopped_at) ? :stopped_at : :stop_at_iteration_fallback
-end
-
-function _solver_stats(
-    model_cost,
-    model_grad,
-    M,
-    p_opt,
-    state,
-    ::Nothing;
-    tol_T,
-    maxiter::Int,
-    solver::Symbol,
-    tiny_grad_tol = nothing,
-    solver_info = (;),
-    use_state_gradient::Bool = true,
-)
-    T = typeof(tol_T)
-    final_cost = model_cost(M, p_opt)
-    cost_for_error = max(T(0), T(2) * final_cost)
-    rel_error = sqrt(cost_for_error)
-    grad_state = use_state_gradient ? _solver_gradient(state) : nothing
-    grad_from_state = !isnothing(grad_state)
-    grad_final =
-        isnothing(grad_state) ? model_grad(M, p_opt) :
-        _align_layout_like_point(p_opt, grad_state)
-    grad_norm = norm(M, p_opt, grad_final)
-    iterations = _solver_iterations(state, maxiter)
-    converged_grad =
-        grad_norm < tol_T || (!isnothing(tiny_grad_tol) && grad_norm < tiny_grad_tol)
-    converged_state = _solver_has_converged(state)
-    solver_info = merge(
-        solver_info,
-        (
-            gradient_source = grad_from_state ? :state : :recomputed,
-            has_converged_state = converged_state,
-            converged_by_gradient_threshold = converged_grad,
-            iteration_source = _solver_iteration_source(),
-        ),
-    )
-    return (
-        point = p_opt,
-        cost = final_cost,
-        rel_error = rel_error,
-        grad_norm = grad_norm,
-        iterations = iterations,
-        converged = converged_state,
-        solver = solver,
-        solver_info = solver_info,
-    )
-end
-
-function _solver_stats(
-    model_cost,
-    model_grad,
-    M,
-    p_opt,
-    state,
-    normA2::Real;
-    tol_T,
-    maxiter::Int,
-    solver::Symbol,
-    tiny_grad_tol = nothing,
-    solver_info = (;),
-    use_state_gradient::Bool = true,
-)
-    T = typeof(tol_T)
-    final_cost = model_cost(M, p_opt)
-    cost_for_error = max(T(0), T(2) * final_cost)
-    rel_error = _relative_error_frob_sq(cost_for_error, T(normA2))
-    grad_state = use_state_gradient ? _solver_gradient(state) : nothing
-    grad_from_state = !isnothing(grad_state)
-    grad_final =
-        isnothing(grad_state) ? model_grad(M, p_opt) :
-        _align_layout_like_point(p_opt, grad_state)
-    grad_norm = norm(M, p_opt, grad_final)
-    iterations = _solver_iterations(state, maxiter)
-    converged_grad =
-        grad_norm < tol_T || (!isnothing(tiny_grad_tol) && grad_norm < tiny_grad_tol)
-    converged_state = _solver_has_converged(state)
-    solver_info = merge(
-        solver_info,
-        (
-            gradient_source = grad_from_state ? :state : :recomputed,
-            has_converged_state = converged_state,
-            converged_by_gradient_threshold = converged_grad,
-            iteration_source = _solver_iteration_source(),
-        ),
-    )
-    return (
-        point = p_opt,
-        cost = final_cost,
-        rel_error = rel_error,
-        grad_norm = grad_norm,
-        iterations = iterations,
-        converged = converged_state,
-        solver = solver,
-        solver_info = solver_info,
-    )
-end
-
-_solver_debug_callbacks(callbacks...) = Any[cb for cb in callbacks if !isnothing(cb)]
-
-# Allow callers to pass `nothing` (e.g., when verbose/debug is omitted)
-_solver_debug_actions(::Nothing, callbacks...) = _solver_debug_callbacks(callbacks...)
-
-function _solver_debug_actions(verbose::Bool, callbacks...)
-    callback_actions = _solver_debug_callbacks(callbacks...)
-    if verbose
-        io = _SOLVER_DEBUG_SINK
-        init_group = Manopt.DebugGroup([
-            Manopt.DebugDivider("Initial "; io, at_init = true),
-            Manopt.DebugCost(; io, format = "f(x): %.6e", at_init = true),
-            Manopt.DebugGradientNorm(; io, format = "|grad f(p)|:%.6e", at_init = true),
-            Manopt.DebugDivider("\n"; io, at_init = true),
-        ])
-        iter_group = Manopt.DebugEvery(
-            Manopt.DebugGroup([
-                Manopt.DebugIteration(; io, format = "# %-6d"),
-                Manopt.DebugDivider(" "; io, at_init = true),
-                Manopt.DebugCost(; io, format = "f(x): %.6e", at_init = true),
-                Manopt.DebugGradientNorm(; io, format = "|grad f(p)|:%.6e", at_init = true),
-                Manopt.DebugDivider("\n"; io, at_init = true),
-            ]),
-            100,
-        )
-        iteration_actions = Any[iter_group]
-        append!(iteration_actions, callback_actions)
-        return Any[:Start=>Any[init_group], :Iteration=>iteration_actions]
-    end
-    return callback_actions
-end
-
-function _solver_progress_callback(
-    progress,
-    model_cost,
-    model_grad,
-    M;
-    normA2 = nothing,
-    diagnostics_recorder = nothing,
-)
-    progress isa NoMethodProgress && return nothing
-    has_relative_scale = !isnothing(normA2) && normA2 > 0
-    target_norm = has_relative_scale ? sqrt(normA2) : nothing
-    return function (problem, state, k)
-        k <= 0 && return nothing
-        p = get_iterate(state)
-        c = model_cost(M, p)
-        g = model_grad(M, p)
-        gnorm = norm(M, p, g)
-        c_display = has_relative_scale ? sqrt(max(2 * c, zero(c))) : c
-        gnorm_display = has_relative_scale ? gnorm * target_norm : gnorm
-        showvalues = Any[("Iter", k), ("Cost", c_display), ("Grad norm", gnorm_display)]
-        if !isnothing(diagnostics_recorder)
-            step = diagnostics_recorder.accepted_stepsize_history
-            trials = diagnostics_recorder.line_search_trial_history
-            !isempty(step) && push!(showvalues, ("Accepted α", step[end]))
-            !isempty(trials) && push!(showvalues, ("Line-search trials", trials[end]))
-        end
-        update_progress!(progress, k; showvalues)
-        return nothing
-    end
-end
-
-mutable struct _SolverDiagnosticsRecorder
-    first_accepted_stepsize::Float64
-    min_accepted_stepsize::Float64
-    first_line_search_trials::Int
-    line_search_trial_count::Int
-    function_evaluations::Int
-    gradient_evaluations::Int
-    prev_function_evaluations::Int
-    prev_gradient_evaluations::Int
-    line_search_enabled::Bool
-    fallback_stepsize::Float64
-    accepted_stepsize_history::Vector{Float64}
-    line_search_trial_history::Vector{Int}
-end
-
-function _SolverDiagnosticsRecorder(;
-    line_search_enabled::Bool,
-    fallback_stepsize::Real = NaN,
-)
-    return _SolverDiagnosticsRecorder(
-        NaN,
-        Inf,
-        0,
-        0,
-        0,
-        0,
-        0,
-        0,
-        line_search_enabled,
-        Float64(fallback_stepsize),
-        Float64[],
-        Int[],
-    )
-end
-
-function _solver_eval_count(problem, sym::Symbol)
-    try
-        count = get_count(get_objective(problem), sym)
-        return count < 0 ? 0 : Int(count)
-    catch
-        return 0
-    end
-end
-
-function _solver_diagnostics_callback(recorder::_SolverDiagnosticsRecorder)
-    return function (problem, state, k)
-        fe = _solver_eval_count(problem, :Cost)
-        ge = _solver_eval_count(problem, :Gradient)
-        if k == 0
-            recorder.function_evaluations = fe
-            recorder.gradient_evaluations = ge
-            recorder.prev_function_evaluations = fe
-            recorder.prev_gradient_evaluations = ge
-            return nothing
-        end
-        step =
-            recorder.line_search_enabled ? get_last_stepsize(problem, state, k) :
-            recorder.fallback_stepsize
-        delta_fe = max(fe - recorder.prev_function_evaluations, 0)
-        step_f = Float64(step)
-        ls_trials = recorder.line_search_enabled ? max(delta_fe - 1, 0) : 0
-        if isnan(recorder.first_accepted_stepsize)
-            recorder.first_accepted_stepsize = step_f
-            recorder.first_line_search_trials = ls_trials
-        end
-        recorder.min_accepted_stepsize = min(recorder.min_accepted_stepsize, step_f)
-        recorder.line_search_trial_count += ls_trials
-        push!(recorder.accepted_stepsize_history, step_f)
-        push!(recorder.line_search_trial_history, ls_trials)
-        recorder.function_evaluations = fe
-        recorder.gradient_evaluations = ge
-        recorder.prev_function_evaluations = fe
-        recorder.prev_gradient_evaluations = ge
-        return nothing
-    end
-end
-
-function _solver_info(recorder::_SolverDiagnosticsRecorder, iterations::Int)
-    return (
-        total_iterations = iterations,
-        first_accepted_stepsize = recorder.first_accepted_stepsize,
-        min_accepted_stepsize = recorder.min_accepted_stepsize,
-        first_line_search_trials = recorder.first_line_search_trials,
-        line_search_trial_count = recorder.line_search_trial_count,
-        function_evaluations = recorder.function_evaluations,
-        gradient_evaluations = recorder.gradient_evaluations,
-        accepted_stepsize_history = recorder.accepted_stepsize_history,
-        line_search_trial_history = recorder.line_search_trial_history,
-    )
-end
-
-function _solver_post_step_callback(
-    model::AbstractDecompositionModel,
-    M,
-    normalization::AbstractNormalizationPolicy,
-    solver_sym::Symbol,
-)
-    normalization isa NoNormalization && return nothing
-    return function (problem, state, k)
-        p_old = get_iterate(state)
-        # Backend postprocessing (for example normalization) runs in canonical
-        # CP coordinates and then packs back into the solver's current layout.
-        p_new = post_step!(
-            model,
-            p_old;
-            normalization,
-            solver = solver_sym,
-            problem,
-            state,
-            iteration = k,
-        )
-        p_new = _align_layout_like_point(p_old, p_new)
-        p_new === p_old && return nothing
-        set_iterate!(state, M, p_new)
-        if solver_sym == :rcg && hasproperty(state, :X)
-            get_gradient!(problem, state.X, get_iterate(state))
-            if hasproperty(state, :δ)
-                state.δ = -copy(M, get_iterate(state), state.X)
-            end
-            hasproperty(state, :β) && (state.β = zero(typeof(state.β)))
-            if hasproperty(state, :coefficient) && hasproperty(state.coefficient, :storage)
-                update_storage!(state.coefficient.storage, problem, state)
-            end
-        end
-        return nothing
-    end
-end
-
 function solve_rgd(
     model_cost,
     model_egrad,

From b1c73cb85445fee64d8f679a0d83b18427e4a8d1 Mon Sep 17 00:00:00 2001
From: PBrdng <pbrdng@gmail.com>
Date: Wed, 17 Jun 2026 20:45:43 +0200
Subject: [PATCH 06/18] add helper documentation

---
 src/solvers/manopt_helpers.jl | 47 ++++++++++++++++++++++++++++++++---
 1 file changed, 43 insertions(+), 4 deletions(-)

diff --git a/src/solvers/manopt_helpers.jl b/src/solvers/manopt_helpers.jl
index a107e27..8c0f7bd 100644
--- a/src/solvers/manopt_helpers.jl
+++ b/src/solvers/manopt_helpers.jl
@@ -1,12 +1,19 @@
 
+######################## 
+# This file contains helper functions for communicating with Manopt
+########################
+
+# Sink Manopt's built-in debug text while still letting debug actions run.
 struct _SolverDebugSink <: IO end
 Base.isopen(::_SolverDebugSink) = true
 Base.write(::_SolverDebugSink, ::UInt8) = 1
 Base.write(::_SolverDebugSink, s::Union{String,SubString{String}}) = sizeof(s)
 Base.unsafe_write(::_SolverDebugSink, ::Ptr{UInt8}, n::UInt) = Int(n)
 
+# Shared no-op IO used by Manopt debug groups.
 const _SOLVER_DEBUG_SINK = _SolverDebugSink()
 
+# Stop when both the relative cost change and Riemannian gradient norm are small.
 mutable struct StopWhenCostRelChangeAndGradientLess{T<:Real} <: Manopt.StoppingCriterion
     tol_cost::T
     tol_grad::T
@@ -16,6 +23,7 @@ mutable struct StopWhenCostRelChangeAndGradientLess{T<:Real} <: Manopt.StoppingC
     at_iteration::Int
 end
 
+# Initialize the dual cost-change/gradient stopping rule with empty history.
 function StopWhenCostRelChangeAndGradientLess(tol_cost::T, tol_grad::T) where {T<:Real}
     return StopWhenCostRelChangeAndGradientLess{T}(
         tol_cost,
@@ -27,6 +35,7 @@ function StopWhenCostRelChangeAndGradientLess(tol_cost::T, tol_grad::T) where {T
     )
 end
 
+# Update the dual stopping rule from the current Manopt problem/state.
 function (c::StopWhenCostRelChangeAndGradientLess)(problem, state, i)
     if i == 0
         c.prev_cost = Manopt.get_cost(problem, Manopt.get_iterate(state))
@@ -51,6 +60,7 @@ function (c::StopWhenCostRelChangeAndGradientLess)(problem, state, i)
     return false
 end
 
+# Explain why the dual stopping rule stopped, for Manopt status reporting.
 function Manopt.get_reason(c::StopWhenCostRelChangeAndGradientLess)
     if c.at_iteration >= 0
         return "At iteration $(c.at_iteration) the relative cost change ($(c.last_cost_rel_change)) " *
@@ -60,14 +70,17 @@ function Manopt.get_reason(c::StopWhenCostRelChangeAndGradientLess)
     return ""
 end
 
+# Summarize the current dual stopping-rule state for Manopt displays.
 function Manopt.status_summary(c::StopWhenCostRelChangeAndGradientLess)
     has_stopped = c.at_iteration >= 0
     status = has_stopped ? "reached" : "not reached"
     return "cost rel change < $(c.tol_cost) and |grad f| < $(c.tol_grad): $status"
 end
 
+# Mark the dual stopping rule as convergence, not failure or exhaustion.
 Manopt.indicates_convergence(::StopWhenCostRelChangeAndGradientLess) = true
 
+# Print the dual stopping rule compactly in Manopt diagnostics.
 function Base.show(io::IO, c::StopWhenCostRelChangeAndGradientLess)
     return print(
         io,
@@ -76,6 +89,7 @@ function Base.show(io::IO, c::StopWhenCostRelChangeAndGradientLess)
 end
 
 
+# Extract the final iterate across Manopt versions and nested state wrappers.
 function _tk_get_solver_result(state)
     try
         return Manopt.get_solver_result(state)
@@ -95,12 +109,14 @@ function _tk_get_solver_result(state)
 end
 
 
+# Match gradients/tangents to the point container Manopt is currently using.
 @inline _align_layout_like_point(p, x) =
     hasproperty(p, :x) ?
     (hasproperty(x, :x) ? x : (x isa Tuple ? ArrayPartition(x...) : x)) :
     (hasproperty(x, :x) ? Tuple(getproperty(x, :x)) : x)
 
 
+# Recursively convert tuple-like product points to ArrayPartition layout.
 function _to_array_partition(x)
     if x isa ArrayPartition
         return ArrayPartition(map(_to_array_partition, x.x)...)
@@ -113,12 +129,14 @@ function _to_array_partition(x)
 end
 
 
+# Adapt an initial point to the layout expected by the solver manifold.
 function _solver_point(M, p0)
     M2 = _unwrap_solver_manifold(M)
     return M2 isa ProductManifold ? _to_array_partition(p0) : p0
 end
 
 
+# Detect pullback nonnegative geometries that need conservative line search.
 function _contains_sqeuclidean_manifold(M)
     M2 = _unwrap_solver_manifold(M)
     if M2 isa SqEuclidean || M2 isa SoftplusEuclidean
@@ -134,6 +152,7 @@ function _contains_sqeuclidean_manifold(M)
     return false
 end
 
+# Detect strict squaring geometries that require extra Armijo safeguards.
 function _contains_strict_sqeuclidean_manifold(M)
     M2 = _unwrap_solver_manifold(M)
     if M2 isa SqEuclidean
@@ -147,6 +166,7 @@ function _contains_strict_sqeuclidean_manifold(M)
 end
 
 
+# Compute how many Armijo contractions are needed before alpha_min is reached.
 function _armijo_max_decreases(initial_stepsize::Real, contraction::Real, alpha_min::Real)
     initial_stepsize <= alpha_min && return 0
     (contraction <= 0 || contraction >= 1) && return 1000
@@ -155,6 +175,7 @@ function _armijo_max_decreases(initial_stepsize::Real, contraction::Real, alpha_
 end
 
 
+# Estimate a guarded initial RGD step from a one-sided curvature probe.
 function _adaptive_initial_stepsize(
     M,
     p0,
@@ -194,6 +215,7 @@ function _adaptive_initial_stepsize(
     return α
 end
 
+# Check recursively that points, tangents, and nested manifold containers are finite.
 function _all_finite(x)
     if x isa Number
         return isfinite(x)
@@ -216,6 +238,7 @@ function _all_finite(x)
 end
 
 
+# Wrap a cost so invalid points return Inf instead of poisoning line search.
 function _safe_cost_function(model_cost)
     return function (M, p)
         _all_finite(p) || return Inf
@@ -225,6 +248,7 @@ function _safe_cost_function(model_cost)
 end
 
 
+# Wrap a gradient so its container layout follows the queried point layout.
 function _layout_adapt_gradient(model_grad)
     return function (M, p)
         g = model_grad(M, p)
@@ -232,6 +256,7 @@ function _layout_adapt_gradient(model_grad)
     end
 end
 
+# Infer the scalar element type from nested points/tangents used by solvers.
 function _scalar_eltype(p)
     if hasproperty(p, :x) || p isa AbstractVector || p isa Tuple
         parts = point_parts(p)
@@ -248,6 +273,7 @@ function _scalar_eltype(p)
     end
 end
 
+# Read Manopt's cached gradient when the current state exposes it.
 @inline function _solver_gradient(state)
     try
         return Manopt.get_gradient(state)
@@ -256,13 +282,13 @@ end
     end
 end
 
+# Scale cost and gradient when using relative error  
 @inline _scale_solver_tangent(x::Number, scale::Real) = x * scale
 _scale_solver_tangent(x::AbstractArray, scale::Real) = x .* scale
 _scale_solver_tangent(x::ArrayPartition, scale::Real) =
     ArrayPartition(map(part -> _scale_solver_tangent(part, scale), x.x)...)
 _scale_solver_tangent(x::Tuple, scale::Real) =
     map(part -> _scale_solver_tangent(part, scale), x)
-
 function _scale_solver_tangent(x, scale::Real)
     try
         return x .* scale
@@ -270,7 +296,6 @@ function _scale_solver_tangent(x, scale::Real)
         return scale * x
     end
 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
@@ -282,9 +307,10 @@ function _relative_solver_functions(model_cost, model_grad, scale::Real)
     )
 end
 
+# Query Manopt's convergence flag through one local compatibility point.
 @inline _solver_has_converged(state) = Manopt.has_converged(state)
 
-
+# Recover the stopping iteration across Manopt versions and wrapped states.
 function _solver_iterations(state, maxiter::Int)
     if isdefined(Manopt, :stopped_at)
         try
@@ -300,10 +326,12 @@ function _solver_iterations(state, maxiter::Int)
            state.stop.at_iteration : maxiter
 end
 
+# Report which Manopt iteration API was used for solver metadata.
 @inline function _solver_iteration_source()
     return isdefined(Manopt, :stopped_at) ? :stopped_at : :stop_at_iteration_fallback
 end
 
+# Build common solver result stats when no target norm is available.
 function _solver_stats(
     model_cost,
     model_grad,
@@ -353,6 +381,7 @@ function _solver_stats(
     )
 end
 
+# Build common solver result stats and relative error when ||A||^2 is available.
 function _solver_stats(
     model_cost,
     model_grad,
@@ -402,11 +431,14 @@ function _solver_stats(
     )
 end
 
+# Collect only active Manopt debug callbacks, dropping omitted hooks.
 _solver_debug_callbacks(callbacks...) = Any[cb for cb in callbacks if !isnothing(cb)]
 
 # Allow callers to pass `nothing` (e.g., when verbose/debug is omitted)
+# Build debug actions when the verbose flag was omitted.
 _solver_debug_actions(::Nothing, callbacks...) = _solver_debug_callbacks(callbacks...)
 
+# Build Manopt debug actions and attach TensorKitchen callback hooks.
 function _solver_debug_actions(verbose::Bool, callbacks...)
     callback_actions = _solver_debug_callbacks(callbacks...)
     if verbose
@@ -434,6 +466,7 @@ function _solver_debug_actions(verbose::Bool, callbacks...)
     return callback_actions
 end
 
+# Create a Manopt iteration callback that updates TensorKitchen progress output.
 function _solver_progress_callback(
     progress,
     model_cost,
@@ -465,6 +498,7 @@ function _solver_progress_callback(
     end
 end
 
+# Accumulate line-search, step-size, and evaluation diagnostics during solving.
 mutable struct _SolverDiagnosticsRecorder
     first_accepted_stepsize::Float64
     min_accepted_stepsize::Float64
@@ -480,6 +514,7 @@ mutable struct _SolverDiagnosticsRecorder
     line_search_trial_history::Vector{Int}
 end
 
+# Initialize a diagnostics recorder for solvers with or without line search.
 function _SolverDiagnosticsRecorder(;
     line_search_enabled::Bool,
     fallback_stepsize::Real = NaN,
@@ -500,6 +535,7 @@ function _SolverDiagnosticsRecorder(;
     )
 end
 
+# Read Manopt objective evaluation counters defensively across configurations.
 function _solver_eval_count(problem, sym::Symbol)
     try
         count = get_count(get_objective(problem), sym)
@@ -509,6 +545,7 @@ function _solver_eval_count(problem, sym::Symbol)
     end
 end
 
+# Create a Manopt callback that records per-iteration solver diagnostics.
 function _solver_diagnostics_callback(recorder::_SolverDiagnosticsRecorder)
     return function (problem, state, k)
         fe = _solver_eval_count(problem, :Cost)
@@ -542,6 +579,7 @@ function _solver_diagnostics_callback(recorder::_SolverDiagnosticsRecorder)
     end
 end
 
+# Convert recorded diagnostics to the public solver_info named tuple.
 function _solver_info(recorder::_SolverDiagnosticsRecorder, iterations::Int)
     return (
         total_iterations = iterations,
@@ -556,6 +594,7 @@ function _solver_info(recorder::_SolverDiagnosticsRecorder, iterations::Int)
     )
 end
 
+# Create a callback that normalizes/postprocesses iterates after each solver step.
 function _solver_post_step_callback(
     model::AbstractDecompositionModel,
     M,
@@ -591,4 +630,4 @@ function _solver_post_step_callback(
         end
         return nothing
     end
-end
\ No newline at end of file
+end

From f76d069608bf6dab442c9391c4b8dc40b680aea7 Mon Sep 17 00:00:00 2001
From: PBrdng <pbrdng@gmail.com>
Date: Wed, 17 Jun 2026 20:49:00 +0200
Subject: [PATCH 07/18] clean code

---
 src/solvers/manopt_helpers.jl | 59 +++++------------------------------
 1 file changed, 7 insertions(+), 52 deletions(-)

diff --git a/src/solvers/manopt_helpers.jl b/src/solvers/manopt_helpers.jl
index 8c0f7bd..b997a24 100644
--- a/src/solvers/manopt_helpers.jl
+++ b/src/solvers/manopt_helpers.jl
@@ -331,64 +331,19 @@ end
     return isdefined(Manopt, :stopped_at) ? :stopped_at : :stop_at_iteration_fallback
 end
 
-# Build common solver result stats when no target norm is available.
-function _solver_stats(
-    model_cost,
-    model_grad,
-    M,
-    p_opt,
-    state,
-    ::Nothing;
-    tol_T,
-    maxiter::Int,
-    solver::Symbol,
-    tiny_grad_tol = nothing,
-    solver_info = (;),
-    use_state_gradient::Bool = true,
-)
-    T = typeof(tol_T)
-    final_cost = model_cost(M, p_opt)
-    cost_for_error = max(T(0), T(2) * final_cost)
-    rel_error = sqrt(cost_for_error)
-    grad_state = use_state_gradient ? _solver_gradient(state) : nothing
-    grad_from_state = !isnothing(grad_state)
-    grad_final =
-        isnothing(grad_state) ? model_grad(M, p_opt) :
-        _align_layout_like_point(p_opt, grad_state)
-    grad_norm = norm(M, p_opt, grad_final)
-    iterations = _solver_iterations(state, maxiter)
-    converged_grad =
-        grad_norm < tol_T || (!isnothing(tiny_grad_tol) && grad_norm < tiny_grad_tol)
-    converged_state = _solver_has_converged(state)
-    solver_info = merge(
-        solver_info,
-        (
-            gradient_source = grad_from_state ? :state : :recomputed,
-            has_converged_state = converged_state,
-            converged_by_gradient_threshold = converged_grad,
-            iteration_source = _solver_iteration_source(),
-        ),
-    )
-    return (
-        point = p_opt,
-        cost = final_cost,
-        rel_error = rel_error,
-        grad_norm = grad_norm,
-        iterations = iterations,
-        converged = converged_state,
-        solver = solver,
-        solver_info = solver_info,
-    )
-end
+# Convert a squared residual value to the public relative-error diagnostic.
+_solver_rel_error(cost_for_error, ::Nothing, ::Type) = sqrt(cost_for_error)
+_solver_rel_error(cost_for_error, normA2::Real, ::Type{T}) where {T} =
+    _relative_error_frob_sq(cost_for_error, T(normA2))
 
-# Build common solver result stats and relative error when ||A||^2 is available.
+# Build common solver result stats, with optional ||A||^2 for relative scaling.
 function _solver_stats(
     model_cost,
     model_grad,
     M,
     p_opt,
     state,
-    normA2::Real;
+    normA2::Union{Nothing,Real};
     tol_T,
     maxiter::Int,
     solver::Symbol,
@@ -399,7 +354,7 @@ function _solver_stats(
     T = typeof(tol_T)
     final_cost = model_cost(M, p_opt)
     cost_for_error = max(T(0), T(2) * final_cost)
-    rel_error = _relative_error_frob_sq(cost_for_error, T(normA2))
+    rel_error = _solver_rel_error(cost_for_error, normA2, T)
     grad_state = use_state_gradient ? _solver_gradient(state) : nothing
     grad_from_state = !isnothing(grad_state)
     grad_final =

From ec80afdc028211015cbb9841a2286c5ba76f8375 Mon Sep 17 00:00:00 2001
From: PBrdng <pbrdng@gmail.com>
Date: Wed, 17 Jun 2026 20:53:32 +0200
Subject: [PATCH 08/18] unify code

---
 src/api/approx.jl             |  3 +--
 src/solvers/lbfgs.jl          | 17 +----------------
 src/solvers/manopt_helpers.jl | 16 ++++++----------
 src/solvers/rcg.jl            | 17 +----------------
 src/solvers/rgd.jl            | 34 ++--------------------------------
 5 files changed, 11 insertions(+), 76 deletions(-)

diff --git a/src/api/approx.jl b/src/api/approx.jl
index 7765c18..3aa4866 100644
--- a/src/api/approx.jl
+++ b/src/api/approx.jl
@@ -324,8 +324,7 @@ function approx(
     return _approx_tucker_rank(approx_dispatch(dispatch), base, r, target; kwargs...)
 end
 
-_reject_generic_rank_dispatch(::AutoApproxDispatch) = nothing
-_reject_generic_rank_dispatch(::GenericApproxDispatch) = nothing
+_reject_generic_rank_dispatch(::Union{AutoApproxDispatch,GenericApproxDispatch}) = nothing
 
 function _reject_generic_rank_dispatch(::CPDApproxDispatch)
     throw(ArgumentError("approx(...; dispatch=:cpd) requires Manifolds.Segre inputs."))
diff --git a/src/solvers/lbfgs.jl b/src/solvers/lbfgs.jl
index fa0009c..626f48f 100644
--- a/src/solvers/lbfgs.jl
+++ b/src/solvers/lbfgs.jl
@@ -172,22 +172,7 @@ function solve_lbfgs(
             uses_nonpositive_curvature_behavior = false,
         ),
     )
-    return isnothing(normA2) ?
-           _solver_stats(
-        model_cost,
-        model_grad_local,
-        M,
-        p_opt,
-        state,
-        nothing;
-        tol_T = T(tol),
-        maxiter,
-        solver = :lbfgs,
-        tiny_grad_tol = tol_g_raw,
-        solver_info,
-        use_state_gradient = !uses_relative_objective,
-    ) :
-           _solver_stats(
+    return _solver_stats(
         model_cost,
         model_grad_local,
         M,
diff --git a/src/solvers/manopt_helpers.jl b/src/solvers/manopt_helpers.jl
index b997a24..ac46e85 100644
--- a/src/solvers/manopt_helpers.jl
+++ b/src/solvers/manopt_helpers.jl
@@ -331,10 +331,10 @@ end
     return isdefined(Manopt, :stopped_at) ? :stopped_at : :stop_at_iteration_fallback
 end
 
-# Convert a squared residual value to the public relative-error diagnostic.
-_solver_rel_error(cost_for_error, ::Nothing, ::Type) = sqrt(cost_for_error)
-_solver_rel_error(cost_for_error, normA2::Real, ::Type{T}) where {T} =
-    _relative_error_frob_sq(cost_for_error, T(normA2))
+function _solver_rel_error(cost_for_error, normA2::Union{Nothing,Real}, ::Type{T}) where {T}
+    return isnothing(normA2) ? sqrt(cost_for_error) :
+           _relative_error_frob_sq(cost_for_error, T(normA2))
+end
 
 # Build common solver result stats, with optional ||A||^2 for relative scaling.
 function _solver_stats(
@@ -389,14 +389,10 @@ end
 # Collect only active Manopt debug callbacks, dropping omitted hooks.
 _solver_debug_callbacks(callbacks...) = Any[cb for cb in callbacks if !isnothing(cb)]
 
-# Allow callers to pass `nothing` (e.g., when verbose/debug is omitted)
-# Build debug actions when the verbose flag was omitted.
-_solver_debug_actions(::Nothing, callbacks...) = _solver_debug_callbacks(callbacks...)
-
 # Build Manopt debug actions and attach TensorKitchen callback hooks.
-function _solver_debug_actions(verbose::Bool, callbacks...)
+function _solver_debug_actions(verbose::Union{Nothing,Bool}, callbacks...)
     callback_actions = _solver_debug_callbacks(callbacks...)
-    if verbose
+    if verbose === true
         io = _SOLVER_DEBUG_SINK
         init_group = Manopt.DebugGroup([
             Manopt.DebugDivider("Initial "; io, at_init = true),
diff --git a/src/solvers/rcg.jl b/src/solvers/rcg.jl
index a1ac2d9..f667379 100644
--- a/src/solvers/rcg.jl
+++ b/src/solvers/rcg.jl
@@ -181,22 +181,7 @@ function solve_rcg(
     if !return_stats
         return p_opt
     end
-    return isnothing(normA2) ?
-           _solver_stats(
-        model_cost,
-        model_grad_local,
-        M,
-        p_opt,
-        state,
-        nothing;
-        tol_T = T(tol),
-        maxiter,
-        solver = :rcg,
-        tiny_grad_tol = tol_g_raw,
-        solver_info,
-        use_state_gradient = !uses_relative_objective,
-    ) :
-           _solver_stats(
+    return _solver_stats(
         model_cost,
         model_grad_local,
         M,
diff --git a/src/solvers/rgd.jl b/src/solvers/rgd.jl
index be850a6..2490ae8 100644
--- a/src/solvers/rgd.jl
+++ b/src/solvers/rgd.jl
@@ -131,22 +131,7 @@ function solve_rgd(
     if !return_stats
         return p_opt
     end
-    return isnothing(normA2) ?
-           _solver_stats(
-        model_cost,
-        model_grad_local,
-        M,
-        p_opt,
-        state,
-        nothing;
-        tol_T = T(tol),
-        maxiter,
-        solver = :rgd,
-        tiny_grad_tol = tol_g_raw,
-        solver_info,
-        use_state_gradient = !uses_relative_objective,
-    ) :
-           _solver_stats(
+    return _solver_stats(
         model_cost,
         model_grad_local,
         M,
@@ -244,22 +229,7 @@ function solve_rgd_fixed(
     if !return_stats
         return p_opt
     end
-    return isnothing(normA2) ?
-           _solver_stats(
-        model_cost,
-        model_grad_local,
-        M,
-        p_opt,
-        state,
-        nothing;
-        tol_T = T(tol),
-        maxiter,
-        solver = :rgd_fixed,
-        tiny_grad_tol = tiny_grad_tol,
-        solver_info,
-        use_state_gradient = !uses_relative_objective,
-    ) :
-           _solver_stats(
+    return _solver_stats(
         model_cost,
         model_grad_local,
         M,

From 861e9e97dffdcaab49bc82ab53caad752144258f Mon Sep 17 00:00:00 2001
From: PBrdng <pbrdng@gmail.com>
Date: Wed, 17 Jun 2026 20:55:56 +0200
Subject: [PATCH 09/18] make normalized_objective=true be the default

---
 src/api/cpd.jl                | 2 --
 src/solvers/abstract.jl       | 4 ++--
 src/solvers/lbfgs.jl          | 4 ++--
 src/solvers/rcg.jl            | 4 ++--
 src/solvers/rgd.jl            | 8 ++++----
 src/solvers/solve_dispatch.jl | 2 +-
 6 files changed, 11 insertions(+), 13 deletions(-)

diff --git a/src/api/cpd.jl b/src/api/cpd.jl
index 2d19de5..c06fcf3 100644
--- a/src/api/cpd.jl
+++ b/src/api/cpd.jl
@@ -862,8 +862,6 @@ function _run_cpd_solver(
         refinement_verbose = verbose,
         vector_transport_method,
         grad_tol = _cpd_manifold_grad_tol(model, solver, tol),
-        normalized_objective = solver isa
-                               Union{RGDSolver,RGDFixedSolver,RCGSolver,LBFGSSolver},
         iteration_callbacks,
         kwargs...,
     )
diff --git a/src/solvers/abstract.jl b/src/solvers/abstract.jl
index dac468d..52afb35 100644
--- a/src/solvers/abstract.jl
+++ b/src/solvers/abstract.jl
@@ -225,7 +225,7 @@ function solve(
     return_stats::Bool = false,
     vector_transport_method::Union{ManifoldsBase.AbstractVectorTransportMethod,Nothing} = nothing,
     grad_tol = nothing,
-    normalized_objective::Bool = false,
+    normalized_objective::Bool = true,
     iteration_callbacks = (),
 ) where {T<:AbstractFloat}
     setup = _prepare_solver_problem(model; init, p0, gradient_mode, verbose)
@@ -268,7 +268,7 @@ function solve(
     return_stats::Bool = false,
     vector_transport_method::Union{ManifoldsBase.AbstractVectorTransportMethod,Nothing} = nothing,
     grad_tol = nothing,
-    normalized_objective::Bool = false,
+    normalized_objective::Bool = true,
     iteration_callbacks = (),
 ) where {T<:AbstractFloat}
     setup = _prepare_solver_problem(model; init, p0, gradient_mode)
diff --git a/src/solvers/lbfgs.jl b/src/solvers/lbfgs.jl
index 626f48f..6c3badf 100644
--- a/src/solvers/lbfgs.jl
+++ b/src/solvers/lbfgs.jl
@@ -75,7 +75,7 @@ function solve_lbfgs(
     linesearch::Symbol = :wolfe,
     preconditioner = nothing,
     grad_tol = nothing,
-    normalized_objective::Bool = false,
+    normalized_objective::Bool = true,
 )
     p0_local = _solver_point(M, p0)
     T = _scalar_eltype(p0_local)
@@ -200,7 +200,7 @@ function run_second_order_solver(
     diagnostics_recorder,
     iteration_callbacks,
     grad_tol = nothing,
-    normalized_objective::Bool = false,
+    normalized_objective::Bool = true,
 )
     return solve_lbfgs(
         setup.model_cost,
diff --git a/src/solvers/rcg.jl b/src/solvers/rcg.jl
index f667379..6d198a1 100644
--- a/src/solvers/rcg.jl
+++ b/src/solvers/rcg.jl
@@ -106,7 +106,7 @@ function solve_rcg(
     diagnostics_recorder = nothing,
     iteration_callbacks = (),
     grad_tol = nothing,
-    normalized_objective::Bool = false,
+    normalized_objective::Bool = true,
 )
     p0_local = _solver_point(M, p0)
     T = _scalar_eltype(p0_local)
@@ -223,7 +223,7 @@ function run_first_order_solver(
     diagnostics_recorder,
     iteration_callbacks,
     grad_tol = nothing,
-    normalized_objective::Bool = false,
+    normalized_objective::Bool = true,
 )
     return solve_rcg(
         setup.model_cost,
diff --git a/src/solvers/rgd.jl b/src/solvers/rgd.jl
index 2490ae8..14b7af7 100644
--- a/src/solvers/rgd.jl
+++ b/src/solvers/rgd.jl
@@ -19,7 +19,7 @@ function solve_rgd(
     diagnostics_recorder = nothing,
     iteration_callbacks = (),
     grad_tol = nothing,
-    normalized_objective::Bool = false,
+    normalized_objective::Bool = true,
 )
     p0_local = _solver_point(M, p0)
     T = _scalar_eltype(p0_local)
@@ -163,7 +163,7 @@ function solve_rgd_fixed(
     diagnostics_recorder = nothing,
     iteration_callbacks = (),
     grad_tol = nothing,
-    normalized_objective::Bool = false,
+    normalized_objective::Bool = true,
 )
     p0_local = _solver_point(M, p0)
     T = _scalar_eltype(p0_local)
@@ -274,7 +274,7 @@ function run_first_order_solver(
     diagnostics_recorder,
     iteration_callbacks,
     grad_tol = nothing,
-    normalized_objective::Bool = false,
+    normalized_objective::Bool = true,
 )
     return solve_rgd(
         setup.model_cost,
@@ -329,7 +329,7 @@ function run_first_order_solver(
     diagnostics_recorder,
     iteration_callbacks,
     grad_tol = nothing,
-    normalized_objective::Bool = false,
+    normalized_objective::Bool = true,
 )
     return solve_rgd_fixed(
         setup.model_cost,
diff --git a/src/solvers/solve_dispatch.jl b/src/solvers/solve_dispatch.jl
index fcce422..0f132bd 100644
--- a/src/solvers/solve_dispatch.jl
+++ b/src/solvers/solve_dispatch.jl
@@ -65,7 +65,7 @@ function _solve_with_solver(
     verbose::Bool,
     vector_transport_method::Union{ManifoldsBase.AbstractVectorTransportMethod,Nothing} = nothing,
     grad_tol = nothing,
-    normalized_objective::Bool = false,
+    normalized_objective::Bool = true,
     iteration_callbacks = (),
     kwargs...,
 )

From abe7003b3c9e2ddabf4ef5cb1008b8edbdb44bfc Mon Sep 17 00:00:00 2001
From: PBrdng <pbrdng@gmail.com>
Date: Wed, 17 Jun 2026 21:15:51 +0200
Subject: [PATCH 10/18] wire in relative norm natively

---
 src/solvers/lbfgs.jl          | 10 ++++----
 src/solvers/manopt_helpers.jl | 27 +++++++++++----------
 src/solvers/rcg.jl            | 10 ++++----
 src/solvers/rgd.jl            | 21 +++++++----------
 test/basic_tests.jl           | 44 +++++++++++++++++++++++++++++++++--
 5 files changed, 73 insertions(+), 39 deletions(-)

diff --git a/src/solvers/lbfgs.jl b/src/solvers/lbfgs.jl
index 6c3badf..24ab188 100644
--- a/src/solvers/lbfgs.jl
+++ b/src/solvers/lbfgs.jl
@@ -92,7 +92,6 @@ function solve_lbfgs(
         vector_transport_method
     grad_stop_tol = isnothing(grad_tol) ? T(tol) : T(grad_tol)
     tol_g = _dual_stop_grad_tol(T, tol, grad_tol)
-    tol_g_raw = uses_relative_objective ? tol_g * objective_scale : tol_g
     dual_stop = StopWhenCostRelChangeAndGradientLess(T(tol), tol_g)
     stopping = StopWhenAny(
         StopAfterIteration(maxiter),
@@ -116,7 +115,6 @@ function solve_lbfgs(
         solver_cost,
         solver_grad,
         M;
-        normA2 = uses_relative_objective ? normA2 : nothing,
         diagnostics_recorder,
     )
 
@@ -173,8 +171,8 @@ function solve_lbfgs(
         ),
     )
     return _solver_stats(
-        model_cost,
-        model_grad_local,
+        solver_cost,
+        solver_grad,
         M,
         p_opt,
         state,
@@ -182,9 +180,9 @@ function solve_lbfgs(
         tol_T = T(tol),
         maxiter,
         solver = :lbfgs,
-        tiny_grad_tol = tol_g_raw,
+        tiny_grad_tol = tol_g,
         solver_info,
-        use_state_gradient = !uses_relative_objective,
+        normalized_objective = uses_relative_objective,
     )
 end
 
diff --git a/src/solvers/manopt_helpers.jl b/src/solvers/manopt_helpers.jl
index ac46e85..559cccc 100644
--- a/src/solvers/manopt_helpers.jl
+++ b/src/solvers/manopt_helpers.jl
@@ -296,6 +296,7 @@ function _scale_solver_tangent(x, scale::Real)
         return scale * x
     end
 end
+# Build the Manopt objective as squared residual cost, optionally divided by ||target||^2.
 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
@@ -331,12 +332,19 @@ end
     return isdefined(Manopt, :stopped_at) ? :stopped_at : :stop_at_iteration_fallback
 end
 
-function _solver_rel_error(cost_for_error, normA2::Union{Nothing,Real}, ::Type{T}) where {T}
-    return isnothing(normA2) ? sqrt(cost_for_error) :
-           _relative_error_frob_sq(cost_for_error, T(normA2))
+function _solver_rel_error(
+    final_cost,
+    normA2::Union{Nothing,Real},
+    normalized_objective::Bool,
+    ::Type{T},
+) where {T}
+    cost_for_error = max(T(0), T(2) * final_cost)
+    normalized_objective && return sqrt(cost_for_error)
+    return isnothing(normA2) || normA2 <= 0 ? sqrt(cost_for_error) :
+           sqrt(cost_for_error / T(normA2))
 end
 
-# Build common solver result stats, with optional ||A||^2 for relative scaling.
+# Build common solver result stats
 function _solver_stats(
     model_cost,
     model_grad,
@@ -350,11 +358,11 @@ function _solver_stats(
     tiny_grad_tol = nothing,
     solver_info = (;),
     use_state_gradient::Bool = true,
+    normalized_objective::Bool = false,
 )
     T = typeof(tol_T)
     final_cost = model_cost(M, p_opt)
-    cost_for_error = max(T(0), T(2) * final_cost)
-    rel_error = _solver_rel_error(cost_for_error, normA2, T)
+    rel_error = _solver_rel_error(final_cost, normA2, normalized_objective, T)
     grad_state = use_state_gradient ? _solver_gradient(state) : nothing
     grad_from_state = !isnothing(grad_state)
     grad_final =
@@ -423,21 +431,16 @@ function _solver_progress_callback(
     model_cost,
     model_grad,
     M;
-    normA2 = nothing,
     diagnostics_recorder = nothing,
 )
     progress isa NoMethodProgress && return nothing
-    has_relative_scale = !isnothing(normA2) && normA2 > 0
-    target_norm = has_relative_scale ? sqrt(normA2) : nothing
     return function (problem, state, k)
         k <= 0 && return nothing
         p = get_iterate(state)
         c = model_cost(M, p)
         g = model_grad(M, p)
         gnorm = norm(M, p, g)
-        c_display = has_relative_scale ? sqrt(max(2 * c, zero(c))) : c
-        gnorm_display = has_relative_scale ? gnorm * target_norm : gnorm
-        showvalues = Any[("Iter", k), ("Cost", c_display), ("Grad norm", gnorm_display)]
+        showvalues = Any[("Iter", k), ("Cost", c), ("Grad norm", gnorm)]
         if !isnothing(diagnostics_recorder)
             step = diagnostics_recorder.accepted_stepsize_history
             trials = diagnostics_recorder.line_search_trial_history
diff --git a/src/solvers/rcg.jl b/src/solvers/rcg.jl
index 6d198a1..afed20e 100644
--- a/src/solvers/rcg.jl
+++ b/src/solvers/rcg.jl
@@ -123,7 +123,6 @@ function solve_rcg(
         vector_transport_method
     grad_stop_tol = isnothing(grad_tol) ? T(tol) : T(grad_tol)
     tol_g = _dual_stop_grad_tol(T, tol, grad_tol)
-    tol_g_raw = uses_relative_objective ? tol_g * objective_scale : tol_g
     dual_stop = StopWhenCostRelChangeAndGradientLess(T(tol), tol_g)
 
     stopping = StopWhenAny(
@@ -143,7 +142,6 @@ function solve_rcg(
         solver_cost,
         solver_grad,
         M;
-        normA2 = uses_relative_objective ? normA2 : nothing,
         diagnostics_recorder,
     )
 
@@ -182,8 +180,8 @@ function solve_rcg(
         return p_opt
     end
     return _solver_stats(
-        model_cost,
-        model_grad_local,
+        solver_cost,
+        solver_grad,
         M,
         p_opt,
         state,
@@ -191,9 +189,9 @@ function solve_rcg(
         tol_T = T(tol),
         maxiter,
         solver = :rcg,
-        tiny_grad_tol = tol_g_raw,
+        tiny_grad_tol = tol_g,
         solver_info,
-        use_state_gradient = !uses_relative_objective,
+        normalized_objective = uses_relative_objective,
     )
 end
 
diff --git a/src/solvers/rgd.jl b/src/solvers/rgd.jl
index 14b7af7..2503f61 100644
--- a/src/solvers/rgd.jl
+++ b/src/solvers/rgd.jl
@@ -34,7 +34,6 @@ function solve_rgd(
     armijo_alpha_min = T(1e-8) * objective_scale
     grad_stop_tol = isnothing(grad_tol) ? T(tol) : T(grad_tol)
     tol_g = _dual_stop_grad_tol(T, tol, grad_tol)
-    tol_g_raw = uses_relative_objective ? tol_g * objective_scale : tol_g
     dual_stop = StopWhenCostRelChangeAndGradientLess(T(tol), tol_g)
 
     stopping = StopWhenAny(
@@ -91,7 +90,6 @@ function solve_rgd(
         solver_cost,
         solver_grad,
         M;
-        normA2 = uses_relative_objective ? normA2 : nothing,
         diagnostics_recorder,
     )
 
@@ -132,8 +130,8 @@ function solve_rgd(
         return p_opt
     end
     return _solver_stats(
-        model_cost,
-        model_grad_local,
+        solver_cost,
+        solver_grad,
         M,
         p_opt,
         state,
@@ -141,9 +139,9 @@ function solve_rgd(
         tol_T = T(tol),
         maxiter,
         solver = :rgd,
-        tiny_grad_tol = tol_g_raw,
+        tiny_grad_tol = tol_g,
         solver_info,
-        use_state_gradient = !uses_relative_objective,
+        normalized_objective = uses_relative_objective,
     )
 end
 
@@ -175,9 +173,7 @@ function solve_rgd_fixed(
         _relative_solver_functions(model_cost, model_grad_local, objective_scale)
     retraction_method = _solver_retraction_method(M, p0_local)
     grad_stop_tol = isnothing(grad_tol) ? T(tol) : T(grad_tol)
-    tiny_grad_tol =
-        isnothing(grad_tol) ? T(1e-5) :
-        (uses_relative_objective ? T(grad_tol) * objective_scale : T(grad_tol))
+    tiny_grad_tol = isnothing(grad_tol) ? T(1e-5) : T(grad_tol)
     stopping =
         StopWhenAny(StopAfterIteration(maxiter), StopWhenGradientNormLess(grad_stop_tol))
     progress =
@@ -192,7 +188,6 @@ function solve_rgd_fixed(
         solver_cost,
         solver_grad,
         M;
-        normA2 = uses_relative_objective ? normA2 : nothing,
         diagnostics_recorder,
     )
     state = gradient_descent(
@@ -230,8 +225,8 @@ function solve_rgd_fixed(
         return p_opt
     end
     return _solver_stats(
-        model_cost,
-        model_grad_local,
+        solver_cost,
+        solver_grad,
         M,
         p_opt,
         state,
@@ -241,7 +236,7 @@ function solve_rgd_fixed(
         solver = :rgd_fixed,
         tiny_grad_tol = tiny_grad_tol,
         solver_info,
-        use_state_gradient = !uses_relative_objective,
+        normalized_objective = uses_relative_objective,
     )
 end
 
diff --git a/test/basic_tests.jl b/test/basic_tests.jl
index 041a458..6a818c6 100644
--- a/test/basic_tests.jl
+++ b/test/basic_tests.jl
@@ -136,6 +136,45 @@ end
     @test out.solver_info.gradient_evaluations >= 1
 end
 
+@testset "Manopt normalized_objective controls objective units" begin
+    A = randn(6, 5, 4)
+    r = 2
+    model = JoinModel(A, r; geometry = :canonical)
+    p0 = TensorKitchen.initial_point(model, TuckerInit(); verbose = false)
+    target_norm = norm(A)
+
+    rel_out = solve(
+        RGDFixedSolver(0.0),
+        model;
+        p0,
+        maxiter = 1,
+        tol = 0.0,
+        verbose = false,
+        return_stats = true,
+        normalized_objective = true,
+    )
+    abs_out = solve(
+        RGDFixedSolver(0.0),
+        model;
+        p0,
+        maxiter = 1,
+        tol = 0.0,
+        verbose = false,
+        return_stats = true,
+        normalized_objective = false,
+    )
+
+    @test isapprox(2 * rel_out.cost, rel_out.rel_error^2; rtol = 1e-12, atol = 1e-12)
+    @test isapprox(abs_out.rel_error, rel_out.rel_error; rtol = 1e-12, atol = 1e-12)
+    @test isapprox(abs_out.cost, rel_out.cost * target_norm^2; rtol = 1e-12, atol = 1e-12)
+    @test isapprox(
+        abs_out.grad_norm,
+        rel_out.grad_norm * target_norm^2;
+        rtol = 1e-10,
+        atol = 1e-10,
+    )
+end
+
 # =========================================================================
 # cpd/cp_rank.jl (cost/egrad functions)
 # =========================================================================
@@ -433,6 +472,7 @@ end
         rel = norm(A) > 0 ? norm(X .- A) / norm(A) : norm(X .- A)
         cost, rel
     end
+    expected_solver_cost(solver, cost, rel) = solver == :als ? cost : 0.5 * rel^2
     public_columns_unit(res) = all(
         isapprox(norm(TensorKitchen.factors(res)[m][:, k]), 1; atol = 1e-8, rtol = 1e-8) for m in eachindex(TensorKitchen.factors(res)) for
         k in eachindex(TensorKitchen.weights(res))
@@ -443,7 +483,7 @@ end
     for solver in (:rgd, :rcg)
         res = cpd(A1, 1; solver = solver, nonnegative = true, maxiter = 4, verbose = false)
         cost, rel = explicit_stats(A1, res)
-        @test res.cost ≈ cost atol = 1e-8 rtol = 1e-8
+        @test res.cost ≈ expected_solver_cost(solver, cost, rel) atol = 1e-8 rtol = 1e-8
         @test res.rel_error ≈ rel atol = 1e-8 rtol = 1e-8
         @test public_columns_unit(res)
     end
@@ -453,7 +493,7 @@ end
     for solver in (:als, :rgd, :rcg)
         res = cpd(Ar, 3; solver = solver, nonnegative = true, maxiter = 4, verbose = false)
         cost, rel = explicit_stats(Ar, res)
-        @test res.cost ≈ cost atol = 1e-8 rtol = 1e-8
+        @test res.cost ≈ expected_solver_cost(solver, cost, rel) atol = 1e-8 rtol = 1e-8
         @test res.rel_error ≈ rel atol = 1e-8 rtol = 1e-8
         @test public_columns_unit(res)
     end

From 6a5aa55530b04916653cd69697c6bba4a383f573 Mon Sep 17 00:00:00 2001
From: PBrdng <pbrdng@gmail.com>
Date: Wed, 17 Jun 2026 21:16:32 +0200
Subject: [PATCH 11/18] compactify progress meter

---
 src/solvers/manopt_helpers.jl | 6 ------
 1 file changed, 6 deletions(-)

diff --git a/src/solvers/manopt_helpers.jl b/src/solvers/manopt_helpers.jl
index 559cccc..1ae3d2e 100644
--- a/src/solvers/manopt_helpers.jl
+++ b/src/solvers/manopt_helpers.jl
@@ -441,12 +441,6 @@ function _solver_progress_callback(
         g = model_grad(M, p)
         gnorm = norm(M, p, g)
         showvalues = Any[("Iter", k), ("Cost", c), ("Grad norm", gnorm)]
-        if !isnothing(diagnostics_recorder)
-            step = diagnostics_recorder.accepted_stepsize_history
-            trials = diagnostics_recorder.line_search_trial_history
-            !isempty(step) && push!(showvalues, ("Accepted α", step[end]))
-            !isempty(trials) && push!(showvalues, ("Line-search trials", trials[end]))
-        end
         update_progress!(progress, k; showvalues)
         return nothing
     end

From 7fb5557937b6c0d8e77fa04a99d4c13f70546bbb Mon Sep 17 00:00:00 2001
From: PBrdng <pbrdng@gmail.com>
Date: Wed, 17 Jun 2026 21:24:36 +0200
Subject: [PATCH 12/18] print progress meter correctly

---
 src/core/progress.jl | 20 +++++++++++++++++---
 1 file changed, 17 insertions(+), 3 deletions(-)

diff --git a/src/core/progress.jl b/src/core/progress.jl
index 30c514a..e77f00c 100644
--- a/src/core/progress.jl
+++ b/src/core/progress.jl
@@ -176,6 +176,20 @@ end
     p.was_rendered = true
     return nothing
 end
+@inline _meter_was_printed(meter) = getproperty(getproperty(meter, :core), :printed)
+@inline _sync_rendered!(::NoMethodProgress) = nothing
+function _sync_rendered!(p::FamilyProgress)
+    meter = _meter(p)
+    if !isnothing(meter) && _meter_was_printed(meter)
+        _mark_rendered!(p)
+    end
+    return nothing
+end
+function _sync_tracker_rendered!(tracker::PhaseProgress)
+    _sync_rendered!(tracker.initialization)
+    _sync_rendered!(tracker.refinement)
+    return nothing
+end
 
 @inline function _force_visible_phase_finish(tracker::PhaseProgress, progress)
     return tracker.phase == :refinement &&
@@ -226,9 +240,7 @@ function update_progress!(
             Any[("Method", _method_name(progress)); showvalues]
         end
         PM.update!(meter, current; showvalues = showvalues_with_method, force)
-        if force || current < meter.n || _was_rendered(progress)
-            _mark_rendered!(progress)
-        end
+        _sync_rendered!(progress)
     end
     return nothing
 end
@@ -244,6 +256,7 @@ function finish_progress!(
     progress isa NoMethodProgress && return nothing
     meter = _meter(progress)
     isnothing(meter) && return nothing
+    _sync_tracker_rendered!(tracker)
 
     showvalues_with_method = if isnothing(showvalues)
         Any[("Method", _method_name(progress))]
@@ -273,6 +286,7 @@ function finish_progress!(
     isnothing(meter) && return nothing
     tracker = _current_phase_tracker()
     if tracker isa PhaseProgress
+        _sync_tracker_rendered!(tracker)
         set_phase!(tracker, progress.phase)
         if active_progress(tracker) === progress
             return finish_progress!(tracker; current, showvalues)

From 9efdcde9882a73ae85af5f969d302f72cd1d0f0c Mon Sep 17 00:00:00 2001
From: PBrdng <pbrdng@gmail.com>
Date: Wed, 17 Jun 2026 21:24:52 +0200
Subject: [PATCH 13/18] run JuliaFormatter

---
 src/core/progress.jl | 8 +-------
 1 file changed, 1 insertion(+), 7 deletions(-)

diff --git a/src/core/progress.jl b/src/core/progress.jl
index e77f00c..4136f1d 100644
--- a/src/core/progress.jl
+++ b/src/core/progress.jl
@@ -200,13 +200,7 @@ end
 function _render_unrendered_completion!(meter, showvalues)
     # ProgressMeter does not render a meter that reaches 100% before its first
     # visible update. Give it one display-only step before completion.
-    PM.update!(
-        meter,
-        meter.n;
-        showvalues,
-        force = true,
-        max_steps = meter.n + 1,
-    )
+    PM.update!(meter, meter.n; showvalues, force = true, max_steps = meter.n + 1)
     return nothing
 end
 

From 98c2d1e2c829741a338de180a46663ac77f555b8 Mon Sep 17 00:00:00 2001
From: Se Eun Choi <seeun.choi@uni-osnabrueck.de>
Date: Wed, 17 Jun 2026 19:24:43 +0000
Subject: [PATCH 14/18] format code with JuliaFormatter

---
 Project.toml             | 2 ++
 src/core/tensor_ops.jl   | 4 ++--
 src/cpd/core/cpd_init.jl | 2 +-
 src/tucker/sthosvd.jl    | 4 ++--
 4 files changed, 7 insertions(+), 5 deletions(-)

diff --git a/Project.toml b/Project.toml
index bd4414c..cd61afc 100644
--- a/Project.toml
+++ b/Project.toml
@@ -3,6 +3,7 @@ uuid = "3630a16b-0f2f-4d88-afbf-c7d59eccf553"
 version = "0.1.0"
 
 [deps]
+JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
 ManifoldsBase = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb"
@@ -13,6 +14,7 @@ RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
 
 [compat]
+JuliaFormatter = "2.8.5"
 Manifolds = "0.11.20"
 ManifoldsBase = "2.3.5"
 Manopt = "0.5.37"
diff --git a/src/core/tensor_ops.jl b/src/core/tensor_ops.jl
index 5ab2f89..e5c1427 100644
--- a/src/core/tensor_ops.jl
+++ b/src/core/tensor_ops.jl
@@ -470,8 +470,8 @@ gradU_column_cp(
 
 cp_reconstruction_norm2(components::Vector{RankOneTensor{T}}) where {T<:AbstractFloat} =
     sum(
-        cross_component(components[i], components[j]) for i in eachindex(components),
-        j in eachindex(components)
+        cross_component(components[i], components[j]) for
+        i in eachindex(components), j in eachindex(components)
     )
 
 function cp_inner_AX(
diff --git a/src/cpd/core/cpd_init.jl b/src/cpd/core/cpd_init.jl
index b6964e7..bff19c2 100644
--- a/src/cpd/core/cpd_init.jl
+++ b/src/cpd/core/cpd_init.jl
@@ -44,7 +44,7 @@ function _cp_core_diag_init(core::AbstractArray{T,N}, r::Int) where {T<:Abstract
             Um[k, k] = one(T)
         end
         if rm > 0 && r > n_eye
-            Um[:, n_eye+1:r] .= random_unit_matrix(rm, r - n_eye, T)
+            Um[:, (n_eye+1):r] .= random_unit_matrix(rm, r - n_eye, T)
         end
         U0[m] = Um
     end
diff --git a/src/tucker/sthosvd.jl b/src/tucker/sthosvd.jl
index 2b0e575..cf69651 100644
--- a/src/tucker/sthosvd.jl
+++ b/src/tucker/sthosvd.jl
@@ -213,7 +213,7 @@ function sthosvd(
         rk = max(rk, 1)  # keep at least rank 1
 
         if verbose
-            discarded = rk < length(sigma) ? sqrt(sum(sigma[rk+1:end] .^ 2)) : 0.0
+            discarded = rk < length(sigma) ? sqrt(sum(sigma[(rk+1):end] .^ 2)) : 0.0
             update_progress!(
                 progress,
                 step;
@@ -322,7 +322,7 @@ function error_bound(td::TuckerResult{T,N}) where {T,N}
         rk = size(td.core, k)
         sigma = td.singular_values[k]
         if rk < length(sigma)
-            sq_error += sum(sigma[rk+1:end] .^ 2)
+            sq_error += sum(sigma[(rk+1):end] .^ 2)
         end
     end
     return sqrt(sq_error)

From a5a65cd176a4e568e179d2235241b8ac69705606 Mon Sep 17 00:00:00 2001
From: PBrdng <pbrdng@gmail.com>
Date: Wed, 17 Jun 2026 21:43:52 +0200
Subject: [PATCH 15/18] unify code

---
 src/solvers/lbfgs.jl          | 113 ++++++++------------
 src/solvers/manopt_helpers.jl | 130 ++++++++++++++++++++++
 src/solvers/rcg.jl            |  93 +++++++---------
 src/solvers/rgd.jl            | 196 ++++++++++++++--------------------
 4 files changed, 292 insertions(+), 240 deletions(-)

diff --git a/src/solvers/lbfgs.jl b/src/solvers/lbfgs.jl
index 24ab188..1d4c17d 100644
--- a/src/solvers/lbfgs.jl
+++ b/src/solvers/lbfgs.jl
@@ -77,51 +77,49 @@ function solve_lbfgs(
     grad_tol = nothing,
     normalized_objective::Bool = true,
 )
-    p0_local = _solver_point(M, p0)
-    T = _scalar_eltype(p0_local)
-    model_grad_raw = isnothing(model_grad) ? grad(model_egrad) : model_grad
-    model_grad_local = _layout_adapt_gradient(model_grad_raw)
-    objective_scale =
-        normalized_objective && !isnothing(normA2) && normA2 > 0 ? T(normA2) : one(T)
-    solver_cost, solver_grad, uses_relative_objective =
-        _relative_solver_functions(model_cost, model_grad_local, objective_scale)
+    setup = _prepare_manopt_solver_functions(
+        model_cost,
+        model_egrad,
+        M,
+        p0;
+        normA2,
+        model_grad,
+        tol,
+        grad_tol,
+        normalized_objective,
+    )
+    p0_local = setup.p0
+    T = setup.T
     retraction_method = _solver_retraction_method(M, p0_local)
     transport =
         isnothing(vector_transport_method) ?
         _default_vector_transport_method(M, p0_local, retraction_method) :
         vector_transport_method
-    grad_stop_tol = isnothing(grad_tol) ? T(tol) : T(grad_tol)
-    tol_g = _dual_stop_grad_tol(T, tol, grad_tol)
+    tol_g = setup.dual_grad_tol
     dual_stop = StopWhenCostRelChangeAndGradientLess(T(tol), tol_g)
-    stopping = StopWhenAny(
-        StopAfterIteration(maxiter),
-        StopWhenGradientNormLess(grad_stop_tol),
-        dual_stop,
-    )
-    progress =
-        maxiter > 0 ?
-        make_manopt_family_progress(
-            maxiter;
+    stopping = _manopt_stopping(maxiter, setup.grad_stop_tol, dual_stop)
+    callbacks = _manopt_callbacks(
+        n -> make_manopt_family_progress(
+            n;
             enabled = verbose,
             phase = :refinement,
             method = "L-BFGS",
             dt = 0.2,
-        ) : NoMethodProgress()
-    diagnostics_callback =
-        isnothing(diagnostics_recorder) ? nothing :
-        _solver_diagnostics_callback(diagnostics_recorder)
-    progress_callback = _solver_progress_callback(
-        progress,
-        solver_cost,
-        solver_grad,
+        ),
+        maxiter,
+        verbose,
+        setup.solver_cost,
+        setup.solver_grad,
         M;
         diagnostics_recorder,
+        post_step_callback,
+        iteration_callbacks,
     )
 
     state = Manopt.quasi_Newton(
         M,
-        solver_cost,
-        solver_grad,
+        setup.solver_cost,
+        setup.solver_grad,
         p0_local;
         cautious_update = cautious_update,
         direction_update = Manopt.InverseBFGS(),
@@ -132,35 +130,28 @@ function solve_lbfgs(
         vector_transport_method = transport,
         stepsize = _lbfgs_linesearch(linesearch),
         stopping_criterion = stopping,
-        debug = _solver_debug_actions(
-            verbose,
-            post_step_callback,
-            diagnostics_callback,
-            progress_callback,
-            iteration_callbacks...,
-        ),
+        debug = callbacks.debug_actions,
         count = [:Cost, :Gradient],
         return_state = true,
     )
 
-    p_opt = _tk_get_solver_result(state)
-    iterations_done = _solver_iterations(state, maxiter)
-    if verbose
-        finish_progress!(
-            progress;
-            current = iterations_done,
-            showvalues = Any[("Status", "Finished"), ("Iterations", iterations_done)],
-        )
-    end
-    solver_info =
-        isnothing(diagnostics_recorder) ? (;) :
-        _solver_info(diagnostics_recorder, iterations_done)
-    if !return_stats
-        return p_opt
-    end
-    solver_info = merge(
-        solver_info,
-        (
+    return _manopt_finish_result(
+        _tk_get_solver_result(state),
+        state,
+        callbacks.progress,
+        diagnostics_recorder,
+        setup.solver_cost,
+        setup.solver_grad,
+        M,
+        normA2;
+        tol_T = T(tol),
+        maxiter,
+        solver = :lbfgs,
+        tiny_grad_tol = tol_g,
+        return_stats,
+        verbose,
+        normalized_objective = setup.uses_relative_objective,
+        solver_info_extra = (
             memory_size = memory_size,
             cautious_update = cautious_update,
             initial_scale = initial_scale,
@@ -170,20 +161,6 @@ function solve_lbfgs(
             uses_nonpositive_curvature_behavior = false,
         ),
     )
-    return _solver_stats(
-        solver_cost,
-        solver_grad,
-        M,
-        p_opt,
-        state,
-        normA2;
-        tol_T = T(tol),
-        maxiter,
-        solver = :lbfgs,
-        tiny_grad_tol = tol_g,
-        solver_info,
-        normalized_objective = uses_relative_objective,
-    )
 end
 
 function run_second_order_solver(
diff --git a/src/solvers/manopt_helpers.jl b/src/solvers/manopt_helpers.jl
index 1ae3d2e..67b8ac8 100644
--- a/src/solvers/manopt_helpers.jl
+++ b/src/solvers/manopt_helpers.jl
@@ -308,6 +308,136 @@ function _relative_solver_functions(model_cost, model_grad, scale::Real)
     )
 end
 
+# Prepare the shared Manopt point, objective, gradient, and tolerance data.
+function _prepare_manopt_solver_functions(
+    model_cost,
+    model_egrad,
+    M,
+    p0;
+    normA2 = nothing,
+    model_grad = nothing,
+    tol,
+    grad_tol = nothing,
+    normalized_objective::Bool,
+)
+    p0_local = _solver_point(M, p0)
+    T = _scalar_eltype(p0_local)
+    model_grad_raw = isnothing(model_grad) ? grad(model_egrad) : model_grad
+    model_grad_local = _layout_adapt_gradient(model_grad_raw)
+    objective_scale =
+        normalized_objective && !isnothing(normA2) && normA2 > 0 ? T(normA2) : one(T)
+    solver_cost, solver_grad, uses_relative_objective =
+        _relative_solver_functions(model_cost, model_grad_local, objective_scale)
+    return (
+        p0 = p0_local,
+        T = T,
+        solver_cost = solver_cost,
+        solver_grad = solver_grad,
+        uses_relative_objective = uses_relative_objective,
+        objective_scale = objective_scale,
+        grad_stop_tol = isnothing(grad_tol) ? T(tol) : T(grad_tol),
+        dual_grad_tol = _dual_stop_grad_tol(T, tol, grad_tol),
+    )
+end
+
+# Build the common Manopt stopping rule for iteration, gradient, and dual criteria.
+function _manopt_stopping(maxiter::Int, grad_stop_tol, dual_stop; extra = ())
+    return StopWhenAny(
+        StopAfterIteration(maxiter),
+        StopWhenGradientNormLess(grad_stop_tol),
+        extra...,
+        dual_stop,
+    )
+end
+
+# Create progress and debug callbacks shared by Manopt-backed solvers.
+function _manopt_callbacks(
+    make_progress::Function,
+    maxiter::Int,
+    verbose::Bool,
+    solver_cost,
+    solver_grad,
+    M;
+    diagnostics_recorder = nothing,
+    post_step_callback = nothing,
+    iteration_callbacks = (),
+)
+    progress = maxiter > 0 ? make_progress(maxiter) : NoMethodProgress()
+    diagnostics_callback =
+        isnothing(diagnostics_recorder) ? nothing :
+        _solver_diagnostics_callback(diagnostics_recorder)
+    progress_callback = _solver_progress_callback(
+        progress,
+        solver_cost,
+        solver_grad,
+        M;
+        diagnostics_recorder,
+    )
+    debug_actions = _solver_debug_actions(
+        verbose,
+        post_step_callback,
+        diagnostics_callback,
+        progress_callback,
+        iteration_callbacks...,
+    )
+    return (
+        progress = progress,
+        diagnostics_callback = diagnostics_callback,
+        progress_callback = progress_callback,
+        debug_actions = debug_actions,
+    )
+end
+
+# Finish progress, collect diagnostics, and return common solver stats.
+function _manopt_finish_result(
+    p_opt,
+    state,
+    progress,
+    diagnostics_recorder,
+    solver_cost,
+    solver_grad,
+    M,
+    normA2;
+    tol_T,
+    maxiter::Int,
+    solver::Symbol,
+    tiny_grad_tol,
+    return_stats::Bool,
+    verbose::Bool,
+    normalized_objective::Bool,
+    solver_info_extra = (;),
+)
+    iterations_done = _solver_iterations(state, maxiter)
+    if verbose
+        finish_progress!(
+            progress;
+            current = iterations_done,
+            showvalues = Any[("Status", "Finished"), ("Iterations", iterations_done)],
+        )
+    end
+    solver_info =
+        isnothing(diagnostics_recorder) ? (;) :
+        _solver_info(diagnostics_recorder, iterations_done)
+    solver_info = merge(solver_info, solver_info_extra)
+    if !return_stats
+        return p_opt
+    end
+    return _solver_stats(
+        solver_cost,
+        solver_grad,
+        M,
+        p_opt,
+        state,
+        normA2;
+        tol_T,
+        maxiter,
+        solver,
+        tiny_grad_tol,
+        solver_info,
+        normalized_objective,
+    )
+end
+
 # Query Manopt's convergence flag through one local compatibility point.
 @inline _solver_has_converged(state) = Manopt.has_converged(state)
 
diff --git a/src/solvers/rcg.jl b/src/solvers/rcg.jl
index afed20e..7204a10 100644
--- a/src/solvers/rcg.jl
+++ b/src/solvers/rcg.jl
@@ -108,90 +108,69 @@ function solve_rcg(
     grad_tol = nothing,
     normalized_objective::Bool = true,
 )
-    p0_local = _solver_point(M, p0)
-    T = _scalar_eltype(p0_local)
-    model_grad_raw = isnothing(model_grad) ? grad(model_egrad) : model_grad
-    model_grad_local = _layout_adapt_gradient(model_grad_raw)
-    objective_scale =
-        normalized_objective && !isnothing(normA2) && normA2 > 0 ? T(normA2) : one(T)
-    solver_cost, solver_grad, uses_relative_objective =
-        _relative_solver_functions(model_cost, model_grad_local, objective_scale)
+    setup = _prepare_manopt_solver_functions(
+        model_cost,
+        model_egrad,
+        M,
+        p0;
+        normA2,
+        model_grad,
+        tol,
+        grad_tol,
+        normalized_objective,
+    )
+    p0_local = setup.p0
+    T = setup.T
     retraction_method = _solver_retraction_method(M, p0_local)
     transport =
         isnothing(vector_transport_method) ?
         _default_vector_transport_method(M, p0_local, retraction_method) :
         vector_transport_method
-    grad_stop_tol = isnothing(grad_tol) ? T(tol) : T(grad_tol)
-    tol_g = _dual_stop_grad_tol(T, tol, grad_tol)
+    tol_g = setup.dual_grad_tol
     dual_stop = StopWhenCostRelChangeAndGradientLess(T(tol), tol_g)
 
-    stopping = StopWhenAny(
-        StopAfterIteration(maxiter),
-        StopWhenGradientNormLess(grad_stop_tol),
-        dual_stop,
-    )
-    progress =
-        maxiter > 0 ?
-        make_rcg_progress(maxiter; enabled = verbose, phase = :refinement, dt = 0.2) :
-        NoMethodProgress()
-    diagnostics_callback =
-        isnothing(diagnostics_recorder) ? nothing :
-        _solver_diagnostics_callback(diagnostics_recorder)
-    progress_callback = _solver_progress_callback(
-        progress,
-        solver_cost,
-        solver_grad,
+    stopping = _manopt_stopping(maxiter, setup.grad_stop_tol, dual_stop)
+    callbacks = _manopt_callbacks(
+        n -> make_rcg_progress(n; enabled = verbose, phase = :refinement, dt = 0.2),
+        maxiter,
+        verbose,
+        setup.solver_cost,
+        setup.solver_grad,
         M;
         diagnostics_recorder,
+        post_step_callback,
+        iteration_callbacks,
     )
 
     state = conjugate_gradient_descent(
         M,
-        solver_cost,
-        solver_grad,
+        setup.solver_cost,
+        setup.solver_grad,
         p0_local;
         retraction_method = retraction_method,
         vector_transport_method = transport,
         stopping_criterion = stopping,
-        debug = _solver_debug_actions(
-            verbose,
-            post_step_callback,
-            diagnostics_callback,
-            progress_callback,
-            iteration_callbacks...,
-        ),
+        debug = callbacks.debug_actions,
         count = [:Cost, :Gradient],
         return_state = true,
     )
 
-    p_opt = get_solver_result(state)
-    iterations_done = _solver_iterations(state, maxiter)
-    if verbose
-        finish_progress!(
-            progress;
-            current = iterations_done,
-            showvalues = Any[("Status", "Finished"), ("Iterations", iterations_done)],
-        )
-    end
-    solver_info =
-        isnothing(diagnostics_recorder) ? (;) :
-        _solver_info(diagnostics_recorder, iterations_done)
-    if !return_stats
-        return p_opt
-    end
-    return _solver_stats(
-        solver_cost,
-        solver_grad,
-        M,
-        p_opt,
+    return _manopt_finish_result(
+        get_solver_result(state),
         state,
+        callbacks.progress,
+        diagnostics_recorder,
+        setup.solver_cost,
+        setup.solver_grad,
+        M,
         normA2;
         tol_T = T(tol),
         maxiter,
         solver = :rcg,
         tiny_grad_tol = tol_g,
-        solver_info,
-        normalized_objective = uses_relative_objective,
+        return_stats,
+        verbose,
+        normalized_objective = setup.uses_relative_objective,
     )
 end
 
diff --git a/src/solvers/rgd.jl b/src/solvers/rgd.jl
index 2503f61..e7517b5 100644
--- a/src/solvers/rgd.jl
+++ b/src/solvers/rgd.jl
@@ -21,26 +21,29 @@ function solve_rgd(
     grad_tol = nothing,
     normalized_objective::Bool = true,
 )
-    p0_local = _solver_point(M, p0)
-    T = _scalar_eltype(p0_local)
-    model_grad_raw = isnothing(model_grad) ? grad(model_egrad) : model_grad
-    model_grad_local = _layout_adapt_gradient(model_grad_raw)
-    objective_scale =
-        normalized_objective && !isnothing(normA2) && normA2 > 0 ? T(normA2) : one(T)
-    solver_cost_base, solver_grad, uses_relative_objective =
-        _relative_solver_functions(model_cost, model_grad_local, objective_scale)
+    setup = _prepare_manopt_solver_functions(
+        model_cost,
+        model_egrad,
+        M,
+        p0;
+        normA2,
+        model_grad,
+        tol,
+        grad_tol,
+        normalized_objective,
+    )
+    p0_local = setup.p0
+    T = setup.T
     retraction_method = _solver_retraction_method(M, p0_local)
-    stepsize_eff_base = T(stepsize) * objective_scale
-    armijo_alpha_min = T(1e-8) * objective_scale
-    grad_stop_tol = isnothing(grad_tol) ? T(tol) : T(grad_tol)
-    tol_g = _dual_stop_grad_tol(T, tol, grad_tol)
+    stepsize_eff_base = T(stepsize) * setup.objective_scale
+    armijo_alpha_min = T(1e-8) * setup.objective_scale
+    tol_g = setup.dual_grad_tol
     dual_stop = StopWhenCostRelChangeAndGradientLess(T(tol), tol_g)
-
-    stopping = StopWhenAny(
-        StopAfterIteration(maxiter),
-        StopWhenGradientNormLess(grad_stop_tol),
-        StopWhenStepsizeLess(armijo_alpha_min),
-        dual_stop,
+    stopping = _manopt_stopping(
+        maxiter,
+        setup.grad_stop_tol,
+        dual_stop;
+        extra = (StopWhenStepsizeLess(armijo_alpha_min),),
     )
 
     use_squaring_armijo = _contains_sqeuclidean_manifold(M)
@@ -50,7 +53,7 @@ function solve_rgd(
         _adaptive_initial_stepsize(
             M,
             p0_local,
-            solver_grad,
+            setup.solver_grad,
             retraction_method,
             stepsize_eff_base;
             alpha_min = armijo_alpha_min,
@@ -64,7 +67,7 @@ function solve_rgd(
     armijo_additional_decrease =
         use_strict_sqeuclidean ? ((M, q) -> _all_finite(q)) : ((M, q) -> true)
     solver_cost =
-        use_strict_sqeuclidean ? _safe_cost_function(solver_cost_base) : solver_cost_base
+        use_strict_sqeuclidean ? _safe_cost_function(setup.solver_cost) : setup.solver_cost
     armijo = Manopt.ArmijoLinesearch(
         M;
         retraction_method = retraction_method,
@@ -78,70 +81,48 @@ function solve_rgd(
         additional_decrease_condition = armijo_additional_decrease,
     )
 
-    progress =
-        maxiter > 0 ?
-        make_rgd_progress(maxiter; enabled = verbose, phase = :refinement, dt = 0.2) :
-        NoMethodProgress()
-    diagnostics_callback =
-        isnothing(diagnostics_recorder) ? nothing :
-        _solver_diagnostics_callback(diagnostics_recorder)
-    progress_callback = _solver_progress_callback(
-        progress,
+    callbacks = _manopt_callbacks(
+        n -> make_rgd_progress(n; enabled = verbose, phase = :refinement, dt = 0.2),
+        maxiter,
+        verbose,
         solver_cost,
-        solver_grad,
+        setup.solver_grad,
         M;
         diagnostics_recorder,
+        post_step_callback,
+        iteration_callbacks,
     )
 
     state = gradient_descent(
         M,
         solver_cost,
-        solver_grad,
+        setup.solver_grad,
         p0_local;
         retraction_method = retraction_method,
         stepsize = armijo,
         stopping_criterion = stopping,
-        debug = _solver_debug_actions(
-            verbose,
-            post_step_callback,
-            diagnostics_callback,
-            progress_callback,
-            iteration_callbacks...,
-        ),
+        debug = callbacks.debug_actions,
         count = [:Cost, :Gradient],
         return_state = true,
     )
 
-    p_opt = get_solver_result(state)
-    iterations_done = _solver_iterations(state, maxiter)
-    if verbose
-        finish_progress!(
-            progress;
-            current = iterations_done,
-            showvalues = Any[("Status", "Finished"), ("Iterations", iterations_done)],
-        )
-    end
-    solver_info =
-        isnothing(diagnostics_recorder) ? (;) :
-        _solver_info(diagnostics_recorder, iterations_done)
-    solver_info =
-        merge(solver_info, (initial_stepsize_eff = Float64(initial_stepsize_eff),))
-    if !return_stats
-        return p_opt
-    end
-    return _solver_stats(
+    return _manopt_finish_result(
+        get_solver_result(state),
+        state,
+        callbacks.progress,
+        diagnostics_recorder,
         solver_cost,
-        solver_grad,
+        setup.solver_grad,
         M,
-        p_opt,
-        state,
         normA2;
         tol_T = T(tol),
         maxiter,
         solver = :rgd,
         tiny_grad_tol = tol_g,
-        solver_info,
-        normalized_objective = uses_relative_objective,
+        return_stats,
+        verbose,
+        normalized_objective = setup.uses_relative_objective,
+        solver_info_extra = (initial_stepsize_eff = Float64(initial_stepsize_eff),),
     )
 end
 
@@ -163,80 +144,65 @@ function solve_rgd_fixed(
     grad_tol = nothing,
     normalized_objective::Bool = true,
 )
-    p0_local = _solver_point(M, p0)
-    T = _scalar_eltype(p0_local)
-    model_grad_raw = isnothing(model_grad) ? grad(model_egrad) : model_grad
-    model_grad_local = _layout_adapt_gradient(model_grad_raw)
-    objective_scale =
-        normalized_objective && !isnothing(normA2) && normA2 > 0 ? T(normA2) : one(T)
-    solver_cost, solver_grad, uses_relative_objective =
-        _relative_solver_functions(model_cost, model_grad_local, objective_scale)
+    setup = _prepare_manopt_solver_functions(
+        model_cost,
+        model_egrad,
+        M,
+        p0;
+        normA2,
+        model_grad,
+        tol,
+        grad_tol,
+        normalized_objective,
+    )
+    p0_local = setup.p0
+    T = setup.T
     retraction_method = _solver_retraction_method(M, p0_local)
-    grad_stop_tol = isnothing(grad_tol) ? T(tol) : T(grad_tol)
     tiny_grad_tol = isnothing(grad_tol) ? T(1e-5) : T(grad_tol)
-    stopping =
-        StopWhenAny(StopAfterIteration(maxiter), StopWhenGradientNormLess(grad_stop_tol))
-    progress =
-        maxiter > 0 ?
-        make_rgd_fixed_progress(maxiter; enabled = verbose, phase = :refinement, dt = 0.2) :
-        NoMethodProgress()
-    diagnostics_callback =
-        isnothing(diagnostics_recorder) ? nothing :
-        _solver_diagnostics_callback(diagnostics_recorder)
-    progress_callback = _solver_progress_callback(
-        progress,
-        solver_cost,
-        solver_grad,
+    stopping = StopWhenAny(
+        StopAfterIteration(maxiter),
+        StopWhenGradientNormLess(setup.grad_stop_tol),
+    )
+    callbacks = _manopt_callbacks(
+        n -> make_rgd_fixed_progress(n; enabled = verbose, phase = :refinement, dt = 0.2),
+        maxiter,
+        verbose,
+        setup.solver_cost,
+        setup.solver_grad,
         M;
         diagnostics_recorder,
+        post_step_callback,
+        iteration_callbacks,
     )
     state = gradient_descent(
         M,
-        solver_cost,
-        solver_grad,
+        setup.solver_cost,
+        setup.solver_grad,
         p0_local;
         retraction_method = retraction_method,
-        stepsize = Manopt.ConstantStepsize(M, T(stepsize) * objective_scale),
+        stepsize = Manopt.ConstantStepsize(M, T(stepsize) * setup.objective_scale),
         stopping_criterion = stopping,
-        debug = _solver_debug_actions(
-            verbose,
-            post_step_callback,
-            diagnostics_callback,
-            progress_callback,
-            iteration_callbacks...,
-        ),
+        debug = callbacks.debug_actions,
         count = [:Cost, :Gradient],
         return_state = true,
     )
 
-    p_opt = get_solver_result(state)
-    iterations_done = _solver_iterations(state, maxiter)
-    if verbose
-        finish_progress!(
-            progress;
-            current = iterations_done,
-            showvalues = Any[("Status", "Finished"), ("Iterations", iterations_done)],
-        )
-    end
-    solver_info =
-        isnothing(diagnostics_recorder) ? (;) :
-        _solver_info(diagnostics_recorder, iterations_done)
-    if !return_stats
-        return p_opt
-    end
-    return _solver_stats(
-        solver_cost,
-        solver_grad,
-        M,
-        p_opt,
+    return _manopt_finish_result(
+        get_solver_result(state),
         state,
+        callbacks.progress,
+        diagnostics_recorder,
+        setup.solver_cost,
+        setup.solver_grad,
+        M,
         normA2;
         tol_T = T(tol),
         maxiter,
         solver = :rgd_fixed,
         tiny_grad_tol = tiny_grad_tol,
-        solver_info,
-        normalized_objective = uses_relative_objective,
+        return_stats,
+        verbose,
+        normalized_objective = setup.uses_relative_objective,
     )
 end
 

From a9656a6a3923ae03a008a50d0eb4532886e31f27 Mon Sep 17 00:00:00 2001
From: PBrdng <pbrdng@gmail.com>
Date: Wed, 17 Jun 2026 22:00:57 +0200
Subject: [PATCH 16/18] unify more code

---
 src/api/btd.jl            |  4 +-
 src/api/cpd.jl            |  2 +-
 src/core/types.jl         |  4 +-
 src/join/cpd_backend.jl   |  2 +-
 src/results/conversion.jl | 32 +++++---------
 src/results/rel_error.jl  | 14 +-----
 src/solvers/abstract.jl   | 90 ++++++++++++++++++++++++++-------------
 7 files changed, 77 insertions(+), 71 deletions(-)

diff --git a/src/api/btd.jl b/src/api/btd.jl
index aa06ded..97d6f02 100644
--- a/src/api/btd.jl
+++ b/src/api/btd.jl
@@ -24,7 +24,7 @@ function _polish_btd_with_als(
         max_stagnation_restarts = 0,
     )
     rel_error(als_res) < rel_error(result) || return result
-    si0 = hasproperty(result, :solver_info) ? solver_info(result) : (;)
+    si0 = _result_solver_info(result)
     si = merge(
         si0,
         (
@@ -127,7 +127,7 @@ function _btd_warm_start_result(
 end
 
 function _merge_btd_solver_info(result, extra::NamedTuple)
-    si0 = hasproperty(result, :solver_info) ? result.solver_info : (;)
+    si0 = _result_solver_info(result)
     return (
         point = result.point,
         cost = result.cost,
diff --git a/src/api/cpd.jl b/src/api/cpd.jl
index c06fcf3..e03ffab 100644
--- a/src/api/cpd.jl
+++ b/src/api/cpd.jl
@@ -9,7 +9,7 @@ function _pullback_eps_value(::Type{T}, pullback_eps) where {T<:AbstractFloat}
 end
 
 function _merge_res_solver_info(res, patch::NamedTuple)
-    si0 = hasproperty(res, :solver_info) ? solver_info(res) : (;)
+    si0 = _result_solver_info(res)
     return (
         point = point(res),
         cost = cost(res),
diff --git a/src/core/types.jl b/src/core/types.jl
index 8c781a3..d39ee6a 100644
--- a/src/core/types.jl
+++ b/src/core/types.jl
@@ -363,9 +363,7 @@ solver_info(r::Union{CPDResult,ApproxResult,BTDResult}) = r.solver_info
 
 Return the decoded components stored in a decomposition result.
 """
-components(r::CPDResult) = r.components
-components(r::ApproxResult) = r.components
-components(r::BTDResult) = r.components
+components(r::Union{CPDResult,ApproxResult,BTDResult}) = r.components
 components(r::NamedTuple) = getproperty(r, :components)
 
 """
diff --git a/src/join/cpd_backend.jl b/src/join/cpd_backend.jl
index 03493f9..123a6b5 100644
--- a/src/join/cpd_backend.jl
+++ b/src/join/cpd_backend.jl
@@ -109,7 +109,7 @@ end
 function _cpd_result(model::JoinModel{<:AbstractFloat,<:CPDBackend}, result, dims, r)
     m = cpd_model(model)
     solver_sym = _result_solver_symbol(solver(result))
-    si = hasproperty(result, :solver_info) ? solver_info(result) : (;)
+    si = _result_solver_info(result)
     als_family = solver_sym in _CP_ALS_FAMILY_SOLVERS
 
     if r == 1
diff --git a/src/results/conversion.jl b/src/results/conversion.jl
index 7c5795c..15755f8 100644
--- a/src/results/conversion.jl
+++ b/src/results/conversion.jl
@@ -2,12 +2,13 @@
 
 _result_solver_symbol(solver::Symbol) = solver
 _result_solver_symbol(solver) = :unknown
+_result_solver_info(result) = hasproperty(result, :solver_info) ? solver_info(result) : (;)
 
-function _to_approx_result(model::JoinModel{T}, result) where {T<:AbstractFloat}
+function _to_join_result(result_type, model::JoinModel{T}, result) where {T<:AbstractFloat}
     comps = extract_components(model, result.point)
     solver_sym = _result_solver_symbol(result.solver)
-    solver_info = hasproperty(result, :solver_info) ? result.solver_info : (;)
-    return ApproxResult(
+    si = _result_solver_info(result)
+    return result_type(
         result.point,
         comps,
         result.cost,
@@ -16,27 +17,16 @@ function _to_approx_result(model::JoinModel{T}, result) where {T<:AbstractFloat}
         result.iterations,
         result.converged,
         solver_sym,
-        solver_info,
+        si,
     )
 end
 
-function _to_btd_result(model::JoinModel{T}, result) where {T<:AbstractFloat}
-    comps = extract_components(model, result.point)
-    solver_sym = _result_solver_symbol(result.solver)
-    solver_info = hasproperty(result, :solver_info) ? result.solver_info : (;)
-    # Reuse solver-reported cost/error instead of reconstructing the full BTD residual again.
-    return BTDResult(
-        result.point,
-        comps,
-        result.cost,
-        result.rel_error,
-        result.grad_norm,
-        result.iterations,
-        result.converged,
-        solver_sym,
-        solver_info,
-    )
-end
+_to_approx_result(model::JoinModel{T}, result) where {T<:AbstractFloat} =
+    _to_join_result(ApproxResult, model, result)
+
+# Reuse solver-reported cost/error instead of reconstructing the full BTD residual again.
+_to_btd_result(model::JoinModel{T}, result) where {T<:AbstractFloat} =
+    _to_join_result(BTDResult, model, result)
 
 _to_cpd_result(model, result, dims, r) =
     throw(ArgumentError("No CPD result converter for model $(typeof(model))."))
diff --git a/src/results/rel_error.jl b/src/results/rel_error.jl
index c6ea060..567be1c 100644
--- a/src/results/rel_error.jl
+++ b/src/results/rel_error.jl
@@ -18,19 +18,7 @@ function rel_error(A::AbstractArray, Ahat::AbstractArray)
     return relative_frobenius_error(A, Ahat)
 end
 
-function rel_error(A::AbstractArray, res::CPDResult)
-    return relative_frobenius_error(A, reconstruct(res))
-end
-
-function rel_error(A::AbstractArray, tucker_res::TuckerResult)
-    return relative_frobenius_error(A, reconstruct(tucker_res))
-end
-
-function rel_error(A::AbstractArray, res::ApproxResult)
-    return relative_frobenius_error(A, reconstruct(res))
-end
-
-function rel_error(A::AbstractArray, res::BTDResult)
+function rel_error(A::AbstractArray, res::Union{CPDResult,TuckerResult,ApproxResult,BTDResult})
     return relative_frobenius_error(A, reconstruct(res))
 end
 
diff --git a/src/solvers/abstract.jl b/src/solvers/abstract.jl
index 52afb35..5aa7608 100644
--- a/src/solvers/abstract.jl
+++ b/src/solvers/abstract.jl
@@ -212,6 +212,50 @@ function solve(solver::AbstractSolver, model::AbstractDecompositionModel; kwargs
     error("solve not implemented for $(typeof(solver))")
 end
 
+function _solve_ro_solver(
+    solver::AbstractROSolver,
+    model::AbstractDecompositionModel;
+    init,
+    p0,
+    maxiter::Int,
+    tol::Real,
+    gradient_mode,
+    normalization::Union{AbstractNormalizationPolicy,Symbol,Nothing},
+    verbose::Bool,
+    return_stats::Bool,
+    vector_transport_method::Union{ManifoldsBase.AbstractVectorTransportMethod,Nothing},
+    grad_tol,
+    normalized_objective::Bool,
+    iteration_callbacks,
+    diagnostics_recorder,
+    run_solver::Function,
+)
+    setup = _prepare_solver_problem(model; init, p0, gradient_mode, verbose)
+    normalization_policy = _normalization_policy(normalization)
+    supports_normalization_policy(model, normalization_policy) || throw(
+        ArgumentError(
+            "Normalization policy $(typeof(normalization_policy)) is not supported for model $(typeof(model)).",
+        ),
+    )
+    solver_sym = solver_symbol(solver)
+    post_step_callback =
+        _solver_post_step_callback(model, setup.M, normalization_policy, solver_sym)
+    return run_solver(
+        solver,
+        setup;
+        maxiter,
+        tol,
+        verbose,
+        return_stats,
+        vector_transport_method,
+        grad_tol,
+        normalized_objective,
+        post_step_callback,
+        diagnostics_recorder,
+        iteration_callbacks,
+    )
+end
+
 function solve(
     solver::AbstractFirstOrderROSolver,
     model::AbstractDecompositionModel{T};
@@ -228,30 +272,23 @@ function solve(
     normalized_objective::Bool = true,
     iteration_callbacks = (),
 ) where {T<:AbstractFloat}
-    setup = _prepare_solver_problem(model; init, p0, gradient_mode, verbose)
-    normalization_policy = _normalization_policy(normalization)
-    supports_normalization_policy(model, normalization_policy) || throw(
-        ArgumentError(
-            "Normalization policy $(typeof(normalization_policy)) is not supported for model $(typeof(model)).",
-        ),
-    )
-    solver_sym = solver_symbol(solver)
-    post_step_callback =
-        _solver_post_step_callback(model, setup.M, normalization_policy, solver_sym)
-    diagnostics_recorder = first_order_diagnostics_recorder(solver)
-    return run_first_order_solver(
+    return _solve_ro_solver(
         solver,
-        setup;
+        model;
+        init,
+        p0,
         maxiter,
         tol,
+        gradient_mode,
+        normalization,
         verbose,
         return_stats,
         vector_transport_method,
         grad_tol,
         normalized_objective,
-        post_step_callback,
-        diagnostics_recorder,
         iteration_callbacks,
+        diagnostics_recorder = first_order_diagnostics_recorder(solver),
+        run_solver = run_first_order_solver,
     )
 end
 
@@ -271,30 +308,23 @@ function solve(
     normalized_objective::Bool = true,
     iteration_callbacks = (),
 ) where {T<:AbstractFloat}
-    setup = _prepare_solver_problem(model; init, p0, gradient_mode)
-    normalization_policy = _normalization_policy(normalization)
-    supports_normalization_policy(model, normalization_policy) || throw(
-        ArgumentError(
-            "Normalization policy $(typeof(normalization_policy)) is not supported for model $(typeof(model)).",
-        ),
-    )
-    solver_sym = solver_symbol(solver)
-    post_step_callback =
-        _solver_post_step_callback(model, setup.M, normalization_policy, solver_sym)
-    diagnostics_recorder = second_order_diagnostics_recorder(solver)
-    return run_second_order_solver(
+    return _solve_ro_solver(
         solver,
-        setup;
+        model;
+        init,
+        p0,
         maxiter,
         tol,
+        gradient_mode,
+        normalization,
         verbose,
         return_stats,
         vector_transport_method,
         grad_tol,
         normalized_objective,
-        post_step_callback,
-        diagnostics_recorder,
         iteration_callbacks,
+        diagnostics_recorder = second_order_diagnostics_recorder(solver),
+        run_solver = run_second_order_solver,
     )
 end
 

From 85a86d9c046580f6379025aed05c0fa969cdf664 Mon Sep 17 00:00:00 2001
From: Se Eun Choi <seeun.choi@uni-osnabrueck.de>
Date: Thu, 18 Jun 2026 00:26:25 +0000
Subject: [PATCH 17/18] run JuliaFormatter

---
 src/results/rel_error.jl | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/results/rel_error.jl b/src/results/rel_error.jl
index 567be1c..b46ed02 100644
--- a/src/results/rel_error.jl
+++ b/src/results/rel_error.jl
@@ -18,7 +18,10 @@ function rel_error(A::AbstractArray, Ahat::AbstractArray)
     return relative_frobenius_error(A, Ahat)
 end
 
-function rel_error(A::AbstractArray, res::Union{CPDResult,TuckerResult,ApproxResult,BTDResult})
+function rel_error(
+    A::AbstractArray,
+    res::Union{CPDResult,TuckerResult,ApproxResult,BTDResult},
+)
     return relative_frobenius_error(A, reconstruct(res))
 end
 

From 718dcbf9fded5b5981f7a9afa168f560e93c3805 Mon Sep 17 00:00:00 2001
From: Se Eun Choi <seeun.choi@uni-osnabrueck.de>
Date: Thu, 18 Jun 2026 00:32:22 +0000
Subject: [PATCH 18/18] Format: apply JuliaFormatter

---
 src/core/tensor_ops.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/core/tensor_ops.jl b/src/core/tensor_ops.jl
index e5c1427..5ab2f89 100644
--- a/src/core/tensor_ops.jl
+++ b/src/core/tensor_ops.jl
@@ -470,8 +470,8 @@ gradU_column_cp(
 
 cp_reconstruction_norm2(components::Vector{RankOneTensor{T}}) where {T<:AbstractFloat} =
     sum(
-        cross_component(components[i], components[j]) for
-        i in eachindex(components), j in eachindex(components)
+        cross_component(components[i], components[j]) for i in eachindex(components),
+        j in eachindex(components)
     )
 
 function cp_inner_AX(