Best regards,
~ Ash
__________________________________________________
Ash Richardson, BSc
Physical Scientist / Spécialiste des science physiques
Natural Resources Canada / Ressources naturelles Canada
Canadian Forest Service / Service canadien des forêts
Pacific Forestry Centre / Centre de foresterie du Pacifique
506 W. Burnside Road / 506 rue Burnside Ouest
Victoria, BC V8Z 1M5 / Victoria, C-B V8Z 1M5
Tel: (250) 363-6018 Facs: (250) 363-0775
mailto:ashri...@nrcan.gc.ca
___________________________________________________
-----Original Message-----
From: r-devel-boun...@r-project.org on behalf of Ravi Varadhan
Sent: Tue 8/4/2009 7:59 AM
To: 'Martin Becker'
Cc: hwborch...@googlemail.com; r-de...@stat.math.ethz.ch
Subject: Re: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate)
Please forgive me for my lack of understanding of IEEE floating-point
arithmetic. I have a hard time undertsanding why "this is not a problem of
R itself", when "ALL" the other well known computing environments including
Matlab, Octave, S+, and Scilab provide accurate results. My concern is not
really about the "overall" accuracy of R, but just about the complex
arithmetic. Is there something idiosyncratic about the complex arithmetic?
I am really hoping that some one from the R core would speak up and address
this issue. It would be a loss to the R users if this fascinating idea of
"complex-step derivative" could not be implemented in R.
Thanks,
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: Tuesday, August 04, 2009 7:34 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,
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/Vara
dhan.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
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
