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 ...
28
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 ...
21
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/...
19
To add to Lutz Lehmann's answer, you can look up the latency for the CPU instructions in this comprehensive table by Agner Fog.
For example, on the Intel Ivy Bridge processors:
FADD / FSUB (floating point add and subtract) both take 3 cycles
FMUL (multiply) takes 5 cycles
FDIV (divide) takes 10-24 cycles
FYL2X ($y \cdot \log_2(x)$) takes 90-106 cycles
F2XM1 ...
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 ...
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
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
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 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-...
12
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 ...
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
$\exp$, $\sin$, $\tan$ and their inverse and otherwise related functions are transcendental, defined by an infinite power series. Meaning it takes some effort to evaluate uniformly good approximations. This is done via argument simplification using the properties of these functions, and polynomial approximations, or using quotients of polynomials à la Padé ...
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
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},
$$
...
8
First, observe that if you have a method that gives a most accurate answer in all cases,
then it will satisfy your required condition.
(Note that I say a most accurate answer rather than the most accurate answer, since there may be two winners.)
Proof: If, to the contrary, you have an accurate-as-possible answer that does not satisfy the
required condition, ...
8
This is trying to compute the eigenvalues by computing the roots of the characteristic polynomial. In this case, the characteristic polynomial is $p(t) = t^3-2t^2x$, $x=1.25\times 10^6$, and zero is a root of multiplicity two.
There is generally no good way of computing double roots of a polynomial to precision better than $O(\sqrt{\epsilon})$, because if ...
8
Non-reproducible behaviors in computing amidst different runs can involve several mechanisms, sometimes mixed. They can be especially sensitive when one iterates calculations on large sets of data, like inverse 3D tomography.
soft errors, caused by defects, interaction with high-energy particles in spacecrafts (let's rule them out for now)
number ...
8
The Short Answer
If the two BLAS implementations are written to carry out the operations in the exact same order, and the libraries were compiled using the same compiler flags and with the same compiler, then they'll give you the same result. Floating point arithmetic is not random, so two identical implementations will give identical results.
However, ...
8
There is a nice discussion on StackOverflow regarding floating point vs integer operations. In short, the performance of the operations depends a lot on
processor architecture
how the data is stored in memory and in which order it is accessed
if (and which) SSE/AVX/AVX2/etc instructions are used (and how efficiently)
This probably provides some insight ...
8
In my opinion, not very uniformly. Low precision arithmetic seems to have gained some traction in machine learning, but there's varying definitions for what people mean by low precision. There's the IEEE-754 half (10 bit mantissa, 5 bit exponent, 1 bit sign) but also bfloat16 (7 bit mantissa, 8 bit exponent, 1 bit sign) which favors dynamic range over ...
7
It's worth a bit of a discussion as to why denormals don't matter in typical scientific computation. Usually, we are given some input $x$ and want to compute a function $f(x)$. Due to numerical errors, we instead compute $\tilde{f}(x) \approx f(x)$. There are two kinds of accuracy typically sought:
Relative: $$|f(x)-\tilde{f}(x)| < \alpha |f(x)|$$
...
Only top voted, non community-wiki answers of a minimum length are eligible
