Smoothing splines on Riemannian manifolds, with applications to 3D shape space

There has been increasing interest in statistical analysis of data lying in manifolds. This paper generalizes a smoothing spline fitting method to Riemannian manifold data based on the technique of unrolling and unwrapping originally proposed in Jupp and Kent (1987) for spherical data. In particular we develop such a fitting procedure for shapes of configurations in general $m$-dimensional Euclidean space, extending our previous work for two dimensional shapes. We show that parallel transport along a geodesic on Kendall shape space is linked to the solution of a homogeneous first-order differential equation, some of whose coefficients are implicitly defined functions. This finding enables us to approximate the procedure of unrolling and unwrapping by simultaneously solving such equations numerically, and so to find numerical solutions for smoothing splines fitted to higher dimensional shape data. This fitting method is applied to the analysis of simulated 3D shapes and to some dynamic 3D peptide data. Model selection procedures are also developed.


Introduction
Analysis of temporal shape data has become increasingly important for applications in many fields, where the shape of an object is what is left after removing the effects of rotation, translation and re-scaling (Kendall, 1984). For an introduction to the statistical analysis of shape see Dryden and Mardia (2016); Kendall et al. (1999); Patrangenaru and Ellingson (2016) and Srivastava and Klassen (2016). Suppose that we are given time series of data which are measurements of a moving object. One may wish to find the overall trend or be interested in the behaviour of the object at unobserved times, including the explanation of how the shapes change between successive times or prediction in the future. For example, in the field of molecular biology, the study of dynamic proteins is of interest.
The main difficulty in applying the classical statistical approach directly to analyse changes of shapes of configurations of k labelled landmarks over time lies in the fact that spaces of such shapes are curved, so that classical methods are not always appropriate. In the case where there is small variability in the set of shape data, one way to overcome this difficulty is to carry out the classical techniques on the projection of the data onto the Procrustes tangent space at their mean. See, for example,  for Procrustes tangent projections, and  and Morris et al. (1999) for their applications. The authors in  applied such projections to two dimensional configuration data to construct models for the growth of shapes using roughness penalties, while the authors in Morris et al. (1999) combined such projections with principal component analysis to study the growth with age of the size-andshapes of two dimensional profiles of human faces.
For data lying on a general Riemannian manifold and with large variability, geodesic based fitting approaches generalize the regression model for Euclidean data. Several different methods for fitting geodesics have been proposed. An early geodesic based model for detecting shape changes over time in Kendall shape space for planar configurations was proposed in Le and Kume (2000). Fletcher et al. (2004); Fletcher (2013) introduced principal geodesic analysis for data lying on a manifold, which uses the first principal component of the image of data under the inverse exponential map at their mean to determine the fitted geodesic for the given data. In Kenobi, Dryden and Le (2010), the authors investigated a method for fitting minimal geodesics in Kendall shape space based on minimizing sums of squares of Procrustes distances. A method of principal geodesic analysis based on intrinsic methods of fitting geodesics to data on Riemannian manifolds was developed in Huckemann and Ziezold (2006) and Huckemann, Hotz and Munk (2010).
It is also possible to generalize Euclidean smoothing spline fitting methods (Green and Silverman, 1994) to manifold-valued data. For example, Su et al. (2012) discussed the problem of fitting smoothing splines to data points on a Riemannian manifold measured over time using a Palais metric-based gradient method. To be able to implement this method in practice, complete knowledge of the geometric structure of manifolds is required. This is usually difficult for Kendall shape spaces of configurations lying in three dimensional Euclidean space.
The authors in Jupp and Kent (1987) proposed a method for fitting spherical smoothing splines to spherical data based on the techniques of unrolling and unwrapping onto an appropriate tangent space. The advantage of this approach is that, on the one hand, it retains certain geometric features of the data, such as lengths and angles, and on the other hand it does not require as much geometric knowledge as that in Su et al. (2012). This smoothing spline fitting method was generalized to two dimensional shape data in Kume, Dryden and Le (2007), using the result obtained in Le (2003) for unrolling procedures in general Kendall shape spaces. However, the study of shapes of configurations in three or more dimensions is much more complicated and, as noted in Kume, Dryden and Le (2007), the method obtained there cannot be generalized to three dimensional shape data.
In this paper, we propose a smoothing spline fitting method on general manifolds, based on the techniques of unrolling and unwrapping, and then focus on fitting smoothing shape splines to m-dimensional time series shape data for m 3. Although we are unable to obtain a closed expression for parallel transport, which underpins the procedure of unrolling and unwrapping, we show that parallel transport along a geodesic on Kendall shape space is linked to the solution of a homogeneous first-order differential equation, some of whose coefficients are implicitly defined functions. This enables us to approximate the procedure of unrolling and unwrapping by simultaneously solving such equations numerically, and hence obtaining numerical solutions for smoothing splines fitted to three dimensional shape data. This paper is organized as follows. In Section 2, we reformulate the procedures of unrolling and unwrapping, and then use them to define smoothing splines fitting on general manifolds. Since, as is clear from Section 2, the procedures of unrolling and unwrapping use the concept of so-called parallel transport, Section 3 provides some necessary basic geometry of shape space of labelled configurations in R m to develop parallel transport for general m for practical use. Numerical simulations and an application to real 3D peptide data are given in Section 4.

