Generalized partially linear models on Riemannian manifolds

The generalized partially linear models on Riemannian manifolds are introduced. These models, like ordinary generalized linear models, are a generalization of partially linear models on Riemannian manifolds that allow for response variables with error distribution models other than a normal distribution. Partially linear models are particularly useful when some of the covariates of the model are elements of a Riemannian manifold, because the curvature of these spaces makes it difficult to define parametric models. The model was developed to address an interesting application, the prediction of children's garment fit based on 3D scanning of their body. For this reason, we focus on logistic and ordinal models and on the important and difficult case where the Riemannian manifold is the three-dimensional case of Kendall's shape space. An experimental study with a well-known 3D database is carried out to check the goodness of the procedure. Finally it is applied to a 3D database obtained from an anthropometric survey of the Spanish child population. A comparative study with related techniques is carried out.


Introduction
Classification problems arise in many real-life situations.A new observation has to be classified on the basis of a training set, which is described by a set of features whose class memberships are known.Supervised learning techniques have been widely studied when the features lie on a vector space (Hastie et al., 2009).When features do not form a vector space, well-known supervised learning techniques are not well suited to the classifiers Tuzel et al. (2008).However, features can also take values on a Riemannian manifold.This is common in fields such as astronomy, geology, meteorology, etc., which include natural distributions on spheres, tangent bundles and Lie groups (González-Manteiga et al., 2012).Another example, this time in the field of computer vision, would be the space of non-singular covariance matrices (Tuzel et al., 2008).One discipline that certainly offers many examples in different fields of applications (biology, medicine, chemistry, etc.) is statistical shape analysis (Dryden and Mardia, 2016).Many problems involve predicting a categorical variable as a function of the shape of an object that lies in a Riemannian manifold.
Although different approaches can be identified in shape analysis based on how the object is treated in mathematical terms (Stoyan and Stoyan, 1995), the majority of research has been restricted to landmark-based analysis, where objects are represented using k labeled points in the Euclidean space R m .These landmarks are required to appear in each data object, and to correspond to each other in a physical sense.Seminal papers on this topic are Bookstein (1978), Kendall (1984), and Goodall (1991).The main references are Dryden and Mardia (2016) and Kendall et al. (2009).In this paper we concentrate on this approach.
In a formal way, shape can be defined as the geometrical information about the object that is invariant under a Euclidean similarity transformation, i. e., location, orientation, and scale.The shape space is the resulting quotient space.When the landmark-based approach is used, the corresponding shape space is a finite-dimensional Riemannian manifold, and statistical methodologies on manifolds must be used.There are several difficulties in generalizing probability distributions and statistical procedures to measurements in a non-vectorial space like a Riemannian manifold, but fortunately, there has been a significant amount of research and activity in this area over recent years.An excellent review can be found in Pennec (2006).
The most immediate approach for solving the classification problem when predictive variables take values on a Riemmanian manifold would be to map the manifold to a Euclidean space, i.e. to flatten the manifold.But, in a general case, mapping that globally preserves the distances between the points on the manifold is not available.As a consequence, the flattened space would not represent the global structure of the points appropriately.Although statistical analysis of manifold-valued data has gained a great deal of attention in recent years, there is little literature on classification.In fact, to our knowledge the only reference is Tuzel et al. (2008), and it is restricted to a binary classification problem.A LogitBoost (Friedman et al., 2000) on Riemannian manifolds is proposed in Tuzel et al. (2008).It is similar to the original LogitBoost, except for differences at the level of weak learners (the regression functions are learned on the tangent space at the weighted mean of the points).
Another related work is González-Manteiga et al. ( 2012), but they studied a regression problem (the predicted variable is real valued) rather than a classification problem.They introduced partially linear models on Riemannian manifolds (robust estimators can be found in Henry and Rodríguez (2014)).Partially linear models were proposed by Engle et al. (1986).Since then, partially linear models have been used in the study of complex nonlinear problems, some recent examples are Zhang et al. (2017), Qian and Wang (2017), Cui et al. (2017) and Hilafu and Wu (2017).In this semiparametric regression method, the dependent variable is modeled with a parametric linear part and a nonparametric part.In González-Manteiga et al. (2012) the variable to be non-parametrically modeled is in a manifold, and they proposed a kernel-type estimator.
Based on this idea, we can generalize partially linear models (GPLM) to solve the classification problem with features in a Riemmanian manifold.We benefit from the flexibility of partially linear models and at the same time we can include features in a Riemmanian manifold.To our knowledge, this is the first time GPLMs have been defined on Riemmanian manifolds.At the same time, we also propose a solution for the classification problem for more than two classes, the only case studied to date.In particular, we also introduce a solution for when the dependent variable is ordinal.Furthermore, unlike the method proposed in Tuzel et al. (2008) where features in a Riemmanian manifold were the only predictive variables, other predictive variables together with those in a Riemmanian manifold are managed jointly by our proposal.
This paper addresses an important current application: size fitting for online garment shops, in particular children's garment size matching.Customers face a challenge when they have to choose the right size of garment without try it on when buying these items both in store and, especially, in online clothes shops (Ding et al., 2011).Although users can base their decision on their previous experience (their virtual closet), children are constantly growing, so this not suitable strategy (Sindicich and Black, 2011).Not only that, but each company also has its own sizing, and what is more, this can change over time (Schofield and LaBat, 2005).As a consequence, size matching in children should be based on their current form.
There is usually a sizing chart that corresponds to several anthropometric measurements, together with their ranges to show the size assignation.Nevertheless, customer's measurements can lead to different size assignations depending on which measurements are considered.Therefore, customers cannot know which size will fit them best (Labat, 2007).As a result, size fitting problems lead to a high percentage of returns, which represents one of the main costs of this sales channel for distributors and manufacturers.The return rates of some e-commerce businesses are between 20 and 50% (Eneh, 2015).This also decreases customer satisfaction (Otieno et al., 2005) and the likelihood that the customer will buy again.Moreover, concern about poor fit is the main obstacle to purchasing clothes online.
To address the child garment size matching problem, a fit assessment study was carried out by the Biomechanics Institute of Valencia.A sample of Spanish children aged between 3 and 12 years were scanned using the Vitus Smart 3D body scanner from Human Solutions.This has a non-intrusive laser system consisting of four columns that house the optic system, which moves from head to feet in ten seconds, performing a sweep of the body.The body shape of each child in our data set was represented by 3075 3D landmarks.Although a 3D body scanner is not usually available for customers, nowadays customers can obtain their detailed body shapes using their own digital cameras or other measuring technologies (Cordier et al., 2003;Ballester et al., 2015).Recently, 3D bodies have been reconstructed from images captured with a smartphone or tablet in Ballester et al. (2016).Furthermore, a subsample of these children tested different garments of different sizes, and their fit was assessed by an expert.This expert labeled the fit as 2 (correct), 1 (if the garment was small for the child) or 3 (if the garment was large for the child) in an ordered factor called Size-fit.
Therefore, finding the garment size that best fits the user is a statistical classification problem (Meunier, 2000).In this problem, the children's body shapes represented by landmarks are predictive variables in a Riemmanian manifold.The proposed method has been applied to the aforementioned database of children with excellent results.
To our knowledge, the only previous reference about the child garment size matching problem is Pierola et al. (2016).However, they used multivariate features, not the complete information about the child's form.In particular, they used the differences between the reference mannequin of the evaluated size and the child for several anthropometric measurements.If the reference mannequin is not available, that methodology cannot be used.
As regards other works that also use variables in a Riemmanian manifold in the context of the apparel industry, in Vinué et al. (2016) women's body shapes represented by landmarks were used to define a new sizing system by adapting clustering algorithms to the shape space.Unlike our supervised learning problem, they dealt with an unsupervised learning problem.Another unsupervised learning problem is faced in Epifanio et al. (2017), where archetypal shapes of children are discovered.
The R language (R Core Team, 2017) was employed in our implementations.We used the shapes package by Ian Dryden (Dryden, 2017).This is a very powerful and complete package for the statistical analysis of shapes.
The article is organized as follows: Section 2 reviews the partially linear models on Riemmanian manifolds, which are generalized in Section 3 to Riemmanian manifolds.Algorithms for their estimation are also given.Their R (R Core Team (2017)) code is available at www3.uji.es/~epifanio/RESEARCH/partly.rar.Section 4 describes the basic concepts of statistical shape analysis, and explains how to estimate generalized partially linear models on the Kendall's 3D Shape Space.The use of the logistic partly linear model on the Kendall's 3D Shape Space is illustrated by a well-known data set in Section 5, while the ordered partially linear model is applied to solve the child garment size matching problem in Section 6.Finally, conclusions are discussed in Section 7.

