|
| 1 | +""" |
| 2 | + HalleyDescent(; linsolve = nothing) |
| 3 | +
|
| 4 | +Improve the NewtonDescent with higher-order terms. First compute the descent direction as ``J a = -fu``. |
| 5 | +Then compute the hessian-vector-vector product and solve for the second-order correction term as ``J b = H a a``. |
| 6 | +Finally, compute the descent direction as ``δu = a * a / (b / 2 - a)``. |
| 7 | +
|
| 8 | +Note that `import TaylorDiff` is required to use this descent algorithm. |
| 9 | +
|
| 10 | +See also [`NewtonDescent`](@ref). |
| 11 | +""" |
| 12 | +@kwdef @concrete struct HalleyDescent <: AbstractDescentDirection |
| 13 | + linsolve = nothing |
| 14 | +end |
| 15 | + |
| 16 | +supports_line_search(::HalleyDescent) = true |
| 17 | + |
| 18 | +@concrete mutable struct HalleyDescentCache <: AbstractDescentCache |
| 19 | + f |
| 20 | + p |
| 21 | + δu |
| 22 | + δus |
| 23 | + b |
| 24 | + fu |
| 25 | + hvvp |
| 26 | + lincache |
| 27 | + timer |
| 28 | + preinverted_jacobian <: Union{Val{false}, Val{true}} |
| 29 | +end |
| 30 | + |
| 31 | +@internal_caches HalleyDescentCache :lincache |
| 32 | + |
| 33 | +function InternalAPI.init( |
| 34 | + prob::NonlinearProblem, alg::HalleyDescent, J, fu, u; stats, |
| 35 | + shared = Val(1), pre_inverted::Val = Val(false), |
| 36 | + linsolve_kwargs = (;), abstol = nothing, reltol = nothing, |
| 37 | + timer = get_timer_output(), kwargs...) |
| 38 | + @bb δu = similar(u) |
| 39 | + @bb b = similar(u) |
| 40 | + @bb fu = similar(fu) |
| 41 | + @bb hvvp = similar(fu) |
| 42 | + δus = Utils.unwrap_val(shared) ≤ 1 ? nothing : map(2:Utils.unwrap_val(shared)) do i |
| 43 | + @bb δu_ = similar(u) |
| 44 | + end |
| 45 | + lincache = Utils.unwrap_val(pre_inverted) ? nothing : |
| 46 | + construct_linear_solver( |
| 47 | + alg, alg.linsolve, J, Utils.safe_vec(fu), Utils.safe_vec(u); |
| 48 | + stats, abstol, reltol, linsolve_kwargs... |
| 49 | + ) |
| 50 | + return HalleyDescentCache( |
| 51 | + prob.f, prob.p, δu, δus, b, fu, hvvp, lincache, timer, pre_inverted) |
| 52 | +end |
| 53 | + |
| 54 | +function InternalAPI.solve!( |
| 55 | + cache::HalleyDescentCache, J, fu, u, idx::Val = Val(1); |
| 56 | + skip_solve::Bool = false, new_jacobian::Bool = true, kwargs...) |
| 57 | + δu = SciMLBase.get_du(cache, idx) |
| 58 | + skip_solve && return DescentResult(; δu) |
| 59 | + if preinverted_jacobian(cache) |
| 60 | + @assert J!==nothing "`J` must be provided when `pre_inverted = Val(true)`." |
| 61 | + @bb δu = J × vec(fu) |
| 62 | + else |
| 63 | + @static_timeit cache.timer "linear solve 1" begin |
| 64 | + linres = cache.lincache(; |
| 65 | + A = J, b = Utils.safe_vec(fu), |
| 66 | + kwargs..., linu = Utils.safe_vec(δu), |
| 67 | + reuse_A_if_factorization = !new_jacobian || (idx !== Val(1))) |
| 68 | + δu = Utils.restructure(SciMLBase.get_du(cache, idx), linres.u) |
| 69 | + if !linres.success |
| 70 | + set_du!(cache, δu, idx) |
| 71 | + return DescentResult(; δu, success = false, linsolve_success = false) |
| 72 | + end |
| 73 | + end |
| 74 | + end |
| 75 | + b = cache.b |
| 76 | + # compute the hessian-vector-vector product |
| 77 | + hvvp = evaluate_hvvp(cache.hvvp, cache, cache.f, cache.p, u, δu) |
| 78 | + # second linear solve, reuse factorization if possible |
| 79 | + if preinverted_jacobian(cache) |
| 80 | + @bb b = J × vec(hvvp) |
| 81 | + else |
| 82 | + @static_timeit cache.timer "linear solve 2" begin |
| 83 | + linres = cache.lincache(; |
| 84 | + A = J, b = Utils.safe_vec(hvvp), |
| 85 | + kwargs..., linu = Utils.safe_vec(b), |
| 86 | + reuse_A_if_factorization = true) |
| 87 | + b = Utils.restructure(cache.b, linres.u) |
| 88 | + if !linres.success |
| 89 | + set_du!(cache, δu, idx) |
| 90 | + return DescentResult(; δu, success = false, linsolve_success = false) |
| 91 | + end |
| 92 | + end |
| 93 | + end |
| 94 | + @bb @. δu = δu * δu / (b / 2 - δu) |
| 95 | + set_du!(cache, δu, idx) |
| 96 | + cache.b = b |
| 97 | + return DescentResult(; δu) |
| 98 | +end |
| 99 | + |
| 100 | +evaluate_hvvp(hvvp, cache, f, p, u, δu) = error("not implemented. please import TaylorDiff") |
0 commit comments