need fzero (root-finding) function

[stevengj]

Julia needs to have an analogue to Matlab's fzero function for finding roots of nonlinear functions (at least scalar functions, though a multidimensional version would be great too).

Not sure what the current best-of-breed free/open-source package is for this. MINPACK is a little long in the tooth, and MINPACK-2 is not free/open-source (it is noncommercial-use only). NL2SOL is from ACM TOMS and hence is under the restrictive ACM copyright. GSL either uses MINPACK or textbook techniques at the Numerical Recipes level of sophistication.

[JeffBezanson]

https://github.com/JeffBezanson/numal/tree/master/CHAPTER5 has some code that I believe is basically equivalent to what's in GSL, but is public domain.

[stevengj]

@JeffBezanson, I don't see why you think it is "public domain". The copyright statement in the original NUMAL code online (also reproduced in your translation) only says:

The Stichting Centrum Wiskunde & Informatica (Stichting CWI) (legal successor of Stichting Mathematisch Centrum) at Amsterdam has granted permission to Paul McJones to attach the integral NUMAL library manual to his software preservation project web page. It may be freely used. It may be copied provided that the name NUMAL and the attribution to the Stichting CWI are retained.

Not only is this not "public domain" (granting permission != relinquishing copyright ownership), it is not even unambiguously open-source/free-software. "Freely used" and "may be copied" do not grant explicit permission to modify, create derived works, or sell copies, and all of these things are prohibited by default unless explicit written permission to the contrary is granted.

This is why it is dangerous to use or rely on nonstandard licenses, and why huge amounts of "freely available" numerical code (e.g. large parts of Netlib) are basically unusable in free/open-source software (unless you and your users like to live dangerously) because their authors didn't understand copyright law well enough to grant the permission they probably(?) intended.

You might be able to get in touch with Stichting CWI and get them to grant explicit permission to distribute NUMAL under a 3-clause BSD or MIT/X11 license. But honestly, for textbook algorithms like this, it is probably easier to just rewrite them from scratch as needed, starting from the English/pseudocode description of the algorithms.

[ViralBShah]

It may just be easier to write this from scratch. I wonder if the Optim package has code for something like fzero already.

https://github.com/Johnmyleswhite/Optim.jl

Ping @johnmyleswhite

[ViralBShah]

Also, ping @homerreid.

[johnmyleswhite]

We didn't build any code like that, but I believe this code base may be useful: https://github.com/matthewclegg/univariate_opt.jl

I feel like I've seen someone else mention doing this on the mailing list at one point.

[stevengj]

Note to implementer: A lot of the recent literature on univariate root-finding seems to revolve around variations on Steffensen's method (e.g. this paper). It might be nice to switch over to one of these newer methods once one is sufficiently close to the root, although I haven't found much info about robustifying these techniques (i.e. falling back to a lower-order method if they are not close enough to the root to converge quickly).

[JeffBezanson]

Steve you're right, I incorrectly remembered the license status of that code. I believe it is all classic textbook algorithms in any case.

[timholy]

For scalar functions at the Numerical Recipes level of sophistication, I can't imagine why we wouldn't just implement something directly in Julia, rather than wrapping a library.

[stevengj]

The PyCall module now supports passing Julia functions as callbacks to Python, so one option for now is simply to call Scipy, e.g.:

julia> @pyimport scipy.optimize as so
julia> so.newton(sin, 3)
3.141592653589793

But a pure Julia fzero would be good to have eventually.

[JeffBezanson]

Very impressive!

[jtravs]

While it probably falls into what you guys call the naive category, I wrote a simple fzero function in julia (I needed it for my own purposes). You can find it at: https://gist.github.com/jtravs/5327056 along with a large number of varied tests.

It is not as modern as the paper linked to by Steven, but it is more modern than the standard Matlab and scipy fzero codes and beats them both in speed, while achieving at least as good accuracy. It is based on:

G. E. Alefeld, F. A. Potra, and Y. Shi, "Algorithm 748: enclosing zeros of continuous functions," ACM Trans. Math. Softw. 21, 327–344 (1995).

But the code is written from scratch (although I peeked at the Boost C++ implementation in math/tools for some tips). I guess it might be useful for julia using people until a bigger fish comes along and implements something better.

P.S. If any of you test this code and find errors, please let me know. Also, I guess you are too busy, but if you read the code and see what you regard as awful style, do tell me, as I'm still learning julia and trying to pick up style/idiom tips.

[ViralBShah]

I guess this should be included in Optim.jl.

[johnmyleswhite]

I'm happy to incorporate code into Optim, but I see two issues: (1) I suspect many people will start to use NLopt rather than Optim now that it exists and (2) the root-finding problem solved by fzero feels more fundamental than finding minima. This is totally subjective, but I would think fzero might even belong in Base.

