Principal arc analysis on direct product manifolds

Introduction

Principal Component Analysis (PCA) has been frequently used as a method of dimension reduction and data visualization for high-dimensional data. For data that naturally lie in a curved manifold, application of PCA is not straightforward since the sample space is not linear. Nevertheless, the need for PCA-like methods is growing as more manifold data sets are encountered and as the dimensions of the manifolds increase.

Our approach to a manifold version of PCA builds upon earlier work, especially the principal geodesic analysis proposed by Fletcher et al. 2004 and the geodesic PCA proposed by Huckemann and Ziezold 2006 and Huckemann, Hotz and Munk 2010. A detailed catalogue of current methodologies can be found in Huckemann, Hotz and Munk 2010. An important approach among these is to approximate the manifold by a linear space. Fletcher et al. 2004 take the tangent space of the manifold at the geodesic mean as the linear space, and work with appropriate mappings between the manifold and the tangent space. This results in finding the best fitting geodesics among those passing through the geodesic mean. This was improved in an important way by Huckemann, Hotz and Munk, who found the best fit over the set of all geodesics. Huckemann, Hotz and Munk went on to propose a new notion of center point, the PCmean, which is an intersection of the first two principal geodesics. This approach gives significant advantages, especially when the curvature of the manifold makes the geodesic mean inadequate, an example of which is depicted in Figure 2b.

Our method inherits advantages of these methods and improves further by effectively capturing more complex nongeodesic modes of variation. Note that the curvature of direct product manifolds is mainly due to the spherical part, which motivates careful investigation of $S^{2}$-valued variables. We point out that (small) circles in $S^{2}$, including geodesics, can be used to capture the nongeodesic variation. We introduce the principal circles and principal circle mean, analogous to, yet more flexible than, the geodesic principal component and PCmean of Huckemann, Hotz and Munk. These become principal arcs when the manifold is indeed $S^{2}$. For more complex direct product manifolds, we suggest transforming the data points in $S^{2}$ into a linear space by a special mapping utilizing the principal circles. For the other components of the manifold, the tangent space mappings can be used to map the data into a linear space as done in Fletcher et al. 2004. Once manifold-valued data are mapped onto the linear space, then the classical linear PCA can be applied to find principal components in the transformed linear space. The estimated principal components in the linear space can be back-transformed to the manifold, which leads to principal arcs.

We illustrate the potential of our method by an example of m-rep data in Figure 1. Here, m-reps with 15 sample points called atoms model the prostate gland (an organ in the male reproductive system) and come from the simulator developed and analyzed in Jeong et al. 2008. Figure 1 shows that the $S^{2}$ components of the data tend to be distributed along small circles, which frequently are not geodesics. We emphasize the curvature of variation along each sphere by fitting a great circle and a small circle (by the method discussed in Section 2). Our method is adapted to capture this nonlinear (nongeodesic) variation of the data. A potential application of our method is to improve accuracy of segmentation of objects from CT images. Detailed description of the data and results of our analysis can be found in Section 5.

Note that the previous approaches [Fletcher et al. 2004, Huckemann and Ziezold 2006] are defined for general manifolds, while our method focuses on these particular direct product manifolds. Although the method is not applicable for general manifolds, it is useful for this common class of manifolds that is often found in applications. Our results inform our belief that focusing on specific types of manifolds allow more precise and informative statistical modeling than methods that attempt to be fully universal. This happens through using special properties (e.g., presence of small circles) that are not available for all other manifolds.

The rest of the article is organized as follows. We begin by introducing a circle class on $S^{2}$ as an alternative to the set of geodesics. Section 2 discusses principal circles in $S^{2}$, which will be the basis of the special transformation. The first principal circle is defined by the least-squares circle, minimizing the sum of squared residuals. In Section 3 we introduce a data-driven method to decide whether the least-squares circle is appropriate. A recipe for principal arc analysis on direct product manifolds is proposed in Section 4 with discussion on the transformations. A detailed introduction of the space of m-reps and the results from applying the proposed method follow. A novel computational algorithm for the least-squares circles is presented in Section 6. In the Appendix we provide some necessary background for treating direct product manifolds as sample spaces, including the notion of geodesic mean, tangent space, exponential map and log map.

Circle class for nongeodesic variation on S2S^{2}

The circle class includes the simple geodesic case.

Each circle can be parameterized, which leads to an easy interpretation.

There is an orthogonal complement of each circle, which gives two important advantages:

Two orthogonal circles can be used as a basis of a further extension to principal arc analysis.

Building a sensible notion of principal components on $S^{2}$ alone is easily done by utilizing the circles.

