{'fraction_non_alphanumeric': 0.03463219776906844, 'fraction_numerical': 0.01043864938197166, 'mean_word_length': 4.443487179487179, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Polynomial Regression on Riemannian Manifolds

In this paper we develop the theory of parametric polynomial regression in Riemannian manifolds and Lie groups. We show application of Riemannian polynomial regression to shape analysis in Kendall shape space. Results are presented, showing the power of polynomial regression on the classic rat skull growth data of Bookstein as well as the analysis of the shape changes associated with aging of the corpus callosum from the OASIS Alzheimer's study.


Introduction
The study of the relationship between measured data and descriptive variables is known as the field of regression analysis. As with most statistical techniques, regression analysis can be broadly divided into two classes: parametric and non-parametric. The most widely known parametric regression methods are linear and polynomial regression in Euclidean space, wherein a linear or polynomial function is fit in a least-squares fashion to observed data. Such methods are the staple of modern data analysis. However, classical regression suffers from the fundamental limitation that the data must lie in a Euclidian vector space. The most common non-parametric regression approaches are kernel-based methods and spline smoothing approaches which provide much more flexibility in the class of regression functions. However, their non-parametric nature presents a challenge to inference problems, for instance if one wishes to perform a hypothesis test to determine whether the trend for one group of data is significantly different from that of another group.
Fundamental to the analysis of anatomical imaging data within the framework of computational anatomy is the analysis of transformations and shape which are best represented as elements of Remaninan Manifolds rather than Euclidian vector spaces. In previous work, non-parametric kernel-based or spline-based methods have been extended to observations that lie on a Remannian manifold with some success [Davis et al., 2010, Jupp andKent, 1987], but intrinsic parametric regression on Riemannian manifolds has received limited attention. Most recently, Fletcher [Fletcher, 2011] and Niethammer et al. [Niethammer et al., 2011] have each independently developed geodesic regression which generalizes the notion of linear regression to Riemannian manifolds.
The goal of the current work is to extend such work in order to accommodate more flexibility in the model while remaining in the parametric setting. The increased flexibility introduced by the methods in this manuscript allow a better description of the variability in the data, and ultimately will allow more powerful statistical inference. Our work builds off that of Jupp & Kent [Jupp and Kent, 1987], whose method for fitting parametric curves to the sphere involved intricate unwrapping and rolling processes. The work presented in this paper allows one to fit regression curves on a general Riemannian manifold, using intrinsic methods and avoiding the need for unwrapping and unrolling.
We demonstrate the usefulness of our algorithm in three studies of shape. By applying our algorithm to Bookstein's classical rat skull growth dataset [Bookstein, 1991], we show that we are able to obtain a parametric regression curve of similar quality to that produced by non-parametric methods [Kent et al., 2001]. We also demonstrate in a 2D corpus callosum aging study that, in addition to providing more flexibility in the traced path, our polynomial model provides information about the optimal parametrization of the time variable.

Preliminaries
Let (M, g) be a Riemannian manifold [do Carmo, 1992]. For each point p ∈ M , the metric g determines an inner product on the tangent space T p M as well as a way to differentiate vector fields X, Y with respect to one another. That derivative is referred to as the covariant derivative and is denoted ∇ is a smooth vector field on M and γ : [0, T ] → M is a smooth curve on M then the covariant derivative of V along γ is the time derivative of V in a reference frame along γ where the intrinsic derivative D dt is determined by the metric. Geodesics γ : [0, T ] → M are characterized by the second-order covariant differential equation (c.f. [do Carmo, 1992]) This equation, called the geodesic equation, uniquely determines geodesics up to a choice of initial conditions (γ(0),γ(0)) ∈ T M . The mapping from the tangent bundle into the manifold is called the exponential map, Exp : T M → M . Fixing the base point p ∈ M , the exponential map is injective on a zero-centered ball B in T p M of some non-zero (possibly infinite) radius. Thus for a point q in the image of B under Exp p there exists a unique vector v ∈ T p M corresponding to a minimal length path under the exponential map from p to q. The mapping of such points q to their associated tangent vectors v at p is called the log map of q at p, denoted v = Log p q. Given a curve γ : [0, T ] → M we'll want to relate tangent vectors at different points along the curve. These relations are governed infinitesimally by the covariant derivative ∇γ. A vector field X : [0, T ] → T γ(t) M along the curve γ is parallel transported along γ if it satisfies the parallel transport equation: for all times t ∈ [0, T ]. Notice that the geodesic equation is a special case of parallel transport, in which we require that the velocity is parallel transported along the curve itself.

Riemannian Polynomials
Given a vector field X along a curve γ, the covariant derivative of X gives us a way to define vector fields which are "constant" along γ, as parallel transported vectors. Geodesics are generalizations to the Riemannian manifold setting of curves with constant first derivative. The covariant derivative of a vector field along a curve gives another vector field along the curve. We will apply that derivative repeatedly to examine curves with constant higher derivatives. For instance, we refer to the vector field ∇γ (t)γ (t) as the acceleration of the curve γ. Curves with constant acceleration are generalizations of quadratic curves in R and satisfy the second order polynomial equation with initial conditions: γ(0) ,γ(0) andγ(0).
Extending this idea, a cubic polynomial is defined as curve having constant jerk (time derivative of acceleration), and so on. Generally, a kth order polynomial in M is defined as a curve γ : with initial conditions: for all time points t ∈ [0, T ]. As with polynomials in Euclidian space, a k th order polynomial is determined by the (k + 1) initial conditions at t = 0.
The covariant differential equation governing the evolution of Riemannian polynomials is linearized in the same way that a Euclidian ordinary differential equation is. Introducing vector fields v 1 (t), . . . , v k (t) ∈ T γ(t) M , we can write the system of covariant differential equations asγ Algorithm 1 Pseudocode for forward integration of k th order Riemannian polynomial γ ← γ(0) Sample polynomial curves emanating from a common basepoint (green) on the sphere (black=geodesic,blue=quadratic,red=cubic).
The Riemannian polynomial equation cannot, in general, be solved in closed form, and must be integrated numerically. In order to discretize this system of covariant differential equations, we implement the covariant integrator depicted in Alg. 1. At each step of the integrator, each vector is incremented within the tangent space at γ(t) and the results are parallel transported infinitesimally along a geodesic from γ(t) to γ(t + ∆t). The only ingredients necessary to integrate a polynomial forward in time are the exponential map and parallel transport on the manifold. Figure 1 shows the result of integrating polynomials of order one, two, and three. Initial velocity, acceleration, and jerk were chosen and a cubic polynomial integrated to obtain the blue curve. Then the initial jerk was set to zero and the blue quadratic curve was integrated, followed by the black geodesic whose acceleration was set to zero.

Estimation of Parameters for Regression
In order to regress polynomials against observed data y i ∈ M, i = 1, . . . , N , we define the following constrained objective function which is minimized in order to find the optimal initial conditions γ(0), v i (0), i = 1, . . . , k, which we will refer to as the parameters of our model. In order to determine the optimal parameters γ(0), v i (0), i = 1, . . . , k, we introduce Lagrange multiplier vector fields λ i ∈ X(M ), i = 0, . . . , k, often called the adjoint variables, and define the unconstrained objective function As is standard practice, the optimality conditions for this equation are obtained by taking variations with respect to all arguments of E, integrating by parts when necessary. The resulting variations with respect the the adjoint variables yield the original dynamic constraints: the polynomial equations. Variations with respect to the primal variables gives rise to the following system of equations, termed the adjoint equations. The adjoint equations take the following form (see appendix for derivation): where R is the Riemannian curvature tensor and the Dirac δ functional indicates that the order zero adjoint variable takes on jump discontinuities at time points where data is present. Gradients of E with respect to initial and final conditions give rise to the terminal endpoint conditions for the adjoint variables, as well as expressions for the gradients with respect to initial conditions.
In order to determine the value of the adjoint vector fields at t = 0, and thus the gradients of our functional E 0 , the adjoint variables are first initialized to zero at time T , then the system Eq. 24 is evolved backward in time to t = 0.
Given the gradients with respect to the initial conditions, a simple steepest descent algorithm is used to optimize the functional. The update to γ(0) is computed using the exponential map, and the vectors v i (0) are updated via parallel translation.
Note that in the special case of a zero-order polynomial (k = 0), the only gradient λ 0 is simply the mean of the log map vectors at the current estimate of the Fréchet mean. So this method generalizes the common method of Fréchet averaging on manifolds [Fletcher et al., 2004].
The curvature term, in the case k = 1, indicates that λ 1 is a sum of adjoint Jacobi fields. So this approach subsumes geodesic regression as presented by Fletcher [Fletcher, 2011]. For higher order polynomials, the adjoint equations represent a generalized Jacobi system.
In practice, it is often not necessary to explicitly compute the curvature terms. In the case that the manifold M embeds into a Hilbert space, the extrinsic adjoint equations can be computed by taking variations in the ambient space, using standard methods. Such an approach gives rise to the regression algorithm found in Niethammer et al. [Niethammer et al., 2011], for example.
Algorithm 2 Pseudocode for reverse integration of adjoint equations for k th order Riemannian polynomial

Time Reparametrization
In the geodesic model, curves propagate at a constant speed, a result of their extremal action property. However, using higher-order polynomials it is possible to generate curves whose images match those of a geodesic, but whose timedependence has been reparametrized. If the initial conditions consist of all v i collinear, then this will necessarily be the case. Polynomials provide flexibility not only in the class of paths that are possible, but in the time dependence of the curves traversing those paths. Regression models could even be imagined in which the operator wishes to estimate geodesic paths, but is unsure of parametrization, and so enforces the estimated parameters to be collinear.

Polynomials on Riemannian Lie Groups
A special case arises when the manifold is a Lie group and the metric is left-invariant. In this case, left-translation by γ and γ −1 , denoted L γ , L γ −1 , allow one to isometrically represent tangent vectors at γ as tangent vectors at the identity, using the pushforward L γ −1 * . The tangent space at the identity is isomorphically identified with the Lie algebra g. The algebraic operation in g is called the Lie bracket [, ] : g × g → g. For X, Y ∈ g, this bracket operation, also called the adjoint action of X on Y and often denoted ad X Y = [X, Y ], is bilinear and alternating in its arguments (i.e. ad X Y = − ad Y X). Fixing X, the operator ad X is linear and its adjoint is computed with respect to the metric: We will refer to the operator ad † X as the adjoint-transpose action of X. Suppose X is a vector field along the curve γ and X c = L γ −1 * X ∈ g is its representative at the identity. Similarly, ω c = L γ −1 * γ is the representative of the curve's velocity in the Lie algebra. Then the covariant derivative of X along γ evolves in the Lie algebra via In the special case when X =γ so that X c = ω c , if we set this covariant derivative equal to zero we have the geodesic equation:ω which in this form is often refered to as the Euler-Poincare equation [Arnol'd, 1989, Holm et al., 1998]. The curvature tensor is given by the standard formula Applying Eq. 29, in a Lie group with left-invariant metric, these terms are computed in the Lie algebra using the formulas: Also note that if the metric is both left-and right-invariant, then the adjoint-transpose action is alternating and we can simplify the covariant derivative further. In particular, in the presence of such a bi-invariant metric, the first few covariant derivatives take the following forms: This allows us to solve the forward evolution in the convenient setting of the vector space g. Parallel transport requires solution of the equation which can be integrated using Euler integration or similar time-stepping algorithms, so long as the adjoint action can be easily computed. Because of the simplified covariant derivative formula for a bi-invariant metric, the curvature takes the simple form [do Carmo, 1992]:

