Use operator norm upper bound in solvers

Project.toml

@@ -26,7 +26,7 @@ SolverCore = "ff4d7338-4cf1-434d-91df-b86cb86fb843"
 [compat]
 ADNLPModels = "0.8.13"
 Arpack = "0.5"
-LinearOperators = "2.10.0"
+LinearOperators = "2.13.0"
 ManualNLPModels = "0.2.0"
 NLPModels = "0.19, 0.20, 0.21"
 NLPModelsModifiers = "0.7"

src/R2N.jl

@@ -147,7 +147,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;
-- `opnorm_maxiter::Int = 5`: how many iterations of the power method to use to compute the operator norm of Bₖ. If a negative number is provided, then Arpack is used instead;
+- `opnorm_maxiter::Int = 5`: how many iterations of the power method to use to compute the operator norm of Bₖ. If a negative number is provided, an upper bound of the operator norm is computed: see `opnorm_upper_bound`.
 - `m_monotone::Int = 1`: monotonicity parameter. By default, R2N is monotone but the non-monotone variant will be used if `m_monotone > 1`;
 - `compute_obj::Bool = true`: (advanced) whether `f(x₀)` should be computed or not. If set to false, then the value is retrieved from `stats.solver_specific[:smooth_obj]`;
 - `compute_grad::Bool = true`: (advanced) whether `∇f(x₀)` should be computed or not. If set to false, then the value is retrieved from `solver.∇fk`;
