This is the mail archive of the gsl-discuss@sources.redhat.com mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

Re: specfunc error analysis - generated and propagated errors.


Hi Gerard,

Firstly, thanks for your considered and detailed reply, which I've been mulling over. What I'd like to do is just make a summary of what I believe is the current model for error analysis in the library, and then discuss the 'bigger picture' issues you've raised.

So, if I understand what you're saying correctly, the current model for error analysis is:

(1) Assume exact input(s). [except for a few functions where input errors are explicitly taken into account]
(2) Assume no errors on constants used in calculations.
(3) Each "primitive" operation {+,-,/,*} between two numbers adds a _relative_ error of 2*GSL_DBL_EPSILON to the (intermediate) result. [This is what I meant by generated error]
(4) Errors at each step are propogated to the next step according to normal error analysis i.e.
(i) Addition and subtraction of A and B with absolute errors a and b respectively gives an absolute propagated error in the result of a+b
(ii) Multiplication and division of A and B with absolute errors a and b respectively gives a relative propagated error in the result of [|a/A| + |b/B|]
(5) Errors from (3) and (4) are added for each step.
(6) At the end, if the error calculated by this analysis means that the computed value is not equal to the known real value within the computed error, add a Jungman relative error factor of a few GSL_DBL_EPSILON until it is :)


I hope I've not missrepresented what you've said here, but I thought it would be useful to nail down the current prescription.

Regarding the more general points about the reasoning behind the error analysis and what it achieves, that all seems reasonable. Essentially, the model is one of forward error analysis, and your arguments about the benefit this has had in the development of the library seem unquestionable, as is your argument that the user needs to be given some idea of the merit/accuracy of the results of any given function. Given though that the current forward error analysis only really gives a worst case scenario, and doesn't take into account non-symmetric errors, or cancellation etc etc, I think it's worth at least considering some alternatives. Certainly the presence of point (6) above indicates that the current approach is somewhat subjective and inprecise in some circumstances.

I wonder if a cleaner design might be to move to a backward error analysis philosophy, whereby the functions are evaluated as they are currently, but without the forward error analysis, and the test suite is used to chacterise/establish the accuracy of the implementation by comparing known values of the functions to computed values. For backwards compatibility, a typical backward-analysed error could be returned by the _e functions. It would make for a much cleaner design of the specfun code, and also remove the runtime computational overhead of (most-times often unused) the forward error analysis. This would of course require further development of the test suite, and relies upon testing over the full range of use, and of course, that the "true" values of the functions are known for comparison and error calculation. I'm not really questioning the value of the forward error analysis during the design and implementation/coding of the functions, but am questioning the worth of it to the library user, and also the cost in terms of code complexity, maintainability, and efficiency. A backward error analysis I think could be (a) more accurate/representative and (b) more efficient (b) lend itself to cleaner, simpler code. Of course, the error may vary substantially over the range/domain of the function, but that shouldn't be too difficult to deal with.

I agree that often, if the bottle kneck in the code is the calculation of eg. bessel functions, then perhaps that's an ill thought out piece of code. Often times though, it's unavoidable. Specifically for the angular momentum code I am working on, 3j/6j and sometimes 9j symbols are calculated in often deeply nested loops, and a factor of 2 or 3 here makes a massive difference. Given that these symbols have 6 or 9 indices, often many of which are summed over in a computation, efficiency of coding becomes important in many applications. This is in fact why, even though I have transplanted the code to GSL, I'll probably still have to use my own (previous) implementation without error analysis.

Just looking through NAG and numerical recipes type equivalents, I see they don't implement any error analysis returned to the user. Many conclusions could be reached from that, though :)

I agree on your point of avoiding having more than a single version of the compiled library entirely, we'll write that off as an idea based on 1 am porridge-brain.

Anyway, I'd be really interested to know what you think, and am also keen to help with the implementation of any changes that these discussions lead to.

Goggling around a bit on these matters is interesting. Many people seem to favour proper interval arithmetic for this sort of thing, but that would be a massive headache to implement. Interestingly people always argue against forward error analysis on the basis of some statement by von Neumann who showed that it was overly pessimistic for some matrix problems, couldn't find the details on that, though.


(a) Multiplication: Looking at coupling.c, it looks like the _relative_ generated error for multplication (A*B) is (2 * GSL_DBL_EPSILON). [See lines 73 and 74 of coupling.c]. However, I then looked at legendre_poly.c, as this contains some very simple algebra. Looking at line 84, I see that my assertion must be wrong. But then, comparing line 122, which is the same algebraic expression as line 84, gives an entirely different generated error analysis. What gives?



Yes, these are inconsistent by a factor of 3. This is a code
duplication problem, and I have very little idea how extensive
it is. Rubber mallet time?


I also noticed something else about this code - although Legendre polynomials are "hard coded" up until P_3, in the general function, only up to P_2 is hard coded. Anyway, if you let me know which of the two error analyses you think is more approriate, I'll happily cook up a patch to sort it out.

What do you mean by truncation error?



er, ignore that, not sure what I was on about there.


By the way, lately I have been influenced by a nice
little book by Muller called "Elementary Functions: Algorithms and
Implementation". I wish I had seen it before I started working on
the special functions for GSL.



thanks for the reference, I'll have a look at that.

Best wishes,

Jonathan


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]