# How do you *accurately* calculate distances in the hyperboloid model of hyperbolic geometry?

In the hyperboloid model, the distance between two points has a straightforward calculation: $$D(\mathbf{u}, \mathbf{v}) = \operatorname{arcosh}(-B(\mathbf{u}, \mathbf{v})) = \operatorname{arcosh}(\mathbf{u}_0 \mathbf{v}_0 - \mathbf{u}_1 \mathbf{v}_1 - \ldots - \mathbf{u}_n \mathbf{v}_n)$$

However, implementing this definition in a straightforward manner leads to terrible numerical issues. For starters, as $$D \rightarrow 0$$, the expression inside goes to 1, losing most of its precision. Here's a straightforward Python implementation showing the issues in the case of 3 dimensions (with the zeroth coordinate implicitly computed from the other 3):

def dist_bad(m, n):
partials = [math.hypot(*m, 1) * math.hypot(*n, 1),
-m[0]*n[0], -m[1]*n[1], -m[2]*n[2]]
return math.acosh(math.fsum(partials))

print(dist_bad([0,0,0], [1.5e-8,0,0])

2.1073424255447017e-08
0.0


(The correct value is ~1.5e-8 and ~1.4e-8, respectively).

It also fails badly for reasonable distances when away from the origin:

print(dist_bad([1e6,0,0], [1e6,0.02,0])

0.022096637370389118


(The correct value is ~0.01999966)

What's a better transformation of the distance formula that preserves precision for small distances?

• it is just the numerical stability issue, you can use an expansion of acosh to compute. Apr 15 at 2:53
• I'm still thinking about this. You can reduce your second problem (large aspect ratio between components with a large component) by scaling the largest component, $\max($*m$,1)$, so that this 4-vector is $m_{\text{max}}(m^1/m_{\text{max}}, m^2/m_{\text{max}}, m^3/m_{\text{max}}, 1/m_{\text{max}})$. Do the same to $n$, then chase the maxs through your distance formula separately from the components. They're repeated common factors, so this isn't so awful. So that reduces several general bad cases to your initial bad case. :-) Apr 15 at 3:19
• When I have worked something out for what to do with these "normalized"-ish 4-vectors, I'll post an answer. If I don't, ... Apr 15 at 3:20

The biggest source of precision loss is how the value inside of $$\operatorname{arcosh}$$ goes to 1. We'll use the definition of $$\operatorname{arcosh}(x) = \log(x + \sqrt{x^2-1})$$. Making the substitution $$t = x - 1$$, this becomes $$\log(1 + t + \sqrt{t(t+2)})$$, and we can compute this expression much more precisely using the standard log1p() function.

There is still the issue of catastrophic cancellation in the calculation of $$B$$. To handle this, note that $$t = x - 1 = \sqrt{(\mathbf{u}_1^2 + \ldots + \mathbf{u}_n^2 + 1)(\mathbf{v}_1^2 + \ldots + \mathbf{v}_n^2 + 1)} - \mathbf{u}_1 \mathbf{v}_1 - \ldots - \mathbf{u}_n \mathbf{v}_n - 1$$. If we multiply and divide by the conjugate, we can get terms to cancel, and also the conjugate will not have any cancellation (for points that are close to each other).

Specifically, the numerator becomes \begin{align} &(\mathbf{u}_1^2 + \ldots + \mathbf{u}_n^2 + 1)(\mathbf{v}_1^2 + \ldots + \mathbf{v}_n^2 + 1) - (\mathbf{u}_1 \mathbf{v}_1 + \ldots + \mathbf{u}_n \mathbf{v}_n + 1)^2\\ =\ &\left((\mathbf{u}_1^2 + \ldots + \mathbf{u}_n^2)(\mathbf{v}_1^2 + \ldots + \mathbf{v}_n^2) + (\mathbf{u}_1^2 + \ldots + \mathbf{u}_n^2) + (\mathbf{v}_1^2 + \ldots + \mathbf{v}_n^2) + 1\right)\\&\ - \left((\mathbf{u}_1 \mathbf{v}_1 + \ldots + \mathbf{u}_n \mathbf{v}_n)^2 + 2(\mathbf{u}_1 \mathbf{v}_1 + \ldots + \mathbf{u}_n \mathbf{v}_n) + 1\right)\\ =\ &(\mathbf{u}_1^2 + \ldots + \mathbf{u}_n^2)(\mathbf{v}_1^2 + \ldots + \mathbf{v}_n^2) - (\mathbf{u}_1 \mathbf{v}_1 + \ldots + \mathbf{u}_n \mathbf{v}_n)^2 + (\mathbf{u}_1 - \mathbf{v}_1)^2 + \ldots + (\mathbf{u}_n - \mathbf{v}_n)^2\\ =\ &\sum_{i,j}\mathbf{u}_i^2 \mathbf{v}_j^2 - \sum_{i}\mathbf{u}_i^2 \mathbf{v}_i^2 - \sum_{i}\sum_{j \ne i}\mathbf{u}_i \mathbf{v}_i \mathbf{u}_j \mathbf{v}_j + \sum_{i}(\mathbf{u}_i - \mathbf{v}_i)^2\\ =\ &\sum_{i}\sum_{j \ne i}\left(\mathbf{u}_i^2 \mathbf{v}_j^2 - \mathbf{u}_i \mathbf{v}_i \mathbf{u}_j \mathbf{v}_j\right) + \sum_{i}(\mathbf{u}_i - \mathbf{v}_i)^2\\ =\ &\sum_{i}\sum_{j > i}(\mathbf{u}_i \mathbf{v}_j - \mathbf{u}_j \mathbf{v}_i)^2 + \sum_{i}(\mathbf{u}_i - \mathbf{v}_i)^2\\ \end{align}

This confines the cancellation to a much smaller domain, dealing with many cross-product terms instead of the entire sum. We can test this implementation in 3 dimensions in Python:

def dist_good(m, n):
c1 = m[0]*n[1] - m[1]*n[0]
c2 = m[0]*n[2] - m[2]*n[0]
c3 = m[1]*n[2] - m[2]*n[1]
s1 = m[0] - n[0]
s2 = m[1] - n[1]
s3 = m[2] - n[2]
numer = [c1*c1, c2*c2, c3*c3, s1*s1, s2*s2, s3*s3]
denom = [math.hypot(*m, 1) * math.hypot(*n, 1),
m[0]*n[0], m[1]*n[1], m[2]*n[2], 1]
v = math.fsum(numer)/math.fsum(denom)
return math.log1p(v + math.sqrt(v * (v+2)))

print(dist_good([0,0,0], [1.5e-8,0,0])
print(dist_good([0,0,0], [1.4e-8,0,0])
print(dist_good([1e6,0,0], [1e6,0.02,0])

1.5e-8
1.4e-8
0.01999966668166577


Edit: Testing with the following bc (the arbitrary-precision calculator) program shows that the results are accurate to machine precision:

#!/usr/bin/bc -l

define dist(m0, m1, m2, n0, n1, n2) {
auto x;
x = sqrt((m0*m0 + m1*m1 + m2*m2 + 1)*(n0*n0 + n1*n1 + n2*n2 + 1)) - m0*n0 - m1*n1 - m2*n2;
return l(x + sqrt(x*x - 1));
}

scale=50
dist(1000000, 0, 0, 1000000, 0.02, 0);
quit;

.01999966668166577287033322057922018961564473589608

• Instead of using $$D=\operatorname{arcosh}(-\mathbf u\cdot\mathbf v)=\operatorname{arcosh}\left(u_0v_0-\sum_iu_iv_i\right),$$ you could try either of these (but I haven't tried them myself): $$D=\operatorname{arsinh}\lVert\mathbf u\wedge\mathbf v\rVert=\operatorname{arsinh}\sqrt{\sum_i(u_0v_i-u_iv_0)^2-\sum_i\sum_{j>i}(u_iv_j-u_jv_i)^2}$$ $$D=2\operatorname{arsinh}\tfrac12\lVert\mathbf u-\mathbf v\rVert=2\operatorname{arsinh}\frac12\sqrt{-(u_0-v_0)^2+\sum_i(u_i-v_i)^2}$$ Apr 21 at 23:54
• The real issue is dealing with the cancellation of $u_0$ and $v_0$ terms. AFAICT, this is easiest to work with in the arcosh form, and (if we want) we can transform to use asinh directly by using $\mathbb{arsinh}(\sqrt{t(t+2)})$ instead of log1p. Apr 22 at 11:28
• No, $u_0-v_0$ cancellation shouldn't be a problem. Geometrically, if $\mathbf u$ and $\mathbf v$ are near the origin $(1,0,0\cdots)$, the distance on the hyperboloid ($u_0=\sqrt{1+u_1^2+u_2^2+\cdots}$) is about the same as the distance on the plane ($u_0=1$); there's no great loss in rounding $u_0-v_0$ to $0$. -- Now I have tested the $\operatorname{arsinh}$ formulas. They do work better than the $\operatorname{arcosh}$ formula. The last one (with the factors of $2$) gives the same results as your dist_good. Apr 25 at 0:35
• (Of course I meant "distance on the hyperboloid" using the Minkowski metric.) Apr 25 at 0:45
• I didn't speak with enough clarity. What I meant was "rearranging the formula so that the $u_0$ and $v_0$ terms cancel out". If you don't do that, then (taking the 2nd asinh formula for example) the $(u_0 - v_0)^2$ term has similar magnitude to the sums term, and you get catastrophic cancellation. I'm not sure how you implemented the formulas, but for me when I tested the straightforward implementation of the 2nd arsinh one, I got 1.4999999999999995e-08 for [1.5e-8,0,0] and 12.206072645562394 for the distance from [0,0,0] to [1e5,0,0], which is off in the last 13 bits. Apr 25 at 21:50

Method I: Prior to converting coordinates to floating point, apply an isometry to move one to the origin, then measure the distance as precisely as you like. Possibly using ideas from your Answer or this one.

Method II: The problem initially is loss of precision for arguments to acosh near $$1$$. This is tricky because Python's floating point is implemented as C99's double which has maximal precision of 53 bits around $$1$$. For instance, $$\frac{\mathrm{arcosh}(1+2^{-53}) - \mathrm{arcosh}(1)}{2^{-53}} = 1.3{\dots}\cdot 10^{8} \text{.}$$ Thinking of this differently, rounding errors in the last bit of a double input to acosh are amplified about $$10^8 \approx 2^{26}$$ in the output. That is, a least significant bit change in the input's mantissa alters half the bits in the output's mantissa.

Next, there is catastrophic cancellation (which I think you have mentioned in the Question or in a comment to it or to your Answer). Suppose $$\vec{m} = \vec{n} + \vec{\varepsilon}$$, where as in your first examples, we think of $$\vec{m}$$ and $$\vec{n}$$ as small and of $$\vec{\varepsilon}$$ as very small. Then the angle, $$\theta$$, between $$\vec{m}$$ and $$\vec{n}$$ ranges from $$\theta = 0$$, when $$\vec{m}$$ and $$\vec{n}$$ are parallel, to $$\theta = \arctan(|\vec{\varepsilon}|/|\vec{m}|) = 1/\sqrt{1+|\vec{\varepsilon}|^2/|\vec{m}|^2}$$, when $$\vec{m}$$ and $$\vec{\varepsilon}$$ are orthogonal (think of a ball centered at $$\vec{n}$$ of radius $$|\vec{\varepsilon}|$$ and that the line spanned by $$\vec{m}$$ opens the widest angle with $$\vec{n}$$ when that line is tangent to the ball). Also $$|\vec{m}|-|\vec{\varepsilon}| \leq |\vec{n}| \leq |\vec{m}|+|\vec{\varepsilon}|$$. In fact, there is a relation between $$\theta$$ and where $$|\vec{n}|$$ is in this interval, but we allow our estimate to be a bit pessimistic by ignoring that relation. With these, we have \begin{align*} \sqrt{1 + \vec{m}\cdot\vec{m}} &= 1 + |\vec{m}|^2/2 - |\vec{m}|^4/8 + \cdots \\ \sqrt{1 + \vec{n}\cdot\vec{n}} &= 1 + |\vec{n}|^2/2 - |\vec{n}|^4/8 + \cdots \\ \vec{m}\cdot \vec{n} &= |\vec{m}|\,|\vec{n}|\cos \theta \end{align*} Putting these together, replacing $$\vec{n}$$ and $$\theta$$, and expanding in $$|\vec{\varepsilon}|$$, we obtain $$1 + |\vec{m}|^4/4 - |\vec{m}|^3|\vec{\varepsilon}|/2 - |\vec{\varepsilon}|^2 + |\vec{m}|^2|\vec{\varepsilon}|^2/4 + \cdots \text{.}$$ In your early examples, $$\vec{m} = \vec{0}$$ and $$|\vec{\varepsilon}| \approx 10^{-8}$$. Since the surviving term is $$|\vec{\varepsilon}|^2 \approx 10^{-16} \approx 2^{-53}$$, it should be unsurprising that you found nonzero inputs that when evaluated passed the representation of $$1$$ to acosh and obtained $$0$$.

To fix this, we need to have better precision resolution near $$1$$. An example of a library function that does this is log1p. We need an acosh1p so that we get the usual high precision resolution of a double near zero, but with semantics of a value near one. Your Answer finds one using log1p. It's the same one I found in the source for glibc's libm in e_acosh.c.

Your large aspect ratio example runs into a different/additional problem. We have $$\theta = \arctan(|\vec{\varepsilon}|/|\vec{m}|) = |\vec{\varepsilon}|/|\vec{m}| - \cdots \text{.}$$ In the example, $$|\vec{n}|\approx |\vec{m}| = 10^6$$ and $$|\vec{\varepsilon}| = 2 \cdot 10^{-2}$$, so $$\theta \approx 10^{-14} \approx 2^{-46}$$ so subtracting the dot product cancels the leading 46 bits of the mantissa, meaning that you are mostly passing guard bits (of unknown quality, possibly zeroes, possibly noise, possibly deterministic noise) to acosh and we should not be surprised that the large derivative of arcosh near $$1$$ amplifies the contribution of these guard bits to the output.

To address this, \begin{align*} \cos \theta &= \cos \arctan \frac{|\vec{\varepsilon}|}{|\vec{m}| } \\ &= 1 - \frac{|\vec{\varepsilon}|^2}{2|\vec{m}|^2} + \frac{3|\vec{\varepsilon}|^4}{8|\vec{m}|^4} + \cdots \end{align*} and expanding $$\vec{n}$$ around $$\vec{m}$$ (and multiplying the terms out), \begin{align*} &\sqrt{(1+\vec{m}\cdot\vec{m})(1+\vec{n}\cdot\vec{n})} - |\vec{m}|\,|\vec{n}| \left( 1 - \frac{|\vec{\varepsilon}|^2}{2|\vec{m}|^2} + \frac{3|\vec{\varepsilon}|^4}{8|\vec{m}|^4} + \cdots \right) + \cdots \\ &\quad {} = 1 + \frac{|\vec{m}|^2}{2(|\vec{m}|^2+1)} - \frac{|\vec{m}|\,|\vec{n}|}{|\vec{m}|^2+1} + \frac{|\vec{n}|^2}{2(|\vec{m}|^2+1)} + \frac{|\vec{n}|\,|\vec{\varepsilon}|^2}{2|\vec{m}|} - \frac{2|\vec{n}|\,|\vec{\varepsilon}|^4}{8|\vec{m}|^3} + \cdots \\ &\quad {} = 1 + \frac{|\vec{m}| - |\vec{n}|}{2(|\vec{m}|^2+1)} + \frac{|\vec{n}|\,|\vec{\varepsilon}|^2}{2|\vec{m}|} - \frac{2|\vec{n}|\,|\vec{\varepsilon}|^4}{8|\vec{m}|^3} + \cdots \end{align*}

In the large aspect ratio example, the second term is much smaller than the third, so the ordering of the terms above may not always be decreasing. Nevertheless, passing the second and third (and maybe fourth, on the off chance it's not too small) to your acosh1p should also address your large aspect ratio examples.

• I'm not sure what you mean by "prior to converting coordinates to floating point.” When dealing with numerical computations, typically all your data is floating point, there is no "before conversion." 2 days ago
• @D0SBoots : One computes for insight, not to randomly bang bits together. Consequently, problem instances come from somewhere that is not of limited precision. 2 days ago
• I feel like you took a tangent off into philosophy, instead of numerical analysis. Yes, the field of reals is not of limited precision, and we use the tools of continuous arithmetic to analyze functions we plan to compute. But this does not change the fact that these functions have a finite and discrete domain, since the inputs are floating-point (or other realizable) values. yesterday