Update criticality measure (#117)

Co-authored-by: Dominique <dominique.orban@gmail.com>

Author: geoffroyleconte

benchmarks/tables/fh-table.jl

@@ -13,9 +13,9 @@ uvar = fill(Inf, 5)
 data, simulate, resid, misfit, x0 = RegularizedProblems.FH_smooth_term()
 model = ADNLPModel(misfit, ones(5), lvar, uvar)
 f = LBFGSModel(model)
-λ = cstr ? 2.0e1 : 1.0e1
+λ = cstr ? 4.0e1 : 1.0e1
 
-h = NormL1(λ)
+h = cstr ? NormL1(λ) : NormL0(λ)
 ν = 1.0e0
 verbose = 0 #10
 maxIter = 500

benchmarks/tables/regulopt-tables.jl

@@ -107,7 +107,7 @@ function benchmark_table(
         "solver",
         L"$f(x)$",
         L"$h(x) / \lambda$",
-        L"$\xi$",
+        L"$\sqrt{\xi / \nu}$",
         L"$\# \ f$",
         L"$\# \ \nabla f$",
         L"$\# \ prox$",
@@ -118,7 +118,7 @@ function benchmark_table(
         "solver",
         "\$f(x)\$",
         L"$h(x)/\lambda$",
-        L"$\xi$",
+        L"$\sqrt{\xi / \nu}$",
         pb_name[1:3] == "SVM" ? L"$(Train, Test)$" : L"$\|x-x_T\|_2$",
         L"$\# \ f$",
         L"$\# \ \nabla f$",
@@ -128,13 +128,13 @@ function benchmark_table(
     end
   else
     if length(sol) == 0
-      header = ["solver", "f(x)", "h(x)/λ", "ξ", "# f", "# ∇f", "# prox", "t (s)"]
+      header = ["solver", "f(x)", "h(x)/λ", "√(ξ/ν)", "# f", "# ∇f", "# prox", "t (s)"]
     else
       header = [
         "solver",
         "f(x)",
         "h(x)/λ",
-        "ξ",
+        "√ξ/√ν",
         pb_name[1:3] == "SVM" ? "(Train, Test)" : "||x-x*||",
         "# f",
         "# ∇f",

benchmarks/tables/svm-table.jl

@@ -15,7 +15,7 @@ verbose = 0 #10
 ϵ = 1.0e-4
 ϵi = 1.0e-3
 ϵri = 1.0e-6
-maxIter = 500
+maxIter = 1000
 maxIter_inner = 100
 options =
   ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true)

src/R2_alg.jl

@@ -117,7 +117,7 @@ function R2(nlp::AbstractNLPModel, args...; kwargs...)
   set_status!(stats, outdict[:status])
   set_solution!(stats, xk)
   set_objective!(stats, outdict[:fk] + outdict[:hk])
-  set_residuals!(stats, zero(eltype(xk)), ξ ≥ 0 ? sqrt(ξ) : ξ)
+  set_residuals!(stats, zero(eltype(xk)), ξ)
   set_iter!(stats, k)
   set_time!(stats, outdict[:elapsed_time])
   set_solver_specific!(stats, :Fhist, outdict[:Fhist])
@@ -258,14 +258,15 @@ function R2!(
 
   if verbose > 0
     #! format: off
-    @info @sprintf "%6s %8s %8s %7s %8s %7s %7s %7s %1s" "iter" "f(x)" "h(x)" "√ξ" "ρ" "σ" "‖x‖" "‖s‖" ""
+    @info @sprintf "%6s %8s %8s %7s %8s %7s %7s %7s %1s" "iter" "f(x)" "h(x)" "√(ξ/ν)" "ρ" "σ" "‖x‖" "‖s‖" ""
     #! format: off
   end
 
   local ξ::R
   k = 0
   σk = max(1 / ν, σmin)
   ν = 1 / σk
+  sqrt_ξ_νInv = one(R)
 
   fk = f(xk)
   ∇f!(∇fk, xk)
@@ -288,12 +289,13 @@ function R2!(
     Complex_hist[k] += 1
     mks = mk(s)
     ξ = hk - mks + max(1, abs(hk)) * 10 * eps()
+    sqrt_ξ_νInv = ξ ≥ 0 ? sqrt(ξ / ν) : sqrt(-ξ / ν)
 
     if ξ ≥ 0 && k == 1
-      ϵ += ϵr * sqrt(ξ/ν)  # make stopping test absolute and relative
+      ϵ += ϵr * sqrt_ξ_νInv # make stopping test absolute and relative
     end
 
-    if (ξ < 0 && sqrt(-ξ/ν) ≤ neg_tol) || (ξ ≥ 0 && sqrt(ξ/ν) ≤ ϵ)
+    if (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv ≤ ϵ)
       optimal = true
       continue
     end
@@ -310,7 +312,7 @@ function R2!(
     if (verbose > 0) && (k % ptf == 0)
       #! format: off
       σ_stat = (η2 ≤ ρk < Inf) ? "↘" : (ρk < η1 ? "↗" : "=")
-      @info @sprintf "%6d %8.1e %8.1e %7.1e %8.1e %7.1e %7.1e %7.1e %1s" k fk hk sqrt(ξ/ν) ρk σk norm(xk) norm(s) σ_stat
+      @info @sprintf "%6d %8.1e %8.1e %7.1e %8.1e %7.1e %7.1e %7.1e %1s" k fk hk sqrt_ξ_νInv ρk σk norm(xk) norm(s) σ_stat
 
       #! format: on
     end
@@ -350,7 +352,7 @@ function R2!(
       #! format: off
       @info @sprintf "%6d %8.1e %8.1e %7.1e %8s %7.1e %7.1e %7.1e" k fk hk sqrt(ξ/ν) "" σk norm(xk) norm(s)
       #! format: on
-      @info "R2: terminating with √ξ/√ν = $(sqrt(ξ/ν))"
+      @info "R2: terminating with √(ξ/ν) = $(sqrt_ξ_νInv)"
     end
   end
 
@@ -364,5 +366,5 @@ function R2!(
     :exception
   end
 
-  return k, status, fk, hk, ξ
+  return k, status, fk, hk, sqrt_ξ_νInv
 end

src/TRDH_alg.jl

@@ -72,7 +72,7 @@ function TRDH(
   set_status!(stats, outdict[:status])
   set_solution!(stats, xk)
   set_objective!(stats, outdict[:fk] + outdict[:hk])
-  set_residuals!(stats, zero(eltype(xk)), ξ ≥ 0 ? sqrt(ξ) : ξ)
+  set_residuals!(stats, zero(eltype(xk)), ξ)
   set_iter!(stats, k)
   set_time!(stats, outdict[:elapsed_time])
   set_solver_specific!(stats, :Fhist, outdict[:Fhist])
@@ -188,9 +188,9 @@ function TRDH(
   if verbose > 0
     #! format: off
     if reduce_TR
-      @info @sprintf "%6s %8s %8s %7s %7s %8s %7s %7s %7s %7s %1s" "outer" "f(x)" "h(x)" "√ξ1" "√ξ" "ρ" "Δ" "‖x‖" "‖s‖" "‖Dₖ‖" "TRDH"
+      @info @sprintf "%6s %8s %8s %7s %7s %8s %7s %7s %7s %7s %1s" "outer" "f(x)" "h(x)" "√(ξ1/ν)" "√ξ" "ρ" "Δ" "‖x‖" "‖s‖" "‖Dₖ‖" "TRDH"
     else
-      @info @sprintf "%6s %8s %8s %7s %8s %7s %7s %7s %7s %1s" "outer" "f(x)" "h(x)" "√ξ" "ρ" "Δ" "‖x‖" "‖s‖" "‖Dₖ‖" "TRDH"
+      @info @sprintf "%6s %8s %8s %7s %8s %7s %7s %7s %7s %1s" "outer" "f(x)" "h(x)" "√(ξ/ν)" "ρ" "Δ" "‖x‖" "‖s‖" "‖Dₖ‖" "TRDH"
     end
     #! format: off
   end
@@ -209,6 +209,7 @@ function TRDH(
   νInv = (DkNorm + one(R) / (α * Δk))
   ν = one(R) / νInv
   mν∇fk = -ν .* ∇fk
+  sqrt_ξ_νInv = one(R)
 
   optimal = false
   tired = maxIter > 0 && k ≥ maxIter || elapsed_time > maxTime
@@ -227,12 +228,13 @@ function TRDH(
       prox!(s, ψ, mν∇fk, ν)
       Complex_hist[k] += 1
       ξ1 = hk - mk1(s) + max(1, abs(hk)) * 10 * eps()
+      sqrt_ξ_νInv = ξ1 ≥ 0 ? sqrt(ξ1 / ν) : sqrt(-ξ1 / ν)
 
       if ξ1 ≥ 0 && k == 1
-        ϵ += ϵr * sqrt(ξ1)  # make stopping test absolute and relative
+        ϵ += ϵr * sqrt_ξ_νInv  # make stopping test absolute and relative
       end
 
-      if (ξ1 < 0 && sqrt(-ξ1) ≤ neg_tol) || (ξ1 ≥ 0 && sqrt(ξ1) < ϵ)
+      if (ξ1 < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ_νInv < ϵ)
         # the current xk is approximately first-order stationary
         optimal = true
         continue
@@ -267,11 +269,13 @@ function TRDH(
     ξ = hk - mk(s) + max(1, abs(hk)) * 10 * eps()
 
     if !reduce_TR
+
+      sqrt_ξ_νInv = ξ ≥ 0 ? sqrt(ξ / ν) : sqrt(-ξ / ν)
       if ξ ≥ 0 && k == 1
-        ϵ += ϵr * sqrt(ξ)  # make stopping test absolute and relative
+        ϵ += ϵr * sqrt_ξ_νInv  # make stopping test absolute and relative
       end
 
-      if (ξ < 0 && sqrt(-ξ) ≤ neg_tol) || (ξ ≥ 0 && sqrt(ξ) < ϵ)
+      if (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv < ϵ)
         # the current xk is approximately first-order stationary
         optimal = true
         continue
@@ -289,9 +293,9 @@ function TRDH(
     if (verbose > 0) && (k % ptf == 0)
       #! format: off
       if reduce_TR
-        @info @sprintf "%6d %8.1e %8.1e %7.1e %7.1e %8.1e %7.1e %7.1e %7.1e %7.1e %1s" k fk hk sqrt(ξ1) sqrt(ξ) ρk Δk χ(xk) sNorm norm(Dk.d) TR_stat
+        @info @sprintf "%6d %8.1e %8.1e %7.1e %7.1e %8.1e %7.1e %7.1e %7.1e %7.1e %1s" k fk hk sqrt_ξ_νInv sqrt(ξ) ρk Δk χ(xk) sNorm norm(Dk.d) TR_stat
       else
-        @info @sprintf "%6d %8.1e %8.1e %7.1e %8.1e %7.1e %7.1e %7.1e %7.1e %1s" k fk hk sqrt(ξ) ρk Δk χ(xk) sNorm norm(Dk.d) TR_stat
+        @info @sprintf "%6d %8.1e %8.1e %7.1e %8.1e %7.1e %7.1e %7.1e %7.1e %1s" k fk hk sqrt_ξ_νInv ρk Δk χ(xk) sNorm norm(Dk.d) TR_stat
       end
       #! format: on
     end
@@ -321,7 +325,7 @@ function TRDH(
       has_bnds ? set_bounds!(ψ, l_bound_k, u_bound_k) : set_radius!(ψ, Δk)
     end
 
-    νInv = (DkNorm + one(R) / (α * Δk))
+    νInv = reduce_TR ? (DkNorm + one(R) / (α * Δk)) : (DkNorm + one(R) / α)
     ν = one(R) / νInv
     mν∇fk .= -ν .* ∇fk
 
@@ -334,14 +338,12 @@ function TRDH(
     elseif optimal
       #! format: off
       if reduce_TR
-        @info @sprintf "%6d %8.1e %8.1e %7.1e %7.1e %8s %7.1e %7.1e %7.1e %7.1e" k fk hk sqrt(ξ1) sqrt(ξ1) "" Δk χ(xk) χ(s) norm(Dk.d)
+        @info @sprintf "%6d %8.1e %8.1e %7.1e %7.1e %8s %7.1e %7.1e %7.1e %7.1e" k fk hk sqrt_ξ_νInv sqrt(ξ1) "" Δk χ(xk) χ(s) norm(Dk.d)
         #! format: on
-        @info "TRDH: terminating with √ξ1 = $(sqrt(ξ1))"
+        @info "TRDH: terminating with √(ξ1/ν) = $(sqrt_ξ_νInv)"
       else
-        @info @sprintf "%6d %8.1e %8.1e %7.1e %8s %7.1e %7.1e %7.1e %7.1e" k fk hk sqrt(ξ) "" Δk χ(
-          xk,
-        ) χ(s) norm(Dk.d)
-        @info "TRDH: terminating with √ξ = $(sqrt(ξ))"
+        @info @sprintf "%6d %8.1e %8.1e %7.1e %8s %7.1e %7.1e %7.1e %7.1e" k fk hk sqrt_ξ_νInv "" Δk χ(xk) χ(s) norm(Dk.d)
+        @info "TRDH: terminating with √(ξ/ν) = $(sqrt_ξ_νInv)"
       end
     end
   end
@@ -365,7 +367,7 @@ function TRDH(
     :status => status,
     :fk => fk,
     :hk => hk,
-    :ξ => ξ1,
+    :ξ => sqrt_ξ_νInv,
     :elapsed_time => elapsed_time,
   )
 

src/TR_alg.jl

@@ -48,13 +48,13 @@ function TR(
   f::AbstractNLPModel,
   h::H,
   χ::X,
-  options::ROSolverOptions;
-  x0::AbstractVector = f.meta.x0,
+  options::ROSolverOptions{R};
+  x0::AbstractVector{R} = f.meta.x0,
   subsolver_logger::Logging.AbstractLogger = Logging.NullLogger(),
   subsolver = R2,
   subsolver_options = ROSolverOptions(ϵa = options.ϵa),
   selected::AbstractVector{<:Integer} = 1:(f.meta.nvar),
-) where {H, X}
+) where {H, X, R}
   start_time = time()
   elapsed_time = 0.0
   # initialize passed options
@@ -117,7 +117,7 @@ function TR(
   Complex_hist = zeros(Int, maxIter)
   if verbose > 0
     #! format: off
-    @info @sprintf "%6s %8s %8s %8s %7s %7s %8s %7s %7s %7s %7s %1s" "outer" "inner" "f(x)" "h(x)" "√ξ1" "√ξ" "ρ" "Δ" "‖x‖" "‖s‖" "‖Bₖ‖" "TR"
+    @info @sprintf "%6s %8s %8s %8s %7s %7s %8s %7s %7s %7s %7s %1s" "outer" "inner" "f(x)" "h(x)" "√(ξ1/ν)" "√ξ" "ρ" "Δ" "‖x‖" "‖s‖" "‖Bₖ‖" "TR"
     #! format: on
   end
 
@@ -133,6 +133,8 @@ function TR(
 
   λmax = opnorm(Bk)
   νInv = (1 + θ) * λmax
+  ν = 1 / (νInv + 1 / (Δk * α))
+  sqrt_ξ1_νInv = one(R)
 
   optimal = false
   tired = k ≥ maxIter || elapsed_time > maxTime
@@ -160,24 +162,24 @@ function TR(
 
     # Take first proximal gradient step s1 and see if current xk is nearly stationary.
     # s1 minimizes φ1(s) + ‖s‖² / 2 / ν + ψ(s) ⟺ s1 ∈ prox{νψ}(-ν∇φ1(0)).
-    ν = 1 / (νInv + 1 / (Δk * α))
     prox!(s, ψ, -ν * ∇fk, ν)
     ξ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 / ν)
 
     if ξ1 ≥ 0 && k == 1
-      ϵ_increment = ϵr * sqrt(ξ1)
+      ϵ_increment = ϵr * sqrt_ξ1_νInv
       ϵ += ϵ_increment  # make stopping test absolute and relative
       ϵ_subsolver += ϵ_increment
     end
 
-    if sqrt(ξ1) < ϵ
+    if sqrt_ξ1_νInv < ϵ
       # the current xk is approximately first-order stationary
       optimal = true
       continue
     end
 
-    subsolver_options.ϵa = k == 1 ? 1.0e-5 : max(ϵ_subsolver, min(1e-2, sqrt(ξ1)) * ξ1)
+    subsolver_options.ϵa = k == 1 ? 1.0e-5 : max(ϵ_subsolver, min(1e-2, sqrt_ξ1_νInv))
     ∆_effective = min(β * χ(s), Δk)
     (has_bounds(f) || subsolver == TRDH) ?
     set_bounds!(ψ, max.(-∆_effective, l_bound - xk), min.(∆_effective, u_bound - xk)) :
@@ -214,7 +216,7 @@ function TR(
 
     if (verbose > 0) && (k % ptf == 0)
       #! format: off
-      @info @sprintf "%6d %8d %8.1e %8.1e %7.1e %7.1e %8.1e %7.1e %7.1e %7.1e %7.1e %1s" k iter fk hk sqrt(ξ1) sqrt(ξ) ρk ∆_effective χ(xk) sNorm νInv TR_stat
+      @info @sprintf "%6d %8d %8.1e %8.1e %7.1e %7.1e %8.1e %7.1e %7.1e %7.1e %7.1e %1s" k iter fk hk sqrt_ξ1_νInv sqrt(ξ) ρk ∆_effective χ(xk) sNorm νInv TR_stat
       #! format: on
     end
 
@@ -248,6 +250,7 @@ function TR(
       (has_bounds(f) || subsolver == TRDH) ?
       set_bounds!(ψ, max.(-Δk, l_bound - xk), min.(Δk, u_bound - xk)) : set_radius!(ψ, Δk)
     end
+    ν = 1 / (νInv + 1 / (Δk * α))
     tired = k ≥ maxIter || elapsed_time > maxTime
   end
 
@@ -256,9 +259,9 @@ function TR(
       @info @sprintf "%6d %8s %8.1e %8.1e" k "" fk hk
     elseif optimal
       #! format: off
-      @info @sprintf "%6d %8d %8.1e %8.1e %7.1e %7.1e %8s %7.1e %7.1e %7.1e %7.1e" k 1 fk hk sqrt(ξ1) sqrt(ξ1) "" Δk χ(xk) χ(s) νInv
+      @info @sprintf "%6d %8d %8.1e %8.1e %7.1e %7.1e %8s %7.1e %7.1e %7.1e %7.1e" k 1 fk hk sqrt_ξ1_νInv sqrt(ξ1) "" Δk χ(xk) χ(s) νInv
       #! format: on
-      @info "TR: terminating with √ξ1 = $(sqrt(ξ1))"
+      @info "TR: terminating with √(ξ1/ν) = $(sqrt_ξ1_νInv)"
     end
   end
 
@@ -276,7 +279,7 @@ function TR(
   set_status!(stats, status)
   set_solution!(stats, xk)
   set_objective!(stats, fk + hk)
-  set_residuals!(stats, zero(eltype(xk)), ξ1 ≥ 0 ? sqrt(ξ1) : ξ1)
+  set_residuals!(stats, zero(eltype(xk)), sqrt_ξ1_νInv)
   set_iter!(stats, k)
   set_time!(stats, elapsed_time)
   set_solver_specific!(stats, :Fhist, Fobj_hist[1:k])

src/input_struct.jl

@@ -23,7 +23,7 @@ mutable struct ROSolverOptions{R}
   function ROSolverOptions{R}(;
     ϵa::R = √eps(R),
     ϵr::R = √eps(R),
-    neg_tol::R = eps(R)^(1 / 3),
+    neg_tol::R = eps(R)^(1 / 4),
     Δk::R = one(R),
     verbose::Int = 0,
     maxIter::Int = 10000,