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

Author: MaxenceGollier

Project.toml

@@ -4,6 +4,7 @@ author = ["Robert Baraldi <[email protected]> and Dominique Orban <dominique.orban
 version = "0.1.0"
 
 [deps]
+Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearOperators = "5c8ed15e-5a4c-59e4-a42b-c7e8811fb125"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
@@ -15,18 +16,19 @@ ProximalOperators = "a725b495-10eb-56fe-b38b-717eba820537"
 RegularizedProblems = "ea076b23-609f-44d2-bb12-a4ae45328278"
 ShiftedProximalOperators = "d4fd37fa-580c-4e43-9b30-361c21aae263"
 SolverCore = "ff4d7338-4cf1-434d-91df-b86cb86fb843"
-Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
+SparseMatricesCOO = "fa32481b-f100-4b48-8dc8-c62f61b13870"
 
 [compat]
+Arpack = "0.5"
 LinearOperators = "2.10.0"
-NLPModels = "0.19, 0.20"
+NLPModels = "0.19, 0.20, 0.21"
 NLPModelsModifiers = "0.7"
 Percival = "0.7.2"
 ProximalOperators = "0.15"
 RegularizedProblems = "0.1.1"
 ShiftedProximalOperators = "0.2"
 SolverCore = "0.3.0"
-Arpack = "0.5"
+SparseMatricesCOO = "0.2.4"
 julia = "^1.6.0"
 
 [extras]

src/R2N.jl

@@ -18,6 +18,7 @@ mutable struct R2NSolver{
   s::V
   s1::V
   has_bnds::Bool
+  store_h::Bool
   l_bound::V
   u_bound::V
   l_bound_m_x::V
@@ -32,6 +33,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
@@ -60,7 +62,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)
@@ -76,6 +86,7 @@ function R2NSolver(
     s,
     s1,
     has_bnds,
+    store_h,
     l_bound,
     u_bound,
     l_bound_m_x,
@@ -215,7 +226,7 @@ function SolverCore.solve!(
   γ::T = T(3),
   β::T = 1 / eps(T),
   θ::T = 1/(1 + eps(T)^(1 / 5)),
-  sub_kwargs::Dict{Symbol} = Dict(),
+  sub_kwargs::Dict{Symbol, T} = Dict{Symbol, T}(),
 ) where {T, V, G}
   reset!(stats)
 
@@ -287,7 +298,11 @@ function SolverCore.solve!(
 
   quasiNewtTest = isa(nlp, QuasiNewtonModel)
   λmax::T = T(1)
-  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
 
   λmax, found_λ = opnorm(solver.subpb.model.B)
   found_λ || error("operator norm computation failed")
@@ -417,7 +432,11 @@ function SolverCore.solve!(
         @. ∇fk⁻ = ∇fk - ∇fk⁻
         push!(nlp, s, ∇fk⁻)
       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
 
       λmax, found_λ = opnorm(solver.subpb.model.B)
       found_λ || error("operator norm computation failed")

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,14 +34,21 @@ 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)
+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)
   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)
@@ -50,3 +57,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,7 +8,7 @@ using Arpack, ProximalOperators
 
 # dependencies from us
 using LinearOperators,
-  NLPModels, NLPModelsModifiers, RegularizedProblems, ShiftedProximalOperators, SolverCore
+  NLPModels, NLPModelsModifiers, RegularizedProblems, ShiftedProximalOperators, SolverCore, SparseMatricesCOO
 using Percival: AugLagModel, update_y!, update_μ!
 
 include("utils.jl")