Inaccurate complex arithmetic of R (Matlab is accurate)
On 04/08/2009 1:41 PM, Richardson, Ash wrote:
Dear Ravi, Here's an idea: if R uses the C++ library complex.h, then the problem may lie in the C++ library. (Or how the C++ library is being called: perhaps it's defaulting to single precision floating point). I recently wrote an eigen-solver for a complex matrix using explicit formula; using the C++ library I got residuals of 10^-5. Using a complex number class I wrote myself using the polar representation, the same program got residuals of 10^-15. The complex functions in standard C don't seem to have the same problem. I think there is something wrong with the C++ implementation of complex numbers. I think matlab uses the standard C implementation of complex numbers, so it doesn't suffer from this problem.
R uses the C runtime, not C++. But the cpow function was not very good.
Blame Microsoft for that on Windows, some other scapegoat on whatever
other platform you use.
Martin M. has now committed some changes to R-devel that give these
results on Windows:
> #-- 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)
[1] 190.3080
> Im (rx)
[1] -4.667764e-15
These seem to agree with the desired numbers.
Duncan Murdoch
Best regards, ~ Ash
__________________________________________________ Ash Richardson, BSc Physical Scientist / Sp?cialiste des science physiques Natural Resources Canada / Ressources naturelles Canada Canadian Forest Service / Service canadien des for?ts Pacific Forestry Centre / Centre de foresterie du Pacifique 506 W. Burnside Road / 506 rue Burnside Ouest Victoria, BC V8Z 1M5 / Victoria, C-B V8Z 1M5 Tel: (250) 363-6018 Facs: (250) 363-0775 mailto:ashricha at nrcan.gc.ca ___________________________________________________ -----Original Message----- From: r-devel-bounces at r-project.org on behalf of Ravi Varadhan Sent: Tue 8/4/2009 7:59 AM To: 'Martin Becker' 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: http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h tml ---------------------------------------------------------------------------- -------- -----Original Message----- From: Martin Becker [mailto:martin.becker at mx.uni-saarland.de] 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: http://r-forge.r-project.org/projects/rmpfr/ , 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: 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 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 -- Dr. Martin Becker Statistics and Econometrics Saarland University Campus C3 1, Room 206 66123 Saarbruecken Germany ______________________________________________ R-devel at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel ______________________________________________ R-devel at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel