Toroidal PCA via density ridges

Principal Component Analysis (PCA) is a well-known linear dimension-reduction technique designed for Euclidean data. In a wide spectrum of applied fields, however, it is common to observe multivariate circular data (also known as toroidal data), rendering spurious the use of PCA on it due to the periodicity of its support. This paper introduces Toroidal Ridge PCA (TR-PCA), a novel construction of PCA for bivariate circular data that leverages the concept of density ridges as a flexible first principal component analog. Two reference bivariate circular distributions, the bivariate sine von Mises and the bivariate wrapped Cauchy, are employed as the parametric distributional basis of TR-PCA. Efficient algorithms are presented to compute density ridges for these two distribution models. A complete PCA methodology adapted to toroidal data (including scores, variance decomposition, and resolution of edge cases) is introduced and implemented in the companion R package ridgetorus. The usefulness of TR-PCA is showcased with a novel case study involving the analysis of ocean currents on the coast of Santa Barbara.


Introduction
Toroidal data is formed by tuples of observations that lie on the p-dimensional torus T p = [−π, π) p , p ≥ 1, where −π and π are identified.This kind of data is present in a variety of applied fields ranging from bioinformatics, when modeling angles in atom bonds (e.g., Boomsma et al., 2008), to environmental sciences, when modeling sea conditions (e.g., Jona-Lasinio et al., 2012), among many others (see, e.g., Ley and Verdebout, 2018).The fact that circular variables present different properties than linear variables renders spurious the application of standard statistical tools, with their adaptations comprising the so-called directional statistics (Mardia and Jupp, 1999); see, e.g., Pewsey and García-Portugués (2021) for a recent review.In particular, Principal Component Analysis (PCA) is a widely used (linear) statistical tool whose extension to circular variables has been challenging due to the ill-posedness of linear relations on T p .On T2 , PCA produces linear subspaces that are either non-periodic or, if extended periodically by "wrapping" at the −π ≡ π boundaries, would almost surely result in useless perfect fits to a sample.Indeed, infinitely-dense line wrappings are produced by irrational slopes for the first sample principal component, which arise with probability one in a sample generated by a continuous random vector.
There have been several attempts to extend PCA to toroidal data, with all of them facing some relevant limitations or foundational issues; see Section 8.1 in Pewsey and García-Portugués (2021) for an overview.These approaches can be divided into two general branches.The first one consists of "geodesic-based" methods, like Fletcher et al. (2004)'s Principal Geodesic Analysis (PGA).PGA defines the first principal geodesic as the one passing through the intrinsic mean that minimizes the sum of squared intrinsic residuals.Nodehi et al. (2015) applied PGA to the torus.Both of these approaches inevitably may lead to the aforementioned infinite wrapping.The other branch of techniques is based on embedding toroidal data in Euclidean space to use classical PCA afterward.For example, Mu et al. (2005) proposed dihedral PCA (dPCA), which maps the angles onto S 1 × S 1 via (ϕ, ψ) → (cos ϕ, sin ϕ, cos ψ, sin ψ) and then applies regular PCA.Although this method transforms the data in a unique way, the geometry of S 1 × S 1 is ignored in the subsequent PCA application, leading to substantial transformation-induced artifacts in the resulting PCA scores.Riccardi et al. (2009) introduced angular PCA (aPCA): a centering of the data with the circular mean followed by an application of PCA that disregards periodicity.Another idea was that of characterizing the covariance matrix, which Kent and Mardia (2009) did on a wrapped normal model using trigonometric moments to facilitate PCA on it.To reduce the distortion generated by transformation-based approaches, Sittel et al. (2017) proposed a modification on aPCA consisting of a shift of the data so that the "periodic region" is located at the lowest-density area (determined by histogram approximation) and most of the data is located in the more "non-periodic region".Like aPCA, the periodicity of the resulting principal components and scores is not enforced.Eltzner et al. (2018) introduced Torus PCA (T-PCA) by mapping T d onto S d , a geometrically-benign space that, unlike geodesic PCA, produces non-dense linear subspaces.Principal Nested Spheres (PNS) (Jung et al., 2012) is then carried out on S d .Although a principled approach, the transformation in T-PCA is not invariant under permutations of variables, which may yield significantly different outcomes, and is prone to create artifacts in, e.g., datasets with marginal uniform-like distributions.More recently, Zoubouloglou et al. (2022) introduced Scaled Torus PCA (ST-PCA), which uses multidimensional scaling maps between T d and S d designed to alleviate the aforementioned drawbacks in T-PCA, and then applies PNS.However, these maps are computationally demanding and difficult to invert.Despite many partial advances, the development of a flexible and successful PCA on the torus that is free of significant drawbacks remains remarkably open.
Differently from the previous PCA approaches, the purpose of this article is to advance a flexible density-driven PCA on T 2 .A landmark in flexible dimension reduction on R p was the principal curves of Hastie and Stuetzle (1989).The linear subspaces spanned by PCA's principal directions are self-consistent when the data is normally distributed, i.e., the expectation of all the points projecting onto a point x ∈ R p of the subspace is also x.This fact motivates the definition of principal curves as those that are self-consistent, although their existence is not guaranteed.Another alternative definition of principal curves was introduced by Delicado (2001Delicado ( , 2003) ) in terms of the total variance, which has the advantage of a guaranteed existence.The total variance of a random variable is minimal when its hyperplane is orthogonal to the first principal component.A different approach to flexible PCA was given by Ozertem and Erdogmus (2011), who defined the concept of density ridges, sets of points that characterize the main features of the density, for locally defined principal curves and surfaces.Notably, Ozertem and Erdogmus (2011) introduced an estimation approach for density ridges based on the mean shift algorithm.Although the first notions of density ridges were already introduced by Hall et al. (1992), it has been in the last decade that the topic has attracted more attention.To review just a few contributions, Genovese et al. (2014) showed Hausdorff-consistency of the mean-shift estimated ridge sets, Chen et al. (2015) gave a Gaussian-process representation of the asymptotic distribution of the Hausdorff distances between estimated and theoretical ridges, and Qiao and Polonik (2016) established the pointwise asymptotic normality of the estimated ridges and the asymptotic distribution of its uniform deviation from theoretical ridges.
This article aims to generalize PCA to bivariate toroidal data by making use of density ridges.The contributions towards this goal are manifold.First, we provide efficient algorithms for computing periodic density ridges for two standard toroidal distributions, the Bivariate Sine von Mises (BSvM) (Singh et al., 2002) and the Bivariate Wrapped Cauchy (BWC) (Kato and Pewsey, 2015).These two algorithms substantially reduce the computational burden of iterative methods, based on integral curves, by (1) targeting the implicit equations that define the density ridge, (2) leveraging specific solutions of these implicit equations, and (3) exploiting the ridge symmetries induced by the BSvM and BWC densities.Second, we present Toroidal Ridge PCA (TR-PCA), which uses the connected component of the density ridge that passes through the mode(s) of the density as a flexible first principal component in the parametric setting of BSvM and BWC distributions.TR-PCA provides PCA-like scores through the obtention of signed distances along the ridge, signed projections onto the ridge, and Fréchet means within the ridge.These operations are facilitated by Fourier approximations of the ridge curve, which provide a fast analytical handle.Using likelihood theory, TR-PCA automatically handles edge cases arising in practice, such as diagonal ridges, and decides the best-fitting underlying distribution.Third, we empirically evidence the advantages of TR-PCA over aPCA, arguably the most readily implementable alternative, in a series of illustrative numerical examples.In particular, unlike transformation-based approaches, TR-PCA does not introduce distortions and, unlike other approaches like aPCA, yields scores that inherit the periodicity of the data.Fourth, the usefulness of TR-PCA is demonstrated with a novel application on the study of ocean currents at the coast of Santa Barbara, where the obtained principal ridge explains the currents' behavior successfully.A final contribution is the companion R package ridgetorus, which provides an implementation of TR-PCA and allows the end-to-end replicability of the data application.
The organization of the rest of this article is as follows: Sections 2-4 study the population case of density ridges, and their sample version appears in Section 5, with the definition of TR-PCA, and in Section 6.More precisely, Section 2 provides the definition and some useful properties of density ridges on R p , as well as the computation methods thereof.Section 3 is devoted to the computation and properties of the density ridge for the BSvM distribution on T 2 .Section 4 gives an analogous construction for the BWC distribution.Then, Section 5 introduces an effective parametrization of the connected component of the BSvM/BWC ridge and explains how to obtain scores from it in the presence of a sample, resulting in the definition of TR-PCA in Algorithm 1. Section 6 shows an application of TR-PCA on real bivariate data.Finally, a discussion of the main conclusions, alternatives, and limitations of the present work is given in Section 7.

Definition
Density ridges are higher-dimensional extensions of the concept of mode that are informative of the main features of a density.A mode is a local maximum of the density function f that, if f ∈ C 2 (R p ), p ≥ 2, can be characterized by a null gradient and negative Hessian eigenvalues, if the degenerate cases are excluded.The eigendecomposition of the Hessian of a density function f ∈ C 2 (R p ) evaluated at a point x ∈ R p is given by Hf (x) = U(x)Λ(x)(U(x)) ′ , where U(x) = (u 1 (x), . . ., u p (x)) is a matrix whose columns are the eigenvectors of Hf (x) and Λ(x) = diag(λ 1 (x), . . ., λ p (x)), λ 1 (x) ≥ . . .≥ λ p (x), contains their eigenvalues.Denoting U (p−1) (x) := (u 2 (x), . . ., u p (x)) and defining the projected gradient on {u 2 (x), . . ., u p (x)} as D (p−1) f (x) := U (p−1) (x)(U (p−1) (x)) ′ Df (x), where Df (x) stands for the (column vector) gradient, the density ridge is defined by Genovese et al. (2014) as follows.
Definition 1 (Density ridge; Genovese et al. 2014) Clearly, two cases imply that x ∈ R p satisfies D (p−1) f (x) = 0.The first is Df (x) = 0.In that case, if λ 2 (x), . . ., λ p (x) < 0, x is a maximum or a saddle point.The second is that Df (x) is perpendicular to U (p−1) (x)(U (p−1) (x)) ′ , that is, the gradient is parallel to u 1 (x).In this case, the directions of maximum ascent (gradient) and maximum signed curvature (u 1 (x)) coincide (see Figure 1), with "signed" emphasizing that the maximum is not in absolute value terms since λ 1 (x) ∈ R.

