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/...


19

(I have tested this approach before, and I remember it worked correctly, but I haven't tested it specifically for this question.) As far as I can tell, both $\|\mathbf{v}_1\times \mathbf{v}_2\|$ and $\mathbf{v}_1\cdot \mathbf{v}_2$ can suffer from catastrophic cancellation if they are almost parallel/perpendicular—atan2 can't give you good accuracy if ...


16

Discretization of the continuum PDEs usually commits much more error than the finite precision. I find that about 90% of the people requesting higher precision have just been lazy with problem formulation and are trying to solve a problem using poor scaling, bad discretizations, or bad continuum modeling. The remaining 10% may have justifiably ill-...


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 ...


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

I am not an expert in machine learning, but I can outline the considerations that are relevant. The numerical calculations in machine learning are generally linear algebra -- either solving linear systems or linear least squares. For both types of problems, there are well-known backward-stable methods, so I will assume you are using a backward-stable ...


10

For $a+b\sqrt{-3}$ you can use the representation $$\begin{pmatrix}a&-3b\\b&a\end{pmatrix}$$ Addition works obviously. For multiplication, you can verify $$\begin{pmatrix}a_1&-3b_1\\b_1&a_1\end{pmatrix}\begin{pmatrix}a_2&-3b_2\\b_2&a_2\end{pmatrix} = \begin{pmatrix}a_1a_2-3b_1b_2&-3(a_1b_2+b_1a_2)\\a_1b_2+b_1a_2&a_1a_2-3b_1b_2\...


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.


10

Jean-Michel Muller, et. al., "Handbook of Floating-Point Arithmetic 2nd ed.", Birkhäuser 2018, gives the following example due to Muller, specifically constructed to deliver incorrect results with floating-point evaluation: $$ {u_{0} = 2,\\ u_{1} = -4,\\ u_{n} = 111 - \frac{1130}{u_{n-1}} + \frac{3000}{u_{n-1}u_{n-2}},\>\>\>\>n \ge 2.} $$ ...


9

In the 15 years that we have provided FEM software in the form of the deal.II project (http://www.dealii.org/), I don't think that we've ever had a genuine request to solve PDEs to higher accuracy than double precision. The reason is as Jed suggests in the other answer: The error one makes discretizing the PDE is much larger than the 16 digits of accuracy ...


9

Overview Good question. There is a paper entitled "Improving the accuracy of the matrix differentiation method for arbitrary collocation points" by R. Baltensperger. It's no big deal in my opinion, but it has a point (that already was known before the appearance in 2000): it stresses the importance of an accurate representation of the fact that the ...


9

I would like to know if there are any practical methods for obtaining an approximation to $\mathbf{x}(t_f)$ (where $t_f \in \mathbb{R}$ is some given final time) which is provably correct to $N$ [sic] digits. That all depends on your opinion of the practicality of interval arithmetic. There are validated integrators available, such as the COSY code out of ...


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

The efficient answer to this question is, not too surprisingly, in another note by Velvel Kahan: $$\alpha=2\arctan\left(\left\|\frac{\mathbf v_1}{\|\mathbf v_1\|}+\frac{\mathbf v_2}{\|\mathbf v_2\|}\right\|,\left\|\frac{\mathbf v_1}{\|\mathbf v_1\|}-\frac{\mathbf v_2}{\|\mathbf v_2\|}\right\|\right)$$ where I use $\arctan(x,y)$ as the angle made by $(x,y)$ ...


7

Calculate the SVD in place of the spectral decomposition. The results are the same in exact arithmetic, as your matrix is symmetric positve definite, but in finite precision arithmetic, you'll get the small eigenvalues with much more accuracy. Edit: See Demmel & Kahan, Accurate Singular Values of Bidiagonal Matrices, SIAM J. Sci. Stat. Comput. 11 (1990)...


7

For all non-trivial problems (i.e., for those where performance matters) almost all of the memory you have will be in the matrix, and relatively little in vectors. For example, for 3d Taylor-Hood elements for the Stokes equation, you have a few hundred elements per row in the matrix, and this vastly outweighs the amount of memory needed for vectors. We have ...


7

It is the difference in scales between terms in your equations that tend to cause numerical difficulties. You may work in any units you like as long as you are consistent. My approach has been to always consistently non-dimensionalize my equations in order to reduce the number of parameters to the minimum required, but this is only for my convenience.


7

At least in MATLAB, I believe abs(z) is implemented as sqrt(z*z'). The extra square-root and squaring operation reduces numerical precision. >> z = randn + randn * i z = 0.5377 + 1.8339i >> abs(z)^2 - z*z' ans = 4.4409e-16 >> abs(z)^2 - sqrt(z*z')^2 ans = 0


7

[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 fast growing function like $\exp$. At one point, $\exp x > x$. If $M$ is the maximum element from $S$, then $f(M)\ge M$. The strict version $f(M)> M$ is ...


7

Added after my initial answer: It appears to me now that the author of the referenced paper is giving condition numbers (apparently 2-norm condition numbers but possibly infinity-norm condition numbers) in the table while giving maximum absolute errors rather than norm relative errors or maximum elementwise relative errors (these are all different ...


7

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 ...


6

This may be a bit outdated, but Hairer and Wanner's book recommends their own radau5 code. In my own experience, this code is both extremely robust and efficient. It is also a rather straight-forward Runge-Kutta scheme, and thus not too difficult to implement. There's a Matlab implementation of radau5 by Christian Engstler somewhere out there.


6

There are many good open-source implementations of root-finding methods out there already. One example is boost, whose implementation of Newton-Raphson and related methods you can find here. If you read its source code, you will be able to see what issues the authors wanted to address, and how they dealt with issues like convergence, user-specified tolerance ...


6

There are three issues that are likely to cause such problems in pseudospectral methods: Gibbs oscillations Aliasing Time step too large In any case you likely develop oscillations in the solution until some point ends up with a negative density, resulting in a NaN when computing the pressure or sound speed or some other term. The solution to 3 is obvious, ...


6

Yes, it is possible to show that the statistical behavior of the approximate system will reach that of the "exact" system. (This is true even though hard-sphere dynamics do not accurately describe molecular systems!) The basic premise underlying molecular dynamics is the ergodic theorem, which states that, in the limit of long times, the time average of a ...


6

In my opinion, an example could be the calculation of the minimal solution of a three term recurrence relation (TTRR) $$y_{n+1} +a_{n}y_{n}+b_{n}y_{n-1} = 0, \quad n=1,2,3,\ldots$$. For example, the Bessel $J_{n}(x)$ function for a fixed $x$ satisfies the TTRR $$y_{n+1}-\frac{2n}{x}y_{n}+y_{n-1} = 0$$ Suppose you know $y_{0} = J_{0}(1)$ and $y_{1} = J_{1}(...


6

The accepted answer provides an overview. I'll add a few more details about support in NVIDIA processors. The support I'm describing here is 16 bit, IEEE 754 compliant, floating point arithmetic support, including add, multiply, multiply-add, and conversions to/from other formats. Maxwell (circa 2015) The earliest IEEE 754 FP16 ("binary16" or "half ...


5

The computation of integer-valued matrix determinants has been a subject of considerable research. Using exact arithmetic the Smith normal form can be computed, and from this diagonal form the determinant is easily found. Saunders and Wan (2004), Smith Normal Form of Dense Integer Matrices, Fast Algorithms into Practice, say "Over the past thirty years, ...


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