Coefficient of Determination for Regression in Metric Spaces
It will be useful to define a statistic which will indicate how well our model fits some set of observed data. As in [Fletcher, 2011], we compute the coefficient of determination of our regression polynomial γ(t), denoted R 2 . The first step to computing R 2 is to compute the variance of the data. As discussed in [Fletcher, 2011], the natural choice of total variance statistic is the Fréchet variance, defined by Note that the Fréchet meanȳ itself is the 0-th order polynomial regression against the data {y i } and the variance is the value of the objective function E 0 at that point. We also define the sum of squared error for a curve γ as the value E 0 (γ): Then the coefficient of determination is defined as Any regressed polynomial will have SSE < Var{y i }, R 2 will have a value between 0 and 1, with 1 indicating a perfect fit and 0 indicating that the curve γ provides no better fit than does the Fréchet mean.

The n-Dimensional Sphere
Suppose M = S n = {p ∈ R n+1 : p = 1}. Then geodesics are great circles and the exponential map is given by The corresponding log map is the inverse of this function: The Riemannian curvature tensor is [do Carmo, 1992] Parallel translation of a vector X along the exponential map of a vector v is performed as follows. The vector X is decomposed into a vector parallel to v, which we denote X v , and a vector orthogonal to v, X ⊥ . X ⊥ is left unchanged by parallel transport along Exp p v, while X v transforms in a similar way as v: Using this exponential map and parallel transport, the integrator depicted in Alg. 1 was implemented. Figure 1 shows the result of integrating polynomials of order one, two, and three, using the equations from this section, for the special case of a two-dimensional sphere.

