Re: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate)
Dear All,
I just tested the complex-step derivative method on the latest windows
binary version: R 2.10.0 pre-release.
The complex-step derivatives are "exact" to within machine epsilon, with
very little additional computational effort compared to the simple,
first-order forward difference scheme.
Thanks very much to all, especially to, Martin Maechler.
Best,
Ravi.
-Original Message-
From: Martin Maechler [mailto:maech...@stat.math.ethz.ch]
Sent: Wednesday, August 05, 2009 3:35 AM
To: John Nolan
Cc: Ravi Varadhan; hwborch...@googlemail.com; r-de...@stat.math.ethz.ch
Subject: Re: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate)
>>>>> "JN" == John Nolan
>>>>> on Tue, 4 Aug 2009 18:05:47 -0400 writes:
JN> Ravi,
JN> There has been a lot of chatter about this, and people don't seem to
be
JN> reading carefully. Perhaps this will help clarify things.
Yes, I hope so; thanks a lot, John!!
{Indeed, I've been appalled too about the amount of proclaiming
without listening, and even more about the no-bug report... }
JN> The problem appears to be that R was evaluating x^2 not as
multiplication
JN> x*x, but as x^2.0=exp(2.0*log(x)), using standard C functions for
the
JN> complex log and exponentiation. Apparently the C library functions
for
JN> this have some inaccuracies in it, which showed up in your
calculations.
JN> One of the other responders said that matlab, octave, etc. detect
the case
JN> when there is a positive integer power and use multiplication to
evaluate.
JN> It appears that Martin Maechler has now submitted a change to R to
detect
JN> integer powers and evaluate them via repeated multiplication, which
should
JN> eliminate your problem.
Yes; note however that -- exactly for the reason John has given
(and Martin Becker and I had tried to explain earlier in this
thread!) --
that your problem (*yours* indeed, not R's !!) will resurface if
you change "^2" by "^2.01" in your original problem.
JN> There is no guarantee that R or matlab or any other
JN> program will evaluate every expression correctly.
JN> matlab has many years of use and evaluation that have
JN> guided them in correct problems. R is catching up, but
JN> not there yet.
hhmm, hmm,... are you sure that the catching up is only on our side?
I dare say that I'm pretty sure R does quite a few computations
better than e.g., matlab... But let's not get into more unproductive
chattering...
JN> We are all part of improving R; your question led to a careful
JN> examination of the issue and a fix within a few days. No commercial
JN> company responds so quickly!
JN> HTH, John
it did, thank you very much!
Martin
JN>
...
JN> John P. Nolan
JN> Math/Stat Department
JN> 227 Gray Hall
JN> American University
JN> 4400 Massachusetts Avenue, NW
JN> Washington, DC 20016-8050
JN> jpno...@american.edu
JN> 202.885.3140 voice
JN> 202.885.3155 fax
JN> http://academic2.american.edu/~jpnolan
JN>
...
JN> -r-devel-boun...@r-project.org wrote: -
JN> To: "'Martin Becker'"
JN> From: "Ravi Varadhan"
JN> Sent by: r-devel-boun...@r-project.org
JN> Date: 08/04/2009 10:59AM
JN> cc: hwborch...@googlemail.com, r-de...@stat.math.ethz.ch
JN> Subject: Re: [Rd] Inaccurate complex arithmetic of R (Matlab is
accurate)
JN> Please forgive me for my lack of understanding of IEEE
floating-point
JN> arithmetic. I have a hard time undertsanding why "this is not a
problem of
JN> R itself", when "ALL" the other well known computing environments
including
JN> Matlab, Octave, S+, and Scilab provide accurate results. My concern
is not
JN> really about the "overall" accuracy of R, but just about the complex
JN> arithmetic. Is there something idiosyncratic about the complex
arithmetic?
JN> I am really hoping that some one from the R core would speak up and
address
JN> this issue. It would be a loss to the R users if this fascinating
idea of
JN> "complex-step derivative" could not be implemented in R.
JN> Thanks,
JN> Ravi.
JN>
JN> ---
JN> Ravi Varadhan, Ph.D.
JN> Assistant Professor, The Center on Aging and Health
JN> Division of Geriatric Medicine and Gerontology
JN> Johns Hopkins University
JN> Ph: (410) 502-2619
JN>
Re: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate)
Dear John, Thank you. I am amazed at the insightful responses of people in the R-devel group and even more so amazed at the rapdity with which the R core group and other knowledgeable R folks respond to and fix the problems. Big thanks to Martin Becker, Ash Richardson, Duncan Murdoch and Martin Maechler! My frustration was only due to my impatience in having to wait a few days (!) to experiment with this marvellous idea of analytic continuation for computing numerical derivatives that are "exact", but with only minimal computation. 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 - Original Message - From: John Nolan Date: Tuesday, August 4, 2009 8:36 pm Subject: Re: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate) To: Ravi Varadhan Cc: 'Martin Becker' , hwborch...@googlemail.com, r-de...@stat.math.ethz.ch > Ravi, > > There has been a lot of chatter about this, and people don't seem to > be > reading carefully. Perhaps this will help clarify things. > > The problem appears to be that R was evaluating x^2 not as multiplication > x*x, but as x^2.0=exp(2.0*log(x)), using standard C functions for the > complex log and exponentiation. Apparently the C library functions for > this have some inaccuracies in it, which showed up in your calculations. > > One of the other responders said that matlab, octave, etc. detect the > case > when there is a positive integer power and use multiplication to evaluate. > It appears that Martin Maechler has now submitted a change to R to detect > integer powers and evaluate them via repeated multiplication, which should > eliminate your problem. > > There is no guarantee that R or matlab or any other program will evaluate > every expression correctly. matlab has many years of use and evaluation > that have guided them in correct problems. R is catching up, but not > there > yet. We are all part of improving R; your question led to a careful > examination of the issue and a fix within a few days. No commercial > company responds so quickly! > > HTH, John > > ... > > John P. Nolan > Math/Stat Department > 227 Gray Hall > American University > 4400 Massachusetts Avenue, NW > Washington, DC 20016-8050 > > jpno...@american.edu > 202.885.3140 voice > 202.885.3155 fax > > ... > > > > -r-devel-boun...@r-project.org wrote: - > > > To: "'Martin Becker'" > From: "Ravi Varadhan" > Sent by: r-devel-boun...@r-project.org > Date: 08/04/2009 10:59AM > cc: hwborch...@googlemail.com, r-de...@stat.math.ethz.ch > Subject: Re: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate) > > Please forgive me for my lack of understanding of IEEE floating-point > arithmetic. I have a hard time undertsanding why "this is not a > problem of > R itself", when "ALL" the other well known computing environments including > Matlab, Octave, S+, and Scilab provide accurate results. My concern > is not > really about the "overall" accuracy of R, but just about the complex > arithmetic. Is there something idiosyncratic about the complex arithmetic? > > > I am really hoping that some one from the R core would speak up and address > this issue. It would be a loss to the R users if this fascinating > idea of > "complex-step derivative" could not be implemented in R. > > Thanks, > 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: > > > tml > > > > > > > > > -Original Message- > From: Martin Becker [ > Sent: Tuesday, August 04, 2009 7:34 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, > > I suspect that, in gener
Re: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate)
>>>>> "JN" == John Nolan
>>>>> on Tue, 4 Aug 2009 18:05:47 -0400 writes:
JN> Ravi,
JN> There has been a lot of chatter about this, and people don't seem to be
JN> reading carefully. Perhaps this will help clarify things.
Yes, I hope so; thanks a lot, John!!
{Indeed, I've been appalled too about the amount of proclaiming
without listening, and even more about the no-bug report... }
JN> The problem appears to be that R was evaluating x^2 not as
multiplication
JN> x*x, but as x^2.0=exp(2.0*log(x)), using standard C functions for the
JN> complex log and exponentiation. Apparently the C library functions for
JN> this have some inaccuracies in it, which showed up in your calculations.
JN> One of the other responders said that matlab, octave, etc. detect the
case
JN> when there is a positive integer power and use multiplication to
evaluate.
JN> It appears that Martin Maechler has now submitted a change to R to
detect
JN> integer powers and evaluate them via repeated multiplication, which
should
JN> eliminate your problem.
Yes; note however that -- exactly for the reason John has given
(and Martin Becker and I had tried to explain earlier in this
thread!) --
that your problem (*yours* indeed, not R's !!) will resurface if
you change "^2" by "^2.01" in your original problem.
JN> There is no guarantee that R or matlab or any other
JN> program will evaluate every expression correctly.
JN> matlab has many years of use and evaluation that have
JN> guided them in correct problems. R is catching up, but
JN> not there yet.
hhmm, hmm,... are you sure that the catching up is only on our side?
I dare say that I'm pretty sure R does quite a few computations
better than e.g., matlab... But let's not get into more unproductive
chattering...
JN> We are all part of improving R; your question led to a careful
JN> examination of the issue and a fix within a few days. No commercial
JN> company responds so quickly!
JN> HTH, John
it did, thank you very much!
Martin
JN>
...
JN> John P. Nolan
JN> Math/Stat Department
JN> 227 Gray Hall
JN> American University
JN> 4400 Massachusetts Avenue, NW
JN> Washington, DC 20016-8050
JN> jpno...@american.edu
JN> 202.885.3140 voice
JN> 202.885.3155 fax
JN> http://academic2.american.edu/~jpnolan
JN>
...
JN> -r-devel-boun...@r-project.org wrote: -
JN> To: "'Martin Becker'"
JN> From: "Ravi Varadhan"
JN> Sent by: r-devel-boun...@r-project.org
JN> Date: 08/04/2009 10:59AM
JN> cc: hwborch...@googlemail.com, r-de...@stat.math.ethz.ch
JN> Subject: Re: [Rd] Inaccurate complex arithmetic of R (Matlab is
accurate)
JN> Please forgive me for my lack of understanding of IEEE floating-point
JN> arithmetic. I have a hard time undertsanding why "this is not a
problem of
JN> R itself", when "ALL" the other well known computing environments
including
JN> Matlab, Octave, S+, and Scilab provide accurate results. My concern is
not
JN> really about the "overall" accuracy of R, but just about the complex
JN> arithmetic. Is there something idiosyncratic about the complex
arithmetic?
JN> I am really hoping that some one from the R core would speak up and
address
JN> this issue. It would be a loss to the R users if this fascinating idea
of
JN> "complex-step derivative" could not be implemented in R.
JN> Thanks,
JN> Ravi.
JN>
JN> ---
JN> Ravi Varadhan, Ph.D.
JN> Assistant Professor, The Center on Aging and Health
JN> Division of Geriatric Medicine and Gerontology
JN> Johns Hopkins University
JN> Ph: (410) 502-2619
JN> Fax: (410) 614-9625
JN> Email: rvarad...@jhmi.edu
JN> Webpage:
JN>
http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h
JN> tml
JN>
------------------------------------
JN>
JN> -Original Message-
JN> From: Martin Becker [mailto:martin.bec...@mx.uni-saarland.de]
JN> Sent: Tuesday, August 04, 2009 7:34 AM
JN> To: Ravi Varadhan
JN> Cc: r-de...@stat.math.ethz.ch; hwborch...@googlemail.com
JN> Subject: Re: [Rd] Inaccurate complex arithmetic of R (Matlab is
accurate)
JN> Dear Ravi,
Re: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate)
Ravi, There has been a lot of chatter about this, and people don't seem to be reading carefully. Perhaps this will help clarify things. The problem appears to be that R was evaluating x^2 not as multiplication x*x, but as x^2.0=exp(2.0*log(x)), using standard C functions for the complex log and exponentiation. Apparently the C library functions for this have some inaccuracies in it, which showed up in your calculations. One of the other responders said that matlab, octave, etc. detect the case when there is a positive integer power and use multiplication to evaluate. It appears that Martin Maechler has now submitted a change to R to detect integer powers and evaluate them via repeated multiplication, which should eliminate your problem. There is no guarantee that R or matlab or any other program will evaluate every expression correctly. matlab has many years of use and evaluation that have guided them in correct problems. R is catching up, but not there yet. We are all part of improving R; your question led to a careful examination of the issue and a fix within a few days. No commercial company responds so quickly! HTH, John ... John P. Nolan Math/Stat Department 227 Gray Hall American University 4400 Massachusetts Avenue, NW Washington, DC 20016-8050 jpno...@american.edu 202.885.3140 voice 202.885.3155 fax http://academic2.american.edu/~jpnolan ... -r-devel-boun...@r-project.org wrote: - To: "'Martin Becker'" From: "Ravi Varadhan" Sent by: r-devel-boun...@r-project.org Date: 08/04/2009 10:59AM cc: hwborch...@googlemail.com, r-de...@stat.math.ethz.ch Subject: Re: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate) Please forgive me for my lack of understanding of IEEE floating-point arithmetic. I have a hard time undertsanding why "this is not a problem of R itself", when "ALL" the other well known computing environments including Matlab, Octave, S+, and Scilab provide accurate results. My concern is not really about the "overall" accuracy of R, but just about the complex arithmetic. Is there something idiosyncratic about the complex arithmetic? I am really hoping that some one from the R core would speak up and address this issue. It would be a loss to the R users if this fascinating idea of "complex-step derivative" could not be implemented in R. Thanks, 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: Tuesday, August 04, 2009 7:34 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, I suspect that, in general, you may be facing the limitations of machine accuracy (more precisely, IEEE 754 arithmetics on [64-bit] doubles) in your application. This is not a problem of R itself, but rather a problem of standard arithmetics provided by underlying C compilers/CPUs. In fact, every operation in IEEE arithmetics (so, this is not really a problem only for complex numbers) may suffer from inexactness, a particularly difficult one is addition/subtraction. Consider the following example for real numbers (I know, it is not a very good one...): The two functions badfn <- function(x) 1-(1+x)*(1-x) goodfn <- function(x) x^2 both calculate x^2 (theoretically, given perfect arithmetic). So, as you want to allow the user to 'specify the mathematical function ... in "any" form the user chooses', both functions should be ok. But, unfortunately: > badfn(1e-8) [1] 2.220446049250313e-16 > goodfn(1e-8) [1] 1e-16 I don't know what happens in matlab/octave/scilab for this example. They may do better, but probably at some cost (dropping IEEE arithmetic/do "clever" calculations should result in massive speed penalties, try evaluating hypergeom([1,-99.9],[-49.9-24.9*I],(1+1.71*I)/2); in Maple...). Now, you have some options: - assume, that the user is aware of the numerical inexactness of ieee arithmetics and that he is able to supply some "robust" version of the mathematical function. - use some other software (eg., matlab) for the critical calculations (there is a R <-> Matlab interface, see package R.matlab on CRAN), if you are sure, that this helps.
Re: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate)
On 04/08/2009 1:41 PM, Richardson, Ash wrote:
Dear Ravi,
Here's an idea: if R uses the C++ library complex.h, then the problem may lie
in the C++ library.
(Or how the C++ library is being called: perhaps it's defaulting to single
precision floating point).
I recently wrote an eigen-solver for a complex matrix using explicit formula;
using the C++ library I got residuals of 10^-5.
Using a complex number class I wrote myself using the polar representation,
the same program got residuals of 10^-15.
The complex functions in standard C don't seem to have the same problem.
I think there is something wrong with the C++ implementation of complex numbers.
I think matlab uses the standard C implementation of complex numbers, so it
doesn't suffer from this problem.
R uses the C runtime, not C++. But the cpow function was not very good.
Blame Microsoft for that on Windows, some other scapegoat on whatever
other platform you use.
Martin M. has now committed some changes to R-devel that give these
results on Windows:
> #-- 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)
[1] 190.3080
> Im (rx)
[1] -4.667764e-15
These seem to agree with the desired numbers.
Duncan Murdoch
Best regards,
~ Ash
__
Ash Richardson, BSc
Physical Scientist / Spécialiste des science physiques
Natural Resources Canada / Ressources naturelles Canada
Canadian Forest Service / Service canadien des forêts
Pacific Forestry Centre / Centre de foresterie du Pacifique
506 W. Burnside Road / 506 rue Burnside Ouest
Victoria, BC V8Z 1M5 / Victoria, C-B V8Z 1M5
Tel: (250) 363-6018 Facs: (250) 363-0775
mailto:ashri...@nrcan.gc.ca
___
-Original Message-
From: r-devel-boun...@r-project.org on behalf of Ravi Varadhan
Sent: Tue 8/4/2009 7:59 AM
To: 'Martin Becker'
Cc: hwborch...@googlemail.com; r-de...@stat.math.ethz.ch
Subject: Re: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate)
Please forgive me for my lack of understanding of IEEE floating-point
arithmetic. I have a hard time undertsanding why "this is not a problem of
R itself", when "ALL" the other well known computing environments including
Matlab, Octave, S+, and Scilab provide accurate results. My concern is not
really about the "overall" accuracy of R, but just about the complex
arithmetic. Is there something idiosyncratic about the complex arithmetic?
I am really hoping that some one from the R core would speak up and address
this issue. It would be a loss to the R users if this fascinating idea of
"complex-step derivative" could not be implemented in R.
Thanks,
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: Tuesday, August 04, 2009 7:34 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,
I suspect that, in general, you may be facing the limitations of machine
accuracy (more precisely, IEEE 754 arithmetics on [64-bit] doubles) in your
application. This is not a problem of R itself, but rather a problem of
standard arithmetics provided by underlying C compilers/CPUs.
In fact, every operation in IEEE arithmetics (so, this is not really a
problem only for complex numbers) may suffer from inexactness, a
particularly difficult one is addition/subtraction. Consider the following
example for real numbers (I know, it is not a very good one...):
The two functions
badfn <- function(x) 1-(1+x)*(1-x)
goodfn <- function(x) x^2
both calculate x^2 (theoretically, given perfect arithmetic). So, as you
want to allow the user to 'specify the mathematical function ... in "any"
form the user chooses', both functions should be ok.
But, unfortunately:
> badfn(1e-8)
[1] 2.220446049250313e-16
> goodfn(1e-8)
[1] 1e-16
I don't know what happens in matlab/octave/scilab for this example. They may
do better, but probably at some cost (dropping IEEE arithmet
Re: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate)
Dear Ravi, Here's an idea: if R uses the C++ library complex.h, then the problem may lie in the C++ library. (Or how the C++ library is being called: perhaps it's defaulting to single precision floating point). I recently wrote an eigen-solver for a complex matrix using explicit formula; using the C++ library I got residuals of 10^-5. Using a complex number class I wrote myself using the polar representation, the same program got residuals of 10^-15. The complex functions in standard C don't seem to have the same problem. I think there is something wrong with the C++ implementation of complex numbers. I think matlab uses the standard C implementation of complex numbers, so it doesn't suffer from this problem. Best regards, ~ Ash __ Ash Richardson, BSc Physical Scientist / Spécialiste des science physiques Natural Resources Canada / Ressources naturelles Canada Canadian Forest Service / Service canadien des forêts Pacific Forestry Centre / Centre de foresterie du Pacifique 506 W. Burnside Road / 506 rue Burnside Ouest Victoria, BC V8Z 1M5 / Victoria, C-B V8Z 1M5 Tel: (250) 363-6018 Facs: (250) 363-0775 mailto:ashri...@nrcan.gc.ca ___ -Original Message- From: r-devel-boun...@r-project.org on behalf of Ravi Varadhan Sent: Tue 8/4/2009 7:59 AM To: 'Martin Becker' Cc: hwborch...@googlemail.com; r-de...@stat.math.ethz.ch Subject: Re: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate) Please forgive me for my lack of understanding of IEEE floating-point arithmetic. I have a hard time undertsanding why "this is not a problem of R itself", when "ALL" the other well known computing environments including Matlab, Octave, S+, and Scilab provide accurate results. My concern is not really about the "overall" accuracy of R, but just about the complex arithmetic. Is there something idiosyncratic about the complex arithmetic? I am really hoping that some one from the R core would speak up and address this issue. It would be a loss to the R users if this fascinating idea of "complex-step derivative" could not be implemented in R. Thanks, 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: Tuesday, August 04, 2009 7:34 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, I suspect that, in general, you may be facing the limitations of machine accuracy (more precisely, IEEE 754 arithmetics on [64-bit] doubles) in your application. This is not a problem of R itself, but rather a problem of standard arithmetics provided by underlying C compilers/CPUs. In fact, every operation in IEEE arithmetics (so, this is not really a problem only for complex numbers) may suffer from inexactness, a particularly difficult one is addition/subtraction. Consider the following example for real numbers (I know, it is not a very good one...): The two functions badfn <- function(x) 1-(1+x)*(1-x) goodfn <- function(x) x^2 both calculate x^2 (theoretically, given perfect arithmetic). So, as you want to allow the user to 'specify the mathematical function ... in "any" form the user chooses', both functions should be ok. But, unfortunately: > badfn(1e-8) [1] 2.220446049250313e-16 > goodfn(1e-8) [1] 1e-16 I don't know what happens in matlab/octave/scilab for this example. They may do better, but probably at some cost (dropping IEEE arithmetic/do "clever" calculations should result in massive speed penalties, try evaluating hypergeom([1,-99.9],[-49.9-24.9*I],(1+1.71*I)/2); in Maple...). Now, you have some options: - assume, that the user is aware of the numerical inexactness of ieee arithmetics and that he is able to supply some "robust" version of the mathematical function. - use some other software (eg., matlab) for the critical calculations (there is a R <-> Matlab interface, see package R.matlab on CRAN), if you are sure, that this helps. - implement/use multiple precision arithmetics within R (Martin Maechler's Rmpfr package may be very useful: http://r-forge.r-project.org/projects/rmpfr/ , but this will slow down calculations considerably) All in all, I think it is unfair just to blame R here. Of course, it woul
Re: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate)
Please forgive me for my lack of understanding of IEEE floating-point arithmetic. I have a hard time undertsanding why "this is not a problem of R itself", when "ALL" the other well known computing environments including Matlab, Octave, S+, and Scilab provide accurate results. My concern is not really about the "overall" accuracy of R, but just about the complex arithmetic. Is there something idiosyncratic about the complex arithmetic? I am really hoping that some one from the R core would speak up and address this issue. It would be a loss to the R users if this fascinating idea of "complex-step derivative" could not be implemented in R. Thanks, 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: Tuesday, August 04, 2009 7:34 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, I suspect that, in general, you may be facing the limitations of machine accuracy (more precisely, IEEE 754 arithmetics on [64-bit] doubles) in your application. This is not a problem of R itself, but rather a problem of standard arithmetics provided by underlying C compilers/CPUs. In fact, every operation in IEEE arithmetics (so, this is not really a problem only for complex numbers) may suffer from inexactness, a particularly difficult one is addition/subtraction. Consider the following example for real numbers (I know, it is not a very good one...): The two functions badfn <- function(x) 1-(1+x)*(1-x) goodfn <- function(x) x^2 both calculate x^2 (theoretically, given perfect arithmetic). So, as you want to allow the user to 'specify the mathematical function ... in "any" form the user chooses', both functions should be ok. But, unfortunately: > badfn(1e-8) [1] 2.220446049250313e-16 > goodfn(1e-8) [1] 1e-16 I don't know what happens in matlab/octave/scilab for this example. They may do better, but probably at some cost (dropping IEEE arithmetic/do "clever" calculations should result in massive speed penalties, try evaluating hypergeom([1,-99.9],[-49.9-24.9*I],(1+1.71*I)/2); in Maple...). Now, you have some options: - assume, that the user is aware of the numerical inexactness of ieee arithmetics and that he is able to supply some "robust" version of the mathematical function. - use some other software (eg., matlab) for the critical calculations (there is a R <-> Matlab interface, see package R.matlab on CRAN), if you are sure, that this helps. - implement/use multiple precision arithmetics within R (Martin Maechler's Rmpfr package may be very useful: http://r-forge.r-project.org/projects/rmpfr/ , but this will slow down calculations considerably) All in all, I think it is unfair just to blame R here. Of course, it would be great if there was a simple trigger to turn on multiple precision arithmetics in R. Packages such as Rmpfr may provide a good step in this direction, since operator overloading via S4 classes allows for easy code adaption. But Rmpfr is still declared "beta", and it relies on some external library, which could be problematic on Windows systems. Maybe someone else has other/better suggestions, but I do not think that there is an easy solution for the "general" problem. Best wishes, Martin Ravi Varadhan wrote: > 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. > > --
Re: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate)
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.6677637664e-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
Re: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate)
> "MM" == Martin Maechler
> on Mon, 3 Aug 2009 19:30:24 +0200 writes:
> "HWB" == Hans W Borchers
> on Mon, 3 Aug 2009 13:15:11 + (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 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@r-project.org mailing list
HWB> https://stat.ethz.ch/mailman/listinfo/r-devel
MM> __
MM> R-devel@r-project.org mailing list
MM> https://stat.ethz.ch/mailman/listinfo/r-devel
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate)
> I suspect that, in general, you may be facing the limitations of machine
> accuracy (more precisely, IEEE 754 arithmetics on [64-bit] doubles) in
Dear Martin,
I definitely do not agree with this. Consider your own proposal of
writing the Rosenbrock function:
rosen2 <- 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))
}
The complex-step derivative approximation suggests to set h <- 1.0e-20
and then to compute the expression Im(rosen(x+hi))/h , so let's do it:
h <- 1.0e-20
xh <- c(0.0094 + h*1i, 0.7146, 0.2179, 0.6883, 0.5757,
0.9549,0.7136, 0.0849, 0.4147, 0.4540)
Im(rosen2(xh)) / h
# [1] -4.6677637664000
which is exactly the correct derivative in the first variable, namely
d/dx1 rosen(x). To verify, you could even calculate this by hand on a
"Bierdeckel".
Now look at the former definition, say rosen1(), using '^' instead:
rosen1 <- function(x) {
n <- length(x)
x1 <- x[2:n]
x2 <- x[1:(n-1)]
sum(100*(x1-x2^2)^2 + (1-x2)^2)
}
Im(rosen1(xh)) / h
# [1] -747143.8793904837
We find that R is "running wild". And this has nothing to do with IEEE
arithmetics or machine accuracy, as we can see from the first example
where R is able to do it right.
And as I said, the second example is working correctly in free software
such as Octave which I guess does not do any "clever" calculations here.
We do _not_ ask for multiple precision arithmetics here !
Regards
Hans Werner
Martin Becker mx.uni-saarland.de> writes:
>
> Dear Ravi,
>
> I suspect that, in general, you may be facing the limitations of machine
> accuracy (more precisely, IEEE 754 arithmetics on [64-bit] doubles) in
> your application. This is not a problem of R itself, but rather a
> problem of standard arithmetics provided by underlying C compilers/CPUs.
> In fact, every operation in IEEE arithmetics (so, this is not really a
> problem only for complex numbers) may suffer from inexactness, a
> particularly difficult one is addition/subtraction. Consider the
> following example for real numbers (I know, it is not a very good one...):
> The two functions
>
> badfn <- function(x) 1-(1+x)*(1-x)
> goodfn <- function(x) x^2
>
> both calculate x^2 (theoretically, given perfect arithmetic). So, as you
> want to allow the user to 'specify the mathematical function ... in
> "any" form the user chooses', both functions should be ok.
> But, unfortunately:
>
> > badfn(1e-8)
> [1] 2.220446049250313e-16
> > goodfn(1e-8)
> [1] 1e-16
>
> I don't know what happens in matlab/octave/scilab for this example. They
> may do better, but probably at some cost (dropping IEEE arithmetic/do
> "clever" calculations should result in massive speed penalties, try
> evaluating hypergeom([1,-99.9],[-49.9-24.9*I],(1+1.71*I)/2); in
> Maple...).
> Now, you have some options:
>
> - assume, that the user is aware of the numerical inexactness of ieee
> arithmetics and that he is able to supply some "robust" version of the
> mathematical function.
> - use some other software (eg., matlab) for the critical calculations
> (there is a R <-> Matlab interface, see package R.matlab on CRAN), if
> you are sure, that this helps.
> - implement/use multiple precision arithmetics within R (Martin
> Maechler's Rmpfr package may be very useful:
> http://r-forge.r-project.org/projects/rmpfr/ , but this will slow down
> calculations considerably)
>
> All in all, I think it is unfair just to blame R here. Of course, it
> would be great if there was a simple trigger to turn on multiple
> precision arithmetics in R. Packages such as Rmpfr may provide a good
> step in this direction, since operator overloading via S4 classes allows
> for easy code adaption. But Rmpfr is still declared "beta", and it
> relies on some external library, which could be problematic on Windows
> systems. Maybe someone else has other/better suggestions, but I do not
> think that there is an easy solution for the "general" problem.
>
> Best wishes,
>
> Martin
>
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate)
I checked, and both octave and yorick use multiplication for z^i where i is an integer, leading to better accuracy. Octave uses an integer power if it's stored as a double if it's close enough to an integer. See: http://hg.savannah.gnu.org/hgweb/octave/file/fb22dd5d6242/src/xpow.cc http://yorick.sourcearchive.com/documentation/2.1.01cvs20060706/ops2_8c-source.html - Tom __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate)
Dear Ravi,
I suspect that, in general, you may be facing the limitations of machine
accuracy (more precisely, IEEE 754 arithmetics on [64-bit] doubles) in
your application. This is not a problem of R itself, but rather a
problem of standard arithmetics provided by underlying C compilers/CPUs.
In fact, every operation in IEEE arithmetics (so, this is not really a
problem only for complex numbers) may suffer from inexactness, a
particularly difficult one is addition/subtraction. Consider the
following example for real numbers (I know, it is not a very good one...):
The two functions
badfn <- function(x) 1-(1+x)*(1-x)
goodfn <- function(x) x^2
both calculate x^2 (theoretically, given perfect arithmetic). So, as you
want to allow the user to 'specify the mathematical function ... in
"any" form the user chooses', both functions should be ok.
But, unfortunately:
> badfn(1e-8)
[1] 2.220446049250313e-16
> goodfn(1e-8)
[1] 1e-16
I don't know what happens in matlab/octave/scilab for this example. They
may do better, but probably at some cost (dropping IEEE arithmetic/do
"clever" calculations should result in massive speed penalties, try
evaluating hypergeom([1,-99.9],[-49.9-24.9*I],(1+1.71*I)/2); in
Maple...).
Now, you have some options:
- assume, that the user is aware of the numerical inexactness of ieee
arithmetics and that he is able to supply some "robust" version of the
mathematical function.
- use some other software (eg., matlab) for the critical calculations
(there is a R <-> Matlab interface, see package R.matlab on CRAN), if
you are sure, that this helps.
- implement/use multiple precision arithmetics within R (Martin
Maechler's Rmpfr package may be very useful:
http://r-forge.r-project.org/projects/rmpfr/ , but this will slow down
calculations considerably)
All in all, I think it is unfair just to blame R here. Of course, it
would be great if there was a simple trigger to turn on multiple
precision arithmetics in R. Packages such as Rmpfr may provide a good
step in this direction, since operator overloading via S4 classes allows
for easy code adaption. But Rmpfr is still declared "beta", and it
relies on some external library, which could be problematic on Windows
systems. Maybe someone else has other/better suggestions, but I do not
think that there is an easy solution for the "general" problem.
Best wishes,
Martin
Ravi Varadhan wrote:
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 <
Re: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate)
> "HWB" == Hans W Borchers
> on Mon, 3 Aug 2009 13:15:11 + (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 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@r-project.org mailing list
HWB> https://stat.ethz.ch/mailman/listinfo/r-devel
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] Inaccurate complex arithmetic of R (Matlab is accurate)
>
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 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
>
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
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.6677637664e-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
