From e8a3001a3883e7cb36064a29a1cab776bfc18552 Mon Sep 17 00:00:00 2001 From: Simon Carlson Date: Wed, 4 Jan 2023 12:15:57 +0100 Subject: [PATCH 1/5] change singularity handeling, and init of J --- src/klement.jl | 29 ++++++++++++++++------------- src/utils.jl | 10 ++++++++++ 2 files changed, 26 insertions(+), 13 deletions(-) diff --git a/src/klement.jl b/src/klement.jl index 6035167..77e0fa1 100644 --- a/src/klement.jl +++ b/src/klement.jl @@ -16,7 +16,7 @@ function SciMLBase.solve(prob::NonlinearProblem, x = float(prob.u0) fₙ = f(x) T = eltype(x) - J = ArrayInterfaceCore.zeromatrix(x) + I + J = init_J(x) if SciMLBase.isinplace(prob) error("Klement currently only supports out-of-place nonlinear problems") @@ -30,8 +30,19 @@ function SciMLBase.solve(prob::NonlinearProblem, xₙ₋₁ = x fₙ₋₁ = fₙ for _ in 1:maxiters - xₙ = xₙ₋₁ - inv(J) * fₙ₋₁ + tmp = J \ fₙ₋₁ + xₙ = xₙ₋₁ - tmp fₙ = f(xₙ) + + iszero(fₙ) && + return SciMLBase.build_solution(prob, alg, xₙ, fₙ; + retcode = ReturnCode.Success) + + if isapprox(xₙ, xₙ₋₁, atol = atol, rtol = rtol) + return SciMLBase.build_solution(prob, alg, xₙ, fₙ; + retcode = ReturnCode.Success) + end + Δxₙ = xₙ - xₙ₋₁ Δfₙ = fₙ - fₙ₋₁ @@ -41,19 +52,11 @@ function SciMLBase.solve(prob::NonlinearProblem, k = (Δfₙ - J * Δxₙ) ./ denominator J += (k * Δxₙ' .* J) * J - # Prevent inverting singular matrix - if det(J) ≈ 0 - J = ArrayInterfaceCore.zeromatrix(x) + I + # Singularity test + if cond(J) > 1e9 + J = init_J(xₙ) end - iszero(fₙ) && - return SciMLBase.build_solution(prob, alg, xₙ, fₙ; - retcode = ReturnCode.Success) - - if isapprox(xₙ, xₙ₋₁, atol = atol, rtol = rtol) - return SciMLBase.build_solution(prob, alg, xₙ, fₙ; - retcode = ReturnCode.Success) - end xₙ₋₁ = xₙ fₙ₋₁ = fₙ end diff --git a/src/utils.jl b/src/utils.jl index 25f6ba3..a4d86e2 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -33,3 +33,13 @@ value_derivative(f::F, x::AbstractArray) where {F} = f(x), ForwardDiff.jacobian( value(x) = x value(x::Dual) = ForwardDiff.value(x) value(x::AbstractArray{<:Dual}) = map(ForwardDiff.value, x) + +function init_J(x) + J = ArrayInterfaceCore.zeromatrix(x) + if ismutable(x) + J[diagind(J)] .= one(eltype(x)) + else + J += I + end + return J +end From b435ae56cd63344833bcc2fb4539d5d356784851 Mon Sep 17 00:00:00 2001 From: Simon Carlson Date: Wed, 4 Jan 2023 18:57:58 +0100 Subject: [PATCH 2/5] using lu-factorization --- src/broyden.jl | 2 +- src/klement.jl | 94 ++++++++++++++++++++++++++++++++++++-------------- 2 files changed, 69 insertions(+), 27 deletions(-) diff --git a/src/broyden.jl b/src/broyden.jl index 49f5395..06f923d 100644 --- a/src/broyden.jl +++ b/src/broyden.jl @@ -16,7 +16,7 @@ function SciMLBase.solve(prob::NonlinearProblem, x = float(prob.u0) fₙ = f(x) T = eltype(x) - J⁻¹ = ArrayInterfaceCore.zeromatrix(x) + I + J⁻¹ = init_J(x) if SciMLBase.isinplace(prob) error("Broyden currently only supports out-of-place nonlinear problems") diff --git a/src/klement.jl b/src/klement.jl index 77e0fa1..1209e95 100644 --- a/src/klement.jl +++ b/src/klement.jl @@ -4,7 +4,7 @@ Klement() ``` A low-overhead implementation of [Klement](https://jatm.com.br/jatm/article/view/373). -This method is non-allocating on scalar and static array problems. +This method is non-allocating on scalar problems. """ struct Klement <: AbstractSimpleNonlinearSolveAlgorithm end @@ -16,7 +16,7 @@ function SciMLBase.solve(prob::NonlinearProblem, x = float(prob.u0) fₙ = f(x) T = eltype(x) - J = init_J(x) + singular_tol = 1e-9 if SciMLBase.isinplace(prob) error("Klement currently only supports out-of-place nonlinear problems") @@ -29,36 +29,78 @@ function SciMLBase.solve(prob::NonlinearProblem, xₙ = x xₙ₋₁ = x fₙ₋₁ = fₙ - for _ in 1:maxiters - tmp = J \ fₙ₋₁ - xₙ = xₙ₋₁ - tmp - fₙ = f(xₙ) - - iszero(fₙ) && - return SciMLBase.build_solution(prob, alg, xₙ, fₙ; - retcode = ReturnCode.Success) - - if isapprox(xₙ, xₙ₋₁, atol = atol, rtol = rtol) - return SciMLBase.build_solution(prob, alg, xₙ, fₙ; - retcode = ReturnCode.Success) - end - Δxₙ = xₙ - xₙ₋₁ - Δfₙ = fₙ - fₙ₋₁ + # x is scalar + if isa(x, Number) + J = 1.0 + for _ in 1:maxiters + + xₙ = xₙ₋₁ - fₙ₋₁/J + fₙ = f(xₙ) + + iszero(fₙ) && + return SciMLBase.build_solution(prob, alg, xₙ, fₙ; + retcode = ReturnCode.Success) + + if isapprox(xₙ, xₙ₋₁, atol = atol, rtol = rtol) + return SciMLBase.build_solution(prob, alg, xₙ, fₙ; + retcode = ReturnCode.Success) + end - # Prevent division by 0 - denominator = max.(J' .^ 2 * Δxₙ .^ 2, 1e-9) + Δxₙ = xₙ - xₙ₋₁ + Δfₙ = fₙ - fₙ₋₁ - k = (Δfₙ - J * Δxₙ) ./ denominator - J += (k * Δxₙ' .* J) * J + # Prevent division by 0 + denominator = max(J ^ 2 * Δxₙ ^ 2, 1e-9) - # Singularity test - if cond(J) > 1e9 - J = init_J(xₙ) + k = (Δfₙ - J * Δxₙ) / denominator + J += (k * Δxₙ * J) * J + + # Singularity test + if J < singular_tol + J = 1.0 + end + + xₙ₋₁ = xₙ + fₙ₋₁ = fₙ end + # x is a vector + else + J = init_J(x) + F = lu(J, check = false) + for _ in 1:maxiters + tmp = F \ fₙ₋₁ + xₙ = xₙ₋₁ - tmp + fₙ = f(xₙ) + + iszero(fₙ) && + return SciMLBase.build_solution(prob, alg, xₙ, fₙ; + retcode = ReturnCode.Success) - xₙ₋₁ = xₙ - fₙ₋₁ = fₙ + if isapprox(xₙ, xₙ₋₁, atol = atol, rtol = rtol) + return SciMLBase.build_solution(prob, alg, xₙ, fₙ; + retcode = ReturnCode.Success) + end + + Δxₙ = xₙ - xₙ₋₁ + Δfₙ = fₙ - fₙ₋₁ + + # Prevent division by 0 + denominator = max.(J' .^ 2 * Δxₙ .^ 2, 1e-9) + + k = (Δfₙ - J * Δxₙ) ./ denominator + J += (k * Δxₙ' .* J) * J + F = lu(J, check = false) + + # Singularity test + if any(abs.(F.U[diagind(F.U)]) .< singular_tol) + J = init_J(xₙ) + F = lu(J, check = false) + end + + xₙ₋₁ = xₙ + fₙ₋₁ = fₙ + end end return SciMLBase.build_solution(prob, alg, xₙ, fₙ; retcode = ReturnCode.MaxIters) From b3488099cea9d4b6e75662222eb0875256d34168 Mon Sep 17 00:00:00 2001 From: Simon Carlson Date: Wed, 4 Jan 2023 19:41:03 +0100 Subject: [PATCH 3/5] formating --- src/klement.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/klement.jl b/src/klement.jl index 1209e95..4804e09 100644 --- a/src/klement.jl +++ b/src/klement.jl @@ -34,8 +34,7 @@ function SciMLBase.solve(prob::NonlinearProblem, if isa(x, Number) J = 1.0 for _ in 1:maxiters - - xₙ = xₙ₋₁ - fₙ₋₁/J + xₙ = xₙ₋₁ - fₙ₋₁ / J fₙ = f(xₙ) iszero(fₙ) && @@ -51,7 +50,7 @@ function SciMLBase.solve(prob::NonlinearProblem, Δfₙ = fₙ - fₙ₋₁ # Prevent division by 0 - denominator = max(J ^ 2 * Δxₙ ^ 2, 1e-9) + denominator = max(J^2 * Δxₙ^2, 1e-9) k = (Δfₙ - J * Δxₙ) / denominator J += (k * Δxₙ * J) * J @@ -64,7 +63,7 @@ function SciMLBase.solve(prob::NonlinearProblem, xₙ₋₁ = xₙ fₙ₋₁ = fₙ end - # x is a vector + # x is a vector else J = init_J(x) F = lu(J, check = false) From 44c42ee7ce452d5f52bb54bb6996032de02fb51c Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 4 Jan 2023 15:25:04 -0500 Subject: [PATCH 4/5] Update src/klement.jl --- src/klement.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/klement.jl b/src/klement.jl index 4804e09..04b5787 100644 --- a/src/klement.jl +++ b/src/klement.jl @@ -31,7 +31,7 @@ function SciMLBase.solve(prob::NonlinearProblem, fₙ₋₁ = fₙ # x is scalar - if isa(x, Number) + if x isa Number J = 1.0 for _ in 1:maxiters xₙ = xₙ₋₁ - fₙ₋₁ / J From 0f42658d158ac35882594a445e304d8991c00ed1 Mon Sep 17 00:00:00 2001 From: Simon Carlson Date: Wed, 4 Jan 2023 22:05:44 +0100 Subject: [PATCH 5/5] moving the rows around => 1 less lu-factorization --- src/klement.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/klement.jl b/src/klement.jl index 4804e09..1d8b8af 100644 --- a/src/klement.jl +++ b/src/klement.jl @@ -66,8 +66,15 @@ function SciMLBase.solve(prob::NonlinearProblem, # x is a vector else J = init_J(x) - F = lu(J, check = false) for _ in 1:maxiters + F = lu(J, check = false) + + # Singularity test + if any(abs.(F.U[diagind(F.U)]) .< singular_tol) + J = init_J(xₙ) + F = lu(J, check = false) + end + tmp = F \ fₙ₋₁ xₙ = xₙ₋₁ - tmp fₙ = f(xₙ) @@ -89,13 +96,6 @@ function SciMLBase.solve(prob::NonlinearProblem, k = (Δfₙ - J * Δxₙ) ./ denominator J += (k * Δxₙ' .* J) * J - F = lu(J, check = false) - - # Singularity test - if any(abs.(F.U[diagind(F.U)]) .< singular_tol) - J = init_J(xₙ) - F = lu(J, check = false) - end xₙ₋₁ = xₙ fₙ₋₁ = fₙ