> 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 > ______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel