add an option to store the hessian in sparse format for R2N

Author: MaxenceGollier

Project.toml

@@ -17,12 +17,11 @@ ProximalOperators = "a725b495-10eb-56fe-b38b-717eba820537"
 RegularizedProblems = "ea076b23-609f-44d2-bb12-a4ae45328278"
 ShiftedProximalOperators = "d4fd37fa-580c-4e43-9b30-361c21aae263"
 SolverCore = "ff4d7338-4cf1-434d-91df-b86cb86fb843"
+SparseMatricesCOO = "fa32481b-f100-4b48-8dc8-c62f61b13870"
 
 [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"
@@ -31,6 +30,7 @@ ProximalOperators = "0.15"
 RegularizedProblems = "0.1.1"
 ShiftedProximalOperators = "0.2"
 SolverCore = "0.3.0"
+SparseMatricesCOO = "0.2.4"
 julia = "^1.6.0"
 
 [extras]

src/R2N.jl

@@ -20,6 +20,7 @@ mutable struct R2NSolver{
   s1::V
   v0::V
   has_bnds::Bool
+  store_h::Bool
   l_bound::V
   u_bound::V
   l_bound_m_x::V
@@ -34,6 +35,7 @@ function R2NSolver(
   reg_nlp::AbstractRegularizedNLPModel{T, V};
   subsolver = R2Solver,
   m_monotone::Int = 1,
+  store_h = false
 ) where {T, V}
   x0 = reg_nlp.model.meta.x0
   l_bound = reg_nlp.model.meta.lvar
@@ -67,7 +69,15 @@ function R2NSolver(
     has_bnds ? shifted(reg_nlp.h, xk, l_bound_m_x, u_bound_m_x, reg_nlp.selected) :
     shifted(reg_nlp.h, xk)
 
-  Bk = hess_op(reg_nlp.model, x0)
+  store_h = isa(reg_nlp.model, QuasiNewtonModel) ? false : store_h
+
+  if !store_h
+    Bk = hess_op(reg_nlp.model, x0) 
+  else
+    rows, cols = hess_structure(reg_nlp.model)
+    vals = hess_coord(reg_nlp.model, x0)
+    Bk = SparseMatrixCOO(reg_nlp.model.meta.nvar, reg_nlp.model.meta.nvar, rows, cols, vals)
+  end
   sub_nlp = R2NModel(Bk, ∇fk, T(1), x0)
   subpb = RegularizedNLPModel(sub_nlp, ψ)
   substats = RegularizedExecutionStats(subpb)
@@ -85,6 +95,7 @@ function R2NSolver(
     s1,
     v0,
     has_bnds,
+    store_h,
     l_bound,
     u_bound,
     l_bound_m_x,
@@ -217,8 +228,7 @@ function SolverCore.solve!(
   γ::T = T(3),
   β::T = 1 / eps(T),
   θ::T = 1/(1 + eps(T)^(1 / 5)),
-  opnorm_maxiter::Int = 5,
-  sub_kwargs::NamedTuple = NamedTuple(),
+  sub_kwargs::Dict{Symbol, T} = Dict{Symbol, T}(),
 ) where {T, V, G}
   reset!(stats)
 
@@ -291,8 +301,11 @@ function SolverCore.solve!(
 
   quasiNewtTest = isa(nlp, QuasiNewtonModel)
   λmax::T = T(1)
-  found_λ = true
-  solver.subpb.model.B = hess_op(nlp, xk)
+  if !solver.store_h 
+    solver.subpb.model.B = hess_op(nlp, xk)
+  else
+    hess_coord!(nlp, xk, solver.subpb.model.B.vals)
+  end
 
   if opnorm_maxiter ≤ 0
     λmax, found_λ = opnorm(solver.subpb.model.B)
@@ -445,7 +458,11 @@ function SolverCore.solve!(
         push!(nlp, s, solver.y)
         qn_copy!(nlp, solver, stats)
       end
-      solver.subpb.model.B = hess_op(nlp, xk)
+      if !solver.store_h 
+        solver.subpb.model.B = hess_op(nlp, xk)
+      else
+        hess_coord!(nlp, xk, solver.subpb.model.B.vals)
+      end
 
       if opnorm_maxiter ≤ 0
         λmax, found_λ = opnorm(solver.subpb.model.B)

src/R2NModel.jl

@@ -14,7 +14,7 @@ this model represents the smooth R2N subproblem:
 where `B` is either an approximation of the Hessian of `f` or the Hessian itself and `∇f` represents the gradient of `f` at `x0`.
 `σ > 0` is a regularization parameter and `v` is a vector of the same size as `x0` used for intermediary computations.
 """
-mutable struct R2NModel{T <: Real, V <: AbstractVector{T}, G <: AbstractLinearOperator{T}} <:
+mutable struct R2NModel{T <: Real, V <: AbstractVector{T}, G <: Union{AbstractLinearOperator{T}, AbstractMatrix{T}}} <:
                AbstractNLPModel{T, V}
   B::G
   ∇f::V
@@ -34,15 +34,22 @@ function R2NModel(B::G, ∇f::V, σ::T, x0::V) where {T, V, G}
   return R2NModel(B::G, ∇f::V, v::V, σ::T, meta, Counters())
 end
 
-function NLPModels.obj(nlp::R2NModel, x::AbstractVector; skip_sigma = false)
+function NLPModels.obj(nlp::R2NModel{T, V, M}, x::V) where{T, V, M <: AbstractLinearOperator{T}}
   @lencheck nlp.meta.nvar x
   increment!(nlp, :neval_obj)
   mul!(nlp.v, nlp.B, x)
   skip_sigma && return dot(nlp.v, x)/2 + dot(nlp.∇f, x)
   return dot(nlp.v, x)/2 + dot(nlp.∇f, x) + nlp.σ * dot(x, x) / 2
 end
 
-function NLPModels.grad!(nlp::R2NModel, x::AbstractVector, g::AbstractVector)
+function NLPModels.obj(nlp::R2NModel{T, V, M}, x::V) where{T, V, M <: AbstractMatrix{T}}
+  @lencheck nlp.meta.nvar x
+  increment!(nlp, :neval_obj)
+  mul!(nlp.v, Symmetric(nlp.B,:L), x)
+  return dot(nlp.v, x)/2 + dot(nlp.∇f, x) + nlp.σ * dot(x, x) / 2
+end
+
+function NLPModels.grad!(nlp::R2NModel{T, V, M}, x::V, g::V) where{T, V, M <: AbstractLinearOperator{T}}
   @lencheck nlp.meta.nvar x
   @lencheck nlp.meta.nvar g
   increment!(nlp, :neval_grad)
@@ -51,3 +58,13 @@ function NLPModels.grad!(nlp::R2NModel, x::AbstractVector, g::AbstractVector)
   g .+= nlp.σ .* x
   return g
 end
+
+function NLPModels.grad!(nlp::R2NModel{T, V, M}, x::V, g::V) where{T, V, M <: AbstractMatrix{T}}
+  @lencheck nlp.meta.nvar x
+  @lencheck nlp.meta.nvar g
+  increment!(nlp, :neval_grad)
+  mul!(g, Symmetric(nlp.B,:L), x)
+  g .+= nlp.∇f
+  g .+= nlp.σ .* x
+  return g
+end

src/RegularizedOptimization.jl

@@ -8,12 +8,7 @@ using Arpack, ProximalOperators
 
 # dependencies from us
 using LinearOperators,
-  ManualNLPModels,
-  NLPModels,
-  NLPModelsModifiers,
-  RegularizedProblems,
-  ShiftedProximalOperators,
-  SolverCore
+  NLPModels, NLPModelsModifiers, RegularizedProblems, ShiftedProximalOperators, SolverCore, SparseMatricesCOO
 using Percival: AugLagModel, update_y!, update_μ!
 
 const callback_docstring = "