Implemented a Klement solver #18
Conversation
Codecov Report
@@ Coverage Diff @@
## main #18 +/- ##
===========================================
- Coverage 91.48% 33.79% -57.70%
===========================================
Files 7 8 +1
Lines 188 216 +28
===========================================
- Hits 172 73 -99
- Misses 16 143 +127
📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more |
| if det(J) ≈ 0 | ||
| J = ArrayInterfaceCore.zeromatrix(x) + I |
There was a problem hiding this comment.
This seems a little hacky. Is it part of the method?
There was a problem hiding this comment.
Ah, yeah, it's hacky...
I'll have a look at it tomorrow!
There was a problem hiding this comment.
Let's merge and improve.
| x = float(prob.u0) | ||
| fₙ = f(x) | ||
| T = eltype(x) | ||
| J = ArrayInterfaceCore.zeromatrix(x) + I |
There was a problem hiding this comment.
This could be a fallback. If x is mutable, then we could do J = ArrayInterfaceCore.zeromatrix(x); J[diagind(J)] .= one(eltype(x)).
| J += (k * Δxₙ' .* J) * J | ||
|
|
||
| # Prevent inverting singular matrix | ||
| if det(J) ≈ 0 |
There was a problem hiding this comment.
det(J) ≈ 0 is a bad check for singularity. Consider det(0.1I(1000)).
There was a problem hiding this comment.
Also, by default, x ≈ 0 is true iff x == 0.
There was a problem hiding this comment.
Yes, this was left for a follow-up: #19
| xₙ₋₁ = x | ||
| fₙ₋₁ = fₙ | ||
| for _ in 1:maxiters | ||
| xₙ = xₙ₋₁ - inv(J) * fₙ₋₁ |
There was a problem hiding this comment.
Always use the matrix factorization object to solve.
There was a problem hiding this comment.
I missed that. Instead, just \ here.
There was a problem hiding this comment.
Or F = lu(J), and the singularity check is then abs(F.factors[end, end]) < singular_tol.
Implementation of a Klement solver according to https://jatm.com.br/jatm/article/view/373
This implementation is very similar to "broyden.jl".
I did some experiments with this solver compared to the Broyden-solver,
and the Klement-solver was able to solve the problems a lot faster (order of magnitudes).
However, I haven't tried the solvers on particularly "hard" problems yet.
(I also found some improvements I think I can do to the Broyden-solver, but that will come in another PR)