-
-
Notifications
You must be signed in to change notification settings - Fork 86
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
New Root Solver: ITP Method #544
base: main
Are you sure you want to change the base?
Conversation
max, | ||
tol, | ||
// kappa1, suggested from paper | ||
float!(0.2) / (max - min), |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Wondering, maybe we want to have this function return Result
and verify max - min
cannot be zero so we don't panic here. Not sure though, maybe not worth it right now.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This sounds like a good addition!
let sol = (self.a + self.b) * float!(0.5); | ||
// TODO: This function evaluation serves no purpose other than to serve argmin's cost | ||
// method on the state. It feels wasteful. | ||
let f_sol = problem.cost(&sol)?; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This bothers me, but maybe someone smarter than me might know a clever way to not do this.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I agree. Do I understand correctly that the algorithm itself does not require to compute the final cost function value? It is only required here because otherwise the cost function value would not be related to the final solution?
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #544 +/- ##
==========================================
+ Coverage 91.92% 91.94% +0.01%
==========================================
Files 177 178 +1
Lines 23724 23940 +216
==========================================
+ Hits 21808 22011 +203
- Misses 1916 1929 +13 ☔ View full report in Codecov by Sentry. |
@stefan-k Does this seem like a worthwhile addition? If so, I did take a peek at CI here and some of these seem like flakes- I could totally be misreading, though. |
Hi @duncanam , this is definitely a worthwhile addition, thank you! Unfortunately it will take a bit until I can give a useful review, hopefully around the upcoming holidays. Apologies for being so unresponsive recently :/ |
No worries! Have a good holiday season! |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thank you for this PR! This is a very valuable addition and looks pretty good to me. :) I only have a few minor points to discuss. Regarding the math I'll have to trust you here as I currently am not able to dive into this in detail.
if self.tol < F::zero() { | ||
return Err(ItpRootError::NegativeTol.into()); | ||
} | ||
// This helps ensure the log evaluation is stable | ||
if self.a > self.b { | ||
return Err(ItpRootError::MinLargerThanMax.into()); | ||
} | ||
// It's important to check this to verify n1o2 doesn't panic | ||
if self.tol.is_zero() { | ||
return Err(ItpRootError::ZeroTol.into()); | ||
} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I know that BrentRoot
also checks these fields in the init
method of Solver
, but I think it may be better to do these checks in the new
method of ItpRoot
. What do you think?
tol, | ||
kappa1, | ||
kappa2, | ||
n0, | ||
a: min, | ||
b: max, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Question out of curiosity: Are there any sane defaults for these parameters (ideally something mentioned in the paper)?
let sol = (self.a + self.b) * float!(0.5); | ||
// TODO: This function evaluation serves no purpose other than to serve argmin's cost | ||
// method on the state. It feels wasteful. | ||
let f_sol = problem.cost(&sol)?; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I agree. Do I understand correctly that the algorithm itself does not require to compute the final cost function value? It is only required here because otherwise the cost function value would not be related to the final solution?
I'll try to fix the CI problems in another PR. |
Intro
Hello! I noticed that this lovely crate does not have an implementation of the ITP method root solver, which is, quoting Wikipedia:
"The ITP method, short for Interpolate Truncate and Project, is the first root-finding algorithm that achieves the superlinear convergence of the secant method while retaining the optimal worst-case performance of the bisection method. It is also the first method with guaranteed average performance strictly better than the bisection method under any continuous distribution. In practice it performs better than traditional interpolation and hybrid based strategies (Brent's Method, Ridders, Illinois), since it not only converges super-linearly over well behaved functions but also guarantees fast performance under ill-behaved functions where interpolations fail."
Practically, it is faster than
BrentRoot
due to fewer function evaluations. I thought this might be a handy contribution towardsargmin
, since I didn't see too many implementations of this online (There's theitp
package in R, and in Julia I think it exists inRoots.jl
. I also found it located within the Rust cratekurbo
here, but this was tuned specifically for curve fitting instead of generic root solving. I templated this implementation based off ofbrentroot.rs
.Verification
I used the quadratic test from
brentroot.rs
, and it achieved a solution in one iteration faster. However, there's an additional function evaluation that occurs to satisfy theargmin
iteration state, so currently comes in equal in terms of function evaluations with BrentRoot for that example specifically.Additionally, I added a test against the polynomial shown in the example on Wikipedia, and stepping through in a debugger I was able to verify that
a
andb
bounds on the bracket (of the solver) indeed match what is in that table there.Let me know what else should be changed, as I'm new to this community. Thanks! 😃