Detect unbounded and infeasible problems

Author: BenjaminPINEAU

src/AL_alg.jl

@@ -102,7 +102,12 @@ If adopted, the Hessian is accessed as an abstract operator and need not be the
 
 - `x::AbstractVector`: a primal initial guess (default: `reg_nlp.model.meta.x0`)
 - `y::AbstractVector`: a dual initial guess (default: `reg_nlp.model.meta.y0`)
-- `atol::T = √eps(T)`: absolute optimality tolerance;
+- `atol::T = eps(T)^(1/3)`: absolute tolerance
+- `diverging_iterates_tol::T = eps(T)^(-1)`: diverging tolerance for the norm of the iterates (the norm should be lower than the tolerance);
+- `diverging_obj_tol::T = -eps(T)^(-1)`: diverging tolerance for the objective function (the objective function should be higher than the tolerance);
+- `cviol_tol::T = eps(T)^(-1)`: tolerance to determine whether the constraints are infeasible
+- `diverging_max_iter::Int = 5`: maximum number of iteration at which `diverging_obj_tol` or `diverging_iterates_tol` is violated;
+- `cviol_max_iter::Int = 5`: maximum number of iteration at which the regularisation parameter is increasing and the constraints are still violated;
 - `ctol::T = atol`: absolute feasibility tolerance;
 - `verbose::Int = 0`: if > 0, display iteration details every `verbose` iteration;
 - `max_iter::Int = 10000`: maximum number of iterations;