Partially linear models on Riemannian manifolds
Partially linear models (PLM) Engle et al. (1986) are regression models in which the response depends on some covariates linearly but on other covariates nonparametrically.PLMs generalize standard linear regression techniques and are special cases of additive models (Hastie and Tibshirani, 1990;Stone, 1985), which makes it easier to interpret the effect of each variable.
Partially linear models when one of the predictive variables takes values on a Riemannian manifold were introduced in González- Manteiga et al. (2012).In this work, they consider a sample {(y i , , where the response variable, Y , is a real valued scalar variable, x t i is a real valued p-dimensional vector and s i is a point of a Riemannian manifold, M , of dimension d.They assume the partially linear model: and and with independent errors i and η ij .Therefore, β, φ 0 (s) and φ(s) are the parameters to estimate.Manteiga et al.González-Manteiga et al. (2012) suggest estimating φ 0 (s) and φ(s) using non-parametric kernel-type estimators on Riemannian manifolds (see section 2.1) and then estimating the parameter β considering the leastsquares estimator obtained by minimizing:

Non-parametric estimators on Riemannian manifolds
Let {(x 1 , s 1 ), . . ., (x n , s n )} be iid random vectors that take values on R × M .Due to the curvature of M , kernel-type estimators of φ(s) = E(x | s) must be adapted to this space.
Pelletier Pelletier (2006) proposes the following non parametric estimator: where θ s (s i ) is the volume density function of M ; ρ is the Riemannian distance on M and K hn is a univariate kernel function with bandwidth h n with lim n→∞ h n = 0 and h n < i M , i M being the injectivity radius of M .In Pelletier (2006) we can find some good properties of this estimator.
3 Generalized Partially Linear Model on Riemannian manifolds Although our proposal can be extended to generalized linear models in general, we will focus on two particular important models that we will use in our applications: logistic and ordered logistic models.

Logistic Partially Linear Model on Riemannian manifolds
Let {(y 1 , x 1 , s 1 ), . . ., (y n , x n , s n )} be a set, where y i are binary variables, x i real valued p-dimensional vectors and s i are points in M , a Riemannian manifold of dimension d.
we can assume the logistic partially linear model: and with g(s) = φ 0 (s) − φ t (s)β; φ(s) = (φ 1 (s), . . ., φ p (s)); and where β, φ 0 (s) and φ(s) are the parameters to estimate.As in González-Manteiga et al. ( 2012), because s is in a Riemannian manifold, the estimation of φ j (s) j = 1, ..., p must be obtained using equation 3.In the next section the expression of this estimator will be given for the particular and difficult case of Kendall's 3D shape space.
The algorithm that we propose follows the ideas of additive generalized linear models (Hastie and Tibshirani, 1990;Hastie et al., 2009): a partially linear model is applied to the adjusted dependent variable at each step of the iteratively reweighted least squares (IRLS) algorithm.It is as follows (the superindex (j) indicates the estimation in the j-th iteration): While (e (j) > specified threshold) and (j < specified maximum number of steps) do Calculate g For i = 1, . . ., n Construct the working target variable

Ordered Partially Linear Model on Riemannian manifolds
The above algorithm can be modified to model an ordinal response, in particular we will assume the cumulative logistic model or proportional odds model McCullagh (1980); Agresti (2010).Let {(y 1 , x 1 , s 1 ), . . ., (y n , x n , s n )} be a set, where y i are response variables, x i real valued p-dimensional vectors and s i are points in M , a Riemannian manifold of dimension d.Suppose that the response variable y has K ordered categories and and Following Walker and Duncan (1967); McCullagh (1980) and Thompson and Baker (1981), we treat the cumulative link model as a multivariate generalized linear model Fahrmeir and Tutz (2013) defining and otherwise as zero.In the multivariate case one merely has to substitute vectors and matrices for the multivariate versions.
We define the total design matrix xi =    x t i . . .
Let D i be the derivative of the link function and the weight matrix which can be considered an approximation of the inverse of the covariance matrix of the transformed response).
i ) Construct the working target variable Apply partly linear model to the targets z = (z i ) i=1,...,n with weight matrix ))

