Hermite interpolation and data processing errors on Riemannian matrix manifolds

TL;DR

This paper develops a general Hermite interpolation framework on Riemannian manifolds and links data processing errors to manifold curvature, with applications to matrix factorizations like SVD and QR.

Contribution

It introduces a versatile interpolation method on Riemannian manifolds and relates data errors to sectional curvature, providing new error bounds for manifold data processing.

Findings

Interpolation framework applicable to matrix manifolds in data analysis and signal processing.

Error bounds depend on the sectional curvature of the manifold.

Numerical experiments demonstrate the method on orthogonal matrix factorizations.

Abstract

The main contribution of this paper is twofold: On the one hand, a general framework for performing Hermite interpolation on Riemannian manifolds is presented. The method is applicable, if algorithms for the associated Riemannian exponential and logarithm mappings are available. This includes many of the matrix manifolds that arise in practical Riemannian computing application such as data analysis and signal processing, computer vision and image processing, structured matrix optimization problems and model reduction. On the other hand, we expose a natural relation between data processing errors and the sectional curvature of the manifold in question. This provides general error bounds for manifold data processing methods that rely on Riemannian normal coordinates. Numerical experiments are conducted for the compact Stiefel manifold of rectangular column-orthogonal matrices. As use cases, we compute Hermite interpolation curves for orthogonal matrix factorizations such as the singular value decomposition and the QR-decomposition.

Equations

$$P(t_i)=f_i,\dot{P}(t_i)={\dot{f}}_i,i=0,\dots ,k.$$

$$f:[a,b]\to \mathcal{M},t\mapsto f(t)$$

$$c(t_i)=p_i\in \mathcal{M},\dot{c}(t_i)=v_{p_i}\in T_{p_i}\mathcal{M},i=0,\dots ,k.$$

$$O(r)=\{\Phi \in {\mathbb{R}}^{r\times r}|{\Phi }^T\Phi =\Phi {\Phi }^T=I_r\}.$$

$${\exp }_m(X):=\sum _{j=0}^{\infty }\frac{X^j}{j!},{\log }_m(I+X):=\sum _{j=1}^{\infty }(-1{)}^{j+1}\frac{X^j}{j}.$$

$${\mathrm{Exp}}_p^{\mathcal{M}}:T_p\mathcal{M}\supset D_0\to {\mathcal{D}}_p\subset \mathcal{M},v\mapsto {\mathrm{Exp}}_p^{\mathcal{M}}(v):=c_{p,v}(1)$$

$${\mathrm{Log}}_p^{\mathcal{M}}:\mathcal{M}\supset {\mathcal{D}}_p\to D_0\subset T_p\mathcal{M},q\mapsto v=({\mathrm{Exp}}_p^{\mathcal{M}}{)}^{-1}(q).$$

$$d{f}_p:T_p\mathcal{M}\to T_{f(p)}\mathcal{N}.$$

$$c(t)=a_0(t)p+a_1(t)q+b_0(t)v_0+b_1(t)v_1.$$

$$f:[t_0,t_1]\to \mathcal{M},t\mapsto f(t).$$

$$c:[t_0,t_1]\to \mathcal{M},c(t)={\mathrm{Exp}}_q^{\mathcal{M}}(\gamma (t)).$$

$$\gamma (t)=a_0(t){\Delta }_p+a_1(t){\Delta }_q+b_0(t){\hat{v}}_p+b_1(t){\hat{v}}_q.$$

$$\dot{c}(t_0)$$

$$\dot{c}(t_1)$$

$$\frac{d}{dt}\big{|}_{t=t_{1}}\operatorname{Exp}^{\mathcal{M}}_{q}(\gamma(t))=\frac{d}{dt}\big{|}_{h=0}\operatorname{Exp}^{\mathcal{M}}_{q}(h\hat{v}_{q}+\mathcal{O}(h^{2}))=d(\operatorname{Exp}^{\mathcal{M}}_{q})_{0}(\hat{v}_{q})=\hat{v}_{q}.$$

$$v_p=d({\mathrm{Exp}}_q^{\mathcal{M}}{)}_{{\Delta }_p}({\hat{v}}_p).$$

$$d({\mathrm{Log}}_q^{\mathcal{M}}{)}_p(v_p)$$

$$f(t_0)=p,f(t_1)=q\in \mathcal{M},\dot{f}(t_0)=v_p\in T_p\mathcal{M},\dot{f}(t_1)=v_q\in T_q\mathcal{M}$$

$${\hat{v}}_p=d({\mathrm{Log}}_q^{\mathcal{M}}{)}_p(v_p)\in T_q\mathcal{M},{\hat{v}}_q=v_q\in T_q\mathcal{M}.$$

$$c:[t_{0},t_{1}]\to\mathcal{M},\quad t\mapsto\operatorname{Exp}_{q}^{\mathcal{M}}\bigl{(}a_{0}(t)\Delta_{p}+b_{0}(t)\hat{v}_{p}+b_{1}(t)\hat{v}_{q}\bigr{)}$$

$$C:[t_0,t_k]\to \mathcal{M}t\mapsto c_{[t_i,t_{i+1}]}(t)\text{ for }t\in [t_i,t_{i+1}]$$

$$d({\mathrm{Log}}_q^{\mathcal{M}}{)}_p:$$

$$d({\mathrm{Exp}}_p^{\mathcal{M}}{)}_0:$$

$$d(\operatorname{Log}^{\mathcal{M}}_{q})_{p}(v_{p})=d(\operatorname{Log}^{\mathcal{M}}_{q})_{\tilde{\gamma}(0)}(\dot{\tilde{\gamma}}(0))=\frac{d}{ds}\big{|}_{s=0}\operatorname{Log}^{\mathcal{M}}_{q}\left(\tilde{\gamma}(s)\right).$$

$$\hat{v}_{p}=\frac{d}{ds}\big{|}_{s=0}\left(\operatorname{Log}^{\mathcal{M}}_{q}\circ\operatorname{Exp}^{\mathcal{M}}_{p}\right)(sv_{p}).$$

$${\hat{v}}_p=\frac{\left({\mathrm{Log}}_q^{\mathcal{M}}\circ {\mathrm{Exp}}_p^{\mathcal{M}}\right)(h{v}_p)-\left({\mathrm{Log}}_q^{\mathcal{M}}\circ {\mathrm{Exp}}_p^{\mathcal{M}}\right)(-h{v}_p)}{2h}+\mathcal{O}(h^2).$$

$${\mathrm{dist}}_{\mathcal{M}}(p,\tilde{p})\le |\delta -\tilde{\delta }|+s_0\delta (1-\frac{K_q(\sigma )}{6}{\delta }^2+o({\delta }^2))+\mathcal{O}(s_0^2),$$

$$w:=\frac{\Delta }{\Vert \Delta \Vert },w^{\perp }=\frac{\tilde{\Delta }-⟨w,\tilde{\Delta }⟩w}{\Vert \tilde{\Delta }-⟨w,\tilde{\Delta }⟩w\Vert }.$$

$$w:[0,\pi /2]\to T_q\mathcal{M},s\mapsto w(s)=\cos (s)w+\sin (s)w^{\perp }.$$

$${\mathrm{dist}}_{\mathcal{M}}(p,\tilde{p})$$

$$\le$$

$${\mathrm{dist}}_{\mathcal{M}}\left(p,{\mathrm{Exp}}_q^{\mathcal{M}}(\delta w(s_0))\right)$$

Full text

On the other hand, we expose a natural relation between data processing errors and the sectional curvature of the manifold in question. This provides general error bounds for manifold data processing methods that rely on Riemannian normal coordinates.

Numerical experiments are conducted for the compact Stiefel manifold of rectangular column-orthogonal matrices. As use cases, we compute Hermite interpolation curves for orthogonal matrix factorizations such as the singular value decomposition and the QR-decomposition.

1 Introduction

Given a data set that consists of locations $t_{0},\ldots,t_{k}\in\mathbb{R}$, function values $f_{0}=f(t_{0}),\ldots,f_{k}=f(t_{k})$ and derivatives $\dot{f}_{0}=\dot{f}(t_{0}),\ldots,\dot{f}_{k}=\dot{f}(t_{k})$, the (first-order) Hermite interpolation problem reads:

Find a polynomial $P$ of suitable degree such that

[TABLE]

Local cubic Hermite interpolation is the special case of Hermite-interpolating a two-points data set $\{f_{i},\dot{f}_{i},f_{i+1},\dot{f}_{i+1}\}$ on $t_{i},t_{i+1}\in\mathbb{R}$. Cubic Hermite interpolation is achieved by joining the local pieces on each sub-interval $[t_{i},t_{i+1}]$. By construction, the derivative at the end point of $[t_{i},t_{i+1}]$ coincides with the derivative of the start point of $[t_{i+1},t_{i+2}]$ so that the resulting curve is globally $C^{1}$, [20, Remark 7.7].

In this paper, we address the Hermite interpolation problem for a function that takes values on a Riemannian manifold $\mathcal{M}$ with tangent bundle $T\mathcal{M}$. More precisely, consider a differentiable function

