Add sigma_cauchy and scp_norm to solver statistics in LMTR, LMSolver, R2DH, R2N, and R2Solver

MohamedLaghdafHABIBOULLAH · MohamedLaghdafHABIBOULLAH · commit 6026bf3965a9 · 2025-10-15T14:35:56.000-04:00
diff --git a/.gitignore b/.gitignore
@@ -0,0 +1 @@
+Manifest.toml
diff --git a/src/LMTR_alg.jl b/src/LMTR_alg.jl
@@ -285,6 +285,7 @@ function SolverCore.solve!(
   set_objective!(stats, fk + hk)
   set_solver_specific!(stats, :smooth_obj, fk)
   set_solver_specific!(stats, :nonsmooth_obj, hk)
+  set_solver_specific!(stats, :sigma_cauchy, 1/ν)
   set_solver_specific!(stats, :prox_evals, prox_evals + 1)
 
   φ1 = let Fk = Fk, ∇fk = ∇fk
@@ -302,6 +303,7 @@ function SolverCore.solve!(
   # Take first proximal gradient step s1 and see if current xk is nearly stationary.
   # s1 minimizes φ1(d) + ‖d‖² / 2 / ν + ψ(d) ⟺ s1 ∈ prox{νψ}(-ν∇φ1(0))
   prox!(s, ψ, mν∇fk, ν)
+  set_solver_specific!(stats, :scp_norm, norm(s))
   ξ1 = fk + hk - mk1(s) + max(1, abs(fk + hk)) * 10 * eps()
   sqrt_ξ1_νInv = ξ1 ≥ 0 ? sqrt(ξ1 / ν) : sqrt(-ξ1 / ν)
   solved = (ξ1 < 0 && sqrt_ξ1_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ1_νInv ≤ atol)
@@ -448,9 +450,11 @@ function SolverCore.solve!(
     set_solver_specific!(stats, :prox_evals, prox_evals + 1)
 
     ν = α * Δk / (1 + σmax^2 * (α * Δk + 1))
+    set_solver_specific!(stats, :sigma_cauchy, 1/ν)
     @. mν∇fk = -∇fk * ν
 
     prox!(s, ψ, mν∇fk, ν)
+    set_solver_specific!(stats, :scp_norm, norm(s))
     mks = mk1(s)
 
     ξ1 = fk + hk - mks + max(1, abs(hk)) * 10 * eps()
diff --git a/src/LM_alg.jl b/src/LM_alg.jl
@@ -290,6 +290,7 @@ function SolverCore.solve!(
   set_objective!(stats, fk + hk)
   set_solver_specific!(stats, :smooth_obj, fk)
   set_solver_specific!(stats, :nonsmooth_obj, hk)
+  set_solver_specific!(stats, :sigma_cauchy, 1/ν)
   set_solver_specific!(stats, :prox_evals, prox_evals + 1)
   m_monotone > 1 && (m_fh_hist[stats.iter % (m_monotone - 1) + 1] = fk + hk)
 
@@ -306,6 +307,7 @@ function SolverCore.solve!(
   end
 
   prox!(s, ψ, mν∇fk, ν)
+  set_solver_specific!(stats, :scp_norm, norm(s))
   ξ1 = fk + hk - mk1(s) + max(1, abs(fk + hk)) * 10 * eps()
   sqrt_ξ1_νInv = ξ1 ≥ 0 ? sqrt(ξ1 / ν) : sqrt(-ξ1 / ν)
   solved = (ξ1 < 0 && sqrt_ξ1_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ1_νInv ≤ atol)
@@ -436,9 +438,11 @@ function SolverCore.solve!(
     set_solver_specific!(stats, :prox_evals, prox_evals + 1)
 
     ν = θ / (σmax^2 + σk) # ‖J'J + σₖ I‖ = ‖J‖² + σₖ
+    set_solver_specific!(stats, :sigma_cauchy, 1/ν)
 
     @. mν∇fk = - ν * ∇fk
     prox!(s, ψ, mν∇fk, ν)
+    set_solver_specific!(stats, :scp_norm, norm(s))
     mks = mk1(s)
 
     ξ1 = fk + hk - mks + max(1, abs(hk)) * 10 * eps()
diff --git a/src/R2DH.jl b/src/R2DH.jl
@@ -328,6 +328,7 @@ function SolverCore.solve!(
   mk(d)::T = φ(d) + ψ(d)::T
 
   spectral_test ? prox!(s, ψ, mν∇fk, ν₁) : iprox!(s, ψ, ∇fk, dkσk)
+  set_solver_specific!(stats, :scp_norm, norm(s))
 
   mks = mk(s)
 
@@ -412,19 +413,20 @@ function SolverCore.solve!(
     set_solver_specific!(stats, :smooth_obj, fk)
     set_solver_specific!(stats, :nonsmooth_obj, hk)
     set_solver_specific!(stats, :sigma, σk)
-    set_solver_specific!(stats, :sigma_cauchy, 1/ν₁)
     set_iter!(stats, stats.iter + 1)
     set_time!(stats, time() - start_time)
 
     @. dkσk = D.d .+ σk
     DNorm = norm(D.d, Inf)
 
     ν₁ = θ / (DNorm + σk)
+    set_solver_specific!(stats, :sigma_cauchy, 1/ν₁)
 
     @. mν∇fk = -ν₁ * ∇fk
     m_monotone > 1 && (m_fh_hist[stats.iter % (m_monotone - 1) + 1] = fk + hk)
 
     spectral_test ? prox!(s, ψ, mν∇fk, ν₁) : iprox!(s, ψ, ∇fk, dkσk)
+    set_solver_specific!(stats, :scp_norm, norm(s))
     mks = mk(s)
 
     ξ = hk - mks + max(1, abs(hk)) * 10 * eps()
diff --git a/src/R2N.jl b/src/R2N.jl
@@ -331,6 +331,7 @@ function SolverCore.solve!(
   end
 
   prox!(s1, ψ, mν∇fk, ν₁)
+  set_solver_specific!(stats, :scp_norm, norm(s1))
   mks = mk1(s1)
 
   ξ1 = hk - mks + max(1, abs(hk)) * 10 * eps()
@@ -477,6 +478,7 @@ function SolverCore.solve!(
 
     @. mν∇fk = - ν₁ * ∇fk
     prox!(s1, ψ, mν∇fk, ν₁)
+    set_solver_specific!(stats, :scp_norm, norm(s1))
     mks = mk1(s1)
 
     ξ1 = hk - mks + max(1, abs(hk)) * 10 * eps()
diff --git a/src/R2_alg.jl b/src/R2_alg.jl
@@ -402,6 +402,7 @@ function SolverCore.solve!(
   mk(d)::T = φk(d) + ψ(d)::T
 
   prox!(s, ψ, mν∇fk, ν)
+  set_solver_specific!(stats, :scp_norm, norm(s))
   mks = mk(s)
 
   ξ = hk - mks + max(1, abs(hk)) * 10 * eps()
@@ -489,6 +490,7 @@ function SolverCore.solve!(
     set_time!(stats, time() - start_time)
 
     prox!(s, ψ, mν∇fk, ν)
+    set_solver_specific!(stats, :scp_norm, norm(s))
     mks = mk(s)
 
     ξ = hk - mks + max(1, abs(hk)) * 10 * eps()
diff --git a/src/TRDH_alg.jl b/src/TRDH_alg.jl
@@ -333,6 +333,7 @@ function SolverCore.solve!(
   set_objective!(stats, fk + hk)
   set_solver_specific!(stats, :smooth_obj, fk)
   set_solver_specific!(stats, :nonsmooth_obj, hk)
+  set_solver_specific!(stats, :sigma_cauchy, 1/ν)
 
   # models
   φ1 = let ∇fk = ∇fk
@@ -358,6 +359,7 @@ function SolverCore.solve!(
 
   if reduce_TR
     prox!(s, ψ, mν∇fk, ν)
+    set_solver_specific!(stats, :scp_norm, norm(s))
     mks = mk1(s)
 
     ξ1 = hk - mks + max(1, abs(hk)) * 10 * eps()
@@ -383,6 +385,7 @@ function SolverCore.solve!(
   sNorm = χ(s)
 
   if !reduce_TR
+    set_solver_specific!(stats, :scp_norm, norm(s))
     sqrt_ξ_νInv = ξ ≥ 0 ? sqrt(ξ / ν) : sqrt(-ξ / ν)
     solved = (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv < atol)
     (ξ < 0 && sqrt_ξ_νInv > neg_tol) &&
@@ -473,10 +476,12 @@ function SolverCore.solve!(
     set_time!(stats, time() - start_time)
 
     ν = reduce_TR ? (α * Δk)/(DNorm + one(T)) : α / (DNorm + one(T))
+    set_solver_specific!(stats, :sigma_cauchy, 1/ν)
     mν∇fk .= -ν .* ∇fk
 
     if reduce_TR
       prox!(s, ψ, mν∇fk, ν)
+      set_solver_specific!(stats, :scp_norm, norm(s))
       ξ1 = hk - mk1(s) + max(1, abs(hk)) * 10 * eps()
       sqrt_ξ_νInv = ξ1 ≥ 0 ? sqrt(ξ1 / ν) : sqrt(-ξ1 / ν)
       solved = (ξ1 < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ_νInv < atol)
@@ -489,6 +494,7 @@ function SolverCore.solve!(
     sNorm = χ(s)
 
     if !reduce_TR
+      set_solver_specific!(stats, :scp_norm, norm(s))
       ξ = hk - mk(s) + max(1, abs(hk)) * 10 * eps()
       sqrt_ξ_νInv = ξ ≥ 0 ? sqrt(ξ / ν) : sqrt(-ξ / ν)
       solved = (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv < atol)
diff --git a/src/TR_alg.jl b/src/TR_alg.jl
@@ -314,6 +314,7 @@ function SolverCore.solve!(
   set_objective!(stats, fk + hk)
   set_solver_specific!(stats, :smooth_obj, fk)
   set_solver_specific!(stats, :nonsmooth_obj, hk)
+  set_solver_specific!(stats, :sigma_cauchy, 1/ν₁)
   set_solver_specific!(stats, :prox_evals, prox_evals + 1)
   m_monotone > 1 && (m_fh_hist[stats.iter % (m_monotone - 1) + 1] = fk + hk)
 
@@ -330,6 +331,7 @@ function SolverCore.solve!(
   end
 
   prox!(s, ψ, mν∇fk, ν₁)
+  set_solver_specific!(stats, :scp_norm, norm(s))
   ξ1 = hk - mk1(s) + max(1, abs(hk)) * 10 * eps()
   ξ1 > 0 || error("TR: first prox-gradient step should produce a decrease but ξ1 = $(ξ1)")
   sqrt_ξ1_νInv = sqrt(ξ1 / ν₁)
@@ -489,9 +491,11 @@ function SolverCore.solve!(
     set_solver_specific!(stats, :prox_evals, prox_evals + 1)
 
     ν₁ = α * Δk / (1 + λmax * (α * Δk + 1))
+    set_solver_specific!(stats, :sigma_cauchy, 1/ν₁)
     @. mν∇fk = -ν₁ * ∇fk
 
     prox!(s, ψ, mν∇fk, ν₁)
+    set_solver_specific!(stats, :scp_norm, norm(s))
     ξ1 = hk - mk1(s) + max(1, abs(hk)) * 10 * eps()
     sqrt_ξ1_νInv = sqrt(ξ1 / ν₁)
 
