Add CaNNOLeS as an alternative feasibility step solver

Project.toml

@@ -14,6 +14,12 @@ SolverCore = "ff4d7338-4cf1-434d-91df-b86cb86fb843"
 SolverTools = "b5612192-2639-5dc1-abfe-fbedd65fab29"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
+[weakdeps]
+CaNNOLeS = "5a1c9e79-9c58-5ec0-afc4-3298fdea2875"
+
+[extensions]
+DCISolverCaNNOLeSExt = ["CaNNOLeS"]
+
 [compat]
 HSL = "0.5"
 Krylov = "0.10.1"

README.md

@@ -21,6 +21,14 @@ In particular, [NLPModels.jl](https://github.com/JuliaSmoothOptimizers/NLPModels
 It uses [LDLFactorizations.jl](https://github.com/JuliaSmoothOptimizers/LDLFactorizations.jl) by default to compute the factorization in the tangent step. [HSL.jl](https://github.com/JuliaSmoothOptimizers/HSL.jl) provides alternative linear solvers if [libHSL](https://licences.stfc.ac.uk/product/libhsl) can be downloaded.
 The feasibility steps are factorization-free and use iterative methods from [Krylov.jl](https://github.com/JuliaSmoothOptimizers/Krylov.jl)
 
+## CaNNOLeS extension
+
+CaNNOLeS is optional and loaded through package extensions. Install it in your environment to enable the CaNNOLeS-based feasibility step:
+
+- `using Pkg; Pkg.add("CaNNOLeS")`
+- call `dci(...; feas_step = :feasibility_step_cannoles, cannoles_options = Dict(...))`
+- extra keywords in `cannoles_options` are forwarded directly to `CaNNOLeS.cannoles`; without CaNNOLeS the default trust-region feasibility step is used.
+
 ## References
 
 > Bielschowsky, R. H., & Gomes, F. A.

ext/DCISolverCaNNOLeSExt.jl

@@ -0,0 +1,88 @@
+module DCISolverCaNNOLeSExt
+
+import DCISolver
+using DCISolver: MetaDCI, cons_norhs!, FeasibilityResidual
+using CaNNOLeS: cannoles
+using LinearAlgebra: norm
+using NLPModels: AbstractNLPModel, jac_op!, neval_obj, neval_cons
+using SolverCore: log_row
+
+function DCISolver.feasibility_step_cannoles(
+  nlp::AbstractNLPModel,
+  x::AbstractVector{T},
+  cx::AbstractVector{T},
+  normcx::T,
+  Jx,
+  ρ::T,
+  ctol::AbstractFloat,
+  meta::MetaDCI,
+  workspace,
+  verbose::Bool;
+  max_eval::Int = 1_000,
+  max_time::AbstractFloat = 60.0,
+  max_iter::Int = typemax(Int64),
+  cannoles_options = Dict{Symbol, Any}(),
+) where {T}
+  nls = FeasibilityResidual(nlp, x)
+
+  default_options = Dict{Symbol, Any}(
+    :atol => ctol,
+    :rtol => ctol,
+    :Fatol => ctol,
+    :Frtol => ctol,
+    :max_eval => max_eval,
+    :max_time => max_time,
+    :max_iter => max_iter,
+    :verbose => verbose ? 1 : 0,
+  )
+
+  options = merge(default_options, cannoles_options)
+
+  start_time = time()
+  stats = cannoles(nls; options...)
+  el_time = time() - start_time
+
+  z = stats.solution
+  DCISolver.cons_norhs!(nlp, z, workspace.cz)
+  cz = workspace.cz
+  normcz = norm(cz)
+  Jz = jac_op!(nlp, z, workspace.Jv, workspace.Jtv)
+
+  status = if stats.status == :first_order && normcz ≤ ρ
+    :success
+  elseif stats.status == :first_order || stats.status == :acceptable
+    normcz ≤ ρ ? :success : :unknown
+  elseif stats.status == :max_eval
+    :max_eval
+  elseif stats.status == :max_time
+    :max_time
+  elseif stats.status == :max_iter
+    :max_iter
+  elseif stats.status == :infeasible
+    :infeasible
+  else
+    :unknown
+  end
+
+  verbose && @info log_row(
+    Any[
+      "F-CaNNOLeS",
+      stats.iter,
+      neval_obj(nlp) + neval_cons(nlp),
+      Float64,
+      Float64,
+      Float64,
+      normcz,
+      Float64,
+      ρ,
+      status,
+      Float64,
+      Float64,
+      el_time,
+    ],
+  )
+
+  return z, cz, normcz, Jz, status
+end
+
+end # module

src/DCISolver.jl

@@ -17,9 +17,11 @@ using LinearAlgebra: LinearAlgebra, I, Symmetric, convert, ldiv!, mul!, norm, tr
 using NLPModels:
   NLPModels,
   AbstractNLPModel,
+  AbstractNLSModel,
   cons!,
   equality_constrained,
   get_lcon,
+  get_ucon,
   hess_coord!,
   hess_op,
   hess_structure!,
@@ -29,7 +31,10 @@ using NLPModels:
   jac_structure!,
   neval_cons,
   neval_obj,
-  unconstrained
+  unconstrained,
+  NLPModelMeta,
+  NLSMeta,
+  NLSCounters
 using SolverCore:
   SolverCore,
   AbstractOptimizationSolver,

src/dci_feasibility.jl

@@ -315,3 +315,173 @@ function TR_lsmr(
 
   return d, Jd, infeasible, solved
 end
+
+"""
+	feasibility_step_cannoles(nlp, x, cx, normcx, Jx, ρ, ctol, meta, workspace, verbose; kwargs...)
+
+Provided by the optional CaNNOLeS extension. Add `CaNNOLeS.jl` (and `NLPModelsModifiers.jl`) to enable.
+"""
+function feasibility_step_cannoles end
+
+"""
+	FeasibilityResidual
+
+A wrapper to convert an NLP with equality constraints into a nonlinear least squares
+problem by treating the constraints as residuals. This allows us to use CaNNOLeS
+to solve the feasibility problem min ||c(x)||^2 / 2.
+"""
+mutable struct FeasibilityResidual{T, S, M <: AbstractNLPModel{T, S}} <: AbstractNLSModel{T, S}
+  nlp::M
+  meta::NLPModelMeta{T, S}
+  nls_meta::NLSMeta{T, S}
+  counters::NLSCounters
+end
+
+function FeasibilityResidual(nlp::AbstractNLPModel{T, S}, x0::S) where {T, S}
+  ncon = nlp.meta.ncon
+  nvar = nlp.meta.nvar
+
+  meta = NLPModelMeta(
+    nvar,
+    x0 = x0,
+    ncon = ncon,
+    lcon = get_lcon(nlp),
+    ucon = get_ucon(nlp),
+    nnzj = nlp.meta.nnzj,
+    nnzh = nlp.meta.nnzh,
+    minimize = true,
+  )
+
+  nls_meta = NLSMeta{T, S}(
+    ncon,  # nequ - number of residuals equals number of constraints
+    nvar,
+    nnzj = nlp.meta.nnzj,
+    nnzh = nlp.meta.nnzh,
+  )
+
+  return FeasibilityResidual(nlp, meta, nls_meta, NLSCounters())
+end
+
+function NLPModels.residual!(nls::FeasibilityResidual, x::AbstractVector, rx::AbstractVector)
+  cons!(nls.nlp, x, rx)
+  if nls.nlp.meta.ncon > 0
+    rx .-= get_lcon(nls.nlp)
+  end
+  nls.counters.neval_residual += 1
+  return rx
+end
+
+function NLPModels.jac_structure_residual!(
+  nls::FeasibilityResidual,
+  rows::AbstractVector{<:Integer},
+  cols::AbstractVector{<:Integer},
+)
+  return jac_structure!(nls.nlp, rows, cols)
+end
+
+function NLPModels.jac_nln_structure!(
+  nls::FeasibilityResidual,
+  rows::AbstractVector{<:Integer},
+  cols::AbstractVector{<:Integer},
+)
+  return jac_structure!(nls.nlp, rows, cols)
+end
+
+function NLPModels.jac_coord_residual!(
+  nls::FeasibilityResidual,
+  x::AbstractVector,
+  vals::AbstractVector,
+)
+  jac_coord!(nls.nlp, x, vals)
+  nls.counters.neval_jac_residual += 1
+  return vals
+end
+
+function NLPModels.jprod_residual!(
+  nls::FeasibilityResidual,
+  x::AbstractVector,
+  v::AbstractVector,
+  Jv::AbstractVector,
+)
+  jprod!(nls.nlp, x, v, Jv)
+  nls.counters.neval_jprod_residual += 1
+  return Jv
+end
+
+function NLPModels.jtprod_residual!(
+  nls::FeasibilityResidual,
+  x::AbstractVector,
+  v::AbstractVector,
+  Jtv::AbstractVector,
+)
+  jtprod!(nls.nlp, x, v, Jtv)
+  nls.counters.neval_jtprod_residual += 1
+  return Jtv
+end
+
+function NLPModels.hess_structure!(
+  nls::FeasibilityResidual,
+  rows::AbstractVector{<:Integer},
+  cols::AbstractVector{<:Integer},
+)
+  return hess_structure!(nls.nlp, rows, cols)
+end
+
+function NLPModels.hess_structure_residual!(
+  nls::FeasibilityResidual,
+  rows::AbstractVector{<:Integer},
+  cols::AbstractVector{<:Integer},
+)
+  return hess_structure!(nls.nlp, rows, cols)
+end
+
+function NLPModels.hess_coord_residual!(
+  nls::FeasibilityResidual,
+  x::AbstractVector,
+  v::AbstractVector,
+  vals::AbstractVector,
+)
+  # For feasibility problems, hessian of residuals is zero since residuals are linear
+  fill!(vals, zero(eltype(vals)))
+  nls.counters.neval_hess_residual += 1
+  return vals
+end
+
+function NLPModels.hess_coord!(
+  nls::FeasibilityResidual,
+  x::AbstractVector,
+  y::AbstractVector,
+  vals::AbstractVector;
+  obj_weight::Real = one(eltype(x)),
+)
+  hess_coord!(nls.nlp, x, y, vals, obj_weight = obj_weight)
+  nls.counters.neval_hess += 1
+  return vals
+end
+
+function NLPModels.hprod!(
+  nls::FeasibilityResidual,
+  x::AbstractVector,
+  y::AbstractVector,
+  v::AbstractVector,
+  Hv::AbstractVector;
+  obj_weight::Real = one(eltype(x)),
+)
+  hprod!(nls.nlp, x, y, v, Hv, obj_weight = obj_weight)
+  nls.counters.neval_hprod += 1
+  return Hv
+end
+
+# Define cons_nln! for FeasibilityResidual since CaNNOLeS needs it
+function NLPModels.cons_nln!(nls::FeasibilityResidual, x::AbstractVector, cx::AbstractVector)
+  cons!(nls.nlp, x, cx)
+  if nls.nlp.meta.ncon > 0
+    cx .-= get_lcon(nls.nlp)
+  end
+  return cx
+end
+
+function NLPModels.jac_nln_coord!(nls::FeasibilityResidual, x::AbstractVector, vals::AbstractVector)
+  jac_coord!(nls.nlp, x, vals)
+  return vals
+end

src/param_struct.jl

@@ -136,7 +136,7 @@ The keyword arguments may include:
 - `decrease_γ::T=T(0.1)`: Regularization for the factorization: reduce `γ` if possible, `> √eps(T)`, between tangent steps.
 - `increase_γ::T=T(100.0)`: Regularization for the factorization: up `γ` if possible, `< 1/√eps(T)`, during the factorization.
 - `δmin::T=√eps(T)`: Regularization for the factorization: smallest value of `δ` used for the regularization.
-- `feas_step::Symbol=:feasibility_step`: Normal step.
+- `feas_step::Symbol=:feasibility_step`: Normal step. Options: `:feasibility_step` (default trust-region method) or `:feasibility_step_cannoles` (CaNNOLeS solver).
 - `feas_η₁::T=T(1e-3)`: Feasibility step: decrease the trust-region radius when `Ared/Pred < η₁`.
 - `feas_η₂::T=T(0.66)`: Feasibility step: increase the trust-region radius when `Ared/Pred > η₂`.
 - `feas_σ₁::T=T(0.25)`: Feasibility step: decrease coefficient of the trust-region radius.
@@ -195,7 +195,7 @@ struct MetaDCI{
   δmin::T
 
   # Normal step
-  feas_step::Symbol #:feasibility_step
+  feas_step::Symbol #:feasibility_step or :feasibility_step_cannoles
   ## Feasibility step (called inside the normal step)
   feas_η₁::T
   feas_η₂::T

test/Project.toml

@@ -1,16 +1,18 @@
 [deps]
 ADNLPModels = "54578032-b7ea-4c30-94aa-7cbd1cce6c9a"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+CaNNOLeS = "5a1c9e79-9c58-5ec0-afc4-3298fdea2875"
+DCISolver = "bee2e536-65f6-11e9-3844-e5bb4c9c55c9"
 HSL = "34c5aeac-e683-54a6-a0e9-6e0fdc586c50"
 Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6"
 NLPModelsTest = "7998695d-6960-4d3a-85c4-e1bceb8cd856"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SolverCore = "ff4d7338-4cf1-434d-91df-b86cb86fb843"
 SolverTest = "4343dc35-3317-4c6e-8877-f0cc8502c90e"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
@@ -23,7 +25,7 @@ Logging = "1.10"
 NLPModels = "0.21"
 NLPModelsTest = "0.10"
 Random = "1.10"
-SparseArrays = "1.10"
 SolverCore = "0.3"
 SolverTest = "0.3"
+SparseArrays = "1.10"
 Test = "1.10"

test/runtests.jl

@@ -17,9 +17,20 @@ for (root, dirs, files) in walkdir(@__DIR__)
     if isnothing(match(r"^test-.*\.jl$", file))
       continue
     end
+    if file == "test-cannoles-feasibility.jl"
+      continue
+    end
     title = titlecase(replace(splitext(file[6:end])[1], "-" => " "))
     @testset "$title" begin
       include(file)
     end
   end
 end
+
+if Base.find_package("CaNNOLeS") !== nothing
+  @testset "Cannoles Feasibility" begin
+    include("test-cannoles-feasibility.jl")
+  end
+else
+  @info "Skipping CaNNOLeS feasibility tests: CaNNOLeS not installed"
+end