End while
In the particular case of the ordinal model with three categories of our application: 4 Kendall's 3D Shape Space In the previous section the logistic and ordered logistic partially linear models were given for a general Riemannian manifold.In this section we give the expressions that we need in order to apply them in the particular and important case of the Kendall 's 3D Shape Space.This manifold has a complicated structure and the calculus of the expressions needed in equation 3 is not trivial.We begin by introducing some basic concepts, a complete introduction to which can be found in Dryden and Mardia (2016).
In our approach to shape analysis each object is identified by a set of landmarks, i.e. a set of points in the Euclidean space R m that identifies each object and match between and within populations.
Definition 1.A configuration matrix X is a k × m matrix with the Cartesian coordinates of the k landmarks of an object.
The shape of an object is all the geometric information that remains invariant with translations, rotations and changes of scale.Thus: Definition 2. The shape space Σ k m is the set of equivalence classes T X of k × m configuration matrices X ∈ R k×m under the action of Euclidean similarity transformations.
As mentioned above, the shape space Σ k m admits a Riemannian manifold structure.The complexity of this Riemannian structure depends on k and m.For example, Σ k 2 is the well-known complex projective space.For m > 2, which is the case of our application, they are not familiar spaces and may have singularities.
A representative of each equivalence class T X can be obtained by removing the similarity transformations one at a time.There are different ways to do that.
Let X be a configuration matrix.One way to remove the location effect consists of multiplying it by the Helmert submatrix, H, i. e., X H = HX.
To filter scale, we can divide X H by the centroid size, which is given by • is the Frobenius norm.So, is called the pre-shape of the configuration matrix X because all information about location and scale is removed, but rotation information remains.
Definition 3. The pre-shape space S k m is the set of all possible pre-shapes.
S k m is a hypersphere of unit radius in R m(k−1) (a Riemannian manifold that is widely studied and known).Σ k m is the quotient space of S k m under rotations.As a result, a shape S X is an orbit associated with the action of the rotation group SO(m) on the pre-shape.
From now on, in order to simplify the notation, we will use S X to denote both, a configuration matrix and its shape, provided that it is understood by context.
For m = 2, this quotient space is isometric with the complex projective space CP k−2 , a familiar Riemannian manifold without singularities.For m > 2, Σ k m is not a familiar space, and it has singularities.The singularities are shapes whose preshapes have rank m − 2 or less.With real world applications we can usually assume that our data are almost certainly in the non-singular part of the shape space and, fortunately, the Riemannian structure of the non-singular part of Σ k m can be obtained taking into account that π : S k m → Σ k m is a Riemannian submersion (Kendall et al., 2009), for any π(Z) ∈ Σ k m the tangent space T π(Z) Σ k m can be identified with the horizontal space H Z of T Z S k m .

