Merge branch 'JuliaSmoothOptimizers:master' into documentation

Author: MaxenceGollier

Project.toml

@@ -19,11 +19,13 @@ ShiftedProximalOperators = "d4fd37fa-580c-4e43-9b30-361c21aae263"
 SolverCore = "ff4d7338-4cf1-434d-91df-b86cb86fb843"
 
 [compat]
+ADNLPModels = "0.8.13"
 Arpack = "0.5"
 LinearOperators = "2.10.0"
 ManualNLPModels = "0.2.0"
 NLPModels = "0.19, 0.20, 0.21"
 NLPModelsModifiers = "0.7"
+OptimizationProblems = "0.9.2"
 Percival = "0.7.2"
 ProximalOperators = "0.15"
 RegularizedProblems = "0.1.1"
@@ -32,10 +34,12 @@ SolverCore = "0.3.0"
 julia = "^1.6.0"
 
 [extras]
+ADNLPModels = "54578032-b7ea-4c30-94aa-7cbd1cce6c9a"
+OptimizationProblems = "5049e819-d29b-5fba-b941-0eee7e64c1c6"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RegularizedProblems = "ea076b23-609f-44d2-bb12-a4ae45328278"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestSetExtensions = "98d24dd4-01ad-11ea-1b02-c9a08f80db04"
 
 [targets]
-test = ["Random", "RegularizedProblems", "Test", "TestSetExtensions"]
+test = ["ADNLPModels", "OptimizationProblems", "Random", "RegularizedProblems", "Test", "TestSetExtensions"]

src/AL_alg.jl

