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

[Martin Maechler]

>>>>> "HWB" == Hans W Borchers <hwborchers at googlemail.com>
>>>>>     on Mon, 3 Aug 2009 13:15:11 +0000 (UTC) writes:

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

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

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

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

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

R's  "a ^ b"  on non-Windows typically uses the C library's
'cpow()'  when that is provided  (HAVE_C99_COMPLEX).

Indeed, this seems to use the "general" complex power function
which loses a few bits -- unavoidably.
Our Windows-version of complex  a ^ b  does so as well.

Here's an R version of what (I believe) once was the C library
cpow(); at least I confirm that it has the same slight
inaccuracy [if you are as this very border line case with '1e-15';
	    use  1e-12  and you have no problems !! ]

Cpow <- function(a,b)
{
  ## Purpose: a ^ b  in complex
  ## Find bug in  complex_pow()
  ## -------------------------------------------------------------------------
  ## Author: Martin Maechler, Date: 15 Jan 2000, 21:33

    a <- as.complex(a)
    b <- as.complex(b)

    hypot <- function(x,y)Mod(x + 1i*y)

    logr <- log(hypot(Re(a), Im(a)) );
    logi <- atan2(Im(a), Re(a));
    x <- exp( logr * Re(b) - logi * Im(b) );
    y <- logr * Im(b) + logi * Re(b);

    x * complex(re= cos(y), im= sin(y))
}

----------------

So, yes we could check for the special case of  "^2" and use
multiplication and then for   " ^ n " and ... ...
and only otherwise call cpow(x,y) {or our own C-version of
that}.

I'm looking into implenting that now.
Expect to hear more about it within 24 hours.

Martin Maechler, ETH Zurich


    HWB> Hans W. Borchers


    HWB> Martin Becker <martin.becker <at> mx.uni-saarland.de>
    HWB> 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,
    HWB> 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
    >> 

______________________________________________
    HWB> R-devel at r-project.org mailing list
    HWB> https://stat.ethz.ch/mailman/listinfo/r-devel