Some useful properties
Two useful properties of density ridges are given below.Their proofs are given in Appendix A.
Figure 1: From left to right, representation of the eigenvector fields u 1 , u 2 : R 2 → R 2 , and the gradient and projected gradient vector fields Df, D 1 f : R 2 → R 2 for a Gaussian density f on R 2 .Red plus (respectively, blue minus) indicate a positive (negative) eigenvalue λ j (x) associated with u j (x), for j = 1, 2, with λ 1 (x) ≥ λ 2 (x) and x ∈ R 2 .For visualization purposes, vector fields are unit-norm standardized and low-density regions are shown in white.
Proposition 1 (Ridge invariance to translations and rotations).Let X be a p-random vector with density f ∈ C 2 (R p ), p ≥ 2. Let Y = µ + RX represent a shifting and rotation family of p-random vectors spanned by X, where µ ∈ R p and R is a p × p orthogonal matrix.If f (•; µ, R) stands for the density function of Y and x ∈ R p , then This result yields two convenient handles to manipulate density ridges.First, it reduces the problem of computing density ridges to those that are centered at a certain origin.Second, it shows how to exploit the symmetries of f to reduce the computational costs of evaluating R(f ) (e.g., halve them if f is reflective symmetric in R 2 ).
Proposition 2 (Ridges for elliptically-symmetric densities).Let f ∈ C 2 (R p ), p ≥ 2, be an ellipticallysymmetric density of the form f (x) = g((x − µ) ′ Σ −1 (x − µ)), with x, µ ∈ R p , Σ a p × p positive definite matrix, and g ∈ C 2 (R + 0 ) a (strictly) decreasing function.Then, the subspace spanned by the first eigenvector v 1 of Σ belongs to the density ridge of f : (1) This proposition has several important consequences.First, it proves that R(f ) is intimately related to the first principal component of PCA for elliptically-symmetric densities, since the subspace generated by the latter direction is included in it.A decreasing g implies that f has µ as a unique mode.This is satisfied by many elliptically-symmetric distributions, with its prime representative being the normal distribution, a class of distributions in which PCA is particularly meaningful.Second, like the first principal component, the density ridge contains the center (µ) in ellipticallysymmetric distributions.Third, the proposition shows that the density ridge is more than just the first principal direction subspace.This observation is key for defining the µ-connected density ridge for a mode or mean µ of f with a view to constructing a density-driven first principal component analog.
We return to the µ-connected density ridges in Sections 3 and 4. Before, we describe in the following subsections two approaches to determining in practice the set of points that conform R(f ) for f ∈ C 2 (R p ).Both approaches are exploited in Sections 3 and 4.

Integral curve approach
Genovese et al. (2014) showed that the vector field of the projected gradient defines a global flow.Hence, the trajectory defined by the projected gradient converges to the points where it is null.This fact allows characterizing the density ridge as an integral curve: and ϕ x 0 (0) = x 0 .This differential equation can be solved numerically using the Euler method until attaining a point where D (p−1) f (x) ≈ 0. To move faster in low-density regions and slower in high-density regions, it is useful to consider the normalized projected gradient η (p−1) (x) := D (p−1) f (x)/f (x): The above scheme is started from an initial point x 0 and is iterated using a step h > 0 until a criterion for convergence is met after N iterations.The resulting x N (approximately) belongs to R(f ).Then, it is possible to approximately determine R(f ) by running (2) from a sufficiently dense grid of initial values x 0 .

