[R] Numerical stability of: 1/(1 - cos(x)) - 2/x^2
@vi@e@gross m@iii@g oii gm@ii@com
@vi@e@gross m@iii@g oii gm@ii@com
Sat Aug 19 01:55:12 CEST 2023
This discussion is sooo familiar.
If you want indefinite precision arithmetic, feel free to use a language and data type that supports it.
Otherwise, only do calculations that fit in a safe zone.
This is not just about this scenario. Floating point can work well when adding (or subtracting) two numbers of about the same size. But if one number is .123456789... and another is the same except raised to the -45th power of ten, then adding them effectively throws away the second number.
This is a well-known problem for any finite binary representation of numbers. In the example given, yes, the smaller the number is, the worse the behavior in this case tends to be.
There are many solutions and some are fairly expensive in terms of computation time and sometimes memory usage.
Are there any good indefinite (or much higher) precision packages out there that would not only support the data type needed but also properly be used and passed along to the functions used to do complex calculations? No, I am not asking for indefinite precision complex numbers, but generally that would be a tuple of such numbers.
-----Original Message-----
From: R-help <r-help-bounces using r-project.org> On Behalf Of Bert Gunter
Sent: Friday, August 18, 2023 7:06 PM
To: Leonard Mada <leo.mada using syonic.eu>
Cc: R-help Mailing List <r-help using r-project.org>; Martin Maechler <maechler using stat.math.ethz.ch>
Subject: Re: [R] Numerical stability of: 1/(1 - cos(x)) - 2/x^2
"The ugly thing is that the error only gets worse as x decreases. The
value neither drops to 0, nor does it blow up to infinity; but it gets
worse in a continuous manner."
If I understand you correctly, this is wrong:
> x <- 2^(-20) ## considerably less then 1e-4 !!
> y <- 1 - x^2/2;
> 1/(1 - y) - 2/x^2
[1] 0
It's all about the accuracy of the binary approximation of floating point
numbers (and their arithmetic)
Cheers,
Bert
On Fri, Aug 18, 2023 at 3:25 PM Leonard Mada via R-help <
r-help using r-project.org> wrote:
> I have added some clarifications below.
>
> On 8/18/2023 10:20 PM, Leonard Mada wrote:
> > [...]
> > After more careful thinking, I believe that it is a limitation due to
> > floating points:
> > [...]
> >
> > The problem really stems from the representation of 1 - x^2/2 as shown
> > below:
> > x = 1E-4
> > print(1 - x^2/2, digits=20)
> > print(0.999999995, digits=20) # fails
> > # 0.99999999500000003039
>
> The floating point representation of 1 - x^2/2 is the real culprit:
> # 0.99999999500000003039
>
> The 3039 at the end is really an error due to the floating point
> representation. However, this error blows up when inverting the value:
> x = 1E-4;
> y = 1 - x^2/2;
> 1/(1 - y) - 2/x^2
> # 1.215494
> # should be 1/(x^2/2) - 2/x^2 = 0
>
>
> The ugly thing is that the error only gets worse as x decreases. The
> value neither drops to 0, nor does it blow up to infinity; but it gets
> worse in a continuous manner. At least the reason has become now clear.
>
>
> >
> > Maybe some functions of type cos1p and cos1n would be handy for such
> > computations (to replace the manual series expansion):
> > cos1p(x) = 1 + cos(x)
> > cos1n(x) = 1 - cos(x)
> > Though, I do not have yet the big picture.
> >
>
> Sincerely,
>
>
> Leonard
>
> >
> >
> > On 8/17/2023 1:57 PM, Martin Maechler wrote:
> >>>>>>> Leonard Mada
> >>>>>>> on Wed, 16 Aug 2023 20:50:52 +0300 writes:
> >> > Dear Iris,
> >> > Dear Martin,
> >>
> >> > Thank you very much for your replies. I add a few comments.
> >>
> >> > 1.) Correct formula
> >> > The formula in the Subject Title was correct. A small glitch
> >> swept into
> >> > the last formula:
> >> > - 1/(cos(x) - 1) - 2/x^2
> >> > or
> >> > 1/(1 - cos(x)) - 2/x^2 # as in the subject title;
> >>
> >> > 2.) log1p
> >> > Actually, the log-part behaves much better. And when it fails,
> >> it fails
> >> > completely (which is easy to spot!).
> >>
> >> > x = 1E-6
> >> > log(x) -log(1 - cos(x))/2
> >> > # 0.3465291
> >>
> >> > x = 1E-8
> >> > log(x) -log(1 - cos(x))/2
> >> > # Inf
> >> > log(x) - log1p(- cos(x))/2
> >> > # Inf => fails as well!
> >> > # although using only log1p(cos(x)) seems to do the trick;
> >> > log1p(cos(x)); log(2)/2;
> >>
> >> > 3.) 1/(1 - cos(x)) - 2/x^2
> >> > It is possible to convert the formula to one which is
> >> numerically more
> >> > stable. It is also possible to compute it manually, but it
> >> involves much
> >> > more work and is also error prone:
> >>
> >> > (x^2 - 2 + 2*cos(x)) / (x^2 * (1 - cos(x)))
> >> > And applying L'Hospital:
> >> > (2*x - 2*sin(x)) / (2*x * (1 - cos(x)) + x^2*sin(x))
> >> > # and a 2nd & 3rd & 4th time
> >> > 1/6
> >>
> >> > The big problem was that I did not expect it to fail for x =
> >> 1E-4. I
> >> > thought it is more robust and works maybe until 1E-5.
> >> > x = 1E-5
> >> > 2/x^2 - 2E+10
> >> > # -3.814697e-06
> >>
> >> > This is the reason why I believe that there is room for
> >> improvement.
> >>
> >> > Sincerely,
> >> > Leonard
> >>
> >> Thank you, Leonard.
> >> Yes, I agree that it is amazing how much your formula suffers from
> >> (a generalization of) "cancellation" --- leading you to think
> >> there was a problem with cos() or log() or .. in R.
> >> But really R uses the system builtin libmath library, and the
> >> problem is really the inherent instability of your formula.
> >>
> >> Indeed your first approximation was not really much more stable:
> >>
> >> ## 3.) 1/(1 - cos(x)) - 2/x^2
> >> ## It is possible to convert the formula to one which is numerically
> >> more
> >> ## stable. It is also possible to compute it manually, but it
> >> involves much
> >> ## more work and is also error prone:
> >> ## (x^2 - 2 + 2*cos(x)) / (x^2 * (1 - cos(x)))
> >> ## MM: but actually, that approximation does not seem better (close
> >> to the breakdown region):
> >> f1 <- \(x) 1/(1 - cos(x)) - 2/x^2
> >> f2 <- \(x) (x^2 - 2 + 2*cos(x)) / (x^2 * (1 - cos(x)))
> >> curve(f1, 1e-8, 1e-1, log="xy" n=2^10)
> >> curve(f2, add = TRUE, col=2, n=2^10)
> >> ## Zoom in:
> >> curve(f1, 1e-4, 1e-1, log="xy",n=2^9)
> >> curve(f2, add = TRUE, col=2, n=2^9)
> >> ## Zoom in much more in y-direction:
> >> yl <- 1/6 + c(-5, 20)/100000
> >> curve(f1, 1e-4, 1e-1, log="x", ylim=yl, n=2^9)
> >> abline(h = 1/6, lty=3, col="gray")
> >> curve(f2, add = TRUE, n=2^9, col=adjustcolor(2, 1/2))
> >>
> >> Now, you can use the Rmpfr package (interface to the GNU MPFR
> >> multiple-precision C library) to find out more :
> >>
> >> if(!requireNamespace("Rmpfr")) install.packages("Rmpfr")
> >> M <- function(x, precBits=128) Rmpfr::mpfr(x, precBits)
> >>
> >> (xM <- M(1e-8))# yes, only ~ 16 dig accurate
> >> ## 1.000000000000000020922560830128472675327e-8
> >> M(10, 128)^-8 # would of course be more accurate,
> >> ## but we want the calculation for the double precision number 1e-8
> >>
> >> ## Now you can draw "the truth" into the above plots:
> >> curve(f1, 1e-4, 1e-1, log="xy",n=2^9)
> >> curve(f2, add = TRUE, col=2, n=2^9)
> >> ## correct:
> >> curve(f1(M(x, 256)), add = TRUE, col=4, lwd=2, n=2^9)
> >> abline(h = 1/6, lty=3, col="gray")
> >>
> >> But, indeed we take note how much it is the formula instability:
> >> Also MPFR needs a lot of extra bits precision before it gets to
> >> the correct numbers:
> >>
> >> xM <- c(M(1e-8, 80), M(1e-8, 96), M(1e-8, 112),
> >> M(1e-8, 128), M(1e-8, 180), M(1e-8, 256))
> >> ## to and round back to 70 bits for display:
> >> R <- \(x) Rmpfr::roundMpfr(x, 70)
> >> R(f1(xM))
> >> R(f2(xM))
> >> ## [1] 0 0
> >> 0.15407439555097885670915
> >> ## [4] 0.16666746653133802175779 0.16666666666666666749979
> >> 0.16666666666666666750001
> >>
> >> ## 1. f1() is even worse than f2() {here at x=1e-8}
> >> ## 2. Indeed, even 96 bits precision is *not* sufficient at all, ...
> >> ## which is amazing to me as well !!
> >>
> >> Best regards,
> >> Martin
>
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
[[alternative HTML version deleted]]
______________________________________________
R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
More information about the R-help
mailing list