32 votes
Accepted

How to avoid catastrophic cancellation in python function?

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}}. $$ ...
user avatar
  • 11.4k
29 votes
Accepted

When should log1p and expm1 be used?

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 ...
user avatar
  • 3,789
25 votes
Accepted

Is half precision supported by modern architecture?

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 ...
user avatar
  • 1,991
22 votes

Meaning of "-0.0" in Python?

Floating point numbers (according to the standard1 nearly all programming languages use) are stored with a certain number of bits in the mantissa, in the exponent, and with a sign bit. As such, ...
user avatar
20 votes
Accepted

Why are log and exp considered 'expensive' computations in ML?

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 ...
user avatar
20 votes
Accepted

What are some good strategies to test a floating point arithmetic implementation for double numbers?

You should test transition points. Floating-point numbers have several distinct "ranges": Standard/Normal arithmetic Subnormal arithmetic Infinite arithmetic NaN arithmetic Zero arithmetic ...
user avatar
  • 3,091
18 votes

Robust computation of the mean of two numbers in floating-point?

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 ...
user avatar
  • 11.4k
15 votes
Accepted

Why should I renormalize physical variables?

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 ...
user avatar
15 votes
Accepted

How to determine the amount of FLOPs my computer is capable of

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 ...
user avatar
15 votes
Accepted

Are BLAS implementations guaranteed to give the exact same result?

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 ...
user avatar
14 votes

Matrix multiplication accuracy Matlab vs Python

Here is R1, as computed in MATLAB: ...
user avatar
14 votes

Matrix multiplication accuracy Matlab vs Python

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 ...
user avatar
12 votes
Accepted

Hardware performance, floating point functions

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 ...
user avatar
12 votes
Accepted

Evaluating $\log(\exp(x)+1)$ for negative $x$

Use the (IEEE standard) library function log1p, which should be present in all programming languages. The function log1p(x) ...
user avatar
11 votes

Why are log and exp considered 'expensive' computations in ML?

$\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....
user avatar
  • 3,511
10 votes

Why is $\exp(\ln(x))-x\neq0$ in floating point arithmetic?

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 ...
user avatar
10 votes
Accepted

Accuracy loss in single-precision Euclidean norm computation

The squares are harmless (as long as they don't overflow/underflow), because the relative perturbation they introduce is of the order of the machine precision $u\approx 10^{-8}$. Your troubles here ...
user avatar
9 votes
Accepted

Accurate Polynomial Evaluation in Floating Point

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 ...
user avatar
  • 6,046
9 votes
Accepted

Small, unpredictable results in runs of a deterministic model

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 ...
user avatar
9 votes
Accepted

Is it possible to proof a-b+b = a for all double floating-point numbers?

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)-...
user avatar
  • 11.4k
8 votes

Robust computation of the mean of two numbers in floating-point?

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 ...
user avatar
  • 181
8 votes
Accepted

Order of operations, numerical algorithms

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 (...
user avatar
  • 8,287
8 votes
Accepted

Stabilizing a 3x3 real symmetric matrix eigenvalue calculation

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 ...
user avatar
  • 11.4k
8 votes
Accepted

Can floating point error (in FFTW3) cause non-deterministic behavior?

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, ...
user avatar
8 votes

Are BLAS implementations guaranteed to give the exact same result?

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 ...
user avatar
  • 1,522
8 votes
Accepted

Integer operations vs floating point operations

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 ...
user avatar
  • 8,287
8 votes

Is half precision supported by modern architecture?

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 ...
user avatar
  • 4,286
8 votes
Accepted

Polynomial approximation for floating-point arithmetic

The sine is an odd function, so you want that also in an approximation. A polynomial with $p(0)=0$ can be factored as $p(x)=xq(x)$, so $q(x)\approx \frac{\sin(x)}{x}$. Each interval $[2^n,2^{n+1})$ ...
user avatar
  • 3,511
7 votes
Accepted

Choosing epsilons

The Wikipedia article cite Numerical Recipes, 3rd edition, Section 5.7, which is pages 229-230 (a limited number of page views is available at http://www.nrbook.com/empanel/). Sure enough, they ...
user avatar
7 votes

Why is $\exp(\ln(x))-x\neq0$ in floating point arithmetic?

[EDIT] An alternate view: 64-bit floating numbers represent a discrete set $S$. For a function $f$ to be exactly invertible, it should be a bijection from $S$ to $S$. Suppose we are interested in a ...
user avatar

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