A circle that best fits the points $\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\in S^{2}$ is found by minimizing the sum of squared residuals. The residual of $\mathbf{x}_{i}$ is defined as the signed geodesic distance from $\mathbf{x}_{i}$ to the circle $\delta(\mathbf{c},r)$. Then the least-squares circle is obtained by

Note that there are always multiple solutions of (1). In particular, whenever $(\mathbf{c},r)$ is a solution, $(-\mathbf{c},\pi-r)$ also solves the problem as $\delta(\mathbf{c},r)=\delta(-\mathbf{c},\pi-r)$. This ambiguity does not affect any essential result in this paper. Our convention is to use the circle with smaller geodesic radius.

The optimization task (1) is a constrained nonlinear least squares problem. We propose an algorithm to solve the problem that features a simplified optimization task and approximation of $S^{2}$ by tangent planes. The algorithm works in a doubly iterative fashion, which has been shown by experience to be stable and fast. Section 6 contains a detailed illustration of the algorithm.

Analogous to principal geodesics in $S^{2}$, we can define principal circles in $S^{2}$ by utilizing the least-squares circle. The principal circles are two orthogonal circles in $S^{2}$ that best fit the data. We require the first principal circle to minimize the variance of the residuals, so it is the least-squares circle (1). The second principal circle is a geodesic which passes through the center of the first circle and thus is orthogonal at the points of intersection. Moreover, the second principal circle is chosen so that one intersection point is the intrinsic mean [defined in (2) later] of the projections of the data onto the first principal circle.

Based on a belief that the intrinsic (or extrinsic) mean defined on a curved manifold may not be a useful notion of the center point of the data [see, e.g., Huckemann, Hotz and Munk 2010 and Figure 2b], the principal circles do not use the pre-determined means. To develop a better notion of center point, we locate the best 0-dimensional representation of the data in a data-driven manner. Inspired by the PCmean idea of Huckemann, Hotz and Munk 2010, given the first principal circle $\delta_{1}$, the principal circle mean $\mathbf{u}\in\delta_{1}$ is defined (in an intrinsic way) as

where $P_{\delta_{1}}\mathbf{x}$ is the projection of $\mathbf{x}$ onto $\delta_{1}$, that is, the point on $\delta_{1}$ of the shortest geodesic distance to $\mathbf{x}$. Then

as in equation (3.3) of Mardia and Gadsden 1977. We assume that $\mathbf{c}$ is the north pole $\mathbf{e}_{3}$, without losing generality since otherwise the sphere can be rotated. Then

where $\mathbf{u}=(u_{1},u_{2},u_{3})^{\prime}$, $\mathbf{x}=(x_{1},x_{2},x_{3})^{\prime}$ and $\rho_{S^{1}}$ is the geodesic (angular) distance function on $S^{1}$. The optimization problem (2) is equivalent to finding the geodesic mean in $S^{1}$. See equation (17) in the Appendix for computation of the geodesic mean in $S^{1}$.

The second principal circle $\delta_{2}$ is then the geodesic passing through the principal circle mean $\mathbf{u}$ and the center $\mathbf{c}$ of $\delta_{1}$. Denote $\bar{\delta}\equiv\bar{\delta}(\mathbf{x}_{1},\ldots,\mathbf{x}_{n})$ as a combined representation of $(\delta_{1},\mathbf{u})$ or, equivalently, $(\delta_{1},\delta_{2})$.

As a special case, we can force the principal circles to be great circles. The best fitting geodesic is obtained as a solution of the problem (1) with $r=\pi/2$ and becomes the first principal circle. The optimization algorithm for this case is slightly modified from the original algorithm for the varying $r$ case by simply setting $r=\pi/2$. The principal circle mean $\mathbf{u}$ and the $\delta_{2}$ for this case are defined in the same way as in the small circle case. Note that the principal circles with $r=\pi/2$ are essentially the same as the method of Huckemann and Ziezold 2006.

Figure 2 illustrates the advantages of using the circle class to efficiently summarize variation. On four different sets of toy data, the first principal circle $\delta_{1}$ is plotted with principal circle mean $\mathbf{u}$. The first principal geodesics from the methods of Fletcher and Huckemann are also plotted with their corresponding mean. Figure 2a illustrates the case where the data were indeed stretched along a geodesic. The solutions from the three methods are similar to one another. The advantage of Huckemann’s method over Fletcher’s can be found in Figure 2b. The geodesic mean is found far from the data, which leads to poor performance of the principal geodesic analysis, because it considers only great circles passing through the geodesic mean. Meanwhile, the principal circle and Huckemann’s method, which do not utilize the geodesic mean, work well. The case where geodesic mean and any geodesic do not fit the data well is illustrated in Figure 2c, which is analogous to the Euclidean case, where a nonlinear fitting may do a better job of capturing the variation than PCA. To this data set, the principal circle fits best, and our definition of mean is more sensible than the geodesic mean and the PCmean. The points in Figure 2d are generated from the von Mises–Fisher distribution with $\kappa=10$, thus having no principal mode of variation. In this case the first principal circle $\delta_{1}$ follows a contour of the apparent density of the points. We shall discuss this phenomenon in detail in the following section.