[TABLE]

and a sample plan $a=t_{0},\ldots,t_{k}=b$. Sampling of the function values and the derivatives of $f$ at the parameter instants $t_{i}$ produces a data set consisting of manifold locations $p_{i}=f(t_{i})\in\mathcal{M}$ and velocity vectors $v_{p_{i}}\in T_{p_{i}}\mathcal{M}$ in the respective tangent spaces of $\mathcal{M}$ at $p_{i}$. The Hermite manifold interpolation problem is:

Find a curve $c:[a,b]\to\mathcal{M}$ of class $C^{1}$ such that

[TABLE]

1.1 Original contributions

(1) We introduce a method to tackle problem (2) that is a direct analogue to Hermite interpolation in Euclidean spaces. The method has the following features:

1. (i)

The approach works on arbitrary Riemannian manifolds, i.e., no special structure (Lie Group, homogeneous space, symmetric space,…) is required.

In order to conduct practical computations, only algorithms for evaluating the Riemannian exponential map and the Riemannian logarithm map must be available.111The Riemannian exp and log maps for some of the most prominent matrix manifolds are collected in [38]. 2. (ii)

The computational effort, in particular, the number of Riemannian exp and log evaluations is lower than that of any other Hermite manifold interpolation method known to the author.

(2) In addition, we expose a natural relation between data processing errors and the sectional curvature of the manifold in question. This provides general error bounds for data processing methods (including but not limited to interpolation) that work via a back-and-forth mapping of data between the manifold and its tangent space, or, more precisely, data processing methods that rely on Riemannian normal coordinates.

For convenience, the exposition will focus on cubic polynomial Hermite interpolation. However, the techniques may be readily combined with any interpolation method that is linear in the sampled locations and derivative values. Apart from polynomial interpolation, this includes radial basis function approaches [4] and gradient-enhanced Kriging [36].

As a use-case, we provide an explicit and efficient method for the cubic Hermite interpolation of column-orthogonal matrices, which form the so-called Stiefel manifold $St(n,r)=\{U\in\mathbb{R}^{n\times r}|U^{T}U=I\}$. Stiefel matrices arise in orthogonal matrix factorizations such as the singular value decomposition and the QR-decomposition.

1.2 Related work

Interpolation problems with manifold-valued sample data and spline-related approaches have triggered an extensive amount of research work.

It is well-known that cubic splines in Euclidean spaces are acceleration-minimizing. This property allows for a generalization to Riemannian manifolds in form of a variational problem for the intrinsic, covariant acceleration of curves, whose solutions can be interpreted as generalized cubic polynomials on Riemannian manifolds. The variational approach to interpolation on manifolds has been investigated e.g. in [28, 12, 11, 33, 9, 31, 22], see also [29] and references therein. While the property of minimal mean-acceleration is certainly desirable in many a context, including automobile, aircraft and ship designs and digital animations, there is no conceptual reason to impose this condition when interpolating general smooth non-linear manifold-valued functions.

A related line of research is the generalization of Bézier curves and the De Casteljau-algorithm [6] to Riemannian manifolds [29, 23, 27, 1, 16, 32]. Bézier curves in Euclidean spaces are polynomial splines that rely on a number of so-called control points. A Bézier curve starts at the first control point and ends at the last control point, the starting velocity is tangent to the line between the first two-pair of control points; the velocity at the endpoint is tangent to the line between the penultimate and the last control point. This is illustrated in Fig. 1. The number of control points determines the degree of the polynomial spline. To obtain the value $B(t)$ of a Bézier curve at time $t$, a recursive sequence of straight-line convex combinations of two locations must be computed. The transition of this technique to Riemannian manifolds is via replacing the inherent straight lines with geodesics [29]. The start and end velocities of the resulting spline are proportional to the velocity vectors of the geodesics that connect the first two and the last two control points, respectively [29, Theorem 1].

Note that the actual applications and use cases featured in the work referenced above are almost exclusively on low-dimensional matrix manifolds like $S^{2},S^{3},SO(3)$ or $SE(3)$.

A Hermite-type method that is specifically tailored for interpolation problems on the Grassmann manifold is sketched in [4, §3.7.4]. General Hermitian manifold interpolation has been considered explicitly in [21]. The idea is as follows: Given two points $p,q\in\mathcal{M}$ on a manifold and two tangent directions $v_{p}\in T_{p}\mathcal{M},v_{q}\in T_{q}\mathcal{M}$, the the authors of [21] approach the task to construct a connecting curve $c:[t_{i},t_{i+1}]\to\mathcal{M}$ such that $c(t_{i})=p,\dot{c}(t_{i})=v_{p},c(t_{i+1})=q,\dot{c}(t_{i+1})=v_{q}$ by constructing a “left” arc $l_{i}$ that starts at $t=t_{i}$ from $p$ with the prescribed velocity $v_{p}$ and a “right” arc $r_{i}$ that ends at $t=t_{i+1}$ at $q$ with the prescribed velocity $v_{q}$. The two arcs are then blended to a single spline arc via a certain geometric convex combination. In Euclidean spaces, this would read $s(t)=(1-\Phi(t))l_{i}(t)+\Phi(t)r_{i}(t)$, where $\Phi$ is a suitable weight function. Because a general Riemannian manifold lacks a vector space structure, the challenge is to construct a manifold analogue of a convex combination and [21] proposes a method that works on compact, connected Lie groups with a bi-invariant metric.

This same idea of blending a left and a right arc has been followed up in [16]. Here, the Euclidean convex combination is replaced with a geodesics average $s(t)=\operatorname{Exp}_{l_{i}(t)}(\Phi(t)\operatorname{Log}_{l_{i}(t)}(r_{i}(t)))$. In combination, this constitutes a valid approach for solving (2) in arbitrary Riemannian manifolds.222In practice, the building arcs $l_{i}(t)$ and $r_{i}(t)$ may be taken to be the geodesics with the prescribed velocities in their respective start and end points.

It should be mentioned that none of the papers on Bézier curves referenced above tackle the Hermite interpolation problem explicitly. However, the Bézier approach can be turned into an Hermite method by choosing the control points such that the sampled start and terminal velocities are met. It is clear that this requires at least $4$ control points in each subinterval $[t_{i},t_{i+1}]$, see Fig. 1.

Interpolation problems on Stiefel Manifolds have been considered in [23], however with using quasi-geodesics rather than geodesics. The work [39] includes preliminary numerical experiments for interpolating orthogonal frames on the Stiefel manifold that relies the canonical Riemannian Stiefel logarithm [30, 37].

Remark: (Hermite) interpolation of curves on Riemannian manifolds, i.e., of manifold-valued functions $f:[a,b]\rightarrow\mathcal{M}$ must not be confused with (Hermite) interpolation of real-valued functions with domain of definition on a manifold, $f:\mathcal{M}\rightarrow\mathbb{R}$. The latter line of research is pursued, e.g., in [26] but is not considered here.

1.3 Organization

The paper is organized as follows: Starting from the classical Euclidean case, Section 2 introduces an elementary approach to Hermite interpolation on general Riemannian manifolds. Section 3 relates the data processing errors of calculations in Riemannian normal coordinates to the curvature of the manifold in question. In Section 4, the specifics of performing Hermite interpolation of column-orthogonal matrices are discussed and Section 5 illustrates the theory by means of numerical examples. Conclusions are given in Section 6.

1.4 Notational specifics

The $(r\times r)$-identity matrix is denoted by $I_{r}\in\mathbb{R}^{r\times r}$, or simply $I$ if the dimensions are clear. The $(r\times r)$-orthogonal group is denoted by

[TABLE]

Throughout, the QR-decomposition $A=QR$ of $A\in\mathbb{R}^{n\times r}$, $n\geq r$, is understood as the ‘economy size’ QR-decomposition with $Q\in\mathbb{R}^{n\times r}$, $R\in\mathbb{R}^{r\times r}$.

The standard matrix exponential and the principal matrix logarithm are defined by

[TABLE]

The latter is well-defined for matrices that have no eigenvalues on $\mathbb{R}^{-}$.

For a Riemannian manifold $\mathcal{M}$, the geodesic that starts from $p\in\mathcal{M}$ with velocity $v\in T_{p}\mathcal{M}$ is denoted by $t\mapsto c_{p,v}(t)$. The Riemannian exponential function at $p$ is

[TABLE]

and maps a small star-shaped domain $D_{0}$ around the origin in $T_{p}\mathcal{M}$ diffeomorphically to a domain $\mathcal{D}_{p}\subset\mathcal{M}$, see Fig. 2. The Riemannian logarithm at $p$ is

[TABLE]

Recall that for a differentiable function $f:\mathcal{M}\to\mathcal{N}$, the differential at $p$ is a linear map between the tangent spaces

[TABLE]

2 Hermite interpolation on Riemannian manifolds