Riemannian distance
The induced Riemannian distance in the shape space is given by the Procrustes distance defined as follows.
Definition 4. Given two configuration matrices X 1 , X 2 , the Procrustes distance of its corresponding shapes, ρ(S X1 , S X2 ), is the closest great-circle distance between Z 1 and Z 2 on the pre-shape hypersphere S k m , where Z j = HXj HXj , j = 1, 2. The minimization is carried out over rotations.
The solution for this optimization problem is: where λ 1 ≥ λ 2 ≥ . . .λ m−1 ≥| λ m | are the square roots of the eigenvalues of Z T 1 Z 2 Z T 2 Z 1 , and the smallest value λ m is the negative square root if and only if det(Z T 1 Z 2 ) < 0 Dryden and Mardia (2016).Note that the range of this distance is [0, π/2].

Volume density function
The volume density function can be obtained taking into account that the mapping that assigns the corresponding element on the shape space to each preshape Y : is a Riemannian submersion.Then the volume density function is (see A) We must stress here the importance of the volume density function in the case of a large number of landmarks k, because in the limit k → ∞, the definition formula (10) becomes

Generalized Partially Linear Models on the Kendall's 3D Shape Space
Once the necessary concepts have been introduced, we turn to the algorithm to fit a generalized partially linear model on Kendall's Shape Space.We focus on the particular case of the ordered partially linear model (for the logistic partially linear model, it is analogous but instead we apply the algorithm 1).
Algorithm 3. Given a sample {(y 1 , X 1 , x 1 ), . . ., (y n , X n , x n )}, where y i is a realization of an ordered variable with K categories, X i are configuration matrices and x i real valued p-dimensional vectors.