Fitting a (small) circle to data on a sphere has been investigated for some time, especially in statistical applications in geology. Those approaches can be distinguished in three different ways, where our choice fits into the first category:

Least-squares of intrinsic residuals: Gray, Geiser and Geiser 1980 formulated the same problem as in (1), finding a circle that minimizes sum of squared residuals, where residuals are defined in a geodesic sense.

Least-squares of extrinsic residuals: A different measure of residual was chosen by Mardia and Gadsden 1977 and Rivest 1999, where the residual of $\mathbf{x}$ from $\delta(c,r)$ is defined by the shortest Euclidean distance between $\mathbf{x}$ and $\delta(\mathbf{c},r)$. Their objective is to find

where $\xi_{i}$ denotes the intrinsic residual. This type of approach can be numerically close to the intrinsic method as $\cos(\xi_{i})=1-\xi_{i}^{2}/2+O(\xi_{i}^{4})$.

Distributional approach: Mardia and Gadsden 1977 and Bingham and Mardia 1978 proposed appropriate distributions to model $S^{2}$-valued data that cluster near a small circle. These models essentially depend on the quantity $\cos(\xi)$, which is easily interpreted in the extrinsic sense but not in the intrinsic sense.

The principal circle and principal circle mean always exist. This is because the objective function (1) is a continuous function of $\mathbf{c}$, with the compact domain $S^{2}$. The minimizer $r$ has a closed-form solution (see Section 6). A similar argument can be made for the existence of $\mathbf{u}$. On the other hand, the uniqueness of the solution is not guaranteed. We conjecture that if the manifold is approximately linear or, equivalently, the data set is well approximated by a linear space, then the principal circle will be unique. However, this does not lead to the uniqueness of $\mathbf{u}$, whose sufficient condition is that the projected data on $\delta_{1}$ is strictly contained in a half-circle [Karcher 1977]. Note that a sufficient condition for the uniqueness of the principal circle is not clear even in the Euclidean case [Chernov 2010].

Suppressing small least-squares circles

When the first principal circle $\delta_{1}$ has a small radius, sometimes it is observed that $\delta_{1}$ does not fit the data in a manner that gives useful decomposition, as shown in Figure 2d. This phenomenon has been also observed for the related principal curve fitting method of Hastie and Stuetzle 1989. We view this as unwanted overfitting, which is indeed a side effect caused by using the full class of circles with free radius parameter instead a class of great circles. In this section a data-driven method to flag this overfitting is discussed. In essence, the fitted small circle is replaced by the best fitting geodesics when the data do not cluster along the circle but instead tend to cluster near the center of the circle.

The problem of suppressing a small circle can be paraphrased as “how to determine whether a nonzero point is a mode of the function $f_{X_{1}|X_{2}=0}$, when observing only a length-biased sample.”

The spectrum from the circle-clustered case (mode at a nonzero point) to the center-clustered case (mode at origin) can be modeled as

where the signal is along a circle with radius $\mu$, and the error accounts for the perpendicular deviation from the circle (see Figure 3). Then, in polar coordinates $(R,\Theta)$, $\Theta$ is uniformly distributed on $(0,2\pi]$ and $R$ is a positive random variable with mean $\mu$. First assume that $R$ follows a truncated Normal distribution with standard deviation $\sigma$, with the marginal p.d.f. proportional to

where $\phi$ is the standard Normal density function. The conditional density $f_{X_{1}|X_{2}=0}$ is then

Nonzero local extrema of $f_{X_{1}|X_{2}=0}$ can be characterized as a function of $(\mu,\sigma)$ in terms of $r_{+},r_{-}=\{\mu\pm\sqrt{(\mu-2\sigma)(\mu+2\sigma)}\}/2$ as follows:

When $\mu>2\sigma$, $f_{X_{1}|X_{2}=0}$ has a local maximum at $r_{+}$, minimum at $r_{-}$.

When $\mu=2\sigma$, $r_{+}=r_{-}=\frac{\mu}{2}.$

When $\mu<2\sigma$, $f_{X_{1}|X_{2}=0}$ is strictly decreasing, for $r\geq 0$.

whenever the ratio $\mu/\sigma>2$, $f_{X_{1}|X_{2}=0}$ has a mode at $r_{+}$.

This idea can be applied for circles in $S^{2}$ with some modification, shown next. We point out that the model (5) is useful for understanding the small circle fitting: signal as a circle with radius $\mu$, and error as the deviation along geodesics perpendicular to the circle. Moreover, a spherically symmetric distribution centered at $\mathbf{c}$ on $S^{2}$ can be mapped to a spherically symmetric distribution on the tangent space at $\mathbf{c}$, preserving the radial distances by the log map (defined in the Appendix). A modification needs to be made on the truncated density $f_{R}$. It is more natural to let the error be so large that the deviation from the great circle is greater than $\mu$. Then the observed value may be found near the opposite side of the true signal, which is illustrated in Figure 3 as the large deviation case. To incorporate this case, we consider a wrapping approach. The distribution of errors (on the real line) is wrapped around the sphere along a great circle through $\mathbf{c}$, and the marginal p.d.f. $f_{R}$ in (6) is modified to

The corresponding conditional p.d.f., $f_{X_{1}|X_{2}=0}^{w}$, is similar to $f_{X_{1}|X_{2}=0}$ and a numerical calculation shows that $f_{X_{1}|X_{2}=0}^{w}$ has a mode at some nonzero point whenever $\mu/\sigma>2.0534$, for $\mu<\pi/2$. In other words, we use the small circle when $\mu/\sigma$ is large. Note that in what follows we only consider the first term ($k=0$) of (7) since other terms are negligible in most situations. We have plotted $f_{X_{1}|X_{2}=0}^{w}$ for some selected values of $\mu$ and $\sigma$ in Figure 4.

With a data set on $S^{2}$, we need to estimate $\mu$ and $\sigma$, or the ratio $\mu/\sigma$. Let $\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\in S^{2}$ and let $\hat{\mathbf{c}}$ be the samples and the center of the fitted circle, respectively. Denote $\xi_{i}$ for the errors of the model (5) such that $\xi_{i}\sim N(0,\sigma^{2})$. Then $r_{i}\equiv\rho(\mathbf{x}_{i},\hat{\mathbf{c}})=|\mu+\xi_{i}|$, which has the folded normal distribution [Leone, Nelson and Nottingham 1961]. Estimation of $\mu$ and $\sigma$ based on unsigned $r_{i}$ is not straightforward. We present two different approaches to this problem.

Robust approach. The observations $r_{1},\ldots,r_{n}$ can be thought of as a set of positive numbers contaminated by the folded negative numbers. Therefore, the left half (near zero) of the data are more contaminated than the right half. We only use the right half of the data, which are less contaminated than the other half. We propose to estimate $\mu$ and $\sigma$ by

where $\mbox{Q}_{3}(\Phi)$ is the third quantile of the standard normal distribution. The ratio can be estimated by $\hat{\mu}/\hat{\sigma}$.

Likelihood approach via EM algorithm. The problem may also be solved by a likelihood approach. Early solutions can be found in Leone, Nelson and Nottingham (Leone, Nelson and Nottingham 1961), Elandt 1961 and Johnson 1962, in which the MLEs were given by numerically solving nonlinear equations based on the sample moments. As those methods were very complicated, we present a simpler approach based on the EM algorithm. Consider unobserved binary variables $s_{i}$ with values $-1$ and $+1$ so that $s_{i}r_{i}\sim N(\mu,\sigma^{2})$. The idea of the EM algorithm is that if we have observed $s_{i}$, then the maximum likelihood estimator of $\vartheta=(\mu,\sigma^{2})$ would be easily obtained. The EM algorithm is an iterative algorithm consisting of two steps. Suppose that the $k$th iteration produced an estimate $\hat{\vartheta}_{k}$ of $\vartheta$. The E-step is to impute $s_{i}$ based on $r_{i}$ and $\hat{\vartheta}_{k}$ by forming a conditional expectation of log-likelihood for $\vartheta$,

where $f$ is understood as an appropriate density function, and $p_{i(k)}$ is easily computed as

The M-step is to maximize $Q(\vartheta)$ whose solution becomes the next estimator $\hat{\vartheta}_{k+1}$. Now the $(k+1)$th estimates are calculated by a simple differentiation and given by

With the sample mean and variance of $r_{1},\ldots,r_{n}$ as an initial estimator $\hat{\vartheta}_{0}$, the algorithm iterates E-steps and M-steps until the iteration changes the estimates less than a predefined criteria (e.g., $10^{-10}$). $\mu/\sigma$ is estimated by the ratio of the solutions.