The Lie Group SO(3)
In order to illustrate the algorithm in a Lie group, we consider the Lie group SO(3) of orthogonal 3-by-3 real matrices with determinant one. We review the basics of SO(3) here, but the reader is encouraged to consult [Arnol'd, 1989] for a full treatment of SO(3) and derivation of the Euler-Poincare equation there. The Lie algebra for the rotation group is so (3), the set of all skew-symmetric 3-by-3 real matrices. Such matrices can be identified with vectors in R 3 by the identification under which the Lie bracket takes the convenient form of the cross-product of vectors: With the R 3 representation of vectors in so (3), an invariant metric is determined by the choice of a symmetric, positivedefinite matrix A: With this metric the adjoint-transpose action can be written as If we plug this into the general Lie group formula we obtain the geodesic equation in the Lie algebra: which must be integrated numerically to find ω c at each time. To compute the exponential map Exp p v, the algorithm depicted in Alg. 1 is employed. The matrix v is converted to a skew-symmetric matrix by multiplication by p −1 , which translates v back into the Lie algebra. The resulting vector is converted to a vector ω c (0). The geodesic equation is then integrated to find ω c at future times. The point p then evolves along the exponential map by time-stepping, where at each step ω c is converted back into a skew-symmetric matrix W , multiplied by ∆t and exponentiated to obtain a rotation matrix. The Euler step is then given by left-multiplication by the resulting matrix. Plugging into Eq. 35, the terms need to compute the curvature of SO (3) are and the parallel transport equation is given bẏ The evolution equations are then integrated numerically as in Alg. 1.

Kendall Shape Space
The analysis of shape is a common challenge in medical imaging. Shape is defined as the property of the data that is invariant to scale or rotation. It was in this setting that Kendall [Kendall, 1984] originally developed his theory of shape space. We describe here Kendall's shape space of m-landmark point sets in R d , denoted Σ m d . For a more complete overview of Kendall's shape space, including computation of the curvature tensor, the reader is encouraged to consult the work of Kendall and Le [Kendall, 1984, Kendall, 1989, Le and Kendall, 1993.
..,m , x i ∈ R d be a labelled pointset in R d . Translate the pointset so that the centroid resides at 0, removing any effects due to translations. Another degree of freedom due to scaling is removed by requiring that m i=1 x i 2 = 1. After this standardization, one is left with a point in S md−1 , commonly referred to as preshape space in this context.
Kendall shape space Σ m d is obtained by taking the quotient of the preshape space by the natural action of the rotation group SO(d) on pointsets. In practice, points in the quotient (refered to as shapes) are represented by a member of their equivalence class in preshape space. The work of O'Neill [O'Neill, 1966] characterizes the link in geometry between the shape and preshape spaces. We describe now how to compute exponential maps, log maps, and parallel transport in shape space, using representatives in S md−1 .
First, notice that there exist tangent vectors in S md−1 corresponding to global rotations of the preshape. At a particular point p in preshape space, these vectors form a linear subspace of the tangent space T p S md−1 , which we define as the vertical subspace. Curves moving along vertical tangent vectors result in rotations of a preshape, and so do not indicate any change in actual shape. Since we are concerned with shape change and wish to ignore effects due to rotation, it will be important to project such tangent vectors to be purely non-vertical, or horizontal. We denote the horizontal projection by H : T S md−1 → T S md−1 .
A vertical vector in preshape space arises as the derivative of a rotation of a preshape. The derivative of such a rotation is a skew-symmetric matrix W , and its action on a preshape x has the form (W x 1 , . . . , W x n ) ∈ T S md−1 . The vertical subspace is then spanned by such tangent vectors arising from any linearly independent set of skew-symmetric matrices. The projection H is performed by taking such a spanning set, performing Gram-Schmidt orthonormalization, and removing each component.
The horizontal projection allows us to relate the covariant derivative on the sphere to that on shape space. Lemma 1 of O'Neill [O'Neill, 1966] states that if X, Y are horizontal vector fields at some point p in preshape space, then where ∇ denotes the covariant derivative on preshape space and ∇ * , X * , and Y * are their counterparts in shape space. For a general shape space Σ m d , the exponential map and parallel translation are performed using representatives in S md−1 . In practice, this usually must be done in a time-stepping algorithm, in which at each time step an infinitesimal spherical parallel transport is performed, followed by the horizontal projection. The resulting algorithm can be used to compute the exponential map as well. Computation of the log map is less trivial, as it requires an iterative optimization routine. The target shape is first aligned to the base shape using Procrustes alignment. Then, at each iteration of the log map routine, the exponential map is integrated forward, compared to the target preshape, then the resulting tangent vector is adjoint parallel transported back to the base point, for update in a steepest descent iteration.
A special case arises in the case when d = 2. In this case the exponential map, parallel transport and log map have closed form. The reader is encouraged to consult [Fletcher, 2011] for more details about the two-dimensional case.
With exponential map, log map, and parallel transport implemented as described here, one can perform polynomial regression on Kendall shape space in any dimension. We demonstrate such regression now on examples in two dimensions.

Rat Calivarium Growth
The first dataset we consider was first analyzed by Bookstein [Bookstein, 1991]. The data is available for download at http://life.bio.sunysb.edu/morph/data/datasets.html and consists of m = 8 landmarks on a midsagittal section of rat calivaria (upper skulls). Landmark positions are available for 18 rats and at 8 ages apiece. Riemannian polynomials of orders k = 0, 1, 2, 3 were fit to this data. The resulting curves are shown in Fig. 2 and in zoomed detail in Fig. 3. Clearly the quadratic and cubic curves differ from that of the geodesic regression. The R 2 values agree with this qualitative difference: the geodesic regression has R 2 = 0.79, while the quadratic and cubic regressions both have R 2 values of 0.85 and 0.87, respectively. While this shows that there is a clear improvement in the fit due to increasing k from one to two, it also shows that little is gained by further increasing the order of the polynomial. Qualitatively, Fig. 3 shows that the slight increase in R 2 obtained by moving from a quadratic to cubic model corresponds to a marked difference in the curves, possibly indicating that the cubic curve is overfitting the data.

Corpus Callosum Aging
The corpus callosum is the major white matter bundle connecting the two hemispheres of the brain. In order to investigate shape change of the corpus callosum during normal aging, polynomial regression was performed on a collection of data from the OASIS brain database (www.oasis-brains.org). Magnetic resonance imaging (MRI) scans from 32 normal subjects with ages between 19 and 90 years were obtained from the database. A midsagittal slice was extracted from each volumetric image and segmented using the ITK-SNAP program (www.itksnap.org). Sets of 64 landmarks were optimized on the contours of these segmentations using the ShapeWorks program [Cates et al., 2007] (www.sci.utah.edu/software.html). The algorithm generates samplings of each shape boundary with optimal correspondences among the populaton. Regression results for geodesic, quadratic, and cubic regression are shown in Fig. 4. At first glance the results appear similar for the three different models, since the motion envelopes show close agreement. However, the R 2 values show an improvement from geodesic to quadratic (from 0.12 to 0.13) and from quadratic to cubic (from 0.13 to 0.21). Inspection of the estimated initial conditions, shown in Fig. 5 reveals that the tangent vectors appear to be rather collinear. For the reasons stated in Sec. 2.4, this suggests that the differences can essentially be described as reparametrization of the time variable, which is only accommodated in the higher order polynomial models.

