diff --git a/Project.toml b/Project.toml index 5a42353b..06c4ddb9 100644 --- a/Project.toml +++ b/Project.toml @@ -3,6 +3,7 @@ uuid = "10dff2fc-5484-5881-a0e0-c90441020f8a" version = "0.14.1" [deps] +HSL = "34c5aeac-e683-54a6-a0e9-6e0fdc586c50" Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearOperators = "5c8ed15e-5a4c-59e4-a42b-c7e8811fb125" @@ -10,18 +11,27 @@ Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6" NLPModelsModifiers = "e01155f1-5c6f-4375-a9d8-616dd036575f" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +QRMumps = "422b30a1-cc69-4d85-abe7-cc07b540c444" SolverCore = "ff4d7338-4cf1-434d-91df-b86cb86fb843" SolverParameters = "d076d87d-d1f9-4ea3-a44b-64b4cdd1e470" SolverTools = "b5612192-2639-5dc1-abfe-fbedd65fab29" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +SparseMatricesCOO = "fa32481b-f100-4b48-8dc8-c62f61b13870" +TSVD = "9449cd9e-2762-5aa3-a617-5413e99d722e" [compat] -Krylov = "0.10.0" -LinearOperators = "2.0" +HSL = "0.5.1" +Krylov = "0.10.1" +LinearOperators = "2.11.0" NLPModels = "0.21" NLPModelsModifiers = "0.7" +QRMumps = "0.3.1" SolverCore = "0.3" SolverParameters = "0.1" SolverTools = "0.9" +SparseArrays = "1.11.0" +SparseMatricesCOO = "0.2.5" +TSVD = "0.4.4" julia = "1.10" [extras] diff --git a/README.md b/README.md index cd21ff94..3969c9de 100644 --- a/README.md +++ b/README.md @@ -26,7 +26,8 @@ This package provides an implementation of four classic algorithms for unconstra > large scale optimization. *Mathematical Programming*, 45(1), 503-528. > DOI: [10.1007/BF01589116](https://doi.org/10.1007/BF01589116) - +- `R2N`: An inexact second-order quadratic regularization method for unconstrained optimization (with shifted L-BFGS or shifted Hessian operator); + - `R2`: a first-order quadratic regularization method for unconstrained optimization; > E. G. Birgin, J. L. Gardenghi, J. M. Martínez, S. A. Santos, Ph. L. Toint. (2017). @@ -34,6 +35,8 @@ This package provides an implementation of four classic algorithms for unconstra > high-order regularized models. *Mathematical Programming*, 163(1), 359-368. > DOI: [10.1007/s10107-016-1065-8](https://doi.org/10.1007/s10107-016-1065-8) +- `R2NLS`: an inexact second-order quadratic regularization method for nonlinear least-squares problems; + - `fomo`: a first-order method with momentum for unconstrained optimization; - `tron`: a pure Julia implementation of TRON, a trust-region solver for bound-constrained optimization described in @@ -63,7 +66,7 @@ using JSOSolvers, ADNLPModels # Rosenbrock nlp = ADNLPModel(x -> 100 * (x[2] - x[1]^2)^2 + (x[1] - 1)^2, [-1.2; 1.0]) -stats = lbfgs(nlp) # or trunk, tron, R2 +stats = lbfgs(nlp) # or trunk, tron, R2, R2N ``` ## How to cite diff --git a/docs/src/solvers.md b/docs/src/solvers.md index 322f7c2e..44e33bcd 100644 --- a/docs/src/solvers.md +++ b/docs/src/solvers.md @@ -7,11 +7,12 @@ - [`trunk`](@ref) - [`R2`](@ref) - [`fomo`](@ref) +- [`R2N`](@ref) | Problem type | Solvers | | --------------------- | -------- | -| Unconstrained NLP | [`lbfgs`](@ref), [`tron`](@ref), [`trunk`](@ref), [`R2`](@ref), [`fomo`](@ref)| -| Unconstrained NLS | [`trunk`](@ref), [`tron`](@ref) | +| Unconstrained NLP | [`lbfgs`](@ref), [`tron`](@ref), [`trunk`](@ref), [`R2`](@ref), [`fomo`](@ref), ['R2N'](@ref)| +| Unconstrained NLS | [`trunk`](@ref), [`tron`](@ref), [`R2NLS`](@ref)| | Bound-constrained NLP | [`tron`](@ref) | | Bound-constrained NLS | [`tron`](@ref) | @@ -22,5 +23,7 @@ lbfgs tron trunk R2 +R2NLS fomo +R2N ``` diff --git a/myR2N/Project.toml b/myR2N/Project.toml new file mode 100644 index 00000000..458aa458 --- /dev/null +++ b/myR2N/Project.toml @@ -0,0 +1,2 @@ +[deps] +HSL_jll = "017b0a0e-03f4-516a-9b91-836bbd1904dd" diff --git a/my_R2N_dev/Project.toml b/my_R2N_dev/Project.toml new file mode 100644 index 00000000..dd952fd3 --- /dev/null +++ b/my_R2N_dev/Project.toml @@ -0,0 +1,13 @@ +[deps] +ADNLPModels = "54578032-b7ea-4c30-94aa-7cbd1cce6c9a" +HSL = "34c5aeac-e683-54a6-a0e9-6e0fdc586c50" +HSL_jll = "017b0a0e-03f4-516a-9b91-836bbd1904dd" +JSOSolvers = "10dff2fc-5484-5881-a0e0-c90441020f8a" +Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" +LinearOperators = "5c8ed15e-5a4c-59e4-a42b-c7e8811fb125" +NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6" +NLPModelsModifiers = "e01155f1-5c6f-4375-a9d8-616dd036575f" +Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" +SolverCore = "ff4d7338-4cf1-434d-91df-b86cb86fb843" +SolverTools = "b5612192-2639-5dc1-abfe-fbedd65fab29" +SparseMatricesCOO = "fa32481b-f100-4b48-8dc8-c62f61b13870" diff --git a/my_R2N_dev/ideas_for_test.jl b/my_R2N_dev/ideas_for_test.jl new file mode 100644 index 00000000..77768c94 --- /dev/null +++ b/my_R2N_dev/ideas_for_test.jl @@ -0,0 +1,64 @@ +# 1. Test error for constrained problems +@testset "Constrained Problems Error" begin + f(x) = (x[1] - 1.0)^2 + 100 * (x[2] - x[1]^2)^2 + x0 = [-1.2; 1.0] + lvar = [-Inf; 0.1] + uvar = [0.5; 0.5] + c(x) = [x[1] + x[2] - 2; x[1]^2 + x[2]^2] + lcon = [0.0; -Inf] + ucon = [Inf; 1.0] + nlp_constrained = ADNLPModel(f, x0, lvar, uvar, c, lcon, ucon) + + solver = R2NSolver(nlp_constrained) + stats = GenericExecutionStats(nlp_constrained) + + @test_throws ErrorException begin + SolverCore.solve!(solver, nlp_constrained, stats) + end +end + +# 2. Test error when non_mono_size < 1 +@testset "non_mono_size < 1 Error" begin + f(x) = (x[1] - 1)^2 + 4 * (x[2] - x[1]^2)^2 + nlp = ADNLPModel(f, [-1.2; 1.0]) + + @test_throws ErrorException begin + R2NSolver(nlp; non_mono_size = 0) + end + @test_throws ErrorException begin + R2N(nlp; non_mono_size = 0) + end + + @test_throws ErrorException begin + R2NSolver(nlp; non_mono_size = -1) + end + @test_throws ErrorException begin + R2N(nlp; non_mono_size = -1) + end +end + +# 3. Test error when subsolver_type is ShiftedLBFGSSolver but nlp is not of type LBFGSModel +@testset "ShiftedLBFGSSolver with wrong nlp type Error" begin + f(x) = (x[1] - 1)^2 + 4 * (x[2] - x[1]^2)^2 + nlp = ADNLPModel(f, [-1.2; 1.0]) + + @test_throws ErrorException begin + R2NSolver(nlp; subsolver= :shifted_lbfgs) + end + @test_throws ErrorException begin + R2N(nlp; subsolver= :shifted_lbfgs) + end +end + +# 4. Test error when subsolver_type is not a subtype of R2N_allowed_subsolvers +@testset "Invalid subsolver type Error" begin + f(x) = (x[1] - 1)^2 + 4 * (x[2] - x[1]^2)^2 + nlp = ADNLPModel(f, [-1.2; 1.0]) + + @test_throws ErrorException begin + R2NSolver(nlp; subsolver_type = CgSolver) + end + @test_throws ErrorException begin + R2N(nlp; subsolver_type = CgSolver) + end +end \ No newline at end of file diff --git a/my_R2N_dev/test_HSL.jl b/my_R2N_dev/test_HSL.jl new file mode 100644 index 00000000..badf6817 --- /dev/null +++ b/my_R2N_dev/test_HSL.jl @@ -0,0 +1,212 @@ +using Revise +# using ADNLPModels, JSOSolvers, LinearAlgebra +using NLPModels, ADNLPModels +using HSL_jll # ] dev C:\Users\Farhad\Documents\Github\openblas_HSL_jll.jl-2023.11.7\HSL_jll.jl-2023.11.7 +using HSL +using SparseMatricesCOO +using SparseArrays, LinearAlgebra + +# TO test HSL_jll +# TODO A version 2.0 of HSL_jll.jl based on a dummy libHSL has been precompiled with Yggdrasil. Therefore HSL_jll.jl is +# a registered Julia package and can be added as a dependency of any Julia package. +using HSL_jll +function LIBHSL_isfunctional() + @ccall libhsl.LIBHSL_isfunctional()::Bool +end +bool = LIBHSL_isfunctional() + +T = Float64 + +n = 4 +x = ones(T, 4) +nlp = ADNLPModel( + x -> sum(100 * (x[i + 1] - x[i]^2)^2 + (x[i] - 1)^2 for i = 1:(n - 1)), + x, + name = "Extended Rosenbrock", +) + +x = nlp.meta.x0 +# 1. get problem dimensions and Hessian structure +meta_nlp = nlp.meta +n = meta_nlp.nvar +nnzh = meta_nlp.nnzh + +# 2. Allocate COO arrays for the augmented matrix [H; sqrt(σ)I] +# Total non-zeros = non-zeros in Hessian (nnzh) + n diagonal entries for the identity block. +rows = Vector{Int}(undef, nnzh + n) +cols = Vector{Int}(undef, nnzh + n) +vals = Vector{T}(undef, nnzh + n) + +# 3. fill in the structure of the Hessian +hess_structure!(nlp, view(rows, 1:nnzh), view(cols, 1:nnzh)) + +# 4. Fill in the sparsity pattern for the √σ·Iₙ block +# This block lives in rows n+1 to n+n and columns n+1 to n+n. +@inbounds for i = 1:n + rows[nnzh + i] = n + i + cols[nnzh + i] = n + i +end + +# 5. Pre-allocate the augmented right-hand-side vector +b_aug = Vector{T}(undef, n + n) + +#testing the solver + +#1. update the hessian values +hess_coord!(nlp, x, view(vals, 1:nnzh)) +σ = 1.0 +@inbounds for i = 1:n + vals[nnzh + i] = σ +end + +b_aug[1:n] .= -1.0 +fill!(view(b_aug, (n + 1):(n + n)), zero(eltype(b_aug))) + +# Hx = SparseMatrixCOO( +# nnzh + n,#nvar, +# nnzh + n, #nvar, +# rows, +# cols, #ls_subsolver.jcn[1:ls_subsolver.nnzj], +# vals, #ls_subsolver.val[1:ls_subsolver.nnzj], +# ) +H = sparse(cols, rows, vals) #TODO add this also to the struct of MA97 + +LBL = Ma97(H) +ma97_factorize!(LBL) +x_sol = ma97_solve(LBL, b_aug) # or x = LBL \ rhs + +acc = norm(H * x_sol - b_aug) / norm(b_aug) +println("The accuracy of the solution is: $acc") + + + + +#### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### +# TEsting to update the values of H +new_vals = ones(nnzh + n) .* 10 + +H.nzval .= new_vals +#### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### + +LBL = Ma97(H) +ma97_factorize!(LBL) +x_sol = ma97_solve(LBL, b_aug) # or x = LBL \ rhs + +acc = norm(H * x_sol - b_aug) / norm(b_aug) +println("The accuracy of the solution is: $acc") + + + + + + + + + +##################OLD way + +abstract type AbstractMA97Solver end + +mutable struct MA97Solver{T} <: AbstractMA97Solver + # HSL MA97 factorization object + ma97_obj::Ma97{T} + + # Sparse coordinate representation for (B + σI) + rows::Vector{Int32} # MA97 prefers Int32 + cols::Vector{Int32} + vals::Vector{T} + + # Keep track of sizes + n::Int + nnzh::Int # Non-zeros in the Hessian B + + function MA97Solver(nlp::AbstractNLPModel{T, V}) where {T, V} + n = nlp.meta.nvar + nnzh = nlp.meta.nnzh + total_nnz = nnzh + n + + # 1. Build the coordinate structure of B + σI + rows = Vector{Int32}(undef, total_nnz) + cols = Vector{Int32}(undef, total_nnz) + hess_structure!(nlp, view(rows, 1:nnzh), view(cols, 1:nnzh)) + for i = 1:n + rows[nnzh + i] = i + cols[nnzh + i] = i + end + + # 2. Create a temporary sparse matrix ONCE to establish the final sparsity pattern. + # The values (ones) are just placeholders. + K_pattern = sparse(rows, cols, one(T), n, n) + + # 3. Initialize the Ma97 object using the final CSC pattern + ma97_obj = ma97_csc( + K_pattern.n, + Int32.(K_pattern.colptr), + Int32.(K_pattern.rowval), + K_pattern.nzval, # Pass initial values + ) + + # 4. Perform the expensive symbolic analysis here, only ONCE. + # ma97_analyze!(ma97_obj) + + # 5. Allocate a buffer for the values that will be updated in each iteration + vals = Vector{T}(undef, total_nnz) + + return new{T}(ma97_obj, rows, cols, vals, n, nnzh) + end +end + + +# Dispatch for MA97Solver +function subsolve!( + r2_subsolver::MA97Solver{T}, + R2N::R2NSolver, + nlp::AbstractNLPModel, + s, + atol, + n, + subsolver_verbose, +) where {T} + + # Unpack for clarity + g = R2N.gx # Note: R2N main loop has g = -∇f + σ = R2N.σ + n = r2_subsolver.n + nnzh = r2_subsolver.nnzh + + # 1. Update the Hessian part of the values array + hess_coord!(nlp, R2N.x, view(r2_subsolver.vals, 1:nnzh)) + + # 2. Update the shift part of the values array + # The last 'n' entries correspond to the diagonal for σI + @inbounds for i = 1:n + r2_subsolver.vals[nnzh + i] = σ + end + + # 3. Create the sparse matrix K = B + σI in CSC format. + # The `sparse` function sums up duplicate entries, which is exactly + # what we need for the diagonal. + #TODO with Prof. Orban, check if this is efficient + # K = sparse(r2_subsolver.rows, r2_subsolver.cols, r2_subsolver.vals, n, n) + + + # 4. Copy the new numerical values into the MA97 object + copyto!(r2_subsolver.ma97_obj.nzval, K.nzval) + + # 5. Factorize the matrix + #TODO Prof Orban, do I need this?# I think we only need to do this once + + ma97_factorize!(r2_subsolver.ma97_obj) + if r2_subsolver.ma97_obj.info.flag != 0 + @warn("MA97 factorization failed with flag = $(r2_subsolver.ma97_obj.info.flag)") + return false, :err, 1, 0 # Indicate failure + end + + # 6. Solve the system (B + σI)s = g, where g = -∇f + # s = ma97_solve(r2_subsolver.ma97_obj, g) # Solves in-place + # 6. Solve the system (B + σI)s = g, where g = -∇f + s .= g # Copy the right-hand-side (g) into the solution buffer (s) #TODO confirm with Prof. Orban + ma97_solve!(r2_subsolver.ma97_obj, s) # Solve in-place, overwriting s + + return true, :first_order, 1, 0 +end diff --git a/my_R2N_dev/test_R2N.jl b/my_R2N_dev/test_R2N.jl new file mode 100644 index 00000000..f1d6ab2c --- /dev/null +++ b/my_R2N_dev/test_R2N.jl @@ -0,0 +1,127 @@ +using Revise +using JSOSolvers +using HSL +using NLPModels, ADNLPModels +using SparseArrays, LinearAlgebra +using ADNLPModels, Krylov, LinearOperators, NLPModels, NLPModelsModifiers, SolverCore, SolverTools + +f(x) = (x[1] - 1)^2 + 4 * (x[2] - x[1]^2)^2 +nlp = ADNLPModel(f, [-1.2; 1.0]) +stats = R2N(nlp, verbose = 1) + println("The status is: ", stats.status == :first_order) + println("The solution is approximately: ", isapprox(stats.solution, [1.0; 1.0], atol = 1e-6)) + +s = :R2NSolver +fun = :R2N +f(x) = (x[1] - 1)^2 + 4 * (x[2] - x[1]^2)^2 +nlp = ADNLPModel(f, [-1.2; 1.0]) +# nlp.meta.x0 .= 2.0 + + +if fun == :R2N_exact + nlp = LBFGSModel(nlp) + solver = eval(s)(nlp,subsolver= :shifted_lbfgs) + else + solver = eval(s)(nlp) + end + + stats = GenericExecutionStats(nlp) + stats = SolverCore.solve!(solver, nlp, stats, verbose = 1) + println("The status is: ", stats.status == :first_order) + println("The solution is approximately: ", isapprox(stats.solution, [1.0; 1.0], atol = 1e-6)) + + nlp.meta.x0 .= 2.0 + SolverCore.reset!(solver) + + stats = GenericExecutionStats(nlp) + stats = SolverCore.solve!(solver, nlp, stats, verbose = 1) + println("The status is: ", stats.status == :first_order) + println("The solution is approximately: ", isapprox(stats.solution, [1.0; 1.0], atol = 1e-6)) + + + +# # # const npc_handler_allowed = [:armijo, :sigma, :prev, :cp] +# for mysub in [:cg, :cr, :minres, :minres_qlp] +# println("\n\n\t\t===================================") +# println("============ Testing subsolver: $mysub ==============\n\n") + +# for (name, mySolver) in [ +# ("R2N_cg_cp", (nlp, ; kwargs...) -> R2N(nlp; subsolver = mysub, npc_handler = :cp, kwargs...)), +# ( +# "R2N_cg_prev", +# (nlp, ; kwargs...) -> R2N(nlp; subsolver = mysub, npc_handler = :prev, kwargs...), +# ), +# ( +# "R2N_cg_sigma", +# (nlp, ; kwargs...) -> R2N(nlp; subsolver = mysub, npc_handler = :sigma, kwargs...), +# ), +# ( +# "R2N_cg_armijo", +# (nlp, ; kwargs...) -> R2N(nlp; subsolver = mysub, npc_handler = :armijo, kwargs...), +# ), +# ] +# println("Testing solver: $name,$mySolver") + +# println("===================================") +# f(x) = (x[1] - 1)^2 + 4 * (x[2] - x[1]^2)^2 +# nlp = ADNLPModel(f, [-1.2; 1.0]) + +# stats = mySolver(nlp, verbose = 1, max_iter = 14) +# println("The status is: ", stats.status == :first_order) +# println("The solution is approximately: ", isapprox(stats.solution, [1.0; 1.0], atol = 1e-6)) +# println("===================================") +# end +# end + +# # testing R2N defualt and testing different subsolvers Krylov +# for (name, mySolver) in [ +# ("R2N", (nlp, ; kwargs...) -> R2N(nlp; kwargs...)), +# ("R2N_cg", (nlp, ; kwargs...) -> R2N(nlp, subsolver = :cg; kwargs...)), +# ("R2N_cr", (nlp, ; kwargs...) -> R2N(nlp, subsolver = :cr; kwargs...)), +# ("R2N_minres", (nlp, ; kwargs...) -> R2N(nlp, subsolver = :minres; kwargs...)), +# ("R2N_minres_qlp", (nlp, ; kwargs...) -> R2N(nlp, subsolver = :minres_qlp; kwargs...)), +# ] +# println("Testing solver: $name,$mySolver") + +# println("===================================") +# f(x) = (x[1] - 1)^2 + 4 * (x[2] - x[1]^2)^2 +# nlp = ADNLPModel(f, [-1.2; 1.0]) + +# stats = mySolver(nlp, verbose = 1) +# println("The status is: ", stats.status == :first_order) +# println("The solution is approximately: ", isapprox(stats.solution, [1.0; 1.0], atol = 1e-6)) +# end + +# # testing LBFGS and Krylov + +# # TESITNG lbfgs EXACT SOLVER AS SUBSOLVER +# for (name, mySolver) in [ +# ("R2N_lbfgs", (nlp, ; kwargs...) -> R2N(LBFGSModel(nlp), subsolver = :shifted_lbfgs; kwargs...)), +# ] +# println("Testing solver: $name,$mySolver") + +# println("===================================") +# f(x) = (x[1] - 1)^2 + 4 * (x[2] - x[1]^2)^2 +# nlp = ADNLPModel(f, [-1.2; 1.0]) + +# stats = mySolver(nlp, verbose = 1) +# println("The status is: ", stats.status == :first_order) +# println("The solution is approximately: ", isapprox(stats.solution, [1.0; 1.0], atol = 1e-6)) +# end + +# # testing MA97 +# if LIBHSL_isfunctional() #TODO +# f(x) = (x[1] - 1)^2 + 4 * (x[2] - x[1]^2)^2 +# nlp = ADNLPModel(f, [-1.2; 1.0]) + +# stats = GenericExecutionStats(nlp) + +# solver = R2NSolver(nlp, subsolver = :ma97) +# stats = SolverCore.solve!(solver, nlp, stats, verbose = 1, max_iter = 50) +# println("The status is: ", stats.status == :first_order) +# println("The solution is approximately: ", isapprox(stats.solution, [1.0; 1.0], atol = 1e-6)) +# else +# @warn("HSL library is not functional. Skipping MA97 tests.") +# end + +# # testing different -ve handling strategies with Krylov solvers diff --git a/src/JSOSolvers.jl b/src/JSOSolvers.jl index 0e86dfd6..afb9c19e 100644 --- a/src/JSOSolvers.jl +++ b/src/JSOSolvers.jl @@ -1,7 +1,7 @@ module JSOSolvers # stdlib -using LinearAlgebra, Logging, Printf +using LinearAlgebra, Logging, Printf, SparseArrays # JSO packages using Krylov, @@ -10,6 +10,7 @@ using Krylov, import SolverCore.solve! export solve! + const Callback_docstring = " The callback is called at each iteration. The expected signature of the callback is `callback(nlp, solver, stats)`, and its output is ignored. @@ -45,9 +46,12 @@ end include("lbfgs.jl") include("trunk.jl") include("fomo.jl") +include("R2N.jl") # Unconstrained solvers for NLS include("trunkls.jl") +include("R2NLS.jl") + # List of keywords accepted by TRONTrustRegion const tron_keys = ( diff --git a/src/R2N.jl b/src/R2N.jl new file mode 100644 index 00000000..c8f1f996 --- /dev/null +++ b/src/R2N.jl @@ -0,0 +1,874 @@ +export R2N, R2NSolver, R2NParameterSet +export ShiftedLBFGSSolver, HSLDirectSolver +export ShiftedOperator + +using LinearOperators, LinearAlgebra +using SparseArrays +using HSL + +#TODO move to LinearOperators +# Define a new mutable operator for A = H + σI +mutable struct ShiftedOperator{T, V, OpH <: Union{AbstractLinearOperator{T}, AbstractMatrix{T}}} <: AbstractLinearOperator{T} + H::OpH + σ::T + n::Int + symmetric::Bool + hermitian::Bool +end + +function ShiftedOperator(H::OpH) where {T, OpH <: Union{AbstractLinearOperator{T}, AbstractMatrix{T}}} + is_sym = isa(H, AbstractLinearOperator) ? H.symmetric : issymmetric(H) + is_herm = isa(H, AbstractLinearOperator) ? H.hermitian : ishermitian(H) + return ShiftedOperator{T, Vector{T}, OpH}(H, zero(T), size(H, 1), is_sym, is_herm) +end + +# Define required properties for AbstractLinearOperator +Base.size(A::ShiftedOperator) = (A.n, A.n) +LinearAlgebra.isreal(A::ShiftedOperator{T}) where {T <: Real} = true +LinearOperators.issymmetric(A::ShiftedOperator) = A.symmetric +LinearOperators.ishermitian(A::ShiftedOperator) = A.hermitian + +# Define the core multiplication rules: y = (H + σI)x +function LinearAlgebra.mul!( + y::AbstractVecOrMat, + A::ShiftedOperator{T, V, OpH}, + x::AbstractVecOrMat{T}, +) where {T, V, OpH} + mul!(y, A.H, x) + LinearAlgebra.axpy!(A.σ, x, y) + return y +end + + +""" + R2NParameterSet([T=Float64]; θ1, θ2, η1, η2, γ1, γ2, γ3, σmin, non_mono_size) + +Parameter set for the R2N solver. Controls algorithmic tolerances and step acceptance. + +# Keyword Arguments +- `θ1 = T(0.5)`: Cauchy step parameter (0 < θ1 < 1). +- `θ2 = eps(T)^(-1)`: Cauchy step parameter. +- `η1 = eps(T)^(1/4)`: Step acceptance parameter (0 < η1 ≤ η2 < 1). +- `η2 = T(0.95)`: Step acceptance parameter (0 < η1 ≤ η2 < 1). +- `γ1 = T(1.5)`: Regularization update parameter (1 < γ1 ≤ γ2). +- `γ2 = T(2.5)`: Regularization update parameter (γ1 ≤ γ2). +- `γ3 = T(0.5)`: Regularization update parameter (0 < γ3 ≤ 1). +- `δ1 = T(0.5)`: Cauchy point calculation parameter. +- `σmin = eps(T)`: Minimum step parameter. +- `non_mono_size = 1`: the size of the non-monotone behaviour. If > 1, the algorithm will use a non-monotone strategy to accept steps. +""" +struct R2NParameterSet{T} <: AbstractParameterSet + θ1::Parameter{T, RealInterval{T}} + θ2::Parameter{T, RealInterval{T}} + η1::Parameter{T, RealInterval{T}} + η2::Parameter{T, RealInterval{T}} + γ1::Parameter{T, RealInterval{T}} + γ2::Parameter{T, RealInterval{T}} + γ3::Parameter{T, RealInterval{T}} + δ1::Parameter{T, RealInterval{T}} + σmin::Parameter{T, RealInterval{T}} + non_mono_size::Parameter{Int, IntegerRange{Int}} +end + +# Default parameter values +const R2N_θ1 = DefaultParameter(nlp -> eltype(nlp.meta.x0)(0.5), "T(0.5)") +const R2N_θ2 = DefaultParameter(nlp -> inv(eps(eltype(nlp.meta.x0))), "eps(T)^(-1)") +const R2N_η1 = DefaultParameter(nlp -> begin + T = eltype(nlp.meta.x0) + T(eps(T))^(T(1) / T(4)) +end, "eps(T)^(1/4)") +const R2N_η2 = DefaultParameter(nlp -> eltype(nlp.meta.x0)(0.95), "T(0.95)") +const R2N_γ1 = DefaultParameter(nlp -> eltype(nlp.meta.x0)(1.5), "T(1.5)") +const R2N_γ2 = DefaultParameter(nlp -> eltype(nlp.meta.x0)(2.5), "T(2.5)") +const R2N_γ3 = DefaultParameter(nlp -> eltype(nlp.meta.x0)(0.5), "T(0.5)") +const R2N_δ1 = DefaultParameter(nlp -> eltype(nlp.meta.x0)(0.5), "T(0.5)") +const R2N_σmin = DefaultParameter(nlp -> begin + T = eltype(nlp.meta.x0) + T(eps(T)) +end, "eps(T)") +const R2N_non_mono_size = DefaultParameter(1) + +function R2NParameterSet( + nlp::AbstractNLPModel; + θ1::T = get(R2N_θ1, nlp), + θ2::T = get(R2N_θ2, nlp), + η1::T = get(R2N_η1, nlp), + η2::T = get(R2N_η2, nlp), + γ1::T = get(R2N_γ1, nlp), + γ2::T = get(R2N_γ2, nlp), + γ3::T = get(R2N_γ3, nlp), + δ1::T = get(R2N_δ1, nlp), + σmin::T = get(R2N_σmin, nlp), + non_mono_size::Int = get(R2N_non_mono_size, nlp), +) where {T} + R2NParameterSet{T}( + Parameter(θ1, RealInterval(zero(T), one(T))), + Parameter(θ2, RealInterval(one(T), T(Inf))), + Parameter(η1, RealInterval(zero(T), one(T))), + Parameter(η2, RealInterval(zero(T), one(T))), + Parameter(γ1, RealInterval(one(T), T(Inf))), + Parameter(γ2, RealInterval(one(T), T(Inf))), + Parameter(γ3, RealInterval(zero(T), one(T))), + Parameter(δ1, RealInterval(zero(T), one(T))), + Parameter(σmin, RealInterval(zero(T), T(Inf))), + Parameter(non_mono_size, IntegerRange(1, typemax(Int))), + ) +end + +abstract type AbstractShiftedLBFGSSolver end + +mutable struct ShiftedLBFGSSolver <: AbstractShiftedLBFGSSolver #TODO Ask what I can do inside here + # Shifted LBFGS-specific fields +end + +abstract type AbstractMASolver end + +""" +HSLDirectSolver: Generic wrapper for HSL direct solvers (e.g., MA97, MA57). +- hsl_obj: The HSL solver object (e.g., Ma97 or Ma57) +- rows, cols, vals: COO format for the Hessian + shift +- n: problem size +- nnzh: number of Hessian nonzeros +""" +mutable struct HSLDirectSolver{T, S} <: AbstractMASolver + hsl_obj::S + rows::Vector{Int} + cols::Vector{Int} + vals::Vector{T} + n::Int + nnzh::Int +end + +""" + HSLDirectSolver(nlp, hsl_constructor) + +Constructs an HSLDirectSolver for the given NLP model and HSL solver constructor (e.g., ma97_coord or ma57_coord). +""" +function HSLDirectSolver(nlp::AbstractNLPModel{T, V}, hsl_constructor) where {T, V} + n = nlp.meta.nvar + nnzh = nlp.meta.nnzh + total_nnz = nnzh + n + + rows = Vector{Int}(undef, total_nnz) + cols = Vector{Int}(undef, total_nnz) + vals = zeros(T, total_nnz) + + hess_structure!(nlp, view(rows, 1:nnzh), view(cols, 1:nnzh)) + hess_coord!(nlp, nlp.meta.x0, view(vals, 1:nnzh)) + + @inbounds for i = 1:n + rows[nnzh + i] = i + cols[nnzh + i] = i + vals[nnzh + i] = one(T) + end + + hsl_obj = hsl_constructor(n, cols, rows, vals) + return HSLDirectSolver{T, typeof(hsl_obj)}(hsl_obj, rows, cols, vals, n, nnzh) +end + + +const npc_handler_allowed = [:armijo, :sigma, :prev, :cp] + +const R2N_allowed_subsolvers = [:cg, :cr, :minres, :minres_qlp, :shifted_lbfgs, :ma97, :ma57] #TODO Prof. Orban check if I can add more + +""" + R2N(nlp; kwargs...) + +An inexact second-order quadratic regularization method for unconstrained optimization (with shifted L-BFGS or shifted Hessian operator). + +For advanced usage, first define a `R2NSolver` to preallocate the memory used in the algorithm, and then call `solve!`: + + solver = R2NSolver(nlp) + solve!(solver, nlp; kwargs...) + +# Arguments +- `nlp::AbstractNLPModel{T, V}` is the model to solve, see `NLPModels.jl`. + + + +# Keyword arguments +- `params::R2NParameterSet = R2NParameterSet(nlp)`: algorithm parameters, see [`R2NParameterSet`](@ref). +- `η1::T = $(R2N_η1)`: step acceptance parameter, see [`R2NParameterSet`](@ref). +- `η2::T = $(R2N_η2)`: step acceptance parameter, see [`R2NParameterSet`](@ref). +- `θ1::T = $(R2N_θ1)`: Cauchy step parameter, see [`R2NParameterSet`](@ref). +- `θ2::T = $(R2N_θ2)`: Cauchy step parameter, see [`R2NParameterSet`](@ref). +- `γ1::T = $(R2N_γ1)`: regularization update parameter, see [`R2NParameterSet`](@ref). +- `γ2::T = $(R2N_γ2)`: regularization update parameter, see [`R2NParameterSet`](@ref). +- `γ3::T = $(R2N_γ3)`: regularization update parameter, see [`R2NParameterSet`](@ref). +- `δ1::T = $(R2N_δ1)`: Cauchy point calculation parameter, see [`R2NParameterSet`](@ref). +- `σmin::T = $(R2N_σmin)`: minimum step parameter, see [`R2NParameterSet`](@ref). +- `non_mono_size::Int = $(R2N_non_mono_size)`: the size of the non-monotone behaviour. If > 1, the algorithm will use a non-monotone strategy to accept steps. +- `x::V = nlp.meta.x0`: the initial guess. +- `atol::T = √eps(T)`: absolute tolerance. +- `rtol::T = √eps(T)`: relative tolerance: algorithm stops when ‖∇f(xᵏ)‖ ≤ atol + rtol * ‖∇f(x⁰)‖. +- `max_eval::Int = -1`: maximum number of evaluations of the objective function. +- `max_time::Float64 = 30.0`: maximum time limit in seconds. +- `max_iter::Int = typemax(Int)`: maximum number of iterations. +- `verbose::Int = 0`: if > 0, display iteration details every `verbose` iteration. +- `subsolver::Symbol = :shifted_lbfgs`: the subsolver to solve the shifted system. The `MinresWorkspace` which solves the shifted linear system exactly at each iteration. Using the exact solver is only possible if `nlp` is an `LBFGSModel`. See `JSOSolvers.R2N_allowed_subsolvers` for a list of available subsolvers. +- `subsolver_verbose::Int = 0`: if > 0, display iteration information every `subsolver_verbose` iteration of the subsolver if KrylovWorkspace type is selected. +- `scp_flag::Bool = true`: if true, we compare the norm of the calculate step with `θ2 * norm(scp)`, each iteration, selecting the smaller step. +- `npc_handler::Symbol = :armijo`: the non_positive_curve handling strategy. + - `:armijo`: uses the Armijo rule to handle non-positive curvature. + - `:sigma`: increase the regularization parameter σ. + - `:prev`: if subsolver return after first iteration, increase the sigma, but if subsolver return after second iteration, set s_k = s_k^(t-1). + - `:cp`: set s_k to Cauchy point. +See `JSOSolvers.npc_handler_allowed` for a list of available `npc_handler` strategies. + +All algorithmic parameters (θ1, θ2, η1, η2, γ1, γ2, γ3, σmin) can be set via the `params` keyword or individually as shown above. If both are provided, individual keywords take precedence. + +# Output +The value returned is a `GenericExecutionStats`, see `SolverCore.jl`. + +# Callback +$(Callback_docstring) + +# Examples #TODO +```jldoctest +using JSOSolvers, ADNLPModels +nlp = ADNLPModel(x -> sum(x.^2), ones(3)) +stats = R2N(nlp) + +# output + +"Execution stats: first-order stationary" +``` + +```jldoctest +using JSOSolvers, ADNLPModels +nlp = ADNLPModel(x -> sum(x.^2), ones(3)) +solver = R2NSolver(nlp); +stats = solve!(solver, nlp) + +# output + +"Execution stats: first-order stationary" +``` +""" + +mutable struct R2NSolver{ + T, + V, + Op <: Union{AbstractLinearOperator{T}, AbstractMatrix{T}}, #TODO confirm with Prof. Orban + ShiftedOp <: Union{ShiftedOperator{T, V, Op}, Nothing}, # for cg and cr solvers + Sub <: Union{KrylovWorkspace{T, T, V}, ShiftedLBFGSSolver, HSLDirectSolver{T, S} where S}, +} <: AbstractOptimizationSolver + x::V + xt::V + gx::V + gn::V + H::Op + A::ShiftedOp + Hs::V + s::V + scp::V + obj_vec::V # used for non-monotone + r2_subsolver::Sub + subtol::T + σ::T + params::R2NParameterSet{T} +end + +function R2NSolver( + nlp::AbstractNLPModel{T, V}; + η1 = get(R2N_η1, nlp), + η2 = get(R2N_η2, nlp), + θ1 = get(R2N_θ1, nlp), + θ2 = get(R2N_θ2, nlp), + γ1 = get(R2N_γ1, nlp), + γ2 = get(R2N_γ2, nlp), + γ3 = get(R2N_γ3, nlp), + δ1 = get(R2N_δ1, nlp), + σmin = get(R2N_σmin, nlp), + non_mono_size = get(R2N_non_mono_size, nlp), + subsolver::Symbol = :minres, +) where {T, V} + params = R2NParameterSet( + nlp; + η1 = η1, + η2 = η2, + θ1 = θ1, + θ2 = θ2, + γ1 = γ1, + γ2 = γ2, + γ3 = γ3, + δ1 = δ1, + σmin = σmin, + non_mono_size = non_mono_size, + ) + subsolver in R2N_allowed_subsolvers || + error("subproblem solver must be one of $(R2N_allowed_subsolvers)") + + value(params.non_mono_size) >= 1 || error("non_mono_size must be greater than or equal to 1") + + !(subsolver == :shifted_lbfgs) || + (nlp isa LBFGSModel) || + error("Unsupported subsolver type, ShiftedLBFGSSolver can only be used by LBFGSModel") + + nvar = nlp.meta.nvar + x = V(undef, nvar) + xt = V(undef, nvar) + gx = V(undef, nvar) + gn = isa(nlp, QuasiNewtonModel) ? V(undef, nvar) : V(undef, 0) + Hs = V(undef, nvar) + A = nothing + + local H, r2_subsolver + if subsolver == :ma97 + LIBHSL_isfunctional() || error("HSL library is not functional") + r2_subsolver = HSLDirectSolver(nlp, ma97_coord) + H = spzeros(T, nvar, nvar) + elseif subsolver == :ma57 + LIBHSL_isfunctional() || error("HSL library is not functional") + r2_subsolver = HSLDirectSolver(nlp, ma57_coord) + H = spzeros(T, nvar, nvar) + + else # not using ma971 + # H = isa(nlp, QuasiNewtonModel) ? nlp.op : hess_op!(nlp, x, Hs) + if subsolver == :shifted_lbfgs + H = nlp.op + r2_subsolver = ShiftedLBFGSSolver() + else + H = hess_op!(nlp, x, Hs) + r2_subsolver = krylov_workspace(Val(subsolver), nvar, nvar, V) + if subsolver in (:cg, :cr) + A = ShiftedOperator(H) + end + end + end + + Op = typeof(H) + ShiftedOp = typeof(A) + σ = zero(T) + s = V(undef, nvar) + scp = V(undef, nvar) + subtol = one(T) + obj_vec = fill(typemin(T), non_mono_size) + Sub = typeof(r2_subsolver) + + return R2NSolver{T, V, Op, ShiftedOp, Sub}( + x, + xt, + gx, + gn, + H, + A, + Hs, + s, + scp, + obj_vec, + r2_subsolver, + subtol, + σ, + params, + ) +end +#TODO Prof. Orban, check if I need to reset H in reset! function +function SolverCore.reset!(solver::R2NSolver{T}) where {T} + fill!(solver.obj_vec, typemin(T)) + reset!(solver.H) + if solver.r2_subsolver isa CgWorkspace || solver.r2_subsolver isa CrWorkspace + solver.A = ShiftedOperator(solver.H) + end + solver +end +function SolverCore.reset!(solver::R2NSolver{T}, nlp::AbstractNLPModel) where {T} + fill!(solver.obj_vec, typemin(T)) + # @assert (length(solver.gn) == 0) || isa(nlp, QuasiNewtonModel) + # solver.H = isa(nlp, QuasiNewtonModel) ? nlp.op : hess_op!(nlp, solver.x, solver.Hs) + solver +end + +@doc (@doc R2NSolver) function R2N( + nlp::AbstractNLPModel{T, V}; + η1::Real = get(R2N_η1, nlp), + η2::Real = get(R2N_η2, nlp), + θ1::Real = get(R2N_θ1, nlp), + θ2::Real = get(R2N_θ2, nlp), + γ1::Real = get(R2N_γ1, nlp), + γ2::Real = get(R2N_γ2, nlp), + γ3::Real = get(R2N_γ3, nlp), + δ1::Real = get(R2N_δ1, nlp), + σmin::Real = get(R2N_σmin, nlp), + non_mono_size::Int = get(R2N_non_mono_size, nlp), + subsolver::Symbol = :minres, + kwargs..., +) where {T, V} + solver = R2NSolver( + nlp; + η1 = convert(T, η1), + η2 = convert(T, η2), + θ1 = convert(T, θ1), + θ2 = convert(T, θ2), + γ1 = convert(T, γ1), + γ2 = convert(T, γ2), + γ3 = convert(T, γ3), + δ1 = convert(T, δ1), + σmin = convert(T, σmin), + non_mono_size = non_mono_size, + subsolver = subsolver, + ) + return solve!(solver, nlp; kwargs...) +end + +function SolverCore.solve!( + solver::R2NSolver{T, V}, + nlp::AbstractNLPModel{T, V}, + stats::GenericExecutionStats{T, V}; + callback = (args...) -> nothing, + x::V = nlp.meta.x0, + atol::T = √eps(T), + rtol::T = √eps(T), + max_time::Float64 = 30.0, + max_eval::Int = -1, + max_iter::Int = typemax(Int), + verbose::Int = 0, + subsolver_verbose::Int = 0, + npc_handler::Symbol = :armijo, + scp_flag::Bool = true, +) where {T, V} + unconstrained(nlp) || error("R2N should only be called on unconstrained problems.") + npc_handler in npc_handler_allowed || error("npc_handler must be one of $(npc_handler_allowed)") + + reset!(stats) + params = solver.params + η1 = value(params.η1) + η2 = value(params.η2) + θ1 = value(params.θ1) + θ2 = value(params.θ2) + γ1 = value(params.γ1) + γ2 = value(params.γ2) + γ3 = value(params.γ3) + δ1 = value(params.δ1) + σmin = value(params.σmin) + non_mono_size = value(params.non_mono_size) + + + @assert(η1 > 0 && η1 < 1) + @assert(θ1 > 0 && θ1 < 1) + @assert(θ2 > 1) + @assert(γ1 >= 1 && γ1 <= γ2 && γ3 <= 1) + @assert(δ1 > 0 && δ1 < 1) + + start_time = time() + set_time!(stats, 0.0) + + n = nlp.meta.nvar + x = solver.x .= x + xt = solver.xt + ∇fk = solver.gx # k-1 + ∇fn = solver.gn #current + s = solver.s + scp = solver.scp + H = solver.H + A = solver.A + Hs = solver.Hs + σk = solver.σ + + subtol = solver.subtol + subsolver_solved = false + + set_iter!(stats, 0) + f0 = obj(nlp, x) + set_objective!(stats, f0) + + grad!(nlp, x, ∇fk) + isa(nlp, QuasiNewtonModel) && (∇fn .= ∇fk) + norm_∇fk = norm(∇fk) + set_dual_residual!(stats, norm_∇fk) + + σk = 2^round(log2(norm_∇fk + 1)) / norm_∇fk + ρk = zero(T) + + # Stopping criterion: + fmin = min(-one(T), f0) / eps(T) + unbounded = f0 < fmin + + ϵ = atol + rtol * norm_∇fk + optimal = norm_∇fk ≤ ϵ + if optimal + @info "Optimal point found at initial point" + @info log_header( + [:iter, :f, :dual, :σ, :ρ], + [Int, Float64, Float64, Float64, Float64], + hdr_override = Dict(:f => "f(x)", :dual => "‖∇f‖"), + ) + @info log_row([stats.iter, stats.objective, norm_∇fk, σk, ρk]) + end + cp_step_log = " " + if verbose > 0 && mod(stats.iter, verbose) == 0 + @info log_header( + [:iter, :f, :dual, :σ, :ρ, :sub_iter, :dir, :cp_step_log, :sub_status], + [Int, Float64, Float64, Float64, Float64, Int, String, String, String], + hdr_override = Dict( + :f => "f(x)", + :dual => "‖∇f‖", + :sub_iter => "subiter", + :dir => "dir", + :cp_step_log => "cp step", + :sub_status => "status", + ), + ) + @info log_row([stats.iter, stats.objective, norm_∇fk, σk, ρk, 0, " ", " ", " "]) + end + + set_status!( + stats, + get_status( + nlp, + elapsed_time = stats.elapsed_time, + optimal = optimal, + unbounded = unbounded, + max_eval = max_eval, + iter = stats.iter, + max_iter = max_iter, + max_time = max_time, + ), + ) + + # subtol initialization for subsolver + subtol = max(rtol, min(T(0.1), √norm_∇fk, T(0.9) * subtol)) + solver.σ = σk + solver.subtol = subtol + r2_subsolver = solver.r2_subsolver + + callback(nlp, solver, stats) + subtol = solver.subtol + σk = solver.σ + + done = stats.status != :unknown + + ν_k = one(T) # used for scp calculation + γ_k = zero(T) + scp_recal = true + + while !done + scp_recal = true # Recalculate the Cauchy point if needed + + # Compute the Cauchy step. + # Note that we use buffer values Hs to avoid reallocating memory. + # mul!(Hs, H, ∇fk) #TODO check for HSL + if r2_subsolver isa HSLDirectSolver + coo_sym_prod!(r2_subsolver.rows, r2_subsolver.cols, r2_subsolver.vals, ∇fk, Hs) + else + mul!(Hs, H, ∇fk) # Use linear operator + end + curv = dot(∇fk, Hs) + slope = σk * norm_∇fk^2 # slope= σ * ||∇f||^2 + γ_k = (curv + slope) / norm_∇fk^2 + if γ_k > 0 + cp_step_log = "α_k" + ν_k = 2*(1-δ1) / (γ_k) + else + # we have to calcualte the scp, since we have encounter a negative curvature + if H isa AbstractLinearOperator + λmax, found_λ = opnorm(H) # This uses iterative methods (p=2) + else + # H is a SparseMatrixCSC (MA97 case), use p=Inf as suggested by the error + λmax = opnorm(H, Inf) + found_λ = true # We assume the Inf-norm was found + end + cp_step_log = "ν_k" + ν_k = θ1 / (λmax + σk) + end + + # Solving for step direction s_k + ∇fk .*= -1 + if r2_subsolver isa CgWorkspace || r2_subsolver isa CrWorkspace + # Update the shift in the operator + solver.A.σ = σk + solver.H = H + end + subsolver_solved, sub_stats, subiter, npcCount = + subsolve!(r2_subsolver, solver, nlp, s, zero(T), n, subsolver_verbose) + + if !subsolver_solved + @warn("Subsolver failed to solve the system") + # It might be better to increase σ and retry instead of stopping, + # but for now, we'll stop. + break + end + if !(r2_subsolver isa ShiftedLBFGSSolver) && !(r2_subsolver isa HSLDirectSolver) + if r2_subsolver.stats.npcCount >= 1 #npc case + if npc_handler == :armijo #TODO check this logic with Prof. Orban + c = one(T) * 1e-4 #TODO do I want user to set this value? + α = one(T) * 1.0 + increase_factor = one(T) * 1.5 + decrease_factor = one(T) * 0.5 + min_alpha = one(T) * 1e-8 + max_alpha = one(T) * 1e2 + dir = r2_subsolver.npc_dir + f0 = stats.objective + grad_dir = dot(∇fk, dir) + # First, try α = 1 + xt .= x .+ α * dir + fck = obj(nlp, xt) + if fck > f0 + c * α * grad_dir + # Backward tracking + while fck > f0 + c * α * grad_dir && α > min_alpha + α *= decrease_factor + xt .= x .+ α * dir + fck = obj(nlp, xt) + end + else + # Forward tracking + while true + α_new = α * increase_factor + xt .= x .+ α_new * dir + fck_new = obj(nlp, xt) + if fck_new > f0 + c * α_new * grad_dir || α_new > max_alpha + break + end + α = α_new + fck = fck_new + end + end + s .= α * dir + elseif npc_handler == :prev #Cr and cg will return the last iteration s + s .= r2_subsolver.x + elseif npc_handler == :cp + # Cauchy point + scp .= ν_k * ∇fk # the ∇fk is already negative + s .= scp + scp_recal = false # we don't need to recalculate the scp later + end + end + end + + if scp_flag && scp_recal # if the case was cp, then s = scp already + # Based on the flag, scp is calcualted + scp .= ν_k * ∇fk # the ∇fk is already negative + if norm(s) > θ2 * norm(scp) + s .= scp + end + end + if npc_handler == :sigma && r2_subsolver.stats.npcCount >= 1 # non-positive curvature case happen and the npc_handler is sigma + step_accepted = false + else + # Correctly compute curvature s' * B * s + #TODO Prof Orban, check if this is correct + if solver.r2_subsolver isa HSLDirectSolver + coo_sym_prod!(solver.r2_subsolver.rows, solver.r2_subsolver.cols, solver.r2_subsolver.vals, s, Hs) + else + mul!(Hs, H, s) # Use linear operator + end + + curv = dot(s, Hs) + slope = dot(s, ∇fk) # = -∇fkᵀ s because we flipped the sign of ∇fk + + ΔTk = slope - curv / 2 + xt .= x .+ s + fck = obj(nlp, xt) + + if non_mono_size > 1 #non-monotone behaviour + k = mod(stats.iter, non_mono_size) + 1 + solver.obj_vec[k] = stats.objective + fck_max = maximum(solver.obj_vec) + ρk = (fck_max - fck) / (fck_max - stats.objective + ΔTk) #TODO Prof. Orban check if this is correct the denominator + else + # Avoid division by zero/negative. If ΔTk <= 0, the model is bad. + ρk = (ΔTk > 10 * eps(T)) ? (stats.objective - fck) / ΔTk : -one(T) + end + + # Update regularization parameters and Acceptance of the new candidate + step_accepted = ρk >= η1 + if step_accepted + x .= xt + grad!(nlp, x, ∇fk) + if isa(nlp, QuasiNewtonModel) + ∇fn .-= ∇fk + ∇fn .*= -1 # = ∇f(xₖ₊₁) - ∇f(xₖ) + push!(nlp, s, ∇fn) + ∇fn .= ∇fk + end + set_objective!(stats, fck) + unbounded = fck < fmin + norm_∇fk = norm(∇fk) + if ρk >= η2 + σk = max(σmin, γ3 * σk) + else # η1 ≤ ρk < η2 + σk = min(σmin, γ1 * σk) + end + # we need to update H if we use Ma97 or ma57 + if solver.r2_subsolver isa HSLDirectSolver + hess_coord!(nlp, x, view(solver.r2_subsolver.vals, 1:solver.r2_subsolver.nnzh)) + end + else # η1 > ρk + σk = max(σmin, γ2 * σk) + ∇fk .*= -1 + end + end + + set_iter!(stats, stats.iter + 1) + set_time!(stats, time() - start_time) + + subtol = max(rtol, min(T(0.1), √norm_∇fk, T(0.9) * subtol)) + set_dual_residual!(stats, norm_∇fk) + + solver.σ = σk + solver.subtol = subtol + + callback(nlp, solver, stats) + + norm_∇fk = stats.dual_feas # if the user change it, they just change the stats.norm , they also have to change subtol + σk = solver.σ + subtol = solver.subtol + optimal = norm_∇fk ≤ ϵ + + if verbose > 0 && mod(stats.iter, verbose) == 0 + dir_stat = step_accepted ? "↘" : "↗" + @info log_row([ + stats.iter, + stats.objective, + norm_∇fk, + σk, + ρk, + subiter, + cp_step_log, + dir_stat, + sub_stats, + ]) + end + + if stats.status == :user + done = true + else + set_status!( + stats, + get_status( + nlp, + elapsed_time = stats.elapsed_time, + optimal = optimal, + unbounded = unbounded, + max_eval = max_eval, + iter = stats.iter, + max_iter = max_iter, + max_time = max_time, + ), + ) + done = stats.status != :unknown + end + end + + set_solution!(stats, x) + return stats +end + +# Dispatch for subsolvers KrylovWorkspace: cg and cr +function subsolve!( + r2_subsolver::KrylovWorkspace{T, T, V}, + R2N::R2NSolver, + nlp::AbstractNLPModel, + s, + atol, + n, + subsolver_verbose, +) where {T, V} + krylov_solve!( + r2_subsolver, + R2N.A, # Use the ShiftedOperator A + R2N.gx, + itmax = max(2 * n, 50), + atol = atol, + rtol = R2N.subtol, + verbose = subsolver_verbose, + linesearch = true, + ) + s .= r2_subsolver.x + return Krylov.issolved(r2_subsolver), + r2_subsolver.stats.status, + r2_subsolver.stats.niter, + r2_subsolver.stats.npcCount +end + +# Dispatch for MinresWorkspace and MinresQlpWorkspace +function subsolve!( + r2_subsolver::Union{MinresWorkspace{T, V}, MinresQlpWorkspace{T, V}}, + R2N::R2NSolver, + nlp::AbstractNLPModel, + s, + atol, + n, + subsolver_verbose, +) where {T, V} + krylov_solve!( + r2_subsolver, + R2N.H, + R2N.gx, + λ = R2N.σ, + itmax = max(2 * n, 50), + atol = atol, + rtol = R2N.subtol, + verbose = subsolver_verbose, + linesearch = true, + ) + s .= r2_subsolver.x + return Krylov.issolved(r2_subsolver), + r2_subsolver.stats.status, + r2_subsolver.stats.niter, + r2_subsolver.stats.npcCount +end + +# Dispatch for ShiftedLBFGSSolver +function subsolve!( + r2_subsolver::ShiftedLBFGSSolver, + R2N::R2NSolver, + nlp::AbstractNLPModel, + s, + atol, + n, + subsolver_verbose, +) + ∇f_neg = R2N.gx + H = R2N.H + σ = R2N.σ + solve_shifted_system!(s, H, ∇f_neg, σ) + return true, :first_order, 1, 0 +end + +# Dispatch for HSLDirectSolver +""" + subsolve!(r2_subsolver::HSLDirectSolver, ...) + +Solves the shifted system using the selected HSL direct solver (MA97 or MA57). +""" +function subsolve!( + r2_subsolver::HSLDirectSolver{T, S}, + R2N::R2NSolver, + nlp::AbstractNLPModel, + s, + atol, + n, + subsolver_verbose, +) where {T, S} + + g = R2N.gx + σ = R2N.σ + + @inbounds for i = 1:n + r2_subsolver.vals[r2_subsolver.nnzh + i] = σ + end + + # Dispatch to correct factorization/solve based on S + if S <: Ma97{T} + ma97_factorize!(r2_subsolver.hsl_obj) + if r2_subsolver.hsl_obj.info.flag != 0 + @warn("MA97 factorization failed with flag = $(r2_subsolver.hsl_obj.info.flag)") + return false, :err, 1, 0 + end + s .= g + ma97_solve!(r2_subsolver.hsl_obj, s) + elseif S <: Ma57{T} + ma57_factorize!(r2_subsolver.hsl_obj) + if r2_subsolver.hsl_obj.info.flag != 0 + @warn("MA57 factorization failed with flag = $(r2_subsolver.hsl_obj.info.flag)") + return false, :err, 1, 0 + end + s .= g + ma57_solve!(r2_subsolver.hsl_obj, s) + else + error("Unknown HSL solver type") + end + + return true, :first_order, 1, 0 +end \ No newline at end of file diff --git a/src/R2NLS.jl b/src/R2NLS.jl new file mode 100644 index 00000000..49f8230c --- /dev/null +++ b/src/R2NLS.jl @@ -0,0 +1,706 @@ +export R2NLS, R2SolverNLS, R2NLSParameterSet +export QRMumpsSolver +using SparseMatricesCOO + +using QRMumps, LinearAlgebra, SparseArrays + +#TODO prof Orban, the name should be R2SolverNLS or R2NSolverNLS +""" + R2NLSParameterSet([T=Float64]; η1, η2, θ1, θ2, γ1, γ2, γ3, δ1, σmin, non_mono_size) + +Parameter set for the R2NLS solver. Controls algorithmic tolerances and step acceptance. + +# Keyword Arguments +- `η1 = eps(T)^(1/4)`: Step acceptance parameter. +- `η2 = T(0.95)`: Step acceptance parameter. +- `θ1 = T(0.5)`: Cauchy step parameter. +- `θ2 = eps(T)^(-1)`: Cauchy step parameter. +- `γ1 = T(1.5)`: Regularization update parameter. +- `γ2 = T(2.5)`: Regularization update parameter. +- `γ3 = T(0.5)`: Regularization update parameter. +- `δ1 = T(0.5)`: Cauchy point calculation parameter. +- `σmin = eps(T)`: Minimum step parameter. +- `non_mono_size = 1`: the size of the non-monotone behaviour. If > 1, the algorithm will use a non-monotone strategy to accept steps. +""" +struct R2NLSParameterSet{T} <: AbstractParameterSet + η1::Parameter{T, RealInterval{T}} + η2::Parameter{T, RealInterval{T}} + θ1::Parameter{T, RealInterval{T}} + θ2::Parameter{T, RealInterval{T}} + γ1::Parameter{T, RealInterval{T}} + γ2::Parameter{T, RealInterval{T}} + γ3::Parameter{T, RealInterval{T}} + δ1::Parameter{T, RealInterval{T}} + σmin::Parameter{T, RealInterval{T}} + non_mono_size::Parameter{Int, IntegerRange{Int}} +end + +# Default parameter values +const R2NLS_η1 = DefaultParameter(nlp -> begin + T = eltype(nlp.meta.x0) + T(eps(T))^(T(1)/T(4)) +end, "eps(T)^(1/4)") +const R2NLS_η2 = DefaultParameter(nlp -> eltype(nlp.meta.x0)(0.95), "T(0.95)") +const R2NLS_θ1 = DefaultParameter(nlp -> eltype(nlp.meta.x0)(0.5), "T(0.5)") +const R2NLS_θ2 = DefaultParameter(nlp -> inv(eps(eltype(nlp.meta.x0))), "eps(T)^(-1)") +const R2NLS_γ1 = DefaultParameter(nlp -> eltype(nlp.meta.x0)(1.5), "T(1.5)") +const R2NLS_γ2 = DefaultParameter(nlp -> eltype(nlp.meta.x0)(2.5), "T(2.5)") +const R2NLS_γ3 = DefaultParameter(nlp -> eltype(nlp.meta.x0)(0.5), "T(0.5)") +const R2NLS_δ1 = DefaultParameter(nlp -> eltype(nlp.meta.x0)(0.5), "T(0.5)") +const R2NLS_σmin = DefaultParameter(nlp -> eps(eltype(nlp.meta.x0)), "eps(T)") +const R2NLS_non_mono_size = DefaultParameter(1) + +function R2NLSParameterSet( + nlp::AbstractNLSModel; + η1::T = get(R2NLS_η1, nlp), + η2::T = get(R2NLS_η2, nlp), + θ1::T = get(R2NLS_θ1, nlp), + θ2::T = get(R2NLS_θ2, nlp), + γ1::T = get(R2NLS_γ1, nlp), + γ2::T = get(R2NLS_γ2, nlp), + γ3::T = get(R2NLS_γ3, nlp), + δ1::T = get(R2NLS_δ1, nlp), + σmin::T = get(R2NLS_σmin, nlp), + non_mono_size::Int = get(R2NLS_non_mono_size, nlp), +) where {T} + R2NLSParameterSet{T}( + Parameter(η1, RealInterval(zero(T), one(T))), + Parameter(η2, RealInterval(zero(T), one(T))), + Parameter(θ1, RealInterval(zero(T), one(T))), + Parameter(θ2, RealInterval(one(T), T(Inf))), + Parameter(γ1, RealInterval(one(T), T(Inf))), + Parameter(γ2, RealInterval(one(T), T(Inf))), + Parameter(γ3, RealInterval(zero(T), one(T))), + Parameter(δ1, RealInterval(zero(T), one(T))), + Parameter(σmin, RealInterval(zero(T), T(Inf))), + Parameter(non_mono_size, IntegerRange(1, typemax(Int))), + ) +end + +abstract type AbstractQRMumpsSolver end + +""" + QRMumpsSolver + +A solver structure for handling the linear least-squares subproblems within R2NLS +using the QRMumps package. This structure pre-allocates all necessary memory +for the sparse matrix representation and the factorization. +""" +mutable struct QRMumpsSolver{T} <: AbstractQRMumpsSolver + # QRMumps structures + spmat::qrm_spmat{T} + spfct::qrm_spfct{T} + + # COO storage for the augmented matrix [J; sqrt(σ) * I] + irn::Vector{Int} + jcn::Vector{Int} + val::Vector{T} + + # Augmented RHS vector + b_aug::Vector{T} + + # Problem dimensions + m::Int + n::Int + nnzj::Int + + closed::Bool # Avoid double-destroy + + function QRMumpsSolver(nlp::AbstractNLSModel{T}) where {T} + # Safely initialize QRMumps library + qrm_init() + + # 1. Get problem dimensions and Jacobian structure + meta_nls = nls_meta(nlp) + n = nlp.nls_meta.nvar + m = nlp.nls_meta.nequ + nnzj = meta_nls.nnzj + + # 2. Allocate COO arrays for the augmented matrix [J; sqrt(σ)I] + # Total non-zeros = non-zeros in Jacobian (nnzj) + n diagonal entries for the identity block. + irn = Vector{Int}(undef, nnzj + n) + jcn = Vector{Int}(undef, nnzj + n) + val = Vector{T}(undef, nnzj + n) + + # 3. Fill in the sparsity pattern of the Jacobian J(x) + jac_structure_residual!(nlp, view(irn, 1:nnzj), view(jcn, 1:nnzj)) + + # 4. Fill in the sparsity pattern for the √σ·Iₙ block + # This block lives in rows m+1 to m+n and columns 1 to n. + @inbounds for i = 1:n + irn[nnzj + i] = m + i + jcn[nnzj + i] = i + end + + # 5. Initialize QRMumps sparse matrix and factorization structures + spmat = qrm_spmat_init(m + n, n, irn, jcn, val; sym = false) + spfct = qrm_spfct_init(spmat) + qrm_analyse!(spmat, spfct; transp = 'n') + + # 6. Pre-allocate the augmented right-hand-side vector + b_aug = Vector{T}(undef, m+n) + + # 7. Create the solver object and set a finalizer for safe cleanup. + solver = new{T}(spmat, spfct, irn, jcn, val, b_aug, m, n, nnzj) + return solver + end +end + +const R2NLS_allowed_subsolvers = (:cgls, :crls, :lsqr, :lsmr, :qrmumps) + +""" + R2NLS(nlp; kwargs...) + +An inexact second-order quadratic regularization method designed specifically for nonlinear least-squares problems: + + min ½‖F(x)‖² + +where `F: ℝⁿ → ℝᵐ` is a vector-valued function defining the least-squares residuals. + +For advanced usage, first create a `R2SolverNLS` to preallocate the necessary memory for the algorithm, and then call `solve!`: + + solver = R2SolverNLS(nlp) + solve!(solver, nlp; kwargs...) + +# Arguments +- `nlp::AbstractNLSModel{T, V}` is the nonlinear least-squares model to solve. See `NLPModels.jl` for additional details. +# Keyword Arguments +- `x::V = nlp.meta.x0`: the initial guess. +- `atol::T = √eps(T)`: absolute tolerance. +- `rtol::T = √eps(T)`: relative tolerance; the algorithm stops when ‖J(x)ᵀF(x)‖ ≤ atol + rtol * ‖J(x₀)ᵀF(x₀)‖. +- `params::R2NLSParameterSet = R2NLSParameterSet()`: algorithm parameters, see [`R2NLSParameterSet`](@ref). +- `η1::T = $(R2NLS_η1)`: step acceptance parameter, see [`R2NLSParameterSet`](@ref). +- `η2::T = $(R2NLS_η2)`: step acceptance parameter, see [`R2NLSParameterSet`](@ref). +- `θ1::T = $(R2NLS_θ1)`: Cauchy step parameter, see [`R2NLSParameterSet`](@ref). +- `θ2::T = $(R2NLS_θ2)`: Cauchy step parameter, see [`R2NLSParameterSet`](@ref). +- `γ1::T = $(R2NLS_γ1)`: regularization update parameter, see [`R2NLSParameterSet`](@ref). +- `γ2::T = $(R2NLS_γ2)`: regularization update parameter, see [`R2NLSParameterSet`](@ref). +- `γ3::T = $(R2NLS_γ3)`: regularization update parameter, see [`R2NLSParameterSet`](@ref). +- `δ1::T = $(R2NLS_δ1)`: Cauchy point calculation parameter, see [`R2NLSParameterSet`](@ref). +- `σmin::T = $(R2NLS_σmin)`: minimum step parameter, see [`R2NLSParameterSet`](@ref). +- `non_mono_size::Int = $(R2NLS_non_mono_size)`: the size of the non-monotone behaviour. If > 1, the algorithm will use a non-monotone strategy to accept steps. +- `scp_flag::Bool = true`: if true, compare the norm of the calculated step with `θ2 * norm(scp)` each iteration, selecting the smaller step. +- `subsolver::Symbol = :lsmr`: method used as subproblem solver, see `JSOSolvers.R2NLS_allowed_subsolvers` for a list. +- `subsolver_verbose::Int = 0`: if > 0, display subsolver iteration details every `subsolver_verbose` iterations when a KrylovWorkspace type is selected. +- `max_eval::Int = -1`: maximum number of objective function evaluations. +- `max_time::Float64 = 30.0`: maximum allowed time in seconds. +- `max_iter::Int = typemax(Int)`: maximum number of iterations. +- `verbose::Int = 0`: if > 0, displays iteration details every `verbose` iterations. + +# Output +Returns a `GenericExecutionStats` object containing statistics and information about the optimization process (see `SolverCore.jl`). + +# Callback +$(Callback_docstring) + +# Examples +```jldoctest +using JSOSolvers, ADNLPModels +F(x) = [x[1] - 1; 2 * (x[2] - x[1]^2)] +model = ADNLSModel(F, [-1.2; 1.0], 2) +stats = R2NLS(model) +# output +"Execution stats: first-order stationary" +``` +```jldoctest +using JSOSolvers, ADNLPModels +F(x) = [x[1] - 1; 2 * (x[2] - x[1]^2)] +model = ADNLSModel(F, [-1.2; 1.0], 2) +solver = R2SolverNLS(model) +stats = solve!(solver, model) +# output +"Execution stats: first-order stationary" +``` +""" +mutable struct R2SolverNLS{ + T, + V, + Op <: Union{AbstractLinearOperator{T}, SparseMatrixCOO{T, Int}}, + Sub <: Union{KrylovWorkspace{T, T, V}, QRMumpsSolver{T}}, +} <: AbstractOptimizationSolver + x::V + xt::V + temp::V + gx::V + Fx::V + rt::V + Jv::V + Jtv::V + Jx::Op + ls_subsolver::Sub + obj_vec::V # used for non-monotone behaviour + subtol::T + s::V + scp::V + σ::T + params::R2NLSParameterSet{T} +end + +function R2SolverNLS( + nlp::AbstractNLSModel{T, V}; + η1::T = get(R2NLS_η1, nlp), + η2::T = get(R2NLS_η2, nlp), + θ1::T = get(R2NLS_θ1, nlp), + θ2::T = get(R2NLS_θ2, nlp), + γ1::T = get(R2NLS_γ1, nlp), + γ2::T = get(R2NLS_γ2, nlp), + γ3::T = get(R2NLS_γ3, nlp), + δ1::T = get(R2NLS_δ1, nlp), + σmin::T = get(R2NLS_σmin, nlp), + non_mono_size::Int = get(R2NLS_non_mono_size, nlp), + subsolver::Symbol = :lsmr, +) where {T, V} + params = R2NLSParameterSet( + nlp; + η1 = η1, + η2 = η2, + θ1 = θ1, + θ2 = θ2, + γ1 = γ1, + γ2 = γ2, + γ3 = γ3, + δ1 = δ1, + σmin = σmin, + non_mono_size = non_mono_size, + ) + subsolver in R2NLS_allowed_subsolvers || + error("subproblem solver must be one of $(R2NLS_allowed_subsolvers)") + value(params.non_mono_size) >= 1 || error("non_mono_size must be greater than or equal to 1") + + nvar = nlp.meta.nvar + nequ = nlp.nls_meta.nequ + x = V(undef, nvar) + xt = V(undef, nvar) + temp = V(undef, nequ) + gx = V(undef, nvar) + Fx = V(undef, nequ) + rt = V(undef, nequ) + Jv = V(undef, nequ) + Jtv = V(undef, nvar) + s = V(undef, nvar) + scp = V(undef, nvar) + σ = eps(T)^(1 / 5) + if subsolver == :qrmumps + Jv = V(undef, 0) + Jtv = V(undef, 0) + ls_subsolver = QRMumpsSolver(nlp) + Jx = SparseMatrixCOO( + nequ, + nvar, + ls_subsolver.irn[1:ls_subsolver.nnzj], + ls_subsolver.jcn[1:ls_subsolver.nnzj], + ls_subsolver.val[1:ls_subsolver.nnzj], + ) + else + Jx = jac_op_residual!(nlp, x, Jv, Jtv) + ls_subsolver = krylov_workspace(Val(subsolver), nequ, nvar, V) + end + Sub = typeof(ls_subsolver) + Op = typeof(Jx) + + subtol = one(T) # must be ≤ 1.0 + obj_vec = fill(typemin(T), value(params.non_mono_size)) + + return R2SolverNLS{T, V, Op, Sub}( + x, + xt, + temp, + gx, + Fx, + rt, + Jv, + Jtv, + Jx, + ls_subsolver, + obj_vec, + subtol, + s, + scp, + σ, + params, + ) +end + +function SolverCore.reset!(solver::R2SolverNLS{T}) where {T} + fill!(solver.obj_vec, typemin(T)) + solver +end +function SolverCore.reset!(solver::R2SolverNLS{T}, nlp::AbstractNLSModel) where {T} + fill!(solver.obj_vec, typemin(T)) + solver +end + +@doc (@doc R2SolverNLS) function R2NLS( + nlp::AbstractNLSModel{T, V}; + η1::Real = get(R2NLS_η1, nlp), + η2::Real = get(R2NLS_η2, nlp), + θ1::Real = get(R2NLS_θ1, nlp), + θ2::Real = get(R2NLS_θ2, nlp), + γ1::Real = get(R2NLS_γ1, nlp), + γ2::Real = get(R2NLS_γ2, nlp), + γ3::Real = get(R2NLS_γ3, nlp), + δ1::Real = get(R2NLS_δ1, nlp), + σmin::Real = get(R2NLS_σmin, nlp), + non_mono_size::Int = get(R2NLS_non_mono_size, nlp), + subsolver::Symbol = :lsmr, + kwargs..., +) where {T, V} + solver = R2SolverNLS( + nlp; + η1 = convert(T, η1), + η2 = convert(T, η2), + θ1 = convert(T, θ1), + θ2 = convert(T, θ2), + γ1 = convert(T, γ1), + γ2 = convert(T, γ2), + γ3 = convert(T, γ3), + δ1 = convert(T, δ1), + σmin = convert(T, σmin), + non_mono_size = non_mono_size, + subsolver = subsolver, + ) + return solve!(solver, nlp; kwargs...) +end + +function SolverCore.solve!( + solver::R2SolverNLS{T, V}, + nlp::AbstractNLSModel{T, V}, + stats::GenericExecutionStats{T, V}; + callback = (args...) -> nothing, + x::V = nlp.meta.x0, + atol::T = √eps(T), + rtol::T = √eps(T), + Fatol::T = zero(T), + Frtol::T = zero(T), + max_time::Float64 = 30.0, + max_eval::Int = -1, + max_iter::Int = typemax(Int), + verbose::Int = 0, + scp_flag::Bool = true, + subsolver_verbose::Int = 0, +) where {T, V} + unconstrained(nlp) || error("R2NLS should only be called on unconstrained problems.") + if !(nlp.meta.minimize) + error("R2NLS only works for minimization problem") + end + + reset!(stats) + params = solver.params + η1 = value(params.η1) + η2 = value(params.η2) + θ1 = value(params.θ1) + θ2 = value(params.θ2) + γ1 = value(params.γ1) + γ2 = value(params.γ2) + γ3 = value(params.γ3) + δ1 = value(params.δ1) + σmin = value(params.σmin) + non_mono_size = value(params.non_mono_size) + + @assert(η1 > 0 && η1 < 1) + @assert(θ1 > 0 && θ1 < 1) + @assert(θ2 > 1) + @assert(γ1 >= 1 && γ1 <= γ2 && γ3 <= 1) + @assert(δ1>0 && δ1<1) + + start_time = time() + set_time!(stats, 0.0) + + n = nlp.nls_meta.nvar + m = nlp.nls_meta.nequ + x = solver.x .= x + xt = solver.xt + ∇f = solver.gx # k-1 + ls_subsolver = solver.ls_subsolver + r, rt = solver.Fx, solver.rt + s = solver.s + scp = solver.scp + subtol = solver.subtol + + σk = solver.σ + # preallocate storage for products with Jx and Jx' + Jx = solver.Jx + if Jx isa SparseMatrixCOO + jac_coord_residual!(nlp, x, view(ls_subsolver.val, 1:ls_subsolver.nnzj)) + Jx.vals .= view(ls_subsolver.val, 1:ls_subsolver.nnzj) + end + + residual!(nlp, x, r) + f = obj(nlp, x, r, recompute = false) + f0 = f + + mul!(∇f, Jx', r) + + norm_∇fk = norm(∇f) + ρk = zero(T) + + # Stopping criterion: + fmin = min(-one(T), f0) / eps(T) + unbounded = f < fmin + + σk = 2^round(log2(norm_∇fk + 1)) / norm_∇fk + ϵ = atol + rtol * norm_∇fk + ϵF = Fatol + Frtol * 2 * √f + + # Preallocate xt. + xt = solver.xt + temp = solver.temp + + optimal = norm_∇fk ≤ ϵ + small_residual = 2 * √f ≤ ϵF + + set_iter!(stats, 0) + set_objective!(stats, f) + set_dual_residual!(stats, norm_∇fk) + + if optimal + @info "Optimal point found at initial point" + @info log_header( + [:iter, :f, :dual, :σ, :ρ], + [Int, Float64, Float64, Float64, Float64], + hdr_override = Dict(:f => "f(x)", :dual => "‖∇f‖"), + ) + @info log_row([stats.iter, stats.objective, norm_∇fk, σk, ρk]) + end + cp_step_log = " " + if verbose > 0 && mod(stats.iter, verbose) == 0 + @info log_header( + [:iter, :f, :dual, :σ, :ρ, :sub_iter, :dir, :cp_step_log, :sub_status], + [Int, Float64, Float64, Float64, Float64, Int, String, String, String], + hdr_override = Dict( + :f => "f(x)", + :dual => "‖∇f‖", + :sub_iter => "subiter", + :dir => "dir", + :cp_step_log => "cp step", + :sub_status => "status", + ), + ) + @info log_row([stats.iter, stats.objective, norm_∇fk, σk, ρk, 0, " ", " ", " "]) + end + + set_status!( + stats, + get_status( + nlp, + elapsed_time = stats.elapsed_time, + optimal = optimal, + unbounded = unbounded, + max_eval = max_eval, + iter = stats.iter, + small_residual = small_residual, + max_iter = max_iter, + max_time = max_time, + ), + ) + + subtol = max(rtol, min(T(0.1), √norm_∇fk, T(0.9) * subtol)) + solver.σ = σk + solver.subtol = subtol + + callback(nlp, solver, stats) + + subtol = solver.subtol + σk = solver.σ + + done = stats.status != :unknown + ν_k = one(T) # used for scp calculation + + while !done + + # Compute the Cauchy step. + mul!(temp, Jx, ∇f) # temp <- Jx'*∇f + curv = dot(temp, temp) # curv = ∇f' Jx'Jx *∇f + slope = σk * norm_∇fk^2 # slope= σ * ||∇f||^2 + γ_k = (curv + slope) / norm_∇fk^2 + temp .= .-r + solver.σ = σk + + if γ_k > 0 + ν_k = 2*(1-δ1) / (γ_k) + cp_step_log = "α_k" + # Compute the step s. + subsolver_solved, sub_stats, subiter = + subsolve!(ls_subsolver, solver, nlp, s, atol, n, m, max_time, subsolver_verbose) + if scp_flag + # Based on the flag, scp is calcualted + scp .= -ν_k * ∇f + if norm(s) > θ2 * norm(scp) + s .= scp + end + end + else # when zero curvature occures + # we have to calcualte the scp, since we have encounter a negative curvature + λmax, found_λ = opnorm(Jx) + found_λ || error("operator norm computation failed") + cp_step_log = "ν_k" + ν_k = θ1 / (λmax + σk) + scp .= -ν_k * ∇f + s .= scp + end + + # Compute actual vs. predicted reduction. + xt .= x .+ s + mul!(temp, Jx, s) + slope = dot(r, temp) + curv = dot(temp, temp) + + residual!(nlp, xt, rt) + fck = obj(nlp, x, rt, recompute = false) + + ΔTk = -slope - curv / 2 + if non_mono_size > 1 #non-monotone behaviour + k = mod(stats.iter, non_mono_size) + 1 + solver.obj_vec[k] = stats.objective + fck_max = maximum(solver.obj_vec) + ρk = (fck_max - fck) / (fck_max - stats.objective + ΔTk) + else + ρk = (stats.objective - fck) / ΔTk + end + + # Update regularization parameters and Acceptance of the new candidate + step_accepted = ρk >= η1 + if step_accepted + if Jx isa SparseMatrixCOO # we need to update the values of Jx in QRMumpsSolver + jac_coord_residual!(nlp, x, view(ls_subsolver.val, 1:ls_subsolver.nnzj)) + Jx.vals .= view(ls_subsolver.val, 1:ls_subsolver.nnzj) + end + # update Jx implicitly for other solvers + x .= xt + r .= rt + f = fck + mul!(∇f, Jx', r) # ∇f = Jx' * r + set_objective!(stats, fck) + unbounded = fck < fmin + norm_∇fk = norm(∇f) + if ρk >= η2 + σk = max(σmin, γ3 * σk) + else # η1 ≤ ρk < η2 + σk = min(σmin, γ1 * σk) + end + else # η1 > ρk + σk = max(σmin, γ2 * σk) + end + + set_iter!(stats, stats.iter + 1) + set_time!(stats, time() - start_time) + + subtol = max(rtol, min(T(0.1), √norm_∇fk, T(0.9) * subtol)) + + solver.σ = σk + solver.subtol = subtol + set_dual_residual!(stats, norm_∇fk) + + callback(nlp, solver, stats) + + σk = solver.σ + subtol = solver.subtol + norm_∇fk = stats.dual_feas + + optimal = norm_∇fk ≤ ϵ + small_residual = 2 * √f ≤ ϵF + + if verbose > 0 && mod(stats.iter, verbose) == 0 + dir_stat = step_accepted ? "↘" : "↗" + @info log_row([ + stats.iter, + stats.objective, + norm_∇fk, + σk, + ρk, + subiter, + cp_step_log, + dir_stat, + sub_stats, + ]) + end + + if stats.status == :user + done = true + else + set_status!( + stats, + get_status( + nlp, + elapsed_time = stats.elapsed_time, + optimal = optimal, + unbounded = unbounded, + small_residual = small_residual, + max_eval = max_eval, + iter = stats.iter, + max_iter = max_iter, + max_time = max_time, + ), + ) + end + + done = stats.status != :unknown + end + + set_solution!(stats, x) + return stats +end + +# Dispatch for KrylovWorkspace +function subsolve!( + ls_subsolver::KrylovWorkspace, + R2NLS::R2SolverNLS, + nlp, + s, + atol, + n, + m, + max_time, + subsolver_verbose, +) + krylov_solve!( + ls_subsolver, + R2NLS.Jx, + R2NLS.temp, + atol = atol, + rtol = R2NLS.subtol, + λ = √(R2NLS.σ), # λ ≥ 0 is a regularization parameter. + itmax = max(2 * (n + m), 50), + # timemax = max_time - R2SolverNLS.stats.elapsed_time, + verbose = subsolver_verbose, + ) + s .= ls_subsolver.x + return Krylov.issolved(ls_subsolver), ls_subsolver.stats.status, ls_subsolver.stats.niter +end + +# Dispatch for QRMumpsSolver +function subsolve!( + ls::QRMumpsSolver, + R2NLS::R2SolverNLS, + nlp, + s, + atol, + n, + m, + max_time, + subsolver_verbose, +) + + # 1. Update Jacobian values at the current point x + # jac_coord_residual!(nlp, R2NLS.x, view(ls.val, 1:ls.nnzj)) + + # 2. Update regularization parameter σ + sqrt_σ = sqrt(R2NLS.σ) + @inbounds for i = 1:n + ls.val[ls.nnzj + i] = sqrt_σ + end + + # 3. Build the augmented right-hand side vector: b_aug = [-F(x); 0] + ls.b_aug[1:m] .= R2NLS.temp # -F(x) + fill!(view(ls.b_aug, (m + 1):(m + n)), zero(eltype(ls.b_aug))) # we have to do this for some reason #Applying all of its Householder (or Givens) transforms to the entire RHS vector b_aug—i.e. computing QTbQTb. + # Update spmat + qrm_update!(ls.spmat, ls.val) + + # 4. Solve the least-squares system + qrm_factorize!(ls.spmat, ls.spfct; transp = 'n') + qrm_apply!(ls.spfct, ls.b_aug; transp = 't') + qrm_solve!(ls.spfct, ls.b_aug, s; transp = 'n') + + # 5. Return status. For a direct solver, we assume success. + return true, "QRMumps", 1 +end diff --git a/test/allocs.jl b/test/allocs.jl index 70fbd214..82bad961 100644 --- a/test/allocs.jl +++ b/test/allocs.jl @@ -30,17 +30,27 @@ end if Sys.isunix() @testset "Allocation tests" begin - @testset "$symsolver" for symsolver in - (:LBFGSSolver, :FoSolver, :FomoSolver, :TrunkSolver, :TronSolver) + @testset "$name" for (name, symsolver) in ( + (:R2N, :R2NSolver), + (:R2N_exact, :R2NSolver), + (:R2, :FoSolver), + (:fomo, :FomoSolver), + (:lbfgs, :LBFGSSolver), + (:tron, :TronSolver), + (:trunk, :TrunkSolver), + ) for model in NLPModelsTest.nlp_problems nlp = eval(Meta.parse(model))() - if unconstrained(nlp) || (bound_constrained(nlp) && (symsolver == :TronSolver)) - if (symsolver == :FoSolver || symsolver == :FomoSolver) + if unconstrained(nlp) || (bound_constrained(nlp) && (name == :TronSolver)) + if (name == :FoSolver || name == :FomoSolver) solver = eval(symsolver)(nlp; M = 2) # nonmonotone configuration allocates extra memory + elseif name == :R2N_exact + solver = eval(symsolver)(LBFGSModel(nlp), subsolver= :shifted_lbfgs) + #TODO MA97 else solver = eval(symsolver)(nlp) end - if symsolver == :FomoSolver + if name == :FomoSolver T = eltype(nlp.meta.x0) stats = GenericExecutionStats(nlp, solver_specific = Dict(:avgβmax => T(0))) else @@ -48,8 +58,8 @@ if Sys.isunix() end with_logger(NullLogger()) do SolverCore.solve!(solver, nlp, stats) - SolverCore.reset!(solver) - NLPModels.reset!(nlp) + reset!(solver) + reset!(nlp) al = @wrappedallocs SolverCore.solve!(solver, nlp, stats) @test al == 0 end @@ -57,11 +67,21 @@ if Sys.isunix() end end - @testset "$symsolver" for symsolver in (:TrunkSolverNLS, :TronSolverNLS) + @testset "$name" for (name, symsolver) in ( + (:TrunkSolverNLS, :TrunkSolverNLS), + (:R2SolverNLS, :R2SolverNLS), + (:R2SolverNLS_QRMumps, :R2SolverNLS), + (:TronSolverNLS, :TronSolverNLS), + ) for model in NLPModelsTest.nls_problems nlp = eval(Meta.parse(model))() if unconstrained(nlp) || (bound_constrained(nlp) && (symsolver == :TronSolverNLS)) - solver = eval(symsolver)(nlp) + if name == :R2SolverNLS_QRMumps + solver = eval(symsolver)(nlp, subsolver = :qrmumps) + else + solver = eval(symsolver)(nlp) + end + stats = GenericExecutionStats(nlp) with_logger(NullLogger()) do SolverCore.solve!(solver, nlp, stats) diff --git a/test/callback.jl b/test/callback.jl index ddadc799..72234965 100644 --- a/test/callback.jl +++ b/test/callback.jl @@ -17,6 +17,11 @@ using ADNLPModels, JSOSolvers, LinearAlgebra, Logging #, Plots end @test stats.iter == 8 + stats = with_logger(NullLogger()) do + R2N(nlp, callback = cb) + end + @test stats.iter == 8 + stats = with_logger(NullLogger()) do lbfgs(nlp, callback = cb) end @@ -56,6 +61,11 @@ end tron(nls, callback = cb) end @test stats.iter == 8 + + stats = with_logger(NullLogger()) do + R2NLS(nls, callback = cb) + end + @test stats.iter == 8 end @testset "Testing Solver Values" begin diff --git a/test/consistency.jl b/test/consistency.jl index 8ec758fa..31e8b267 100644 --- a/test/consistency.jl +++ b/test/consistency.jl @@ -10,58 +10,86 @@ function consistency() @testset "Consistency" begin args = Pair{Symbol, Number}[:atol => 1e-6, :rtol => 1e-6, :max_eval => 20000, :max_time => 60.0] - @testset "NLP with $mtd" for mtd in [trunk, lbfgs, tron, R2, fomo] + @testset "NLP with $mtd" for (mtd, solver) in [ + ("trunk", trunk), + ("lbfgs", lbfgs), + ("tron", tron), + ("R2", R2), + ("R2N", R2N), + ( + "R2N_exact", + (nlp; kwargs...) -> + R2N(LBFGSModel(nlp), subsolver= :shifted_lbfgs; kwargs...), + ), + ("fomo", fomo), + ] with_logger(NullLogger()) do - NLPModels.reset!(unlp) - stats = mtd(unlp; args...) + reset!(unlp) + stats = solver(unlp; args...) @test stats isa GenericExecutionStats @test stats.status == :first_order - NLPModels.reset!(unlp) - stats = mtd(unlp; max_eval = 1) + reset!(unlp) + stats = solver(unlp; max_eval = 1) @test stats.status == :max_eval slow_nlp = ADNLPModel(x -> begin sleep(0.1) f(x) end, unlp.meta.x0) - stats = mtd(slow_nlp; max_time = 0.0) + stats = solver(slow_nlp; max_time = 0.0) @test stats.status == :max_time end end - @testset "Quasi-Newton NLP with $mtd" for mtd in [trunk, lbfgs, tron, R2, fomo] + @testset "Quasi-Newton NLP with $mtd" for (mtd, solver) in [ + ("trunk", trunk), + ("lbfgs", lbfgs), + ("tron", tron), + ("R2", R2), + ("R2N", R2N), + ( + "R2N_exact", + (nlp; kwargs...) -> + R2N(LBFGSModel(nlp), subsolver= :shifted_lbfgs; kwargs...), + ), + ("fomo", fomo), + ] with_logger(NullLogger()) do - NLPModels.reset!(qnlp) - stats = mtd(qnlp; args...) + reset!(qnlp) + stats = solver(qnlp; args...) @test stats isa GenericExecutionStats @test stats.status == :first_order end end - @testset "NLS with $mtd" for mtd in [trunk] + @testset "NLS with $mtd" for (mtd, solver) in [ + ("trunk", trunk), + ("R2NLS", (unls; kwargs...) -> R2NLS(unls; kwargs...)), + ("R2NLS_CGLS", (unls; kwargs...) -> R2NLS(unls, subsolver = :cgls; kwargs...)), + ("R2NLS_QRMumps", (unls; kwargs...) -> R2NLS(unls, subsolver = :qrmumps; kwargs...)), + ] with_logger(NullLogger()) do - stats = mtd(unls; args...) + stats = solver(unls; args...) @test stats isa GenericExecutionStats @test stats.status == :first_order NLPModels.reset!(unls) - stats = mtd(unls; max_eval = 1) + stats = solver(unls; max_eval = 1) @test stats.status == :max_eval slow_nls = ADNLSModel(x -> begin sleep(0.1) F(x) end, unls.meta.x0, nls_meta(unls).nequ) - stats = mtd(slow_nls; max_time = 0.0) + stats = solver(slow_nls; max_time = 0.0) @test stats.status == :max_time end end - @testset "Quasi-Newton NLS with $mtd" for mtd in [trunk] + @testset "Quasi-Newton NLS with $mtd" for (mtd, solver) in [("trunk", trunk)] with_logger(NullLogger()) do - stats = mtd(qnls; args...) + stats = solver(qnls; args...) @test stats isa GenericExecutionStats @test stats.status == :first_order end end end end - consistency() diff --git a/test/restart.jl b/test/restart.jl index 38765465..915f6719 100644 --- a/test/restart.jl +++ b/test/restart.jl @@ -1,4 +1,6 @@ @testset "Test restart with a different initial guess: $fun" for (fun, s) in ( + (:R2N, :R2NSolver), + (:R2N_exact, :R2NSolver), (:R2, :FoSolver), (:fomo, :FomoSolver), (:lbfgs, :LBFGSSolver), @@ -7,9 +9,14 @@ ) f(x) = (x[1] - 1)^2 + 4 * (x[2] - x[1]^2)^2 nlp = ADNLPModel(f, [-1.2; 1.0]) + if fun == :R2N_exact + nlp = LBFGSModel(nlp) + solver = eval(s)(nlp,subsolver= :shifted_lbfgs) + else + solver = eval(s)(nlp) + end stats = GenericExecutionStats(nlp) - solver = eval(s)(nlp) stats = SolverCore.solve!(solver, nlp, stats) @test stats.status == :first_order @test isapprox(stats.solution, [1.0; 1.0], atol = 1e-6) @@ -25,12 +32,30 @@ end @testset "Test restart NLS with a different initial guess: $fun" for (fun, s) in ( (:tron, :TronSolverNLS), (:trunk, :TrunkSolverNLS), + (:R2SolverNLS, :R2SolverNLS), + (:R2SolverNLS_CG, :R2SolverNLS), + (:R2SolverNLS_LSQR, :R2SolverNLS), + (:R2SolverNLS_CR, :R2SolverNLS), + (:R2SolverNLS_LSMR, :R2SolverNLS), + (:R2SolverNLS_QRMumps, :R2SolverNLS), ) F(x) = [x[1] - 1; 2 * (x[2] - x[1]^2)] nlp = ADNLSModel(F, [-1.2; 1.0], 2) stats = GenericExecutionStats(nlp) - solver = eval(s)(nlp) + if fun == :R2SolverNLS_CG + solver = eval(s)(nlp, subsolver = :cgls) + elseif fun == :R2SolverNLS_LSQR + solver = eval(s)(nlp, subsolver = :lsqr) + elseif fun == :R2SolverNLS_CR + solver = eval(s)(nlp, subsolver = :crls) + elseif fun == :R2SolverNLS_LSMR + solver = eval(s)(nlp, subsolver = :lsmr) + elseif fun == :R2SolverNLS_QRMumps + solver = eval(s)(nlp, subsolver = :qrmumps) + else + solver = eval(s)(nlp) + end stats = SolverCore.solve!(solver, nlp, stats) @test stats.status == :first_order @test isapprox(stats.solution, [1.0; 1.0], atol = 1e-6) @@ -44,6 +69,8 @@ end end @testset "Test restart with a different problem: $fun" for (fun, s) in ( + (:R2N, :R2NSolver), + (:R2N_exact, :R2NSolver), (:R2, :FoSolver), (:fomo, :FomoSolver), (:lbfgs, :LBFGSSolver), @@ -52,15 +79,25 @@ end ) f(x) = (x[1] - 1)^2 + 4 * (x[2] - x[1]^2)^2 nlp = ADNLPModel(f, [-1.2; 1.0]) + if fun == :R2N_exact + nlp = LBFGSModel(nlp) + solver = eval(s)(nlp,subsolver= :shifted_lbfgs) + else + solver = eval(s)(nlp) + end stats = GenericExecutionStats(nlp) - solver = eval(s)(nlp) stats = SolverCore.solve!(solver, nlp, stats) @test stats.status == :first_order @test isapprox(stats.solution, [1.0; 1.0], atol = 1e-6) f2(x) = (x[1])^2 + 4 * (x[2] - x[1]^2)^2 nlp = ADNLPModel(f2, [-1.2; 1.0]) + if fun == :R2N_exact + nlp = LBFGSModel(nlp) + else + solver = eval(s)(nlp) + end SolverCore.reset!(solver, nlp) stats = SolverCore.solve!(solver, nlp, stats, atol = 1e-10, rtol = 1e-10) @@ -68,15 +105,34 @@ end @test isapprox(stats.solution, [0.0; 0.0], atol = 1e-6) end + @testset "Test restart NLS with a different problem: $fun" for (fun, s) in ( (:tron, :TronSolverNLS), (:trunk, :TrunkSolverNLS), + (:R2SolverNLS, :R2SolverNLS), + (:R2SolverNLS_CG, :R2SolverNLS), + (:R2SolverNLS_LSQR, :R2SolverNLS), + (:R2SolverNLS_CR, :R2SolverNLS), + (:R2SolverNLS_LSMR, :R2SolverNLS), + (:R2SolverNLS_QRMumps, :R2SolverNLS), ) F(x) = [x[1] - 1; 2 * (x[2] - x[1]^2)] nlp = ADNLSModel(F, [-1.2; 1.0], 2) stats = GenericExecutionStats(nlp) - solver = eval(s)(nlp) + if fun == :R2SolverNLS_CG + solver = eval(s)(nlp, subsolver = :cgls) + elseif fun == :R2SolverNLS_LSQR + solver = eval(s)(nlp, subsolver = :lsqr) + elseif fun == :R2SolverNLS_CR + solver = eval(s)(nlp, subsolver = :crls) + elseif fun == :R2SolverNLS_LSMR + solver = eval(s)(nlp, subsolver = :lsmr) + elseif fun == :R2SolverNLS_QRMumps + solver = eval(s)(nlp, subsolver = :qrmumps) + else + solver = eval(s)(nlp) + end stats = SolverCore.solve!(solver, nlp, stats) @test stats.status == :first_order @test isapprox(stats.solution, [1.0; 1.0], atol = 1e-6) diff --git a/test/runtests.jl b/test/runtests.jl index 22461248..07f49dcf 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,7 +2,7 @@ using Printf, LinearAlgebra, Logging, SparseArrays, Test # additional packages -using ADNLPModels, LinearOperators, NLPModels, NLPModelsModifiers, SolverCore, SolverTools +using ADNLPModels, Krylov, LinearOperators, NLPModels, NLPModelsModifiers, SolverCore, SolverTools using NLPModelsTest, SolverParameters # this package @@ -14,6 +14,7 @@ using JSOSolvers (TRONParameterSet, tron), (TRUNKParameterSet, trunk), (FOMOParameterSet, fomo), + (R2NParameterSet, R2N), ) nlp = BROWNDEN() params = eval(paramset)(nlp) @@ -48,7 +49,7 @@ end end @testset "Test iteration limit" begin - @testset "$fun" for fun in (R2, fomo, lbfgs, tron, trunk) + @testset "$fun" for fun in (R2, fomo, lbfgs, tron, trunk, R2N) f(x) = (x[1] - 1)^2 + 4 * (x[2] - x[1]^2)^2 nlp = ADNLPModel(f, [-1.2; 1.0]) @@ -56,7 +57,7 @@ end @test stats.status == :max_iter end - @testset "$(fun)-NLS" for fun in (tron, trunk) + @testset "$(fun)-NLS" for fun in (tron, trunk, R2NLS) f(x) = [x[1] - 1; 2 * (x[2] - x[1]^2)] nlp = ADNLSModel(f, [-1.2; 1.0], 2) @@ -66,13 +67,24 @@ end end @testset "Test unbounded below" begin - @testset "$fun" for fun in (R2, fomo, lbfgs, tron, trunk) + @testset "$name" for (name, solver) in [ + ("trunk", trunk), + ("lbfgs", lbfgs), + ("tron", tron), + ("R2", R2), + # ("R2N", R2N), + ( + "R2N_exact", + (nlp; kwargs...) -> + R2N(LBFGSModel(nlp), subsolver= :shifted_lbfgs; kwargs...), + ), + ("fomo", fomo), + ] T = Float64 x0 = [T(0)] f(x) = -exp(x[1]) nlp = ADNLPModel(f, x0) - - stats = eval(fun)(nlp) + stats = solver(nlp) @test stats.status == :unbounded @test stats.objective < -one(T) / eps(T) end diff --git a/test/test_solvers.jl b/test/test_solvers.jl index 30216e4e..9cd913b2 100644 --- a/test/test_solvers.jl +++ b/test/test_solvers.jl @@ -8,6 +8,8 @@ function tests() ("lbfgs", lbfgs), ("tron", tron), ("R2", R2), + ("R2N", R2N), + ("R2N_exact", (nlp; kwargs...) -> R2N(LBFGSModel(nlp), subsolver= :shifted_lbfgs; kwargs...)), ("fomo_r2", fomo), ("fomo_tr", (nlp; kwargs...) -> fomo(nlp, step_backend = JSOSolvers.tr_step(); kwargs...)), ] @@ -41,6 +43,11 @@ function tests() ("trunk full Hessian", (nls; kwargs...) -> trunk(nls, variant = :Newton; kwargs...)), ("tron+cgls", (nls; kwargs...) -> tron(nls, subsolver = :cgls; kwargs...)), ("tron full Hessian", (nls; kwargs...) -> tron(nls, variant = :Newton; kwargs...)), + ("R2NLS", (unls; kwargs...) -> R2NLS(unls; kwargs...)), + ("R2NLS_CGLS", (unls; kwargs...) -> R2NLS(unls, subsolver = :cgls; kwargs...)), + ("R2NLS_LSQR", (unls; kwargs...) -> R2NLS(unls, subsolver = :lsqr; kwargs...)), + ("R2NLS_CRLS", (unls; kwargs...) -> R2NLS(unls, subsolver = :crls; kwargs...)), + ("R2NLS_LSMR", (unls; kwargs...) -> R2NLS(unls, subsolver = :lsmr; kwargs...)), ] unconstrained_nls(solver) multiprecision_nls(solver, :unc)