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 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, Kent and Mardia (2001) for Procrustes tangent projections, and Kent et al. (2001) and Morris et al. (1999) for their applications. The authors in Kent et al. (2001) 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-and-shapes 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).
Various methods for fitting regression models on manifolds with Euclidean covariates have been recently introduced, including regression on Riemannian symmetric spaces (Cornea et al., 2017), intrinsic polynomial regression (Hinkle, Fletcher and Joshi, 2014), extrinsic local regression (Lin et al., 2017), global and local Fréchet regression (Petersen and Müller, 2017)
, and intrinsic and varying coefficient regression models for diffusion tensors(Yuan et al., 2013; Zhu et al., 2009).
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 -dimensional time series shape data for . 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 to develop parallel transport for general for practical use. Numerical simulations and an application to real 3D peptide data are given in Section 4.
2 Smoothing splines on manifolds
2.1 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 unrolling a curve in the relevant space to the tangent space at its starting point; and 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 be a complete and connected Riemannian manifold with induced Riemannian distance and denote by the tangent space to at . Then, the exponential map on at can be expressed by
where is the unit-speed geodesic such that and . Thus, the inverse of the exponential map can be expressed as
when there is a unique shortest unit-speed geodesic from to .
. 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, is a smooth curve in and that , where is a fixed constant. Then, the parallel transport of along is the extension of to a vector field on such that
where denotes the covariant derivative on . 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 .
Then, as explained in Kume, Dryden and Le (2007), in terms of the parallel transport ,
the unrolling of onto is the curve on such that and
the unwrapping at , with respect to , of a point into is the sum of the two tangent vectors in : and the parallel transport of the tangent vector along to , i.e. the unwrapping at of , with respect to , is
the wrapping at , with respect to , of a tangent vector back into 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 to .
2.2 Smoothing spline fitting on manifolds
Turning to smoothing splines on manifolds, we recall that, for observed at times,
respectively, the cubic smoothing spline estimator for this dataset inis the function that minimizes
among all -functions, where and denotes a smoothing parameter. This minimization problem can be solved on each component of separately to reduce the -dimensional minimization problem to one dimensional ones, so that
for . The solution to such a minimization problem is a cubic spline with knots at the unique values of , while the smoothing parameter controls the trade-off of the model complexity. When tends to zero,
becomes any function interpolating the data points. On the other hand,
becomes a classical simple linear regression line whentends 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.
For a given dataset in , where is observed at time , and smoothing parameter , we define the -valued smoothing spline fitted to with parameter to be the -function
such that its unrolling onto is the cubic smoothing spline fitted to the data obtained by unwrapping at times , with respect to , into .
To find the smoothing spline fitted to a given dataset , 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 ; at each step , unwrap with respect to to get in ; fit a Euclidean smoothing spline to in ; then is the wrapped path, with respect to , of in . The iterative procedure stops when and 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 , , on , which is a geodesic between and , where , the above unrolling, unwrapping and wrapping procedures along can be formulated as follows.
The unrolling of onto is the piecewise linear segment in the tangent space joining the following successive points:
The unwrapping at time of , with respect to , into is the tangent vector given by
as long as is not in the cut locus of .
The wrapping at time , along , of a tangent vector back into is the point given by
where is the tangent vector in specified 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 :
a unit speed geodesic starting from takes the form
where such that and ;
the exponential map at has the expression
for any such that ;
the inverse exponential at is given by
for all ;
the parallel transport of the tangent vector along the geodesic defined by (2) is the vector field along given by
where , as and as the covariant derivative of along is .
3 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 () 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 , 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 for . In this section, we summarize the facts required for the exponential maps on the shape space of labelled configurations in and develop the result in Le (2003) further for general to the extent that it is usable for practical purposes.
Recall that, for a configuration in with 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 matrix with . The space of such pre-shapes is known as the pre-shape space of configurations of labelled landmarks in and is the unit sphere of dimension , i.e. . The Kendall shape space of configurations with labelled landmarks in is then the quotient space of by the rotation group acting on the left (cf. Dryden and Mardia (2016) & Kendall et al. (1999)). We shall use to denote the shape of the pre-shape .
Writing for the space of real matrices, then the tangent space to at is
The horizontal subspace, which is the same as the Procrustean tangent space, of , with respect to the quotient map from to , can be expressed as
The horizontal subspace is isometric to the tangent space to at the shape of , and this gives a useful isometric representation of the tangent space of at for the statistical analysis of shapes (cf. Kendall et al., 1999).
Similarly, a unit-speed geodesic in the shape space starting from can be isometrically represented by a so-called horizontal geodesic in of the form
where and . The path is usually referred to as the horizontal lift of . To obtain a representation of a shortest geodesic between two shapes and , we take the two pre-shapes and of and respectively such that
is symmetric and all its eigenvalues are non-negative except possibly for, the smallest one, where . Then, a unit-speed shortest geodesic from to can be isometrically represented by the horizontal geodesic from to :
and is the Riemannian shape distance between and . Note that .
Thus, it follows from the expression (3) for the horizontal geodesic from to that, if there is a unique geodesic between and , the inverse exponential map on the shape space can be isometrically represented by its horizontal lift on given by
3.2 Parallel transport
Turning to parallel transport, it is shown in Le (2003) that, along a horizontal geodesic in , a vector field is horizontal and its projection to is the parallel transport, along the shape geodesic , of the projection of to if and only if 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 matrixin (6), which is answered by the following result.
Let , , be a given horizontal -curve in and be a given horizontal tangent vector in . Assume that the rank of is at least , except for at most finitely many . Then, the vector field along is horizontal and the projection of to is the parallel transport, along the shape curve , of the projection of to if and only if is the solution of
where is skew-symmetric and is the unique solution to
Note that given in Proposition 2 is generally not the parallel transport of along in .
Note first that implies that and is a symmetric matrix.
which is equivalent to
Hence, if (6) holds for some such that , must satisfy
The uniqueness of the solution to (8) follows from the fact that, for a given symmetric non-negative definite matrix of rank at least , where and , and a given skew-symmetric matrix , there is a unique skew-symmetric matrix satisfying the equation . To see the latter, write and . Then, if and only if , which is if and only if for , i.e. if and only if .
The result of Proposition 2 provides a numerical method to approximate the horizontal vector field , along a horizontal geodesic , whose projection to is the parallel transport, along the shape geodesic , of the projection of to in the usual way: with step size ,
where is updated at each step using (8). However, obtained in this way is generally neither tangent to nor horizontal. Hence, this requires finding, at each step, the re-normalized horizontal component of . That is, at each step, we need to project to and then to the horizontal subtangent space of at , followed by normalizing it so that its length remains equal to that of .
If is a pre-shape matrix with rank, the projection of onto the tangent space is given by
To obtain the projection of onto , we note that the orthogonal complementary subspace to in is . Since forms an orthonormal basis for the space of skew-symmetric matrices where and is the square matrix with all its elements zero except for the th element which is one, forms an orthonormal basis for . Thus, the projection of onto is given by
3.3 Shape spline fitting algorithm
We are now in a position to construct an algorithm to calculate the shape-space spline, or simply the shape spline, for a given set of configurations in , observed at times , whose pre-shapes are denoted by , in a similar fashion to that given in Kume, Dryden and Le (2007), as follows. Let and be two given small positive numbers.
Rotate the successive pre-shapes, , such that the resulting is the Procrustes fit of the original one onto . That is, rotate the successive pre-shapes to satisfy the conditions that 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 . Denote the set of the resulting pre-shapes by and let be the initial piecewise horizontal geodesic of such that .
Construct grid points between successive two times such that
which gives for and for where the difference between successive is less than or equal to .
Set the base path to be the piecewise horizontal geodesic passing through .
Using the numerical procedure, including the required necessary projections specified following the proof of Proposition 2, unwrap the data into with respect to the base path which gives
Fit the smoothing spline (1) to and find fitted values at the times of the grid points given in Step 2, giving
Using the numerical procedure, including required necessary projections specified following the proof of Proposition 2, wrap back into the pre-shape sphere with respect to the base path , giving .
Successively rotate in such that the resulting is the Procrustes fit of the original one onto , and obtain the piecewise horizontal geodesic passing through the resulting . Then, becomes the base path in the next iteration.
If , replace by and go to Step 9. Otherwise, stop the iterations.
Successively rotate in such that the resulting is the Procrustes fit of the original one onto . Update by the resulting and go to Step 4.
3.4 Size-and-shape splines
All these results have analogues in the size-and-shape setting. For example, for a configuration in with 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 matrix and the space of the pre-size-and-shapes, , is identical with . For , the tangent space is then also and the horizontal subspace of , with respect to the quotient map from to the Kendall size-and-shape space , is given by
Horizontal geodesics in starting from take the form
where . For two given size-and-shapes and , let and be the pre-size-and-shapes of and respectively such that is symmetric and all its eigenvalues are non-negative except possibly for , the smallest one, where . Then, and a shortest geodesic from to can be represented by the horizontal geodesic that connects and :
Thus, the inverse exponential map on the size-and-shape space , using its horizontal lift on , is isometrically represented by
as long as the geodesic between and is unique.
Moreover, a modification of the argument given in Le (2003) shows that, along a horizontal size-and-shape geodesic in , the vector field is horizontal and its projection to is the parallel translation, along the size-and-shape geodesic , of the projection of onto if and only if satisfies the conditions (5) and
Let , , be a given horizontal -curve in and be a given horizontal tangent vector in . Assume that , except for at most finitely many . Then, the vector field along is horizontal and the projection of to is the parallel transport, along the size-and-shape curve , of the projection of onto if and only if is the solution of
where is skew-symmetric and is the unique solution to (8) with being replaced by .
4 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 data points are measured at times .
Our simulation and real data analyses have been implemented in R. The built-in 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 for a smooth base path and for convergence. The choice of the smoothing parameter , adapted to the data, is determined by applying the (Euclidean) leave-one-out 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 is the unrolling of to the tangent space at its starting point, is the smoothing shape spline with the th observation excluded and is the unwrapped th observed shape data with respect to . The candidates for in the cross-validation were taken from the set . Consequently in the following, the chosen smoothing parameters for the original and perturbed data in the simulation study were and , respectively, and that for the peptide data analysis was . The step size for the approximation in parallel transport was calculated by allowing two equally spaced steps between every two successive shapes.
We consider configurations in with 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 , 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 Figure2.
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