R2NLS solver for unconstrained nonlinear least-squares (NLS) problems

Project.toml

@@ -3,25 +3,35 @@ 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"
 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]

README.md

@@ -26,14 +26,17 @@ 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).
     > Worst-case evaluation complexity for unconstrained nonlinear optimization using
     > 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

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
 ```

myR2N/Project.toml

@@ -0,0 +1,2 @@
+[deps]
+HSL_jll = "017b0a0e-03f4-516a-9b91-836bbd1904dd"

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"

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

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