Appendix: Derivation of Adjoint Equations
In this section we derive the adjoint system for the polynomial regression problem. The approach to calculus of variations followed here is outlined by, for example, Noakes et al. [Noakes et al., 1989]. Consider a simplified objective function containing only a single data term, at time T : The first term is integrated by parts to yield The variation with respect to v i for i = 1, . . . , k is then given by In order to determine the differential equation for λ 0 , the variation with respect to γ must be computed. Let W again denote a test vector field along γ. For some > 0, let {γ s : s ∈ (− , )} be a differentiable family of curves satisfying If is chosen small enough, the vector field W can be extended to a neighborhood of γ such that [W,γ s ] = 0, where a dot indicates the derivative in the ∂ ∂t direction. The vanishing Lie bracket implies the following identities Finally, the vector fields v i , λ i are extended along γ s via parallel translation, so that The variation of E with respect to γ satisfies As the λ i are extended via parallel translation, their inner products satisfy Then applying this to each term in the previous equation, Then by construction, since ∇ W v i = 0, Then using the Lie bracket and curvature identities, this is written as which is further simplified, again using the identity ∇ W v i = 0: Using the Bianchi identities, it can be demonstrated that the curvature tensor satisfies the identity [do Carmo, 1992]: A, R(B, C)D = − B, R(D, A)C , for any vectors A, B, C, D. The covariant derivative along γ is also integrated by parts to arrive at T 0 δ γ E, W dt = − Log γ(T ) y, W (T ) Finally, gathering terms, the adjoint equation for λ 0 and its gradients are obtained: δ γ(T ) E = 0 = − Log γ(T ) y + λ 0 (86) Along with the variations with respect to v i , this constitutes the full adjoint system. Extension to the case of multiple data at multiple time points is trivial, and results in the adjoint system presented in Sec. 2.3.