From RVaradhan at jhmi.edu  Thu Aug  6 21:04:58 2009
From: RVaradhan at jhmi.edu (Ravi Varadhan)
Date: Thu, 6 Aug 2009 15:04:58 -0400
Subject: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate)
In-Reply-To: <19065.13838.892579.558430@lynne.math.ethz.ch>
References: <f5139d6efa66.4a75bdf7@johnshopkins.edu>
	<4A76B2BE.1060301@mx.uni-saarland.de>
	<004001ca1453$929c0520$7c94100a@win.ad.jhu.edu>
	<4A781CB9.8030505@mx.uni-saarland.de>
	<007701ca1514$33169e50$7c94100a@win.ad.jhu.edu>
	<OF99A54C31.A5EDECCA-ON85257608.0079614A-85257608.0079616B@american.edu>
	<19065.13838.892579.558430@lynne.math.ethz.ch>
Message-ID: <000d01ca16c8$cb1bfdc0$6153f940$@edu>

Dear All,

I just tested the complex-step derivative method on the latest windows
binary version: R 2.10.0 pre-release.

The complex-step derivatives are "exact" to within machine epsilon, with
very little additional computational effort compared to the simple,
first-order forward difference scheme.  

Thanks very much to all, especially to, Martin Maechler.

Best,
Ravi.

-----Original Message-----
From: Martin Maechler [mailto:maechler at stat.math.ethz.ch] 
Sent: Wednesday, August 05, 2009 3:35 AM
To: John Nolan
Cc: Ravi Varadhan; hwborchers at googlemail.com; r-devel at stat.math.ethz.ch
Subject: Re: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate)

>>>>> "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



[Rd] Inaccurate complex arithmetic of R (Matlab is accurate)
Ravi Varadhan