R devel - Aug 2009 - Inaccurate complex arithmetic of R (Matlab is accurate)

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

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