Implicit equation approach
When p = 2, from Definition 1 it can be seen that R(f , with the subscript indicating the first/second component of the gradient and the first/second component of the second eigenvector u 2 (x) of the Hessian at x. Furthermore, still when p = 2, the eigenvectors and eigenvalues of the Hessian admit a closed form.For example, Qiao and Polonik (2016) give: . This means that the implicit equation in terms of the derivatives of f is given by and that the eigenvalue condition reads as Therefore, it is possible to obtain R(f ) for a bivariate density f by solving Equation (3), restricted to Equation ( 4), along a grid on one of its variables.This approach is much faster than that in Section 2.3, yet it is restricted to p = 2.
3 Density ridges for bivariate sine von Mises

Bivariate sine von Mises
Let Θ 1 and Θ 2 be two circular random variables.Singh et al. (2002)'s Bivariate Sine von Mises (BSvM) distribution has density given by , where I m is the modified Bessel function of order m.The distribution is pointwise symmetric about (µ 1 , µ 2 ), with these two parameters representing the marginal circular means.The parameters κ 1 and κ 2 measure the marginal concentrations of the distribution about µ 1 and µ 2 , respectively.The parameter λ measures dependence: positive/negative values of λ correspond to positive/negative correlation between Θ 1 and Θ 2 .If λ = 0, then Θ 1 and Θ 2 are independent, with each variable having a (univariate) von Mises distribution.Density (5) can be bimodal.A sufficient unimodality condition is given by κ 1 κ 2 > λ 2 (Mardia et al., 2007).In the bimodal case where µ 1 , µ 2 ∈ T 2 are the two modes of ( 5), due to the symmetry of the BSvM, it is easy to see that R µ 1 (f BSvM ) = R µ 2 (f BSvM ), thus ensuring that a connected component is well defined.
The BSvM has some properties that make it a candidate toroidal analog to the bivariate normal.Among these are the facts that it is also part of the exponential family and that it has asymptotic bivariate normal distributions under high concentrations (Singh et al., 2002).Furthermore, the marginals of the BSvM are not von Mises, but the conditionals are. Figure 2 shows different forms of the BSvM density.For simplicity, and without loss of generality due to Proposition 1, we will work with (µ 1 , µ 2 ) = (0, 0).
No closed expressions for the maximum likelihood estimators of the BSvM parameters exist (Mardia et al., 2008), but these can be computed numerically using moment-based estimators as starting values in the optimization of the log-likelihood.In particular, this enables testing the homogeneity hypothesis H 0 : κ 1 = κ 2 vs. H 1 : κ 1 ̸ = κ 2 using the Likelihood Ratio Test (LRT) that asymptotically rejects H 0 at the significance level α whenever −2( l − l0 ) > χ 2 α;1 , where l and l0 are the maximum likelihoods of (5), the latter under H 0 , and χ 2 α;1 is the α-upper quantile of the chi-square distribution with one degree of freedom.This homogeneity LRT is convenient for distinguishing in practice the edge case κ 1 = κ 2 , which has a simple diagonal ridge associated.An analogous LRT can be constructed for testing H 0 : λ = 0 vs. H 1 : λ ̸ = 0, a null hypothesis of independence associated with horizontal/vertical ridges if κ 1 is smaller/larger than κ 2 .

Integral curve approach
The first alternative to compute the R(f BSvM ) of a BSvM(0, 0, κ 1 , κ 2 , λ) uses the integral curve approach.With the Euler algorithm, an initial grid of points converges to the density ridge by following the trajectory defined by the projected gradient, as Figure 2 shows.It can be seen that the main part of the density ridge captures the main features of the distribution, shape, and correlation, while also being periodic.The figure also shows that there is a significant number of "secondary" ridges.

Implicit ridge equation approach and connected ridges
The implicit equation offers a faster alternative to the computationally-demanding Euler algorithm.In the following, the normalizing constant of the BSvM and a common positive factor will be ignored, as these do not affect the direction of η (p−1) .It is easy to check that D These expressions are readily usable in (3)-(4).Figure 3 shows the density ridges obtained with this method, which allow a clear identification of R 0 (f BSvM ).
However, in general, there is no explicit parametrization of R µ (f BSvM ) with µ ∈ T 2 , so a method is needed to filter it from the full R(f BSvM ) obtained from the implicit equation.The symmetry of f BSvM about µ and Proposition 1 imply (see Figure 3 for graphical insight) that: (1) to the first or second quadrant, depending on sign(λ).Finally, since R µ (f BSvM ) passes through µ by definition, a connected component can be obtained by iteratively adding the sufficiently close ridge points to the last point added to the set, starting from µ.The following result presents limit cases for which the density ridge R 0 (f BSvM ) is explicit.
Proposition 3. The 0-connected density ridge R 0 (f BSvM ) of a BSvM(0, 0, κ 1 , κ 2 , λ) admits the following representations: Remark 1.Note that (ii) does not characterize R 0 , unlike (i).Although in view of the third leftmost panel in Figure 3 it seems possible to achieve this characterization with an additional restriction in Definition 2, we do not pursue this further due to its limited practical interest: given will not always wrap periodically at −π ≡ π.To solve this edge-case issue in practice, and in coherence with the other cases, we simply extend R 0 (f BSvM ) according to the diagonal.

Bivariate wrapped Cauchy
The density of a random vector (Θ 1 , Θ 2 ) following the Bivariate Wrapped Cauchy (BWC) distribution proposed by Kato and Pewsey (2015) is The several constants are given as follows: . Analogously to the BSvM, µ 1 and µ 2 represent the marginal circular means of the density.The parameters ξ 1 and ξ 2 regulate the concentrations of the marginal distributions, that of Θ j being circular uniform when ξ j = 0 (j = 1, 2).When ξ 1 , ξ 2 > 0, the density is unimodal and pointwise symmetric about (µ 1 , µ 2 ).As ξ j → 1, the marginal distribution of Θ j tends to a point mass at µ j .The association between Θ 1 and Θ 2 is controlled by the parameter ρ, ρ = 0 corresponding to independence.Positive/negative values of ρ correspond to positive/negative correlation between Θ 1 and Θ 2 .Figure 4 shows different forms of the BWC density.This distribution is always unimodal (Kato and Pewsey, 2015), thus ensuring that R µ is well-defined.
The BWC is closed under marginalization and conditioning, meaning that the resulting marginals and conditional densities belong to the wrapped Cauchy family.This distinctive closedness property is shared with the bivariate normal distribution, whose marginals and conditionals are also normal.The functional form of marginals is immediate from the BWC's construction as a particular Wehrly and Johnson (1980) model, but the conditionals are not (despite being wrapped Cauchys).
As in the BSvM case, there are no closed expressions for the maximum likelihood estimators for the BWC parameters (Kato and Jones, 2015), so numerical routines are needed.Analogously to the BSvM case, the LRTs for H 0 : ξ 1 = ξ 2 vs. H 1 : ξ 1 ̸ = ξ 2 (homogeneity) and H 0 : ρ = 0 vs. H 1 : ρ ̸ = 0 (independence) are also relevant to distinguish diagonal and horizontal/vertical ridges in a practical and principled manner.
The integral curve approach for the BWC is completely analogous to the BSvM case.

Implicit ridge equation approach and connected ridges
It is simple to check that the derivatives of the BWC density, excluding a common positive factor, are , where f * = f /c. Figure 4 shows the density ridges obtained with the implicit equation approach.This figure illustrates that, differently to the sinusoidal shapes of R 0 (f BSvM ), the shapes of R 0 (f BWC ) involve ridges that connect (0, 0) with (±π, ±π) for unequal marginal concentrations.

Ridge parametrization
Due to the difficulty in explicitly solving the implicit equation (3), we have not found any simple parametric form for the curve defined by R µ (f ) beyond limit cases.This poses an important hindrance to the tractability of the distance operation along R µ (f ) and the projection operation that maps a point on T 2 to its closest point on R µ (f ), both being core mechanisms for the definition of the forthcoming scores.To prevent this issue from draining the performance of TR-PCA, we consider a Fourier-type parametrization of R µ (f ) given by , depending on which is the index coordinate j (see the next paragraph).The consideration of only the cosine/sine parts in ( 6) is justified by the pointwise symmetry of R 0 (f BSvM ) and R 0 (f BWC ), which renders null Fourier cosine/sine coefficients for sine/cosine components.In practice, {(a k , b k )} m k=0 is approximated by Gaussian quadrature using a grid of points in R µ (f ) determined using the implicit equation method.In numerical experiments, the truncation of the Fourier series (6) to m = 15 was found to be a sensible choice for a large number of parameter specifications of the BSvM and BWC densities.With m = 15, the maximum distance between the implicit-equation computation of R µ (f ) and the Fourier-parametrized approximation was found to be smaller than 10 −2 .Therefore, we consider m = 15 to Fourier-parametrize R µ (f BSvM ) and R µ (f BWC ).
To obtain a one-to-one parametrization for the BSvM and BWC densities, in practice r f,µ,j is indexed along the variable with the smallest concentration (e.g., θ 2 and θ 1 for the first and second leftmost panels in Figure 2, respectively), which is straightforward to identify.For the sake of notational simplicity, henceforth we assume without loss of generality that j = 1.
To define the first scores, it is fundamental to parametrize the curve r f,µ,1 through its arc length with id denoting the identity function.Two further tweaks are required on r f,µ to: (1) scale the parametrization to s ∈ [0, 2π); (2) center the parametrization to s ∈ [−π, π).The first step sets the range of the scores to 2π, just as the original data in T, while the second step is crucial for assigning signs.The definition below summarizes the construction.

Scores
Projections on R µ (f ) are defined via the fast handle rf,µ .They involve the toroidal distance Definition 4 (Ridge projections).For θ ∈ T 2 , its projection to the Fourier-parametrized where the projection argument is The TR-PCA scores for an arbitrary point θ ∈ T 2 are defined as an analogy to ordinary PCA.The first score is the signed distance along rf,µ and between proj f,µ (θ) and µ.The second score is the signed distance between θ and proj f,µ (θ).The sign is set according to the relative position of the tangent and normal vectors, both at the projection point.

Proportion of variance explained
Given a sample Θ 1 , . . ., Θ n in T p , n, p ≥ 1, its Fréchet mean (or intrinsic mean) is defined as The Fréchet variance of the sample is the minimum of the previous objective function: F stands for the Fréchet variance of the jth scores {s j (Θ i )} n i=1 , j = 1, 2.

Complete TR-PCA procedure
The complete TR-PCA procedure for a given toroidal sample involves all the concepts introduced so far.It is divided into three main stages.
The first scenario (leftmost column) shows the almost equivalence between both methods when dealing with highly-concentrated samples, an expected consequence of the torus being locally Euclidean.Here, both the first scores of TR-PCA and aPCA explain 81% of the variance.When the distribution is more dispersed, as in the other three samples, periodicity becomes relevant and aPCA begins introducing artifacts.In the second scenario, TR-PCA (83%) does not introduce any distortion on the scores thanks to honoring the toroidal geometry, while aPCA (75%) creates some spurious clusters with scores that exit T 2 (observe the vertical scale) and induces a slight rotation on the scores.These artifacts are magnified in the third scenario (TR-PCA: 88%; aPCA: 74%).Finally, the fourth scenario represents a Simpson's paradox case in which TR-PCA (78%) is able to successfully separate the two clusters along the first scores, while aPCA (38%) fails to do so.In all scenarios, TR-PCA yields periodic scores, unlike aPCA.The PVEs of both methods were computed according to (10).

Data application
This section illustrates the application of TR-PCA to the study of currents in the Santa Barbara Channel (California, USA).Studying the direction of currents is of crucial importance to understand the supply of nutrients in marine habitats (Allen et al., 2012) and the genetic propagation of marine fauna and flora (White et al., 2010), as well as is important for environmental purposes and prevention of contamination (DiGiacomo et al., 2004).In the present context of increased contamination and climate change, the Santa Barbara Channel is an area well known for the confluence of several important ocean currents and vast marine biodiversity.The complexity of local currents (Winant et al., 2003) makes the use of standard statistical tools not directly applicable and motivates the search for new methodologies able to explain their variability.In their Figure 10, Winant et al. (2003) explain that there exists a counterclockwise vortex in the Santa Barbara Channel, which is represented by a general westward flux on its northern coast.In addition, studies such as Auad et al. (1998) show that there exists a net influx of water in the Santa Barbara Channel that exits heading south through the Santa Cruz Channel.These two facts are taken as a guideline to validate the ability of TR-PCA on indexing the main variability of the data in a fully data-driven way.
The data was obtained from the High-Frequency Radar Database (https://hfradar.ndbc.noaa.gov/), which maps surface currents and wave fields over wide areas.In particular, the currents' direction, d = atan2(u, v), is hourly available through the measured eastward and northward surface velocities (u and v, respectively) of the water body.We smoothed this direction by taking the daily speed-weighted circular mean in a given zone Z, obtaining the circular variable of the study, θ Z .Since large-scale currents involve timespans longer than hours, daily averages allow smoothing of the data while still keeping the variability associated with these ocean currents.We used the data from the three-year period 2019-2021, using complete years so that seasonality events were neither over-nor underrepresented.Four locations were selected: zones A and B, located along the northern coast of the Santa Barbara Channel; and zones C and D, corresponding to the north and south of the Santa Cruz Channel, respectively.These areas are shown in Figure 6.The study focuses on the dependency between θ A and θ B to further analyze the top part of the vortex shown in (Winant et al., 2003, Figure 10) and the water flux parallel to the coast, as well as on the dependency between θ C and θ D for the water output through the inter-island strait (Auad et al., 1998).The third analysis on (θ A , θ C ) is performed to investigate the dependence between both regions.A first analysis of the data's density in the three zones (A-B, A-C, C-D) was performed to check the adequacy of the BSvM and BWC distributions to model the current directions.Toroidal-adapted kernel density estimators for estimating the unknown bivariate densities were compared with the estimated parametric densities of the BSvM and BWC.These estimates are shown in Figure 7, where it can be seen that the BWC and BSvM give similar shapes, although simplifying some of the asymmetries of the kernel density estimators.To formally assess the existence of dependence in the four scenarios, we performed the ϕ (n) dc test described in García-Portugués et al. (2023, Section 3), emphatically rejecting the null hypothesis of independence in all cases.
The BSvM and BWC models were fitted to the samples of (θ A , θ B ), (θ A , θ C ), and (θ C , θ D ) using maximum likelihood.The estimated parameters and log-likelihoods are shown in Table 1.In terms of log-likelihoods, it can be seen that the BSvM estimate has better performance in two of the analyses.As expected, the BWC distribution is more concentrated, leading to higher values of the density around the modes.The estimated ridges R μ( f ) for the three bivariate analyses are shown in Figure 7.It can be seen how R μ( f ) indexes the main variability of the sample and informs on the positioning of high-density regions.On the one hand, the best-fitting models capture the local correlations of the currents' directions about their modes, which are seen to be positive in the A-B case, and negative in the A-C and C-D cases.On the other hand, the curve R μ( f ) allows synthesizing the behavior of the currents in a sensible way.For example, in the A-B case, the BSvM model not only reproduces the W-W mode (western direction both at zone A and B) found in previous studies (Winant et al., 2003), but also provides additional information in terms of the sinusoidal-like dependence between A-B.A final insight revealed by the estimated ridges is that there is a high negative correlation in C-D.Due to the geographical shape of the strait, water can flow through it in two main directions: SE-SE and NW-NW, as previous studies show (Auad et al., 1998).As seen in Figure 7, both main modes are captured fairly well: there is a high concentration toward SE-S and a secondary concentration at NW-W.Furthermore, it can also be seen that the higher-density transition region between the two modes has a negative correlation, a result that cannot be easily obtained by only analyzing the modes.
The summarizing capability of R μ( f ) can be further exploited by visualizing a march along it.Figure 8 shows this march in the case C-D using the BSvM fit shown previously in Figure 7. Unlike previous studies that were limited to finding the main modes of direction (Auad et al., 1998;Winant et al., 2003), this parametrization also allows an estimation of the variability when the flow has other directions.For example, the net influx of water outside the strait is reproduced in the second leftmost plot in Figure 8.
TR-PCA scores were computed in the rightmost column of Figure 7, resulting in narrower, more concentrated distributions that capture a large part of the total variance, with the PVEs collected in Table 1 being greater than 74%.These scores could be used for clustering purposes to find particular meteorological events or outliers.

Discussion and conclusions
In this work we have advocated the use of density ridges for constructing a well-grounded bivariate toroidal PCA.The construction is based on using the implicit equation approach to determine ridges, which has proven to be more efficient and robust than the Euler algorithm.By tailoring this procedure to bivariate data, we have corroborated empirically that two reference toroidal distributions, the BWC and the BSvM, present stable connected components that go through the distributions' modes.We have proposed a practical way to parametrize R µ (f ) to yield a tractable computation of PCA-like scores that allows for a full dimension-reduction analysis.The real data application has shown that TR-PCA explains 75% − 80% of the variability of the ocean currents and gives interpretations that are consistent with previous explanations of large-scale water movements in the area.
An important takeaway of our investigations is that the BWC seems to be, in general, a more robust parametric distribution than the BSvM for TR-PCA.Although the BSvM is recognized as a somewhat canonical choice for "the normal density on T 2 ", for this application, the squarish form of its density contours may introduce "elbows" on its associated density ridges, these in turn introducing artifacts on the resulting scores.In comparison, the BWC does not present this problem and tends to yield more flexible and descriptive principal ridges.Nevertheless, the choice of the reference parametric distribution is subject to improvement, as any sufficiently well-behaved density in C 2 (T 2 ) could be considered within our methodology.In this regard, we highlight that the ridge parametrization introduced in Section 5.1 could be used in other densities apart from BSvM/BWC.Deriving theoretical results giving the conditions under which the existence and well-definedness of a connected ridge for parametric families of densities would be a useful endeavor to be addressed in the future.In this respect, we experienced with another normal analog on the torus, the wrapped normal, for the analyses in Sections 3.1, 3.2, and 5.5.We found numerically that some ridges could be multivalued functions in both variables θ 1 and θ 2 , giving convoluted ridge curves.When this does not happen, somehow unpleasant Z-shaped ridges with rough corners arise, thus making this distribution less appropriate for TR-PCA than BSvM and BWC.Although not pursued in this paper, we note that asymptotic inference for R µ (f ) is readily addressable when using a density model f with tractable maximum likelihood.
TR-PCA has been constructed as a toroidal analog of PCA on R 2 under the following optic.PCA can be seen as a parametric dimension-reduction method driven by Gaussianity: in an arbitrary sample, it fits a normal distribution N (µ, Σ) by maximum likelihood and extracts the eigendecomposition of Σ; the subspace spanned by the first eigenvector coincides with the μ-connected ridge.TR-PCA follows this view of PCA, replacing the normal distribution with the BSvM/BWC model.It might be argued that PCA is a model-free technique since ( μ, Σ) estimate the mean and covariance matrix of a general population.In the torus, these two descriptive summaries are less canonical; e.g., there exist two definitions of circular means (extrinsic or intrinsic).The maximum likelihood estimators of the BSvM and BWC distributions do not necessarily coincide with extrinsic/intrinsic means, and so TR-PCA does not mimic this Gaussian-specific aspect of PCA.However, the maximum likelihood estimators of the BSvM/BWC model f ξ estimate ξ 0 = arg min ξ D KL (f 0 ∥f ξ ), the parameter that minimizes the Kullback-Leibler divergence of a general population f 0 from f ξ .In the Gaussian model f µ,Σ , (µ 0 , Σ 0 ) = arg min µ,Σ D KL (f 0 ∥f µ,Σ ) are the mean and covariance matrix of the population f 0 , highlighting that both TR-PCA and PCA estimate a general population descriptor under a parametric model.
A clear limitation of the present work is its restriction to p = 2.The concept of ridges can be extended to higher dimensions by redefining the eigenvector matrix, the projected gradient, and the density ridge of a multivariate toroidal density f ∈ C 2 (R p ) as U (p−q) (x) := (u q+1 (x), . . ., u p (x)), D (p−q) f (x) := U (p−q) (x)(U (p−q) (x)) ′ Df (x), and R q (f ) := {x ∈ R p : ∥D (p−q) f (x)∥ = 0, λ q+1 (x), . . ., λ p (x) < 0}, with 1 ≤ q ≤ p − 1.However, in this multivariate case, obtaining an equivalent of the implicit equation is more challenging, so one may need to rely on the (computationally expensive) Euler algorithm.Still, important concepts such as iterative projections on nested spaces and scores computation would need to be carefully defined for the multivariate case.The multivariate extension of TR-PCA is also currently limited by the scarcity of multivariate toroidal models beyond the multivariate von Mises distribution (Mardia et al., 2008) (which extends (5)).For example, currently there does not exist a multivariate extension of the BWC.Therefore, the larger endeavor of extending TR-PCA to higher dimensions is open for future research.
The statement (ii) follows due to the arc-length parametrization of (7) and the scaling in (8).

Figure 3 :
Figure 3: Catalog of the implicit-equation-approximated ridges R(f BSvM ) (black) and R 0 (f BSvM ) (red) of a zero-centered BSvM for the same parameter values as in Figure 2. The background shows the density contour of the BSvM.All the contourplots share the same color scale.

F
Due to the product structure of T p , , where the superscript denotes a marginal Fréchet mean/ variance.For p = 2, this variance decomposition facilitates the definition of the Proportion of Variance Explained (PVE) in TR- (a) Fit the BSvM and/or BWC models (Sections 3.1 and 4.1) with maximum likelihood estimation.If both models are fit, select the one with the smallest Bayesian Information Criterion (BIC).(b) Inspect edge cases using LRTs (Sections 3.1 and 4.1) at 5% significance level.(b.i) Test diagonal vs. non-diagonal ridges with the independence LRT.(b.ii) Test straight vs. non-straight ridges with the homogeneity LRT.(c) If any of the LRTs does not reject, refit the model with maximum likelihood estimation restricted to the decisions of (b.i)-(b.ii).(d) Retrieve f , μ (location parameter), and ĵ (index of the lowest concentration).(ii) Ridge computation.(a) Determine a grid of R 0 ( f ) with the implicit equation approach (Sections 3.3 and 4.2).(b) Construct r f , θ, ĵ in (6) with {(â k , bk )} m k=0 computed with Gauss-Legendre quadrature on the previous grid.(c) Obtain the arc-length parametrized ridge curve r f , θ from Definition 3. (iii) Scores and PVE computation.

Figure 7 :
Figure 7: From left to right, the columns represent the density contours of kernel density estimates, fitted BSvM densities, fitted BWC densities, and kernel density estimates of the TR-PCA scores of the bestperforming model.The second and third columns show the estimated ridge curves R μ( f ) (red curves) from the sample (black points in the first column) for the four bivariate analyses and two models.From top to bottom, rows represent zones A-B, A-C, and C-D.The contourplots within the same row share a common color scale in the first three columns.The parameters of the fits are shown in Table1.

Figure 8 :
Figure 8: Four snapshots of the march along the ridge curve R μ( f ) for C-D.The arrows at C (top) and D (bottom) are colored according to their position on the ridge curve.Past directions are also shown with transparencies in order to visualize the movement's direction.The main variability of the currents on C-D is manifested on a pendular-like variation of the current direction at D that is aligned with a counterclockwise variation of the current direction at C. Video available at https://github.com/egarpor/ridgetorus/tree/main/application.

Table 1 :
Estimated parameters, log-likelihood ( l), and Proportion of Variance Explained (PVE) by TR-PCA for the different analyses.Left and right blocks give the BSvM and BWC fits, respectively.The bold font indicates the best-performing model in terms of the log-likelihood.