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

[Martin Maechler]

>>>>> "MM" == Martin Maechler <maechler at stat.math.ethz.ch>
>>>>>     on Mon, 3 Aug 2009 19:30:24 +0200 writes:

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

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

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

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

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

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

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

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

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

    MM> ----------------

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

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

I have now committed a change to R-devel (not sure if to be
back-ported to R-patched)
which uses  ~ log2(n) multiplications  for  z^n  when
n is integer (and  |n| <= 2^16 ; that cut-off being somewhat arbitrary).

Along the same line, I've looked what we do for  x^y  when x,y
are double. Till now, we had only special cased the cases
y == 0 (, 1), 2; and after some simple tests, I've decided to
use the log(n) #{multiplications} algorithm, whenever |n| <= 256.

Thanks also to Tom Short for investigating what other free
packages use.

Martin Maechler, ETH Zurich




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

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

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