Skip to content
New issue

Have a question about this project? # for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “#”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? # to your account

Implementation of a Levenberg–Marquardt algorithm #131

Merged
merged 11 commits into from
Jan 23, 2023

Conversation

Deltadahl
Copy link
Contributor

@Deltadahl Deltadahl commented Jan 22, 2023

Implementation of a Levenberg–Marquardt algorithm #119, according to the paper Improvements to the Levenberg-Marquardt algorithm for nonlinear least-squares minimization.
The solver in raphson.jl was used as a template for this solver.

@codecov
Copy link

codecov bot commented Jan 22, 2023

Codecov Report

Merging #131 (518747a) into master (37034a5) will increase coverage by 2.09%.
The diff coverage is 95.94%.

@@            Coverage Diff             @@
##           master     #131      +/-   ##
==========================================
+ Coverage   88.79%   90.89%   +2.09%     
==========================================
  Files           6        7       +1     
  Lines         348      494     +146     
==========================================
+ Hits          309      449     +140     
- Misses         39       45       +6     
Impacted Files Coverage Δ
src/NonlinearSolve.jl 60.00% <ø> (ø)
src/ad.jl 100.00% <ø> (ø)
src/levenberg.jl 95.68% <95.68%> (ø)
src/trustRegion.jl 98.40% <100.00%> (+0.09%) ⬆️

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

src/ad.jl Outdated
@@ -26,7 +26,7 @@ end
function SciMLBase.solve(prob::NonlinearProblem{<:Union{Number, StaticArraysCore.SVector},
iip,
<:Dual{T, V, P}},
alg::Union{NewtonRaphson, TrustRegion},
alg::Union{NewtonRaphson, TrustRegion, LevenbergMarquardt},
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
alg::Union{NewtonRaphson, TrustRegion, LevenbergMarquardt},
alg::AbstractNewtonAlgorithm,

src/ad.jl Outdated
@@ -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...;
alg::Union{NewtonRaphson, TrustRegion, LevenbergMarquardt},
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
alg::Union{NewtonRaphson, TrustRegion, LevenbergMarquardt},
alg::AbstractNewtonAlgorithm,

src/levenberg.jl Outdated
Comment on lines 300 to 306
cache.v = -(JᵀJ .+ λ .* DᵀD) \ mul!(cache.du_tmp, J', fu)

# Geodesic acceleration (step_size = v + a / 2).
@unpack v, alg = cache
h = alg.finite_diff_step_geodesic
f(cache.fu_tmp, u .+ h .* v, p)
cache.a = -J \ ((2 / h) .* ((cache.fu_tmp .- fu) ./ h .- mul!(cache.du_tmp, J, v)))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The two \'s should be changed to linsolve handlings.

src/levenberg.jl Outdated
@unpack u, p, λ, JᵀJ, DᵀD, J = cache

# Usual Levenberg-Marquardt step ("velocity").
cache.v = -(JᵀJ .+ λ .* DᵀD) \ mul!(cache.du_tmp, J', fu)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

-(JᵀJ .+ λ .* DᵀD) should be put into a temporary

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

mul!(cache.du_tmp, J', fu) I'd move that to a different line for clarity. Though it returns cache.du_tmp in most dispatches, it's usually best to mutate in one line and then use in the next

@Deltadahl
Copy link
Contributor Author

Edit:
I will include this algorithm in precompile_algs. I naively thought I would have the same loading time as the 13 sec stated in #128, and when I got 20 sec for all three algorithms I thought the last 7 sec was due to the newly added LevenbergMarquardt.
And when I removed both LevenbergMarquardt and TrustRegion the loading time vanished.
But I now tested with only LevenbergMarquardt and NewtonRaphson and the loading time was good. While for TrustRegion and NewtonRaphson Is was 20 sec. Therefore I'll now include LevenbergMarquardt in the precompile_algs.

@Deltadahl
Copy link
Contributor Author

Testing to remove the precompilation of LevenbergMarquardt to see if it fixes the error in the check "IntegrationTest / ModelingToolkit.jl/All/1 (pull_request)"

@Deltadahl
Copy link
Contributor Author

Think I've found the issue from #128.
Have to go now, unfortunately, but I'll be back in 12h again.

@Deltadahl Deltadahl mentioned this pull request Jan 23, 2023
Comment on lines +306 to +320
linres = dolinsolve(alg.precs, linsolve, A = cache.mat_tmp, b = _vec(cache.fu_tmp),
linu = _vec(cache.du_tmp), p = p, reltol = cache.abstol)
cache.linsolve = linres.cache
@. cache.v = -cache.du_tmp

# Geodesic acceleration (step_size = v + a / 2).
@unpack v, alg = cache
h = alg.finite_diff_step_geodesic
f(cache.fu_tmp, u .+ h .* v, p)

# The following lines do: cache.a = -J \ cache.fu_tmp
mul!(cache.du_tmp, J, v)
@. cache.fu_tmp = (2 / h) * ((cache.fu_tmp - fu) / h - cache.du_tmp)
linres = dolinsolve(alg.precs, linsolve, A = J, b = _vec(cache.fu_tmp),
linu = _vec(cache.du_tmp), p = p, reltol = cache.abstol)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We may want to split the preconditioners between these two. But that's for follow up.

@ChrisRackauckas
Copy link
Member

This looks good for a first pass. More to be done, but we should split correctness from improvements.

@ChrisRackauckas ChrisRackauckas merged commit e839519 into SciML:master Jan 23, 2023
avik-pal added a commit that referenced this pull request Nov 1, 2024
Forward Mode overloads for Least Squares Problem
# for free to join this conversation on GitHub. Already have an account? # to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants