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

[Hans W. Borchers]

>

Thanks for pointing out the weak point in this computation. I tried out your
suggestions and they both deliver the correct and accurate result.

But as a general solution this approach is not feasible. We want to provide
"complex-step derivatives" as a new method for computing exact gradients, for
example in 'numDeriv::grad' as

    grad(fun, x, method="complex-step")

and we can not reasonably assume that a user prepares his function specially 
for calling this method.

I tried out other numerical math software besides Matlab, that is Octave,
Scilab, Euler and others, and they all return the same correct value up to 15
digits. Should we not expect that R is capable of doing the same?

Hans W. Borchers


Martin Becker <martin.becker <at> mx.uni-saarland.de> writes:

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