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: rvarad...@jhmi.edu Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.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 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 ______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
