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

[Martin Becker]

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


Ravi Varadhan wrote:
Dear Martin,

Thank you for this useful trick.  However, we are interested in a "general"
approach for exact derivative computation.  This approach should allow the
user to specify the mathematical function that needs to be differentiated in
"any" form that the user chooses.  So, your trick will be difficult to
implement there.  Furthermore, do we know for sure that `exponentiation' is
the only operation that results in inaccuracy?  Are there other operations
that also yield inaccurate results for complex arithmetic?  

Hans Borchers also checked the computations with other free numerical
software, such as Octave, Scilab, Euler, and they all return exactly the
same results as Matlab.  It would be a shame if R could not do the same.  

It would be great if the R core could address the "fundamental" issue. 

Thank you.

Best regards,
Ravi.

----------------------------------------------------------------------------
-------

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology 

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvarad...@jhmi.edu

Webpage:
http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h
tml

 

----------------------------------------------------------------------------
--------


-----Original Message-----
From: Martin Becker [mailto:martin.bec...@mx.uni-saarland.de] 
Sent: Monday, August 03, 2009 5:50 AM
To: Ravi Varadhan
Cc: r-de...@stat.math.ethz.ch; hwborch...@googlemail.com
Subject: Re: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate)

Dear Ravi,

the inaccuracy seems to creep in when powers are calculated. Apparently,
some quite general function is called to calculate the squares, and one can
avoid the error by reformulating the example as follows:

rosen <- 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)) }

x0 <- c(0.0094, 0.7146, 0.2179, 0.6883, 0.5757, 0.9549, 0.7136, 0.0849,
0.4147, 0.4540) h <- c(1.e-15*1i, 0, 0, 0, 0, 0, 0, 0, 0, 0) xh <- x0 + h

rx <- rosen(xh)
Re(rx)
Im (rx)


I don't know which arithmetics are involved in the application you
mentioned, but writing some auxiliary function for the calculation of x^n
when x is complex and n is (a not too large) integer may solve some of the
numerical issues. A simple version is:

powN <- function(x,n) sapply(x,function(x) prod(rep(x,n)))

The corresponding summation in 'rosen' would then read:

sum(100*powN(x1-powN(x2,2),2) + powN(1-x2,2))


HTH,

  Martin


Ravi Varadhan wrote:
  
Dear All,

Hans Borchers and I have been trying to compute "exact" derivatives in R
    
using the idea of complex-step derivatives that Hans has proposed.  This is
a really, really cool idea.  It gives "exact" derivatives with only a
minimal effort (same as that involved in computing first-order
forward-difference derivative).  
  
Unfortunately, we cannot implement this in R as the "complex arithmetic"
    
in R appears to be inaccurate.
  
Here is an example:

#-- Classical Rosenbrock function in n variables rosen <- function(x) 
{ n <- length(x)
x1 <- x[2:n]
x2 <- x[1:(n-1)]
sum(100*(x1-x2^2)^2 + (1-x2)^2)
}


x0 <- c(0.0094, 0.7146, 0.2179, 0.6883, 0.5757, 0.9549, 0.7136, 
0.0849, 0.4147, 0.4540) h <- c(1.e-15*1i, 0, 0, 0, 0, 0, 0, 0, 0, 0) 
xh <- x0 + h

rx <- rosen(xh)
Re(rx)
Im (rx)

#  rx = 190.3079796814885 - 12.13915588266717e-15 i  # incorrect 
imaginary part in R

However, the imaginary part of the above answer is inaccurate.  The
    
correct imaginary part (from Matlab) is:
  
190.3079796814886 - 4.66776376640000e-15 i  # correct imaginary part 
from Matlab

This inaccuracy is serious enough to affect the acuracy of the compex-step
    
gradient drastically.
  
Hans and I were wondering if there is a way to obtain accurate "small"
    
imaginary part for complex arithmetic.  
  
I am using Windows XP operating system.

Thanks for taking a look at this.

Best regards,
Ravi.


____________________________________________________________________

Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology School of Medicine 
Johns Hopkins University

Ph. (410) 502-2619
email: rvarad...@jhmi.edu

______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
  
    


--
Dr. Martin Becker
Statistics and Econometrics
Saarland University
Campus C3 1, Room 206
66123 Saarbruecken
Germany

  


--
Dr. Martin Becker
Statistics and Econometrics
Saarland University
Campus C3 1, Room 206
66123 Saarbruecken
Germany

______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel