From 51206a55a8374ab54eababd56be73e30d77f3315 Mon Sep 17 00:00:00 2001 From: Simon Carlson Date: Tue, 17 Jan 2023 11:00:03 +0100 Subject: [PATCH 01/24] Implemented a TrustRegion-solver --- README.md | 2 +- src/NonlinearSolve.jl | 11 +- src/trustRegion.jl | 322 ++++++++++++++++++++++++++++++++++++++++++ src/utils.jl | 28 ++++ test/basictests.jl | 167 ++++++++++++++++++++++ 5 files changed, 524 insertions(+), 6 deletions(-) create mode 100644 src/trustRegion.jl diff --git a/README.md b/README.md index 207bad382..68673a077 100644 --- a/README.md +++ b/README.md @@ -28,7 +28,7 @@ using NonlinearSolve, StaticArrays f(u,p) = u .* u .- 2 u0 = @SVector[1.0, 1.0] probN = NonlinearProblem(f, u0) -solver = solve(probN, NewtonRaphson(), tol = 1e-9) +solver = solve(probN, NewtonRaphson(), abstol = 1e-9) ## Bracketing Methods diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl index 99fc4d61c..7f16bdb4a 100644 --- a/src/NonlinearSolve.jl +++ b/src/NonlinearSolve.jl @@ -31,22 +31,23 @@ end include("utils.jl") include("jacobian.jl") include("raphson.jl") +include("trustRegion.jl") include("ad.jl") import SnoopPrecompile SnoopPrecompile.@precompile_all_calls begin for T in (Float32, Float64) prob = NonlinearProblem{false}((u, p) -> u .* u .- p, T(0.1), T(2)) - for alg in (NewtonRaphson,) - solve(prob, alg(), abstol = T(1e-2)) + for alg in (NewtonRaphson(), TrustRegion(T(10.0))) + solve(prob, alg, abstol = T(1e-2)) end prob = NonlinearProblem{true}((du, u, p) -> du[1] = u[1] * u[1] - p[1], T[0.1], T[2]) - for alg in (NewtonRaphson,) - solve(prob, alg(), abstol = T(1e-2)) + for alg in (NewtonRaphson(), TrustRegion(T(10.0))) + solve(prob, alg, abstol = T(1e-2)) end end end -export NewtonRaphson +export NewtonRaphson, TrustRegion end # module diff --git a/src/trustRegion.jl b/src/trustRegion.jl new file mode 100644 index 000000000..9c75815b6 --- /dev/null +++ b/src/trustRegion.jl @@ -0,0 +1,322 @@ +""" +```julia +TrustRegion(max_trust_radius::Number; chunk_size = Val{0}(), autodiff = Val{true}(), + standardtag = Val{true}(), concrete_jac = nothing, + diff_type = Val{:forward}, linsolve = nothing, precs = DEFAULT_PRECS, + initial_trust_radius::Number = max_trust_radius / 11, + step_threshold::Number = 0.1, + shrink_threshold::Number = 0.25, + expand_threshold::Number = 0.75, + shrink_factor::Number = 0.25, + expand_factor::Number = 2.0, + max_shrink_times::Int = 32) +``` + +An advanced TrustRegion implementation with support for efficient handling of sparse +matrices via colored automatic differentiation and preconditioned linear solvers. Designed +for large-scale and numerically-difficult nonlinear systems. + +### Keyword Arguments + +- `chunk_size`: the chunk size used by the internal ForwardDiff.jl automatic differentiation + system. This allows for multiple derivative columns to be computed simultaneously, + improving performance. Defaults to `0`, which is equivalent to using ForwardDiff.jl's + default chunk size mechanism. For more details, see the documentation for + [ForwardDiff.jl](https://juliadiff.org/ForwardDiff.jl/stable/). +- `autodiff`: whether to use forward-mode automatic differentiation for the Jacobian. + Note that this argument is ignored if an analytical Jacobian is passed, as that will be + used instead. Defaults to `Val{true}`, which means ForwardDiff.jl via + SparseDiffTools.jl is used by default. If `Val{false}`, then FiniteDiff.jl is used for + finite differencing. +- `standardtag`: whether to use a standardized tag definition for the purposes of automatic + differentiation. Defaults to true, which thus uses the `NonlinearSolveTag`. If `Val{false}`, + then ForwardDiff's default function naming tag is used, which results in larger stack + traces. +- `concrete_jac`: whether to build a concrete Jacobian. If a Krylov-subspace method is used, + then the Jacobian will not be constructed and instead direct Jacobian-vector products + `J*v` are computed using forward-mode automatic differentiation or finite differencing + tricks (without ever constructing the Jacobian). However, if the Jacobian is still needed, + for example for a preconditioner, `concrete_jac = true` can be passed in order to force + the construction of the Jacobian. +- `diff_type`: the type of finite differencing used if `autodiff = false`. Defaults to + `Val{:forward}` for forward finite differences. For more details on the choices, see the + [FiniteDiff.jl](https://github.com/JuliaDiff/FiniteDiff.jl) documentation. +- `linsolve`: the [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl) used for the + linear solves within the Newton method. Defaults to `nothing`, which means it uses the + LinearSolve.jl default algorithm choice. For more information on available algorithm + choices, see the [LinearSolve.jl documentation](https://docs.sciml.ai/LinearSolve/stable/). +- `precs`: the choice of preconditioners for the linear solver. Defaults to using no + preconditioners. For more information on specifying preconditioners for LinearSolve + algorithms, consult the + [LinearSolve.jl documentation](https://docs.sciml.ai/LinearSolve/stable/). +- `initial_trust_radius`: the initial trust region radius. Defaults to + `max_trust_radius / 11`. +- `step_threshold`: the threshold for taking a step. In every iteration, the threshold is + compared with a value `r`, which is the actual reduction in the objective function divided + by the predicted reduction. If `step_threshold > r` the model is not a good approximation, + and the step is rejected. Defaults to `0.1`. For more details, see + [Trust-region methods](https://optimization.mccormick.northwestern.edu/index.php/Trust-region_methods) +- `shrink_threshold`: the threshold for shrinking the trust region radius. In every + iteration, the threshold is compared with a value `r` which is the actual reduction in the + objective function divided by the predicted reduction. If `shrink_threshold > r` the trust + region radius is shrunk by `shrink_factor`. Defaults to `0.25`. For more details, see + [Trust-region methods](https://optimization.mccormick.northwestern.edu/index.php/Trust-region_methods) +- `expand_threshold`: the threshold for expanding the trust region radius. If a step is + taken, i.e `step_threshold < r` (with `r` defined in `shrink_threshold`), a check is also + made to see if `expand_threshold < r`. If that is true, the trust region radius is + expanded by `expand_factor`. Defaults to `0.75`. +- `shrink_factor`: the factor to shrink the trust region radius with if + `shrink_threshold > r` (with `r` defined in `shrink_threshold`). Defaults to `0.25`. +- `expand_factor`: the factor to expand the trust region radius with if + `expand_threshold < r` (with `r` defined in `shrink_threshold`). Defaults to `2.0`. +- `max_shrink_times`: the maximum number of times to shrink the trust region radius in a + row, `max_shrink_times` is exceeded, the algorithm returns. Defaults to `32`. + +!!! note + + Currently, the linear solver and chunk size choice only applies to in-place defined + `NonlinearProblem`s. That is expected to change in the future. +""" +struct TrustRegion{CS, AD, FDT, L, P, ST, CJ} <: + AbstractNewtonAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::L + precs::P + max_trust_radius::Number + initial_trust_radius::Number + step_threshold::Number + shrink_threshold::Number + expand_threshold::Number + shrink_factor::Number + expand_factor::Number + max_shrink_times::Int +end + +function TrustRegion(max_trust_radius::Number; chunk_size = Val{0}(), + autodiff = Val{true}(), + standardtag = Val{true}(), concrete_jac = nothing, + diff_type = Val{:forward}, linsolve = nothing, precs = DEFAULT_PRECS, + initial_trust_radius::Number = max_trust_radius / 11, + step_threshold::Number = 0.1, + shrink_threshold::Number = 0.25, + expand_threshold::Number = 0.75, + shrink_factor::Number = 0.25, + expand_factor::Number = 2.0, + max_shrink_times::Int = 32) + TrustRegion{_unwrap_val(chunk_size), _unwrap_val(autodiff), diff_type, + typeof(linsolve), typeof(precs), _unwrap_val(standardtag), + _unwrap_val(concrete_jac)}(linsolve, precs, max_trust_radius::Number, + initial_trust_radius::Number, + step_threshold::Number, + shrink_threshold::Number, + expand_threshold::Number, + shrink_factor::Number, + expand_factor::Number, + max_shrink_times::Int) +end + +mutable struct TrustRegionCache{iip, fType, algType, uType, duType, resType, pType, + INType, tolType, probType, ufType, L, jType, JC + } + f::fType + alg::algType + u::uType + fu::resType + p::pType + uf::ufType + linsolve::L + J::jType + du1::duType + jac_config::JC + iter::Int + force_stop::Bool + maxiters::Int + internalnorm::INType + retcode::SciMLBase.ReturnCode.T + abstol::tolType + prob::probType + trust_r::Number + loss::Number + loss_new::Number + H::jType + g::resType + shrink_counter::Int + step_size::uType + u_new::uType + fu_new::resType + make_new_J::Bool + + function TrustRegionCache{iip}(f::fType, alg::algType, u::uType, fu::resType, p::pType, + uf::ufType, linsolve::L, J::jType, du1::duType, + jac_config::JC, iter::Int, + force_stop::Bool, maxiters::Int, internalnorm::INType, + retcode::SciMLBase.ReturnCode.T, abstol::tolType, + prob::probType, trust_r::Number, loss::Number, + loss_new::Number, H::jType, g::resType, + shrink_counter::Int, step_size::uType, u_new::uType, + fu_new::resType, + make_new_J::Bool) where {iip, fType, algType, uType, + duType, resType, pType, INType, + tolType, probType, ufType, L, + jType, JC} + new{iip, fType, algType, uType, duType, resType, pType, + INType, tolType, probType, ufType, L, jType, JC + }(f, alg, u, fu, p, uf, linsolve, J, + du1, jac_config, iter, force_stop, + maxiters, internalnorm, retcode, + abstol, prob, trust_r, loss, + loss_new, H, g, shrink_counter, + step_size, u_new, fu_new, + make_new_J) + end +end + +function jacobian_caches(alg::TrustRegion, f, u, p, ::Val{true}) + uf = JacobianWrapper(f, p) + J = ArrayInterfaceCore.undefmatrix(u) + + linprob = LinearProblem(J, _vec(zero(u)); u0 = _vec(zero(u))) + weight = similar(u) + recursivefill!(weight, false) + + Pl, Pr = wrapprecs(alg.precs(J, nothing, u, p, nothing, nothing, nothing, nothing, + nothing)..., weight) + linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, + Pl = Pl, Pr = Pr) + + du1 = zero(u) + du2 = zero(u) + tmp = zero(u) + jac_config = build_jac_config(alg, f, uf, du1, u, tmp, du2) + + uf, linsolve, J, du1, jac_config +end + +function jacobian_caches(alg::TrustRegion, f, u, p, ::Val{false}) + J = ArrayInterfaceCore.undefmatrix(u) + JacobianWrapper(f, p), nothing, J, nothing, nothing +end + +function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::TrustRegion, + args...; + alias_u0 = false, + maxiters = 1000, + abstol = 1e-6, + internalnorm = DEFAULT_NORM, + kwargs...) where {uType, iip} + if alias_u0 + u = prob.u0 + else + u = deepcopy(prob.u0) + end + f = prob.f + p = prob.p + if iip + fu = zero(u) + f(fu, u, p) + else + fu = f(u, p) + end + + loss = get_loss(fu) + uf, linsolve, J, du1, jac_config = jacobian_caches(alg, f, u, p, Val(iip)) + + return TrustRegionCache{iip}(f, alg, u, fu, p, uf, linsolve, J, du1, jac_config, + 1, false, maxiters, internalnorm, + ReturnCode.Default, abstol, prob, alg.initial_trust_radius, + loss, loss, J, fu, 0, u, u, fu, true) +end + +function perform_step!(cache::TrustRegionCache{true}) + @unpack make_new_J, J, fu, f, u, p = cache + if cache.make_new_J + jacobian!(J, cache) + cache.H = J * J + cache.g = J * fu + end + + dogleg!(cache) + + # Compute the potentially new u + cache.u_new = u .+ cache.step_size + f(cache.fu_new, cache.u_new, p) + + trust_region_step!(cache) + return nothing +end + +function perform_step!(cache::TrustRegionCache{false}) + @unpack make_new_J, fu, f, u, p = cache + + if make_new_J + J = jacobian(cache, f) + cache.H = J * J + cache.g = J * fu + end + dogleg!(cache) + + # Compute the potentially new u + cache.u_new = u .+ cache.step_size + cache.fu_new = f(cache.u_new, p) + + trust_region_step!(cache) + return nothing +end + +function trust_region_step!(cache::TrustRegionCache) + @unpack fu_new, u_new, step_size, g, H, loss, alg = cache + cache.loss_new = get_loss(fu_new) + + # Compute the ratio of the actual reduction to the predicted reduction. + model = -(step_size' * g + 0.5 * step_size' * H * step_size) + r = (loss - cache.loss_new) / model + + # Update the trust region radius. + if r < alg.shrink_threshold + cache.trust_r *= alg.shrink_factor + cache.shrink_counter += 1 + else + cache.shrink_counter = 0 + end + if r > alg.step_threshold + + # Take the step. + cache.u = u_new + cache.fu = fu_new + cache.loss = cache.loss_new + + # Update the trust region radius. + if r > alg.expand_threshold + cache.trust_r = min(alg.expand_factor * cache.trust_r, alg.max_trust_radius) + end + + cache.make_new_J = true + else + # No need to make a new J, no step was taken, so we try again with a smaller trust_r + cache.make_new_J = false + end + + if iszero(cache.fu) || cache.internalnorm(cache.step_size) < cache.abstol + cache.force_stop = true + end +end + +function SciMLBase.solve!(cache::TrustRegionCache) + while !cache.force_stop && cache.iter < cache.maxiters && + cache.shrink_counter < cache.alg.max_shrink_times + perform_step!(cache) + cache.iter += 1 + end + + if cache.iter == cache.maxiters + cache.retcode = ReturnCode.MaxIters + else + cache.retcode = ReturnCode.Success + end + + SciMLBase.build_solution(cache.prob, cache.alg, cache.u, cache.fu; + retcode = cache.retcode) +end + +function get_loss(fu) + return 0.5 * norm(fu)^2 +end diff --git a/src/utils.jl b/src/utils.jl index f4946091b..4d0429c00 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -124,3 +124,31 @@ function _nfcount(N, ::Type{diff_type}) where {diff_type} end tmp end + +function dogleg!(cache) + @unpack g, H, trust_r = cache + # Compute the Newton step. + δN = -H \ g + # Test if the full step is within the trust region. + if norm(δN) ≤ trust_r + cache.step_size = δN + return + end + + # Calcualte Cauchy point, optimum along the steepest descent direction. + δsd = -g + norm_δsd = norm(δsd) + if norm_δsd ≥ trust_r + cache.step_size = δsd .* trust_r / norm_δsd + return + end + + # Find the intersection point on the boundary. + N_sd = δN - δsd + dot_N_sd = dot(N_sd, N_sd) + dot_sd_N_sd = dot(δsd, N_sd) + dot_sd = dot(δsd, δsd) + fact = dot_sd_N_sd^2 - dot_N_sd * (dot_sd - trust_r^2) + τ = (-dot_sd_N_sd + sqrt(fact)) / dot_N_sd + cache.step_size = δsd + τ * N_sd +end diff --git a/test/basictests.jl b/test/basictests.jl index 85b093c1c..17876d722 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -3,6 +3,8 @@ using StaticArrays using BenchmarkTools using Test +# --- NewtonRaphson tests --- + function benchmark_immutable(f, u0) probN = NonlinearProblem{false}(f, u0) solver = init(probN, NewtonRaphson(), abstol = 1e-9) @@ -121,3 +123,168 @@ for u0 in [1.0, [1, 1.0]] @test solve(probN, NewtonRaphson()).u ≈ sol @test solve(probN, NewtonRaphson(; autodiff = false)).u ≈ sol end + +# --- TrustRegion tests --- + +function benchmark_immutable(f, u0) + probN = NonlinearProblem(f, u0) + solver = init(probN, TrustRegion(10.0), abstol = 1e-9) + sol = solve!(solver) +end + +function benchmark_mutable(f, u0) + probN = NonlinearProblem(f, u0) + solver = init(probN, TrustRegion(10.0), abstol = 1e-9) + sol = solve!(solver) +end + +function benchmark_scalar(f, u0) + probN = NonlinearProblem(f, u0) + sol = (solve(probN, TrustRegion(10.0))) +end + +function ff(u, p) + u .* u .- 2 +end + +function sf(u, p) + u * u - 2 +end + +sol = benchmark_immutable(ff, cu0) +@test sol.retcode === ReturnCode.Success +@test all(sol.u .* sol.u .- 2 .< 1e-9) +sol = benchmark_mutable(ff, cu0) +@test sol.retcode === ReturnCode.Success +@test all(sol.u .* sol.u .- 2 .< 1e-9) +sol = benchmark_scalar(sf, csu0) +@test sol.retcode === ReturnCode.Success +@test sol.u * sol.u - 2 < 1e-9 + +function benchmark_inplace(f, u0) + probN = NonlinearProblem(f, u0) + solver = init(probN, TrustRegion(10.0), abstol = 1e-9) + sol = solve!(solver) +end + +function ffiip(du, u, p) + du .= u .* u .- 2 +end +u0 = [1.0, 1.0] + +sol = benchmark_inplace(ffiip, u0) +@test sol.retcode === ReturnCode.Success +@test all(sol.u .* sol.u .- 2 .< 1e-9) + +# AD Tests +using ForwardDiff + +# Immutable +f, u0 = (u, p) -> u .* u .- p, @SVector[1.0, 1.0] + +g = function (p) + probN = NonlinearProblem(f, csu0, p) + sol = solve(probN, TrustRegion(10.0), abstol = 1e-9) + return sol.u[end] +end + +for p in 1.1:0.1:100.0 + @test g(p) ≈ sqrt(p) + @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) +end + +# Scalar +f, u0 = (u, p) -> u * u - p, 1.0 + +g = function (p) + probN = NonlinearProblem(f, oftype(p, u0), p) + sol = solve(probN, TrustRegion(10.0), abstol = 1e-10) + return sol.u +end + +@test ForwardDiff.derivative(g, 3.0) ≈ 1 / (2 * sqrt(3.0)) + +for p in 1.1:0.1:100.0 + @test g(p) ≈ sqrt(p) + @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) +end + +f = (u, p) -> p[1] * u * u - p[2] +t = (p) -> [sqrt(p[2] / p[1])] +p = [0.9, 50.0] +gnewton = function (p) + probN = NonlinearProblem(f, 0.5, p) + sol = solve(probN, TrustRegion(10.0)) + return [sol.u] +end +@test gnewton(p) ≈ [sqrt(p[2] / p[1])] +@test ForwardDiff.jacobian(gnewton, p) ≈ ForwardDiff.jacobian(t, p) + +# Error Checks + +f, u0 = (u, p) -> u .* u .- 2.0, @SVector[1.0, 1.0] +probN = NonlinearProblem(f, u0) + +@test solve(probN, TrustRegion(10.0)).u[end] ≈ sqrt(2.0) +@test solve(probN, TrustRegion(10.0; autodiff = false)).u[end] ≈ sqrt(2.0) + +for u0 in [1.0, [1, 1.0]] + local f, probN, sol + f = (u, p) -> u .* u .- 2.0 + probN = NonlinearProblem(f, u0) + sol = sqrt(2) * u0 + + @test solve(probN, TrustRegion(10.0)).u ≈ sol + @test solve(probN, TrustRegion(10.0)).u ≈ sol + @test solve(probN, TrustRegion(10.0; autodiff = false)).u ≈ sol +end + +# Test that `TrustRegion` passes a test that `NewtonRaphson` fails on. +u0 = [-10.0, -1.0, 1.0, 2.0, 3.0, 4.0, 10.0] +global g, f +f = (u, p) -> 0.010000000000000002 .+ + 10.000000000000002 ./ (1 .+ + (0.21640425613334457 .+ + 216.40425613334457 ./ (1 .+ + (0.21640425613334457 .+ + 216.40425613334457 ./ + (1 .+ 0.0006250000000000001(u .^ 2.0))) .^ 2.0)) .^ 2.0) .- + 0.0011552453009332421u .- p +g = function (p) + probN = NonlinearProblem{false}(f, u0, p) + sol = solve(probN, TrustRegion(100.0)) + return sol.u +end +p = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] +u = g(p) +f(u, p) +@test all(abs.(f(u, p)) .< 1e-10) + +# Test kwars in `TrustRegion` +max_trust_radius = [10.0, 100.0, 1000.0] +initial_trust_radius = [10.0, 1.0, 0.1] +step_threshold = [0.0, 0.01, 0.25] +shrink_threshold = [0.25, 0.3, 0.5] +expand_threshold = [0.5, 0.8, 0.9] +shrink_factor = [0.1, 0.3, 0.5] +expand_factor = [1.5, 2.0, 3.0] +max_shrink_times = [10, 20, 30] + +list_of_options = zip(max_trust_radius, initial_trust_radius, step_threshold, + shrink_threshold, expand_threshold, shrink_factor, + expand_factor, max_shrink_times) +for options in list_of_options + local probN, sol, alg + alg = TrustRegion(options[1]; + initial_trust_radius = options[2], + step_threshold = options[3], + shrink_threshold = options[4], + expand_threshold = options[5], + shrink_factor = options[6], + expand_factor = options[7], + max_shrink_times = options[8]) + + probN = NonlinearProblem(f, u0, p) + sol = solve(probN, alg) + @test all(f(u, p) .< 1e-10) +end From 3f434c59368663721db765fb469ddc1ae2e59a5a Mon Sep 17 00:00:00 2001 From: Simon Carlson Date: Tue, 17 Jan 2023 11:14:45 +0100 Subject: [PATCH 02/24] Moving functions. --- src/trustRegion.jl | 32 ++++++++++++++++++++++++++++---- src/utils.jl | 28 ++-------------------------- 2 files changed, 30 insertions(+), 30 deletions(-) diff --git a/src/trustRegion.jl b/src/trustRegion.jl index 9c75815b6..b12ef7c92 100644 --- a/src/trustRegion.jl +++ b/src/trustRegion.jl @@ -300,6 +300,34 @@ function trust_region_step!(cache::TrustRegionCache) end end +function dogleg!(cache::TrustRegionCache) + @unpack g, H, trust_r = cache + # Compute the Newton step. + δN = -H \ g + # Test if the full step is within the trust region. + if norm(δN) ≤ trust_r + cache.step_size = δN + return + end + + # Calcualte Cauchy point, optimum along the steepest descent direction. + δsd = -g + norm_δsd = norm(δsd) + if norm_δsd ≥ trust_r + cache.step_size = δsd .* trust_r / norm_δsd + return + end + + # Find the intersection point on the boundary. + N_sd = δN - δsd + dot_N_sd = dot(N_sd, N_sd) + dot_sd_N_sd = dot(δsd, N_sd) + dot_sd = dot(δsd, δsd) + fact = dot_sd_N_sd^2 - dot_N_sd * (dot_sd - trust_r^2) + τ = (-dot_sd_N_sd + sqrt(fact)) / dot_N_sd + cache.step_size = δsd + τ * N_sd +end + function SciMLBase.solve!(cache::TrustRegionCache) while !cache.force_stop && cache.iter < cache.maxiters && cache.shrink_counter < cache.alg.max_shrink_times @@ -316,7 +344,3 @@ function SciMLBase.solve!(cache::TrustRegionCache) SciMLBase.build_solution(cache.prob, cache.alg, cache.u, cache.fu; retcode = cache.retcode) end - -function get_loss(fu) - return 0.5 * norm(fu)^2 -end diff --git a/src/utils.jl b/src/utils.jl index 4d0429c00..9cb414f34 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -125,30 +125,6 @@ function _nfcount(N, ::Type{diff_type}) where {diff_type} tmp end -function dogleg!(cache) - @unpack g, H, trust_r = cache - # Compute the Newton step. - δN = -H \ g - # Test if the full step is within the trust region. - if norm(δN) ≤ trust_r - cache.step_size = δN - return - end - - # Calcualte Cauchy point, optimum along the steepest descent direction. - δsd = -g - norm_δsd = norm(δsd) - if norm_δsd ≥ trust_r - cache.step_size = δsd .* trust_r / norm_δsd - return - end - - # Find the intersection point on the boundary. - N_sd = δN - δsd - dot_N_sd = dot(N_sd, N_sd) - dot_sd_N_sd = dot(δsd, N_sd) - dot_sd = dot(δsd, δsd) - fact = dot_sd_N_sd^2 - dot_N_sd * (dot_sd - trust_r^2) - τ = (-dot_sd_N_sd + sqrt(fact)) / dot_N_sd - cache.step_size = δsd + τ * N_sd +function get_loss(fu) + return 0.5 * norm(fu)^2 end From 76edc1fdff0372389b140a62863d6886b7237e65 Mon Sep 17 00:00:00 2001 From: Simon Carlson Date: Tue, 17 Jan 2023 11:19:16 +0100 Subject: [PATCH 03/24] Format --- test/basictests.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/basictests.jl b/test/basictests.jl index 17876d722..4eba12f83 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -221,7 +221,6 @@ end @test ForwardDiff.jacobian(gnewton, p) ≈ ForwardDiff.jacobian(t, p) # Error Checks - f, u0 = (u, p) -> u .* u .- 2.0, @SVector[1.0, 1.0] probN = NonlinearProblem(f, u0) From 569599a656a2890d1bb7876c50d612162a0f7770 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 17 Jan 2023 08:20:14 -0500 Subject: [PATCH 04/24] Update src/trustRegion.jl --- src/trustRegion.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/trustRegion.jl b/src/trustRegion.jl index b12ef7c92..70bb5a295 100644 --- a/src/trustRegion.jl +++ b/src/trustRegion.jl @@ -249,8 +249,8 @@ function perform_step!(cache::TrustRegionCache{false}) if make_new_J J = jacobian(cache, f) - cache.H = J * J - cache.g = J * fu + mul!(cache.H,J,J) + mul!(cache.g,J,fu) end dogleg!(cache) From 18abfba80eda8714e4c2d92e678f4cbe88fb896c Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 17 Jan 2023 08:21:11 -0500 Subject: [PATCH 05/24] Update src/trustRegion.jl --- src/trustRegion.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/trustRegion.jl b/src/trustRegion.jl index 70bb5a295..b12ef7c92 100644 --- a/src/trustRegion.jl +++ b/src/trustRegion.jl @@ -249,8 +249,8 @@ function perform_step!(cache::TrustRegionCache{false}) if make_new_J J = jacobian(cache, f) - mul!(cache.H,J,J) - mul!(cache.g,J,fu) + cache.H = J * J + cache.g = J * fu end dogleg!(cache) From 7caf196fa19699432dd0c2f9f86d7c2860518fd1 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 17 Jan 2023 08:21:18 -0500 Subject: [PATCH 06/24] Update src/trustRegion.jl --- src/trustRegion.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/trustRegion.jl b/src/trustRegion.jl index b12ef7c92..0fc84d25d 100644 --- a/src/trustRegion.jl +++ b/src/trustRegion.jl @@ -230,8 +230,8 @@ function perform_step!(cache::TrustRegionCache{true}) @unpack make_new_J, J, fu, f, u, p = cache if cache.make_new_J jacobian!(J, cache) - cache.H = J * J - cache.g = J * fu + mul!(cache.H,J,J) + mul!(cache.g,J,fu) end dogleg!(cache) From c532be2c39a8091c421692212a9ab8d21d1c0472 Mon Sep 17 00:00:00 2001 From: Simon Carlson Date: Tue, 17 Jan 2023 19:11:23 +0100 Subject: [PATCH 07/24] fixing comments in #116 PR --- src/NonlinearSolve.jl | 8 +-- src/trustRegion.jl | 120 +++++++++++++++++++++++++----------------- src/utils.jl | 2 +- test/basictests.jl | 34 ++++++------ 4 files changed, 96 insertions(+), 68 deletions(-) diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl index 7f16bdb4a..ceec75389 100644 --- a/src/NonlinearSolve.jl +++ b/src/NonlinearSolve.jl @@ -38,13 +38,13 @@ import SnoopPrecompile SnoopPrecompile.@precompile_all_calls begin for T in (Float32, Float64) prob = NonlinearProblem{false}((u, p) -> u .* u .- p, T(0.1), T(2)) - for alg in (NewtonRaphson(), TrustRegion(T(10.0))) - solve(prob, alg, abstol = T(1e-2)) + for alg in (NewtonRaphson, TrustRegion) + solve(prob, alg(), abstol = T(1e-2)) end prob = NonlinearProblem{true}((du, u, p) -> du[1] = u[1] * u[1] - p[1], T[0.1], T[2]) - for alg in (NewtonRaphson(), TrustRegion(T(10.0))) - solve(prob, alg, abstol = T(1e-2)) + for alg in (NewtonRaphson, TrustRegion) + solve(prob, alg(), abstol = T(1e-2)) end end end diff --git a/src/trustRegion.jl b/src/trustRegion.jl index b12ef7c92..044940d9b 100644 --- a/src/trustRegion.jl +++ b/src/trustRegion.jl @@ -1,8 +1,9 @@ """ ```julia -TrustRegion(max_trust_radius::Number; chunk_size = Val{0}(), autodiff = Val{true}(), +TrustRegion(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(), concrete_jac = nothing, diff_type = Val{:forward}, linsolve = nothing, precs = DEFAULT_PRECS, + max_trust_radius::Number = nothing, initial_trust_radius::Number = max_trust_radius / 11, step_threshold::Number = 0.1, shrink_threshold::Number = 0.25, @@ -49,6 +50,7 @@ for large-scale and numerically-difficult nonlinear systems. preconditioners. For more information on specifying preconditioners for LinearSolve algorithms, consult the [LinearSolve.jl documentation](https://docs.sciml.ai/LinearSolve/stable/). +- `max_trust_radius`: the maximal trust region radius. Defaults to nothing. - `initial_trust_radius`: the initial trust region radius. Defaults to `max_trust_radius / 11`. - `step_threshold`: the threshold for taking a step. In every iteration, the threshold is @@ -77,41 +79,43 @@ for large-scale and numerically-difficult nonlinear systems. Currently, the linear solver and chunk size choice only applies to in-place defined `NonlinearProblem`s. That is expected to change in the future. """ -struct TrustRegion{CS, AD, FDT, L, P, ST, CJ} <: +struct TrustRegion{CS, AD, FDT, L, P, ST, CJ, MTR} <: AbstractNewtonAlgorithm{CS, AD, FDT, ST, CJ} linsolve::L precs::P - max_trust_radius::Number - initial_trust_radius::Number - step_threshold::Number - shrink_threshold::Number - expand_threshold::Number - shrink_factor::Number - expand_factor::Number + max_trust_radius::MTR + initial_trust_radius::MTR + step_threshold::MTR + shrink_threshold::MTR + expand_threshold::MTR + shrink_factor::MTR + expand_factor::MTR max_shrink_times::Int end -function TrustRegion(max_trust_radius::Number; chunk_size = Val{0}(), +function TrustRegion(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(), concrete_jac = nothing, diff_type = Val{:forward}, linsolve = nothing, precs = DEFAULT_PRECS, - initial_trust_radius::Number = max_trust_radius / 11, - step_threshold::Number = 0.1, - shrink_threshold::Number = 0.25, - expand_threshold::Number = 0.75, - shrink_factor::Number = 0.25, - expand_factor::Number = 2.0, + max_trust_radius::Real = 0.0, + initial_trust_radius::Real = 0.0, + step_threshold::Real = 0.1, + shrink_threshold::Real = 0.25, + expand_threshold::Real = 0.75, + shrink_factor::Real = 0.25, + expand_factor::Real = 2.0, max_shrink_times::Int = 32) TrustRegion{_unwrap_val(chunk_size), _unwrap_val(autodiff), diff_type, typeof(linsolve), typeof(precs), _unwrap_val(standardtag), - _unwrap_val(concrete_jac)}(linsolve, precs, max_trust_radius::Number, - initial_trust_radius::Number, - step_threshold::Number, - shrink_threshold::Number, - expand_threshold::Number, - shrink_factor::Number, - expand_factor::Number, - max_shrink_times::Int) + _unwrap_val(concrete_jac), typeof(max_trust_radius) + }(linsolve, precs, max_trust_radius, + initial_trust_radius, + step_threshold, + shrink_threshold, + expand_threshold, + shrink_factor, + expand_factor, + max_shrink_times) end mutable struct TrustRegionCache{iip, fType, algType, uType, duType, resType, pType, @@ -135,6 +139,7 @@ mutable struct TrustRegionCache{iip, fType, algType, uType, duType, resType, pTy abstol::tolType prob::probType trust_r::Number + max_trust_r::Number loss::Number loss_new::Number H::jType @@ -144,29 +149,32 @@ mutable struct TrustRegionCache{iip, fType, algType, uType, duType, resType, pTy u_new::uType fu_new::resType make_new_J::Bool + r::Real function TrustRegionCache{iip}(f::fType, alg::algType, u::uType, fu::resType, p::pType, uf::ufType, linsolve::L, J::jType, du1::duType, jac_config::JC, iter::Int, force_stop::Bool, maxiters::Int, internalnorm::INType, retcode::SciMLBase.ReturnCode.T, abstol::tolType, - prob::probType, trust_r::Number, loss::Number, + prob::probType, trust_r::Number, max_trust_r::Number, + loss::Number, loss_new::Number, H::jType, g::resType, shrink_counter::Int, step_size::uType, u_new::uType, fu_new::resType, - make_new_J::Bool) where {iip, fType, algType, uType, - duType, resType, pType, INType, - tolType, probType, ufType, L, - jType, JC} + make_new_J::Bool, + r::Real) where {iip, fType, algType, uType, + duType, resType, pType, INType, + tolType, probType, ufType, L, + jType, JC} new{iip, fType, algType, uType, duType, resType, pType, INType, tolType, probType, ufType, L, jType, JC }(f, alg, u, fu, p, uf, linsolve, J, du1, jac_config, iter, force_stop, maxiters, internalnorm, retcode, - abstol, prob, trust_r, loss, + abstol, prob, trust_r, max_trust_r, loss, loss_new, H, g, shrink_counter, step_size, u_new, fu_new, - make_new_J) + make_new_J, r) end end @@ -220,24 +228,40 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::TrustRegion, loss = get_loss(fu) uf, linsolve, J, du1, jac_config = jacobian_caches(alg, f, u, p, Val(iip)) + # Set default trust region radius if not specified + max_trust_radius = alg.max_trust_radius + initial_trust_radius = alg.initial_trust_radius + if max_trust_radius == 0.0 + max_trust_radius = max(norm(fu), maximum(u) - minimum(u)) + end + if initial_trust_radius == 0.0 + initial_trust_radius = max_trust_radius / 11 + end + return TrustRegionCache{iip}(f, alg, u, fu, p, uf, linsolve, J, du1, jac_config, 1, false, maxiters, internalnorm, - ReturnCode.Default, abstol, prob, alg.initial_trust_radius, - loss, loss, J, fu, 0, u, u, fu, true) + ReturnCode.Default, abstol, prob, initial_trust_radius, + max_trust_radius, loss, loss, J, fu, 0, u, u, fu, true, + 0.0) end function perform_step!(cache::TrustRegionCache{true}) - @unpack make_new_J, J, fu, f, u, p = cache + @unpack make_new_J, J, fu, f, u, p, trust_r, du1, alg, linsolve = cache if cache.make_new_J jacobian!(J, cache) cache.H = J * J cache.g = J * fu end - dogleg!(cache) + # u = u - J \ fu + linres = dolinsolve(alg.precs, linsolve, A = cache.H, b = _vec(cache.g), + linu = _vec(du1), + p = p, reltol = cache.abstol) + cache.linsolve = linres.cache + dogleg!(-du1, trust_r, cache) # Compute the potentially new u - cache.u_new = u .+ cache.step_size + cache.u_new .= u .+ cache.step_size f(cache.fu_new, cache.u_new, p) trust_region_step!(cache) @@ -245,14 +269,18 @@ function perform_step!(cache::TrustRegionCache{true}) end function perform_step!(cache::TrustRegionCache{false}) - @unpack make_new_J, fu, f, u, p = cache + @unpack make_new_J, fu, f, u, p, trust_r = cache if make_new_J J = jacobian(cache, f) cache.H = J * J cache.g = J * fu end - dogleg!(cache) + + @unpack g, H = cache + # Compute the Newton step. + δN = -H \ g + dogleg!(δN, trust_r, cache) # Compute the potentially new u cache.u_new = u .+ cache.step_size @@ -263,12 +291,12 @@ function perform_step!(cache::TrustRegionCache{false}) end function trust_region_step!(cache::TrustRegionCache) - @unpack fu_new, u_new, step_size, g, H, loss, alg = cache + @unpack fu_new, u_new, step_size, g, H, loss, alg, max_trust_r = cache cache.loss_new = get_loss(fu_new) # Compute the ratio of the actual reduction to the predicted reduction. - model = -(step_size' * g + 0.5 * step_size' * H * step_size) - r = (loss - cache.loss_new) / model + cache.r = -(loss - cache.loss_new) / (step_size' * g + step_size' * H * step_size / 2) + @unpack r = cache # Update the trust region radius. if r < alg.shrink_threshold @@ -286,7 +314,7 @@ function trust_region_step!(cache::TrustRegionCache) # Update the trust region radius. if r > alg.expand_threshold - cache.trust_r = min(alg.expand_factor * cache.trust_r, alg.max_trust_radius) + cache.trust_r = min(alg.expand_factor * cache.trust_r, max_trust_r) end cache.make_new_J = true @@ -300,10 +328,8 @@ function trust_region_step!(cache::TrustRegionCache) end end -function dogleg!(cache::TrustRegionCache) - @unpack g, H, trust_r = cache - # Compute the Newton step. - δN = -H \ g +function dogleg!(δN, trust_r, cache::TrustRegionCache) + # Test if the full step is within the trust region. if norm(δN) ≤ trust_r cache.step_size = δN @@ -311,7 +337,7 @@ function dogleg!(cache::TrustRegionCache) end # Calcualte Cauchy point, optimum along the steepest descent direction. - δsd = -g + δsd = -cache.g norm_δsd = norm(δsd) if norm_δsd ≥ trust_r cache.step_size = δsd .* trust_r / norm_δsd diff --git a/src/utils.jl b/src/utils.jl index 9cb414f34..4118b525c 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -126,5 +126,5 @@ function _nfcount(N, ::Type{diff_type}) where {diff_type} end function get_loss(fu) - return 0.5 * norm(fu)^2 + return norm(fu)^2 / 2 end diff --git a/test/basictests.jl b/test/basictests.jl index 4eba12f83..878aa48e3 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -128,19 +128,19 @@ end function benchmark_immutable(f, u0) probN = NonlinearProblem(f, u0) - solver = init(probN, TrustRegion(10.0), abstol = 1e-9) + solver = init(probN, TrustRegion(), abstol = 1e-9) sol = solve!(solver) end function benchmark_mutable(f, u0) probN = NonlinearProblem(f, u0) - solver = init(probN, TrustRegion(10.0), abstol = 1e-9) + solver = init(probN, TrustRegion(), abstol = 1e-9) sol = solve!(solver) end function benchmark_scalar(f, u0) probN = NonlinearProblem(f, u0) - sol = (solve(probN, TrustRegion(10.0))) + sol = (solve(probN, TrustRegion())) end function ff(u, p) @@ -163,7 +163,7 @@ sol = benchmark_scalar(sf, csu0) function benchmark_inplace(f, u0) probN = NonlinearProblem(f, u0) - solver = init(probN, TrustRegion(10.0), abstol = 1e-9) + solver = init(probN, TrustRegion(), abstol = 1e-9) sol = solve!(solver) end @@ -184,11 +184,12 @@ f, u0 = (u, p) -> u .* u .- p, @SVector[1.0, 1.0] g = function (p) probN = NonlinearProblem(f, csu0, p) - sol = solve(probN, TrustRegion(10.0), abstol = 1e-9) + sol = solve(probN, TrustRegion(), abstol = 1e-9) return sol.u[end] end -for p in 1.1:0.1:100.0 +# TODO look in the the forward diff failure when p = 100.0 +for p in 1.1:0.1:99.9 @test g(p) ≈ sqrt(p) @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) end @@ -198,13 +199,14 @@ f, u0 = (u, p) -> u * u - p, 1.0 g = function (p) probN = NonlinearProblem(f, oftype(p, u0), p) - sol = solve(probN, TrustRegion(10.0), abstol = 1e-10) + sol = solve(probN, TrustRegion(), abstol = 1e-10) return sol.u end @test ForwardDiff.derivative(g, 3.0) ≈ 1 / (2 * sqrt(3.0)) -for p in 1.1:0.1:100.0 +# TODO look in the the forward diff failure when p = 100.0 +for p in 1.1:0.1:99.9 @test g(p) ≈ sqrt(p) @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) end @@ -214,7 +216,7 @@ t = (p) -> [sqrt(p[2] / p[1])] p = [0.9, 50.0] gnewton = function (p) probN = NonlinearProblem(f, 0.5, p) - sol = solve(probN, TrustRegion(10.0)) + sol = solve(probN, TrustRegion()) return [sol.u] end @test gnewton(p) ≈ [sqrt(p[2] / p[1])] @@ -224,8 +226,8 @@ end f, u0 = (u, p) -> u .* u .- 2.0, @SVector[1.0, 1.0] probN = NonlinearProblem(f, u0) -@test solve(probN, TrustRegion(10.0)).u[end] ≈ sqrt(2.0) -@test solve(probN, TrustRegion(10.0; autodiff = false)).u[end] ≈ sqrt(2.0) +@test solve(probN, TrustRegion()).u[end] ≈ sqrt(2.0) +@test solve(probN, TrustRegion(; autodiff = false)).u[end] ≈ sqrt(2.0) for u0 in [1.0, [1, 1.0]] local f, probN, sol @@ -233,9 +235,9 @@ for u0 in [1.0, [1, 1.0]] probN = NonlinearProblem(f, u0) sol = sqrt(2) * u0 - @test solve(probN, TrustRegion(10.0)).u ≈ sol - @test solve(probN, TrustRegion(10.0)).u ≈ sol - @test solve(probN, TrustRegion(10.0; autodiff = false)).u ≈ sol + @test solve(probN, TrustRegion()).u ≈ sol + @test solve(probN, TrustRegion()).u ≈ sol + @test solve(probN, TrustRegion(; autodiff = false)).u ≈ sol end # Test that `TrustRegion` passes a test that `NewtonRaphson` fails on. @@ -251,7 +253,7 @@ f = (u, p) -> 0.010000000000000002 .+ 0.0011552453009332421u .- p g = function (p) probN = NonlinearProblem{false}(f, u0, p) - sol = solve(probN, TrustRegion(100.0)) + sol = solve(probN, TrustRegion()) return sol.u end p = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] @@ -274,7 +276,7 @@ list_of_options = zip(max_trust_radius, initial_trust_radius, step_threshold, expand_factor, max_shrink_times) for options in list_of_options local probN, sol, alg - alg = TrustRegion(options[1]; + alg = TrustRegion(max_trust_radius = options[1], initial_trust_radius = options[2], step_threshold = options[3], shrink_threshold = options[4], From e4c91f1d884ac4fcedcfb53ee9d5e8bd7a4c0337 Mon Sep 17 00:00:00 2001 From: Simon Carlson Date: Tue, 17 Jan 2023 19:14:32 +0100 Subject: [PATCH 08/24] fix misspell --- test/basictests.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/basictests.jl b/test/basictests.jl index 878aa48e3..157b42a40 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -188,7 +188,7 @@ g = function (p) return sol.u[end] end -# TODO look in the the forward diff failure when p = 100.0 +# TODO look in to the the forward diff failure when p = 100.0 for p in 1.1:0.1:99.9 @test g(p) ≈ sqrt(p) @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) @@ -205,7 +205,7 @@ end @test ForwardDiff.derivative(g, 3.0) ≈ 1 / (2 * sqrt(3.0)) -# TODO look in the the forward diff failure when p = 100.0 +# TODO look in to the the forward diff failure when p = 100.0 for p in 1.1:0.1:99.9 @test g(p) ≈ sqrt(p) @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) From 2043c9520da62bcde09ff2e14f3ba3569533d195 Mon Sep 17 00:00:00 2001 From: Simon Carlson Date: Tue, 17 Jan 2023 19:15:02 +0100 Subject: [PATCH 09/24] fix misspell --- test/basictests.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/basictests.jl b/test/basictests.jl index 157b42a40..bc8a480da 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -188,7 +188,7 @@ g = function (p) return sol.u[end] end -# TODO look in to the the forward diff failure when p = 100.0 +# TODO look into the the forward diff failure when p = 100.0 for p in 1.1:0.1:99.9 @test g(p) ≈ sqrt(p) @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) @@ -205,7 +205,7 @@ end @test ForwardDiff.derivative(g, 3.0) ≈ 1 / (2 * sqrt(3.0)) -# TODO look in to the the forward diff failure when p = 100.0 +# TODO look into the the forward diff failure when p = 100.0 for p in 1.1:0.1:99.9 @test g(p) ≈ sqrt(p) @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) From 8e89003d23df598dcf0730980a6189c3e2f408df Mon Sep 17 00:00:00 2001 From: Simon Carlson Date: Tue, 17 Jan 2023 19:20:54 +0100 Subject: [PATCH 10/24] fix docs --- src/trustRegion.jl | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/src/trustRegion.jl b/src/trustRegion.jl index 1c423bec3..e745159d7 100644 --- a/src/trustRegion.jl +++ b/src/trustRegion.jl @@ -3,13 +3,13 @@ TrustRegion(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(), concrete_jac = nothing, diff_type = Val{:forward}, linsolve = nothing, precs = DEFAULT_PRECS, - max_trust_radius::Number = nothing, - initial_trust_radius::Number = max_trust_radius / 11, - step_threshold::Number = 0.1, - shrink_threshold::Number = 0.25, - expand_threshold::Number = 0.75, - shrink_factor::Number = 0.25, - expand_factor::Number = 2.0, + max_trust_radius::Real = 0.0, + initial_trust_radius::Real = 0.0, + step_threshold::Real = 0.1, + shrink_threshold::Real = 0.25, + expand_threshold::Real = 0.75, + shrink_factor::Real = 0.25, + expand_factor::Real = 2.0, max_shrink_times::Int = 32) ``` @@ -50,7 +50,8 @@ for large-scale and numerically-difficult nonlinear systems. preconditioners. For more information on specifying preconditioners for LinearSolve algorithms, consult the [LinearSolve.jl documentation](https://docs.sciml.ai/LinearSolve/stable/). -- `max_trust_radius`: the maximal trust region radius. Defaults to nothing. +- `max_trust_radius`: the maximal trust region radius. + Defaults to `max(norm(fu), maximum(u) - minimum(u))`. - `initial_trust_radius`: the initial trust region radius. Defaults to `max_trust_radius / 11`. - `step_threshold`: the threshold for taking a step. In every iteration, the threshold is @@ -163,9 +164,9 @@ mutable struct TrustRegionCache{iip, fType, algType, uType, duType, resType, pTy fu_new::resType, make_new_J::Bool, r::Real) where {iip, fType, algType, uType, - duType, resType, pType, INType, - tolType, probType, ufType, L, - jType, JC} + duType, resType, pType, INType, + tolType, probType, ufType, L, + jType, JC} new{iip, fType, algType, uType, duType, resType, pType, INType, tolType, probType, ufType, L, jType, JC }(f, alg, u, fu, p, uf, linsolve, J, @@ -249,8 +250,8 @@ function perform_step!(cache::TrustRegionCache{true}) @unpack make_new_J, J, fu, f, u, p, trust_r, du1, alg, linsolve = cache if cache.make_new_J jacobian!(J, cache) - mul!(cache.H,J,J) - mul!(cache.g,J,fu) + mul!(cache.H, J, J) + mul!(cache.g, J, fu) end # u = u - J \ fu From cae5f58f5630357466e241e19970afa632efb261 Mon Sep 17 00:00:00 2001 From: Simon Carlson Date: Tue, 17 Jan 2023 19:37:34 +0100 Subject: [PATCH 11/24] fixing allocation of du1 --- src/trustRegion.jl | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/src/trustRegion.jl b/src/trustRegion.jl index e745159d7..62d751121 100644 --- a/src/trustRegion.jl +++ b/src/trustRegion.jl @@ -151,6 +151,7 @@ mutable struct TrustRegionCache{iip, fType, algType, uType, duType, resType, pTy fu_new::resType make_new_J::Bool r::Real + δN::uType function TrustRegionCache{iip}(f::fType, alg::algType, u::uType, fu::resType, p::pType, uf::ufType, linsolve::L, J::jType, du1::duType, @@ -163,7 +164,7 @@ mutable struct TrustRegionCache{iip, fType, algType, uType, duType, resType, pTy shrink_counter::Int, step_size::uType, u_new::uType, fu_new::resType, make_new_J::Bool, - r::Real) where {iip, fType, algType, uType, + r::Real, δN::uType) where {iip, fType, algType, uType, duType, resType, pType, INType, tolType, probType, ufType, L, jType, JC} @@ -175,7 +176,7 @@ mutable struct TrustRegionCache{iip, fType, algType, uType, duType, resType, pTy abstol, prob, trust_r, max_trust_r, loss, loss_new, H, g, shrink_counter, step_size, u_new, fu_new, - make_new_J, r) + make_new_J, r, δN) end end @@ -242,12 +243,12 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::TrustRegion, return TrustRegionCache{iip}(f, alg, u, fu, p, uf, linsolve, J, du1, jac_config, 1, false, maxiters, internalnorm, ReturnCode.Default, abstol, prob, initial_trust_radius, - max_trust_radius, loss, loss, J, fu, 0, u, u, fu, true, - 0.0) + max_trust_radius, loss, loss, copy(J), fu, 0, u, u, fu, true, + 0.0, u) end function perform_step!(cache::TrustRegionCache{true}) - @unpack make_new_J, J, fu, f, u, p, trust_r, du1, alg, linsolve = cache + @unpack make_new_J, J, fu, f, u, p, du1, alg, linsolve = cache if cache.make_new_J jacobian!(J, cache) mul!(cache.H, J, J) @@ -259,7 +260,8 @@ function perform_step!(cache::TrustRegionCache{true}) linu = _vec(du1), p = p, reltol = cache.abstol) cache.linsolve = linres.cache - dogleg!(-du1, trust_r, cache) + cache.δN = -du1 + dogleg!(cache) # Compute the potentially new u cache.u_new .= u .+ cache.step_size @@ -270,7 +272,7 @@ function perform_step!(cache::TrustRegionCache{true}) end function perform_step!(cache::TrustRegionCache{false}) - @unpack make_new_J, fu, f, u, p, trust_r = cache + @unpack make_new_J, fu, f, u, p = cache if make_new_J J = jacobian(cache, f) @@ -280,8 +282,8 @@ function perform_step!(cache::TrustRegionCache{false}) @unpack g, H = cache # Compute the Newton step. - δN = -H \ g - dogleg!(δN, trust_r, cache) + cache.δN = -H \ g + dogleg!(cache) # Compute the potentially new u cache.u_new = u .+ cache.step_size @@ -329,7 +331,8 @@ function trust_region_step!(cache::TrustRegionCache) end end -function dogleg!(δN, trust_r, cache::TrustRegionCache) +function dogleg!(cache::TrustRegionCache) + @unpack δN, trust_r = cache # Test if the full step is within the trust region. if norm(δN) ≤ trust_r From 65458fca18827f0bd52ae4e0b31ee3f9bf4d1d6d Mon Sep 17 00:00:00 2001 From: Simon Carlson Date: Tue, 17 Jan 2023 19:48:30 +0100 Subject: [PATCH 12/24] fixing Copy of J -> new J --- src/trustRegion.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/trustRegion.jl b/src/trustRegion.jl index 62d751121..cff85da95 100644 --- a/src/trustRegion.jl +++ b/src/trustRegion.jl @@ -239,11 +239,12 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::TrustRegion, if initial_trust_radius == 0.0 initial_trust_radius = max_trust_radius / 11 end - + H = ArrayInterfaceCore.undefmatrix(u) + return TrustRegionCache{iip}(f, alg, u, fu, p, uf, linsolve, J, du1, jac_config, 1, false, maxiters, internalnorm, ReturnCode.Default, abstol, prob, initial_trust_radius, - max_trust_radius, loss, loss, copy(J), fu, 0, u, u, fu, true, + max_trust_radius, loss, loss, H, fu, 0, u, u, fu, true, 0.0, u) end From 47ed7280719477adcfa42d4b316793b7d16e29b2 Mon Sep 17 00:00:00 2001 From: Simon Carlson Date: Tue, 17 Jan 2023 20:08:18 +0100 Subject: [PATCH 13/24] replaced $delta N$ with du1 --- src/trustRegion.jl | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/src/trustRegion.jl b/src/trustRegion.jl index cff85da95..7736706bc 100644 --- a/src/trustRegion.jl +++ b/src/trustRegion.jl @@ -151,7 +151,6 @@ mutable struct TrustRegionCache{iip, fType, algType, uType, duType, resType, pTy fu_new::resType make_new_J::Bool r::Real - δN::uType function TrustRegionCache{iip}(f::fType, alg::algType, u::uType, fu::resType, p::pType, uf::ufType, linsolve::L, J::jType, du1::duType, @@ -176,7 +175,7 @@ mutable struct TrustRegionCache{iip, fType, algType, uType, duType, resType, pTy abstol, prob, trust_r, max_trust_r, loss, loss_new, H, g, shrink_counter, step_size, u_new, fu_new, - make_new_J, r, δN) + make_new_J, r) end end @@ -240,12 +239,12 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::TrustRegion, initial_trust_radius = max_trust_radius / 11 end H = ArrayInterfaceCore.undefmatrix(u) - + return TrustRegionCache{iip}(f, alg, u, fu, p, uf, linsolve, J, du1, jac_config, 1, false, maxiters, internalnorm, ReturnCode.Default, abstol, prob, initial_trust_radius, max_trust_radius, loss, loss, H, fu, 0, u, u, fu, true, - 0.0, u) + 0.0) end function perform_step!(cache::TrustRegionCache{true}) @@ -261,7 +260,7 @@ function perform_step!(cache::TrustRegionCache{true}) linu = _vec(du1), p = p, reltol = cache.abstol) cache.linsolve = linres.cache - cache.δN = -du1 + cache.du1 .= -1 .* du1 dogleg!(cache) # Compute the potentially new u @@ -283,7 +282,7 @@ function perform_step!(cache::TrustRegionCache{false}) @unpack g, H = cache # Compute the Newton step. - cache.δN = -H \ g + cache.du1 = -H \ g dogleg!(cache) # Compute the potentially new u @@ -333,11 +332,11 @@ function trust_region_step!(cache::TrustRegionCache) end function dogleg!(cache::TrustRegionCache) - @unpack δN, trust_r = cache + @unpack du1, trust_r = cache # Test if the full step is within the trust region. - if norm(δN) ≤ trust_r - cache.step_size = δN + if norm(du1) ≤ trust_r + cache.step_size = du1 return end @@ -350,7 +349,7 @@ function dogleg!(cache::TrustRegionCache) end # Find the intersection point on the boundary. - N_sd = δN - δsd + N_sd = du1 - δsd dot_N_sd = dot(N_sd, N_sd) dot_sd_N_sd = dot(δsd, N_sd) dot_sd = dot(δsd, δsd) From 2751b1618c36e368d9d039834ede7c26af2df171 Mon Sep 17 00:00:00 2001 From: Simon Carlson Date: Tue, 17 Jan 2023 20:15:53 +0100 Subject: [PATCH 14/24] undo last commit... --- src/trustRegion.jl | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/src/trustRegion.jl b/src/trustRegion.jl index cff85da95..7736706bc 100644 --- a/src/trustRegion.jl +++ b/src/trustRegion.jl @@ -151,7 +151,6 @@ mutable struct TrustRegionCache{iip, fType, algType, uType, duType, resType, pTy fu_new::resType make_new_J::Bool r::Real - δN::uType function TrustRegionCache{iip}(f::fType, alg::algType, u::uType, fu::resType, p::pType, uf::ufType, linsolve::L, J::jType, du1::duType, @@ -176,7 +175,7 @@ mutable struct TrustRegionCache{iip, fType, algType, uType, duType, resType, pTy abstol, prob, trust_r, max_trust_r, loss, loss_new, H, g, shrink_counter, step_size, u_new, fu_new, - make_new_J, r, δN) + make_new_J, r) end end @@ -240,12 +239,12 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::TrustRegion, initial_trust_radius = max_trust_radius / 11 end H = ArrayInterfaceCore.undefmatrix(u) - + return TrustRegionCache{iip}(f, alg, u, fu, p, uf, linsolve, J, du1, jac_config, 1, false, maxiters, internalnorm, ReturnCode.Default, abstol, prob, initial_trust_radius, max_trust_radius, loss, loss, H, fu, 0, u, u, fu, true, - 0.0, u) + 0.0) end function perform_step!(cache::TrustRegionCache{true}) @@ -261,7 +260,7 @@ function perform_step!(cache::TrustRegionCache{true}) linu = _vec(du1), p = p, reltol = cache.abstol) cache.linsolve = linres.cache - cache.δN = -du1 + cache.du1 .= -1 .* du1 dogleg!(cache) # Compute the potentially new u @@ -283,7 +282,7 @@ function perform_step!(cache::TrustRegionCache{false}) @unpack g, H = cache # Compute the Newton step. - cache.δN = -H \ g + cache.du1 = -H \ g dogleg!(cache) # Compute the potentially new u @@ -333,11 +332,11 @@ function trust_region_step!(cache::TrustRegionCache) end function dogleg!(cache::TrustRegionCache) - @unpack δN, trust_r = cache + @unpack du1, trust_r = cache # Test if the full step is within the trust region. - if norm(δN) ≤ trust_r - cache.step_size = δN + if norm(du1) ≤ trust_r + cache.step_size = du1 return end @@ -350,7 +349,7 @@ function dogleg!(cache::TrustRegionCache) end # Find the intersection point on the boundary. - N_sd = δN - δsd + N_sd = du1 - δsd dot_N_sd = dot(N_sd, N_sd) dot_sd_N_sd = dot(δsd, N_sd) dot_sd = dot(δsd, δsd) From 12a26cd7bded38d6ee3c698f8e263dc277245275 Mon Sep 17 00:00:00 2001 From: Simon Carlson Date: Tue, 17 Jan 2023 20:18:50 +0100 Subject: [PATCH 15/24] undo last commit that gave error --- src/trustRegion.jl | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/trustRegion.jl b/src/trustRegion.jl index 7736706bc..ee2d25750 100644 --- a/src/trustRegion.jl +++ b/src/trustRegion.jl @@ -151,6 +151,7 @@ mutable struct TrustRegionCache{iip, fType, algType, uType, duType, resType, pTy fu_new::resType make_new_J::Bool r::Real + δN::uType function TrustRegionCache{iip}(f::fType, alg::algType, u::uType, fu::resType, p::pType, uf::ufType, linsolve::L, J::jType, du1::duType, @@ -175,7 +176,7 @@ mutable struct TrustRegionCache{iip, fType, algType, uType, duType, resType, pTy abstol, prob, trust_r, max_trust_r, loss, loss_new, H, g, shrink_counter, step_size, u_new, fu_new, - make_new_J, r) + make_new_J, r, δN) end end @@ -244,7 +245,7 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::TrustRegion, 1, false, maxiters, internalnorm, ReturnCode.Default, abstol, prob, initial_trust_radius, max_trust_radius, loss, loss, H, fu, 0, u, u, fu, true, - 0.0) + 0.0, u) end function perform_step!(cache::TrustRegionCache{true}) @@ -260,7 +261,7 @@ function perform_step!(cache::TrustRegionCache{true}) linu = _vec(du1), p = p, reltol = cache.abstol) cache.linsolve = linres.cache - cache.du1 .= -1 .* du1 + cache.δN .= -1 .* du1 dogleg!(cache) # Compute the potentially new u @@ -282,7 +283,7 @@ function perform_step!(cache::TrustRegionCache{false}) @unpack g, H = cache # Compute the Newton step. - cache.du1 = -H \ g + cache.δN = -H \ g dogleg!(cache) # Compute the potentially new u @@ -332,11 +333,11 @@ function trust_region_step!(cache::TrustRegionCache) end function dogleg!(cache::TrustRegionCache) - @unpack du1, trust_r = cache + @unpack δN, trust_r = cache # Test if the full step is within the trust region. - if norm(du1) ≤ trust_r - cache.step_size = du1 + if norm(δN) ≤ trust_r + cache.step_size = δN return end @@ -349,7 +350,7 @@ function dogleg!(cache::TrustRegionCache) end # Find the intersection point on the boundary. - N_sd = du1 - δsd + N_sd = δN - δsd dot_N_sd = dot(N_sd, N_sd) dot_sd_N_sd = dot(δsd, N_sd) dot_sd = dot(δsd, δsd) From fd6c9a3d8ffa458fcc3f4c4f25aacf9d0e745f52 Mon Sep 17 00:00:00 2001 From: Simon Carlson Date: Tue, 17 Jan 2023 20:31:33 +0100 Subject: [PATCH 16/24] replaced $delta N$ with du1 --- src/trustRegion.jl | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/src/trustRegion.jl b/src/trustRegion.jl index f3309edbe..b4abf8658 100644 --- a/src/trustRegion.jl +++ b/src/trustRegion.jl @@ -151,7 +151,6 @@ mutable struct TrustRegionCache{iip, fType, algType, uType, duType, resType, pTy fu_new::resType make_new_J::Bool r::Real - δN::uType function TrustRegionCache{iip}(f::fType, alg::algType, u::uType, fu::resType, p::pType, uf::ufType, linsolve::L, J::jType, du1::duType, @@ -164,7 +163,7 @@ mutable struct TrustRegionCache{iip, fType, algType, uType, duType, resType, pTy shrink_counter::Int, step_size::uType, u_new::uType, fu_new::resType, make_new_J::Bool, - r::Real, δN::uType) where {iip, fType, algType, uType, + r::Real) where {iip, fType, algType, uType, duType, resType, pType, INType, tolType, probType, ufType, L, jType, JC} @@ -176,7 +175,7 @@ mutable struct TrustRegionCache{iip, fType, algType, uType, duType, resType, pTy abstol, prob, trust_r, max_trust_r, loss, loss_new, H, g, shrink_counter, step_size, u_new, fu_new, - make_new_J, r, δN) + make_new_J, r) end end @@ -203,7 +202,7 @@ end function jacobian_caches(alg::TrustRegion, f, u, p, ::Val{false}) J = ArrayInterfaceCore.undefmatrix(u) - JacobianWrapper(f, p), nothing, J, nothing, nothing + JacobianWrapper(f, p), nothing, J, zero(u), nothing end function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::TrustRegion, @@ -241,12 +240,11 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::TrustRegion, end H = ArrayInterfaceCore.undefmatrix(u) - return TrustRegionCache{iip}(f, alg, u, fu, p, uf, linsolve, J, du1, jac_config, 1, false, maxiters, internalnorm, ReturnCode.Default, abstol, prob, initial_trust_radius, max_trust_radius, loss, loss, H, fu, 0, u, u, fu, true, - 0.0, copy(u)) + 0.0) end function perform_step!(cache::TrustRegionCache{true}) @@ -262,7 +260,7 @@ function perform_step!(cache::TrustRegionCache{true}) linu = _vec(du1), p = p, reltol = cache.abstol) cache.linsolve = linres.cache - cache.δN .= -1 .* du1 + cache.du1 .= -1 .* du1 dogleg!(cache) # Compute the potentially new u @@ -284,7 +282,7 @@ function perform_step!(cache::TrustRegionCache{false}) @unpack g, H = cache # Compute the Newton step. - cache.δN = -H \ g + cache.du1 = -H \ g dogleg!(cache) # Compute the potentially new u @@ -334,11 +332,11 @@ function trust_region_step!(cache::TrustRegionCache) end function dogleg!(cache::TrustRegionCache) - @unpack δN, trust_r = cache + @unpack du1, trust_r = cache # Test if the full step is within the trust region. - if norm(δN) ≤ trust_r - cache.step_size = δN + if norm(du1) ≤ trust_r + cache.step_size = du1 return end @@ -351,7 +349,7 @@ function dogleg!(cache::TrustRegionCache) end # Find the intersection point on the boundary. - N_sd = δN - δsd + N_sd = du1 - δsd dot_N_sd = dot(N_sd, N_sd) dot_sd_N_sd = dot(δsd, N_sd) dot_sd = dot(δsd, δsd) From b9ab3262fe732bf9a7f91bc7168319254cc94d55 Mon Sep 17 00:00:00 2001 From: Simon Carlson Date: Tue, 17 Jan 2023 20:53:42 +0100 Subject: [PATCH 17/24] replaced du1 with u_tmp --- src/trustRegion.jl | 52 ++++++++++++++++++++++------------------------ 1 file changed, 25 insertions(+), 27 deletions(-) diff --git a/src/trustRegion.jl b/src/trustRegion.jl index b4abf8658..f176096f7 100644 --- a/src/trustRegion.jl +++ b/src/trustRegion.jl @@ -119,7 +119,7 @@ function TrustRegion(; chunk_size = Val{0}(), max_shrink_times) end -mutable struct TrustRegionCache{iip, fType, algType, uType, duType, resType, pType, +mutable struct TrustRegionCache{iip, fType, algType, uType, resType, pType, INType, tolType, probType, ufType, L, jType, JC } f::fType @@ -130,7 +130,6 @@ mutable struct TrustRegionCache{iip, fType, algType, uType, duType, resType, pTy uf::ufType linsolve::L J::jType - du1::duType jac_config::JC iter::Int force_stop::Bool @@ -147,34 +146,34 @@ mutable struct TrustRegionCache{iip, fType, algType, uType, duType, resType, pTy g::resType shrink_counter::Int step_size::uType - u_new::uType + u_tmp::uType fu_new::resType make_new_J::Bool r::Real function TrustRegionCache{iip}(f::fType, alg::algType, u::uType, fu::resType, p::pType, - uf::ufType, linsolve::L, J::jType, du1::duType, + uf::ufType, linsolve::L, J::jType, jac_config::JC, iter::Int, force_stop::Bool, maxiters::Int, internalnorm::INType, retcode::SciMLBase.ReturnCode.T, abstol::tolType, prob::probType, trust_r::Number, max_trust_r::Number, loss::Number, loss_new::Number, H::jType, g::resType, - shrink_counter::Int, step_size::uType, u_new::uType, + shrink_counter::Int, step_size::uType, u_tmp::uType, fu_new::resType, make_new_J::Bool, r::Real) where {iip, fType, algType, uType, - duType, resType, pType, INType, + resType, pType, INType, tolType, probType, ufType, L, jType, JC} - new{iip, fType, algType, uType, duType, resType, pType, + new{iip, fType, algType, uType, resType, pType, INType, tolType, probType, ufType, L, jType, JC }(f, alg, u, fu, p, uf, linsolve, J, - du1, jac_config, iter, force_stop, + jac_config, iter, force_stop, maxiters, internalnorm, retcode, abstol, prob, trust_r, max_trust_r, loss, loss_new, H, g, shrink_counter, - step_size, u_new, fu_new, + step_size, u_tmp, fu_new, make_new_J, r) end end @@ -227,7 +226,7 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::TrustRegion, end loss = get_loss(fu) - uf, linsolve, J, du1, jac_config = jacobian_caches(alg, f, u, p, Val(iip)) + uf, linsolve, J, u_tmp, jac_config = jacobian_caches(alg, f, u, p, Val(iip)) # Set default trust region radius if not specified max_trust_radius = alg.max_trust_radius @@ -240,32 +239,31 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::TrustRegion, end H = ArrayInterfaceCore.undefmatrix(u) - return TrustRegionCache{iip}(f, alg, u, fu, p, uf, linsolve, J, du1, jac_config, + return TrustRegionCache{iip}(f, alg, u, fu, p, uf, linsolve, J, jac_config, 1, false, maxiters, internalnorm, ReturnCode.Default, abstol, prob, initial_trust_radius, - max_trust_radius, loss, loss, H, fu, 0, u, u, fu, true, + max_trust_radius, loss, loss, H, fu, 0, u, u_tmp, fu, true, 0.0) end function perform_step!(cache::TrustRegionCache{true}) - @unpack make_new_J, J, fu, f, u, p, du1, alg, linsolve = cache + @unpack make_new_J, J, fu, f, u, p, u_tmp, alg, linsolve = cache if cache.make_new_J jacobian!(J, cache) mul!(cache.H, J, J) mul!(cache.g, J, fu) end - # u = u - J \ fu linres = dolinsolve(alg.precs, linsolve, A = cache.H, b = _vec(cache.g), - linu = _vec(du1), + linu = _vec(u_tmp), p = p, reltol = cache.abstol) cache.linsolve = linres.cache - cache.du1 .= -1 .* du1 + cache.u_tmp .= -1 .* u_tmp dogleg!(cache) # Compute the potentially new u - cache.u_new .= u .+ cache.step_size - f(cache.fu_new, cache.u_new, p) + cache.u_tmp .= u .+ cache.step_size + f(cache.fu_new, cache.u_tmp, p) trust_region_step!(cache) return nothing @@ -282,19 +280,19 @@ function perform_step!(cache::TrustRegionCache{false}) @unpack g, H = cache # Compute the Newton step. - cache.du1 = -H \ g + cache.u_tmp = -H \ g dogleg!(cache) # Compute the potentially new u - cache.u_new = u .+ cache.step_size - cache.fu_new = f(cache.u_new, p) + cache.u_tmp = u .+ cache.step_size + cache.fu_new = f(cache.u_tmp, p) trust_region_step!(cache) return nothing end function trust_region_step!(cache::TrustRegionCache) - @unpack fu_new, u_new, step_size, g, H, loss, alg, max_trust_r = cache + @unpack fu_new, u_tmp, step_size, g, H, loss, alg, max_trust_r = cache cache.loss_new = get_loss(fu_new) # Compute the ratio of the actual reduction to the predicted reduction. @@ -311,7 +309,7 @@ function trust_region_step!(cache::TrustRegionCache) if r > alg.step_threshold # Take the step. - cache.u = u_new + cache.u = u_tmp cache.fu = fu_new cache.loss = cache.loss_new @@ -332,11 +330,11 @@ function trust_region_step!(cache::TrustRegionCache) end function dogleg!(cache::TrustRegionCache) - @unpack du1, trust_r = cache + @unpack u_tmp, trust_r = cache # Test if the full step is within the trust region. - if norm(du1) ≤ trust_r - cache.step_size = du1 + if norm(u_tmp) ≤ trust_r + cache.step_size = u_tmp return end @@ -349,7 +347,7 @@ function dogleg!(cache::TrustRegionCache) end # Find the intersection point on the boundary. - N_sd = du1 - δsd + N_sd = u_tmp - δsd dot_N_sd = dot(N_sd, N_sd) dot_sd_N_sd = dot(δsd, N_sd) dot_sd = dot(δsd, δsd) From 6dbae4ee2f3382e8bc94cc43511d0fa442817132 Mon Sep 17 00:00:00 2001 From: Simon Carlson Date: Tue, 17 Jan 2023 21:28:02 +0100 Subject: [PATCH 18/24] solved ForwardDif problem --- src/ad.jl | 5 +++-- test/basictests.jl | 6 ++---- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/src/ad.jl b/src/ad.jl index a2b09d569..eb9ecd94e 100644 --- a/src/ad.jl +++ b/src/ad.jl @@ -25,7 +25,8 @@ end function SciMLBase.solve(prob::NonlinearProblem{<:Union{Number, StaticArraysCore.SVector}, iip, - <:Dual{T, V, P}}, alg::NewtonRaphson, + <:Dual{T, V, P}}, + alg::Union{NewtonRaphson, TrustRegion}, args...; kwargs...) where {iip, T, V, P} sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...) return SciMLBase.build_solution(prob, alg, Dual{T, V, P}(sol.u, partials), sol.resid; @@ -34,7 +35,7 @@ end function SciMLBase.solve(prob::NonlinearProblem{<:Union{Number, StaticArraysCore.SVector}, iip, <:AbstractArray{<:Dual{T, V, P}}}, - alg::NewtonRaphson, args...; kwargs...) where {iip, T, V, P} + alg::Union{NewtonRaphson, TrustRegion}, args...; kwargs...) where {iip, T, V, P} sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...) return SciMLBase.build_solution(prob, alg, Dual{T, V, P}(sol.u, partials), sol.resid; retcode = sol.retcode) diff --git a/test/basictests.jl b/test/basictests.jl index bc8a480da..7a7237b46 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -188,8 +188,7 @@ g = function (p) return sol.u[end] end -# TODO look into the the forward diff failure when p = 100.0 -for p in 1.1:0.1:99.9 +for p in 1.1:0.1:100.0 @test g(p) ≈ sqrt(p) @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) end @@ -205,8 +204,7 @@ end @test ForwardDiff.derivative(g, 3.0) ≈ 1 / (2 * sqrt(3.0)) -# TODO look into the the forward diff failure when p = 100.0 -for p in 1.1:0.1:99.9 +for p in 1.1:0.1:100.0 @test g(p) ≈ sqrt(p) @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) end From 15fa80ef5e0af76fce29512d45ff33fa976fb035 Mon Sep 17 00:00:00 2001 From: Simon Carlson Date: Wed, 18 Jan 2023 07:05:48 +0100 Subject: [PATCH 19/24] Trying to fix allocations --- src/ad.jl | 3 ++- src/trustRegion.jl | 34 +++++++++++++++++----------------- test/basictests.jl | 16 ++++++++++------ 3 files changed, 29 insertions(+), 24 deletions(-) diff --git a/src/ad.jl b/src/ad.jl index eb9ecd94e..b88f23751 100644 --- a/src/ad.jl +++ b/src/ad.jl @@ -35,7 +35,8 @@ end function SciMLBase.solve(prob::NonlinearProblem{<:Union{Number, StaticArraysCore.SVector}, iip, <:AbstractArray{<:Dual{T, V, P}}}, - alg::Union{NewtonRaphson, TrustRegion}, args...; kwargs...) where {iip, T, V, P} + alg::Union{NewtonRaphson, TrustRegion}, args...; + kwargs...) where {iip, T, V, P} sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...) return SciMLBase.build_solution(prob, alg, Dual{T, V, P}(sol.u, partials), sol.resid; retcode = sol.retcode) diff --git a/src/trustRegion.jl b/src/trustRegion.jl index f176096f7..8c0c8982e 100644 --- a/src/trustRegion.jl +++ b/src/trustRegion.jl @@ -120,7 +120,8 @@ function TrustRegion(; chunk_size = Val{0}(), end mutable struct TrustRegionCache{iip, fType, algType, uType, resType, pType, - INType, tolType, probType, ufType, L, jType, JC + INType, tolType, probType, ufType, L, jType, JC, floatType, + trustType } f::fType alg::algType @@ -138,10 +139,10 @@ mutable struct TrustRegionCache{iip, fType, algType, uType, resType, pType, retcode::SciMLBase.ReturnCode.T abstol::tolType prob::probType - trust_r::Number - max_trust_r::Number - loss::Number - loss_new::Number + trust_r::trustType + max_trust_r::trustType + loss::floatType + loss_new::floatType H::jType g::resType shrink_counter::Int @@ -149,25 +150,24 @@ mutable struct TrustRegionCache{iip, fType, algType, uType, resType, pType, u_tmp::uType fu_new::resType make_new_J::Bool - r::Real + r::floatType function TrustRegionCache{iip}(f::fType, alg::algType, u::uType, fu::resType, p::pType, uf::ufType, linsolve::L, J::jType, jac_config::JC, iter::Int, force_stop::Bool, maxiters::Int, internalnorm::INType, retcode::SciMLBase.ReturnCode.T, abstol::tolType, - prob::probType, trust_r::Number, max_trust_r::Number, - loss::Number, - loss_new::Number, H::jType, g::resType, + prob::probType, trust_r::trustType, + max_trust_r::trustType, loss::floatType, + loss_new::floatType, H::jType, g::resType, shrink_counter::Int, step_size::uType, u_tmp::uType, - fu_new::resType, - make_new_J::Bool, - r::Real) where {iip, fType, algType, uType, - resType, pType, INType, - tolType, probType, ufType, L, - jType, JC} + fu_new::resType, make_new_J::Bool, + r::floatType) where {iip, fType, algType, uType, + resType, pType, INType, + tolType, probType, ufType, L, + jType, JC, floatType, trustType} new{iip, fType, algType, uType, resType, pType, - INType, tolType, probType, ufType, L, jType, JC + INType, tolType, probType, ufType, L, jType, JC, floatType, trustType }(f, alg, u, fu, p, uf, linsolve, J, jac_config, iter, force_stop, maxiters, internalnorm, retcode, @@ -243,7 +243,7 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::TrustRegion, 1, false, maxiters, internalnorm, ReturnCode.Default, abstol, prob, initial_trust_radius, max_trust_radius, loss, loss, H, fu, 0, u, u_tmp, fu, true, - 0.0) + loss) end function perform_step!(cache::TrustRegionCache{true}) diff --git a/test/basictests.jl b/test/basictests.jl index 7a7237b46..68f394d72 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -41,9 +41,9 @@ sol = benchmark_scalar(sf, csu0) @test sol.retcode === ReturnCode.Success @test sol.u * sol.u - 2 < 1e-9 -#@test (@ballocated benchmark_immutable(ff, cu0)) < 200 -#@test (@ballocated benchmark_mutable(ff, cu0)) < 200 -#@test (@ballocated benchmark_scalar(sf, csu0)) < 400 +# @test (@ballocated benchmark_immutable(ff, cu0)) < 200 +# @test (@ballocated benchmark_mutable(ff, cu0)) < 200 +# @test (@ballocated benchmark_scalar(sf, csu0)) < 400 function benchmark_inplace(f, u0) probN = NonlinearProblem{true}(f, u0) @@ -127,19 +127,19 @@ end # --- TrustRegion tests --- function benchmark_immutable(f, u0) - probN = NonlinearProblem(f, u0) + probN = NonlinearProblem{false}(f, u0) solver = init(probN, TrustRegion(), abstol = 1e-9) sol = solve!(solver) end function benchmark_mutable(f, u0) - probN = NonlinearProblem(f, u0) + probN = NonlinearProblem{false}(f, u0) solver = init(probN, TrustRegion(), abstol = 1e-9) sol = solve!(solver) end function benchmark_scalar(f, u0) - probN = NonlinearProblem(f, u0) + probN = NonlinearProblem{false}(f, u0) sol = (solve(probN, TrustRegion())) end @@ -161,6 +161,10 @@ sol = benchmark_scalar(sf, csu0) @test sol.retcode === ReturnCode.Success @test sol.u * sol.u - 2 < 1e-9 +@test (@ballocated benchmark_immutable(ff, cu0)) < 400 +@test (@ballocated benchmark_mutable(ff, cu0)) < 400 +@test (@ballocated benchmark_scalar(sf, csu0)) < 400 + function benchmark_inplace(f, u0) probN = NonlinearProblem(f, u0) solver = init(probN, TrustRegion(), abstol = 1e-9) From c518e9b1c46af8328a4236472ecc5534ea1fdbf8 Mon Sep 17 00:00:00 2001 From: Simon Carlson Date: Wed, 18 Jan 2023 08:18:23 +0100 Subject: [PATCH 20/24] testing for allocations --- test/basictests.jl | 28 ++++++++++++++++++---------- 1 file changed, 18 insertions(+), 10 deletions(-) diff --git a/test/basictests.jl b/test/basictests.jl index 68f394d72..f1c249d13 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -30,11 +30,12 @@ function sf(u, p) u * u - 2 end const csu0 = 1.0 +u0 = [1.0, 1.0] sol = benchmark_immutable(ff, cu0) @test sol.retcode === ReturnCode.Success @test all(sol.u .* sol.u .- 2 .< 1e-9) -sol = benchmark_mutable(ff, cu0) +sol = benchmark_mutable(ff, u0) @test sol.retcode === ReturnCode.Success @test all(sol.u .* sol.u .- 2 .< 1e-9) sol = benchmark_scalar(sf, csu0) @@ -60,6 +61,11 @@ sol = benchmark_inplace(ffiip, u0) @test sol.retcode === ReturnCode.Success @test all(sol.u .* sol.u .- 2 .< 1e-9) +u0 = [1.0, 1.0] +probN = NonlinearProblem{true}(ffiip, u0) +solver = init(probN, NewtonRaphson(), abstol = 1e-9) +@test (@ballocated solve!(solver)) < 50 + # AD Tests using ForwardDiff @@ -150,23 +156,20 @@ end function sf(u, p) u * u - 2 end +u0 = [1.0, 1.0] sol = benchmark_immutable(ff, cu0) @test sol.retcode === ReturnCode.Success @test all(sol.u .* sol.u .- 2 .< 1e-9) -sol = benchmark_mutable(ff, cu0) +sol = benchmark_mutable(ff, u0) @test sol.retcode === ReturnCode.Success @test all(sol.u .* sol.u .- 2 .< 1e-9) sol = benchmark_scalar(sf, csu0) @test sol.retcode === ReturnCode.Success @test sol.u * sol.u - 2 < 1e-9 -@test (@ballocated benchmark_immutable(ff, cu0)) < 400 -@test (@ballocated benchmark_mutable(ff, cu0)) < 400 -@test (@ballocated benchmark_scalar(sf, csu0)) < 400 - function benchmark_inplace(f, u0) - probN = NonlinearProblem(f, u0) + probN = NonlinearProblem{true}(f, u0) solver = init(probN, TrustRegion(), abstol = 1e-9) sol = solve!(solver) end @@ -180,6 +183,11 @@ sol = benchmark_inplace(ffiip, u0) @test sol.retcode === ReturnCode.Success @test all(sol.u .* sol.u .- 2 .< 1e-9) +u0 = [1.0, 1.0] +probN = NonlinearProblem{true}(ffiip, u0) +solver = init(probN, TrustRegion(), abstol = 1e-9) +@test (@ballocated solve!(solver)) < 120 + # AD Tests using ForwardDiff @@ -187,7 +195,7 @@ using ForwardDiff f, u0 = (u, p) -> u .* u .- p, @SVector[1.0, 1.0] g = function (p) - probN = NonlinearProblem(f, csu0, p) + probN = NonlinearProblem{false}(f, csu0, p) sol = solve(probN, TrustRegion(), abstol = 1e-9) return sol.u[end] end @@ -201,7 +209,7 @@ end f, u0 = (u, p) -> u * u - p, 1.0 g = function (p) - probN = NonlinearProblem(f, oftype(p, u0), p) + probN = NonlinearProblem{false}(f, oftype(p, u0), p) sol = solve(probN, TrustRegion(), abstol = 1e-10) return sol.u end @@ -217,7 +225,7 @@ f = (u, p) -> p[1] * u * u - p[2] t = (p) -> [sqrt(p[2] / p[1])] p = [0.9, 50.0] gnewton = function (p) - probN = NonlinearProblem(f, 0.5, p) + probN = NonlinearProblem{false}(f, 0.5, p) sol = solve(probN, TrustRegion()) return [sol.u] end From 25de883bc2e3c4f2e55c7de8551e8d4f8fdebb1e Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 18 Jan 2023 04:21:20 -0500 Subject: [PATCH 21/24] Update NonlinearSolve.jl --- src/NonlinearSolve.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl index ceec75389..b33a9a507 100644 --- a/src/NonlinearSolve.jl +++ b/src/NonlinearSolve.jl @@ -38,12 +38,12 @@ import SnoopPrecompile SnoopPrecompile.@precompile_all_calls begin for T in (Float32, Float64) prob = NonlinearProblem{false}((u, p) -> u .* u .- p, T(0.1), T(2)) - for alg in (NewtonRaphson, TrustRegion) + for alg in (NewtonRaphson) solve(prob, alg(), abstol = T(1e-2)) end prob = NonlinearProblem{true}((du, u, p) -> du[1] = u[1] * u[1] - p[1], T[0.1], T[2]) - for alg in (NewtonRaphson, TrustRegion) + for alg in (NewtonRaphson) solve(prob, alg(), abstol = T(1e-2)) end end end From 3348ae039989ec3f51ef79f7e8bbf75447d4108b Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 18 Jan 2023 04:35:28 -0500 Subject: [PATCH 22/24] Update NonlinearSolve.jl --- src/NonlinearSolve.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl index b33a9a507..305d9e5aa 100644 --- a/src/NonlinearSolve.jl +++ b/src/NonlinearSolve.jl @@ -38,12 +38,12 @@ import SnoopPrecompile SnoopPrecompile.@precompile_all_calls begin for T in (Float32, Float64) prob = NonlinearProblem{false}((u, p) -> u .* u .- p, T(0.1), T(2)) - for alg in (NewtonRaphson) + for alg in (NewtonRaphson,) solve(prob, alg(), abstol = T(1e-2)) end prob = NonlinearProblem{true}((du, u, p) -> du[1] = u[1] * u[1] - p[1], T[0.1], T[2]) - for alg in (NewtonRaphson) + for alg in (NewtonRaphson,) solve(prob, alg(), abstol = T(1e-2)) end end end From 1758ecae59d684cf59e96411515595ace2d1972e Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 18 Jan 2023 04:54:54 -0500 Subject: [PATCH 23/24] Update NonlinearSolve.jl --- src/NonlinearSolve.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl index 305d9e5aa..f9c06c9eb 100644 --- a/src/NonlinearSolve.jl +++ b/src/NonlinearSolve.jl @@ -38,12 +38,12 @@ import SnoopPrecompile SnoopPrecompile.@precompile_all_calls begin for T in (Float32, Float64) prob = NonlinearProblem{false}((u, p) -> u .* u .- p, T(0.1), T(2)) - for alg in (NewtonRaphson,) + for alg in (NewtonRaphson,TrustRegion) solve(prob, alg(), abstol = T(1e-2)) end prob = NonlinearProblem{true}((du, u, p) -> du[1] = u[1] * u[1] - p[1], T[0.1], T[2]) - for alg in (NewtonRaphson,) + for alg in (NewtonRaphson,TrustRegion) solve(prob, alg(), abstol = T(1e-2)) end end end From 45d61951e70a7b8b4002628a18a535139db791dd Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 18 Jan 2023 05:08:47 -0500 Subject: [PATCH 24/24] Update NonlinearSolve.jl --- src/NonlinearSolve.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl index f9c06c9eb..305d9e5aa 100644 --- a/src/NonlinearSolve.jl +++ b/src/NonlinearSolve.jl @@ -38,12 +38,12 @@ import SnoopPrecompile SnoopPrecompile.@precompile_all_calls begin for T in (Float32, Float64) prob = NonlinearProblem{false}((u, p) -> u .* u .- p, T(0.1), T(2)) - for alg in (NewtonRaphson,TrustRegion) + for alg in (NewtonRaphson,) solve(prob, alg(), abstol = T(1e-2)) end prob = NonlinearProblem{true}((du, u, p) -> du[1] = u[1] * u[1] - p[1], T[0.1], T[2]) - for alg in (NewtonRaphson,TrustRegion) + for alg in (NewtonRaphson,) solve(prob, alg(), abstol = T(1e-2)) end end end