commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From "Phil Steitz" <>
Subject Re: [math] Priority
Date Fri, 30 May 2003 16:34:16 GMT
Brent Worden wrote:
> "J.Pietschmann" <> wrote in message
>>Phil Steitz wrote:
>>>Can you provide a math reference desribing this?
>>H.M.Antia: Num. Methods for Scientists and Engineers.
>>> I have been referring
>>>to Atkinson and I am having a hard time following the implementation.
>>>What exactly do you mean by "only inverse quadratic interpolation"?
>>Brent's method requires a bracketed root as a start. The primary
>>method for shrinking the bracket is inverse quadratic interpolation.
>>This means you get three points
>>  x0,y0 x1,y1 x2,y2 x0<x2<x1
>>and interpolate a parabola
>>  x=a*y^2+b*y+c
>>and because your goal is to find the x for y=0, your next estimate
>>for the root is c (set y=0). The convergence is of the order 1.8,
>>which is better than the secant method (~1.4 if non bracketing, 1..1.2
>>on average if bracketing).
>>The full Brent algorithm measures how well the interval shrinks, and
>>resorts to bisection if progress is too slow. I didn't implement this
>>because I forgot to take the book with me and had only a hazy memory
>>of the control parameters. All in all the algorithm combines
>>near-guaranteed convergence (occasionally problems might falsely
>>detected as ill-conditioned) with near Newton-speed.
>>>What do you have in mind re: plausibility?
>>If the interval to search is (x0,x1), then absAccuracy should be of
>>the order of max(abs(x0),abs(x1))*relAccuracy. Achievable relative
>>accuracy depends on underlying hardware, although nowadays basically
>>everyone uses IEEE 754 (means, uh, 52 bit for double? Damn, have to
>>look it up!). Both relative and absolute accuracy of function
>>evaluation are also important, which is the reason to let the user
>>override defaults.
>>This means if relAcc is given then reject absAcc if
>>   max(abs(x0),abs(x1))*relAcc+c*absAcc == max(abs(x0),abs(x1))*relAcc
>>for some predermined constant c in 0.2..1.
>>>I guess I like your rootfinding framework better than Brent's (Here I
>>>mean Brent Worden, the famous commons-math contributor, not the
>>>numerical analyst).  I suggest that we agree to use your interfaces and
>>>UnivariateRealSolverImpl base,include Brent's bisection implementation
>>>strategy and modify his CDF inversion to use the new rootfinding
> framework.
>>No argument against :-) I took special care for the JavaDocs, which
>>seems to pay off...
> I agree.  The this looks like a very solid framework.  One suggestion I
> would like to make, is instead of a both a firstDirevative and
> secondDerivate method to evaluate the derivates.  Create a single
> getDerivate() method that returns a UnivariateRealFunction.  That way if a
> user needs the n-th derivate, they just call the getDerivate() method n
> times, once on the original function and once on each of the returned
> functions.  That way, common-math supports creating whatever degree derivate
> a
> method might need without ever having to change the framework.  We provide
> the maximum amout of derivate creation support to the user while providing
> us with the minimual amount of maintenance.

Mathematically, that is certainly a natural suggestion.  I am not sure, 
however, how easy it will be to implement in UnivariateRealFunction 
implementations (contract will be tricky regarding existence) or how 
often it will actually be used. Given that the rootfinding solutions 
that we seem to be converging on (pun intended) for the first release do 
not require derivative computations,
one thing we might consider would be to remove derivatives entirely from 
the UnivariateRealFunction interface for now.

It may be more natural for an R->R function to mark itself (infinitely?) 
differentiable by implementing a DifferentiableUnivariateRealFunction 
(getting a bit long, but you get the idea) interface that extends 
UnivariateRealFunction than to not throw exceptions or return null or 
"NaF" (not a function, of course) when you ask it for its derivative.
>>>I do have a few small questions/comments on the framework:
>>>1. Having solve() *require* brackets makes it impossible (or at least,
>>>un-natural) to add Newton's method or any other method that just starts
>>>with one initial value.
>>Why? Start with -4000000,+4000000000 or so. The algorithm is not
>>*required* to start with a bracket, just with an interval. Individual
>>solver implementations should document whether they require a bracket.
>>Apart from this, there is always the possiblity to add another method.
> Why not allow the user to supply any number of initial values and design the
> solvers to compensate as best they can when not enough values are provided.
> Each algorithm has criteria about the initial values that must be met before
> root finding can be attempted.  If the user provided initial values do not
> satisfy the criteria, the solver should first attempt to morph the initial
> values into a set of values that do satisfy the criteria.  If this can not
> be accomplish, then the solver would raise an exception.  This would take
> away some of the user's responsibility of knowing how these algorithms
> actually work and place it on us to create more robust components.
>>>2. I am curious as to why the "result" property is there.  How exactly
>>>do we expect clients to make use of this?
>>The client can ask for the last result as long as no further attempt
>>to solve for another root was made. I found this handy. I don't insist
>>on keeping it.
>>>3. What were you thinking about putting into the base solve()
>>>implementation as diagnostics?  Do you mean things like checking to make
>>>sure the bracket is good?\
>>No, just adding a message indicating that an unimplemented method
>>was called. Some frameworks don't display the type of the exception
>>to the user and rely on the message.
>>>4. Thinking of the developers who will use this stuff, I feel like we
>>>need a way to insulate them from having to think about rootfinding
>>>strategy choice.  As Al has pointed out, we are not (at least yet ;-))
>>>in the AI business, so we can't automagically make the best choice for
>>>them; but it might be a good idea to somehow specify a default strategy.
> A simple and flexible way to make the "best" choice for the user, is to
> control the creation of solver instances.
> This could be done via a factory or factory method.  The creation could be
> based on the class of function the user want to evaluate (real, polynomial,
> rational, etc.)  Also, we could create a solver chain that provides the
> ability to link various solvers together so if a solver fails to find a
> root, the next solver in the chain is given a chance to solve it.  Either,
> eventually one of the solvers finds a root or all solvers fail and an
> exception is raised.  The complexity of creating the "best chain" would be
> contained in the factory or factory method, totally hidden from the user.
> This again would relieve users of having some numerical method knowledge
> about root finding methods.

Nice idea; but do we really need to do this?  See below.

>>>    Unfortunately, the "safe" choice would likely be bisection.
>>Brent's algorithm is quite safe in this respect.
> I'm by no means beholden to bisection.  I just used it because I needed
> something to use for the distribution work.  Brent's algorithm is
> guarranteed to converge to a root, provided one is bracketed by the initial
> values.  So, if a bracket can be provided or generated, it is just as good
> as bisection as far as safety and it should almost always beat it in speed
> of convergence.

I agree. I suggest that we make the necessary adjustments, test the heck 
out of it, beg Al, J. and whatever other numerical analysts we can find 
to review it, and annoint Richard Brent's algorithm as our Solver and, 
while leaving the extensibility in the framework, we do not add any 
other strategies to the initial release.

>>I'll see whether I can complete the implementation near term.
> Unfortunately
>>I'll go on vacation on saturday and will be offline until  2003-06-09.
> Brent Worden
> ---------------------------------------------------------------------
> To unsubscribe, e-mail:
> For additional commands, e-mail:

To unsubscribe, e-mail:
For additional commands, e-mail:

View raw message