diff --git a/Project.toml b/Project.toml
index 5928bad..ad24f56 100644
--- a/Project.toml
+++ b/Project.toml
@@ -10,15 +10,15 @@ Manopt = "0fc0a36d-df90-57f3-8f93-d78a9fc72bb5"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-LinearAlgebra = "1.12.0"
 Manifolds = "0.11.20"
 ManifoldsBase = "2.3.5"
-Manopt = "0.5.36"
+Manopt = "0.5.37"
 ProgressMeter = "1.11.0"
-Random = "1.11.0"
 RecursiveArrayTools = "4.3"
 TensorOperations = "5.6"
-julia = "1.10"
\ No newline at end of file
+julia = "1.10"
diff --git a/src/api/btd.jl b/src/api/btd.jl
index 476cf91..94c1b02 100644
--- a/src/api/btd.jl
+++ b/src/api/btd.jl
@@ -262,44 +262,51 @@ function btd(
     solve_init = init_eff
     solve_p0 = nothing
     warm_info = (;)
-    if solver ∉ (:als,) && init_eff isa BTDALSWarmStartInit
-        warm = _btd_warm_start_result(
-            model,
-            b,
-            init_eff,
-            solver;
+    short_circuited = Ref(false)
+    result = with_phase_progress() do
+        if solver ∉ (:als,) && init_eff isa BTDALSWarmStartInit
+            warm = _btd_warm_start_result(
+                model,
+                b,
+                init_eff,
+                solver;
+                verbose = verbose,
+                warm_rel_error_gate,
+            )
+            solve_init = warm.init
+            solve_p0 = warm.p0
+            warm_info = warm.warm_info
+            if !isnothing(warm.short_circuit)
+                short_circuited[] = true
+                return warm.short_circuit
+            end
+        end
+        _solve_model(
+            model;
+            init = solve_init,
+            p0 = solve_p0,
+            solver = solver,
+            maxiter = maxiter,
+            stepsize = stepsize,
+            tol = tol,
+            gradient_mode = gradient_mode,
+            normalization = NoNormalization(),
             verbose = verbose,
-            warm_rel_error_gate,
+            vector_transport_method = vector_transport_method,
+            block_method = block_method,
+            block_maxiter = block_maxiter,
+            max_stagnation_restarts = max_stagnation_restarts,
+            stagnation_rel_error = stagnation_rel_error,
+            restart_candidates = restart_candidates,
+            restart_screening_steps = restart_screening_steps,
+            restart_block_maxiter = restart_block_maxiter,
+            restart_seed = restart_seed,
+            kwargs...,
         )
-        solve_init = warm.init
-        solve_p0 = warm.p0
-        warm_info = warm.warm_info
-        if !isnothing(warm.short_circuit)
-            return _to_btd_result(model, warm.short_circuit)
-        end
     end
-    result = _solve_model(
-        model;
-        init = solve_init,
-        p0 = solve_p0,
-        solver = solver,
-        maxiter = maxiter,
-        stepsize = stepsize,
-        tol = tol,
-        gradient_mode = gradient_mode,
-        normalization = NoNormalization(),
-        verbose = verbose,
-        vector_transport_method = vector_transport_method,
-        block_method = block_method,
-        block_maxiter = block_maxiter,
-        max_stagnation_restarts = max_stagnation_restarts,
-        stagnation_rel_error = stagnation_rel_error,
-        restart_candidates = restart_candidates,
-        restart_screening_steps = restart_screening_steps,
-        restart_block_maxiter = restart_block_maxiter,
-        restart_seed = restart_seed,
-        kwargs...,
-    )
+    if short_circuited[]
+        return _to_btd_result(model, result)
+    end
     result =
         isempty(propertynames(warm_info)) ? result :
         _merge_btd_solver_info(result, warm_info)
diff --git a/src/api/cpd.jl b/src/api/cpd.jl
index 59ea874..119cd09 100644
--- a/src/api/cpd.jl
+++ b/src/api/cpd.jl
@@ -87,6 +87,7 @@ function _cpd_als_warm_then_pack(
             nonnegative = inner.nonnegative,
             verbose = verbose,
             return_stats = true,
+            progress_phase = :initialization,
         )
         return pack_cpd_point(target, CPDPoint(warm_out.weights, warm_out.factors))
     end
@@ -112,6 +113,7 @@ function _cpd_als_warm_then_pack(
         verbose = verbose,
         vector_transport_method = nothing,
         nonnegative = inner.nonnegative,
+        progress_phase = :initialization,
         kwargs...,
     )
     warm_cpd = _to_cpd_result(warm_model, warm_result, size(A), r)
@@ -269,20 +271,21 @@ function _cpd_impl(
         use_pullback_metric = (geometry_eff == :squaring_metric),
         pullback_eps = pullback_eps_eff,
     )
-    p_solve = if init_eff isa ALSWarmStartInit && isnothing(p0) && solver ∈ manifold_solvers
-        _cpd_als_warm_then_pack(
-            model,
-            init_eff;
-            tol,
-            normalization = als_normalization_eff,
-            verbose,
-            pullback_eps = pullback_eps_eff,
-            kwargs...,
-        )
-    else
-        _pack_cpd_explicit_p0(model, p0)
-    end
     result_rgd = with_phase_progress() do
+        p_solve =
+            if init_eff isa ALSWarmStartInit && isnothing(p0) && solver ∈ manifold_solvers
+                _cpd_als_warm_then_pack(
+                    model,
+                    init_eff;
+                    tol,
+                    normalization = als_normalization_eff,
+                    verbose,
+                    pullback_eps = pullback_eps_eff,
+                    kwargs...,
+                )
+            else
+                _pack_cpd_explicit_p0(model, p0)
+            end
         _solve_model(
             model;
             init = init_eff,
diff --git a/src/solvers/cp_als.jl b/src/solvers/cp_als.jl
index dc797ad..b908dea 100644
--- a/src/solvers/cp_als.jl
+++ b/src/solvers/cp_als.jl
@@ -478,6 +478,7 @@ function solve(
     nonnegative::Bool = false,
     verbose::Bool = true,
     return_stats::Bool = false,
+    progress_phase::Symbol = :refinement,
     kwargs...,
 ) where {T<:AbstractFloat}
     A, r = cp_als_data(model)
@@ -497,7 +498,7 @@ function solve(
         normalization = normalization,
         verbose = verbose,
         return_stats = true,
-        progress_phase = :refinement,
+        progress_phase = progress_phase,
         mttkrp_method = get(kwargs, :mttkrp_method, :auto),
         nonnegative,
     )
diff --git a/src/solvers/rgd.jl b/src/solvers/rgd.jl
index 0c627e6..e458f63 100644
--- a/src/solvers/rgd.jl
+++ b/src/solvers/rgd.jl
@@ -2,6 +2,13 @@
 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
@@ -339,36 +346,36 @@ function _solver_stats(
     )
 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...) =
-    Any[cb for cb in callbacks if !isnothing(cb)]
+_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 "; at_init = true),
-            Manopt.DebugCost(; format = "f(x): %.6e", at_init = true),
-            Manopt.DebugGradientNorm(; format = "|grad f(p)|:%.6e", at_init = true),
-            Manopt.DebugDivider("\n"; at_init = true),
+            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(; format = "# %-6d"),
-                Manopt.DebugDivider(" "; at_init = true),
-                Manopt.DebugCost(; format = "f(x): %.6e", at_init = true),
-                Manopt.DebugGradientNorm(; format = "|grad f(p)|:%.6e", at_init = true),
-                Manopt.DebugDivider("\n"; at_init = true),
+                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,
         )
-        actions = Any[init_group, iter_group]
-    else
-        actions = Any[]
-    end
-    for callback in callbacks
-        !isnothing(callback) && push!(actions, callback)
+        iteration_actions = Any[iter_group]
+        append!(iteration_actions, callback_actions)
+        return Any[:Start=>Any[init_group], :Iteration=>iteration_actions]
     end
-    return actions
+    return callback_actions
 end
 
 function _solver_progress_callback(