31
This is indeed called catastrophic cancellation. In fact, this particular case is very easy: rewrite the function using the equivalent, numerically stable expression $$ \frac{t}{1+\sqrt{1-t^2}}. $$
Since you probably need a reference, this is discussed in most numerical methods textbooks in relation to the formula for solving quadratic equations (that ...
27
We all know that
\begin{equation}
\exp(x) = \sum_{n=0}^\infty \frac{x^n}{n!} = 1 + x + \frac12 x^2 + \dots
\end{equation}
implies that for $|x| \ll 1$, we have $\exp(x) \approx 1 + x$. This means that if we have to evaluate in floating point $\exp(x) -1$, for $|x| \ll 1$ catastrophic cancellation can occur.
This can be easily demonstrated in python:
>&...
26
Are these integers or floating point numbers? Assuming it's floating point, I would go with the first option. It's better to add the smaller numbers to each other, then add the bigger numbers later. With the second option, you'll end up adding a small number to a big number as i increases, which can lead to problems. Here's a good resource on floating point ...
24
animal_magic's answer is correct that you should add the numbers from smallest to largest, however I want to give an example to show why.
Assume we are working in a floating point format that gives us a staggering 3 digits of accuracy. Now we want to add ten numbers:
[1000, 1, 1, 1, 1, 1, 1, 1, 1, 1]
Of course the exact answer is 1009, but we can't get ...
20
Intel support for IEEE float16 storage format
Intel supports IEEE half as a storage type in processors since Ivy Bridge (2013). Storage type means you can get a memory/cache capacity/bandwidth advantage but the compute is done with single precision after converting to and from the IEEE half precision format.
https://software.intel.com/content/www/us/en/...
18
Adding arbitrary floating point numbers will usually give some rounding error, and the rounding error will be proportional to the size of the result. If you calculate a single sum and start by adding the largest numbers first, the average result will be larger. So you would start adding with the smallest numbers.
But you get better result (and it runs ...
18
I think Higham's Accuracy and Stability of Numerical Algorithms addresses how one can analyze these types of problems. See Chapter 2, especially exercise 2.8.
In this answer I'd like to point out something that isn't really addressed in Higham's book (it doesn't seem to be very widely known, for that matter). If you are interested in proving properties of ...
17
Double precision is fairly common on newer GPUs. For instance I own a NVIDIA GTX560 Ti (fairly low end when it comes to computing) that has no issue running ViennaCL in double precision. From here (section 4) it appears all NVIDIA cards from GTX4xx onward support double precision natively.
I would guess that the GROMACS information is simply outdated.
15
For the general case, I'd use compensated summation (or Kahan summation). Unless the numbers are already sorted, sorting them will be much more expensive than adding them. Compensated summation is also more accurate than sorted summation or naive summation (see the previous link).
As for references, What every programmer should know about floating-point ...
15
This question is very closely related to (and possibly a duplicate of) Is variable scaling essential when solving some PDE problems numerically?.
There are still good practical reasons to nondimensionalize equations, if possible:
It reduces the number of independent parameters for parametric studies (which was one of the original reasons for ...
15
No, that is not guaranteed. If you are using a NETLIB BLAS without any optimizations, it it mostly true that the results are the same. But for any practical usage of BLAS and LAPACK one uses a highly optimized an parallel BLAS. The parallelization causes, even if it only works in parallel inside the vector registers of a CPU, that the order how the single ...
14
Here is R1, as computed in MATLAB:
1.0e+07 *
-7.382605957465515 -9.599867106092937 -2.830412177259742 -0.000000000002830 -0.000000000002830
-1.230434326244253 -1.599977851015490 -0.471735362876624 -0.000000000000472 -0.000000000000472
3.691302978732758 4.799933553046468 1.415206088629871 0.000000000001415 0.000000000001415
-5....
14
First, see Mark L. Stone's answers, which is completely correct. Second, realize that this is the reason why people told you to use relative errors in your numerical analysis class. :)
Third, the real question here is why the results do not coincide exactly, since both languages call some BLAS library functions for their computations. There are several very ...
13
The formula
$$\mathrm{logsum}(x,y)=\max(x,y)+\mathrm{log1p}(\exp(-\operatorname{abs}(x-y))$$
should be numerically stable. It generalizes to a numerically stable computation of
$$\log \sum_i e^{x_i} = \xi+ \log\sum_i e^{x_i-\xi},~~~\xi=\max_i x_i$$
In case the logsum is very close to zero and you want high relative accuracy, you can probably use
$$\...
12
Use the (IEEE standard) library function log1p, which should be present in all programming languages. The function log1p(x) returns $\log(1+x)$, and is implemented with particular attention to accuracy when $x$ is small. It is designed to solve exactly this kind of problem.
11
You can directly check for denormals using isnormal() from math.h (C99 or later, POSIX.1 or later). In Fortran, if the module ieee_arithmetic is available, you can use ieee_is_normal(). To be more precise about fuzzy equality, you have to consider the floating point representation of denormals and decide what you mean for the results to be good enough.
More ...
11
$p \log p$ won't suffer from precision loss anywhere in $[0,1]$, and won't suffer from exponent overflow near $0$ either. Thus, the fast, accurate way is
p ? p * -log(p) : 0
11
The theoretical peak FLOP/s is given by:
$$ \text{Number of Cores} * \text{Average frequency} * \text{Operations per cycle} $$
The number of cores is easy. Average frequency should, in theory, factor in some amount of Turbo Boost (Intel) or Turbo Core (AMD), but the operating frequency is a good lower bound. The operations per cycle is architecture-...
11
In the 1980's era Intel 80x86 architecture, there was a scalar floating point unit that had instructions like FSIN, FCOS, etc. for computing functions like sin and cos. These functions were implemented in microcode and might take 100's of CPU cycles to execute.
Later, Intel added Streaming SIMD Extensions (SSE) which gave the processor parallel floating ...
10
The point they are making is that if you compute the difference
123457.1467 - 123456.659 = 0.4877
then you get an answer that seems to be quite different than the difference
123457.1 - 123456.7 = 0.4
that you get by subtracting the rounded numbers. The ~20% figure is because 0.4877 is more than 20% larger than 0.4.
This example shows that if you subtract ...
10
Some references on rounding error analysis of Krylov methods:
Gutknecht, Martin H., and Zdenvek Strakos. "Accuracy of two three-term and three two-term recurrences for Krylov space solvers." SIAM Journal on Matrix Analysis and Applications 22.1 (2000): 213-229.
Paige, Christopher C., and Zdenvek Strakos. "Residual and backward error bounds in minimum ...
10
If you're looking for a good introduction, I would suggest David Goldberg's What Every Computer Scientist Should Know About Floating Point Arithmetic. It may be a bit too detailed, but it's available online for free.
If you've got a good library, I'd suggest Michael Overton's Numerical Computing with IEEE Floating Point Arithmetic, or the first few chapters ...
10
The numeric precision is not perfect. You get rounding errors during your computation. When working with floats, don't check if they are = 0, but check if their absolute distance to 0 is smaller than some epsilon.
9
Their observation ''the form of the algorithm may needlessly introduce some numerical imprecision'' is correct. But their explanation
''This arises because, in eqn (3.14), a small term ($O(δt^2)$) is added to a difference of large terms ($O(δt^0)$), in order to generate the trajectory.''
is spurious.
The true reason for the slight numerical instability of ...
9
Horner is indeed the most stable way to evaluate a polynomial (and you get the bonus of evaluating its derivatives with not too much extra cost). Higham presents a nice error analysis of the algorithm (Accuracy and stability of numerical algorithms, 2nd edition, p.94). He also presents an algorithm that includes a running error bound so you have an idea on ...
9
There are aspects of modern computing systems that are inherently non-deterministic that can cause these kinds of differences. As long as the differences are very small in comparison with the required accuracy of your solutions, there probably isn't any reason to worry about this.
An example of what can go wrong based on my own experience. Consider the ...
9
You can sometimes prove such results (or get counterexamples) using an SMT solver such as Z3 that supports floating point arithmetic. Here is a proof of a version of your theorem that says $|((x+y)-y)-x| \leq 2^{-23}|x|$ when $x>y>1$ and $x+y\neq\infty_{32}$ in 32-bit floating point arithmetic:
λ> import Data.SBV
λ> :set -XScopedTypeVariables
λ&...
8
The Matlab command help eps says the following:
D = EPS(X), is the positive distance from ABS(X) to the next larger in
magnitude floating point number of the same precision as X. X may be
either double precision or single precision.
In other words, if $\varepsilon_\mathsf{mach}$ is the relative error due to floating point, as defined in the Wikipedia ...
8
Let's denote by $\otimes,\oplus,\ominus$ (I was lazy trying to get circled version of division operator) the floating-point analogs of exact multiplication ($\times$), addition ($+$), and subtraction ($-$), respectively. We'll assume (IEEE-754) that for all of them
$$
[x\oplus y]=(x+ y)(1+\delta_\oplus),\quad |\delta_\oplus|\le\epsilon_\mathrm{mach},
$$
...
Only top voted, non community-wiki answers of a minimum length are eligible
