Inaccurate complex arithmetic of R (Matlab is accurate)

Tue, Aug 4, 2009 7:35 AM

Dear Martin,

I definitely do not agree with this. Consider your own proposal of
writing the Rosenbrock function:

    rosen2 <- function(x) {
        n <- length(x)
        x1 <- x[2:n]
        x2 <- x[1:(n-1)]
        sum(100*(x1-x2*x2)*(x1-x2*x2) + (1-x2)*(1-x2))
    }

The complex-step derivative approximation suggests to set h <- 1.0e-20
and then to compute the expression  Im(rosen(x+hi))/h , so let's do it:

    h  <- 1.0e-20
    xh <- c(0.0094 + h*1i, 0.7146, 0.2179, 0.6883, 0.5757,
            0.9549,        0.7136, 0.0849, 0.4147, 0.4540)

    Im(rosen2(xh)) / h
    # [1] -4.6677637664000

which is exactly the correct derivative in the first variable, namely
d/dx1 rosen(x). To verify, you could even calculate this by hand on a
"Bierdeckel".

Now look at the former definition, say rosen1(), using '^' instead:

    rosen1 <- function(x) {
        n <- length(x)
        x1 <- x[2:n]
        x2 <- x[1:(n-1)]
        sum(100*(x1-x2^2)^2 + (1-x2)^2)
    }

    Im(rosen1(xh)) / h
    # [1] -747143.8793904837

We find that R is "running wild". And this has nothing to do with IEEE
arithmetics or machine accuracy, as we can see from the first example
where R is able to do it right.

And as I said, the second example is working correctly in free software
such as Octave which I guess does not do any "clever" calculations here.

We do _not_ ask for multiple precision arithmetics here !

Regards
Hans Werner


Martin Becker <martin.becker <at> mx.uni-saarland.de> writes:

Dear Ravi,

I suspect that, in general, you may be facing the limitations of machine 
accuracy (more precisely, IEEE 754 arithmetics on [64-bit] doubles) in 
your application. This is not a problem of R itself, but rather a 
problem of standard arithmetics provided by underlying C compilers/CPUs.
In fact, every operation in IEEE arithmetics (so, this is not really a 
problem only for complex numbers) may suffer from inexactness, a 
particularly difficult one is addition/subtraction. Consider the 
following example for real numbers (I know, it is not a very good one...):
The two functions

badfn <- function(x) 1-(1+x)*(1-x)
goodfn <- function(x) x^2

both calculate x^2 (theoretically, given perfect arithmetic). So, as you 
want to allow the user to 'specify the mathematical function ... in 
"any" form the user chooses', both functions should be ok.
But, unfortunately:

 > badfn(1e-8)

[1] 2.220446049250313e-16

 > goodfn(1e-8)

[1] 1e-16

I don't know what happens in matlab/octave/scilab for this example. They 
may do better, but probably at some cost (dropping IEEE arithmetic/do 
"clever" calculations should result in massive speed penalties, try 
evaluating   hypergeom([1,-99.9],[-49.9-24.9*I],(1+1.71*I)/2);   in 
Maple...).
Now, you have some options:

- assume, that the user is aware of the numerical inexactness of ieee 
arithmetics and that he is able to supply some "robust" version of the 
mathematical function.
- use some other software (eg., matlab) for the critical calculations 
(there is a R <-> Matlab interface, see package R.matlab on CRAN), if 
you are sure, that this helps.
- implement/use multiple precision arithmetics within R (Martin 
Maechler's Rmpfr package may be very useful: 
http://r-forge.r-project.org/projects/rmpfr/ , but this will slow down 
calculations considerably)

All in all, I think it is unfair just to blame R here. Of course, it 
would be great if there was a simple trigger to turn on multiple 
precision arithmetics in R. Packages such as Rmpfr may provide a good 
step in this direction, since operator overloading via S4 classes allows 
for easy code adaption. But Rmpfr is still declared "beta", and it 
relies on some external library, which could be problematic on Windows 
systems. Maybe someone else has other/better suggestions, but I do not 
think that there is an easy solution for the "general" problem.

Best wishes,

  Martin

Inaccurate complex arithmetic of R (Matlab is accurate)

Thread (16 messages)