1 Introduction
This paper aims to introduce an R package spherepc that considers several dimension reduction techniques on a sphere, which encompass recently developed approaches such as SPC and LPG as well as some existing methods, and discuss how to implement these methods through spherepc.
Dimension reduction methods are widely used in various fields, including statistics and machine learning, by efficiently compressing data and removing noise (Benner et al. 2005). As one of the dimension reduction methods, the principal curves of Hastie and Stuetzle (1989) are suitable for fitting a curve or a surface of data in Euclidean space, which go through the middle of the data. Hauberg (2016) proposed an algorithm to find the principal curves in Riemannian manifolds based on the concept of the original principal curves. However, the principal curves proposed by Hauberg (2016) no longer represent the data continuously because of the approximation of the projection step required to fit the curves.
Recently, Lee et al. (2021a) proposed a new method, termed spherical principal curves (SPC), that constructs principal curves, ensuring a stationary property on spheres. SPC is useful for representing circular or waveform data with smaller reconstruction errors than conventional methods, including principal geodesic analysis (Fletcher et al. 2004), exact principal circle (Lee et al. 2021a), and principal curves proposed by Hauberg (2016). However, SPC has the disadvantage of being sensitive to initialization. As a result, there are some data structures that SPC does not apply to, for example, data with spirals, zigzags, or branches like tree-shape. A localized version of SPC called local principal geodesics (LPG) is being developed to resolve such a problem. A function for LPG is also provided in the package spherepc. Research on the LPG is underway in progress.
To the best of our knowledge, no available R packages offer the methods of dimension reduction and principal curves on a sphere. The existing R packages providing principal curves, such as princurve (Hastie and Weingessel 2015) and LPCM (Einbeck et al. 2015), are available only on Euclidean space, not on a sphere or (Riemannian) manifold. In addition, most dimension reduction methods on manifolds (Huckemann et al. 2010; Panaretos et al. 2014; Liu et al. 2017) involve somewhat complex optimizations. The proposed package spherepc for R provides the state-of-the-art principal curve technique on the sphere (Lee et al. 2021a) and comprehensively collects and implements the existing methods (Fletcher et al. 2004; Hauberg 2016).
The rest of this paper is organized as follows. The following section
introduces the existing methods for dimension reduction on the sphere
and relevant functions covered in the package
spherepc, which is
available on CRAN. Furthermore, their usages are discussed with examples
in detail. Then, the spherical principal curves proposed by (Lee et al. 2021a)
and principal curves of (Hauberg 2016) are briefly described. In
addition, implementations of the SPC() and SPC.Hauberg() functions
in the spherepc are
presented. The subsequent section discusses the local principal
geodesics (LPG) with the implementation of various simulated data,
demonstrating its promising usability. In the application session, all
the mentioned methods are performed to analyze real seismological data.
Finally, conclusions are given in the last section.
2 Existing methods
Principal geodesic analysis
Principal geodesic analysis (PGA) proposed by Fletcher et al. (2004) can be regarded as a generalization of principal component analysis (PCA) to Riemannian manifolds. In particular, Fletcher et al. (2004) performed dimension reduction of data on the Cartesian product space of the manifolds. In detail, the data are projected onto the tangent spaces at the intrinsic means of each component of the manifolds; thus, the given data are approximated as points on Euclidean vector space, and subsequently, PCA is applied to the points. As a result, the dimension reduction can be performed through the inverse of the tangent projections.
The principal geodesic analysis can be implemented by the PGA()
function available in the
spherepc. The detailed
usage of the PGA() function is described as follows.
PGA(data, col1 = "blue", col2 = "red")Before using the PGA() function, it requires loading the packages
rgl (Adler and Murdoch 2020),
sphereplot
(Robotham 2013), and
geosphere
(Hijmans et al. 2017). The following codes yield an implementation of the
PGA() function.
#### for all simulated datasets, longitude and latitude are expressed in degrees
#### example 1: half-great circle data
> circle <- GenerateCircle(c(150, 60), radius = pi/2, T = 1000)
> sigma <- 2 # noise level
> half.circle <- circle[circle[, 1] < 0, , drop = FALSE]
> half.circle <- half.circle + sigma * rnorm(nrow(half.circle))
> PGA(half.circle)
#### example 2: S-shaped data
# the dataset consists of two parts: lon ~ Uniform[0, 20],
# lat = sqrt(20 * lon - lon^2) + N(0, sigma^2),
# lon ~ Uniform[-20, 0], lat = -sqrt(-20 * lon - lon^2) + N(0, sigma^2)
> n <- 500
> sigma <- 1 # noise level
> lon <- 60 * runif(n)
> lat <- (60 * lon - lon^2)^(1/2) + sigma * rnorm(n)
> simul.S1 <- cbind(lon, lat)
> lon2 <- -60 * runif(n)
> lat2 <- -(-60 * lon2 - lon2^2)^(1/2) + sigma * rnorm(n)
> simul.S2 <- cbind(lon2, lat2)
> simul.S <- rbind(simul.S1, simul.S2)
> PGA(simul.S)Because a principal geodesic is always a great circle, the PGA()
function is suitable for identifying the global data trend. The
implementations of half-circle and S-shaped data are displayed in
Figure 1, where the principal geodesic properly extracts the
global trends in the half-great circle and S-shaped data, while it
cannot identify the circular variations in the S-shaped case. In
addition, the arguments and outputs of the PGA() function are
described in Tables 1 and 2.
| Argument | Description |
|---|---|
data |
matrix or data frame consisting of spatial locations with two columns. Each row represents longitude and latitude (denoted by degrees). |
col1 |
color of data. The default is blue. |
col2 |
color of the principal geodesic line. The default is red. |
| Output | Description |
|---|---|
plot |
plotting of the result in 3D graphics. |
line |
spatial locations (longitude and latitude by degrees) of points in the principal geodesic line. |
Principal circle
In a spherical surface, as shown in Figure 1, the principal
geodesic analysis always results in a great circle, which cannot be
sufficient to identify the non-geodesic structure of data. The circle on
a sphere that minimizes a reconstruction error is called a principal
circle, where the reconstruction error is defined as the total sum of
squares of geodesic distances between the circle and data points.
However, the existing method for generating the principal circle is
still based on the tangent space approximation and its inverse process,
thereby leading to numerical errors. Lee et al. (2021a) have proposed an exact
principal circle in an intrinsic way and its practical algorithm based
on gradient descent. The details are described in Section 3 of Kim et al. (2020)
and Appendix B of Lee et al. (2021b). The
spherepc package
provides the PrincipalCircle() function to implement the intrinsic
principal circle. Its usage is followed by
PrincipalCircle(data, step.size = 1e-3, thres = 1e-5, maxit = 10000).| Argument | Description |
|---|---|
data |
matrix or data frame consisting of spatial locations (longitude and latitude denoted by degrees) with two columns. |
step.size |
step size of gradient descent algorithm. For convergence of the algorithm, step.size is recommended to be below 0.01. The default is 1e-3. |
thres |
threshold of the stopping condition. The default is 1e-5. |
maxit |
maximum number of iterations. The default is 10000. |
The arguments of the PrincipalCircle() are described in
Table 3, and its output is a three-dimensional
vector, where the first and second components are longitude and latitude
(represented by degrees), respectively. The last one is the radius of
the principal circle. To display the circle, the GenerateCircle()
function should be implemented. Its usage is followed by
GenerateCircle(center, radius, T = 1000).The output of the GenerateCircle() function is a matrix consisting of
spatial locations (longitude and latitude by degrees) with two columns,
which can be plotted by the sphereplot::rgl.sphgrid() and
sphereplot::rgl.sphpoints() functions from the
sphereplot package
(Robotham 2013). Note that the
sphereplot package
depends on the rgl package
(Adler and Murdoch 2020). The detailed arguments of the GenerateCircle() function
are described in Table 4.
| Argument | Description |
|---|---|
center |
center of circle with spatial locations (longitude and latitude denoted by degrees). |
radius |
radius of circle. It should be range from 0 to \(\pi\). |
T |
the number of points that make up a circle. The points in a circle are equally spaced. The default is 1000. |
The following codes implement principal circles by the
PrincipalCircle() and GenerateCircle() functions.
## for all the following examples, longitude and latitude are denoted by degrees
#### example 1: half-great circle data
> circle <- GenerateCircle(c(150, 60), radius = pi/2, T = 1000)
> half.great.circle <- circle[circle[, 1] < 0, , drop = FALSE]
> sigma <- 2 # noise level
> half.great.circle <- half.great.circle + sigma * rnorm(nrow(half.great.circle))
## find a principal circle
> PC <- PrincipalCircle(half.great.circle)
> result <- GenerateCircle(PC[1:2], PC[3], T = 1000)
## plot the half-great circle data and principal circle
> sphereplot::rgl.sphgrid(col.lat = "black", col.long = "black")
> sphereplot::rgl.sphpoints(half.great.circle, radius = 1, col = "blue", size = 9)
> sphereplot::rgl.sphpoints(result, radius = 1, col = "red", size = 6)
#### example 2: circular data
> n <- 700 # the number of samples
> sigma <- 5 # noise level
> x <- seq(-180, 180, length.out = n)
> y <- 45 + sigma * rnorm(n)
> simul.circle <- cbind(x, y)
## find a principal circle
> PC <- PrincipalCircle(simul.circle)
> result <- GenerateCircle(PC[1:2], PC[3], T = 1000)
## plot the circular data and principal circle
> sphereplot::rgl.sphgrid(col.lat = "black", col.long = "black")
> sphereplot::rgl.sphpoints(simul.circle, radius = 1, col = "blue", size = 9)
> sphereplot::rgl.sphpoints(result, radius = 1, col = "red", size = 6)The results of the principal circle are shown in Figure 2. As one can see, the principal circle identifies the circular patterns of the noisy half-great circle and circular dataset well.
3 Spherical principal curves
Principal curves proposed by Hastie and Stuetzle (1989) can be considered as a nonlinear generalization of the principal component analysis in the sense that the principal curves pass through the middle of given data and reserve a stationary property. The curve is a smooth function from a one-dimensional closed interval to a given space; then, a curve \(f\) is said to be a principal curve of \(X\) or self-consistent if the curve satisfies \[f(\lambda) = \mathbb{E}\big(X \ | \ \lambda_{f}(X)=\lambda\big),\] where \(f(\lambda_f(x))\) is the closest (projection) point in the curve \(f\) from the point \(x\).
Hauberg (2016) provided an algorithm for principal curves on Riemannian manifolds. However, Hauberg (2016) used approximations for finding the closest point of each data point, which may lead to numerical errors. Recently, Lee et al. (2021a) presented theoretical results of principal curves on spheres and a practical algorithm for constructing principal curves without any approximations, called spherical principal curves (SPC), thereby causing the given data to be represented more precisely and smoothly compared to principal curves of Hauberg (2016). In the both ways of extrinsic and intrinsic approaches, the method of SPC updates curves on the spherical surfaces to represent the given data and fits curves that satisfy the stationary conditions. For more details, refer to Kim et al. (2020) or Lee et al. (2021a).
The package spherepc
provides the SPC() function for implementing spherical principal
curves and the SPC.Hauberg() function for principal curves of
Hauberg (2016). The usage of the SPC() function is as follows.
SPC(data, q = 0.05, T = nrow(data), step.size = 1e-3, maxit = 30,
type = "Intrinsic", thres = 1e-2, deletePoints = FALSE,
plot.proj = FALSE, kernel = "quartic", col1 = "blue",
col2 = "green", col3 = "red").The usage of the SPC.Hauberg() function is the same as that of the
SPC() function. Before implementing the SPC() and SPC.Hauberg()
functions, it requires loading the
rgl (Adler and Murdoch 2020),
sphereplot
(Robotham 2013), and
geosphere
(Hijmans et al. 2017) packages. To implement the SPC() and SPC.Hauberg()
functions, we consider the waveform data used in Liu et al. (2017),
Kim et al. (2020), and Lee et al. (2021a). The generating equation of waveform is
\[\phi = \alpha\cdot \sin(f \theta\cdot \pi/180) + 10,\]
where \(\phi\), \(\theta\), \(\alpha\), and \(f\) denote the longitude, latitude
in degrees, amplitude and frequency of the waveform, respectively.
\(\theta\) is uniformly sampled from the interval \([-180,\, 180]\) and a
Gaussian random noise from \(N(0,\, \sigma^2)\) is added on each \(\phi\)
where \(\sigma = 2,\, 10\). The generating waveform data and
implementations of the SPC() and SPC.Hauberg() functions are as
follows.
#### longitude and latitude are expressed in degrees
#### example: waveform data
> n <- 200
> alpha <- 1/3; freq <- 4 # amplitude and frequency of wave
> sigma1 <- 2; sigma2 <- 10 # noise levels
> lon <- seq(-180, 180, length.out = n) # uniformly sampled longitude
> lat <- alpha * 180/pi * sin(freq * lon * pi/180) + 10 # latitude vector
## add Gaussian noises on the latitude vector
> lat1 <- lat + sigma1 * rnorm(length(lon))
> lat2 <- lat + sigma2 * rnorm(length(lon))
> wave1 <- cbind(lon, lat1); wave2 <- cbind(lon, lat2)
## implement Hauberg's principal curves to the waveform data
> SPC.Hauberg(wave1, q = 0.05)
## implement SPC to the (noisy) waveform data
> SPC(wave1, q = 0.05)
> SPC(wave2, q = 0.05)The above codes generate the results in Figure 3. As one
can see, the SPC() and SPC.Hauberg() functions identify the waveform
pattern of the simulated data. Especially, the SPC() generates a
smoother curve. The detailed arguments and outputs of the SPC() are
described in Tables 5 and 6,
respectively, which are the same for the SPC.Hauberg().
| Argument | Description |
|---|---|
data |
matrix or data frame consisting of spatial locations with two columns. Each row represents longitude and latitude (denoted by degrees). |
q |
numeric value of the smoothing parameter. The parameter plays the same role, as the bandwidth does in kernel regression, in the SPC function. The value should be a numeric value between 0.01 and 0.5. The default is 0.1. |
T |
the number of points making up the resulting curve. The default is 1000. |
step.size |
step size of the PrincipalCircle function. The default is 0.001. The resulting principal circle is used as an initialization of the SPC function. |
maxit |
maximum number of iterations. The default is 30. |
type |
type of mean on the sphere. The default is "Intrinsic" and the other choice is "Extrinsic". |
thres |
threshold of the stopping condition. The default is 0.01. |
deletePoints |
logical value. The argument is an option of whether to delete points or not. If deletePoints is FALSE, this function leaves the points in curves that do not have adjacent data for each expectation step. As a result, the function usually returns a closed curve, i.e. a curve without endpoints. If deletePoints is TRUE, this function deletes the points in curves that do not have adjacent data for each expectation step. As a result, The SPC function usually returns an open curve, i.e. a curve with endpoints. The default is FALSE. |
plot.proj |
logical value. If the argument is TRUE, the projection line for each data point is plotted. The default is FALSE. |
kernel |
kind of kernel function. The default is the quartic kernel, and the alternative is indicator or Gaussian. |
col1 |
color of data. The default is blue. |
col2 |
color of points in principal curves. The default is green. |
col3 |
color of resulting principal curves. The default is red. |
| Output | Description |
|---|---|
plot |
plotting of the result in 3D graphics. |
prin.curves |
spatial locations (denoted by degrees) of points in the resulting principal curves. |
line |
connecting lines between points in prin.curves. |
converged |
whether or not the algorithm converged. |
iteration |
the number of iterations of the algorithm. |
recon.error |
sum of squared distances between the data and their projections. |
num.dist.pt |
the number of distinct projections. |
Options for spherical principal curves
There are some options for the SPC() and SPC.Hauberg() functions. In
particular, we implement using the arguments plot.proj and
deletePoints, described in Table 5. If
plot.proj = TRUE is used, then the projection line for each data point
is plotted. If the argument deletePoints = TRUE is performed, the
SPC() function deletes the points in curves that do not have adjacent
data for each expectation step required to fit the principal curves,
returning an open curve, i.e., a curve with endpoints. As a result,
the principal curves are more parsimonious since a redundant part of the
resulting curves is removed. The SPC.Hauberg() function also contains
the same options. For implementing these two arguments, the following
codes are performed through real earthquake data.
> data(Earthquake)
# collect spatial locations (longitude and latitude denoted by degrees) of data
> earthquake <- cbind(Earthquake$longitude, Earthquake$latitude)
#### example 1: plot the projection lines (option of plot.proj)
> SPC(earthquake, q = 0.1, plot.proj = TRUE)
#### example 2: open principal curves (option of deletePoints)
> SPC(earthquake, q = 0.04, deletePoints = TRUE)The results are illustrated in Figure 4. The left
panel shows a closed principal curve (red) with projection lines (black)
of each data point onto the curve, and the right panel displays an open
principal curve due to the option deletePoints = TRUE. It is a
parsimonious result because the redundant part on the upper right side
of the sphere is removed.
deletePoints=TRUE. The less q is, the more the curve overfits
the data.
4 Local principal geodesics
Suppose that observations have a non-geodesic structure. Then the PGA may not be beneficial to represent such data because PGA always results in a geodesic line. To overcome this problem, we consider performing PGA locally and repeatedly to detect the non-geodesic and complex structures of data, which can be interpreted as a localized version of the PGA and SPC. The newly proposed method is called local principal geodesics (LPG). The main idea behind the LPG is that non-geodesic structures can be regarded as a part of geodesic when viewed locally. Although there is no reference to the LPG because research on LPG is underway, there is a localized principal curve method on Euclidean space (Einbeck et al. 2005), which is similar to LPG and may share some motivation with the LPG. For more details, refer to (Einbeck et al. 2005).
The package spherepc
offers the LPG() function to recognize various data structures, such
as spirals, zigzag, and tree data. The usage of the function is
LPG(data, scale = 0.04, tau = scale/3, nu = 0, maxpt = 500,
seed = NULL, kernel = "indicator", thres = 1e-4, col1 = "blue",
col2 = "green", col3 = "red").Like the previous functions, before the LPG() function is implemented,
it requires to load the rgl
(Adler and Murdoch 2020),
sphereplot
(Robotham 2013), and
geosphere
(Hijmans et al. 2017) packages. The detailed arguments and outputs of this
function are described in Tables 7 and 8.
We implement the following code to apply the LPG() function to the
spiral, zigzag, and tree simulated data illustrated in Figures
5, 6, and 7.
## longitude and latitude are expressed in degrees
#### example 1: spiral data
> set.seed(40)
> n <- 900 # the number of samples
> sigma1 <- 1; sigma2 <- 2.5; # noise levels
> radius <- 73; slope <- pi/16 # radius and slope of the spiral
## polar coordinate of (longitude, latitude)
> r <- runif(n)^(2/3) * radius; theta <- -slope * r + 3
## transform to (longitude, latitude)
> correction <- (0.5 * r/radius + 0.3) # correction of noise level
> lon1 <- r * cos(theta) + correction * sigma1 * rnorm(n)
> lat1 <- r * sin(theta) + correction * sigma1 * rnorm(n)
> lon2 <- r * cos(theta) + correction * sigma2 * rnorm(n)
> lat2 <- r * sin(theta) + correction * sigma2 * rnorm(n)
> spiral1 <- cbind(lon1, lat1); spiral2 <- cbind(lon2, lat2)
## plot the spiral data
> rgl.sphgrid(col.lat = 'black', col.long = 'black')
> rgl.sphpoints(spiral1, radius = 1, col = 'blue', size = 12)
## implement the LPG to (noisy) spiral data
> LPG(spiral1, scale = 0.06, nu = 0.1, seed = 100)
> LPG(spiral2, scale = 0.12, nu = 0.1, seed = 100)
