LMModel → LLSModels and R2NModel → QuadraticModels

Project.toml

@@ -4,16 +4,20 @@ author = ["Robert Baraldi <[email protected]> and Dominique Orban <dominique.orban
 version = "0.1.0"
 
 [deps]
+ADNLPModels = "54578032-b7ea-4c30-94aa-7cbd1cce6c9a"
 Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
+LLSModels = "39f5bc3e-5160-4bf8-ac48-504fd2534d24"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearOperators = "5c8ed15e-5a4c-59e4-a42b-c7e8811fb125"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 ManualNLPModels = "30dfa513-9b2f-4fb3-9796-781eabac1617"
 NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6"
 NLPModelsModifiers = "e01155f1-5c6f-4375-a9d8-616dd036575f"
+OptimizationProblems = "5049e819-d29b-5fba-b941-0eee7e64c1c6"
 Percival = "01435c0c-c90d-11e9-3788-63660f8fbccc"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 ProximalOperators = "a725b495-10eb-56fe-b38b-717eba820537"
+QuadraticModels = "f468eda6-eac5-11e8-05a5-ff9e497bcd19"
 RegularizedProblems = "ea076b23-609f-44d2-bb12-a4ae45328278"
 ShiftedProximalOperators = "d4fd37fa-580c-4e43-9b30-361c21aae263"
 SolverCore = "ff4d7338-4cf1-434d-91df-b86cb86fb843"

src/LMTR_alg.jl

@@ -41,8 +41,9 @@ function LMTRSolver(
   xk = similar(x0)
   ∇fk = similar(x0)
   mν∇fk = similar(x0)
-  Fk = similar(x0, reg_nls.model.nls_meta.nequ)
-  Fkn = similar(Fk)
+  # residuals will be allocated after creating Jacobian operator
+  Fk = similar(x0, 0)
+  Fkn = similar(x0, 0)
   xkn = similar(x0)
   s = similar(x0)
   has_bnds = any(l_bound .!= T(-Inf)) || any(u_bound .!= T(Inf)) || subsolver == TRDHSolver
@@ -67,7 +68,24 @@ function LMTRSolver(
     ) : shifted(reg_nls.h, xk, one(T), χ)
 
   Jk = jac_op_residual(reg_nls.model, xk)
-  sub_nlp = LMModel(Jk, Fk, T(0), xk)
+  m = size(Jk, 1)
+  Fk = zeros(T, m)
+  Fkn = similar(Fk)
+  residual!(reg_nls.model, xk, Fk)
+
+  # Matrix-free Gram operator for J'J
+  gram_op = JacobianGram{T}(Jk, zeros(T, m))
+  # Or, if JacobianGram is not available, use LinearOperator:
+  # gram_op = LinearOperator{T}(size(Jk,2), size(Jk,2), false, false,
+  #   (y, x) -> mul!(y, Jk', Jk * x),
+  #   (y, x) -> mul!(y, Jk, Jk' * x),
+  #   (y, x) -> mul!(y, Jk, Jk' * x),
+  # )
+
+  JtF = similar(xk)
+  mul!(JtF, Jk', Fk)
+  # Don't include the constant term to avoid numerical overflow in obj evaluation
+  sub_nlp = QuadraticModel(JtF, gram_op, c0 = zero(T), x0 = zeros(T, length(xk)), name = "LMTR-subproblem")
   subpb = RegularizedNLPModel(sub_nlp, ψ)
   substats = RegularizedExecutionStats(subpb)
   subsolver = subsolver(subpb)
@@ -273,8 +291,20 @@ function SolverCore.solve!(
   jtprod_residual!(nls, xk, Fk, ∇fk)
   fk = dot(Fk, Fk) / 2
 
-  σmax, found_σ = opnorm(solver.subpb.model.J)
+  # Get Jacobian operator from original NLS model since QuadraticModel doesn't store it
+  Jk = jac_op_residual(nls, xk)
+  σmax, found_σ = opnorm(Jk)
   found_σ || error("operator norm computation failed")
+
+  # Update the QuadraticModel with current Jacobian and gradient
+  JtF = Jk' * Fk
+  JtJ = Jk' * Jk  
+  # In-place update of QuadraticModel fields to avoid allocation
+  solver.subpb.model.data.c = JtF
+  solver.subpb.model.data.H = JtJ
+  solver.subpb.model.data.c0 = zero(T)
+  solver.subpb.model.x0 .= 0
+  solver.subpb.model.name = "LMTR-subproblem"
   ν = α * Δk / (1 + σmax^2 * (α * Δk + 1))
   @. mν∇fk = -∇fk * ν
   sqrt_ξ1_νInv = one(T)
@@ -416,8 +446,15 @@ function SolverCore.solve!(
       shift!(ψ, xk)
       jtprod_residual!(nls, xk, Fk, ∇fk)
 
-      σmax, found_σ = opnorm(solver.subpb.model.J)
+      # Get Jacobian operator from original NLS model
+      Jk_update = jac_op_residual(nls, xk)
+      σmax, found_σ = opnorm(Jk_update)
       found_σ || error("operator norm computation failed")
+
+      # Update the QuadraticModel with new Jacobian and gradient
+      JtF = Jk_update' * Fk
+      JtJ = Jk_update' * Jk_update  
+      update_quadratic_model!(solver.subpb.model, JtF, JtJ)
     end
 
     if η2 ≤ ρk < Inf