From b65336d1f359e8305d8e85919d57e451650415a9 Mon Sep 17 00:00:00 2001 From: tmigot Date: Thu, 28 Jul 2022 22:35:58 +0200 Subject: [PATCH 1/5] add cglsshifts and lsqrshifts --- src/PreProcess/PreProcessLSKARC.jl | 47 ++++++++++++++++++ src/SolveModel/SolveModelLSKARC.jl | 23 +++++++++ src/solvers.jl | 26 +++++----- src/utils/pdata_struct.jl | 76 ++++++++++++++++++++++++++++++ 4 files changed, 160 insertions(+), 12 deletions(-) create mode 100644 src/PreProcess/PreProcessLSKARC.jl create mode 100644 src/SolveModel/SolveModelLSKARC.jl diff --git a/src/PreProcess/PreProcessLSKARC.jl b/src/PreProcess/PreProcessLSKARC.jl new file mode 100644 index 00000000..7f4d0fa5 --- /dev/null +++ b/src/PreProcess/PreProcessLSKARC.jl @@ -0,0 +1,47 @@ +function preprocess(PData::PDataLSKARC, Jx, Fx, gNorm2, calls, max_calls, α) + ζ, ξ, maxtol, mintol = PData.ζ, PData.ξ, PData.maxtol, PData.mintol + nshifts = PData.nshifts + shifts = PData.shifts + + m, n = size(Jx) + # Tolerance used in Assumption 2.6b in the paper ( ξ > 0, 0 < ζ ≤ 1 ) + atol = PData.cgatol(ζ, ξ, maxtol, mintol, gNorm2) + rtol = PData.cgrtol(ζ, ξ, maxtol, mintol, gNorm2) + + nshifts = length(shifts) + cb = (slv) -> begin + ind = setdiff(1:length(shifts), findall(.!slv.converged)) + if length(ind) > 1 + for i in ind + if (norm(slv.x[i]) / shifts[i] - α > 0) + return true + end + end + end + return false + end + solver = PData.solver + solve!( + solver, + Jx, + Fx, + shifts, + itmax = min(max_calls - sum(calls), 2 * (m + n)), + atol = atol, + rtol = rtol, + verbose = 0, + callback = cb, + ) + + PData.indmin = 0 + PData.positives .= solver.converged + for i = 1:nshifts + @. PData.xShift[i] = -solver.x[i] + PData.norm_dirs[i] = norm(solver.x[i]) # Get it from cg_lanczos? + end + PData.shifts .= shifts + PData.nshifts = nshifts + PData.OK = sum(solver.converged) != 0 # at least one system was solved + + return PData +end diff --git a/src/SolveModel/SolveModelLSKARC.jl b/src/SolveModel/SolveModelLSKARC.jl new file mode 100644 index 00000000..9ac94f2f --- /dev/null +++ b/src/SolveModel/SolveModelLSKARC.jl @@ -0,0 +1,23 @@ +function solve_modelLSKARC(X::PDataLSKARC, Jx, Fx, gNorm2, calls, max_calls, α::T) where {T} + # target value should be close to satisfy αλ=||d|| + start = findfirst(X.positives) + if isnothing(start) + start = 1 + end + if VERSION < v"1.7.0" + positives = collect(start:length(X.positives)) + target = [(abs(α * X.shifts[i] - X.norm_dirs[i])) for i in positives] + else + positives = start:length(X.positives) + target = ((abs(α * X.shifts[i] - X.norm_dirs[i])) for i in positives) + end + + # pick the closest shift to the target within positive definite H+λI + indmin = argmin(target) + X.indmin = start + indmin - 1 + p_imin = X.indmin + X.d .= X.xShift[p_imin] + X.λ = X.shifts[p_imin] + + return X.d, X.λ +end diff --git a/src/solvers.jl b/src/solvers.jl index dbd4c6dd..43489d2c 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -16,16 +16,18 @@ const solvers_const = Dict( ) const solvers_nls_const = Dict( - :ARCqKOpGN => ( - HessGaussNewtonOp, - PDataKARC, - solve_modelKARC, - [:shifts => 10.0 .^ (collect(-10.0:0.5:20.0))], - ), - :ST_TROpGN => (HessGaussNewtonOp, PDataST, solve_modelST_TR, ()), - :ST_TROpGNLSCgls => - (HessGaussNewtonOp, PDataNLSST, solve_modelNLSST_TR, [:solver_method => :cgls]), - :ST_TROpGNLSLsqr => - (HessGaussNewtonOp, PDataNLSST, solve_modelNLSST_TR, [:solver_method => :lsqr]), - :ST_TROpLS => (HessOp, PDataNLSST, solve_modelNLSST_TR, ()), + :ARCqKOpGN => ( + HessGaussNewtonOp, + PDataKARC, + solve_modelKARC, + [:shifts => 10.0 .^ (collect(-10.0:0.5:20.0))], + ), + :ST_TROpGN => (HessGaussNewtonOp, PDataST, solve_modelST_TR, ()), + :ST_TROpGNLSCgls => + (HessGaussNewtonOp, PDataNLSST, solve_modelNLSST_TR, [:solver_method => :cgls]), + :ST_TROpGNLSLsqr => + (HessGaussNewtonOp, PDataNLSST, solve_modelNLSST_TR, [:solver_method => :lsqr]), + :ST_TROpLS => (HessOp, PDataNLSST, solve_modelNLSST_TR, ()), + :LSARCqKOpCgls => (HessGaussNewtonOp, PDataLSKARC, solve_modelLSKARC, [:shifts => 10.0 .^ (collect(-10.0:0.5:20.0)), :solver_method => :cgls]), + :LSARCqKOpLsqr => (HessGaussNewtonOp, PDataLSKARC, solve_modelLSKARC, [:shifts => 10.0 .^ (collect(-10.0:0.5:20.0)), :solver_method => :lsqr]), ) diff --git a/src/utils/pdata_struct.jl b/src/utils/pdata_struct.jl index 6964c86f..453dc0ab 100644 --- a/src/utils/pdata_struct.jl +++ b/src/utils/pdata_struct.jl @@ -115,6 +115,82 @@ function PDataKARC( ) end +""" + PDataLSKARC(::Type{S}, ::Type{T}, n) +Return a structure used for the preprocessing of ARCqK methods. +""" +mutable struct PDataLSKARC{T} <: PDataIter{T} + d::Array{T,1} # (H+λI)\g ; on first call = g + λ::T # "active" value of λ; on first call = 0 + ζ::T # Inexact Newton order parameter: stop when ||∇q|| < ξ * ||g||^(1+ζ) + ξ::T # Inexact Newton order parameter: stop when ||∇q|| < ξ * ||g||^(1+ζ) + maxtol::T # Largest tolerance for Inexact Newton + mintol::T # Smallest tolerance for Inexact Newton + cgatol::Any + cgrtol::Any + + indmin::Int # index of best shift value within "positive". On first call = 0 + + positives::Array{Bool,1} # indices of the shift values yielding (H+λI)⪰0 + xShift::Array{Array{T,1},1} # solutions for each shifted system + shifts::Array{T,1} # values of the shifts + nshifts::Int # number of shifts + norm_dirs::Array{T,1} # norms of xShifts + OK::Bool # preprocess success + + solver::Union{CglsLanczosShiftSolver,LsqrShiftSolver} +end + +function PDataLSKARC( + ::Type{S}, + ::Type{T}, + n; + ζ = T(0.5), + ξ = T(0.01), + maxtol = T(0.01), + mintol = sqrt(eps(T)), + cgatol = (ζ, ξ, maxtol, mintol, gNorm2) -> max(mintol, min(maxtol, ξ * gNorm2^(1 + ζ))), + cgrtol = (ζ, ξ, maxtol, mintol, gNorm2) -> max(mintol, min(maxtol, ξ * gNorm2^ζ)), + shifts = 10.0 .^ collect(-20.0:1.0:20.0), + solver_method = :cgls, + kwargs..., +) where {S,T} + d = S(undef, n) + λ = zero(T) + indmin = 1 + nshifts = length(shifts) + positives = Array{Bool,1}(undef, nshifts) + xShift = Array{S,1}(undef, nshifts) + for i = 1:nshifts + xShift[i] = S(undef, n) + end + norm_dirs = S(undef, nshifts) + OK = true + solver = if solver_method == :cgls + CglsLanczosShiftSolver(n, n, nshifts, S) + else + LsqrShiftSolver(n, n, nshifts, S) + end + return PDataLSKARC( + d, + λ, + ζ, + ξ, + maxtol, + mintol, + cgatol, + cgrtol, + indmin, + positives, + xShift, + T.(shifts), + nshifts, + norm_dirs, + OK, + solver, + ) +end + """ PDataTRK(::Type{S}, ::Type{T}, n) Return a structure used for the preprocessing of TRK methods. From 657fdb790dec2cb8c250fc8a3657e47a10ab233e Mon Sep 17 00:00:00 2001 From: tmigot Date: Thu, 28 Jul 2022 22:50:18 +0200 Subject: [PATCH 2/5] update comments --- src/PreProcess/PreProcessLSKARC.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/PreProcess/PreProcessLSKARC.jl b/src/PreProcess/PreProcessLSKARC.jl index 7f4d0fa5..c9f337ff 100644 --- a/src/PreProcess/PreProcessLSKARC.jl +++ b/src/PreProcess/PreProcessLSKARC.jl @@ -13,7 +13,7 @@ function preprocess(PData::PDataLSKARC, Jx, Fx, gNorm2, calls, max_calls, α) ind = setdiff(1:length(shifts), findall(.!slv.converged)) if length(ind) > 1 for i in ind - if (norm(slv.x[i]) / shifts[i] - α > 0) + if (norm(slv.x[i]) / shifts[i] - α > 0) # for lsqr: sqrt(slv.xNorm²[i]) return true end end @@ -37,7 +37,7 @@ function preprocess(PData::PDataLSKARC, Jx, Fx, gNorm2, calls, max_calls, α) PData.positives .= solver.converged for i = 1:nshifts @. PData.xShift[i] = -solver.x[i] - PData.norm_dirs[i] = norm(solver.x[i]) # Get it from cg_lanczos? + PData.norm_dirs[i] = norm(solver.x[i]) end PData.shifts .= shifts PData.nshifts = nshifts From 3094afa32f1952ca4141bcc0a80a241e1b8aa145 Mon Sep 17 00:00:00 2001 From: tmigot Date: Fri, 29 Jul 2022 18:00:41 +0200 Subject: [PATCH 3/5] try fix using pdata with shifts for least-squares (segfault locally) --- src/main.jl | 45 ++++++++++++++++++++++++++------------- src/utils/pdata_struct.jl | 9 ++++---- 2 files changed, 35 insertions(+), 19 deletions(-) diff --git a/src/main.jl b/src/main.jl index 1a3e85b5..a68b0d70 100644 --- a/src/main.jl +++ b/src/main.jl @@ -83,6 +83,21 @@ function preprocess(stp::NLPStopping, PData::TPData, workspace::TRARCWorkspace, return PData end +function preprocess( + stp::NLPStopping, + PData::PDataIterLS, + workspace::TRARCWorkspace{T,S,Hess}, + ∇f, + norm_∇f, + α, +) where {T,S,Hess<:HessGaussNewtonOp} + max_hprod = stp.meta.max_cntrs[:neval_jprod_residual] + Fx = workspace.Fx + Jx = jac_op_residual!(stp.pb, workspace.xt, workspace.Hstruct.Jv, workspace.Hstruct.Jtv) + PData = ARCTR.preprocess(PData, Jx, Fx, norm_∇f, neval_jprod_residual(stp.pb), max_hprod, α) + return PData +end + function compute_direction( stp::NLPStopping, PData::TPData, @@ -132,21 +147,21 @@ function hessian!(workspace::TRARCWorkspace, nlp, x) end function TRARC( - nlp_stop::NLPStopping{Pb, M, SRC, NLPAtX{Score, T, S}, MStp, LoS}; - TR::TrustRegion = TrustRegion(T(10.0)), - hess_type::Type{Hess} = HessOp, - pdata_type::Type{ParamData} = PDataKARC, - kwargs..., -) where {Pb, M, SRC, MStp, LoS, Score, S, T, Hess, ParamData} - nlp = nlp_stop.pb - - if ParamData == PDataNLSST - PData = PDataNLSST(S, T, nlp.meta.nvar, nlp.nls_meta.nequ; kwargs...) - else - PData = ParamData(S, T, nlp.meta.nvar; kwargs...) - end - workspace = TRARCWorkspace(nlp, Hess) - return TRARC(nlp_stop, PData, workspace, TR; kwargs...) + nlp_stop::NLPStopping{Pb,M,SRC,NLPAtX{Score,T,S},MStp,LoS}; + TR::TrustRegion = TrustRegion(T(10.0)), + hess_type::Type{Hess} = HessOp, + pdata_type::Type{ParamData} = PDataKARC, + kwargs..., +) where {Pb,M,SRC,MStp,LoS,Score,S,T,Hess,ParamData} + nlp = nlp_stop.pb + + if ParamData in (PDataNLSST, PDataLSKARC) + PData = ParamData(S, T, nlp.meta.nvar, nlp.nls_meta.nequ; kwargs...) + else + PData = ParamData(S, T, nlp.meta.nvar; kwargs...) + end + workspace = TRARCWorkspace(nlp, Hess) + return TRARC(nlp_stop, PData, workspace, TR; kwargs...) end """ diff --git a/src/utils/pdata_struct.jl b/src/utils/pdata_struct.jl index 453dc0ab..49139a15 100644 --- a/src/utils/pdata_struct.jl +++ b/src/utils/pdata_struct.jl @@ -119,7 +119,7 @@ end PDataLSKARC(::Type{S}, ::Type{T}, n) Return a structure used for the preprocessing of ARCqK methods. """ -mutable struct PDataLSKARC{T} <: PDataIter{T} +mutable struct PDataLSKARC{T} <: PDataIterLS{T} d::Array{T,1} # (H+λI)\g ; on first call = g λ::T # "active" value of λ; on first call = 0 ζ::T # Inexact Newton order parameter: stop when ||∇q|| < ξ * ||g||^(1+ζ) @@ -144,7 +144,8 @@ end function PDataLSKARC( ::Type{S}, ::Type{T}, - n; + n, + m; ζ = T(0.5), ξ = T(0.01), maxtol = T(0.01), @@ -167,9 +168,9 @@ function PDataLSKARC( norm_dirs = S(undef, nshifts) OK = true solver = if solver_method == :cgls - CglsLanczosShiftSolver(n, n, nshifts, S) + CglsLanczosShiftSolver(m, n, nshifts, S) else - LsqrShiftSolver(n, n, nshifts, S) + LsqrShiftSolver(m, n, nshifts, S) end return PDataLSKARC( d, From 446dadbbd514097c890f98acbd1b585ba45c36d1 Mon Sep 17 00:00:00 2001 From: tmigot Date: Wed, 26 Oct 2022 17:13:36 +0200 Subject: [PATCH 4/5] specialize `solve!` --- src/PreProcess/PreProcessLSKARC.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/PreProcess/PreProcessLSKARC.jl b/src/PreProcess/PreProcessLSKARC.jl index c9f337ff..90d880ae 100644 --- a/src/PreProcess/PreProcessLSKARC.jl +++ b/src/PreProcess/PreProcessLSKARC.jl @@ -21,7 +21,7 @@ function preprocess(PData::PDataLSKARC, Jx, Fx, gNorm2, calls, max_calls, α) return false end solver = PData.solver - solve!( + Krylov.solve!( solver, Jx, Fx, From 3926b8fcdbc7a6073a5ad7e9046debad5640c331 Mon Sep 17 00:00:00 2001 From: tmigot Date: Wed, 9 Nov 2022 12:48:00 -0500 Subject: [PATCH 5/5] remove unnecessary import --- src/main.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main.jl b/src/main.jl index a68b0d70..26a2a770 100644 --- a/src/main.jl +++ b/src/main.jl @@ -94,7 +94,7 @@ function preprocess( max_hprod = stp.meta.max_cntrs[:neval_jprod_residual] Fx = workspace.Fx Jx = jac_op_residual!(stp.pb, workspace.xt, workspace.Hstruct.Jv, workspace.Hstruct.Jtv) - PData = ARCTR.preprocess(PData, Jx, Fx, norm_∇f, neval_jprod_residual(stp.pb), max_hprod, α) + PData = preprocess(PData, Jx, Fx, norm_∇f, neval_jprod_residual(stp.pb), max_hprod, α) return PData end