In this section, we construct a quasi-cubic spline between two data points $p_{0},p_{1}\in\mathcal{M}$ on a manifold with prescribed velocities $v_{p_{0}}\in T_{p_{0}}\mathcal{M}$ and $v_{p_{1}}\in T_{p_{1}}\mathcal{M}$. To this end, we develop a manifold equivalent to the classical local cubic Hermite interpolation in Euclidean spaces [20, §7].

2.1 The Euclidean case

We start with a short recap of Hermite cubic space curve interpolation, where the following setting is of special interest to our considerations. Let $\mathbb{V}$ be a real vector space and let $f:[t_{0},t_{1}]\to\mathbb{V}$ be differentiable with $f(t_{0})=p\in\mathbb{V}$, $f(t_{1})=q=0\in\mathbb{V}$, and derivative data $v_{0}=\dot{f}(t_{0}),v_{1}=\dot{f}(t_{1})$. When applied to vector-valued functions, the classical local cubic Hermite interpolating spline is the space curve $c(t)$ that is obtained via a linear combination of the sampled data,333It is an elementary, yet often overlooked fact that for functions $t\mapsto f(t)=(f_{1}(t),\ldots,f_{n}(t))^{T}\in\mathbb{R}^{n}$, component-wise polynomial interpolation of the coordinate functions $f_{l}(t)$ is equivalent to interpolating the coefficients in a linear combination of the sampled data vectors.

[TABLE]

For the reader’s convenience, the basic cubic Hermite polynomials coefficient polynomials $a_{0}(t),a_{1}(t),b_{0}(t),b_{1}(t)$ are listed in Appendix A.

2.2 Transfer to the manifold setting

Let $\mathcal{M}$ be a Riemannian manifold and consider a differentiable function

[TABLE]

Suppose that $f(t_{0})=p,f(t_{1})=q\in\mathcal{M}$ and $\dot{f}(t_{0})=v_{p}\in T_{p}\mathcal{M},\dot{f}(t_{1})=v_{q}\in T_{q}\mathcal{M}$ and assume further that $\operatorname{dist}(p,q)<r_{\mathcal{M}}(q)$, where $r_{\mathcal{M}}(q)$ is the injectivity radius of $\mathcal{M}$ at $q$. The latter condition ensures that the sample data lies within a domain, where the Riemannian normal coordinates are one-to-one, [13, p. 271].

Our approach is to express the interpolating curve in terms of normal coordinates centered at $q=f(t_{1})\in\mathcal{M}$,

[TABLE]

Hence, the task is transferred to constructing a curve $t\mapsto\gamma(t)\subset T_{q}\mathcal{M}$ such that the image curve $c$ under the exponential function solves the Hermite interpolation problem (2). Because $T_{q}\mathcal{M}$ is a vector space, we can utilize the ansatz of (6) but for $\mathbb{V}=T_{q}\mathcal{M}$,

[TABLE]

Here, $\Delta_{p}=\operatorname{Log}_{q}^{\mathcal{M}}(p),\Delta_{q}=\operatorname{Log}_{q}^{\mathcal{M}}(q)=0\in T_{q}\mathcal{M}$ are the normal coordinate images of the locations $p$ and $q$. The tangent vectors $\hat{v}_{p},\hat{v}_{q}\in T_{q}\mathcal{M}$ play the role of the velocity vectors and must be chosen such that

[TABLE]

Since the interpolating curve $c(t)=\operatorname{Exp}^{\mathcal{M}}_{q}(\gamma(t))$ is expressed in normal coordinates centered at $q=f(t_{1})$, condition (8) is readily fulfilled by selecting $\hat{v}_{q}=v_{q}$: According to the properties of the cubic Hermite coefficient functions $a_{0}(t),b_{0}(t),b_{1}(t)$, the Taylor expansion of $\gamma(t)$ around $t_{1}$ is $\gamma(t_{1}+h)=h\hat{v}_{q}+\mathcal{O}(h^{2})$. Therefore, up to first order, $\gamma(t)$ is a ray emerging from the origin $0\in T_{q}\mathcal{M}$ with direction $\hat{v}_{q}\in T_{q}\mathcal{M}$. Hence, the directional derivative of the exponential function is

[TABLE]

The latter equation holds, because $d(\operatorname{Exp}^{\mathcal{M}}_{q})_{0}=\operatorname{id}_{T_{q}\mathcal{M}}$, [13, §3, Prop. 2.9].

The condition (7) is more challenging, because the Taylor expansion of $\gamma(t)$ around $t_{0}$ is $\gamma(t_{0}+h)=\Delta_{p}+h\hat{v}_{p}+\mathcal{O}(h^{2})$ and is not a ray emerging from the origin $0\in T_{q}\mathcal{M}$. As the differential $d(\operatorname{Exp}^{\mathcal{M}}_{q})_{v}$ is not the identity for $v\neq 0$, the computation of \frac{d}{dt}\big{|}_{t=t_{1}}\operatorname{Exp}^{\mathcal{M}}_{q}(\Delta_{p}+h\hat{v}_{p}+\mathcal{O}(h^{2}))=d(\operatorname{Exp}^{\mathcal{M}}_{q})_{\Delta_{p}}(\hat{v}_{p}) is more involved. In fact, it is related to the Jacobi fields on a Riemannian manifold, see [13, §5], [24, §10] and the upcoming Section 3. Yet, for our purposes, it is sufficient to determine the tangent vector $\hat{v}_{p}$ such that

[TABLE]

As long as the sample points $p=\operatorname{Exp}^{\mathcal{M}}_{q}(\Delta_{p})$ and $q$ are not conjugate, we can make use of the fact that $\operatorname{Exp}^{\mathcal{M}}_{q}$ is a local diffeomorphism around $\Delta_{p}$, [24, Prop. 10.11]. Hence, under this assumption, (9) is equivalent to

[TABLE]

Recall that $v_{p}=\dot{f}(t_{0})$ is the given sample data. In summary, we have proved the following theorem.

Theorem 1**.**

Let $t_{0}<t_{1}\in\mathbb{R}$ and let $f:[t_{0},t_{1}]\rightarrow\mathcal{M}$ be a differentiable function on a Riemannian manifold $\mathcal{M}$. Suppose that

[TABLE]

and assume that $p$ and $q$ are not conjugate along the geodesic $t\mapsto\operatorname{Exp}^{\mathcal{M}}_{q}(t\operatorname{Log}^{\mathcal{M}}_{q}(p))$ that connects $p$ and $q$. Set $\Delta_{p}=\operatorname{Log}^{\mathcal{M}}_{q}(p)$ and

[TABLE]

Then

[TABLE]

with the cubic Hermite coefficient functions $a_{0}(t),b_{0}(t),b_{1}(t)$ as in (30)– (33) is a differentiable curve that solves the Hermite interpolation problem (2).

Remark 1**.**

Consider an Hermite sample set $p_{i}=f(t_{i})\in\mathcal{M}$, $v_{p_{i}}\in T_{p_{i}}\mathcal{M}$ $i=0,\ldots,k$. Then, by construction and in complete analogy to the Euclidean case, the composite curve

[TABLE]

that combines the local quasi-cubic spline arcs $c=c_{[t_{i},t_{i+1}]}$ of Theorem 1 is of class $C^{1}$ and solves the Hermite manifold interpolation problem (2).

Practical computation of $\hat{v}_{p}$

In cases, where an explicit formula for the Riemannian logarithm is at hand, the directional derivative $\hat{v}_{p}=d(\operatorname{Log}^{\mathcal{M}}_{q})_{p}(v_{p})$ can be directly computed. For general nonlinear manifolds $\mathcal{M}$, computing the differentials of the Riemannian exponential and logarithm is rather involved. According to (3), (4), (5), it holds

[TABLE]

with the usual identification a linear space with its tangent space.

In order to evaluate $d(\operatorname{Log}^{\mathcal{M}}_{q})_{p}(v_{p})$, we can take any differentiable curve $\tilde{\gamma}(s)\subset\mathcal{M}$ that satisfies $\tilde{\gamma}(0)=p$ and $\dot{\tilde{\gamma}}(0)=v_{p}$. Then,

[TABLE]

An obvious choice is $\tilde{\gamma}(s)=\operatorname{Exp}_{p}(sv_{p})\subset\mathcal{M}$. The final equation for computing $\hat{v}_{p}$ as required by (7) is

[TABLE]

The composite map $\operatorname{Log}^{\mathcal{M}}_{q}\circ\operatorname{Exp}^{\mathcal{M}}_{p}:T_{p}\mathcal{M}\supset D_{0}\to T_{q}\mathcal{M}$ is in fact a transition function for the normal coordinate charts. It is defined on an open subset of a Hilbert space and maps to a Hilbert space, see [25, Fig. 1.6, p. 12] for an illustration. Hence, we can approximate the directional derivative \hat{v}_{p}=\frac{d}{ds}\big{|}_{s=0}\left(\operatorname{Log}^{\mathcal{M}}_{q}\circ\operatorname{Exp}^{\mathcal{M}}_{p}\right)(sv_{p}) via finite difference approaches:

[TABLE]

2.3 Computational effort and preliminary comparison to other methods