@@ -306,9 +306,9 @@ function SolverCore.solve!(
   found_λ = true
 
   if opnorm_maxiter ≤ 0
-    λmax, found_λ = opnorm(solver.subpb.model.B)
+    λmax, found_λ = opnorm_upper_bound(solver.subpb.model.B)
   else
-    λmax = power_method!(solver.subpb.model.B, solver.v0, solver.subpb.model.v, opnorm_maxiter)
+    λmax, found_λ = power_method!(solver.subpb.model.B, solver.v0, solver.subpb.model.v, opnorm_maxiter)
   end
   found_λ || error("operator norm computation failed")
 
@@ -458,9 +458,9 @@ function SolverCore.solve!(
       end
 
       if opnorm_maxiter ≤ 0
-        λmax, found_λ = opnorm(solver.subpb.model.B)
+        λmax, found_λ = opnorm_upper_bound(solver.subpb.model.B)
       else
-        λmax = power_method!(solver.subpb.model.B, solver.v0, solver.subpb.model.v, opnorm_maxiter)
+        λmax, found_λ = power_method!(solver.subpb.model.B, solver.v0, solver.subpb.model.v, opnorm_maxiter)
       end
       found_λ || error("operator norm computation failed")
       set_step_status!(stats, :accepted)

src/TR_alg.jl

@@ -307,9 +307,9 @@ function SolverCore.solve!(
   found_λ = true
 
   if opnorm_maxiter ≤ 0
-    λmax, found_λ = opnorm(solver.subpb.model.B)
+    λmax, found_λ = opnorm_upper_bound(solver.subpb.model.B)
   else
-    λmax = power_method!(solver.subpb.model.B, solver.v0, solver.subpb.model.v, opnorm_maxiter)
+    λmax, found_λ = power_method!(solver.subpb.model.B, solver.v0, solver.subpb.model.v, opnorm_maxiter)
   end
   found_λ || error("operator norm computation failed")
 
@@ -468,9 +468,9 @@ function SolverCore.solve!(
       end
 
       if opnorm_maxiter ≤ 0
-        λmax, found_λ = opnorm(solver.subpb.model.B)
+        λmax, found_λ = opnorm_upper_bound(solver.subpb.model.B)
       else
-        λmax = power_method!(solver.subpb.model.B, solver.v0, solver.subpb.model.v, opnorm_maxiter)
+        λmax, found_λ = power_method!(solver.subpb.model.B, solver.v0, solver.subpb.model.v, opnorm_maxiter)
       end
       found_λ || error("operator norm computation failed")
 

src/utils.jl

@@ -17,7 +17,42 @@ function power_method!(B::M, v₀::S, v₁::S, max_iter::Int = 1) where {M, S}
     normalize!(v₁)
   end
   mul!(v₁, B, v₀)
-  return abs(dot(v₀, v₁))
+  approx = abs(dot(v₀, v₁))
+  return approx, !isnan(approx)
+end
+
+# Compute upper bounds for ‖B‖₂.
+
+# For matrices, we compute the Frobenius norm.
+function opnorm_upper_bound(B::AbstractMatrix)
+  bound = norm(B, 2)
+  return bound, !isnan(bound)
+end
+
+# For LBFGS, using the formula Bₖ = B\_{k-1} - aₖaₖᵀ + bₖbₖᵀ, we compute
+# ‖Bₖ‖₂ ≤ ‖B₀‖₂ + ∑ᵢ ‖bᵢ‖₂².
+function opnorm_upper_bound(B::LBFGSOperator{T}) where {T}
+  upper_bound = B.data.opnorm_upper_bound
+  return upper_bound, !isnan(upper_bound)
+end
+
+# For LSR1, we use the formula Bₖ = B\_{k-1} + σₖaₖaₖᵀ, we compute
+# ‖Bₖ‖₂ ≤ ‖B₀‖₂ + ∑ᵢ |σᵢ|‖aᵢ‖₂².
+function opnorm_upper_bound(B::LSR1Operator{T}) where {T}
+  upper_bound = B.data.opnorm_upper_bound
+  return upper_bound, !isnan(upper_bound)
+end
+
+# For diagonal operators, we compute the exact operator norm.
+function opnorm_upper_bound(B::AbstractDiagonalQuasiNewtonOperator)
+  bound = norm(B.d, Inf)
+  return bound, !isnan(bound)
+end
+
+# In the general case, we use Arpack.
+# Note: Arpack allocates.
+function opnorm_upper_bound(B::AbstractLinearOperator)
+  return opnorm(B)
 end
 
 # use Arpack to obtain largest eigenvalue in magnitude with a minimum of robustness

test/runtests.jl

@@ -2,6 +2,7 @@ using LinearAlgebra: length
 using LinearAlgebra, Random, Test
 using ProximalOperators
 using ADNLPModels,
+  LinearOperators,
   OptimizationProblems,
   OptimizationProblems.ADNLPProblems,
   NLPModels,
@@ -19,6 +20,7 @@ const global bpdn2, bpdn_nls2, sol2 = bpdn_model(compound, bounds = true)
 const global λ = norm(grad(bpdn, zeros(bpdn.meta.nvar)), Inf) / 10
 
 include("test_AL.jl")
+include("test-utils.jl")
 
 for (mod, mod_name) ∈ ((x -> x, "exact"), (LSR1Model, "lsr1"), (LBFGSModel, "lbfgs"))
   for (h, h_name) ∈ ((NormL0(λ), "l0"), (NormL1(λ), "l1"), (IndBallL0(10 * compound), "B0"))

test/test-utils.jl

@@ -0,0 +1,100 @@
+function test_opnorm_upper_bound(B, m, n)
+  T = eltype(B)
+  is_upper_bound = true
+  for _ = 1:m
+    upper_bound, found =  RegularizedOptimization.opnorm_upper_bound(B)
+    if opnorm(Matrix(B)) > upper_bound || !found
+      is_upper_bound = false
+      break
+    end
+
+    push!(B, randn(T, n), randn(T, n))
+  end
+
+  nallocs = @allocated RegularizedOptimization.opnorm_upper_bound(B)
+  @test nallocs == 0
+  @test is_upper_bound == true
+end
+
+# Test norm functions
+
+@testset "Test opnorm upper bound functions" begin
+  n = 10
+  m = 40
+  @testset "LBFGS" begin
+    B = LBFGSOperator(Float64, n, scaling = false)
+    test_opnorm_upper_bound(B, m, n)
+
+    B = LBFGSOperator(Float64, n, scaling = true)
+    test_opnorm_upper_bound(B, m, n)
+
+    B = LBFGSOperator(Float32, n, scaling = false)
+    test_opnorm_upper_bound(B, m, n)
+
+    B = LBFGSOperator(Float32, n, scaling = true)
+    test_opnorm_upper_bound(B, m, n)
+  end
+
+  @testset "LSR1" begin
+    B = LSR1Operator(Float64, n, scaling = false)
+    test_opnorm_upper_bound(B, m, n)
+
+    B = LSR1Operator(Float64, n, scaling = true)
+    test_opnorm_upper_bound(B, m, n)
+
+    B = LSR1Operator(Float32, n, scaling = false)
+    test_opnorm_upper_bound(B, m, n)
+
+    B = LSR1Operator(Float32, n, scaling = true)
+    test_opnorm_upper_bound(B, m, n)
+  end
+
+  @testset "Diagonal" begin
+    B = SpectralGradient(randn(Float64), n)
+    test_opnorm_upper_bound(B, m, n)
+
+    B = SpectralGradient(randn(Float32), n)
+    test_opnorm_upper_bound(B, m, n)
+
+    B = DiagonalPSB(randn(Float64, n))
+    test_opnorm_upper_bound(B, m, n)
+
+    B = DiagonalPSB(randn(Float32, n))
+    test_opnorm_upper_bound(B, m, n)
+
+    B = DiagonalAndrei(randn(Float64, n))
+    test_opnorm_upper_bound(B, m, n)
+
+    B = DiagonalAndrei(randn(Float32, n))
+    test_opnorm_upper_bound(B, m, n)
+  end
+
+  @testset "Matrix" begin
+    is_upper_bound = true
+    for _ = 1:m
+      B = randn(n, n)
+      upper_bound, found =  RegularizedOptimization.opnorm_upper_bound(B)
+      if opnorm(Matrix(B)) > upper_bound || !found
+        is_upper_bound = false
+        break
+      end
+    end
+
+    @test is_upper_bound == true
+  end
+
+  @testset "Arpack" begin
+    is_upper_bound = true
+    for _ = 1:m
+      B = randn(n, n)
+      B = LinearOperator((B + B')/2, symmetric = true)
+      upper_bound, found =  RegularizedOptimization.opnorm_upper_bound(B)
+
+      if opnorm(Matrix(B)) > upper_bound + 1e-12 || !found
+        is_upper_bound = false
+        break
+      end
+    end
+    @test is_upper_bound == true
+  end
+end