Comparison. Performance of these estimators are now examined by a simulation study. Normal random samples are generated with ratios $\mu/\sigma$ being 0, 1, 2 or 3, representing the transition from the center-clustered to circle-clustered case. For each ratio, $n=50$ samples are generated, from which $\hat{\mu}/\hat{\sigma}$ is estimated. These steps are repeated 1000 times to obtain the sampling variation of the estimates. We also study the $n=1000$ case in order to investigate the consistency of the estimators. The results are summarized in Figure 5 and Table 1.

The distribution of estimators are shown for $n=50,1000$ in Figure 5 and the proportion of estimators greater than 2 is summarized in Table 1. When $n=1000$, both estimators are good in terms of the proportion of correct answers. In the following, the proportions of correct answers are corresponding to $n=50$ case. The top left panel in Figure 5 illustrates the circle-centered case with ratio 3. The estimated ratios from the robust approach give correct solutions (greater than 2) 95% of the time (98.5% for likelihood approach). For the borderline case (ratio 2, top right), the small circle will be used about half the time. The center-clustered case is demonstrated with the true ratio 1, that also gives a reasonable answer (proportion of correct answers 95.3% and 94.8% for the robust and likelihood answers respectively). It can be observed that when the true ratio is zero, the robust estimates are far from 0 (the bottom right in Figure 5). However, this is expected to occur because the proportion of uncontaminated data is low when the ratio is too small. However, those ‘inaccurate’ estimates are around 1 and less than 2 most of the time, which leads to ‘correct’ answers. The likelihood approach looks somewhat better with more hits near zero, but an asymptotic study [Johnson 1962] showed that the variance of the maximum likelihood estimator converges to infinity when the ratio tends to zero, as glimpsed in the long right tail of the simulated distribution.

In summary, we recommend use of the robust estimators (8), which are computationally light, straightforward and stable for all cases.

In addition, we point out that Gray, Geiser and Geiser 1980 and Rivest 1999 proposed to use a goodness-of-fit statistic to test whether the small circle fit is better than a geodesic fit. Let $r_{g}$ and $r_{c}$ be the sums of squares of the residuals from great and small circle fits. They claimed that $V=(n-3)(r_{g}-r_{c})/r_{c}$ is approximately distributed as $F_{1,n-3}$ for a large $n$ if the great circle was true. However, this test does not detect the case depicted in Figure 2d. The following numerical example shows the distinction between our approach and the goodness-of-fit approach.

Consider the sets of data depicted in Figure 2. The goodness-of-fit test gives $p$-values of $0.51$, $0.11359$, $0$ and $0.0008$ for (a)–(d), respectively. The estimated ratios $\mu/\sigma$ are $14.92,16.89$, $14.52$ and $1.55$. Note that for (d), when the least-squares circle is too small, our method suggests to use a geodesic fit over a small circle while the goodness-of-fit test gives significance of the small circle. The goodness-of-fit method is not adequate to suppress the overfitting small circle in a way we desire.

A referee pointed out that the transition of the principal circle between great circle and small circle is not continuous. Specifically, when the data set is perturbed so that the principal circle becomes too small, then the principal circle and principal circle mean are abruptly replaced by a great circle and geodesic mean. As an example, we have generated a toy data set spread along a circle with some radial perturbation. The perturbation is continuously inflated, so that with large inflation, the data are no longer circle-clustered. In Figure 6 the $\widehat{\mu/\sigma}$ changes smoothly, but once the estimate hits 2 (our criterion), there is a sharp transition between small and great circles. Sharp transitions do naturally occur in the statistics of manifold data. For example, even the simple geodesic mean can exhibit a major discontinuous transition resulting from an arbitrarily small perturbation of the data. However, the discontinuity between small and great circles does seem more arbitrary and thus may be worth addressing. An interesting open problem is to develop a blended version of our two solutions, for values of $\widehat{\mu/\sigma}$ near 2, which could be done by fitting circles with radii that are smoothly blended between the small circle radius and $\pi/2$.

Principal arc analysis on direct product manifolds

Consider a data set $x_{1},\ldots,x_{n}\in M$, where $x_{i}\equiv(x_{i}^{1},\ldots,x_{i}^{d})$ such that $x_{i}^{j}\in M_{j}$. Denote $d_{0}\geq d$ for the intrinsic dimension of $M$. The geodesic mean $\bar{x}$ of the data is defined component-wise for each simple manifold $M_{j}$. Similarly, the tangent plane at $\bar{x}$, $T_{\bar{x}}M$, is also defined marginally, that is, $T_{\bar{x}}M$ is a direct product of tangent spaces of the simple manifolds. This tangent space gives a way of applying Euclidean space-based statistical methods by mapping the data onto $T_{\bar{x}}M$. We can manipulate this approximation of the data component-wise. In particular, the marginal data on the $S^{2}$ components can be represented in a linear space by a transformation $h_{\bar{\delta}}$, depending on the principal circles, that differs from the tangent space approximation.

