From 9d81e461a77650122fe3ebaf0e52e7a5875a4ac8 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Mon, 23 Oct 2023 15:42:16 -0400 Subject: [PATCH] preserve LM linearsolve structure --- docs/src/solvers/NonlinearSystemSolvers.md | 2 ++ src/jacobian.jl | 35 ++++++++++++---------- src/klement.jl | 16 +++++----- src/levenberg.jl | 6 ++-- 4 files changed, 33 insertions(+), 26 deletions(-) diff --git a/docs/src/solvers/NonlinearSystemSolvers.md b/docs/src/solvers/NonlinearSystemSolvers.md index 89cde3513..5bcfc8c99 100644 --- a/docs/src/solvers/NonlinearSystemSolvers.md +++ b/docs/src/solvers/NonlinearSystemSolvers.md @@ -69,6 +69,8 @@ features, but have a bit of overhead on very small problems. Automatic Jacobian Resetting. This is a fast method but unstable for most problems! - `GeneralKlement()`: Generalization of Klement's Quasi-Newton Method with Line Search and Automatic Jacobian Resetting. This is a fast method but unstable for most problems! + - `LimitedMemoryBroyden()`: An advanced version of `LBroyden` which uses a limited memory + Broyden method. This is a fast method but unstable for most problems! ### SimpleNonlinearSolve.jl diff --git a/src/jacobian.jl b/src/jacobian.jl index e4c2ce011..676678bb2 100644 --- a/src/jacobian.jl +++ b/src/jacobian.jl @@ -50,8 +50,8 @@ jacobian!!(::Number, cache) = last(value_derivative(cache.uf, cache.u)) # Build Jacobian Caches function jacobian_caches(alg::AbstractNonlinearSolveAlgorithm, f, u, p, ::Val{iip}; - linsolve_kwargs = (;), - linsolve_with_JᵀJ::Val{needsJᵀJ} = Val(false)) where {iip, needsJᵀJ} + linsolve_kwargs = (;), lininit::Val{linsolve_init} = Val(true), + linsolve_with_JᵀJ::Val{needsJᵀJ} = Val(false)) where {iip, needsJᵀJ, linsolve_init} uf = JacobianWrapper{iip}(f, p) haslinsolve = hasfield(typeof(alg), :linsolve) @@ -95,25 +95,28 @@ function jacobian_caches(alg::AbstractNonlinearSolveAlgorithm, f, u, p, ::Val{ii Jᵀfu = J' * _vec(fu) end - linprob = LinearProblem(needsJᵀJ ? __maybe_symmetric(JᵀJ) : J, - needsJᵀJ ? _vec(Jᵀfu) : _vec(fu); u0 = _vec(du)) - - if alg isa PseudoTransient - alpha = convert(eltype(u), alg.alpha_initial) - J_new = J - (1 / alpha) * I - linprob = LinearProblem(J_new, _vec(fu); u0 = _vec(du)) + if linsolve_init + linprob_A = alg isa PseudoTransient ? + (J - (1 / (convert(eltype(u), alg.alpha_initial))) * I) : + (needsJᵀJ ? __maybe_symmetric(JᵀJ) : J) + linsolve = __setup_linsolve(linprob_A, needsJᵀJ ? Jᵀfu : fu, du, p, alg) + else + linsolve = nothing end + needsJᵀJ && return uf, linsolve, J, fu, jac_cache, du, JᵀJ, Jᵀfu + return uf, linsolve, J, fu, jac_cache, du +end + +function __setup_linsolve(A, b, u, p, alg) + linprob = LinearProblem(A, _vec(b); u0 = _vec(u)) + weight = similar(u) recursivefill!(weight, true) - Pl, Pr = wrapprecs(alg.precs(needsJᵀJ ? __maybe_symmetric(JᵀJ) : J, nothing, u, p, - nothing, nothing, nothing, nothing, nothing)..., weight) - linsolve = init(linprob, alg.linsolve; alias_A = true, alias_b = true, Pl, Pr, - linsolve_kwargs...) - - needsJᵀJ && return uf, linsolve, J, fu, jac_cache, du, JᵀJ, Jᵀfu - return uf, linsolve, J, fu, jac_cache, du + Pl, Pr = wrapprecs(alg.precs(A, nothing, u, p, nothing, nothing, nothing, nothing, + nothing)..., weight) + return init(linprob, alg.linsolve; alias_A = true, alias_b = true, Pl, Pr) end __get_nonsparse_ad(::AutoSparseForwardDiff) = AutoForwardDiff() diff --git a/src/klement.jl b/src/klement.jl index e4d398445..848370c5a 100644 --- a/src/klement.jl +++ b/src/klement.jl @@ -27,6 +27,10 @@ solves. linesearch end +function set_linsolve(alg::GeneralKlement, linsolve) + return GeneralKlement(alg.max_resets, linsolve, alg.precs, alg.linesearch) +end + function GeneralKlement(; max_resets::Int = 5, linsolve = nothing, linesearch = LineSearch(), precs = DEFAULT_PRECS) linesearch = linesearch isa LineSearch ? linesearch : LineSearch(; method = linesearch) @@ -60,7 +64,7 @@ end get_fu(cache::GeneralKlementCache) = cache.fu -function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::GeneralKlement, args...; +function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg_::GeneralKlement, args...; alias_u0 = false, maxiters = 1000, abstol = 1e-6, internalnorm = DEFAULT_NORM, linsolve_kwargs = (;), kwargs...) where {uType, iip} @unpack f, u0, p = prob @@ -70,17 +74,13 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::GeneralKlemen if u isa Number linsolve = nothing + alg = alg_ else # For General Julia Arrays default to LU Factorization linsolve_alg = alg.linsolve === nothing && u isa Array ? LUFactorization() : nothing - weight = similar(u) - recursivefill!(weight, true) - Pl, Pr = wrapprecs(alg.precs(J, nothing, u, p, nothing, nothing, nothing, nothing, - nothing)..., weight) - linprob = LinearProblem(J, _vec(fu); u0 = _vec(fu)) - linsolve = init(linprob, linsolve_alg; alias_A = true, alias_b = true, Pl, Pr, - linsolve_kwargs...) + alg = set_linsolve(alg_, linsolve_alg) + linsolve = __setup_linsolve(J, _vec(fu), _vec(u), p, alg) end return GeneralKlementCache{iip}(f, alg, u, fu, zero(fu), _mutable_zero(u), p, linsolve, diff --git a/src/levenberg.jl b/src/levenberg.jl index 311f4acdf..3480e6c63 100644 --- a/src/levenberg.jl +++ b/src/levenberg.jl @@ -213,10 +213,12 @@ function SciMLBase.__init(prob::Union{NonlinearProblem{uType, iip}, mat_tmp = zero(JᵀJ) rhs_tmp = nothing else - mat_tmp = similar(JᵀJ, length(fu1) + length(u), length(u)) + # Preserve Types + mat_tmp = vcat(J, DᵀD) fill!(mat_tmp, zero(eltype(u))) - rhs_tmp = similar(mat_tmp, length(fu1) + length(u)) + rhs_tmp = vcat(fu1, u) fill!(rhs_tmp, zero(eltype(u))) + linsolve = __setup_linsolve(mat_tmp, rhs_tmp, u, p, alg) end return LevenbergMarquardtCache{iip, !_unwrap_val(linsolve_with_JᵀJ)}(f, alg, u, fu1,