Computationally, the most involved numerical operations are the evaluations of Riemannian $\operatorname{Log}$- and $\operatorname{Exp}$-mappings. Therefore, as in [16], we measure the computational effort associated with the Hermite interpolation method as the number of such function evaluations.

Constructing a quasi-cubic Hermite interpolant as in Remark 1 requires on each subinterval $[t_{i},t_{i+1}]$

• •

one Riemannian logarithm to compute $\Delta_{p}=\operatorname{Log}^{\mathcal{M}}_{q}(p)$,

• •

two Riemannian $\operatorname{Log}$- and $\operatorname{Exp}$-evaluations for the central difference approximation of (16),

which results in a total of $3k$ Riemannian $\operatorname{Log}$-evaluations and $2k$ Riemannian $\operatorname{Exp}$-evaluations for the whole composite curve. The data to represent the curve (11) can be precomputed and stored.

Evaluating a quasi-cubic Hermite interpolant at time $t$ requires a single Riemannian $\operatorname{Exp}$-evaluation.

As mentioned in the introduction, Bézier-like approaches may be used to tackle the the Hermite interpolation problem (2). This requires a cubic degree and at least four control points on each sub-interval $[t_{i},t_{i+1}]$ to impose the derivative constraints, see Fig. 1. The most efficient of such methods in [16] requires $\mathcal{O}(k^{2})$ Riemannian $\operatorname{Log}$-evaluations for constructing the curve data. Evaluating the curve at time $t$ requires $3$ Riemannian $\operatorname{Exp}$-evaluations plus $1$ Riemannian $\operatorname{Log}$-evaluation [16, Prop. 5.10].

3 Error propagation

The approach introduced in Section 2.2 follows the standard principle of (1) mapping the sampled data onto the tangent space, (2) performing data processing (in this case, interpolation) in the tangent space, (3) mapping the result back to the curved manifold. In this section, we perform a general qualitative analysis of the behavior of the actual errors on the manifold in question in relation to the data processing errors that accumulate in the tangent space. In particular, this allows to obtain error estimates for any manifold interpolation procedure based on the above standard principle and also applies to other data processing operations that subordinate to this pattern. In essence, the error propagation is related to the manifold’s curvature via a standard result from differential geometry on the spreading of geodesics [13, Chapter 5, §2].

Theorem 2**.**

Let $\mathcal{M}$ be a Riemannian manifold, let $q\in\mathcal{M}$ and consider tangent vectors $\Delta,\tilde{\Delta}\in T_{q}\mathcal{M}$, which are to be interpreted as exact datum and associated approximation. Write $\delta=\|\Delta\|$, $\tilde{\delta}=\|\tilde{\Delta}\|$, where it is understood that the norm is that of $T_{q}\mathcal{M}$. Assume that $\delta,\tilde{\delta}<1$. Let $\sigma=\text{span}(\Delta,\tilde{\Delta})\subset T_{q}\mathcal{M}$ and let $K_{q}(\sigma)$ be the sectional curvature at $q$ with respect to the $2$-plane $\sigma$.

If $s_{0}=\angle(\tilde{\Delta},\Delta)$ is the angle between $\tilde{\Delta}$ and $\Delta$, then the Riemannian distance between the manifold locations $p=\operatorname{Exp}_{q}^{\mathcal{M}}(\Delta)$ and $\tilde{p}=\operatorname{Exp}_{q}^{\mathcal{M}}(\tilde{\Delta})$ is

[TABLE]

with the underlying assumption that all data is within the injectivity radius at $q$.

Proof.

Formally, it holds $\operatorname{dist}_{\mathcal{M}}(p,\tilde{p})=\|\operatorname{Log}_{p}^{\mathcal{M}}(\tilde{p})\|_{T_{p}\mathcal{M}}$. However, the data is given in normal coordinates around $q\in\mathcal{M}$ and not around $p$ (nor $\tilde{p}$). The Riemannian exponential is a radial isometry (lengths of rays starting from the origin of the tangent space equal the lengths of the corresponding geodesics). Yet, it is not an isometry so that $\operatorname{dist}_{\mathcal{M}}(p,\tilde{p})\neq\|E\|$, unless $\mathcal{M}$ is flat. Therefore, we will estimate the distance $\operatorname{dist}_{\mathcal{M}}(p,\tilde{p})$ against a component that corresponds to the length of a ray in $T_{q}\mathcal{M}$ and a circular segment in $T_{q}\mathcal{M}$.

To this end, introduce an orthonormal basis ${w,w^{\bot}}$ for the plane $\sigma=\text{span}(\Delta,\tilde{\Delta})\subset T_{q}\mathcal{M}$ via

[TABLE]

The circular segment $w(s)$ of the unit circle in the $\sigma$-plane that starts from $w=w(0)$ and ends in $w^{\bot}=w(\pi/2)$ can be parameterized via the curve

[TABLE]

Let $s_{0}\in[0,\pi/2]$ be the angle such that $\tilde{\Delta}=\tilde{\delta}w(s_{0})$. This setup is illustrated in Figure 3, where the outer dashed circular arc indicates the unit circle and the solid circular arcs are the circles of radius $\delta$ and $\tilde{\delta}$, respectively. By the triangle inequality,

[TABLE]

Since the points $\delta w(s_{0})$ and $\tilde{\delta}w(s_{0})=\tilde{\Delta}$ are on a ray that emerges from the origin in $T_{q}\mathcal{M}$, the distance term in line (20) is exactly $|\delta-\tilde{\delta}|$, see Figure 3. Note that $\operatorname{Exp}_{q}^{\mathcal{M}}(\delta w(0))=\operatorname{Exp}_{q}^{\mathcal{M}}(\Delta)=p$. Hence, the distance term in line (20) is

[TABLE]

A Taylor expansion centered at $s=0$ of the transition function along the circular segment $s\mapsto\left(\operatorname{Log}_{p}^{\mathcal{M}}\circ\operatorname{Exp}_{q}^{\mathcal{M}}\right)(\delta w(s)))$ gives

[TABLE]

To arrive at the last line, $d(\operatorname{Log}_{p}^{\mathcal{M}})_{p}=\operatorname{id}_{T_{p}\mathcal{M}}$ was used, which follows from the standard result $d(\operatorname{Exp}_{p}^{\mathcal{M}})_{0}=\operatorname{id}_{T_{p}\mathcal{M}}$ [13, §3, Prop. 2.9, p. 65] with the usual identification of $T_{p}\mathcal{M}\cong T_{p}(T_{p}\mathcal{M})$, cf. (12), (13).

By construction, $\delta\mapsto d(\operatorname{Exp}_{q}^{\mathcal{M}})_{\delta w}(\delta w^{\bot})$ is a Jacobi field along the geodesic ray that starts from $q=\operatorname{Exp}_{q}^{\mathcal{M}}(0)$ with unit velocity $w\in T_{q}\mathcal{M}$, see [13, §5, Prop. 2.7, p. 114]. Moreover, $w=(0),w^{\bot}=w(\pi/2)$ constitute an orthonormal basis of the $\sigma$-plane in $T_{q}\mathcal{M}$. Therefore, the results [13, §5, Cor. 2.9, Cor. 2.10, p. 115] apply and give

[TABLE]

In summary,

[TABLE]

which establishes the theorem. ∎

Remark 2**.**

1. (i)

The the approximation error in the tangent space $\epsilon:=\|E\|=\|\tilde{\Delta}-\Delta\|$ can be related by elementary trigonometry to the angle $s_{0}=\angle(\tilde{\Delta},\Delta)$. It holds $\epsilon\geq\delta\|w(0)-w(s_{0})\|$, see Fig. 3. Moreover, $s_{0}=2\arcsin\left(\frac{\|w(0)-w(s_{0})\|}{2}\right)$. Thus,

[TABLE]

In regards of practical applications, it is safe to assume $\epsilon<\delta$. Then, in terms of the error $\epsilon$, the distance estimate (17) reads

[TABLE] 2. (ii)

If we travel from $\Delta$ to $\tilde{\Delta}$ in the tangent space on the corresponding curves as in the proof of Theorem 2, i.e. first along the circular arc from $\Delta=\delta w(0)$ to $\delta w(s_{0})$ and than along the ray from $\delta w(s_{0})$to $\tilde{\delta}w(s_{0})$, then we cover precisely a distance of $|\delta-\tilde{\delta}|+\delta s_{0}$. Comparing this with (17), we see that the corresponding distances of the manifold images are (asymptotically) \left\{\begin{array}[]{l}\text{shorter,}\\ \text{longer,}\end{array}\right\} if $\mathcal{M}$ features \left\{\begin{array}[]{l}\text{positive}\\ \text{negative}\\ \end{array}\right\} sectional curvatures. The underlying principle is the well-known effect that geodesics on positively curved spaces spread apart less than straight rays in a flat space, while they spread apart more on negatively curved spaces, see **[13, §5, Remark 2.11, p. 115/116]**.

From the numerical point of view, this means that data processing operations that work in the tangent space followed by a transition to the manifold are rather well-behaved on manifolds of positive curvature, while the opposite holds on negatively curved manifolds. In Section 4, we will show an illustration of Theorem 2 on an interpolation problem on the compact Stiefel manifold.