Application to anatomical data
In an investigation into sex differences in the crania of a species of macaque, random samples of 9 male and 9 female skulls were obtained by Paul OHiggins (Hull-York Medical School) (Dryden andMardia, 2016, 1993).A subset of seven anatomical landmarks was located on each cranium and the three-dimensional (3D) coordinates of each point were recorded.The aim of the study was to assess whether there were any size and shape differences between sexes.The data are available in the R shapes package (Dryden, 2017).Fig. 1 (a) shows the distribution of all the landmarks of the 18 subjects.From the configuration matrices {X i } i=1,••• ,18 ∈ M 7×3 , with the coordinates of the landmarks of the 18 macaques, the full Procrustes mean shapes are computed for males and females separately (see Fig. 1 (b)), and their preshapes, {s i } i=1,...,18 , (eq. 9) and sizes {x i } i=1,...,18 (eq.8), are computed (see Fig. 2 and Table 1).If we define Y i = 1 if the i − th cranium belongs to a female and Y i = 0 if it belongs to a male, then we can model: and algorithm 1 can be used to fit a logistic partially linear model.Table 1: Sizes of the 18 crania, computed from the landmark configurations.m1, m2,..., m9 denote the 9 males and f1, f2,..., f9 the 9 females.
In the smoothing procedure, we considered a Gaussian kernel and the bandwidth parameter h was fixed as h = π/100 by using a leave-one-out crossvalidation (CV) procedure.With this value, a 0% CV error was obtained.With threshold =0.0002, the algorithm stops at 7 iterations.(see Table 2).
The estimation procedure provides a β1 = −6.02,signaling that the probability of being female decreases as skull size increases, and the values presented in Table 3 were obtained for the nonparametric part of the model.

Application to children's body shapes
The aim of this section is to show how the aforementioned algorithm can be used to predict the goodness of fit of a given garment size, i.e. small (Y i = 1), good fit (Y i = 2) or large (Y i = 3), as a function of the garment size, the size of h π/100 π/50 π/25 π/10 CV (% of correct classifications) 100% 88.89% 88.89% 88.89% the child and his/her shape.
There are multiple ways to choose the most suitable size in a potential online sales application, and all of them depend on the manufacturer.
A randomly selected sample of 739 Spanish children aged 3 to 12 was scanned using a Vitus Smart 3D body scanner from Human Solutions.The children were asked to wear a standard white garment in order to standardize the measurements.Several cameras capture images and associated software provided by the scanner manufacturers detects the brightest points and uses them to create a triangulation that provides information about the 3D spatial location of 3075 points on the body's surface.The 3D scan data are processed to create of posture-harmonized homologous models to obtain a database of individual 3D homologous avatars with one-toone anatomical vertex correspondence between them.As a result, each child's body shape is represented by 3075 3D landmarks.Because the children's head, hands, legs and feet are not involved in the shirt size selection, these parts were discarded from the scans, and a total of 1423 3D landmarks per child were considered, i.e. each child's the body shape was represented by a 1423 × 3 configuration matrix.Two of them are shown in Figure 3.
Seventy eight of these children performed an additional fit test.All of them tried on the same shirt model in different sizes: the supposedly correct size, the size above and the size below.Then, an expert in clothing and design qualitatively evaluated the fit in each case (as small, correct fit or large).Due to lack of cooperation by some of the children, not all the children tried on all the three sizes, but only two sizes or even one.In 24 cases, only two sizes were evaluated, and 9 children tested just one shirt size.There were 7 shirt sizes available, supposedly corresponding to ages 3, 4, 5, 6, 8, 10 and 12.
Two subsamples are considered, the sample consisting of the 37 boys and that consisting of the 41 girls in the data base.
Algorithm 2 is applied to fit the model:    x ij = φ j (s i ) + η ij , i = 1, . . ., n, j = 1, 2, to each subsample, s i being the body shape of the i-th child, x i1 his/her body size and x i2 the size evaluated.We consider a Gaussian kernel again and in order to choose the value of the bandwidth parameter h and, at the same time, perform a quantitative analysis of the effectiveness of the method, a leave-one-out cross-validation study is conducted.At each step of this study, a child is left out, and his/her fit predicted for the supposedly correct size, the size above and the size below.In Table 4 we can see the percentage of correct classifications in each case.