Unrolling, unwrapping and wrapping
Recall that two of the main ingredients of the spline fitting technique used in Jupp and Kent (1987) and Kume, Dryden and Le (2007) are respectively (i) unrolling a curve in the relevant space to the tangent space at its starting point; and (ii) unwrapping points in the space at known times, with respect to a base path, onto the tangent space at the starting point of the path. Thus, to generalize spline fitting in these two papers to general manifolds, we first reformulate these two procedures in a general manifold setting.
For this, let M be a complete and connected Riemannian manifold with induced Riemannian distance d and denote by τ x (M ) the tangent space to M at x ∈ M . Then, the exponential map on M at x can be expressed by where γ is the unit-speed geodesic such that γ(0) = x andγ(0) = v/ v . Thus, the inverse of the exponential map can be expressed as when there is a unique shortest unit-speed geodesic γ from x to x .
Intuitive descriptions of the unrolling, unwrapping and wrapping procedures that we shall define can be found in Jupp and Kent (1987) and Kume, Dryden and Le (2007). Formally, they are closely linked to the concept of parallel transport along a curve. The latter is a method of transporting tangent vectors along smooth curves in a manifold and, in some sense, provides a method of isometrically moving the local geometry of a manifold along curves. More precisely, suppose that γ(t), 0 t t * is a smooth curve in M and that v 0 ∈ τ γ(0) (M ), where t * > 0 is a fixed constant. Then, the parallel transport of v 0 along γ is the extension of v 0 to a vector field v on γ such that where ∇ denotes the covariant derivative on M . This provides linear isomorphisms between the tangent spaces at points along the curve γ: This isomorphism is known as the parallel transport map associated with the curve. Note that P . Then, as explained in Kume, Dryden and Le (2007), in terms of the parallel transport P γ(t) γ(s) , • the unrolling of γ onto τ γ(0) (M ) is the curve γ † on τ γ(0) (M ) such that γ † (0) = 0 and • the unwrapping at t ∈ [0, t * ], with respect to γ, of a point x ∈ M into τ γ(0) (M ) is the sum of the two tangent vectors in τ γ(0) (M ): γ † (t) and the parallel transport of the tangent vector exp −1 γ(t) (x) ∈ τ γ(t) (M ) along γ to τ γ(0) (M ), i.e. the unwrapping at t of x, with respect to γ, is • the wrapping at t, with respect to γ, of a tangent vector v ∈ τ γ(0) (M ) back into M is the reverse of the unwrapping procedure.
In particular, if γ is a geodesic, then γ † is the linear segment of the same length as γ and in the same direction as the initial tangent vector of γ, that is, it is the linear segment from the origin of τ γ(0) (M ) to exp −1 γ(0) (γ(t * )).

