Add R2N and R2NLS solvers with HSL and QRMumps support #347

farhadrclass · 2026-01-13T12:46:24Z

Introduces new second-order quadratic regularization solvers R2N and R2NLS for unconstrained and nonlinear least-squares optimization. Adds support for HSL (MA97, MA57) and QRMumps direct solvers, updates documentation and README, and extends test coverage for the new solvers. Updates dependencies and compat entries in Project.toml.

src/R2N.jl

d-monnet · 2026-01-21T17:05:42Z

src/R2N.jl

+
+    if !subsolver_solved && npcCount == 0
+      @warn("Subsolver failed to solve the system")
+      # TODO exit cleaning, update stats


there should be a (existing?) flag to set stats.status with for this early stop.

yes but I think we need a new one and we can decide on it Tuesday?

d-monnet · 2026-01-21T17:06:56Z

src/R2N.jl

+
+      # Check Indefinite (Negative Eigenvalues)
+      if !force_sigma_increase && num_neg > 0
+        # Curvature = s' (H+σI) s


inconsistent comment

d-monnet · 2026-01-21T17:09:53Z

src/R2N.jl

+            # A: Armijo (Upper Bound) - Is step too big?
+            armijo_pass = fck <= f0_val + c_goldstein * α * grad_dir
+
+            # B: Curvature (Lower Bound) - Is step too small?


this is not a condition on the curvature, this is the Goldstein condition. Update comments and variable names accordingly.

src/R2NLS.jl

src/R2N.jl

src/R2NLS.jl

d-monnet · 2026-01-26T19:51:15Z

src/R2N.jl

+        α_min = zero(T)      # Lower bound of valid interval
+        α_max = T(Inf)       # Upper bound of valid interval
+        iter_ls = 0
+        max_iter_ls = 100    # Safety break #TODO Prof Orban?


This should be a kwarg

d-monnet · 2026-01-26T19:58:29Z

src/R2N.jl

+        ρ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)


Delta Tk theoretically cannot be negative. An error should be thrown in that case, see how this case is handled in the other solver.
Keep rho_k as usual.

d-monnet · 2026-01-26T19:58:56Z

src/R2N.jl

+        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   


denominator should be Delta Tk.

d-monnet · 2026-01-26T20:05:05Z

src/R2N.jl

+      if npc_handler == :gs # Goldstein Line Search
+        npcCount = 0
+
+        # --- Goldstein Line Search --- 


put the line search algo in a separate function.

farhadrclass · 2026-01-27T16:05:29Z

src/R2N.jl

+        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   


if very first iteration and subsolver very first iteration return npc, and we selected npc_handler = prev, then we need a better way, maybe use a flag to tell sigma to increase and try agian

farhadrclass · 2026-01-27T16:06:18Z

src/R2N.jl

+        ρ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)


this is the issue , iteration should be failed then move to next

Co-authored-by: Dominique <[email protected]>

dpo

Lots more comments on R2NLS.

src/R2NLS.jl

Co-authored-by: Dominique <[email protected]>

Rename model argument from `nlp` to `nls` throughout and update DefaultParameter lambdas to accept the new name. Rename solver residual field `Fx` to `r`. Improve QRMumps solver initialization and memory handling: use nls_meta, use views for COO arrays/values, zero the augmented RHS with an in-place loop, add a user-facing destroy! wrapper and keep a finalizer. Use views when constructing SparseMatrixCOO from QRMumps buffers and update jacobian/value update calls to the new `nls` API. Introduce `scp_est` option and refactor Cauchy/step safeguard logic (compute Cauchy step via curvature or operator norm), change initial σk heuristic to use Jacobian scale, and correct predicted/actual reduction and non-monotone bookkeeping to use residual/objective consistently. Tidy up logging columns and remove stale cp_step logging. Fix several small bugs in σ update rules and function argument names; general cleanup and minor performance improvements (views, fewer allocations).

Bump LinearOperators compat to 2.12.0 and refactor R2N internals. - Remove local ShiftedOperator implementation (use library/operator-level solutions). - Add LineModel field (h) and model type parameter to R2NSolver, initialize it in the constructor and recreate it when resetting with a new NLP. - Replace large custom Goldstein line-search block with SolverTools.redirect! + armijo_goldstein; compute step (s) and objective (ft) via that helper. - Replace several uses of fck with ft (clarify naming), adjust ρk and σ update logic, and ensure H updates and objective/stats use the computed ft. - Minor logging/formatting, whitespace, and small correctness tweaks (vectorized assignments, return values in reset!, etc.). These changes centralize line-search behavior, simplify solver code, and prepare for compatibility with newer LinearOperators behavior.

dpo

I propose a cleaner subsolver API.

src/R2NLS.jl

dpo · 2026-02-10T13:27:44Z

src/R2NLS.jl

+
+    function destroy!(s::QRMumpsSolver) #for user use, in case they want to free memory before the finalizer runs
+      free_qrm(s)
+    end


I don't think it is possible to call this function since it is declared inside the constructor. At any rate, I don't see what its purpose would be.

src/R2NLS.jl

Co-authored-by: Dominique <[email protected]>

Add an AbstractR2NLSSubsolver interface and unify subproblem solvers behind initialize!, update_subsolver!, solve_subproblem! and get_jacobian. Implement QRMumpsSubsolver (replacing the old QRMumpsSolver) and a GenericKrylovSubsolver with convenience constructors LSMRSubsolver, LSQRSubsolver and CGLSSubsolver. Extend R2NLSParameterSet with compute_cauchy_point and inexact_cauchy_point flags (and validation) and expose new subsolver exports. Refactor R2NLSSolver to accept either a subsolver Type or a preallocated AbstractR2NLSSubsolver instance, update internal fields, instantiate/update subsolvers, and drive the main solve loop through the new subsolver API. Update the trust-region / Levenberg–Marquardt loop to compute the Cauchy point conditionally, update σ initialization heuristic, streamline acceptance logic, and adjust logging field names. Remove legacy QRMumpsSolver and old subsolve dispatches; add resource cleanup (finalizer) for QRMumpsSubsolver.

farhadrclass · 2026-02-11T00:17:24Z

@dpo I have refactor the R2NLS,

farhadrclass commented Jan 20, 2026

View reviewed changes