$\delta_{1}$ is mapped onto the $x$-axis, and

$\delta_{2}$ is mapped onto the $y$-axis.

Two reasonable choices of the mapping $h_{\bar{\delta}}$ will be discussed in Section 4.1 in detail.

The mapping $h_{\bar{\delta}}$ and the tangent space projection together give a linear space representation of the data where the Euclidean PCA is applicable. The line segments corresponding to the sample principal component direction of the transformed data can be mapped back to $M$, and become the principal arcs.

A procedure for principal arc analysis is as follows:

For each $j$ such that $M_{j}$ is $S^{2}$, compute principal circles $\bar{\delta}=\bar{\delta}(x_{1}^{j},\ldots,x_{2}^{j})$ and the ratio $\widehat{\mu/\sigma}$. If the ratio is greater than the predetermined value $\varepsilon=2$, then $\bar{\delta}$ is adjusted to be great circles as explained in Section 2.

where $\operatorname{Log}_{\bar{x}^{j}}$ and $h_{\bar{\delta}}$ are defined in the Appendix and Section 4.1, respectively.

The $k$th principal arc is obtained by mapping the direction vectors $\mathbf{v}_{k}$ onto $M$ by the inverse of $h$, which can be computed component-wise.

Principal arc analysis for data on direct product manifolds often results in a concise summary of the data. When we observe a significant variation along a small circle of a marginal $S^{2}$, that is most likely not a random artifact but, instead, the result of a signal driving the circular variation. Nongeodesic variation of this type is well captured by our method.

Principal arcs can be used to reduce the intrinsic dimensionality of $M$. Suppose we want to reduce the dimension by $k$, where $k$ can be chosen by inspection of the scree plot. Then each data point $x$ is projected to a $k$-dimensional submanifold $M_{0}$ of $M$ in such a way that

Projection. The first approach is based on the projection of $\mathbf{x}$ onto $\delta_{1}$, defined in (3), and a residual $\xi$. The signed distance from $\mathbf{u}$ to $P_{\delta_{1}}x$, whose unsigned version is defined in (4), becomes the $x$-coordinate, while the residual $\xi$ becomes the $y$-coordinate. This approach has the same spirit as the model for the circle class (5), since the direction of the signal is mapped to the $x$-axis, with the perpendicular axis for errors.

The projection $h_{\bar{\delta}}(\mathbf{x})$ that we define here is closely related to the spherical coordinate system. Assume $\mathbf{c}=\mathbf{e}_{3}$, and $\mathbf{u}$ is at the Prime meridian (i.e., on the $x-z$ plane). For $\mathbf{x}$ and its spherical coordinates $(\phi,\theta)$ such that $\mathbf{x}=(x_{1},x_{2},x_{3})=(\cos\phi\sin\theta,\cos\phi\sin\theta,\cos\theta)$,

where $\theta_{\mathbf{u}}=\cos^{-1}(u_{3})$ is the latitude of $\mathbf{u}$. The set of $h_{\bar{\delta}}(\mathbf{x}_{i})$ has mean zero because the principal circle mean $\mathbf{u}$ has been subtracted.

In many cases, both projection and conformal $h_{\bar{\delta}}$ give better representations than just using the tangent space. Figure 7 illustrates the image of $h_{\bar{\delta}}$ with the toy data set depicted in Figure 2c. The tangent space mapping is also plotted for comparison. The tangent space mapping leaves the curvy form of variation, while both $h_{\bar{\delta}}$’s capture the variation and lead to an elliptical distribution of the transformed data.

Application to m-rep data

In this section an application of Principal Arc Analysis to the medial representation (m-rep) data example, introduced below in more detail, is described.

The m-rep gives an efficient way of representing 2- or 3-dimensional objects. The m-rep is constructed from the medial axis, which is a means of representing the middle of geometric objects. The medial axis of a 3-dimensional object is formed by the centers of all spheres that are interior to objects and tangent to the object boundary at two or more points. In addition, the medial description is defined by the centers of the inscribed spheres and by the associated vectors, called spokes, from the sphere center to the two respective tangent points on the object boundary. The medial axis is sampled over an approximately regular lattice and the elements of the lattice are called medial atoms. A medial atom consists of the location of the atom combined with two equal-length spokes, defined as a 4-tuple:

An important topic in medical imaging is developing segmentation methods of 3D objects from CT images; see Cootes and Taylor 2001 and Pizer et al. 2007. A popular approach is similar to a Bayesian estimation scheme, where the knowledge of anatomic geometries is used (as a prior) together with a measure of how the segmentation matches the image (as a likelihood). A prior probability distribution is modeled using m-reps as a means of measuring geometric atypicality of a segmented object. PCA-like methods (including PAA) can be used to reduce the dimensionality of such a model. A detailed description can be found in Pizer et al. 2007.

2 Simulated m-rep object

The data set partly plotted in Figure 1 is from the generator discussed in Jeong et al. 2008. It generates random samples of objects whose shape changes and motions are physically modeled (with some randomness) by anatomical knowledge of the bladder, prostate and rectum in the male pelvis. Jeong et al. have proposed and used the generator to estimate the probability distribution model of shapes of human organs.

In the data set of 60 samples of prostate m-reps we studied, the major motion of prostate is a rotation. In some $S^{2}$ components, the variation corresponding to the rotation is along a small circle. Therefore, PAA should fit better for this type of data than principal geodesics. To make this advantage more clear, we also show results from a data set by removing the location and the spoke length information from the m-reps, the sample space of which is then $\{S^{2}\}^{30}$.

We have applied PAA as described in the previous section. The ratios $\mu/\sigma$, estimated for the 30 $S^{2}$ components, are in general large (with minimum 21.2, median 44.1 and maximum 118), which suggests use of small circles to capture the variation.

Figure 9 shows the proportion of the cumulative variances, as a function of number of components, from the Principal Geodesic Analysis (PGA) of Fletcher et al. 2004 and PAA. In both cases, the first principal arc leaves smaller residuals than the first principal geodesic. What is more important is illustrated in the scatterplots of the data projected onto the first two principal components. The quadratic form of variation that requires two PGA components is captured by a single PAA component.

The probability distribution model estimated by principal geodesics is qualitatively different from the distribution estimated by PAA. Although the difference in the proportion of variance captured is small, the resulting distribution from PAA is no longer elliptical. In this sense, PAA gives a convenient way to describe a nonelliptical distribution by, for example, a Normal density.

3 Prostate m-reps from real patients

We also have applied PAA to a prostate m-rep data set from real CT images. Our data consist of five patients’ image sets, each of which is a series of CT scans containing prostate taken during a series of radiotherapy treatments [Merck 2008]. The prostate in each image is manually segmented by experts and an m-rep model is fitted. The patients, coded as 3106, 3107, 3109, 3112 and 3115, have different numbers of CT scans (17, 12, 18, 16 and 15, respectively). We have in total 78 m-reps.

The proportion of variation captured in the first principal arc is 40.89%, slightly higher than the 40.53% of the first principal geodesic. Also note that the estimated probability distribution model from PAA is different from that of PGA. In particular, PAA gives a better separation of patients in the first two components, as depicted in the scatter plots (Figure 10).

Doubly iterative algorithm to find the least-squares small circle

The rotation operator $q_{\mathbf{c}}$ can be represented by a rotation matrix. For $\mathbf{c}=(c_{x},c_{y},c_{z})^{\prime}$, the rotation $q_{\mathbf{c}}$ is equivalent to rotation through the angle $\theta=\cos^{-1}(c_{z})$ about the axis $\mathbf{u}=(c_{y},-c_{x},0)^{\prime}/\sqrt{1-c_{z}^{2}}$, whenever $\mathbf{c}\neq\pm\mathbf{e}_{3}$. When $\mathbf{c}=\pm\mathbf{e}_{3}$, $\mathbf{u}$ is set to be $\mathbf{e}_{1}$. It is well known that a rotation matrix with axis $\mathbf{u}=(u_{x},u_{y},u_{z})^{\prime}$ and angle $\theta$ in radians is, for $c=\cos(\theta)$, $s=\sin(\theta)$ and $v=1-\cos(\theta)$,

for $\theta=\cos^{-1}(x_{3})$. See (18)–(19) in the Appendix for $\operatorname{Exp}_{\mathbf{e}_{3}}$ and $\operatorname{Log}_{\mathbf{e}_{3}}$. Note that $\operatorname{Log}_{\mathbf{c}}(\mathbf{c})=\mathbf{0}$ and $\operatorname{Log}_{\mathbf{c}}$ is not defined for the antipodal point of $\mathbf{c}$.

which is to find the least-squares circle centered at $\mathbf{v}$ with radius $r$. The general circle fitting problem is discussed in, for example, Umbach and Jones 2003 and Chernov 2010. This problem is much simpler than (1) because it is an unconstrained problem and the number of parameters to optimize is decreased by 1. Moreover, optimal solution of $r$ is easily found as