With the help of Theorem 2, explicit error bounds for manifold interpolation methods can be obtained. For example, cubic Hermite interpolation comes with a standard error bound [20, Thm 7.16] that applies to the interpolant in the tangent space. This can be forwarded to a manifold error via Theorem 2.

4 Cubic Hermite interpolation of column-orthogonal matrices

The set of column-orthogonal matrices

[TABLE]

is the compact homogeneous matrix manifold known as the (compact) Stiefel manifold. This section reviews the essential aspects of the numerical treatment of Stiefel manifolds. For more details, see [2, 14, 38].

The tangent space $T_{U}St(n,r)$ at a point $U\in St(n,r)$ can be thought of as the space of velocity vectors of differentiable curves on $St(n,r)$ passing through $U$:

[TABLE]

For any matrix representative $U\in St(n,r)$, the tangent space of $St(n,r)$ at $U$ is

[TABLE]

Every tangent vector $\Delta\in T_{U}St(n,r)$ may be written as

[TABLE]

The dimension of both $T_{U}St(n,r)$ and $St(n,r)$ is $nr-\frac{1}{2}r(r+1)$.

Each tangent space carries an inner product $\langle\Delta,\tilde{\Delta}\rangle_{U}=tr\left(\Delta^{T}(I-\frac{1}{2}UU^{T})\tilde{\Delta}\right)$ with corresponding norm $\|\Delta\|_{U}=\sqrt{\langle\Delta,\Delta\rangle_{U}}$. This is called the canonical metric on $T_{U}St(n,r)$. It is derived from the quotient space representation $St(n,r)=O(n)/O(n-r)$ that identifies two square orthogonal matrices in $O(n)$ as the same point on $St(n,r)$, if their first $r$ columns coincide [14, §2.4]. For a condensed introduction to quotient spaces, see [38, §2.5]. Endowing each tangent space with this metric (that varies differentiably in $U$) turns $St(n,r)$ into a Riemannian manifold. The associated sectional curvature is non-negative and is bounded by $0\leq K_{U}(\sigma)\leq\frac{5}{4}$ for all $U\in St(n,r)$ and all two-plans $\sigma=\text{span}(\Delta,\tilde{\Delta})\subset T_{U}St(n,r)$, [30, §5].

Given a start point $U\in St(n,r)$ and an initial velocity $\Delta\in T_{U}St(n,r)$ the Stiefel geodesic $c_{U,\Delta}$ (and thus the Riemannian exponential) is

[TABLE]

where

[TABLE]

is the decomposition of the tangent velocity into its horizontal and vertical component with respect to the base point $U$, [14]. Because $\Delta$ is tangent, $A=U^{T}\Delta\in\mathbb{R}^{r\times r}$ is skew. The Riemannian Stiefel logarithm can be computed with the algorithm of [37].

4.1 Differentiating the Stiefel exponential

In this section, we compute the directional derivative of the Stiefel exponential

[TABLE]

This is important for two reasons.

Differentiable gluing of interpolation curves. Consider a manifold data set $t_{i},p_{i}=f(t_{i})$, $i=0,\ldots,j,j+1,\ldots,k$, where the Riemannian distance, say, of the sample points $p_{j}$ and $p_{0}$ and $p_{j}$ and $p_{k}$ exceeds the injectivity radius of $\mathcal{M}$ at $p_{j}$. Then, simple tangent space interpolation with mapping the data set to $T_{p_{j}}\mathcal{M}$ is not possible. A remedy is to split the data set at $p_{j}$ and to compute two interpolation curves, one for the sample set $t_{i},p_{i}=f(t_{i})$, $i=0,\ldots,j$ and one for the sample set $t_{i},p_{i}=f(t_{i})$, $i=j,j+1,\ldots,k$. With the canonical method of tangent space interpolation, the curves have the expressions $c_{1}(t)=\operatorname{Exp}_{p_{\lfloor j/2\rfloor}}(\sum_{i=0}^{j}a_{i}(t)\Delta_{i})$ and $c_{2}(t)=\operatorname{Exp}_{p_{j+\lfloor j/2\rfloor}}(\sum_{i=j}^{k}a_{i}(t)\Delta_{i})$, where $\Delta_{i}=\operatorname{Log}^{St}_{p_{\lfloor j/2\rfloor}}(p_{i})$ for $c_{1}$ and $\Delta_{i}=\operatorname{Log}^{\mathcal{M}}_{p_{j+\lfloor j/2\rfloor}}(p_{i})$ for $c_{2}$. Concatenating the curves $c_{1}$, $c_{2}$ will result in a non-differentiable kink at the intersection location $p_{j}$, where $c_{1}$ ends and $c_{2}$ starts. In order to avoid this, one can compute the derivative $\dot{c}_{1}(t_{j})$ and use $\dot{c}_{1}(t_{j})=\dot{c}_{2}(t_{j})$ as an Hermitian derivative sample when constructing $c_{2}$. For obtaining $\dot{c}_{1}(t_{j})$, a derivative of the form of (24) must be computed. 2. 2.

Method validation. The cubic Hermite manifold interpolation method of Theorem 1 requires the computation of $\hat{v}_{p}=d(\operatorname{Log}^{\mathcal{M}}_{q})_{p}(v_{p})$. As was mentioned in Section 2.2, the differential of the $\operatorname{Log}$-mapping cannot be computed explicitly for general manifolds $\mathcal{M}$. In order to assess the numerical quality of a finite-differences approximation, we can first compute $\hat{v}_{p}$ by (16) and then recompute (9)

[TABLE]

The numerical accuracy is assessed via the error

[TABLE]

Again, a derivative of the form of (24) must be computed.

Now, let us address the derivative (24) for $\Delta_{0},V\in T_{U}St(n,r)$. The underlying computational obstacle is that the exponential law does not hold for the matrix exponential and two non-commuting matrices $\exp_{m}(X+tY)\neq\exp_{m}(X)\exp_{m}(tY)$. Write $\Delta(t)=\Delta_{0}+tV$ and let $Q(t)R(t)=(I-UU^{T})\Delta(t)$ be the $t$-dependent QR-decomposition of the tangent space curve. Moreover, $A(t):=U^{T}\Delta(t)$ and $\dot{A}(0)=U^{T}V$. Then, by the product rule,

[TABLE]

Introduce the matrix function $M(t)=\begin{pmatrix}A(t)&-R^{T}(t)\\ R(t)&0\end{pmatrix}$. It is sufficient to compute d(\exp_{m})_{M(0)}(\dot{M}(0))=\frac{d}{dt}\big{|}_{t=0}\exp_{m}(M(0)+t\dot{M}(0)).555This is a common problem in Lie group theory, see [17, §5.4]. The solution is formally an infinite sequence of nested commutator products in $[M,\dot{M}]=M\dot{M}-\dot{M}M$,

\frac{d}{dt}\big{|}_{t=0}\exp_{m}(M+t\dot{M})=\exp_{m}(M)\left(\dot{M}-\frac{1}{2!}[M,\dot{M}]+\frac{1}{3!}[M,[M,\dot{M}]]-\cdots\right)

In the following, we often omit the parameter $t$ with the implicit understanding that all quantities are evaluated at $t=0$. By Mathias’ Theorem [19, Thm 3.6, p. 58], it holds

[TABLE]

Hence, for data stemming from $St(n,r)$, a $(4r\times 4r)$-matrix exponential must be computed. However, the advantage is that $\exp_{m}(M)$ and \frac{d}{dt}\big{|}_{t=0}\exp_{m}(M+t\dot{M}) are obtained in one go and both are needed for evaluating (24). Moreover, usually $r\ll n$ in practical applications. For details and alternative algorithms for computing \frac{d}{dt}\big{|}_{t=0}\exp_{m}(M+t\dot{M}), see [19, §10.6]. In summary:

Lemma 3**.**

With all quantities as introduced above, let

[TABLE]

be written in terms of subblocks of size $r\times r$. Then

[TABLE]

The derivatives of the QR-factors of the decomposition $Q(t)R(t)=(I-UU^{T})\Delta(t)$ that are required to compute $\dot{Q}$ and $\dot{M}=\begin{pmatrix}\dot{A}&-\dot{R}^{T}\\ \dot{R}&0\end{pmatrix}$ can be obtained from Alg. 1.

4.2 Alternative options for Hermite data preprocessing on St(n,r)

As outlined in Section 2, the Hermite interpolation problem with local Stiefel sample data $f(t_{0})=U,\dot{f}(t_{0})=\Delta$, $f(t_{1})=\tilde{U},\dot{f}(t_{1})=\tilde{\Delta}$, requires us to translate the derivative samples to a common tangent space. On $St(n,r)$, this amounts to compute

[TABLE]