Comparison with other methods
In Pierola et al. (2016), the authors used ordered logistic regression and random forest methodologies to predict a garment's goodness of fit from the differences between the measurements of the reference mannequin for the evaluated size and the child's anthropometric measurements.We could also have used different children's anthropometric measurements to fit a classical proportional odds model McCullagh (1980).So given the response variable Y with K = 3 ordered categories, and given X a vector of explicative variables formed by the garment size to evaluate and the 27 children's anthropometric measurements considered by Pierola et al. (2016), we could have fitted: Performing a leave-one-out cross-validation study, choosing the model on each step by a forward stepwise model selection based on likelihood ratio tests Christensen (2015), we obtain worse results than those obtained with our methodology.The percentages of correct classifications are now 61.76% for boys and 66.19% for girls.
On the other hand, as stated in the introduction, the shape space and sizeand-shape space are not flat Euclidean spaces, so classical statistical methods cannot be directly applied to the manifold valued data.However, if the sample has little variability, the problem can be transferred to a tangent space (at the Procrustes mean of these shapes or size-and-shapes, for example) and then standard multivariate procedures can be performed in this space Dryden and Mardia (2016), such as Principal Component Analysis (PCA).With this approach, in order to reduce the dimensionality of the data set, the first p PC scores, which summarize most of the variability in the tangent plane data, are usually chosen.
The tangent space is defined from a point called pole, so the distance from the shape to the pole is preserved, i.e. the distance from a point in the manifold to the pole is equal to the Euclidean distance between its projections in the tangent space.As one moves away from the pole, the Euclidean distances between some pairs of points in the tangent space are smaller than their corresponding shape distances.This distortion becomes larger as one considers points further from it.For this reason, the pole should be taken close to all of the points and the mean of the observed shapes is the best choice (Dryden and Mardia, 2016).
So, given the configuration matrices X i ∈ M 1423×3 , the size s i of each child is obtained and the full Procrustes mean shapes are computed for boys and girls separately.Then, the coordinates of the projection of X i ∈ M 1423×3 onto the tangent plane defined at their corresponding mean shape are obtained.The first PC scores of these coordinates are calculated and they will be used as covariates in our predictive model.The first PC components that explain 98% of variability are considered.
So, given the response variable Y with K = 3 ordered categories and given a vector X with the garment size to evaluate, the child's size and the first PC scores of his/her coordinates in his/her corresponding tangent space, we can fit the model given by Eq. 11.
Once again, performing a leave-one-out cross-validation study using this model, we obtain worse results than those obtained with our methodology.The percentages of correct classifications are now 61.43% for boys and 67.12% for girls.

