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:maech...@stat.math.ethz.ch] Sent: Wednesday, August 05, 2009 3:35 AM To: John Nolan Cc: Ravi Varadhan; hwborch...@googlemail.com; r-de...@stat.math.ethz.ch Subject: Re: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate) >>>>> "JN" == John Nolan <jpno...@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> jpno...@american.edu JN> 202.885.3140 voice JN> 202.885.3155 fax JN> http://academic2.american.edu/~jpnolan JN> ........................................................................... JN> -----r-devel-boun...@r-project.org wrote: ----- JN> To: "'Martin Becker'" <martin.bec...@mx.uni-saarland.de> JN> From: "Ravi Varadhan" <rvarad...@jhmi.edu> JN> Sent by: r-devel-boun...@r-project.org JN> Date: 08/04/2009 10:59AM JN> cc: hwborch...@googlemail.com, r-de...@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: rvarad...@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.bec...@mx.uni-saarland.de] JN> Sent: Tuesday, August 04, 2009 7:34 AM JN> To: Ravi Varadhan JN> Cc: r-de...@stat.math.ethz.ch; hwborch...@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: 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 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: 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 >> >> 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@r-project.org mailing list JN> https://stat.ethz.ch/mailman/listinfo/r-devel JN> [[alternative HTML version deleted]] JN> ______________________________________________ JN> R-devel@r-project.org mailing list JN> https://stat.ethz.ch/mailman/listinfo/r-devel ______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