Smoothing spline fitting on manifolds
Turning to smoothing splines on manifolds, we recall that, for v 0 , v 1 , . . . , v n ∈ R d observed at times, t 0 , t 1 , . . . , t n respectively, the cubic smoothing spline estima- among all C 2 -functions, where T = [t 0 , t n ] and λ > 0 denotes a smoothing parameter. This minimization problem can be solved on each component f † (j) of f † separately to reduce the d-dimensional minimization problem to d one dimensional ones, so that The solution to such a minimization problem is a cubic spline with knots at the unique values of v i , while the smoothing parameter λ controls the trade-off of the model complexity. When λ tends to zero, f † λ becomes any function interpolating the data points. On the other hand, f † λ becomes a classical simple linear regression line when λ tends to infinity. Although we could use a penalty with other orders of derivatives, throughout this paper we use the integrated squared second derivative as the penalty function which leads to the cubic smoothing spline which often works well in practice.
In a similar fashion to that in Jupp and Kent (1987) and Kume, Dryden and Le (2007), we use the concept of unrolling and unwrapping to define M -valued smooth splines as follows.
Definition 1. For a given dataset D = {x j : 0 j n} in M , where x j is observed at time t j , and smoothing parameter λ, we define the M -valued smoothing spline fitted to D with parameter λ to be the C 2 -function To find the smoothing spline fitted to a given dataset D, the iterative scheme for approximation given in Jupp and Kent (1987) can be applied to general manifolds: take the piecewise geodesic passing through the data points as the initial curve f 0 ; at each step n, unwrap D with respect to f n to get D † n in τ fn(t0) (M ); fit a Euclidean smoothing spline f † n+1 to D † n in τ fn(t0) (M ); then f n+1 is the wrapped path, with respect to f n , of f † n+1 in M . The iterative procedure stops when f n and f n+1 are sufficiently close.
In such an iterative scheme, smooth curves are approximated by piecewise geodesics. Thus, for the implementation in practice, it is sufficient to restrict ourselves to the procedures of unrolling a piecewise geodesic and of unwrapping and wrapping with respect to a piecewise geodesic. Then for the given curve γ(t), t 0 t t n , on M , which is a geodesic between t i and t i+1 , where t 0 < t 1 < · · · < t n , the above unrolling, unwrapping and wrapping procedures along γ can be formulated as follows.
• The unrolling of γ onto τ γ(t0) (M ) is the piecewise linear segment in the tangent space τ γ(t0) (M ) joining the following successive points: given by as long as x is not in the cut locus of γ(t).
• The wrapping at time t ∈ [t 0 , t n ], along γ, of a tangent vector v ∈ τ γ(t0) (M ) back into M is the point x t ∈ M given by It is clear that, to be able to carry out the unrolling, unwrapping and wrapping procedures in practice, the crucial steps are to find the expressions for the exponential map, as well as its inverse, and for the parallel transport along a geodesic. For example, these expressions are known in the case of the unit sphere S d : • a unit speed geodesic starting from x ∈ S d takes the form where v ∈ R d+1 such that v = 1 and x, v = 0; • the exponential map at x has the expression for all x ∈ S d \ {−x}; • the parallel transport of the tangent vector w ∈ τ γ(t0) (S d ) along the geodesic γ defined by (2) is the vector field w(t) along γ(t) given by

Smoothing splines in shape space
3.1. Shape space, its tangent space and shape geodesics Due to their non-trivial geometric structure, the direct implementation on the Kendall shape spaces for configurations in R m (m 3) of the exponential map and, in particular, of the parallel transport turns out to be challenging. One possible way to overcome this difficulty is to explore the fact that Kendall shape spaces are the quotient spaces of Euclidean spheres and to use the much simpler structure of spheres, as has previously been done by many in various different statistical investigations in shape analysis. For example, Le (2003) obtained an equivalent, qualitative description of the parallel transport on shape spaces via the pre-shape sphere. In the case of shapes of configurations in R 2 , this description leads successfully to a closed expression for such an equivalent notion on the pre-shape sphere, which has then been used for statistical analysis of shape in Kume, Dryden and Le (2007). Unfortunately, as pointed out in Le (2003), such a closed expression is generally unavailable on the shape space of configurations in R m for m 3. In this section, we summarize the facts required for the exponential maps on the shape space of labelled configurations in R m and develop the result in Le (2003) further for general m to the extent that it is usable for practical purposes.
Recall that, for a configuration in R m with k(> m) labelled landmarks, its pre-shape is what is left after the effects of translation and scaling are removed and that this pre-shape can be represented by an m × (k − 1) matrix X with X 2 = tr(XX ) = 1. The space S k m of such pre-shapes is known as the preshape space of configurations of k labelled landmarks in R m and is the unit sphere of dimension m(k −1)−1, i.e. S k m = S m(k−1)−1 . The Kendall shape space Σ k m of configurations with k labelled landmarks in R m is then the quotient space of S k m by the rotation group SO(m) acting on the left (cf. Dryden andMardia (2016) &Kendall et al. (1999)). We shall use [X] to denote the shape of the pre-shape X. Writing . The horizontal subspace, which is the same as the Procrustean tangent space, of τ X (S k m ), with respect to the quotient map from S k m to Σ k m , can be expressed as of X, and this gives a useful isometric representation of the tangent space of Σ k m at [X] for the statistical analysis of shapes (cf. Kendall et al., 1999).
Similarly, a unit-speed geodesic γ in the shape space starting from [X] can be isometrically represented by a so-called horizontal geodesic in S k m of the form Γ(s) = X cos s + V sin s, s ∈ [0, π/2), where V ∈ H X (S k m ) and V = 1. The path Γ is usually referred to as the horizontal lift of γ. To obtain a representation of a shortest geodesic between two shapes [X 1 ] and [X 2 ], we take the two pre-shapes X 1 and X 2 of [X 1 ] and [X 2 ] respectively such that X 1 X 2 is symmetric and all its eigenvalues are non-negative except possibly for λ m , the smallest one, where sign(λ m ) = sign(det(X 1 X 2 )). Then, a unit-speed shortest geodesic from [X 1 ] to [X 2 ] can be isometrically represented by the horizontal geodesic from X 1 to X 2 : where Note that V X1,X2 = 1. Thus, it follows from the expression (3) for the horizontal geodesic from X 1 to X 2 that, if there is a unique geodesic between [X 1 ] and [X 2 ], the inverse exponential map exp −1 [X1] ([X 2 ]) on the shape space Σ k m can be isometrically represented by its horizontal lift on S k m given by exp −1 X1 (X 2 ) = s 0 V X1,X2 .

