[Rd] Inaccurate complex arithmetic of R (Matlab is accurate)
Ravi Varadhan
rvaradhan at jhmi.edu
Wed Aug 5 05:20:14 CEST 2009
Dear John,
Thank you.
I am amazed at the insightful responses of people in the R-devel group and even more so amazed at the rapdity with which the R core group and other knowledgeable R folks respond to and fix the problems.
Big thanks to Martin Becker, Ash Richardson, Duncan Murdoch and Martin Maechler!
My frustration was only due to my impatience in having to wait a few days (!) to experiment with this marvellous idea of analytic continuation for computing numerical derivatives that are "exact", but with only minimal computation.
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: rvaradhan at jhmi.edu
----- Original Message -----
From: John Nolan <jpnolan at american.edu>
Date: Tuesday, August 4, 2009 8:36 pm
Subject: Re: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate)
To: Ravi Varadhan <rvaradhan at jhmi.edu>
Cc: 'Martin Becker' <martin.becker at mx.uni-saarland.de>, hwborchers at googlemail.com, r-devel at stat.math.ethz.ch
> Ravi,
>
> There has been a lot of chatter about this, and people don't seem to
> be
> reading carefully. Perhaps this will help clarify things.
>
> The problem appears to be that R was evaluating x^2 not as multiplication
> x*x, but as x^2.0=exp(2.0*log(x)), using standard C functions for the
> complex log and exponentiation. Apparently the C library functions for
> this have some inaccuracies in it, which showed up in your calculations.
>
> One of the other responders said that matlab, octave, etc. detect the
> case
> when there is a positive integer power and use multiplication to evaluate.
> It appears that Martin Maechler has now submitted a change to R to detect
> integer powers and evaluate them via repeated multiplication, which should
> eliminate your problem.
>
> There is no guarantee that R or matlab or any other program will evaluate
> every expression correctly. matlab has many years of use and evaluation
> that have guided them in correct problems. R is catching up, but not
> there
> yet. We are all part of improving R; your question led to a careful
> examination of the issue and a fix within a few days. No commercial
> company responds so quickly!
>
> HTH, John
>
> ...........................................................................
>
> John P. Nolan
> Math/Stat Department
> 227 Gray Hall
> American University
> 4400 Massachusetts Avenue, NW
> Washington, DC 20016-8050
>
> jpnolan at american.edu
> 202.885.3140 voice
> 202.885.3155 fax
>
> ...........................................................................
>
>
>
> -----r-devel-bounces at r-project.org wrote: -----
>
>
> To: "'Martin Becker'" <martin.becker at mx.uni-saarland.de>
> From: "Ravi Varadhan" <RVaradhan at jhmi.edu>
> Sent by: r-devel-bounces at r-project.org
> Date: 08/04/2009 10:59AM
> cc: hwborchers at googlemail.com, r-devel at 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: rvaradhan at jhmi.edu
>
> Webpage:
>
>
> tml
>
>
>
> ----------------------------------------------------------------------------
>
> --------
>
>
> -----Original Message-----
> From: Martin Becker [
> Sent: Tuesday, August 04, 2009 7:34 AM
> To: Ravi Varadhan
> Cc: r-devel at stat.math.ethz.ch; hwborchers at 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:
> , 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: rvaradhan at jhmi.edu
> >
> > Webpage:
> >
> > dhan.h
> > tml
> >
> >
> >
> > ----------------------------------------------------------------------
> > ------
> > --------
> >
> >
> > -----Original Message-----
> > From: Martin Becker [
> > Sent: Monday, August 03, 2009 5:50 AM
> > To: Ravi Varadhan
> > Cc: r-devel at stat.math.ethz.ch; hwborchers at 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: rvaradhan at jhmi.edu
> >>
> >> ______________________________________________
> >> R-devel at r-project.org mailing list
> >>
> >>
> >>
> >
> >
> > --
> > 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 at r-project.org mailing list
More information about the R-devel
mailing list