How to prevent overflow and underflow in the Euclidean distance and Mahalanobis distance

Question

I was working in my project when I was struck by the question of whether it would be necessary, or at least cautious, prevent overflow and underflow in the calculation of these two distances.

I remembered that there is an implementation of the calculation of the hypotenuse to prevent this. Most languages ​​implementers, and is known for Hypot

The calculation of the Euclidean distance remains the same "pattern" and I thought that if Hypot() controls the overflow and underflow should also beware of the Euclidean distance. I've disappointed to note that the language we use, and others, do not control the overflow and underflow for calculating distance. Will not worth spend this "additional effort"?

I did a searchs and came to a question in Math.StackExchange

There is no definitive answer to this issue and is somewhat old. The first thing I wondered is: Will okay? I think that yes, seeing that is a generalization of the same procedure that performs Hypot().

I decided to extrapolate this concept to the Mahalanobis distance. The original is as follows:

$$D_M(X,Y,L) = \sqrt{\sum_{i=1}^{n}\left(\frac{X_i-Y_i}{L_i}\right)^2}$$

Since $L$ is the vector of eigenvalues.

And my proposal is this:

$$D_M(X,Y,L) = C\sqrt{\sum_{i=1}^{n}\left(\frac{X_i-Y_i}{L_i }\frac{1}{C}\right)^2}$$

That is the same to:

$$D_M(X,Y,L) = C\sqrt{\sum_{i=1}^{n}\left(\frac{X_i-Y_i}{L_i C}\right)^2}$$

And $C$ is the max value from the $|(X_i-Y_i)/L_i|$:

$$C = \max_{i}\left(|\frac{X_i-Y_i}{L_i}|\right)$$

Is it okay?

Answer 1

That is a long solved problem. First, if you are using double precision floating point numbers in IEEE754 format (which is most common), that's what extended precision was invented for. Even in the extreme case, $X_i$ and $Y_i$ huge, $L_i$ tiny, the numbers involved are less than $10^{1300}$ and extended precision goes up to about $10^{4000}$.

Without extended precision available, or if your numbers are already in extended precision, you just calculate the result and check whether you had an overflow or if the result is tiny. Tiny results are critical as well: If you use double precision floating point numbers, and your numbers are around $10^-200$, then squaring them gives you zero. There will be a lower bound c where you can prove "if the result is greater than c, then any underflow in the calculation didn't change the final result".

If overflow happens, or underflow happens, then you scale the numbers by a suitable power of two. Do NOT divide by the maximum value, because that introduces rounding errors, and there are cases with n = 2, $L_i = 1$, where $|X_1 - Y_1| ≤ |x_1 - y_1|$, $|X_2 - Y_2| ≤ |x_2 - y_2|$, yet the calculated distance of X, Y is greater than the calculated distance of x, y, due to rounding errors. Scaling by a power of two avoids the problem because it doesn't create any rounding errors.

And you want to make sure that the order of numbers doesn't make a difference. For n = 2 you must make sure that no fused multiply/add operations are used. For n ≥ 3 you probably want to sort the numbers so the ones with smallest absolute difference come first. Otherwise you can get into situations where the distance between $(X_1, X_2)$ and $(Y_1, Y_2)$ is not the same as the distance between $(X_2, X_1)$ and $(Y_2, Y_1)$.

Answer 2

I don't know whether your approach is correct, but to answer the title, one way to avoid overflow is using approximations, for example (copied from here):

$$\frac{1007}{1024}\max(|x|,|y|) + \frac{441}{1204}\min(|x|,|y|)$$

Note that this should work fine to avoid overflow from large values. However, I think the approximation will be quite bad for high-dimensional cases (10 dimensions or more), which is another cause of overflows.