Parallel transport
Turning to parallel transport, it is shown in Le (2003) that, along a horizontal geodesic Γ(s) in S k m , a vector field V (s) is horizontal and its projection to τ [Γ(s)] (Σ k m ) is the parallel transport, along the shape geodesic [Γ(s)], of the projection of V (0) to τ [Γ(0)] (Σ k m ) if and only if V (s) satisfies the following three conditions: Clearly, the usefulness of this result for practical implementation lies crucially in obtaining a more quantitative description of the skew-symmetric matrix A(s) in (6), which is answered by the following result.
where A(s) is skew-symmetric and is the unique solution to Note that V (s) given in Proposition 2 is generally not the parallel transport of V along Γ(s) in S k m . Proof. Note first that V ∈ H Γ(0) (S k m ) implies that tr(V Γ(0) ) = 0 and V Γ(0) is a symmetric matrix.
The uniqueness of the solution to (8) follows from the fact that, for a given m × m symmetric non-negative definite matrix S = RΛR of rank at least m−1, where R ∈ O(m) and Λ = diag{λ 1 , . . . , λ m }, and a given skew-symmetric matrix B, there is a unique skew-symmetric matrix A satisfying the equation AS + SA = B. To see the latter, writeÃ = (ã ij ) = R AR andB = (b ij ) = R BR. Then, AS + SA = B if and only ifÃΛ + ΛÃ =B, which is if and only if λ iãij + λ jãij =b ij for i < j, i.e. if and only ifã ij =b ij /(λ i + λ j ).
The result of Proposition 2 provides a numerical method to approximate the horizontal vector field V (s), along a horizontal geodesic Γ(s), whose projection to τ [Γ(s)] (Σ k m ) is the parallel transport, along the shape geodesic [Γ(s)], of the projection of V (0) to τ [Γ(0)] (Σ k m ) in the usual way: with step size δ, where A(iδ) is updated at each step using (8). However, V ((i + 1)δ) obtained in this way is generally neither tangent to S k m nor horizontal. Hence, this requires finding, at each step, the re-normalized horizontal component of V ((i + 1)δ). That is, at each step, we need to project V ((i + 1)δ) to τ Γ((i+1)δ) (S k m ) and then to the horizontal subtangent space of S k m at Γ((i + 1)δ), followed by normalizing it so that its length remains equal to that of V (0).
If X is a pre-shape matrix with rank(X) m − 1, the projection of V ∈ M (m, k − 1) onto the tangent space τ X (S k m ) is given by To obtain the projection of V ∈ τ X (S k m ) onto H X (S k m ), we note that the orthogonal complementary subspace to H square matrix with all its elements zero except for the (i, j)th element which is one, {A ij X/ A ij X | 1 i < j m} forms an orthonormal basis for V X . Thus, the projection of V ∈ τ X (S k m ) onto H X (S k m ) is given by

