specfunc error analysis - generated and propagated errors.

[Jonathan G. Underwood]

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