[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

    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!

    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