Shape spline fitting algorithm
We are now in a position to construct an algorithm to calculate the shapespace spline, or simply the shape spline, for a given set of configurations in R m , observed at times t 0 , . . . , t n , whose pre-shapes are denoted by X 0 , . . . , X n , in a similar fashion to that given in Kume, Dryden and Le (2007), as follows. Let δ and ε be two given small positive numbers.
1. Rotate the successive pre-shapes, X i+1 , such that the resulting X i+1 is the Procrustes fit of the original one onto X i . That is, rotate the successive pre-shapes X i+1 to satisfy the conditions that X i X i+1 is symmetric and that all its eigenvalues are non-negative except possibly for the smallest one whose sign is the same as the sign of det X i X i+1 . Denote the set of the resulting pre-shapes by D = {X i : 0 i n} and let Γ D be the initial piecewise horizontal geodesic of D such that Γ D (t i ) = X i . 2. Construct G grid points between successive two times such that t 0 = t 00 < t 01 < t 02 < . . . < t 0G < t 1 = t 10 < . . . < t ij < . . . < t n = t n0 which gives t ij , i = 0, 1, . . . , n, j = 0, 1, . . . , G for i n − 1 and j = 0 for i = n where the difference between successive t ij is less than or equal to δ. 3. Set the base path Γ D1 to be the piecewise horizontal geodesic passing through D 1 = {Γ D (t ij ) : ∀i, j}. 4. Using the numerical procedure, including the required necessary projections specified following the proof of Proposition 2, unwrap the data D into τ Γ D 1 (t0) (S k m ) with respect to the base path Γ D1 which gives 5. Fit the smoothing spline (1) to D † and find fitted values at the times of the grid points given in Step 2, giving 6. Using the numerical procedure, including required necessary projections specified following the proof of Proposition 2, wrap D † 2 back into the preshape sphere with respect to the base path Γ D1 , giving D 2 := {Z ij : i, j} ⊂ S k m . 7. Successively rotate Z ij in D 2 such that the resulting Z i j+1 is the Procrustes fit of the original one onto Z ij , and obtain the piecewise horizontal geodesic Γ D2 passing through the resulting Z ij . Then, Γ D2 becomes the base path in the next iteration.
ε, replace D 1 by D 2 and go to Step 9. Otherwise, stop the iterations. 9. Successively rotate X i in D such that the resulting X i is the Procrustes fit of the original one onto Z i0 . Update D by the resulting X i and go to Step 4.

Size-and-shape splines
All these results have analogues in the size-and-shape setting. For example, for a configuration in R m with k(> m) labelled landmarks, its pre-size-and-shape is what is left after the effects of translation are removed. This pre-size-and-shape can be represented by an m × (k − 1) matrixX ∈ M (m, k − 1) and the space of the pre-size-and-shapes, SS k m , is identical with M (m, k − 1). ForX = 0, the tangent space τX (SS k m ) is then also M (m, k − 1) and the horizontal subspace of τX (SS k m ), with respect to the quotient map from SS k m to the Kendall sizeand-shape space SΣ k m , is given by Horizontal geodesics in SS k m starting fromX take the form where V ∈ HX (SS k m ). For two given size-and-shapes [X 1 ] and [X 2 ], letX 1 and X 2 be the pre-size-and-shapes of [X 1 ] and [X 2 ] respectively such thatX 1X 2 is symmetric and all its eigenvalues are non-negative except possibly for λ m , the smallest one, where sign(λ m ) = sign(det(X 1X 2 )). Then,X 2 −X 1 ∈ HX (SS k m ) and a shortest geodesic from [X 1 ] to [X 2 ] can be represented by the horizontal geodesic that connectsX 1 andX 2 : Thus, the inverse exponential map exp −1 [X1] ([X 2 ]) on the size-and-shape space SΣ k m , using its horizontal lift on SS k m , is isometrically represented by as long as the geodesic between [X 1 ] and [X 2 ] is unique. Moreover, a modification of the argument given in Le (2003) shows that, along a horizontal size-and-shape geodesicΓ(s) in SS k m , the vector field V (s) is horizontal and its projection to τ Thus, it follows from the proof for Proposition 2 that the result of Proposition 2 can also be generalized to size-and-shape space.
where A(s) is skew-symmetric and is the unique solution to (8) with Γ being replaced byΓ.