for some differentiable curve $\tilde{\gamma}(s)\subset St(n,r)$ that satisfies $\tilde{\gamma}(0)=U$ and $\dot{\tilde{\gamma}}(0)=\Delta$. There are other option than $\tilde{\gamma}(s)=\operatorname{Exp}_{U}^{St}(s\Delta)$ with this property, which might be cheaper to evaluate, depending on the context: For a skew-symmetric $M_{0}$, the Cayley transformation [35, eq. (7)], [7, p. 284],

[TABLE]

produces a curve of orthogonal matrices that matches the matrix exponential $s\mapsto\exp_{m}(sM_{0})$ up terms of order $\mathcal{O}(s^{2})$. As a consequence,

[TABLE]

and the matrix curve $s\mapsto(U,Q)M(s)\begin{pmatrix}I_{r}\\ 0\end{pmatrix}$ with $M(s)$ from (28) based on $M_{0}=\begin{pmatrix}A&-R^{T}\\ R&0\end{pmatrix}$ as in (23) may be used as the curve $\tilde{\gamma}(s)$ in (14), (15) instead of the Stiefel exponential (23).

Another option is to use retractions as a replacement for the Riemannian exponential [2, §4.1]. By definition, the differential of a retraction map at the origin of the tangent space is the identity map and thus coincides with the differential of the Riemannian exponential at the origin. Suitable matrix curves $s\mapsto\tilde{\gamma}(s)$ that match $\operatorname{Exp}_{U}^{St}(s\Delta)$ up to terms of first order based on Stiefel retractions are

[TABLE]

see [2, Example 4.1.3].

A word of caution: With the QR-based retraction, there is the challenge of computing a differentiable QR-path. Numerical QR-algorithms in high-level programming environments like MATLAB or SciPy might provide discontinuous matrix paths, e.g., because of different internal pivoting strategies.

5 Examples and experimental results

In this section, we conduct various numerical experiments that put the theoretical findings in perspective. All examples are coded and performed in the SciPy programming environment [15].

5.1 The numerical accuracy of the derivative translates

Before we start with the actual interpolation problems, we assess the numerical accuracy of the process of mapping a velocity sample $v_{U}\in T_{U}St(n,r)$ to a tangent velocity $\hat{v}_{U}\in T_{\tilde{U}}St(n,r)$ for two different Stiefel locations $U,\tilde{U}\in St(n,r)$. This requires the numerical computation of $\hat{v}_{U}=d(\operatorname{Log}^{St}_{\tilde{U}})_{U}(v_{U})$ with the help of central finite differences as in (16).

Then, we reverse this process to recover the original input via

[TABLE]

To this end, we utilize formula (27) of Lemma 3. Then we compute the error $\mathcal{E}=\frac{\|v_{U,rec}-v_{U}\|_{F}}{\|v_{U}\|_{F}}.$ As data points, we use samples of the Stiefel function $\mu\mapsto U(\mu)\in St(n,r)$, $n=1001,r=6$ that is featured in the upcoming Section 5.4: More precisely, $U=U(0.9)$, $\tilde{U}=U(1.4)$. The tangent direction to be translated is chosen as $v_{U}=Log_{U}^{St}(U(1.9))\in T_{U}St(n,r)$. All Riemannian log computations are performed with the Algorithm of [37] and a numerical convergence threshold of $\tau=10^{-14}$. The next table shows the reconstruction error $\mathcal{E}$ versus the finite difference step size $h$ used in (16).

[TABLE]

.

Even though there are various numerical processes involved (matrix exp, matrix log, numerical QR-differentiation, iterative Stiefel logarithm etc.) the accuracy of the finite difference approach is surprisingly high. In the following experiments, a step size of $h=10^{-4}$ is used to calculate (16).

5.2 Hermite interpolation of the Q-factor of a QR-decomposition

As a first example, consider a cubic matrix polynomial

[TABLE]

The matrices $Y_{i}$ were produced as random matrices with entries uniformly sampled from $[0,1]$ for $Y_{0}$, entries uniformly sampled from $[0,0.5]$ for $Y_{1},Y_{2}$ and from $[0,0.2]$ for $Y_{3}$. The $t$-dependent $QR$-decomposition is

[TABLE]

The matrix curve $Y(t)$ is sampled at $6$ Chebychev roots in $[-1.1,1.1]$.666 The Chebychev locations read $t_{0}\approx-1.0625$, $t_{1}\approx-0.7778$, $t_{2}\approx-0.2847$, $t_{3}\approx 0.2847$, $t_{4}\approx 0.7778$, $t_{5}\approx 1.0625$. At each sample point $t_{i}$ the Q-factor $Q(t_{i})$ of the QR-decomposition computed. The corresponding derivative $\dot{Q}(t_{i})$ is obtained from Alg. 1 in Appendix B. This constitutes the Hermite sample data set

[TABLE]

For comparison, the following interpolation schemes are conducted.

• •

Quasi-linear interpolation: In this case, the Stiefel samples are connected by geodesics as described in [38, §3.1]. No derivative information is used. This is the manifold version of linear interpolation.

• •

Tangent space interpolation: In this case, all data is mapped to single tangent space attached at $Q(t_{j})$, $j=\lfloor k/2\rfloor$, where $k$ is the number of sample points. Then, RBF interpolation is performed on tangent vectors as described in [5], [38, §3.1]. As an RBF, the inverse multiquadric is selected. No derivative information is used.

• •

Hermite quasi-cubic interpolation, as introduced in Section 2.2.

Since the quasi-linear and the quasi-cubic approach rely on piece-wise splines, it is only the ‘global’ tangent space interpolation that benefits from the choice of Chebychev samples.

For $t\in[-1.1,1.1]$, the relative interpolation errors are computed in the matrix Frobenius norm as $\frac{\|Q^{*}(t)-Q(t)\|_{F}}{\|Q(t)\|_{F}}$, where $Q^{*}(t)$ denotes the manifold interpolant and $Q(t)$ is the reference solution. The error curves are displayed in Fig. 4.

The relative Frobenius errors are

[TABLE]

.

The (discrete) $L_{2}$-norm gives the integrated squared errors on a discrete representation of the interval $[t_{0},t_{k}]$ with a resolution of $100$ points.

5.3 Hermite interpolation of a low-rank SVD

Next, we consider an academic example of a non-linear matrix function with fixed low rank. The goal is to perform a quasi-cubic interpolation of the associated SVD. As above, we construct a cubic matrix polynomial

[TABLE]

with random matrices $Y_{i}$ with entries uniformly sampled from $[0,1]$ for $Y_{0}$ and from $[0,0.5]$ for $Y_{1},Y_{2},Y_{3}$. Then, a second matrix polynomial is considered

[TABLE]

Here, the entries of $Z_{0}$ are sampled uniformly from $[0,1]$ while the entries of $Z_{1},Z_{2}$ are sampled uniformly from $[0,0.5]$. The nonlinear low-rank matrix function is set as

[TABLE]

By construction, $W(t)$ is of fixed $\text{rank}(W(t))\equiv r=10$ $\forall t$. The low rank SVD

[TABLE]

is sampled at the two Chebychev nodes $t_{0}=0.0732,t_{1}=0.4268$ in the interval $[0.0,0.5]$.777The Chebychev sampling is implemented as a standard in the program code that was written for the numerical experiments and is not of importance in this case. The Hermite sample data set

[TABLE]

is computed with Alg. 3 of Appendix C. (To this end $V(t_{i})=(V_{r}(t_{i}),V_{m-r}(t_{i}))\in O(m)$ is required.)

Remark 3**.**

Computing an analytic path of an SVD and thus a proper sample data set is a challenge in its own right, see [10]. This is in part because of the inherent ambiguity of the SVD even in the case of mutually distinct singular values, where $W=U\Sigma V^{T}=(US)\Sigma(SV^{T})$ for any orthogonal and diagonal matrix $S=\mathop{\operator@font diag}\nolimits(\pm 1,\ldots,\pm 1)$, [19, B.11, p. 334]. SVD algorithms from numerical linear algebra packages may return a different ‘sign-matrix’ $S$ for the SVD of $W(t)$ and $W(s)$, even when $t$ and $s$ are close to each other. This introduces discontinuities in the sampled $U$ and $V$ matrices. In the experiments performed in this work, we normalize the SVD as follows. A reference SVD $U_{0}\Sigma_{0}V_{0}^{T}=W(t_{0})$ is computed. At each $t$, we compute an SVD $U_{t}\Sigma_{t}V_{t}^{T}$ and determine $S=\text{sign}(\mathop{\operator@font diag}\nolimits(U_{t}^{T}U_{0}))$, where the sign-function is understood to be applied entry-wise on the diagonal elements. Then, we replace $U_{t}\leftarrow U_{t}S$, $V_{t}\leftarrow V_{t}S$. In the test cases considered here, this hands-on approach is sufficient to ensure a differentiable SVD computation. In general, one has to allow for negative singular values to ensure differentiability, [10].

For $t$ in the sampled range, the relative interpolation errors are computed in the Frobenius norm as $\frac{\|U^{*}(t)\Sigma^{*}(t)(V^{*}(t))^{T}-W(t)\|_{F}}{\|W(t)\|_{F}}$, where $U^{*}(t)$, $\Sigma^{*}(t)$, $V^{*}(t)$ are the interpolants of the matrix factors of the low-rank SVD of $W(t)$ and $W(t)$ is the reference solution. The relative errors are