Conclusions
We define GPLMs on Riemmanian manifolds for the first time.Due the application that we address, our GPLMs have focused on Kendall's 3D Shape Space.Although it is an important and common problem in real applications, this problem has not been addressed until now, to the best of our knowledge.We have developed and illustrated the algorithms for estimating the GPLM in two different applications.We have also compared the results with other simpler approximations in the case of the children's garment size matching problem.
Although we have focused on children's shapes in the application, the methodology can also be used to select the right size for adults, men and women.Furthermore, as pointed out in Sect. 1, the proposed methodologies have great potential in all the fields where statistical shape analysis is used, including biological and medical applications.
Besides opening the door to applications in different fields, other future work could focus on other Riemannian manifolds.Moreover, all the work carried out on GPLMs for multivariate data could be extended to the case where variables also take values on Riemannian manifolds.space at s 1 , let B r (s 1 ) be an open geodesic ball of radius r centered at s 1 , let B r (0 s1 ) be an open ball of radius r centered at 0 s1 ∈ T s1 M , and let inj(s 1 ) be the injectivity radius at s 1 , (i.e. the maximum radius r such that the exponential map exp : B r (0 s1 ) ⊂ T s1 M → M is a diffeomorphism).The volume density function θ s1 : B inj(p) (p) → R + is a function defined for any point s 2 of the normal ball B inj(s1) (s 1 ) by where g = exp * s1 (g) is the pullback of g by the exponential map, g is the canonical metric induced by g in B r (0 s1 ), and (U , ψ) is any chart of B r (0 s1 ) that contains w = exp −1 (s 2 ).
In this work, following Kendall et al. (2009), we identify a point Z in the pre-shape sphere S k m as a matrix.The explicit formula of the volume density function in the shape space Σ k m is given in the following theorem.where here ρ(π(Z 1 ), π(Z 2 ))) is the distance in Σ m k from π(Z 1 ) to π(Z 2 ).
Proof.Since π : S k m → Σ k m is a Riemannian submersion, for any π(p) ∈ Σ k m , the tangent space T π(p) Σ k m can be identified with the horizontal space H p of T p S k m .Moreover since dπ : ker(dπ) ⊥ → T S k m is an isometry, for any v, w ∈ H p , it is the case that g (v, w) = g Σ k m (dπ(v), dπ(w)).On the other hand, if q = exp p (w), then d exp p (w) : H p → H q .Therefore, in order to compute θ π(p) we can make use of the restriction of the exponential map exp p | Hp to the horizontal space H p .
In our particular setting, given Z ∈ S k m , the tangent space is the horizontal space is and the exponential map is given by exp We are now going to obtain θ π(Z1) (π(Z 2 )) using the properties of the exponential map in S k m .Given Z 1 ∈ S k m , Z 2 = exp Z1 (W ) for some W ∈ H Z1 , suppose π(Z 1 ) = π(Z 2 ) and suppose moreover that π(Z 2 ) is in a normal ball of π(Z 1 ).Let us now choose an orthonormal basis {V i } d i=1 of H Z1 with V 1 = W W (and with d = dim(H p )). Then Consequently, Hence, since π(Z 2 ) is in a normal ball of π(Z 1 ), there is a minimal and horizontal geodesic segment starting at Z 1 and ending in Z 2 with initial velocity W such that ρ(π(Z 1 ), π(Z 2 )) = W . Taking into account that since V i , W = V i , V j = 0 for any i = j (i, j > 1) and that V i , Z = W, Z 1 = 0 because W and V i are tangent vectors to T Z1 S m k we conclude that θ π(Z1) (π(Z 2 )) = det d exp Z1 (W )(V i ), d exp Z1 (W )(V j ) = sin ρ(π(Z 1 ), π(Z 2 )) ρ(π(Z 1 ), π(Z 2 ))) where d is the dimension of H Z1 .Now, we are going to compute the dimension of H Z for any Z ∈ S k m .Since the tangent space T Z S k m can be decomposed as where dim(V Z ) is the dimension of the fiber π −1 (π(Z)).The dimension of the fiber π −1 (π(Z)) depends on the rank of Z (see Kendall et al. (2009)) but if rank(Z) ≥ m − 1 (i.e. it is a non singular point), π −1 (π(Z)) is homeomorphic to SO(m) (and hence with dimension m(m−1)

Figure 1 :
Figure 1: (a) Landmarks corresponding to the 18 configurations.(b) Landmarks corresponding to the mean shapes of males and females.The different colors represent the different landmarks; the symbol o is used for the landmarks of males and the symbol * is used for females.

Figure 3 :
Figure 3: 3D landmarks of (a) a girl and (b) a boy in the data set.

Table 2 :
Results of the CV analysis to choose the value of the bandwidth parameter for the crania problem.

Table 3 :
Estimation of the nonparametric part of the logistic PLM for each macaque in the data set.

Table 4 :
Results of the CV analysis for the children's body shape problem.