Numerical study
In this section, we illustrate our proposed method by applying it to simulated data as well as to real data. For notational convenience we designate the starting index by 1 instead of 0, i.e. the n data points are measured at times t 1 , . . . , t n .
Our simulation and real data analyses have been implemented in R. The builtin function smooth.spline() was used for smoothing splines and the package shapes was used for the Procrustes analysis and the relevant shape functions (Dryden, 2017). In both simulated and real data, we used G = 2 for a smooth base path and ε = 10 −5 for convergence. The choice of the smoothing parameter λ, adapted to the data, is determined by applying the (Euclidean) leave-oneout cross-validation method (Efron and Tibshirani (1993), Chapter 17) to the unrolled shape spline and the unwrapped shape data. That is, we chose the optimal λ which minimizes where f † λ,−i is the unrolling of f λ,−i to the tangent space at its starting point, f λ,−i is the smoothing shape spline with the ith observation excluded and v i is the unwrapped ith observed shape data with respect to f λ,−i . The candidates for λ in the cross-validation were taken from the set {10 −9 , 10 −7 , 10 −5 , 10 −3 , 10 −1 }. Consequently in the following, the chosen smoothing parameters for the original and perturbed data in the simulation study were 10 −5 and 10 −3 , respectively, and that for the peptide data analysis was 10 −3 . The step size for the approximation in parallel transport was calculated by allowing two equally spaced steps between every two successive shapes.

Simulation
We consider configurations in R 3 with k = 8 landmarks. Four initial vertex configurations are chosen such that the Riemannian shape distances between successive configurations are 0.47, 0.75 and 0.54 respectively. Note that the Riemannian shape distance ranges over [0, π/2], so that the shapes of these consecutive configurations are not relatively close. Then four equally spaced shapes on each of the three geodesic segments connecting the successive shapes of the vertex configurations are generated; the configurations of these 16 shapes are shown in Figure 1, where the 1st, 6th, 11th and 16th are the four initial configurations. Then, we add Gaussian noise with standard deviation 0.05 independently to each landmark of the original configurations to generate perturbed data which is shown in Figure 2. Figure 3 shows the configurations of the fitted shapes on the smoothing shape spline to the perturbed data. The Riemannian shape distances between the original data and its fitted data (not displayed) are displayed in Figure 4(a), while the corresponding distances between the perturbed data and its fitted values in Figure 3 are displayed in Figure 4(b). It shows that there are extremely small residuals for the original data while, on the other hand, for the perturbed data there are relatively large residuals at some points, such as at points 5, 6, 9, 15. To visualize the fitted smoothing shape splines to the original and perturbed data respectively, in each case we unroll the given 16 shapes onto the tangent space to the first shape along the piecewise geodesic segments through them and unroll the fitted smoothing shape spline onto the tangent space to its starting point. Then, we carry out the principal component analysis on the coordinates of these four unrolled data sets respectively. Figure 5 shows the first two principal component scores of the unrolled original data and the corresponding fitted shape spline, in (a), and the unrolled perturbed data and the corresponding fitted shape spline, in (b). The given data are displayed as squares and the dotted lines are for the fitted shape spline with 's' indicating the starting point. For the fitted shape spline to the original data set, the proportion of variance explained by PC1 and PC2 is 80.2% and 14.1% respectively and, for the fitted shape spline to the perturbed data, 79.5% and 12.7% are explained by PC1 and PC2. As the original data yields three consecutive geodesic segments, the principal component scores are allocated in a corresponding geometrical manner as expected. Noting that the two turning points of these three geodesic segments are at the 6th and 11th configurations, Figure 5(a) clearly shows the structure of the three geodesic segments which is in accord with the simulation setting. For the perturbed data set, the overall pattern is analogous to that of the original    data even though the noise level of standard deviation 0.05 dictates the local curves in the fitted shape spline. Note that the results at the 5th and 6th points for the perturbed data are consistent with the large residuals in Figure 4.