[mlubin]

Root finding and optimization are really equally fundamental. Newton's method is equivalent to root finding on the gradient.

[timholy]

Nice to have, @jtravs!

I don't have strong feelings about where it ends up. Certainly, Optim seems a bit misnamed (root-finding is not optimization, though obviously related), although otherwise that seems a fine place for it to be. Optim could become a place for Julia-implementations of routines, whereas NLopt uses library code.

[jtravs]

Well, to be honest I was always confused by the way root-finding and optimization where combined in scipy - it somehow blurred two subtly different issues and I was never sure if the code I was using was designed for the job at hand. But I don't have a strong opinion either way (and also haven't yet earned the right to a strong one here :)

I also have some implementations of gradient based routines (newton-raphson, halley's method etc.), which optionally use your marvelous Calculus package for the derivatives if a function is not available. In that mode they are slower than my fzero (for simple functions, and probably much more so for expensive ones), and I guess they are less robust (but haven't tested), but they don't require a bracket, just an initial guess. On the other hand, I'm trying to write a gradient-free bracketing algorithm for fzero, but can't convince myself it is that useful/robust. Maybe using the Calculus package to just get the initial direction is enough.

[stevengj]

One basic thing to remember is that it is not possible in general to devise a completely foolproof algorithm to find roots of arbitrary functions from arbitrary starting points, because the root might be hidden in an arbitrarily small portion of the real line. You can only guarantee that you find a root if the user supplies a bracket (assuming continuity).

This is not to say that algorithms that try to find a root given only a starting guess are not useful, but trying to make it completely bulletproof is a fool's errand. The most you can ask is that you ensure that it does not loop forever or jump to Inf or NaN values. (You have to be careful about this with Newton.)

Note that Steffenson's algorithm(s) are an alternative to Newton that don't require derivatives and don't require a bracket (with the usual caveat about possible failure to converge).

[tshort]

KINSOL from the Sundials library solves scalar and multidimensional cases. The Sundials.jl package provides a low-level interface with many options as well as a high-level interface. Eventually, I'd like to use the Calculus package to generate a Jacobian from expressions.

[jtravs]

OK, I tried to package up my simple root finding functions. Package is here: https://github.com/jtravs/Zeros.jl

Pull request here: https://github.com/JuliaLang/METADATA.jl/pull/201

I won't be offended if you don't want it! Or if you think it makes more sense to for these to be merged elsewhere.

[aviks]

I certainly want it available, but not sure if it goes into Optim or a separate package?

[timholy]

Definitely worth having. An independent package seems fine. Have you considered naming it Roots? Or is Zeros more internationally-accepted (or something)?

[milktrader]

I like Roots as well, particularly given we have methods called zeros

[diegozea]

I like Roots too.

[johnmyleswhite]

If the package is called Roots, could we maybe name the function findroot?

[ViralBShah]

fzero is a well-known name, and it would be nice to keep it.

A Roots package seems quite reasonable. Eventually, I would prefer to have a slightly coarser granularity of the packages. Currently, one has to get scores of packages after an initial install of julia to get some basic functionality.

Just for reference, Matlab classifies fzero as optimization: http://www.mathworks.in/help/matlab/ref/fzero.html

[timholy]

Eventually, I would prefer to have a slightly coarser granularity of the packages. Currently, one has to get scores of packages after an initial install of julia to get some basic functionality.

We might need some metapackages.

Just for reference, Matlab classifies fzero as optimization

Indeed. But I still think they are different:

• Another algorithm that would make sense for this package would be to calculate the roots of polynomials, for which one can use eig. To me, this is clearly not optimization.
• Finding the roots of a multivariable, multi-valued function is a much more poorly-defined operation than optimization. In particular, there is no analog of "descent guaranteed to find a solution." That said, I don't know a substantially better approach than using descent-type ideas and hoping for the best.

[jtravs]

OK, thanks for the feedback. I'll rename it to Roots and try and sort the packaging. Again, I really think it is a different problem than optimisation. One reason for separating it is that sometimes it can be confusing which routine to use to minimise, or to zero a function. Having separate packages for each problem removes this ambiguity. I do however, like the idea of meta-packages. Or even of Standard Packages, which are less fundamental than the standard library, but installed by default anyway, to make it easier for people to use.

[ViralBShah]

I agree that this is different problem than optimization - I was just pointing out Matlab's organization.

Meta-packages have come up in a number of other contexts too, and IIRC, there is even an issue open about having those. We can probably already have meta packages by creating an empty package that lists a bunch of packages as dependencies. It would be nice for meta-packages to be maintained by a group of people, curating the dependencies and pushing for higher quality through docs, tests, etc.

[simonster]

Is this still relevant given the current state of Roots.jl?

[vtjnash]

probably not