Approximation of Riemannian Distances and Applications to Distance-Based Learning on Manifolds

TL;DR

This paper introduces a novel method to efficiently approximate Riemannian distances on manifolds, enabling scalable distance-based learning by solving a linear number of boundary value problems, demonstrated in shape analysis applications.

Contribution

The paper presents a new approximation technique for Riemannian distances that reduces computational cost, suitable for large datasets in machine learning and data analysis.

Findings

Significantly reduces computation time for Riemannian distances

Effective in shape analysis on landmark spaces

Applicable to large-scale distance-based learning

Abstract

Several important algorithms for machine learning and data analysis use pairwise distances as input. On Riemannian manifolds these distances may be prohibitively costly to compute, in particular for large datasets. To tackle this problem, we propose a distance approximation which requires only a linear number of geodesic boundary value problems to be solved. The approximation is constructed by fitting a two-dimensional model space with constant curvature to each pair of samples. We demonstrate the usefulness of our approach in the context of shape analysis on landmarks spaces.

Figures

Equations

$$\operatorname{dist}\big{(}\exp_{x}(u),\exp_{x}(v)\big{)}^{2}=\|u-v\|_{g_{x}}^{2}-\frac{1}{3}R_{x}(u,v,v,u)+O(\|u\|+\|v\|)^{6}.$$

$$g_x(u,u)={\tilde{g}}_{\tilde{x}}(\tilde{u},\tilde{u}),g_x(v,v)={\tilde{g}}_{\tilde{x}}(\tilde{v},\tilde{v}),g_x(u,v)={\tilde{g}}_{\tilde{x}}(\tilde{u},\tilde{v}),k_x(u,v)={\tilde{k}}_{\tilde{x}}(\tilde{u},\tilde{v}).$$

$$\operatorname{dist}\big{(}\exp_{x}(u),\exp_{x}(v)\big{)}^{2}=\widetilde{\operatorname{dist}}\big{(}\exp_{\tilde{x}}(\tilde{u}),\exp_{\tilde{x}}(\tilde{v})\big{)}^{2}+O(\|u\|+\|v\|)^{6}.$$

$$(u,v)$$

$$\displaystyle\operatorname{dist}\big{(}\exp_{x}(u),\exp_{x}(v)\big{)}$$

