[Rd] Inaccurate complex arithmetic of R (Matlab is accurate)

[Hans W Borchers]

> I suspect that, in general, you may be facing the limitations of machine 
> accuracy (more precisely, IEEE 754 arithmetics on [64-bit] doubles) in 

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
>