@@ -138,14 +138,14 @@ Notably, you can access, and modify, the following:
   - `stats.solver_specific[:smooth_obj]`: current value of the smooth part of the objective function;
   - `stats.solver_specific[:nonsmooth_obj]`: current value of the nonsmooth part of the objective function.
 """
-mutable struct ALSolver{T, V, M, ST} <: AbstractOptimizationSolver
+mutable struct ALSolver{T, V, M, Pb, ST} <: AbstractOptimizationSolver
   x::V
   cx::V
   y::V
   has_bnds::Bool
-  sub_model::AugLagModel{M, T, V}
+  sub_problem::Pb
   sub_solver::ST
-  sub_stats::GenericExecutionStats{T, V, V, Any}
+  sub_stats::GenericExecutionStats{T, V, V, T}
 end
 
 function ALSolver(reg_nlp::AbstractRegularizedNLPModel{T, V}; kwargs...) where {T, V}
@@ -155,12 +155,21 @@ function ALSolver(reg_nlp::AbstractRegularizedNLPModel{T, V}; kwargs...) where {
   cx = V(undef, ncon)
   y = V(undef, ncon)
   has_bnds = has_bounds(nlp)
-  sub_model = AugLagModel(nlp, V(undef, ncon), T(0), x, T(0), V(undef, ncon))
+  sub_model = AugLagModel(nlp, V(undef, ncon), T(0), x, T(0), cx)
+  sub_problem = RegularizedNLPModel(sub_model, reg_nlp.h, reg_nlp.selected)
   sub_solver = R2Solver(reg_nlp; kwargs...)
-  sub_stats = GenericExecutionStats(sub_model)
+  sub_stats = RegularizedExecutionStats(sub_problem)
   M = typeof(nlp)
   ST = typeof(sub_solver)
-  return ALSolver{T, V, M, ST}(x, cx, y, has_bnds, sub_model, sub_solver, sub_stats)
+  return ALSolver{T, V, M, typeof(sub_problem), ST}(
+    x,
+    cx,
+    y,
+    has_bnds,
+    sub_problem,
+    sub_solver,
+    sub_stats,
+  )
 end
 
 @doc (@doc ALSolver) function AL(::Val{:equ}, reg_nlp::AbstractRegularizedNLPModel; kwargs...)
@@ -182,7 +191,7 @@ function SolverCore.solve!(
   model::AbstractRegularizedNLPModel;
   kwargs...,
 )
-  stats = GenericExecutionStats(model.model)
+  stats = RegularizedExecutionStats(model)
   solve!(solver, model, stats; kwargs...)
 end
 
@@ -209,6 +218,7 @@ function SolverCore.solve!(
   factor_decrease_subtol::T = T(1 // 4),
   dual_safeguard = project_y!,
 ) where {T, V}
+  reset!(stats)
 
   # Retrieve workspace
   nlp = reg_nlp.model
@@ -254,8 +264,8 @@ function SolverCore.solve!(
   set_solver_specific!(stats, :nonsmooth_obj, hx)
 
   mu = init_penalty
-  solver.sub_model.y .= solver.y
-  update_μ!(solver.sub_model, mu)
+  solver.sub_problem.model.y .= solver.y
+  update_μ!(solver.sub_problem.model, mu)
 
   cviol = norm(solver.cx, Inf)
   cviol_old = Inf
@@ -279,15 +289,13 @@ function SolverCore.solve!(
     iter += 1
 
     # dual safeguard
-    dual_safeguard(solver.sub_model)
+    dual_safeguard(solver.sub_problem.model)
 
-    # AL subproblem
-    sub_reg_nlp = RegularizedNLPModel(solver.sub_model, h, selected)
     subtol = max(subtol, atol)
     reset!(subout)
     solve!(
       solver.sub_solver,
-      sub_reg_nlp,
+      solver.sub_problem,
       subout,
       x = solver.x,
       atol = subtol,
@@ -298,8 +306,8 @@ function SolverCore.solve!(
       verbose = subsolver_verbose,
     )
     solver.x .= subout.solution
-    solver.cx .= solver.sub_model.cx
-    subiters += subout.iter
+    solver.cx .= solver.sub_problem.model.cx
+    subiters = subout.iter
 
     # objective
     fx = obj(nlp, solver.x)
@@ -310,8 +318,8 @@ function SolverCore.solve!(
     set_solver_specific!(stats, :nonsmooth_obj, hx)
 
     # dual estimate
-    update_y!(solver.sub_model)
-    solver.y .= solver.sub_model.y
+    update_y!(solver.sub_problem.model)
+    solver.y .= solver.sub_problem.model.y
     set_constraint_multipliers!(stats, solver.y)
 
     # stationarity measure
@@ -362,7 +370,7 @@ function SolverCore.solve!(
       if cviol > max(ctol, factor_primal_linear_improvement * cviol_old)
         mu *= factor_penalty_up
       end
-      update_μ!(solver.sub_model, mu)
+      update_μ!(solver.sub_problem.model, mu)
       cviol_old = cviol
       subtol *= factor_decrease_subtol
       rem_eval = max_eval < 0 ? max_eval : max_eval - neval_obj(nlp)

src/LMModel.jl

@@ -0,0 +1,59 @@
+export LMModel
+
+@doc raw"""
+    LMModel(j_prod!, jt_prod, F, v, σ, xk)
+
+Given the unconstrained optimization problem:
+```math
+\min \tfrac{1}{2} \| F(x) \|^2,
+```
+this model represents the smooth LM subproblem:
+```math
+\min_s \ \tfrac{1}{2} \| F(x) + J(x)s \|^2 + \tfrac{1}{2} σ \|s\|^2
+```
+where `J` is the Jacobian of `F` at `xk`, represented via matrix-free operations.
+`j_prod!(xk, s, out)` computes `J(xk) * s`, and `jt_prod!(xk, r, out)` computes `J(xk)' * r`.
+
+`σ > 0` is a regularization parameter and `v` is a vector of the same size as `F(xk)` used for intermediary computations.
+"""
+mutable struct LMModel{
+  T <: Real,
+  V <: AbstractVector{T},
+  Jac <: Union{AbstractMatrix, AbstractLinearOperator},
+} <: AbstractNLPModel{T, V}
+  J::Jac
+  F::V
+  v::V
+  xk::V
+  σ::T
+  meta::NLPModelMeta{T, V}
+  counters::Counters
+end
+
+function LMModel(J::Jac, F::V, σ::T, xk::V) where {T, V, Jac}
+  meta = NLPModelMeta(
+    length(xk),
+    x0 = xk, # Perhaps we should add lvar and uvar as well here.
+  )
+  v = similar(F)
+  return LMModel(J, F, v, xk, σ, meta, Counters())
+end
+
+function NLPModels.obj(nlp::LMModel, x::AbstractVector{T}) where {T}
+  @lencheck nlp.meta.nvar x
+  increment!(nlp, :neval_obj)
+  mul!(nlp.v, nlp.J, x)
+  nlp.v .+= nlp.F
+  return (dot(nlp.v, nlp.v) + nlp.σ * dot(x, x)) / 2
+end
+
+function NLPModels.grad!(nlp::LMModel, x::AbstractVector{T}, g::AbstractVector{T}) where {T}
+  @lencheck nlp.meta.nvar x
+  @lencheck nlp.meta.nvar g
+  increment!(nlp, :neval_grad)
+  mul!(nlp.v, nlp.J, x)
+  nlp.v .+= nlp.F
+  mul!(g, nlp.J', nlp.v)
+  @. g += nlp.σ .* x
+  return g
+end