[Rd] Inaccurate complex arithmetic of R (Matlab is accurate)
Martin Maechler
maechler at stat.math.ethz.ch
Wed Aug 5 09:34:38 CEST 2009
>>>>> "JN" == John Nolan <jpnolan at american.edu>
>>>>> on Tue, 4 Aug 2009 18:05:47 -0400 writes:
JN> Ravi,
JN> There has been a lot of chatter about this, and people don't seem to be
JN> reading carefully. Perhaps this will help clarify things.
Yes, I hope so; thanks a lot, John!!
{Indeed, I've been appalled too about the amount of proclaiming
without listening, and even more about the no-bug report... }
JN> The problem appears to be that R was evaluating x^2 not as multiplication
JN> x*x, but as x^2.0=exp(2.0*log(x)), using standard C functions for the
JN> complex log and exponentiation. Apparently the C library functions for
JN> this have some inaccuracies in it, which showed up in your calculations.
JN> One of the other responders said that matlab, octave, etc. detect the case
JN> when there is a positive integer power and use multiplication to evaluate.
JN> It appears that Martin Maechler has now submitted a change to R to detect
JN> integer powers and evaluate them via repeated multiplication, which should
JN> eliminate your problem.
Yes; note however that -- exactly for the reason John has given
(and Martin Becker and I had tried to explain earlier in this
thread!) --
that your problem (*yours* indeed, not R's !!) will resurface if
you change "^2" by "^2.01" in your original problem.
JN> There is no guarantee that R or matlab or any other
JN> program will evaluate every expression correctly.
JN> matlab has many years of use and evaluation that have
JN> guided them in correct problems. R is catching up, but
JN> not there yet.
hhmm, hmm,... are you sure that the catching up is only on our side?
I dare say that I'm pretty sure R does quite a few computations
better than e.g., matlab... But let's not get into more unproductive
chattering...
JN> We are all part of improving R; your question led to a careful
JN> examination of the issue and a fix within a few days. No commercial
JN> company responds so quickly!
JN> HTH, John
it did, thank you very much!
Martin
JN> ...........................................................................
JN> John P. Nolan
JN> Math/Stat Department
JN> 227 Gray Hall
JN> American University
JN> 4400 Massachusetts Avenue, NW
JN> Washington, DC 20016-8050
JN> jpnolan at american.edu
JN> 202.885.3140 voice
JN> 202.885.3155 fax
JN> http://academic2.american.edu/~jpnolan
JN> ...........................................................................
JN> -----r-devel-bounces at r-project.org wrote: -----
JN> To: "'Martin Becker'" <martin.becker at mx.uni-saarland.de>
JN> From: "Ravi Varadhan" <RVaradhan at jhmi.edu>
JN> Sent by: r-devel-bounces at r-project.org
JN> Date: 08/04/2009 10:59AM
JN> cc: hwborchers at googlemail.com, r-devel at stat.math.ethz.ch
JN> Subject: Re: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate)
JN> Please forgive me for my lack of understanding of IEEE floating-point
JN> arithmetic. I have a hard time undertsanding why "this is not a problem of
JN> R itself", when "ALL" the other well known computing environments including
JN> Matlab, Octave, S+, and Scilab provide accurate results. My concern is not
JN> really about the "overall" accuracy of R, but just about the complex
JN> arithmetic. Is there something idiosyncratic about the complex arithmetic?
JN> I am really hoping that some one from the R core would speak up and address
JN> this issue. It would be a loss to the R users if this fascinating idea of
JN> "complex-step derivative" could not be implemented in R.
JN> Thanks,
JN> Ravi.
JN> ----------------------------------------------------------------------------
JN> -------
JN> Ravi Varadhan, Ph.D.
JN> Assistant Professor, The Center on Aging and Health
JN> Division of Geriatric Medicine and Gerontology
JN> Johns Hopkins University
JN> Ph: (410) 502-2619
JN> Fax: (410) 614-9625
JN> Email: rvaradhan at jhmi.edu
JN> Webpage:
JN> http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h
JN> tml
JN> ----------------------------------------------------------------------------
JN> --------
JN> -----Original Message-----
JN> From: Martin Becker [mailto:martin.becker at mx.uni-saarland.de]
JN> Sent: Tuesday, August 04, 2009 7:34 AM
JN> To: Ravi Varadhan
JN> Cc: r-devel at stat.math.ethz.ch; hwborchers at googlemail.com
JN> Subject: Re: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate)
JN> Dear Ravi,
JN> I suspect that, in general, you may be facing the limitations of machine
JN> accuracy (more precisely, IEEE 754 arithmetics on [64-bit] doubles) in your
JN> application. This is not a problem of R itself, but rather a problem of
JN> standard arithmetics provided by underlying C compilers/CPUs.
JN> In fact, every operation in IEEE arithmetics (so, this is not really a
JN> problem only for complex numbers) may suffer from inexactness, a
JN> particularly difficult one is addition/subtraction. Consider the following
JN> example for real numbers (I know, it is not a very good one...):
JN> The two functions
JN> badfn <- function(x) 1-(1+x)*(1-x)
JN> goodfn <- function(x) x^2
JN> both calculate x^2 (theoretically, given perfect arithmetic). So, as you
JN> want to allow the user to 'specify the mathematical function ... in "any"
JN> form the user chooses', both functions should be ok.
JN> But, unfortunately:
>> badfn(1e-8)
JN> [1] 2.220446049250313e-16
>> goodfn(1e-8)
JN> [1] 1e-16
JN> I don't know what happens in matlab/octave/scilab for this example. They
JN> may
JN> do better, but probably at some cost (dropping IEEE arithmetic/do "clever"
JN> calculations should result in massive speed penalties, try
JN> evaluating hypergeom([1,-99.9],[-49.9-24.9*I],(1+1.71*I)/2); in
JN> Maple...).
JN> Now, you have some options:
JN> - assume, that the user is aware of the numerical inexactness of ieee
JN> arithmetics and that he is able to supply some "robust" version of the
JN> mathematical function.
JN> - use some other software (eg., matlab) for the critical calculations
JN> (there
JN> is a R <-> Matlab interface, see package R.matlab on CRAN), if you are
JN> sure,
JN> that this helps.
JN> - implement/use multiple precision arithmetics within R (Martin Maechler's
JN> Rmpfr package may be very useful:
JN> http://r-forge.r-project.org/projects/rmpfr/ , but this will slow down
JN> calculations considerably)
JN> All in all, I think it is unfair just to blame R here. Of course, it would
JN> be great if there was a simple trigger to turn on multiple precision
JN> arithmetics in R. Packages such as Rmpfr may provide a good step in this
JN> direction, since operator overloading via S4 classes allows for easy code
JN> adaption. But Rmpfr is still declared "beta", and it relies on some
JN> external
JN> library, which could be problematic on Windows systems. Maybe someone else
JN> has other/better suggestions, but I do not think that there is an easy
JN> solution for the "general" problem.
JN> Best wishes,
JN> Martin
JN> Ravi Varadhan wrote:
>> Dear Martin,
>>
>> Thank you for this useful trick. However, we are interested in a
JN> "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
JN> results
JN> 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
JN> 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:
>> http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Vara
>> dhan.h
>> tml
>>
>>
>>
>> ----------------------------------------------------------------------
>> ------
>> --------
>>
>>
>> -----Original Message-----
>> From: Martin Becker [mailto:martin.becker at mx.uni-saarland.de]
>> 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
JN> 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
>>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>>
>>>
>>
>>
>> --
>> Dr. Martin Becker
>> Statistics and Econometrics
>> Saarland University
>> Campus C3 1, Room 206
>> 66123 Saarbruecken
>> Germany
>>
>>
JN> --
JN> Dr. Martin Becker
JN> Statistics and Econometrics
JN> Saarland University
JN> Campus C3 1, Room 206
JN> 66123 Saarbruecken
JN> Germany
JN> ______________________________________________
JN> R-devel at r-project.org mailing list
JN> https://stat.ethz.ch/mailman/listinfo/r-devel
JN> [[alternative HTML version deleted]]
JN> ______________________________________________
JN> R-devel at r-project.org mailing list
JN> https://stat.ethz.ch/mailman/listinfo/r-devel
More information about the R-devel
mailing list