25

Take \begin{align} 1-\sqrt{ 1-x^2} &= (1-\sqrt{ 1-x^2})\frac{1+\sqrt{ 1-x^2}}{1+\sqrt{ 1-x^2}}\\ &= \frac{x^2}{1+\sqrt{ 1-x^2}} \end{align} So \begin{align} y = x\sqrt{\frac{1}{2+2\sqrt{1-x^2}}} \end{align}


21

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


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


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

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


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

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


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

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

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


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

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


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


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

What is typically done for ODEs is to solve at each timestep using two methods of different order. The discrepancies between the results are then used to estimate the error and then decide whether to decrease or increase the timestep length. You can find more information on the subject here. Two examples of these types of schemes are RKF45 and Dormand-...


5

for 1-D transport problem, under implicit method, we use courant number, a dimensionless number to choose appropriate time steps: $\frac{u\Delta t}{\Delta x}\le C_{max}$ $C_{max}$ should be less than 1, if we make it as 1, we get $\Delta t\le \frac{\Delta x}{u}$ The basic ideal is that, at one time step, a particle placed in the problem domain should ...


5

In contrast to Bill Barth, I usually try to keep things in dimensional form. Within a single equation, this of course does not change the relative scaling of terms. However, not doing the scaling requires that one pays attention to the relative scaling between different equations of a system of equations. A discussion of one case we have documented can be ...


5

I think there's a simple way to do this. You have a rational function of identical cosh/sinh terms, where every expression is a homogeneous polynomial in cosh/sinh, and the only problem is that these exponential terms overflow. The function does not diverge as these terms approach infinity, so if you divide every numerator and denominator by the same power ...


5

Orthogonal matrices are about as well-conditioned as you can get, but numerical errors still occur. One common error is loss of orthogonality. A fix for this could be to re-orthogonalize your columns after some number of multiplications. You can do this by just taking the QR decomposition of your matrix after some number of products and taking the orthogonal ...


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