@@ -209,7 +214,12 @@ function SolverCore.solve!(
   callback = (args...) -> nothing,
   x::V = reg_nlp.model.meta.x0,
   y::V = reg_nlp.model.meta.y0,
-  atol::T = √eps(T),
+  atol::T = eps(T)^(1/3),
+  diverging_iterates_tol::T = eps(T)^(-1),
+  diverging_obj_tol::T = -eps(T)^(-1),
+  cviol_tol::T = eps(T)^(-1),
+  diverging_max_iter::Int = 5,
+  cviol_max_iter::Int = 5,
   verbose::Int = 0,
   max_iter::Int = 10000,
   max_time::Float64 = 30.0,
@@ -225,6 +235,9 @@ function SolverCore.solve!(
 ) where {T, V}
   reset!(stats)
 
+  local diverging_iter::Int = zero(Int)
+  local cviol_iter::Int = zero(Int)
+
   # Retrieve workspace
   nlp = reg_nlp.model
   h = reg_nlp.h
@@ -315,6 +328,8 @@ function SolverCore.solve!(
     # objective
     fx = obj(nlp, solver.x)
     hx = @views h(solver.x[selected])
+    improper = (hx == -Inf)
+
     objx = fx + hx
     set_objective!(stats, objx)
     set_solver_specific!(stats, :smooth_obj, fx)
@@ -345,19 +360,19 @@ function SolverCore.solve!(
     set_time!(stats, time() - start_time)
     set_status!(
       stats,
-      SolverCore.get_status(
-        nlp,
+      get_status(
+        reg_nlp;
         elapsed_time = stats.elapsed_time,
         iter = stats.iter,
         optimal = optimal,
-        infeasible = false,
-        parameter_too_large = false,
-        unbounded = false,
-        stalled = false,
-        exception = false,
+        improper = improper,
+        diverging_iter = diverging_iter,
+        cviol_iter = cviol_iter,
         max_eval = max_eval,
         max_time = max_time,
         max_iter = max_iter,
+        diverging_max_iter = diverging_max_iter,
+        cviol_max_iter = cviol_max_iter,
       ),
     )
 
@@ -372,6 +387,13 @@ function SolverCore.solve!(
     if !done
       if cviol > max(ctol, factor_primal_linear_improvement * cviol_old)
         mu *= factor_penalty_up
+        if cviol > cviol_tol
+          cviol_iter += 1
+        end
+      end
+      if (fx + hx < diverging_obj_tol) || (norm(solver.x) > diverging_iterates_tol)
+        mu *= factor_penalty_up
+        diverging_iter = diverging_iter + 1
       end
       update_μ!(solver.sub_problem.model, mu)
       cviol_old = cviol

src/R2N.jl

@@ -137,8 +137,11 @@ For advanced usage, first define a solver "R2NSolver" to preallocate the memory
 
 # Keyword arguments 
 - `x::V = nlp.meta.x0`: the initial guess;
-- `atol::T = √eps(T)`: absolute tolerance;
-- `rtol::T = √eps(T)`: relative tolerance;
+- `atol::T = eps(T)^(1/3)`: absolute tolerance
+- `diverging_iterates_tol::T = eps(T)^(-1)`: diverging tolerance for the norm of the iterates (the norm should be lower than the tolerance);
+- `diverging_obj_tol::T = -eps(T)^(-1)`: diverging tolerance for the objective function (the objective function should be higher than the tolerance);
+- `diverging_max_iter::Int = 5`: maximum number of iteration at which `diverging_obj_tol` or `diverging_iterates_tol` is violated;
+- `rtol::T = eps(T)^(1/3)`: relative tolerance;
 - `neg_tol::T = eps(T)^(1 / 4)`: negative tolerance;
 - `max_eval::Int = -1`: maximum number of evaluation of the objective function (negative number means unlimited);
 - `max_time::Float64 = 30.0`: maximum time limit in seconds;
@@ -216,8 +219,11 @@ function SolverCore.solve!(
   qn_update_y!::Function = _qn_grad_update_y!,
   qn_copy!::Function = _qn_grad_copy!,
   x::V = reg_nlp.model.meta.x0,
-  atol::T = √eps(T),
-  rtol::T = √eps(T),
+  atol::T = eps(T)^(1/3),
+  diverging_iterates_tol::T = eps(T)^(-1),
+  diverging_obj_tol::T = -eps(T)^(-1),
+  diverging_max_iter::Int = 5,
+  rtol::T = eps(T)^(1/3),
   neg_tol::T = eps(T)^(1 / 4),
   verbose::Int = 0,
   max_iter::Int = 10000,
@@ -299,6 +305,7 @@ function SolverCore.solve!(
   local ξ1::T
   local ρk::T = zero(T)
   local prox_evals::Int = 0
+  local diverging_iter::Int = zero(Int)
 
   fk = compute_obj ? obj(nlp, xk) : stats.solver_specific[:smooth_obj]
   compute_grad && grad!(nlp, xk, ∇fk)
@@ -368,9 +375,11 @@ function SolverCore.solve!(
       iter = stats.iter,
       optimal = solved,
       improper = improper,
+      diverging_iter = diverging_iter,
       max_eval = max_eval,
       max_time = max_time,
       max_iter = max_iter,
+      diverging_max_iter = diverging_max_iter,
     ),
   )
 
@@ -475,13 +484,18 @@ function SolverCore.solve!(
       set_step_status!(stats, :accepted)
     end
 
-    if η2 ≤ ρk < Inf
-      σk = max(σk/γ, σmin)
-    end
-
-    if ρk < η1 || ρk == Inf
+    if (fk + hk < diverging_obj_tol) || (norm(xk) > diverging_iterates_tol)
       σk = σk * γ
-      set_step_status!(stats, :rejected)
+      diverging_iter = diverging_iter + 1
+    else 
+      diverging_iter = 0
+      if η2 ≤ ρk < Inf
+        σk = max(σk / γ, σmin)
+      end
+      if ρk < η1 || ρk == Inf
+        σk = σk * γ
+        set_step_status!(stats, :rejected)
+      end
     end
 
     ν₁ = θ / (λmax + σk)
@@ -515,9 +529,11 @@ function SolverCore.solve!(
         iter = stats.iter,
         optimal = solved,
         improper = improper,
+        diverging_iter = diverging_iter,
         max_eval = max_eval,
         max_time = max_time,
         max_iter = max_iter,
+        diverging_max_iter = diverging_max_iter,
       ),
     )
 

src/R2_alg.jl

@@ -141,9 +141,12 @@ For advanced usage, first define a solver "R2Solver" to preallocate the memory u
 
 # Keyword arguments 
 - `x::V = nlp.meta.x0`: the initial guess;
-- `atol::T = √eps(T)`: absolute tolerance;
-- `rtol::T = √eps(T)`: relative tolerance;
-- `neg_tol::T = eps(T)^(1 / 4)`: negative tolerance
+- `atol::T = eps(T)^(1/3)`: absolute tolerance
+- `diverging_iterates_tol::T = eps(T)^(-1)`: diverging tolerance for the norm of the iterates (the norm should be lower than the tolerance);
+- `diverging_obj_tol::T = -eps(T)^(-1)`: diverging tolerance for the objective function (the objective function should be higher than the tolerance);
+- `diverging_max_iter::Int = 5`: maximum number of iteration at which `diverging_obj_tol` or `diverging_iterates_tol` is violated;
+- `rtol::T = eps(T)^(1/3)`: relative tolerance;
+- `neg_tol::T = eps(T)^(1 / 4)`: negative tolerance;
 - `max_eval::Int = -1`: maximum number of evaluation of the objective function (negative number means unlimited);
 - `max_time::Float64 = 30.0`: maximum time limit in seconds;
 - `max_iter::Int = 10000`: maximum number of iterations;
@@ -313,8 +316,11 @@ function SolverCore.solve!(
   stats::GenericExecutionStats{T, V};
   callback = (args...) -> nothing,
   x::V = reg_nlp.model.meta.x0,
-  atol::T = √eps(T),
-  rtol::T = √eps(T),
+  atol::T = eps(T)^(1/3),
+  diverging_iterates_tol::T = eps(T)^(-1),
+  diverging_obj_tol::T = -eps(T)^(-1),
+  diverging_max_iter::Int = 5,
+  rtol::T = eps(T)^(1/3),
   neg_tol::T = eps(T)^(1 / 4),
   verbose::Int = 0,
   max_iter::Int = 10000,
@@ -386,6 +392,7 @@ function SolverCore.solve!(
 
   local ξ::T
   local ρk::T = zero(T)
+  local diverging_iter::Int = zero(Int)
   σk = max(1 / ν, σmin)
   ν = 1 / σk
   sqrt_ξ_νInv = one(T)
@@ -420,14 +427,16 @@ function SolverCore.solve!(
   set_status!(
     stats,
     get_status(
-      reg_nlp,
+      reg_nlp;
       elapsed_time = stats.elapsed_time,
       iter = stats.iter,
       optimal = solved,
       improper = improper,
+      diverging_iter = diverging_iter,
       max_eval = max_eval,
       max_time = max_time,
       max_iter = max_iter,
+      diverging_max_iter = diverging_max_iter,
     ),
   )
 
@@ -476,12 +485,18 @@ function SolverCore.solve!(
       set_step_status!(stats, :accepted)
     end
 
-    if η2 ≤ ρk < Inf
-      σk = max(σk / γ, σmin)
-    end
-    if ρk < η1 || ρk == Inf
+    if (fk + hk < diverging_obj_tol) || (norm(xk) > diverging_iterates_tol)
       σk = σk * γ
-      set_step_status!(stats, :rejected)
+      diverging_iter = diverging_iter + 1
+    else
+      diverging_iter = 0
+      if η2 ≤ ρk < Inf
+        σk = max(σk / γ, σmin)
+      end
+      if ρk < η1 || ρk == Inf
+        σk = σk * γ
+        set_step_status!(stats, :rejected)
+      end
     end
 
     ν = 1 / σk
@@ -506,14 +521,16 @@ function SolverCore.solve!(
     set_status!(
       stats,
       get_status(
-        reg_nlp,
+        reg_nlp;
         elapsed_time = stats.elapsed_time,
         iter = stats.iter,
         optimal = solved,
         improper = improper,
+        diverging_iter = diverging_iter,
         max_eval = max_eval,
         max_time = max_time,
         max_iter = max_iter,
+        diverging_max_iter = diverging_max_iter,
       ),
     )
 

src/utils.jl

@@ -142,9 +142,13 @@ function get_status(
   iter = 0,
   optimal = false,
   improper = false,
+  diverging_iter = 0,
+  cviol_iter = 0,
   max_eval = Inf,
   max_time = Inf,
   max_iter = Inf,
+  diverging_max_iter = Inf,
+  cviol_max_iter = Inf,
 ) where {M <: AbstractRegularizedNLPModel}
   if optimal
     :first_order
@@ -156,6 +160,10 @@ function get_status(
     :max_time
   elseif neval_obj(reg_nlp.model) >= max_eval && max_eval >= 0
     :max_eval
+  elseif diverging_max_iter < diverging_iter
+    :unbounded
+  elseif cviol_max_iter < cviol_iter
+    :infeasible
   else
     :unknown
   end