Peptide data
In molecular biology, it is of great interest to study the changes in shape of a protein, as shape is an important component of a protein's function. Studies of this type are often approached through molecular dynamics simulations, which are large deterministic simulations using Newtonian mechanics to model the movement of a protein in a box of water molecules (e.g. Salomon-Ferrer, Case and Walker, 2013). We consider a dataset of in the study of the small alanine pentapeptide (Ala 5 ) which consists of k = 29 atoms in R 3 at 10000 equal picosecond (10 −12 s) time intervals. The data were provided by Professor Charles Laughton, School of Pharmacy, University of Nottingham, and further detail on this particular peptide is given by Margulis, Stern and Berne (2002). We apply the smoothing shape spline fitting to a small temporal subsequence of 10 equally spaced configurations, which are each 1 nanosecond (10 −9 s) apart, and we consider several different models. Planar projections of these configurations are presented in Figure 6 as rainbow coloured dots where red and violet points indicate the first and last landmarks, respectively. The peptide configuration at the start of this sequence is very irregular in all three dimensions at t = 1, then it gradually straightens over time so that we see a more smoothly curved form at t = 10. For the value 10 −3 of λ chosen by cross validation, in Figure 6 the fitted configurations to the moving configurations are indicated by '×' with connecting solid lines, where the maximum Riemannian shape distance between base shapepath for the last two iterations for finding the smoothing shape spline is less than 10 −5 , which is small enough to assume convergence.
We consider three different models for fitting the same data, one of which is a smoothing shape spline model with a sufficiently small smoothing parameter. The others involve strict shape-geodesic fitting, assuming one geodesic and two geodesics respectively, the first by setting the smoothing parameter for the shape spline as infinity which in effect becomes a simple linear regression model in the base tangent space. For the two geodesics model, we use the total sum of squares of errors (SSE) as the criterion to determine the change point for the two geodesic segments. In our example, the resulting change point is t = 5. To investigate which of these three models gives the best fit, we test between the (j) t = 10 hypotheses: H 0 : One geodesic model γ 0 (t i ) † , i = 1, . . . , n, with p 0 = 2(3k − 7) parameters; H 1 : Two geodesics model γ 11 (t i ) † , i = 1, . . . , n 1 , γ 12 (t i ) † , i = n 1 + 1, . . . , n, with p 1 = 4(3k − 7) + 1 parameters; parameters.
Let X i be the pre-shape of the moving peptide configuration at time i and X † 0i , X † 1i and X † 2i be respectively the unwrapped pre-shapes of X i , with respect to the horizontal lifts of the fitted paths for these three models, onto the tangent spaces at the starting points of the paths. We perform F -tests for adequateness of models by letting whereΓ i are the horizontal lifts ofγ i . For H 0 versus H 1 , we consider the test statistic for H 0 versus H 2 , and, for H 1 versus H 2 , where X · ∼ F α,β denotes that X is approximately distributed as the distribution F with degrees of freedom α and β (assuming an independent isotropic Gaussian model for the errors). In our peptide data, k = 29, n = 10, p 0 = 160, p 1 = 321, p 2 = 800, and SSE j , j = 0, 1, 2, are given in Table 1. Hence, the test statistics are F 01 = 3.7043, F 02 = 3.6268 and F 12 = 0.6098, respectively. For H 0 versus H 1 , since the p-value is P(F ≥ 3.7043) < 0.001, F · ∼ F 161,321 , we see that the two geodesics model is a better fitting method for the data than the one geodesic model. For H 0 versus H 2 , the p-value is P(F ≥ 3.6268) < 0.001, F · ∼ F 640,800 , so that the smoothing spline model outperforms the one geodesic model. On the other hand, for H 1 versus H 2 , since p-value = P(F ≥ 0.6098) ≈ 1, F · ∼ F 479,800 , the data do not support the smoothing spline model against the two geodesics model. This partially indicates that the estimated change point works well as a breaking point for two separate geodesics. Figure 7 shows the first three shape principal component scores of the ten peptide configurations in the Procrustean tangent space at their mean shape. The connected dots indicate the PC scores of the data points and the red lines on the other hand are the PC scores of the fitted paths, where their starting points are indicated by 's' in all panels. The curved PC score curves in (a) and (b) are as expected as the mean shape of these peptides does not lie in either of these fitted geodesic segments. As shown in (a), the one geodesic model does not explain the data points well in terms of the PC scores. On the other hand, the others show better results.
Note that there are many other applications of smoothing splines on Riemannian manifolds (e.g., see Su et al., 2012). The main advantage of the unwrapping and unrolling method over a simple tangent space method is that much larger variations in data can be handled, such as is present in the peptide application. Also, unlike Su et al. (2012) we do not need to know the complete geometrical structure of the manifold, and hence the range of potential applications is very broad indeed and will be explored in further work.