$$\displaystyle=\left\{\begin{aligned} &r\arccos\big{(}\cos\|u/r\|\cos\|v/r\|+\sin\|u/r\|\sin\|v/r\|\cos\phi\big{)},&&k>0,\\ &\|u-v\|,&&k=0,\\ &2r\operatorname{arctanh}\sqrt{1-\frac{2}{1+\cosh\|u/r\|\cosh\|v/r\|-\sinh\|u/r\|\sinh\|v/r\|\cos\phi}},&&k<0,\end{aligned}\right.$$

$$\displaystyle\lim_{t\to\infty}\frac{\operatorname{dist}\big{(}\exp_{x}(tu),\exp_{x}(tv)\big{)}}{t}=\left\{\begin{aligned} &0,&&k>0,\\ &\|u-v\|,&&k=0,\\ &\|u\|+\|v\|,&&k<0,\cos(\phi)\neq 1,\\ &\big{|}\|u\|-\|v\|\big{|},&&k<0,\cos(\phi)=1.\end{aligned}\right.$$

$$D(t)$$

$$\displaystyle=1+\cosh\big{(}\|tu/r\|-\|tv/r\|\big{)}+(1-\cos\phi)\sinh\|tu/r\|\sinh\|tv/r\|.$$

$$\displaystyle\lim_{t\to\infty}\frac{1}{t}\operatorname{dist}\big{(}\exp_{x}(tu),\exp_{x}(tv)\big{)}=\lim_{t\to\infty}\frac{d}{dt}2r\operatorname{arctanh}\sqrt{1-\frac{2}{D(t)}}=\lim_{t\to\infty}\frac{rD^{\prime}(t)}{D(t)},$$

$$r{D}^{\prime }(t)=(\Vert u\Vert -\Vert v\Vert \cos \varphi )\cosh \Vert tv/r\Vert \sinh \Vert tu/r\Vert +(\Vert v\Vert -\Vert u\Vert \cos \varphi )\cosh \Vert tu/r\Vert \sinh \Vert tv/r\Vert .$$

$$\underset{t\to \infty }{\lim }\frac{r{D}^{\prime }(t)}{D(t)}$$

$$K_q(p,p)=\sum _{i,j=1}^m(p^i{)}^{\top }k(q^i,q^j)p^j\in \mathbb{R},q\in {\mathrm{Land}}^m({\mathbb{R}}^d),p\in T_q^{\ast }{\mathrm{Land}}^m({\mathbb{R}}^d).$$

$$(\operatorname{Diff}(\mathbb{R}^{d}),k)\ni\phi\mapsto\big{(}\phi(\bar{q}^{1}),\dots,\phi(\bar{q}^{m})\big{)}\in(\operatorname{Land}^{m}(\mathbb{R}^{d}),K).$$

$$\forall x,y\in {\mathbb{R}}^d:k(x,y)=\exp \left(\frac{\Vert x-y{\Vert }^2}{2{\sigma }^2}\right)I^{d\times d},$$

Full text

We thank Peter W. Michor, Xavier Pennec, and Stefan Sommer for helpful discussions and gratefully acknowledge support by FRIAS and Penn State University via the Freiburg–Penn State Collaboration Development Program.

1. 1 Introduction

2. 1.1 Context

3. 1.2 Relation to previous work

4. 1.3 Contribution

5. 1.4 Outlook

6. 2 Approximations of pairwise distances

7. 2.1 Taylor approximation of squared distances

8. 2.2 Approximation by constant curvature spaces

9. 3 Applications to shape analysis

10. 3.1 Landmark spaces with kernel metrics

11. 3.2 Computation of sectional curvatures

12. 3.3 Numerical experiments

1. Introduction

1.1. Context

Several important algorithms for machine learning and data analysis on manifolds use pairwise distances as an input. For example, this is the case for multi-dimensional scaling and agglomerative clustering. The main computational burden is typically not the algorithm per se, but the calculation of the pairwise distance matrix. Indeed, each individual distance might be the solution of a costly optimization problem, as e.g. in the case of Riemannian shape analysis, and the size of the distance matrix grows quadratically in the number of samples. This can render distance-based learning prohibitively slow compared to alternative non-geometric methods.

1.2. Relation to previous work

This paper builds on the idea of [5] to use distance approximations for reducing the number of boundary value problems from quadratic down to linear in the number of samples. The approximations in [5] are based on the Baker–Campbell–Hausdorff formula and are therefore restricted to (quotients of) Lie groups with bi-invariant Riemannian metrics. In particular, they do not apply to landmark spaces with kernel metrics. Indeed, the numerical experiments in [5] do not show any improvement of second over first order approximations. This leads [6] to question if such improvements are possible at all. The present paper answers this question affirmatively.

1.3. Contribution

We propose a new second order approximation of the distance function. The approximation is obtained by fitting a two-dimensional model space with constant curvature to each pair of samples, represented by tangent vectors to the mean. The reduction to two dimensions is crucial because higher-dimensional model spaces with prescribed sectional curvature and closed-form expressions for geodesic distances are not available, to the best of our knowledge. As expected, the second-order approximation is more accurate than the first-order one in small-distance regimes, which cover typical applications in shape analysis. Moreover, in contrast to second-order Taylor polynomials, our approximate distances are always non-negative and have more realistic large-distance asymptotics. As an application, we present a numerical implementation of our distance approximation on landmark manifolds with kernel metrics and demonstrate its performance on some simple toy data.

1.4. Outlook

In the special case of landmark manifolds with kernel metrics, our sectional curvature computations could be sped up using Mario’s formula [2]. This will allow us to test our algorithm on high-dimensional data from real-world applications. Pseudo-landmark data will be of particular interest because the small distances of adjacent pseudo-landmarks lead to large sectional curvatures, which highlights the importance of the curvature correction in the second-order distance approximations.

2. Approximations of pairwise distances

This aim of this section is to construct approximations of the Riemannian distance function, which are fast to compute and accurate for sufficiently concentrated sample points. As an auxiliary tool we first describe Taylor expansions of geodesic distances. Subsequently, we develop an alternative and better approximation using model spaces with constant curvature.

2.1. Taylor approximation of squared distances

Recall that the squared Riemannian distance function on any Riemannian manifold is smooth away from the cut locus, including on a neighborhood of the diagonal. This allows one to use Taylor approximations of squared distances between sufficiently close points. More specifically, for any two points $y,z$ near $x$, one may calculate the Riemannian logarithms $u=\log_{x}(y)$ and $v=\log_{x}(z)$ and approximate the squared distance $\operatorname{dist}(y,z)^{2}$ by a Taylor polynomial of the function $(u,v)\mapsto\operatorname{dist}(\exp_{x}(u),\exp_{x}(v))^{2}$. This is made precise in the following lemma, which also provides an explicit expression of the Taylor series up to order five.

Lemma 1**.**

Let $(M,g)$ be a Riemannian manifold with exponential map $\exp$, distance function $\operatorname{dist}$, and curvature tensor $R$, let $x\in M$, and let $u,v\in T_{x}M$. Then it holds for sufficiently small $u$ and $v$ that

[TABLE]

Proof.

This expansion is well-known even for higher orders; see e.g. [3] or [4]. ∎

Note that the lowest-order term in the Taylor series ((1)) is the Euclidean distance in normal coordinates. The next term corrects for the influence of curvature and, by definition, improves the accuracy for small distances. However, for large distances, it may worsen the accuracy, can lead to negative signs, and has unrealistic quartic growth; cf. Figure 1. These problems are even worse for higher order Taylor expansions, which involve covariant derivatives of the curvature tensor. To address these problems, we propose an alternative approximation in the next section.

2.2. Approximation by constant curvature spaces

This section corrects some shortcomings of the Taylor approximation in the previous section. The idea is to replace the Taylor polynomial ((1)) by the squared distance function of a suitably selected constant curvature space. This guarantees non-negativity and leads to more realistic large-distance asymptotics while retaining the same order of accuracy. A visual comparison of the two approximations is presented in Figure 1.

In more details, the proposed approximation works as follows. For any points $y,z$ near $x\in M$, calculate the Riemannian logarithms $u=\log_{x}(y)$ and $v=\log_{x}(z)$ and sectional curvature $k_{x}(u,v)$, find a surface with constant curvature equal to $k_{x}(u,v)$, and approximate the squared distance between $y$ and $z$ by the squared distance function of the constant curvature surface. This is made precise in the following theorem.

Theorem 1**.**

Let $(M,g)$ be a Riemannian manifold with Riemannian distance function $\operatorname{dist}$ and sectional curvature $k$, let $x\in M$, let $u,v\in T_{x}M$, and let $(\tilde{M},\tilde{g})$ be a 2-dimensional constant curvature space with $\tilde{x}\in\tilde{M}$ and $\tilde{u},\tilde{v}\in T_{\tilde{x}}\tilde{M}$ such that

[TABLE]

Then it holds for sufficiently small $u$ and $v$ that

[TABLE]

Proof.

The conditions ensure that $R_{x}(u,v,v,u)=\tilde{R}_{\tilde{x}}(\tilde{u},\tilde{v},\tilde{v},\tilde{u})$. Therefore, the functions

[TABLE]

have the same Taylor expansion up to order 5 by Lemma 1. ∎

The approximation in Theorem 1 can be computed efficiently thanks to the following closed-form expression of the geodesic distance on constant curvature surfaces (i.e., the sphere, plane, and hyperbolic plane).

Lemma 2**.**

Let $(M,g)$ be a 2-dimensional manifold with constant sectional curvature $k\in\mathbb{R}$, let $r=\mathbbm{1}_{\{k\neq 0\}}|k|^{-1/2}$, let $x\in M$, let $\langle\cdot,\cdot\rangle$ and $\|\cdot\|$ denote to the Riemannian metric and norm on $T_{x}M$, respectively, let $u,v\in T_{x}M$, and let $\phi\in\mathbb{R}$ satisfy $\langle u,v\rangle=\|u\|\|v\|\cos\phi$. Then

[TABLE]

In contrast to the Taylor approximations in Lemma 1, the constant-curvature approximations in Theorem 1 behave well for large distances: they have asymptotically linear growth along geodesics emanating from the same point. This is made precise in the following lemma.

Lemma 3**.**

The distances in Lemma 2 have the following (sub)-linear growth rates:

[TABLE]

Proof.

As the statement is clear for spherical and Euclidean distances, as well as for points moving along the same geodesic, the only non-trivial case is $k<0$, $u\neq 0$, $v\neq 0$, and $\cos(\phi)\neq 1$. Let $D(t)$ be the denominator in the expression for the hyperbolic distance, i.e.,

[TABLE]

The latter formula shows that $D(t)$ tends to infinity for large $t$. Accordingly, the distance \operatorname{dist}\big{(}\exp_{x}(tu),\exp_{x}(tv)\big{)} tends to infinity for large $t$. This justifies the use of l’Hôpital’s rule, and one obtains

[TABLE]

provided that the limit on the right-hand side exists. One easily calculates

[TABLE]

Approximating all $\cosh$ and $\sinh$ by $\exp$, one obtains

[TABLE]

As an alternative to Theorem 1, it would be tempting to use $m$-dimensional instead of two-dimensional model spaces. This would lead to approximate distances of the same order of accuracy with the additional benefit that the triangle inequality holds. However, it is difficult to find model spaces with generically prescribed sectional curvatures at a point. Moreover, the geodesic distance on such model spaces, if they exist, might not have a closed-form expression. This points to the advantage of fitting model spaces separately for each pair of samples, as done in Theorem 1.

3. Applications to shape analysis

This section demonstrates the usefulness of distance approximations in the context of shape analysis. More specifically, we consider shape analysis on landmark spaces with kernel metrics. These metrics are widely used and provide an intuitive notion of similarity. As the calculation of the Riemannian distance on these spaces is computationally intensive, there is a high potential for significant speed-ups via distance approximations.

3.1. Landmark spaces with kernel metrics

For any $d,m\in\mathbb{N}_{>0}$, landmark space $\operatorname{Land}^{m}(\mathbb{R}^{d})$ is the set of all configurations of $m$ distinct points in $\mathbb{R}^{d}$. Landmark space is an open subset of $\mathbb{R}^{d\times m}$ and therefore a manifold. Any kernel $k\colon\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}^{d\times d}$ defines a Riemannian co-metric on landmark space via

[TABLE]

The corresponding Riemannian metric is denoted by $G$. If $k$ is the kernel of a reproducing Hilbert space of vector fields on $\mathbb{R}^{d}$, then $k$ can be interpreted as a right-invariant co-metric on the diffeomorphism group $\operatorname{Diff}(\mathbb{R}^{d})$ such that for any landmark configuration $\bar{q}\in\operatorname{Land}^{m}(\mathbb{R}^{d})$, the following map is a Riemannian submersion:

[TABLE]

In this sense the co-metric on landmark space is induced by a co-metric on the diffeomorphism group on ambient space. A typical choice of kernel is the Gaussian kernel, which is given by

[TABLE]

where $\sigma\in(0,\infty)$ and $I^{d\times d}$ denotes the $d\times d$ identity matrix.

3.2. Computation of sectional curvatures

The sectional curvatures on general Riemannian manifolds can be computed numerically by taking advantage of the automatic differentiation features of modern machine learning software, as demonstrated by [1]. This is also the approach we follow in our numerics. Alternatively, in the specific case of landmark manifolds with kernel metrics, the sectional curvatures could also be computed directly thanks to Mario’s formula, which provides closed-form expressions via the relation to the sectional curvatures on $\operatorname{Diff}(\mathbb{R}^{d})$, as described by [2]. If the kernel metric on $\operatorname{Diff}(\mathbb{R}^{d})$ was bi-invariant, which it is not, then the Riemannian exponential would coincide with the Lie group exponential, and the sectional curvatures could be computed using the Baker–Campbell–Hausdorff formula as suggested by [5].

The computational complexity of the curvature computations is as follows. Consider a dataset of $n$ landmark configurations $\{q_{1},\dots,q_{n}\}$ in $\operatorname{Land}^{m}(\mathbb{R}^{d})$, where the dimension $d$ is treated as a constant. Then the initial registration to some fixed template configuration $\bar{q}$, i.e., the computation of $v_{i}=\log_{\bar{q}}(q_{i})$ for each $i\in\{1,\dots,n\}$, has complexity $O(m^{2}n)$. Most importantly, only a linear number $n$ of boundary value problems has to be solved. Subsequently, the approximate pairwise distances can be computed at complexity $O(m^{4}n^{2})$ from the full Riemannian curvature tensor or at complexity $O(m^{2}n^{2})$ using Mario’s formula . While this quadratic-in-$n$ complexity is unavoidable when fully sampled distance matrices are required, the constants in the complexity bound can be very good. For example, if $m$ is sufficiently small to allow the computation of the full Riemannian curvature tensor at $\bar{q}$, then each sectional curvature can be computed by a few matrix-times-vector operations, and the corresponding approximate distances ((2)) by some additional trigonometric function evaluations.

3.3. Numerical experiments

We tested our distance approximation on some toy datasets of landmark configurations. More extensive tests on larger and higher-dimensional real-world datasets are intended in future work.

A first observation of our numerical experiments is that the second-order approximation via constant-curvature spaces outperforms the first-order one, as shown in Figure 2 for standard normally distributed landmark data. In this example the first-order approximations have a negative bias because the sectional curvatures are on average negative, but this bias is corrected by the second-order approximation. This confirms our theoretical predictions and stands in contrast to the results of [5], where an alternative second-order term did not lead to improved accuracy. Our findings are robust with respect to the variance of the noise and the width of the kernel: in all cases considered, the second-order approximation errors were smaller in mean and variance than the first-order ones.

A second observation is that even for moderately large datasets of approximately 15–20 shapes, the distance approximations are faster to compute than the true distances, despite the overhead of the curvature computations. This holds true despite the fact that our numerics do not yet rely on the explicit curvature formulas of [2], which provide potential for further improvement.

A third observation is that the second-order correction can be highly beneficial to distance-based learning. Indeed, 2(f) shows that multi-dimensional scaling based on true distances is nearly identical to multi-dimensional scaling based on second-order approximations, whereas the result based on first-order approximations is significantly off.