10 votes

Why is computation of this function numerically unstable?

The question has been changed since this answer. Another alternative is $$\begin{align}f(x)&=\frac{\sqrt{x+h}-\sqrt{x-h}}{2h}\\ &=\frac{\sqrt{x+h}-\sqrt{x-h}}{2h}\cdot\frac{\sqrt{x+h}+\sqrt{x-...
GoodDeeds's user avatar
  • 851
9 votes

Why is computation of this function numerically unstable?

What you are experiencing is called 'loss of significant digits'. As h gets smaller, the two square-root terms become closer and closer together, in the sense that their values start with more and ...
PMar's user avatar
  • 131
6 votes
Accepted

numerically stable log1pexp calculation

Let $0 < \varepsilon \lll 1$ be the relative error bound of the floating-point system—$2^{-53}$ in IEEE 754 binary64 arithmetic. First, the naive formula ...
Пафнутий Чебышев's user avatar
6 votes
Accepted

The stability of log(1+x)

Consider the case that x is small. (1 + x) has a rounding error; the result that you get is not (1 + x) but (1 + x') for some x' close to x. If x is very small, the relative difference between x' and ...
gnasher729's user avatar
6 votes

fast and stable x * tanh(log1pexp(x)) computation

With some algebraic manipulation (as pointed out in @orlp's answer), we can deduce the following: $$f(x) = x \tanh(\log(1+e^x)) \tag{1}$$ $$ = x\frac{(1+e^x)^2 - 1}{(1+e^x)^2 + 1} = x\frac{e^{2x} + 2e^...
Yashas's user avatar
  • 275
4 votes
Accepted

Do formulas involving fewer repetitions of variables give higher numerical precision?

First, I want to say that it is not the case in general that an algorithm that minimizes the number of uses of the inputs is more accurate, at least for IEEE 754 floating point. For example, ...
Derek Elkins left SE's user avatar
4 votes
Accepted

Avoiding overflows while computing $e^x$ by Taylor series

Of course there are better numerical ways to compute exponential, but if you want to use Taylor expansion only, the better approach is to reformulate the expansion to avoid computing large nominators ...
user172818's user avatar
4 votes
Accepted

fast and stable x * tanh(log1pexp(x)) computation

OP points to a particular implementation of the mish activation function for accuracy specifications, so I had to characterize this first. That implementation uses ...
njuffa's user avatar
  • 520
4 votes
Accepted

How can I compute logarithm when comparison is undecidable?

Even though absolute comparisons may not converge, you should be able to narrow the argument into at least one of several partially overlapping ranges, such that you have a technique that works in ...
Pseudonym's user avatar
  • 22.1k
4 votes
Accepted

How to represent zero as floating point number?

Note: In the interest of making this somewhat self-contained, I am using terminology from the most recent versions of the IEEE-754 standard. Prior to 2008, "subnormal numbers" were called &...
Pseudonym's user avatar
  • 22.1k
3 votes

Alternatives to SVD for rank factorization

The proper search term in scientific journals is "Rank-Revealing Decomposition". If You want some theoretic guarantees on numeric accuracy/stability, the search term would be "Strong ...
DirkT's user avatar
  • 991
3 votes
Accepted

Why do I get different results from two calculation methods?

For $\exp(-10)$, more terms have to be computed since the 30-th term $\dfrac{(-10)^{30}}{30!}\approx 0.00377$ is several times bigger than the partial sum so far, 0.0009703416. For $\exp(10)$, all ...
John L.'s user avatar
  • 39k
3 votes

Proof that (x-y)(x+y) is more accurate than x²-y²

I'm not 100% sure of everything, but here are some elements that are implicit in the proof. Note that I reuse the notations and common knowledge of the paper, therefore this answer isn't self-...
Celelibi's user avatar
  • 453
3 votes

Proof that (x-y)(x+y) is more accurate than x²-y²

Suppose $y = x + \Delta$. Then $(x^2 - y^2) = x^2 - (x^2 + 2\Delta x + \Delta^2) = - (2\Delta x + \Delta^2)$ with leading term on the order of $2 \Delta x$. Multiply by $1 + \delta_1$ and that's still ...
Peter Taylor's user avatar
  • 2,082
3 votes

Why is computation of this function numerically unstable?

We can write it as $$f(x) = \dfrac{\sqrt{-x+a} - \sqrt{2x+a}}{ 4a}= \dfrac{\sqrt{-\frac xa+1} - \sqrt{\frac{2x}a+1}}{ 4\sqrt a}$$ If $x \ll a$ both of the square roots will be just about $1$. When ...
Ross Millikan's user avatar
3 votes
Accepted

Reception of numerical infinities

So the paper's "Numerical infinities and infinitesimals: Methodology, applications, and repercussions on two Hilbert problems", Yaroslav D. Sergeyev (2017), and it's basically a discussion on the ...
Nat's user avatar
  • 1,351
3 votes
Accepted

Avoiding overflow in computing the ratio of two large numbers

Yes, when $x\ge355$, computing $e^{2x}$ as a double-precision float in IEEE 754-1985 standard, leads to an overflow. $$x + \frac{e^x-x}{e^{2x} + 1}= x + e^{-x} - \frac{x+e^{-x}}{e^{2x} + 1}\approx x+ ...
John L.'s user avatar
  • 39k
3 votes

fast and stable x * tanh(log1pexp(x)) computation

There's no need to perform the logarithm. If you let $p = 1+\exp(x)$ then we have $f(x) = x\cdot\dfrac{p^2-1}{p^2+1}$ or alternatively $f(x) = x - \dfrac{2x}{p^2+1}$.
orlp's user avatar
  • 13.4k
3 votes
Accepted

Precise algorithm for finding higher order derivatives

The first thing you should understand is why central differencing gives you a more precise solution. Consider the Taylor expansion of $f$ around $x$: $$f(x + h) = f(x) + h f'(x) + \frac{1}{2} h^2 f''(...
Pseudonym's user avatar
  • 22.1k
2 votes

Big O notation: removing big O from denominator

Let us compute the difference. $$\begin{aligned} &\frac{\xi+\xi^2}{\xi-\frac{1}{2}\xi^2+\frac{1}{3}\xi^3+\mathcal{O}(\xi^{4})} -\frac{1+\xi}{1-\frac{1}{2}\xi+\frac{1}{3}\xi^2}\\ =&\frac{1+\xi}...
John L.'s user avatar
  • 39k
2 votes

Big O notation: removing big O from denominator

Let's start out simple: $${1 \over 1 + x} = 1 - x + x^2 - \cdots$$ so it follows that $${1 \over 1 + O(x)} = 1 + O(x)$$ as $x \to 0$. (Background: I interpret $O(x)$ as representing any function $...
D.W.'s user avatar
  • 159k
2 votes

Proof that a guard digit bound the error of subtraction

I'll post some elements of answer to my own question from what I understand. Anyone feel free to make a better, less sloppy answer. First, there are several versions of this document online. They all ...
Celelibi's user avatar
  • 453
2 votes

Rigorous error bounds for eigenvalue solvers

It might suffice to ask the solver to give you both the eigenvalue $\lambda$ and the corresponding eigenvector $v$. Then you can verify for yourself how much error there is. Note that if $\lambda$ ...
D.W.'s user avatar
  • 159k
2 votes

How much can we trust mathematical software when working with large numbers, and how much memory it needs to work with these numbers?

Question a): Here is the output from python console. ...
John L.'s user avatar
  • 39k
2 votes

Is order of matrix multiplication affecting numerical accuracy of the result?

Yes, there are differences in accuracy since with machine numbers the usual properties of arithmetics don't hold. Machine numbers are defined as $$ F(\beta,t,m,M)= \{ 0 \} \cup \{ x \in \mathbb{R} : ...
Gerardo Zinno's user avatar
2 votes

Robust two lines/segments intersection point in 2D

If you choose a formula with parameter $t$ you'll get the $P_x$ value as a result of a number of operations - additions, multiplications and one division: $$t = \frac{(x_1-x_3)(y_3-y_4)-(y_1-y_3)(x_3-...
HEKTO's user avatar
  • 3,088
2 votes

fast and stable x * tanh(log1pexp(x)) computation

My impression is that someone wanted to multiply x by a function f(x) that goes smoothly from 0 to 1, and experimented until they found an expression using elementary functions that did this, with no ...
gnasher729's user avatar
2 votes

fast and stable x * tanh(log1pexp(x)) computation

The context here is computer vision and the activation function for training neural nets. Chances are this code is going to be executed on a GPU. While performance is going to depend on the ...
Armadillo Jim's user avatar
2 votes

Validity of Algorithm for Testing Two Floating Point Numbers

There is one method to compare floating point numbers for equality, which is both very simple and correct: You use the equality (==) operator. There is another method to compare whether floating point ...
gnasher729's user avatar
2 votes
Accepted

Cancellation in C++

Let me summarize what I wrote in the comments. It is not a complete answer, since the intervals on which to apply each formula still need to be investigated. It is enough to assume $x\geq0$, since $\...
plop's user avatar
  • 1,189

Only top scored, non community-wiki answers of a minimum length are eligible