root: add Brent's method for finding function roots

[fumin]

This change also introduces a new package: root.
The root package is intended to contain algorithms for function root finding.
Aside from Brent's method, another candidate algorithm to be included in the root package is Newton's method.

For gonum/gonum#483.

Please take a look.

[soypat]

Please consider adding the solver as a method on a type. I suggest using the following interface:

type RootFinder interface {
    // x0 is starting x. f is the function we seek root for.
    Root(x0 float64, f func(float64) float64) (float64, error)
}

You may find an adaptive newton's method solver I've used over the years here: https://github.com/soypat/glgl/blob/aae7eebe6ecdcdde211850267709f8e65e4ac642/math/md1/md1.go#L54

I've succesfully used it in many different applications thanks to it's flexible nature.

[fumin]

@soypat the RootFinder interface you suggested seems to support only gradient based methods, which require only a single starting x.
However, for bracket based methods including Brent, the starting object is an interval.
Moreover, among bracket based methods, there are some such as the Secant method that can work with any interval, while others such as Brent require the interval to be a "straddle" ([a, b] where f(a)*f(b) < 0).
Taking into account these possibilities, the RootFinder interface might look like:

type RootFinder interface {
    // RootGradient finds a root using a gradient based method.
    // x0 is starting x. f is the function we seek root for.
    RootGradient(x0 float64, f func(float64) float64) (float64, error)

    // RootBracket finds a root using a bracket based method.
    // [a, b] is the starting interval.
    RootBracket(a, b float64, f func(float64)float64) (float64, error)

    // RootStraddle finds a root using a straddle based method.
    // [a, b] is the starting interval.
    RootStraddle(a, b float64, f func(float64)float64) (float64, error)
}

I feel that downstream users might find this interface difficult to use, due to the bifurcation of parameters and their implications.

Nonetheless, thanks for sharing your implementation of Newton's method.
That should help a lot in making the package feature complete, in the sense of having at least one gradient based method.

[kortschak]

I'm not convinced we need to explicitly satisfy an interface at this stage. Providing the primitive seems like a good start to me. If we need to having interface satisfaction to help that, then that should be added, and if in the future there is a need to allow polymorphic behaviour, that can be added based on these primitives, either here or by client code.

[sbinet]

The package doc says it's finding roots but the package name is 'root'.
Perhaps it's a clue the package should be named 'roots'?

FYI, "my" package is also named 'roots':
https://pkg.go.dev/go-hep.org/x/exp/roots

(Apologies for the bikeshedding)

[soypat]

I agree with Dan that we can worry about designing interfaces in the future. this is experimental after all.

That said I would like to note that the NewtonRaphson root finder is also implements a bracket mode, but the bracket is part of the type :)

[kortschak]

I'm not saying things shouldn't be encapsulated into types, just that the interface satisfaction can happen later.

[fumin]

@sbinet Thanks for chiming in on a better package name.

After a survey of the standard library, I think our situation is more akin to the encoding, hash, and signal packages.
Similar to the way our package provides root-finding algorithms, these packages provide a variety of encoding and hashing methods.
On the other hand, packages that have pluralistic names such as "bytes", "slices", and "maps" tend to refer to truly enumerable collections, which I think our root package doesn't.
That said, I am happy to follow our community's decision with respect to naming.

On the topic on interfaces, perhaps we can defer designing those after we have Newton's method added into the root package?
Without a concrete implementation to tinker with, it might be a bit hard to come up with the right abstraction.

[sbinet]

BTW, I knew this discussion rang a very distant bell...
and indeed :)

gonum/gonum#896

[fumin]

BTW, I knew this discussion rang a very distant bell... and indeed :)

Yes, I am aware of gonum/gonum#896 , and found it a bit of a pity that it wasn't merged.
I can understand it's not easy to design a good interface for a generic root-finder that incorporates gradient-based, bracket-based, and even multivariate functions.

Hopefully, by starting from here in the exp repo, we can achieve the above goal by iterating bit by bit with more freedom to try out different designs.

[fumin]

I think all comments have been addressed, I wonder could anyone give this pull request another look? Thanks.

[soypat]

LGTM! We can add root functions to the package and eventually we'll find a common pattern for an interface, if that ever makes sense.

[kortschak]

This is outside my area of expertise, so leaving to others, but for light review.

[fumin]

@sbinet I wonder if you can take another look at this change? All previous comments have been addressed.

[soypat]

Is there anything blocking the merge?

[sbinet]

I'll try to carve out some time before next week to review the PR.

(apologies for the belated answer.)

[sbinet]

apologies for the belated answer.

here are a couple of comments.
@kortschak (and others, @chewxy ?) may want to chime in.

[fumin]

@sbinet thanks for the review, all comments have been addressed and incorporated.

[sbinet]

last nitpicks from my end.

[fumin]

@sbinet thanks for the comments, they have been incorporated.

[sbinet]

LGTM

Thanks for sticking with us.
I'll wait for a couple of days to let others chime in.

(and feel free to ping us if we let the ball fall on the floor)

[fumin]

@sbinet thanks for the review, it's great to have feedback from you guys experienced in numerical computing.

[kortschak]

Thanks. I'll leave to @sbinet to merge.

[sbinet]

Still LGTM.

@fumin could you send another PR against the main gonum/gonum repository, adding yourself to the AUTHORS and/or CONTRIBUTORS files ?
(then, I'll merge this PR)

thanks again.

[fumin]

@sbinet yes of course, I think gonum/gonum#2028 does this. Thanks for the review.