when $\mathbf{v}$ is given. Note that for great circle fitting, we can simply put $\hat{r}=\pi/2$. Although the problem is still nonlinear, one can use any optimization method that solves nonlinear least squares problems. We use the Levenberg–Marquardt algorithm, modified by Fletcher 1971 [see Chapter 4 of Scales 1985 and Chapter 3 of Bates and Watts 1988], to minimize (13) with $r$ replaced by $\hat{r}$. One can always use $\mathbf{v}=\mathbf{0}$ as an initial guess since $\mathbf{0}=\operatorname{Log}_{\mathbf{c}}(\mathbf{c})$ is the solution from the previous (outer) iteration.

The algorithm is now summarized as follows:

Given $\{\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\}$, $\mathbf{c}_{0}=\mathbf{x}_{1}$.

If $\|\mathbf{v}\|<\varepsilon$, then iteration stops with the solution $\hat{\mathbf{c}}=\mathbf{c}_{k}$, $r=\hat{r}$ as in (14). Otherwise, $\mathbf{c}_{k+1}=\operatorname{Exp}_{\mathbf{c}_{k}}(\mathbf{v})$ and go to step 2.

Note that the radius of the fitted circle in $T_{\mathbf{c}}$ is the same as the radius of the resulting small circle. There could be many variations of this algorithm: as an instance, one can elaborate the initial value selection by using the eigenvector of the sample covariance matrix of $\mathbf{x}_{i}$’s, corresponding to the smallest eigenvalue as done in Gray, Geiser and Geiser 1980. Experience has shown that the proposed algorithm is stable and speedy enough. Gray, Geiser and Geiser proposed to solve (1) directly, which seems to be unstable in some cases.

The idea of the doubly iterative algorithm can be applied to other optimization problems on manifolds. For example, the geodesic mean is also a solution of a nonlinear minimization, where the nonlinearity comes from the use of the geodesic distance. This can be easily solved by an iterative approximation of the manifold to a linear space [See Chapter 4 of Fletcher 2004], which is the same as the gradient descent algorithms [Pennec 2006, Le 2001]. Note that the proposed algorithm, like other iterative algorithms, only finds one solution even if there are multiple solutions.

Appendix: Some background on direct product manifold

We give some necessary geometric background on direct product manifolds. Precise definitions and geometric discussions on a richer class of manifold, Riemannian manifold, can be found in Boothby 1986 and Helgason 2001.

A $d$-dimensional manifold can be thought of as a curved surface embedded in a Euclidean space of higher dimension $d^{\prime}$ (${\geq}d$). The manifold is required to be smooth, that is, infinitely differentiable, so that a sufficiently small neighborhood of any point on the manifold can be well approximated by a linear space. The tangent space at a point $p$ of a manifold $M$, $T_{p}M$, is defined as a linear space of dimension $d$ which is tangent to $M$ at $p$. The notion of distance on $M$ is handled by a Riemannian metric, which is a metric of tangent spaces. In particular, the geodesic distance function $\rho_{M}(p,q)$ is roughly defined as the length of the shortest curve joining $p,q\in M$.

Geodesic distance function. The geodesic distance between $p\equiv(p^{1},\ldots,p^{m})$ and $q\equiv(q^{1},\ldots,q^{m})$ is defined by

Geodesic mean and variance. The geodesic mean of a set of points in $M$, also referred to as the intrinsic mean, is also calculated component-wise. The geodesic mean of $x_{1},\ldots,x_{n}\in M$ is the minimizer in $M$ of the sum of squared geodesic distances to the data. Thus, the geodesic mean is defined as

In fact, each $\bar{x}^{i}$ of $\bar{x}=(\bar{x}^{1},\ldots,\bar{x}^{m})$ is the geodesic mean of $x_{1}^{i},\ldots,x_{n}^{i}\in M_{i}$. The geodesic mean of $\theta_{1},\ldots,\theta_{n}\in[0,2\pi)\cong S^{1}$ is found by examining a candidate set consisting of

A related notion is geodesic variance. A sample geodesic variance is defined by the average squared geodesic distances to the geodesic mean, that is, $\frac{1}{n}\sum_{i=1}^{n}\rho^{2}_{M}(\bar{x},x_{i})$. When $M$ is indeed the Euclidean space, the geodesic variance is the same as the total variance (the trace of the variance–covariance matrix).

This equation can be understood as a rotation of the base point $\mathbf{p}$ to the direction of $\mathbf{v}$ with angle $\|\mathbf{v}\|$. The corresponding log map for a point $\mathbf{x}=(x_{1},x_{2},x_{3})\in S^{2}$ is given by

References