[TABLE]

.

Fig. 5 displays the error curves for the quasi-linear and the quasi-cubic Hermite interpolation approaches.

For the sake of completeness, we repeat this experiment but with selecting $p=U(\mu_{i})$ as the center for the Riemannian normal coordinates. Hence, the derivative data is mapped to $T_{U(\mu_{i})}St(n,r)$ instead of $T_{U(\mu_{i+1})}St(n,r)$ and the tangent space interpolation curve is of the form

[TABLE]

where $p=U(\mu_{i}),q=U(\mu_{i+1})$. This leads to virtually indistinguishable plots. The maximum relative errors are $6.3023\cdot 10^{-4}$ ($q$-centered) vs. $6.3374\cdot 10^{-4}$ ($p$-centered). The $L_{2}$-norms of the relative errors are $2.3819\cdot 10^{-4}$ ($q$-centered) vs. $2.3954\cdot 10^{-4}$ ($p$-centered).

Recall that the local cubic Hermite interpolation scheme works in essence by performing Hermite interpolation in a selected tangent space and subsequently mapping the result to the manifold. For the $U$-factor interpolation featured in the above example, Fig. 6 shows the absolute interpolation errors of the tangent space data in the canonical Riemannian metric together with interpolation errors of final manifold data in terms of the Riemannian distance. The manifold errors are very close to the tangent space errors but are actually slightly smaller, inspite of the additional downstream translation of the tangent space interpolants to the manifold via the Riemannian exponential, which is an additional source of numerical errors. This is in line with Theorem 2, since the Stiefel manifold features positive sectional curvature.

5.4 Hermite interpolation of the left singular values of non-linear function snapshots

In the next experiment, we consider the SVD of discrete snapshots of a nonlinear multi-parameter function. To this end, define

[TABLE]

where $\langle f_{1},f_{2}\rangle_{L_{2}}=\int_{0}^{1}{f_{1}(x)f_{2}(x)dx}$ and $\|\cdot\|_{L_{2}}=\sqrt{\langle\cdot,\cdot\rangle_{L_{2}}}$ on $L_{2}([0,1])$. We will discretize $F$ in $x$, take function ‘snapshots’ at selected values of $t$ and eventually Hermite interpolate the left singular vectors of the discrete snapshot matrices with respect to $\mu$. The partial derivative of $F$ by $\mu$ is

[TABLE]

For the spatial discretization, we use an equidistant decomposition of the unit interval, $0=x_{1},x_{2},\ldots,x_{n}=1$, $n=1001$. Then, we take $r=6$ function snapshots in $t$ at time instants $t=1.0,\hskip 2.84544pt1.6,\hskip 2.84544pt2.2,\hskip 2.84544pt2.8,\hskip 2.84544pt3.4,\hskip 2.84544pt4.0$. In this way, a $\mu$-dependent snapshot matrix function $Y(\mu):=(F(x,t_{1},\mu),\ldots,F(x,t_{6},\mu))\in\mathbb{R}^{n\times r}=\mathbb{R}^{1001\times 6}$ with SVD

[TABLE]

is obtained.

Fig. 7 displays the snapshot matrices at some selected parameter values, together with the associated left singular value matrices. For $\mu\in[0,2]$, the values of $f$ are non-negative and so are all entries in the corresponding snapshot matrices. Beyond $\mu=2.0$, negative entries arise in the snapshot vectors. Fig. 8 tracks the smallest singular value $\sigma_{r}(\mu)$ of the snapshot matrices $Y(\mu)$ for $\mu\in[1.7,2.3]$. A substantial non-linear change in $\sigma_{r}(\mu)$ is apparent around the value of $\mu=2.0$. This makes SVD interpolation beyond the parameter location $\mu=2.0$ a challenging problem.

We sample the left singular value matrix $U(\mu)$ together with the derivative $\dot{U}(\mu)$ at 6 Chebychev samples in the interval $[1.7,2.3]$.

