add rayliegh quotient

Author: MaxenceGollier

src/R2N.jl

@@ -17,6 +17,7 @@ mutable struct R2NSolver{
   xkn::V
   s::V
   s1::V
+  v0::V
   has_bnds::Bool
   l_bound::V
   u_bound::V
@@ -44,6 +45,10 @@ function R2NSolver(
   xkn = similar(x0)
   s = similar(x0)
   s1 = similar(x0)
+
+  v0 = randn(T, length(x0))
+  v0 /= norm(v0)
+
   has_bnds = any(l_bound .!= T(-Inf)) || any(u_bound .!= T(Inf))
   if has_bnds
     l_bound_m_x = similar(xk)
@@ -75,6 +80,7 @@ function R2NSolver(
     xkn,
     s,
     s1,
+    v0,
     has_bnds,
     l_bound,
     u_bound,
@@ -130,6 +136,7 @@ For advanced usage, first define a solver "R2NSolver" to preallocate the memory
 - `η2::T = T(0.9)`: very successful iteration threshold;
 - `γ::T = T(3)`: regularization parameter multiplier, σ := σ/γ when the iteration is very successful and σ := σγ when the iteration is unsuccessful;
 - `θ::T = 1/(1 + eps(T)^(1 / 5))`: is the model decrease fraction with respect to the decrease of the Cauchy model;
+- `compute_opnorm::Bool = false`: whether the operator norm of Bₖ should be computed at each iteration. If false, a Rayleigh quotient is computed instead. The first option causes the solver to converge in fewer iterations but the computational cost per iteration is larger;
 - `m_monotone::Int = 1`: monotonicity parameter. By default, R2N is monotone but the non-monotone variant will be used if `m_monotone > 1`;
 - `sub_kwargs::Dict{Symbol}`: a dictionary containing the keyword arguments to be sent to the subsolver. The solver will fail if invalid keyword arguments are provided to the subsolver.
 
@@ -215,6 +222,7 @@ function SolverCore.solve!(
   γ::T = T(3),
   β::T = 1 / eps(T),
   θ::T = 1/(1 + eps(T)^(1 / 5)),
+  compute_opnorm::Bool = false,
   sub_kwargs::Dict{Symbol} = Dict(),
 ) where {T, V, G}
   reset!(stats)
@@ -287,9 +295,15 @@ function SolverCore.solve!(
 
   quasiNewtTest = isa(nlp, QuasiNewtonModel)
   λmax::T = T(1)
+  found_λ = true
   solver.subpb.model.B = hess_op(nlp, xk)
 
-  λmax, found_λ = opnorm(solver.subpb.model.B)
+  if !compute_opnorm
+    mul!(solver.subpb.model.v, solver.subpb.model.B, solver.v0)
+    λmax = dot(solver.v0, solver.subpb.model.v)
+  else
+    λmax, found_λ = opnorm(solver.subpb.model.B)
+  end
   found_λ || error("operator norm computation failed")
 
   ν₁ = θ / (λmax + σk)
@@ -419,7 +433,13 @@ function SolverCore.solve!(
       end
       solver.subpb.model.B = hess_op(nlp, xk)
 
-      λmax, found_λ = opnorm(solver.subpb.model.B)
+      if !compute_opnorm
+        mul!(solver.subpb.model.v, solver.subpb.model.B, solver.v0)
+        λmax = dot(solver.v0, solver.subpb.model.v)
+      else
+        λmax, found_λ = opnorm(solver.subpb.model.B)
+      end
+      
       found_λ || error("operator norm computation failed")
 
       ∇fk⁻ .= ∇fk

test/test_allocs.jl

@@ -44,14 +44,13 @@ end
   for (h, h_name) ∈ ((NormL0(λ), "l0"),)
     for (solver, solver_name) ∈ ((:R2Solver, "R2"), (:R2DHSolver, "R2DH"), (:R2NSolver, "R2N"))
       @testset "$(solver_name)" begin
-        solver_name == "R2N" && continue #FIXME
         reg_nlp = RegularizedNLPModel(LBFGSModel(bpdn), h)
         solver = eval(solver)(reg_nlp)
         stats = RegularizedExecutionStats(reg_nlp)
         solver_name == "R2" &&
           @test @wrappedallocs(solve!(solver, reg_nlp, stats, ν = 1.0, atol = 1e-6, rtol = 1e-6)) ==
                 0
-        solver_name == "R2DH" && @test @wrappedallocs(
+        (solver_name == "R2DH" || solver_name == "R2N") && @test @wrappedallocs(
           solve!(solver, reg_nlp, stats, σk = 1.0, atol = 1e-6, rtol = 1e-6)
         ) == 0
         @test stats.status == :first_order