As in Section 5.2, we juxtapose the results of quasi-linear interpolation, tangent space interpolation and quasi-cubic Hermite interpolation. For $\mu$ in the sampled range, the relative interpolation errors are computed in the Frobenius norm as $\frac{\|U^{*}(\mu)-U(\mu)\|_{F}}{\|U(\mu)\|_{F}}$, where $U^{*}(\mu)$ denotes the manifold interpolant and $U(\mu)$ is the reference solution. The error curves are displayed in Fig. 9. According to the figure, the tangent space interpolation method fails to interpolate the samples at the first two parameter locations $\mu_{0}\approx 1.7102,\mu_{1}\approx 1.7879$. This is explained as follows. In the tangent space interpolation method, all the Stiefel samples $U(\mu_{i})$ are mapped to the tangent space attached at $U(\mu_{3})$, $\mu_{3}\approx 2.0776$ via $\Delta(\mu_{i})=\operatorname{Log}_{U(\mu_{3})}^{St}(U(\mu_{i})$.888Note that the base point $\mu_{3}$ happens to lie beyond the ‘$2.0$-threshold value’, after which negative function values appear, while $\mu_{0},\mu_{1}<2.0$. It turns out that the Riemannian Stiefel logarithm is not well-defined for $i=0,1$. Put in different words, $U(\mu_{0})$ and $U(\mu_{1})$ are too far from $U(\mu_{3})$ to be mapped to $T_{U(\mu_{3})}St(n,r)$ by the the Stiefel $\log$-algorithm.

The relative Frobenius errors are

[TABLE]

.

5.5 Parametric dimension reduction

SVD interpolation may be used for parametric dimension reduction. Consider samples $Y_{0},\ldots Y_{k}$ of a matrix curve $\mu\mapsto Y(\mu)\subset\mathbb{R}^{n\times m}$. Suppose that instead of the original data, only a low-rank SVD approximation is stored $Y_{i}\approx U_{i}\Sigma_{i}V_{i}^{T}$, where $U_{i}\in\mathbb{R}^{n\times r},\Sigma_{i}\in\mathbb{R}^{r\times r},V_{i}\in\mathbb{R}^{m\times r}$ for a fixed $r\ll\min\{n,m\}$. Then a low-rank approximant to any $Y(\mu)$ in the sampled range can be obtained via interpolating the SVD data.

In this section, we apply this approach to an application from computational option pricing. The value function $y(T,S;K,r,\sigma)$ that gives the fair price for a European call option is determined via the Black-Scholes-equation [8],

[TABLE]

This is a parabolic PDE that depends on time $t$, the stock value $S$, and a number of additional system parameters, namely the strike price $K$, the interest rate $r$, the volatility $\sigma$ and the exercise time $T$. In this experiment, we consider a fixed interest rate of $r=0.01$ and an exercise time of $T=2$ units. The dependency on the underlying $S\in[50,150]$ is resolved via a discretization of the interval by equidistant steps of $\Delta S=0.01$, while the strike price $K\in[30,170]$ is discretized in steps of $\Delta K=1$. Eventually, the volatility $\sigma$ will act as the interpolation parameter. Hermite interpolation requires the option price $y$ as well as its derivative $\partial_{\sigma}y$, in economics referred to as the ‘vega’ of set of the ‘greeks’. A similar test case was considered in [39].

The Black-Scholes equation for a single underlying has a closed-form solution. Yet, here, we will approach it via a numerical scheme in order to mimic the corresponding procedure for real-life problems. Application of a finite volume scheme to the Black-Scholes PDE yields snapshot matrices

[TABLE]

for $\sigma\in[0.1,0.2,\dots,1.0]$. The computation time for each data pair $Y(\sigma),\partial_{\sigma}Y(\sigma)$ is ca. 17min on a standard laptop computer. For each sampled snapshot matrix $Y(\sigma)\in\mathbb{R}^{n\times m}$, $n=10001$, $m=141$, an SVD is performed and is truncated to the $r=5$ dominant singular values/singular vector triples. This yields a compressed representation $Y(\sigma)=U(\sigma)\Sigma(\sigma)V^{T}(\sigma)$, with $U(\sigma)\in St(n,r)$, $\Sigma(\sigma)\in diag(r,r)$, $V(\sigma)\in St(m,r)$ and consumes ca. $0.1s$ on a laptop computer. The relative information content is ric$(r)=\frac{\sum_{j=1}^{r}\sigma_{j}^{2}}{\sum_{k=1}^{m}\sigma_{k}^{2}}\geq 0.99999$. The storage requirements for a $(10001\times 141)$-matrix are 11.3MB, all the low-rank SVD factors truncated to $r=5$ require a total of 0.4MB of disk space, which is ca. $3.5\%$ of the uncompressed representation. We sample full solution data sets $Y(\sigma)$, $\partial_{\sigma}Y(\sigma)$ at $\sigma\in\{0.1,0.4,0.9\}$. The sample data sets are displayed in Fig. 10.

In order to assess the approximation accuracy, we interpolate at $\sigma^{*}\in\{0.2,0.3,\hskip 2.84544pt0.5,0.6,0.7,0.8\}$ and compute the relative Frobenius norm errors of the interpolated low-rank SVD representation $U(\sigma^{*})\Sigma(\sigma^{*})V(\sigma^{*})^{T}$ with respect to the exact full rank data matrix $Y(\sigma^{*})$. We compare the quasi-linear geodesic interpolation (w/o derivative data) to the cubic Hermite approach. The errors are displayed in the bar plot Fig. 11. For completeness, we also include the results of standard linear and Hermite interpolation on the data set of the full, uncompressed matrices $Y_{i}$, where no special geometric structure needs to be addressed.

The values underlying the bar plot are

sigma 0.2 0.3 0.5 0.6 0.7 0.8 Geo. [0.031178 0.028025 0.025577 0.035379 0.032636 0.0201280] Linear full [0.027382 0.022640 0.010562 0.013015 0.010908 0.0062528] Hermite [0.005791 0.004969 0.001416 0.002799 0.002482 0.0009872] Hermite full [0.004298 0.003142 0.000734 0.001374 0.001156 0.0004374]

Mind that the errors for the geodesic and the Hermite low rank interpolation include both the effects of interpolation and data reduction. Even though the data sets underwent a substantial reduction in dimension, the relative errors are of a comparable order of magnitude. For better judging the results, we note that the relative error between the two consecutive samples $Y(0.1)$, $Y(0.4)$ and $Y(0.4)$, $Y(0.9)$ are $27\%$ and $39\%$, respectively.

The technique could be used as a reduced online storage scheme: store the truncated (Hermite) SVD data at some selected sample locations; interpolate, when a prediction at any in-between location is required online.

6 Conclusions and final remarks

We have presented an elementary, general approach to Hermite interpolation on Riemannian manifolds that is applicable to practical problems, whenever algorithms to compute the Riemannian exp and log mappings are available. While our focus was on the manifold counterpart of local cubic Hermite interpolation, the method is flexible and may be combined with any Hermite method that is linear in the sample data. In fact, only the coefficient functions $a_{0},b_{0},b_{1}$ in (10) need to be replaced, no additional changes are necessary. Moreover, combinations of Hermite and Lagrange methods are straightforward generalizations.

In addition, we have exposed a relation between the sectional curvature of the manifold in question the data processing errors, that arise for computations in Riemannian normal coordinates.

As an example, Hermite interpolation of Stiefel data was discussed in more detail. From the observations in the numerical experiments, the main practical constraint on the sampled data is that two consecutive samples be close enough so that the Riemannian Stiefel logarithm is well-defined. As a rule of thumb, if the data points are close enough so that the Riemannian log algorithm converges, then the Hermite interpolation method provides already quite accurate results.

The method constructs piece-wise cubic manifold splines between data points $p_{i}$ and $p_{i+1}$ in terms of normal coordinates centered at $p_{i+1}$. Thus, it is not symmetric in the sense that computations in normal coordinates centered at $p_{i}$ might lead to different results. Yet, in the numerical experiments, these effects prove to be negligible.

Acknowledgments

The data set featured in Section 5.5 was kindly provided by my colleague Kristian Debrabant from the Department for Mathematics and Computer Science (IMADA), SDU Odense.

Appendix A The basic cubic Hermite coefficient polynomials

The coefficient functions in (6) are the cubic Hermite polynomials that are uniquely defined by

[TABLE]

The explicit cubic coefficient functions are

[TABLE]

and are displayed in Fig. 12 for $t_{0}=0,t_{1}=1$. Since on manifolds, we work exclusively in the setting, where $q=0$, the coefficient $a_{1}(t)$ drops out in (6).

Fig. 13 shows the spatial cubic Hermite spline (6) that connects the points $p=(1,0,0),q=(0,0,0)\in\mathbb{R}^{3}$ with a prescribed start and terminal velocity of $v_{p}=(0.5,0.5,0)$ and $v_{q}=(0,0,1)$, respectively.

Appendix B Differentiating the QR-decomposition

Let $t\mapsto T(t)\in\mathbb{R}^{n\times r}$ be a differentiable matrix function with Taylor expansion $T(t_{0}+h)=T(t_{0})+h\dot{T}(t)$. Following [34, Proposition 2.2], the QR-decomposition is characterized via the following set of matrix equations.

[TABLE]

In the latter, $P_{L}=\left(\begin{smallmatrix}0&0&\cdots&0\\ 1&\ddots&\ddots&0\\ \vdots&\ddots&\ddots&\vdots\\ 1&\cdots&1&0\end{smallmatrix}\right)$ and ‘$\odot$’ is the element-wise matrix product so that $P_{L}\odot R$ selects the lower triangle of the square matrix $R$. For brevity, we write T=T(t_{0}),\dot{T}=\frac{d}{dt}\big{|}_{t=t_{0}}T(t), likewise for $Q(t)$, $R(t)$. By the product rule

[TABLE]

According to [34, Proposition 2.2], the derivatives $\dot{Q},\dot{R}$ can be obtained from Alg. 1. The trick is to compute $X=Q^{T}\dot{Q}$ first and then use this to compute $\dot{Q}=QQ^{T}\dot{Q}+(I-QQ^{T})\dot{Q}$ by exploiting that $Q^{T}\dot{Q}$ is skew-symmetric and that $\dot{R}R^{-1}$ is upper triangular.

Appendix C Differentiating the singular value decomposition

Let $m\leq n\in\mathbb{N}$ and suppose that $t\mapsto Y(t)\in\mathbb{R}^{n\times m}$ is a differentiable matrix curve around $t_{0}\in\mathbb{R}$. If the singular values of $Y(t_{0})$ are mutually distinct, then the singular values and both the left and the right singular vectors depend differentiable on $t\in[t_{0}-\delta t,t_{0}+\delta t]$ for $\delta t$ small enough. This is because the associated symmetric eigenvalue problem $Y^{T}(t)Y(t)=V(t)\Lambda(t)V^{T}(t)$ is differentiable under these (and more relaxed) conditions, [3].

Let $t\mapsto Y(t)=U(t)\Sigma(t)V(t)^{T}\in\mathbb{R}^{n\times m}$, where $U(t)\in St(n,m)$, $V(t)\in O(m)=St(m,m)$ and $\Sigma(t)\in\mathbb{R}^{m\times m}$ diagonal and positive definite. Let $u_{j}$ and $v_{j}$, $j=1,\ldots,m$ denote the columns of $U(t_{0})$ and $V(t_{0})$, respectively. For brevity, write Y=Y(t_{0}),\dot{Y}=\frac{d}{dt}\big{|}_{t=t_{0}}Y(t), likewise for the other matrices that feature in the SVD.

The above algorithm is mathematical ‘folklore’, a proof can be found in, e.g., [18]. Note that $U^{T}\dot{U}$ with $\dot{U}$ as above is skew-symmetric, so that indeed $\dot{U}\in T_{U}St(n,m)$. The above equations make use of the inverse $\Sigma^{-1}$ and are therefore unstable, if the singular values are small. This effect can be alleviated by truncating the SVD to the $r\leq m$ dominant singular values. The derivative matrices for the truncated SVD are stated in Alg. 3.

Since this algorithm is based on representing the derivative vectors $\dot{v}_{j}$ in terms of an eigenvector ONB $V=(V_{r},V_{m-r})=(v_{1},\ldots,v_{r},v_{r+1},v_{m})$, a full square orthogonal $V$ is required also in the truncated case. Yet, note that the columns of $V_{m-r}$ feature only in the equation for $\dot{V}_{r}=V\Gamma=V_{r}\Gamma_{r}+V_{m-r}\Gamma_{m-r}$ while all other quantities are independent of $V_{m-r}$.

If the rank of $Y\in\mathbb{R}^{n\times m}$ is exactly $r\leq m$ and is fixed for all $t$, then the computation of the entries of the lower block $\Gamma_{m-r}$ reduces to $\Gamma_{ij}=\frac{u_{j}^{T}\dot{Y}v_{i}}{\sigma_{j}},i=r+1,\ldots,m,j=1,\ldots,r$. In this case, the singular value matrix features a lower-right zero diagonal block $\Sigma=\mathop{\operator@font diag}\nolimits(\sigma_{1},\ldots,\sigma_{r},\sigma_{r+1},\ldots,\sigma_{m})$. In general, computing the derivatives in the presence of multiple singular values/eigenvalues is sophisticated [3]. Here, however, it is sufficient to compute the singular vectors $V_{r}(t)=(v_{1}(t),\ldots,v_{r}(t))$ associated with the pairwise distinct singular values and to perform a $t$-dependent orthogonal completion $V(t)=(v_{1}(t),\ldots,v_{r}(t),v_{r+1}(t),\ldots,v_{m}(t))$ via the modified Gram-Schmidt process, which is differentiable.

Appendix D The Riemannian Stiefel log algorithm

All numerical experiments featured in this work are performed with a SciPy[15] implementation of the following algorithm, for the details, see [37].