Probabilistic Geometric Principal Component Analysis with application to neural data
Authors: Han-Lin Hsieh
Maryam Shanechi
Venue: ICLR 2025 Spotlight (2025)
Keywords: geometry
nonlinear manifold
factor analysis
dimensionality reduction
neural population activity
Primary area: applications to neuroscience & cognitive science
TL;DR: PGPCA generalizes PPCA by incorporating any nonlinear manifold with various distribution coordinates in its probabilistic model for high-dimensional brain data that exhibit noise.

Abstract: Dimensionality reduction is critical across various domains of science including neuroscience.  Probabilistic Principal Component Analysis (PPCA) is a prominent dimensionality reduction method that provides a probabilistic approach unlike the deterministic approach of PCA and serves as a connection between PCA and Factor Analysis (FA). Despite their power, PPCA and its extensions are mainly based on linear models and can only describe the data in a Euclidean coordinate system around the mean of data. However, in many neuroscience applications, data may be distributed around a nonlinear geometry (i.e., manifold) rather than lying in the Euclidean space around the mean. We develop Probabilistic Geometric Principal Component Analysis (PGPCA) for such datasets as a new dimensionality reduction algorithm that can explicitly incorporate knowledge about a given nonlinear manifold that is first fitted from these data. Further, we show how in addition to the Euclidean coordinate system, a geometric coordinate system can be derived for the manifold to capture the deviations of data from the manifold and noise. We also derive a data-driven EM algorithm for learning the PGPCA model parameters. As such, PGPCA generalizes PPCA to better describe data distributions by incorporating a nonlinear manifold geometry. In simulations and brain data analyses, we show that PGPCA can effectively model the data distribution around various given manifolds and outperforms PPCA for such data. Moreover, PGPCA provides the capability to test whether the new geometric coordinate system better describes the data than the Euclidean one. Finally, PGPCA can perform dimensionality reduction and learn the data distribution both around and on the manifold. These capabilities make PGPCA valuable for enhancing the efficacy of dimensionality reduction for analysis of high-dimensional data that exhibit noise and are distributed around a nonlinear manifold, especially for neural data.

Published as a conference paper at ICLR 2025
PROBABILISTIC GEOMETRIC PRINCIPAL COMPONENT ANALYSIS WITH APPLICATION TO NEURAL DATA
Han-Lin Hsieh
Ming Hsich Department of Electrical and Computer Engineering Viterbi School of Engineering University of Southern California
Los Angeles, CA, U.S.A hanlinhs@usc.edu
Maryam M. Shanechi
Ming Hsieh Department of Electrical and Computer Engineering Thomas Lord Department of Computer Science
Alfred E. Mann Department of Biomedical Engineering
Viterbi School of Engineering
University of Southern California
Los Angeles, CA, U.S.A shanechi cerecchi@uceeeii
ABSTRACT
Dimensionality reduction is critical across various domains of science including neuroscience, Probabilistic Principal Component Analysis (PPCA) a prominent dimensionality reduction method that provides probabilistic approach unlike the deterministic approach of PCA and serves as connection between PCA and Factor Analysis (FA). Despite their power, PPCA and its extensions are mainly based on linear models and can only describe the data in a Euclidean coordinate system around the mean of data. However, in many neuroscience applications, data may be distributed around nonlinear geometry (i.e., manifold) rather than lying in the Euclidean space around the mean. develop Probabilistic Geometrie Principal Component Analysis (PGPCA) for such datasets new dimensionality reduc- tion algorithm that can explicitly incorporate knowledge about a given nonlinear manifold that is first fitted from these data. Further, show how in addition to the Euclidean coordinate system, a geometric coordinate system can be derived for the manifold to capture the deviations of data from the manifold and noise. We also derive data-driven EM algorithm for learning the PGPCA model parameters. As such, PGPCA generalizes PPCA better describe data distributions by incorpo- rating a nonlinear manifold geometry. In simulations and brain data analyses, we show that PGPCA can effectively model the data distribution around various given manifolds and outperforms PPCA for such data. Moreover, PGPCA provides the capability to test whether the new geometric coordinate system better descrihes the data than the Euclidean one. Finally, PGPCA can perform dimensionality reduction and learn the data distribution both around and on the manifold. These capabilities make PGPCA valuable for enhancing the efficacy of dimensionality reduction for analysis of high-ding digh-dimensional data that exhibit noise and are distributed around a nonlinear manifold, especially for neural data.
INTRODUCTION
There exist numerous well-established algorithms for dimensionality reduction designed to efficiently
identify principal components that explain crucial features in high-dimensional (high-D) data in IR...
Distinct features necessitate different algorithms. Among them, Principal Component Analysis (PCA) (Greenacre et al., 2022) and maximum likelihood Factor Analysis (FA) (Bartholomew et al., 2011) are
"Corresponding Author.

---

Published as a conference paper at ICLR 2025
widely recognized. PCA, grounded in deterministic model, maximizes a critical feature-th data variance- explainedl by its principal components. In contrast, FA, rooted in a probabilistic model, efficiently captures another important feature the correlation between elements in the data-via its loading matrix (analogous to principal components in PCA) by maximizing the data log- ogglikelihood. This fundamental difference results in PCA and FA being employed in deterministic and probabilistic analyses separately.
The gap stemming from the fundamental difference in the model assumptions of PCA and FA is bridged by Probabilistic Principal Component Analysis (PPCA) (Tipping & Bishop, 1999b). PPCA, a special type of FA with a uniform diagonal noise covariance, proves that this new condition renders FA's loading matrix equal to PCA's principal components (Tipping & Bishop, 1999b). This insight enhances PCA's value both theoretically and practically. Now, PCA's principal components maximize not only the data variance but also the data log-likelihood given an underlying FA probabilistic model. Moreover, these principal components (equivalent to the loading matrix in PPCA) can be computed analytically, unlike FA's loading matrices, which generally require numerical solutions (De Winter & Dodou, 2012).
All of these methods are widely used in real-world applications including in neuroscience for analyses of neural population activity. However, PCA, FA, and PPCA are all grounded in linear models composed of Euclidean coordinate system around the mean of data. As such, they do not capture nonlinear structure in data. Further, this linear assumption explains why these algorithms do not require the selection of coordinate system, all linear bases of R are equivalent under linear transformations, However, in many applications including in neuroscience, data may be distributed around a nonlinear manifold rather than lying in the Euclidean space around the mean. For example, neural population activity has been shown to be distributed around a ring manifold in the head direction system (Chaudhuri et al., 2019; Jensen et al., 2020) or around a Torus manifold in the hippocampus (Gardner et al., 2022). Indeed, knowledge about the manifold can be provided through various existing methods based on data. For example, the type of manifold underlying noisy data can be identified by the topological data analysis (TDA) (Singh et al., 2008), and that manifold can be fitted by splines (Zheng ct al., 2012; He & Shi, 1996; Bojanov et al., 2013) or other graph-base based methods (Yin et al., 2008; Fefferman et al., 2018; 2023). However, incorporating such knowledge of a nonlinear manifold that is first fitted from data within the PPCA framework remains challenging to date.
While a few studies have explored nonlinear extensions of PPCA, they have assumed that data must lie precisely on top of specific manifold without deviations from it (Zhang & Fletcher, 2013; Zhang et al., 2019; Nodehi et al., 2020). However, in many applications including in neuroscience, rather than being on top of a manifold, data is distributed around it with some deviation and with noise. Thus, these prior PPCA extensions have not found application for these datasets such as neural activity. To model these datasets, we need to find a way to not only incorporate a given manifold, but also derive a coordinate system which we term distribution coordinates to capture the deviation outside of the manifold. Indeed, it may be possible to compute a coordinate systems attached to this nonlinear manifold that is not equivalent to the Euclidean coordinate system and that better describes the data; this is different from the case of the linear model in PPCA in which all coordinate systems are equivalent under linear transformations.
Here we address the above challenges by introducing Probabilistic Geometric Principal Component Analysis (PGPCA). PGPCA generalizes PPCA. Given a nonlinear manifold that is first fitted from data, PGPCA can incorporate this manifold with distribution coordinates that are computed for this manifold in its probabilistic model. PGPCA achieves dimensionality reduction by maximizing the data tog-likelihood Due to the nonlinear manifold, the Singular Value Decomposition (SVD) used in PPCA cannot be used to find the loading matrix in PGPCA. Thus, we Maximizatiioin (EM) algorithm to compute the PGPCA loading matrix. Further,
Contributions
derive an we show how in addition to the Euclidean distribution coordinate, a geometric distribution coordinate can be derived for the manifold to capture the deviations of data from the manifold and noise. Due to the nonlinear manifold/geometry, the geometric and Euclidean distribution coordinates yield different data log-like likelihood values. As such, we show how we can compute this log- ligeeiildoood and use it as a metric for distinguishing the distribution coordinates in a data- manner.
We structure this paper as follows. In Section 3, we first provide detailed mathematical derivation of PGPCA, including its prohabilistic model and the corresponding EM learning algorithm. In Section
2

---

Published as a conference paper at ICLR 2025
4, demonstrate the success of PGPCA with simulations of multiple manifolds and analyses neural population data from the mouse head direction system (Peyrache et al., 2015; Chaudhuri et al., 2019). We also show that PGPCA outperforms the existing PPCA framework by capturing the geometry in both simulations and real data. Finally, we illustrate PGPCA's ability to distinguish between geometric and Euclidean distribution coordinates in simulations and real data, In Section 5, we present a summary and discuss limitations.
2 RELATED WORK
Various extensions have been developed based on PPCA. Ahn & Oh (2003) modifies the PPCA EM algorithm to more efficiently compute the PCA principal components in order. To improve the interpretation of PPCA, prior studies have made its loading matrix sparse by, for example, restricting the domain of the probabilistic distribution in the E-step of PPCA EM (Khanna et al., 2015) or adding penalty terms in the cost function (Park et al., 2017). Penalizing the PPCA EM cost function has also been used in finding the efficient PPCA model dimension (Deng & Craiu, 2023). A supervised version of PPCA (Yu et al., 2006) has also been developed for labeled data. Zhang et al. (2017) has focused on using the mixture PPCA (Tipping & Bishop, 1999a) to integrate two monitoring statistics in order to address a fault diagnosis problem. However, all of the above extensions are based on the PPCA linear model lying in the Euclidean space around the mean. As such, these works cannot incorporate the nonlinear manifold underlying the data for dimensionality reduction and modeling, which is what we enable here.
In addition to the above, few studies have explored extending PPCA to include specific nonlinear manifolds. Probabilistic principal geodesic analysis (PPGA) (Zhang & Fletcher, 2013; Fletcher & Zhang, 2016) extends principal geodesic analysis (PGA) (Fletcher et al., 2003) into probabilistic framework. Mixture PPGA (Zhang et al., 2019) combines multiple PPGA models. Nodehi et al. (2020) develops the PPCA linear model within the Torus T space, as opposed the R² space, thereby extending torus PCA (Eltzner et al., 2018) to a probabilistic context. However, all these approaches require data lie precisely on top of specific manifold without any deviation from This assumption is not the case in many applications such as neuroscience, where neural activity data are distributed around manifolds with deviation and also exhibit noise. As such these prior methods have not found application to such datasets such as neural activity. Our method PGPCA is designed for such datasets and models observations that are probabilistically distributed around a given manifold that is first fitted from data. Unlike the above studies, Lawrence & Hyvärinen (2005) develops the Gaussian process latent variable model (GP-LVM), another nonlinear probabilistic model inspired by PPCA. The nonlinearity is encoded by kernel function between the latent states in GP-LVM. As these latent states are treated as parameters rather than random variables, GP-LVM is typically used for categorization tasks rather than for distribution modeling, which is goal. Given these disparate assumptions about the distribution of observations relative to the manifold and the properties/roles of the latent states, PGPCA addresses distinct application and thus serves complementary role compared with these prior studies.
3 METHODOLOGY
We first define the notations and the probabilistic model of PGPCA. Then we derive its log-likelibood and evidence lower bound (ELBO) for the EM algorithm. Finally, we summarize PGPCA EM by providing pseudo code (Algorithm 1).
3.1 PGPCA PROBABILISTIC MODEL.
We define the PGPCA model as
y: = ((2)) (((2))+ + K(Z)) X C x x
where all notations are listed in Table 1. Briefly, we have T observations yıtr E R", Each y: is composed of three parts. The first part is the !-dimensional manifold M = (0(x) Vz E 22 C R where 2 ~ p(z) 2) is the manifold state and random variable (R. in set Ωgg Essentially, 2, specifies the location on top of the manifold. Conditioned on 2,, the second part is a zero-mean normal distribution K (zz; Ca, where C is the loading matrix and K (z,) is the coordinate system for the
3

---

Published as a conference paper at ICLR 2025
Table I: PGPCA model notations
Notation
Description
Notation
Description
Yt R
observation at time [1, an iid random manifold state ~ p(z)
Φ(z+) RR
I-dim manifold R"
2, E Q2 C R
((,,) E ER
distribution coordinate at 2; loading matrix
mt E R
iid normal R.V. (00, I m
C RR
T1 E IR"
an iid normal R. ~ N(0, II.)
data distribution around the manifold, termed distribution coordinate, with orthonormal property (i.e., (z:)) (zz))((z,) (x,) = K (z,, (z:) I,, an identity matrix in Rnxx). Thus, C' follows the distribution coordinate K and determines the principal directions that cover most of the yır distribution. The third part, Te, with its isotropic variance o², captures any residual in Vı.r that is not already covered. We define the dimension of PGPCA model
the dimension of 2 equivalently the rank of the loading matrix (0 & m < n)). When m 0, C' is set to 0. Finally, our PGPCA model covers the PPCA model (Tipping & Bishop, 1999b) this case, the model (1) reduces to y Cæ. + T2 and the linear hyperplanes/subspaces are modeled by Ca., which is the same as in PPCA. Thus, PGPCA is generalization of PPCA and extends from the case where data is assumed to lie around the mean of data which can be considered as the central manifold PPCA the case where data can lie around nonlinear manifolds.
special case by setting ((zz) and K ((z)) I, In
3.2 PGPCA EM: E-STEP
We need to Iearn a PGPCA model (1) that describes the data the best. We formalize this learning problem as follows: given data y1:7, the manifold function & (that first fitted from data), and the distribution coordinate function K (either Euclidean geometric as we derive later in section 4), find the model parameters C, o², and p(z) in (1) by maximizing the data log-likelihood In p(B:.r). Since y: 's for different l's are iid (1 = 1 T), we can write the logllikel as
T
T
L In Inp(y)
Σ In p(u)lz)p(z)) dz
(2)
x=1 21
where p(g (z) is normal distribution from (1) such that
p(y.|=) = (d)22), (z)) (2))
o()))
= (())CKK() +
where' indicates the matrix transpose operation. To find the maximun-likeliood parameter estimates, we need to partial differentiate w.r.t. model parameters to maximize it; but this differentiation tricky because the integration in (2) is inside the In function. We address this challenge by deriving the ELBO _ of L following the standard EM procedure as
p(u(2)p((2) dz (z(2)
L Σ In 9x(2)
2 >Σ (12) In (p(),1=)((2)) (onnnI((pppp)p)l
In (x) dz E
(5)
i=1
where qx) is any probability distribution on Ωg From the standard EM procedure (Beal, 2003; McI.achlan & Krishnan, 2007), we know L = L if and only qı(2) p(2))r) for Vi E [1, T). Therefore, given the model parameters 02 and p(z). the E-step of PGPCA EM is derived
p p(z)
gi (2) p(zlur)
(6)
Ja_ y dz
p(y. 1- p()))
= E (ztar)
otherwise
Note that equation (7) follows after discretizing z() which is provided later in (9). This discretization necessary in fitting p(2) in the M-step and for numerical computations detailed next.
4

---

Published as a conference paper at ICLR 2025
Algorithm PGPCA EM
Input:
M1.r, model dimension landmark 21:M, manifold Φ(-), distribution coordinate K (-). Output: probability wI.M, parameters and o².
Initialize WI.M. C, and o² randomly.
repeat
[E-step)
Compute (e,;;) by (7) for Vi [1,7]& E [M-step)
Compute Wj by (11) for Vj E [1, M].
Compute (()) by (14) and then cig [ (T in descending order.
Compute a² by (16) and then C by (15).
until ELBO _ in (5) converges.
3.3 PGPCA EM: M-STEP TO FIND p(z)
Given (1: from the E-step, the M-step finds the optimal model parameters C' and o² in addition to p(z) maximize the ELBO L Only the first part of equation (5), q((z) In (p(y))))p(s x)))))) relates to these parameters, so we define
T
CM Σ ( (=)2 (p(uilz)(z)) dz Σ q: (x)(11(()))2)) q:(z) (mp()) dz (8)
Parameters C and o² only in the first in (8), which is defined Ltdl in (12), and the distribution p(z) is only in the second term in (8), which is defined as L2 in (10), respectively. But a challenge here is that we must first parameterize p(z) to learn To do so, we select M landmarks
{2 C 22 with nonnegative weights WI:M such that ΣØ = 1, and discretize p(z) as
p(z)
Σ Wy 8(z 2z))
(9)
where 8: RR - 10 is the Dirac delta function. This is how we discretize (z) E-step by substituting (9) into (6) to get (7). So the new M-step goal is: given (1)r(=), find parameters C, o² and WI.M to maximize L To find the optimized wI:M, we define LM the second term (8), as
C2 = -E (=))(n((z)d
ΣΣ (z)z Inwy
(10)
Using Lagrange multipliers (Bertsekas, 2014), the optimal wj to maximize L is found as
T for Vj [1.
W/ T
Σ ((2))
(11)
3.4 PGPCA EM: COMPUTING THE FIRST TERM L. TO DERIVE THE M-STEP FOR C AND o²
Next, we solve for parameters C and ²² that maximize LMM. the first term CM in (8). Due to the nonlinear manifold and the distribution coordinate K(z), finding the model parameters is more challenging than PPCA, which assumes a linear model. We first derive formula for LM (c.f. (13)), and then optimize it to find C' and o² in the next section. First, we expand L using (3) as follows:
T
CM =
Σ
(<(=) Inp(y)lz)
1=1
0
(x(z)
In 2 + In 12l + (v: o..' " (y) (.) dz (12)
1=1
5

---

Published as a conference paper at ICLR 2025
where the simplified notations ((z) Φ from (3) and (z) v from (4) for ease of exposition. The right-hamt side consists three parts that are added together. We compute these three parts of (12) one by one in appendix A. From there, we have
C1 =
T 2
x
n
2n + (n-m) In In In + trrlo (eC,,,
(13)
T M
T = T ==1==1
Ti,,, x (2₁) where T = - ( -.
(14)
Critically, compared to (12) where we started from, our derivations (appendix A) lead to all summa- tions and integrations being captured in [(g), which is interestingly independent of parameters C and a². This derivation makes partial differentiating L w.r.t. C and o² much easier and tractable, thus solving the major M-step challenge for learning the model parameters in the general case that includes nonlinear manifolds. Moreover, formula (13) is the same as PPCA log- log-likelihood (Tipping & Bishop, 1999a), except for the matrix (g)). This makes solving for C' and o² easy. We show this in detail in section 3.5.
3.5 PGPCA EM: M-STEP FOR C AND o²
Now we are ready to find the optimal C' and o² by maximizing LV/ in (13), Critically, our derivation showed that we can summarize all the nonlinear manifold and distribution coordinate information in one term (q)) within the L. As such, interestingly. (13) becomes a generalization of the PPCA log-likelihood in Tipping & Bishop (1999a) in that they have the same formula except that our I()) considers the manifold and the distribution coordinate on it, while PPCA's matrix S in Tipping & Bishop (1999a) does not. Therefore, we can solve for our optimal C' and o ousing the PPCA formula, and all the established guarantees in the PPCA theory also apply to this nonlinear manifold case. Here we write the optimal solution of C' and o² directly. The detailed derivation is in appendix B.
Define 7im the eigenvalues of eig(T (q))) in descending order. The optimal C is derived as
(4)) u; = 7, u
C UD where
Vi [1.
(15)
d: V
where D = diag(d1mm) and u, is the jth column of U and the 1th eigenvector of [(g). The optimal o2 is
12
(16)
71
m
Our pseudo code summarizes PGPCA EM in Algorithm The intuitive explanation behind our solution is that the loading matrix C captures the dominant directions in data y1.r distribution around the manifold Φ, and x² captures the residual directions with their average variance as a noise term. Since all steps in Algorithm are analytical, PGPCA EM is efficient in terms of training time (appendix C), similar to classical EM for linear state-space models (Roweis & Ghahramani, 1999).
4 EXPERIMENTS
We show that our PGPCA model plus its EM algorithm can solve four problems: (1) Given data Fı.r, PGPCA EM can learn an m-dimenssiona probabilistic model that includes a given underlying nonlinear manifold & and a distribution coordinate K K(z), which can be either Euclidean or computed according to our geometric distribution coordinate. (2) It allows us to perform hypothesis testing to select the Euclidean geometric distribution coordinate by fitting alternative PGPCA models with two different K (z)'s, 's, and selecting the one with the higher data likeliiol L in (2). (3) We can perform dimensionality reduction by fitting
PGPCA model with any dimension m E [0, n) that as low as the user desires. (4) PGPCA EM can not only learn the data distribution around the manifold but also the distribution the manifold; indeed, show that the weights of

---

Published as a conference paper at ICLR 2025
(A) K(z) = EuCOV
K(z) = GeCOV
(B)
K(t) = EuCOV
K(z) GeCOV
1st distribution coordinate
1st distribution coordinate
2nd distribution coordinate
2nd distribution coordinate
3rd distribution coordinate
Y2
V1
Figure 1: Distribution coordinate (22) be Euclidean (EuCOV) geometric (GeCOV) loop and a torus. (A) When (2) = EuCOV on a loop C R22 always aligns with the embedding coordinate R matter where is on the loop. In contrast, if K (=) GeCOV, the distribution coordinate follows the tangent vector and the normal vector. (B) Again, K (z) can always align to the axes of R (BuCOV) or composed of two langent vectors plus another vector perpendicular the torus surface (GeCOV). The top-left inset figure shows the PPCA case whose manifold is the mean of data (the black dot) with its only distribution coordinate system, which equal EuCOV.
Table PGPCA (GeCOV/EuCOV), PPCA, and FA log- log-likeli full-rank models (m = n)
loop in IR²
loop in 20
torus in R
data analysis
True
GeCOV
EuCOV
GeCOV
EuCOV
GeCOv
EuCOV
Mouse12
Mouse28
GeCOV
-2.931
-2.725
-27.211
-28.484
-5.626
-5.560
-31.758
-24.752
EuCOV
-2.939
-2.698
-27.993
-28.356
-5.631
-5.523
-31.908
-25.089
PPCA
-3.048
-2.991
-31.945
-31.677
-5.862
-5.907
-34.622
-29.316
FA
-3.048
-2.991
-31.945
-31.677
-5.862
-5.907
-34.615
-29.310
manifold latent state distribution WI.M (from discretizing p(z) in (9)) can be jointly learned with parameters C' and o² in (1) and result in a similar log-likelihood as the true model.
We show that PGPCA can solve the above four problems using neural data analyses and extensive simulations covering various nonlinear manifolds, distribution coordinates, and manifold latent state distributions p()). The nonlinear manifolds include loop (in R2 or R 10) and torus. The distribution coordinate (2) can be Euclidean (EuCOV) or geometrie (GeCOV) (see Figure and appendix C). For the torus, its p(z) has two options: a uniform distribution
the angular space (0,2m] x [0, 2T]] (uniAng), or a uniform distribution on the torus surface (uniTorus) (Figure in the appendix). The real dataset includes neural spike firing rates recorded from anterodorsal thalamic nucleus (ADn) of mice, part of the thalamo-cortical head-direction (HD) circuit, while animals were exploring an open environment (Peyrache al., 2015; Peyrache & Buzsáki, 2015). The firing rates are projected to R 10 following the same preprocessing as that in prior work (Chaudhuri et al., 2019). Details of neural data analyses and simulations are in appendix C. summary of the log- likelihooos for each model, based on the neural data analyses and simulations, is presented in Table 2. In this table, all models are set full rank (m n) to maximize their log-likelihoods, which makes PPCA mathematically equivalent FA. However, slight differences in the log-like likeliooss between PPCA and FA are observed since their models are learned numerically.
4.1 PGPCA FINDS THE CORRECT DISTRIBUTION COORDINATES ON A 1D LOOP
We first show that PGPCA EM can learn a nonlinear probabilistic model from data and distinguish different distribution coordinates in hypothesis testing. To show that our method succeeds in incorpo- rating the manifold, we compare with the widely used PPCA, which linear. Figure 2A shows the 2D probability distribution from the true models and from the learned models by PGPCA/PPCA. First, we see that the learned PGPCA model's 's distributions are closer to the true model's 's distribution compared to PPCA's distribution, matter what the distribution coordinates (GeCOV/EuCOV) in the true or learned PGPCA models are. This demonstrates the importance of modeling data proba- bilistically with underlying nonlinear manifold enabled by PGPCA. Moreover, the true model's

---

Published as a conference paper at ICLR 2025
Manifiid
loop
[A]
Trus modal, GecOy
(B) PGPPA over dimensons
[C]
PGPCA convergence Iove motel 6e00V)
(01 PSPCA over cimensions
modell
PGPCC (aeCDv)
PSPCA jludoy
PRGA
PORCA POPCA
HOPEA HUCIN
2042
28.
NCNCA
POTCA (Caco1)
POTCA (0a000)
PONCA dimaniion
True mode
PSPCC (iecCv)
POPCCA StaceNt
PPCA
model EuCOV)
PORCA
POPCA
%
10 20
35
40
03 01434 4 5 6 1
NOYGA
Asesston DO
smowien
Figure 2: PGPCA can recover the true model distribution, distinguish different distribution coordinates (2), and perform dimensionality reduction simultaneously. Across all panels (A-D), the true model's K (=) are GeCOV and EuCOV in the top and the bottom row, respectively. (A) The PGPCA (GeCOV) model learned by EM recovers the true model distribution while PPCA does not. Further, PGPCA can do so only with the correct K (z), showing that PGPCA can distinguish the correct coordinate. (B) PGPCA with the correct K (2) always has higher trial-average log-likelihood (paired t-test: top and bottom < 1.7 x 10-12). (C) Learned PGPCA models with different dimensions
m E [0, 10] (color bar) and with different (2) (EuCOV//eCOV) converge within 40 EM iterations.
(D) The same conclusion in (B) also holds here when the loop C R (paired t-test: top and bottom
< <.11 x 10-4).
distribution is only recovered by the learned PGPCA model when their distribution coordinates match. Figure 2B and Table 2 further confirm that the learned PGPCA model with the correct distribution coordinate K (z) has higher log- ligeeliodhood than the learned PGPCA model with the incorrect one. As such, fitting the two alternative PGPCA models and comparing their log- ligllikll can successfully distinguish the true distribution coordinate underlying the data. This shows PGPCA's ability to solve problems (1) and (2) listed at the beginning of section 4.
4.2 PGPCA CAN PERFORM DIMENSIONALITY REDUCTION.
For dimensionality reduction, we simulate a true model built on a loop embedded in IR_ so we have enough dimensions for this application. Figure 2C shows that regardless of whether the true model is GeCOV or EuCOV, PGPCA EM can converge under any PGPCA model dimension m E [0, 10]. Thus, this EM method can robustly learn PGPCA model with any dimension. In Figure 2D, the learned PGPCA model with the correct K (x) always has a higher log-likelihood compared to the alternative, even when the PGPCA dimension is selected to be low. Therefore, PGPCA distinguish the correct distribution coordinate K(z) even when its dimension is chosen low. This result shows that PGPCA can simultaneously perform both dimensionality reduction and distribution coordinate selection, solving problem (3) stated at the beginning of section 4.
still
4.3 PGPCA CAN RECOVER THE TRUE MODEL'S 'S DISTRIBUTION EVEN WHILE LEARNING p(z).
Figure 3A shows that PGPCA EM can recover the true model's distribution when its K K(2) matches the true one, regardless of whether p(z) is given or learned. We also find that PGPCA can again distinguish the correct coordinate system K (2) even when p(z) is being jointly learned. This shows that PGPCA EM can learn not only the distribution around the manifold, but also the distribution p(z) on the manifold. This solves problem (4) stated at the beginning of section 4.
Furthermore, Figure 3B shows that the learned PGPCA model with its K (z) matched to the true one always has a higher log-like likelihood than the learned PGPCA model with the unmatched (xz), whether p(z) is learned or not. Thus, the hypothesis testing ability of PGPCA EM in distinguishing different distribution coordinates also holds even when simultancously learning p(z). The average performance across uniA mmiAng/uniToon and given/learned p(=) is provided in Table 2.
8

---

Published as a conference paper at ICLR 2025
learned
Twirreee
FSPCA
PGPCA
FONCA p0l1
XPCA
K10"
541
5.07
112
prodetiiyy
EuCO BUCTI
Figure 3: PGPCA EM can recover the true model's distribution even when simultaneously learning the manifold state probability p(z), while PPCA does not. (A) The first row shows that when the true model's p(z) = uniAng and (()) = GeCOV, the learned PGPCA model's distribution is similar to the true one, regardless of whether p(z) is given (column 2) learned (column 3). Also, this is the case only if PGPCA's ((z) is GeCOV (true coordinate). showing its ability to identify the coordinate. Rows 2-4 show the same conclusion for alternative true models having different p(z) and (=)). (B) The trial-average log-likelihood of the four learned PGPCA models (columns 2-5 in (A)). Again, the learned PGPCA model whose K (z) matches the one always has higher log-likelihood than the unmatched PGPCA model, regardless of whether p( (() is given or learned, showing hypothesis testing capability. For all rows with given or learned p(2), paired
t-test < 2.4 x 10
Mouse
Manifiitt toop
renal readivo
teyelvon
Hag Ecathroot
Bickoood
-30
-20
kam
am
025
25.7
50 10
Figure 4: PGPCA (GeCOV) better captures the distribution of neural firing rates in mice head direction circuit. (A) The fitted loop manifold in IR-0 and the neural data distributed around it. (B) PGPCA (GeCOV) model consistently has higher lig-likelioood than PGPCA (EuCOV) model across all dimensions. (C) and (D) are the same as (A) and (B) for second mouse, with conclusions being the same.
4.4 PGPCA CAN DISTINGUISH DISTRIBUTION COORDINATES ON REAL DATA
We applied PGPCA to neural firing rates recorded from the thalamo-cortical head direction circuit of six mice, and we select two mice examples here. First, we found that the main manifold structure was a loop, consistent with prior work, and so fitted this loop using a cubic spline with 10 knots selected by K-means (appendix C). Figures 4A and 4C display the projected neural firing rates along with the fitted manifolds, which are the ID loops embedded in RR 10 The neural firing rates distributed not precisely on, but around, the manifold, indicating that the main manifold alone

---

Published as a conference paper at ICLR 2025
is insufficient for completely modeling noisy data. This observation underscores the necessity of PGPCA, which captures the deviation outside of the manifold through distribution coordinates and noise. So we constructed the distribution coordinate and ran PGPCA EM, and compared with PPCA and FA.
Table 2 and Figure 6 in the appendix demonstrate that PGPCA significantly outperforms PPCA and FA. Further and interestingly. Figures 4B and 4D show that PGPCA GeCOV more accurately captures the firing rates than PGPCA EuCOV across both mice. This latter result suggests that the noise not accounted for by the main loop manifold also originates from the same geometric structure rather than being in the Euclidean space. This inference can only be made using a model with a coordinate system around the main manifold, which is major capability provided by PGPCA. This again shows that PGPCA can also perform hypothesis testing about the coordinate system in which data is distributed. These conclusions again held on the other mice (Table 3 in the appendix).
CONCLUSION
We developed PGPCA, a method that generalizes the widely-used PPCA for analyses of data that are distributed around given nonlinear manifold that is fitted from data. Unlike PPCA, which assumes that data lies around the mean in Euclidean space, PGPCA incorporates the nonlinear manifold as well as distribution coordinates attached to this manifold to capture deviations from it and noise. Also, in addition to being able to use the Euclidean coordinate around the manifold, PGPCA can also compute a geometric coordinate system around the manifold, which we derived here. Finally, PGPCA can perform hypothesis testing to pick between the Euclidean and the geometric distribution coordinates based on which can better describe the data distribution. In this paper, we focused on the Euclidean (EuCOV) and our geometric (GeCOV) K (z) because they naturally arise from the linear embedding space R" and the underlying nonlinear manifold, respectively. If prior knowledge about the data suggests the hypothesis that another new form of (x) is needed (assuming it can be derived), PGPCA can serve as a tool to validate or reject this hypothesis by comparing its tog-likelihood with that of other ((2) options, such EuCOV and GeCOV. Our PGPCA can accommodate new K(2) because in deriving the PGPCA EM algorithm, we did not impose any specific assumptions on K K(z) beyond the basic orthonormal property.
We demonstrated the success of PGPCA and its efficient analytical EM learning algorithm on real neural firing rate data and over three types of simulated manifolds with different manifold state distribution p(z) (uniAng/uniTorus) and different distribution coordinates K K(*) (EuCOV/GeCOV), Our results show that PGPCA can correctly i) fit the nonlinear probabilistic model, ii) distinguish between Fuclidean and geometric distribution coordinates, iii) perform dimensionality reduction, and iv) learn the manifold state distribution on top and around the manifold. Further, in both simulations and real neural data, PGPCA outperformed PPCA by capturing the manifold. One major application of PGPCA is for modeling of neural data time-series in the fields of neuroscience and neurotechnology given the evidence that neural data distribute around nonlinear manifolds (section 1). However, PGPCA is not limited to neural data time-series and can in principle be applied to any time-series dataset with data-fitted underlying manifold. A limitation of PGPCA, similarly to PPCA and PCA, is that it is a static dimensionality reduction method and thus does not explicitly model the auto-correlations in data. Further, similar to these methods, PGPCA assumes that the data distribution is stable over time. Further work can extend PGPCA to enable manifold -bsedd dynamical analyses Abbaspourazad al., 2024; Sani et al., 2024) adaptive modeling track non-stationarity in the data distribution (Ahmadipour et al., 2021; Yang et al., 2021; Degenhart et al., 2020). Finally, PGPCA allows for incorporation of manifold knowledge, which should first be obtained based on data using existing manifold identification and fitting methods (e.g.. TDA and splines). As we show in our real neural data analyses, incorporating this knowledge can substantially improve dimensionality reduction and distribution modeling.
ACKNOWLEDGMENTT
This work was partly supported by the National Institutes of Health (NIH) grants DP2-MH126378 and RFIDA056402, and by the Army Research Office (ARO) under contract W911NF- 16- 1-0368 as part of the collaboration between the US DOD, the UK MOD and the UK Engineering and Physical Research Council (EPSRC) under the Multidisciplinary University Research Initiative (MURI).
10

---

Published as a conference paper at ICLR 2025
REFERENCES
Hamidreza Abbaspourazad, Eray Erturk, Bijan Pesaran, and Maryam M Shanechi. Dynamical flexible inference of nonlinear latent factors and structures neural population activity. Nature Biomedical Engineering. 8(1):85-108, 2024.
Parima Ahmadipour, Yuxiao Yang, Edward F Chang, and Maryam M Shanechi. Adaptive tracking of human ecog network dynamics. Journal of Neural Engineering, 18(1):01601,, 2021.
Long-Hoon Ahn and Jong-Hoon Oh. A constrained algorithm for principal component analysis. Neural Computation, 15(1):57-65, 2003.
David Bartholomew, Martin Knott, and Irini Moustaki. Latent variable models and factor analysis: A unified approach, volume 904. John Wiley & Sons, 2011.
Matthew James Beal. Variational algorithms for approximate Bayesian inference. University of London, University College London (United Kingdom), 2003.
Dimitri P Bertsekas. Constrained optimization and Lagrange multiplier methods. Academic press, 2014.
Borislav D Bojanov, H Hakopian, and B Sahakian. Spline functions and multivariate interpolations, volume 248. Springer Science & Business Media, 2013.
Rishidev Chaudhuri, Berk Gerçek, Biraj Pandey, Adrien Peyrache, and lla Fiete. The intrinsic attractor manifold and population dynamics of canonical cognitive circuit across waking and sleep. Nature neuroscience, 22(9):1512-1520 2019.
Joost CF De Winter and Dimitra Dodou. Factor recovery by principal axis factoring and maximum likelihood factor analysis as a function of factor pattern and sample size. Journal of applied statistics, 39(4):695-710, 2012.
Alan D Degenhart, William E Bishop, Emily R Oby, Elizabeth C Tlyler-Kabare Steven M Chase, Aaron P Batista, and Byron M Yu. Stabilization of a brain-computer interface via the alignment of low-dimensional spaces of neural activity. Nature biomedical engineering, 4(7):672-665, 2020.
Wei Deng and Radu V Craiu. Exploring dimension learning via a penalized probabilistic principal component analysis. Journal of Statistical Computation and Simulation, 9211iin66292)2226-99999999 2023.
Manfredo P Do Carmo. Differential geomelry of curves and surfaces: revised and updated second edition. Courier Dover Publications, 2016.
Benjamin Eltzner, Stephan Huckemann, and Kanti V. Mardia. Torus principal component analysis with applications to ma structure. The Annals of Applied Statistics, 12(2):1332 - 1359, 2018. doi: 10.1214/1-A0AASI11 URL https: doi org/ 10. 1214 AAAAAA1111
Charles Fefferman, Sergei Ivanov, Yaroslav Kurylev, Matti Lassas, and Hariharan Narayanan. Fitting putative manifold to noisy data. In Conference On Learning Theory, pp. 688-720. PMLR, 2018.
Charles Fefferman, Sergei Ivanov, Matti Lassas, and Hariharan Narayanan. Fitting manifold of large reach to noisy data. Journal of Topology and Analysis. pp. 1-82, 2023.
Thomas Fletcher and Miaomiao Zhang. Probabilistic geodesic models for regression and dimen- sionality reduction on riemannian manifolds. In Riemannian Computing in Computer Vision, pp. 101-121. Springer, 2016.
Thomas Fletcher, Conglin Lu, and Sarang Joshi. Statistics of shape via principal geodesic analysis lie groups In 2003 IEEE Computer Society Conference Computer Vision and Patlern Recognition, 2003. Proceedings, volume 1, pp. 1-1. IEEE, 2003.
Richard Gardner, Erik Hermansen, Marius Pachitariu, Yoram Burak, Nils A Baas, Benjamin A Dunn, May-Brit Moser, and Edvard I Moser. Toroidal topology of population activity grid cells. Nature, 123 2022.
11

---

Published as a conference paper at ICLR 2025
Michael Greenacre, Patrick JF Groenen, Trevor Hastie, Alfonso lodice d' d'Enza, Angelos Markos, and Elena Tuzhilina. Principal component analysis. Nature Reviews Methods Primers, 2(1):100, 2022.
Michael Hahsler and Kurt Hornik. Isp-infrastructure for the traveling salesperson problem. Journal of Statistical Software, 23(2):1-21, 2007.
David A Harville. Matrix from statistician's perspective. Taylor & Francis, 1998.
Trevor Hastie, Robert Tibshirani, Jerome H Friedman, and Jerome Friedman. The elements of statistical learning: data mining, inference, and prediction, volume 2. Springer, 2009.
Xuming He and Peide Shi. Bivariate tensor-product b-splines in partly linear model. Journal of Multivariatte Analysis, 18(2):62--111, 1996.
Johan Ludwig William Valdemar Jensen. Sur les fonctions convexes et les inégalités entre les valeurs moyennes, Acta mathematica, 30(1):175-193, 1906.
Kristopher Jensen, Ta-Chu Kao, Marco Tripodi, and Guillaume Hennequin. Manifold gplvms for discovering non-euclidean latent structure in neural data. Advances in Neural Information Processing Systems, 33:22580-22 225922, 2020.
Rajiv Khanna, Joydeep Ghosh, Russell Poldrack, and Oluwasanmi Koyejo. Sparse submodular probabiliistic pca. In Artificial Intelligence and Statistics, pp. 453-461. PMLR, 2015.
Neil Lawrence and Aapo Hyvärinen. Probabilistic non-linear principal component analysis with gaussian process latent variable models. Journal of machine learning research, 6(1 1), 2005.
Steven J Leon, Åke Björck, and Walter Gander. Gram-schmidt orthogonalization: 100 years and more. Numerical Linear Algebra with Applications, 20(3):492-532 2013.
Geoffrey J McLachlan and Thriyambakam Krishnan. The EM algorithm and extensions. John Wiley & Sons, 2007.
Tristan Needham. visual explanation of jensen's inequality. The American mathematical monthly, 100(8):768-7777 771, 1993.
Anahita Nodehi, Mousa Golalizadeh, Mehdi Maadooliat, and Claudio Agostinelli. Torus probabilistic principal component analysis. arXiv preprint arXiv:2008. 10725, 2020.
Chongsun Park, Morgan C Wang, and Eun Bi Mo. Probabilistic penalized principal component 2017.
analysis. Communications for Statistical Applications and Methods,
Beresford N Parlett. The qr algorithm. Computing science & engineering, 2(1):38-42, 2000.
Kaare Brandt Petersen, Michael Syskind Pedersen, et al. The matrix cookbook. Technical University of Denmark, 7(15):510, 2008
Adrien Peyrache and György Buzsáki. Extracellular recordings from multi-site silicon probes in the anterior thalamus and subicular formation of freely moving mice. CRCNS. org, 2015.
Adrien Peyrache, Marie M Lacroix, Peter Petersen, and György Buzsáki. Internally organized mechanisms of the head direction sense. Nature neuroscience, 18(4):569-57) 575, 2015.
Sam Roweis and Zoubin Ghahramani. A unifying review of linear gaussian models. Neural computation, 11()):05-3455 1999.
Walter Rudin et al. Principles of mathematical analysis, volume 3. McGraw-hill New York, 1964.
Omid G Sani, Bijan Pesaran, and Maryam M Shanechi. Dissociative and prioritized modeling of behaviorally relevant neural dynamics using recurrent neural networks. Nature neuroscience, 27 (10):2033-204) 2024.
Gurjeet Singh, Facundo Memoli, Tigran Ishkhanov, Guillermo Sapiro, Gunnar Carlsson, and Dario L Ringach. Topological analysis of population activity in visual cortex. Journal of vision, 8(8): 11-11, 2008
12

---

Published as a conference paper at ICLR 2025
Joshua B Tenenbaum, Vin de Silva, and John C Langford. A global geometric framework for nonlinear dimensionality reduction, science, 290(5500):2319-2333 2323, 2000.
Michael E Tipping and Christopher M Bishop. Mixtures of probabilistic principal component analyzers. Neural compulation, 11(2):43-482, 1999a.
Michael E Tipping and Christopher Bishop. Probabilistic principal component analysis. Journal of the Royal Statistical Society Series B: Statistical Methodology. 61(3):611-62) 1999b.
Yuxiao Yang, Parima Ahmadipour, and Maryam M Shanechi. Adaptive latent state modeling of brain network dynamics with real-time learning rate optimization. Journal of Neural Engineering, 18(3): 036013, 2021.
Junsong Yin, Dewen Hu, and Zongtan Zhou. Noisy manifold learning using neighborhood smoothing embedding. Pattern Recognition Letters, 29(1 1011131:11116,0 2008.
Shipeng Yu, Kai Yu, Volker Tresp, Hans-Peter Kriegel, and Mingrui Wu. Supervised probabilistic principal component analysis. In Proceedings of the 12th ACM SIGKDD international conference on Knowledge discovery and data mining, pp. 464-473, 2006.
Jingxin Zhang, Hao Chen, Songhang Chen, and Xia Hong. An improved mixture of probabilistic pea for nonlinear data-driven process monitoring. IEEE transactions on cybernetics, 49(1):198-21) 2017.
Miaomiao Zhang and Tom Fletcher. Probabilistic principal geodesic analysis. Advances in neural information processing systems, 26, 2013.
Youshan Zhang, Jiarui Xing, and Miaomiao Zhang, Mixture probabilistic principal geodesic analysis. In Multimodal Brain Image Analysis and Mathematical Foundations of Computational Anatomy: 4th International Workshop, MBIA 2019, and 7th International Workshop, MFCA 2019, Held in Conjunction with MICCAI 2019, Shenzhen, China, October 17, 2019, Proceedings 4, pp. 196-208. Springer, 2019.
Wenni Zheng, Pengbo Bo, Yang Liu, and Wenping Wang. Fast b-spline curve fitting by I-bfgs. Computer Aided Geometric Design, 29(7)4448-462, 2012.
13

---

Published as a conference paper at ICLR 2025
A DERIVE THE CONCISE FORM OF C1
In this appendix, we compute three parts of (12) one by one below to transform (12) into (13). The first part L Since (z(2) is probability distribution Ω2 and In 2x is constant, we have
SM
Σ R2 (((z) x 2 dz T x n 2
(17)
The second part 2 Recall that for the coordinate system, K,K, = I,, Following (4) and the matrix determinant lemma (Harville, 1998), have
V KK,CCK I'
= =|02In|
Im C'K. x (o² x K.C
(o?)" x 1 Im + CCC)
= (o2)m I C'Cl
(18)
Interestingly, this derivation shows that | a constant independent of z. Therefore, the second part simplifies to
- T' [ m)
+ +In 1 c'c ln (m-m)
(19)
The third part C/ The key idea here is to rewrite the vector norm weighted by - 1 i...
|/ &=1l=:1 (y: " (y:
Φz),
using the identity ir( (r(BB) (A) (when AB and A are well-defined) (Petersen et al., 2008),
which gives the following
y: &=11=.' == (y: Φ2)
(yi x.) ur T
y (5) Φ2)
(20)
define as IT.,
Remember K.K- In. so inverse of = in (4) is
W V = (K,CC' + o² K, K. x (0'," K
(21)
Now we transform (20) into
tr
= tr
K (o CC'
= =tr((0²,, + C K
(22)
define T..:
From (20) and (22), since the trace, summation, and integral operators are linear and can be swapped, the third part L be written
Σ 1=1
ur
v 11 dz = T x CCC--
+ x r(w]
(23)
where,
[(g)
Σ
P g:(2) dz
(24)
S2
Finally, since CM - x (cM, + LMM + CM.) LM (12) is equal (13) by combining the derived forms above for the three parts in (17), (19), and (23).
The last thing to notice is that (2) is discretized in (7) for numerical computations. So in practice, can compute [g)) defined in (24) numerically (14). This completes the derivation.
14

---

Published as a conference paper at ICLR 2025
B DERIVE OPTIMAL C AND o² IN PGPCA MODEL
In this appendix, we derive the optimal and x² by maximizing L in (13). Our (13) and the PPCA ongllikelhhood in Tipping & Bishop (1999a) have the same formula except that our (g)) considers the manifoid and the distribution coordinate on it, while PPCA's matrix S in Tipping & Bishop (1999a) does not. Therefore, we solve for our optimal C' and o² using the PPCA formula. We rewrite the derivation for completeness and notation consistency below. For more details, please refer to Tipping Bishop (1999a).
We first optimize C. Using matrix calculus operations, we have
OC
T x 2 2A - =0
(25)
ac
So the optimal satisfies the following condition
(26)
We define the SVD of C = UDV' where U E Raxm D E mR is diagonal, and V E m=
By Woodbury matrix identity (Petersen et al., 2008), we have
A A = (o'1, CCC
=Cx I CCCC--
= =UDV' x (V + D"
UDD (0 Im + DD)-LV'
(27)
Substituting (27) in (26) and multiplying Vio D2 D on both sides, we have
F(g)U C V (02 m + DD
=UD (o I D-)
= =U (o D
(28)
Note that D and 0² Im + D2 can be swapped because they are diagonal. From (28) and C UDV' we conclude that
can be any orthonormal matrix in mmmxm For convenience, we set = Im-
2. Columns of U are eigenvectors of (g). Define eigenvalues eig (T)())))) {imm without order (e.g., axcending/descending Then U = [ui]... [mm]] (u; is the 1 column of U) and diag(d (dim) such that u, is the eigenvector w.r.t. eigenvalue % σ² + d? from (28)
In summary, given o², the optimal C is
((g) 7/s d = V -2
C = UD where
vi [1.m]
(29)
where ui 7h, and d, are defined above. The next step is optimizing o². To do so, substitute C from (29) into (13), and then rewrite LM as
CM = 2 x n + (n - mn +
In 7s + σ²
x
Σ
% +
m
(30)
i=m+1
ac 0, the optimal o² is Setting 3.
n
02
%
(31)
n
m
$=m+1
The remaining challenge now is that the optimal and o² in (29) and (31) do not complete the answer yet because also have to select 71:m. from (i)( The power of our derivation for LM
15

---

Published as a conference paper at ICLR 2025
formula in (13) is that because the manifold and the distribution coordinate summarized in the (4)) term, we can use the results in Tipping & Bishop (1999a) directly to select Yimm- Briefly, to do this "/.m selection, substitute (31) into (30) to rewrite L again as
T CM = - x { Inln2 2
Σ la
Ia% + n + (n m) x In
m
i=+++
1=m+1
(32)
Note that == In; is a constant because {mmm) = eig( Therefore, maximizing C in (32) equivalent to minimizing
In
l1%
(33)
(=++1
(=++1
which a Jensen's ' inequality (Needham, 1993; Jensen, 1906). It's proved in Tipping & Bishop (1999a) that the optimal 7mm+I:on I:n for (33) must be consecutive series in (i)(T More precisely, defining Tınn as the descending series of eig(T we have that 7m is a consecutive series in Tim-
Finally, from (29) and (31), we see that
vj E m1. 7j 0²
%
(34)
Therefore, Tt.m cannot include 7n. the smallest eigenvalue of ( 7n 7m+tmn Combined with the consecutive condition on 7m+ lin we can conclude that 7/ 7/ for Vi = n]. Then (29)
and (31) become (15) and (16), respectively. This completes the M-step of PGPCA EM.
C SETTING OF ALL SIMULATION CASES AND DATA ANALYSIS.
We describe the details of simulations and data analysis below. For simulations, we simulate 3 kinds of manifolds: a ID loop in IR2, a ID loop in 10 and a 2D torus in R3, First, we rewrite the PGPCA model (I) as
y: = ((zz) + K(z)) x (Cx, + +r)
(35)
because K(z;) x I²I, , x K Kzz)) = I²I,, so the covariance of K K(z,) and are the same. Note that
Cov(Ca: + CC +0-L =
(36)
Therefore, every simulation case is specified by the manifold function & with the manifold latent state distribution p(=) top of it, distribution coordinate K(z). and basic covariance A. We describe all simulation cases for the above three manifolds below.
A ID loop embedded in IR?. We define [0, 2n] with p(z) (0, 2T) where U(a b) is continuous uniform distribution within [a, bj. The nonlinear manifold is an ellipse with function
((z) (cos(2), 2), sin)))). The basic covariance is taken as A diag (I0.1, The distribution coordinate K (=) for the simulated data can be Euclidean or geometric. Euclidean means K K(z) 12. and the distribution coordinate system follows the embedded Euclidean coordinate at all points on the manifold. We refer to this scenario as the Euclidean covariance (EuCOV) since A follows the Euclidean coordinate (Figure IA, left). On the contrary, the geometric case refers when K()) is composed of the tangent and normal vectors at each 2 along the manifold, which is an ellipse here. We refer this alternative scenario as the geometric covariance (GeCOV, Figure 1A, right). For convenience, we also say that K (2) is EuCOV or GeCOV when it's the Euclidean geometric coordinate, respectively. Because there are two options for K(z) corresponding to the EuCOV and GeCOV scenarios respectively, we will fit two types of PGPCA models with PGPCA using either Euclidean or a geometric K(2), which we term PGPCA EuCOV and PGPCA GeCOV, respectively; this thus leads to simulation cases for this ellipse in R For both cases, generate 5000 training samples from the true model and learn every PGPCA model (EuCOv/Gecov) with 500 landmarks 21:500 in (9) using 20 EM iterations. The PGPCA model dimension can be E [0, 2].
16

---

Published as a conference paper at ICLR 2025
A ID loop embedded in R To simulate this loop in higher dimensional space, define the manifold points 2 E [o. L] with p(z) = U(0, L) where L is the length of the loop. We form the manifold Φ(2)
cubic spline with 6 knots, and with length L. The basic covariance is A = diag( (20, 2, 18, 4,.... 12, 10)). The PGPCA model can be either FuCOV or GeCOV. In this case, the (2) in GeCOV model is computed using the Gram-Schmidt process (Leon et al., 2013) with the tangent vector as the first vector, and the Euclidean axes e1:10 as the other independent vectors be orthogonalized by the Gram-Schmidt process one by one in sequence. Again, there are two simulation cases w.r.t. this spline in I 1 corresponding to EuCOV or GeCOV being the true distribution coordinate, respectively. For both cases, we generate 5000 samples from the true model, and learn every PGPCA model (BaCOV/GeCOv) with 500 landmarks 21:500 using 40 EM iterations. The dimension of PGPCA model can be m E [0, 10].
A 2D torus embedded in R. Defining [0, 2m) x [0. 2T), the torus manifold is given by Do Carmo (2016)
(x)) [33 + cos 22) cos 21, (3 + cos 22) sin 21, sin 2g]
Here the basic covariance A- = diag( (0.1, 1, 0.3, 0.5)) and K (x)) can be either EuCOV or GeCOV (Figure IB). GeCOV K(2) is composed of two tangent vectors (Js)
@ and dsy and their cross product.
We also give p(z) two options: a uniform distribution on the angular space (0, 2n] [0, 2x] (uniAng), or a uniform distribution on the torus surface (uniTorus). Because we have two options for ((z) and p(=), there are 2 x 2 4 simulation cases w.r.t this torus in R Figure 5 shows the true model distributions under the 4 cases. For all four cases, we generate 50000 samples from the true model, and every PGPCA model (uniAng/uni/Torus x FuCOV/GeCOV) with 1000 landmarks 21:1000 is learned using 40 EM iterations. We increase the number of training samples because we need to fit 47111000 = 100000000) in all four cases. The dimension of PGPCA model can be m E [0.
Performance measures in simulations. After learning the PGPCA model from the training samples using the PGPCA EM algorithm with one of three manifolds above, for each simulation case, we generate 20 test trials from the true model. Each trial includes 2000 samples. For each of the 20 trials, we measure the performance of learned PGPCA model with the average liceliicl defined as f/T with 7 2000 being the trial length. In the figures, all logllikels for the learned PGPCA models are the average of these 20 trial-average log-likelihoods, and comparisons are done with paired t-tests between the raveraee log-likelihood groups from two different learned PGPCA models (i.e., 20 trials in the paired t-test comparisons).
For data analysis, we utilized the neural firing rates recorded from mice's brains. This dataset is publicly available (Peyrache & Buzsáki, 2015), and further details can be found Peyrache et al. (2015). The preprocessing steps prior to applying PGPCA are primarily based Chaudhuri et al. (2019). These steps are summarized below for completeness.
Data. For all 6 mice, spikes were recorded from intracortical shanks implanted in the anterodorsal thalamic nucleus (ADn), a part of the thalamo-cortical head-direction (HD) circuit, while the mice were exploring an open environment. There are 8 shanks with 50 cells for Mouse 12 and 4 shanks with 22 cells for Mouse 28. The sampling rate is 20 kHz. The numbers of shanks and cells of other mice are listed in Table 3. All mouse data have the same preprocessing and PGPCA training and testing procedures.
Preprocessing. We followed the preprocessing steps in prior work. For both mice, we first computed the firing rates by smoothing the spike time-series with a Gaussian kernel with a standard deviation of 100 ms. The firing rates were then down-sampled to 15000 samples with 100 ms step size (equivalent to data-duration of 25 minutes in total). As preprocessing following prior work, we first applied a square root on the firing rates to stabilize the variance (Chaudhuri et al., 2019), and then projected the data using Isomap (Tenenbaum et al., 2000) from 10 (Mouse 12) and IR22 (Mouse 28) to 10 This 10D space the space in which PGPCA operates, as shown Figure 4.
PGPCA training and testing. We split the 15000 samples equally into 5 trials for 5-fold cross- validation. In each fold, we concatenated trials to form a training set. Similar to what has been observed previously (Chaudhuri et al., 2019), we found that neural data was distributed around loop manifold, but had both noise and deviations from it. We thus fitted a 1D loop in 10 using K-means (Hastic et al., 2009) with 10 clusters. The means of these clusters served as the knots of a closed cubic spline. We determined the order for connecting these knots by solving the traveling salesman problem (Hahsler & Hornik, 2007). This resulted in manifold model di constructed by
17

---

Published as a conference paper at ICLR 2025
Table 3: PGPCA (GeCOV/EuCOV) and PPCA log-likelihood of full-rank models (m n)
data analysis
True
Mouse12
Mouse17
Mouse20
Mouse24
Mouse25
Mouse28
shanks cells
8
4
50
29
10
GeCOV
-31.758
-30.407
-19.751
-18.687
-20.298
-24.752
EuCOV
-31,908
-30.595
-19.768
-18,702
-20.356
-25.089
PPCA
-34.622
-32.668
-21.560
-21.134
-24.977
-29.316
Table 4: PGPCA computational complexity for every iteration
computational step
time complexity
E-step
Compute (=;) for Vi E [1. and Vj E M).
O(TM²)
M-step
Compute wj for Yjr [1,M] M]. Compute (4)).
((TM)
O(TMa²)
Compute (q)).
Onn") (Parlett, 2000)
Compute x² and C.
O m) and (nm)
closed cubic spline. We built the distribution coordinates (EuCOV/GeCOV) in the same manner as the simulation of a loop in IR 10 (appendix C). We then trained the model using the PGPCA EM algorithm. After model training, in each cross validation fold, assessed performance using the 3000 samples in the test set. As our performance measure, we used the data log-likeling bog-likelihood. The average performance was computed the average of log-likeell liglliielssoods over the 15000 samples across the 5 test sets in the 5 cross-validation folds. We compared PGPCA models with two distinct distribution coordinates, Euclidean (EuCOV) and Geometric (GeCOV), We also compared with PPCA. Comparisons between two different learned PGPCA or PPCA models were conducted using paired t-tests on the liglliielsooods of the 15000 test samples.
PGPCA training time. It takes about 13 minutes a regular desktop computer to learn a 10D PGPCA (GeCOV) model with 12000 training samples. This shows that because all steps in Algorithm 1 are analytical, PGPCA EM is efficient in terms of training time, similar to classical EM for linear state-space models (Roweis Ghahramani, 1999). The theoretical computational complexity cach step in one PGPCA EM iteration (Algorithm 1) is listed in Table 4, and every PGPCA EM iteration's computational complexity (TM2) + MTTM²²) Note that PPCA is not iterative (unlike EM), and its computational complexity is (TT?)). Additionally, because the E-step of our PGPCA EM can find the posterior distribution q, (Zz) and our M-step can analytically find the parameters that maximize the ELBO, together our E-step and M-step ensure a monotonic increase ELBO C with each iteration. Consequently, the sequence of L values is monotonically increasing and bounded by L from above, ensuring PGPCA EM convergence by the completeness property real numbers (Rudin al., 1964).
18

---

Published as a conference paper at ICLR 2025
3D probability slices
plane (Axis: y, x vy)
plane (Axis: y, x v,)
XZ plane (Axis: y, x v,)
x10
covariance: GeCOv
p(z) = uniAng
probability
y.
-6
ys
covariance: EUCOV
1.2
plz) uniAng
probability
y.
-6
covariance: GecOV
p(z) = uniTorus
probability
Y.
-6
covariance: EuCOV
piz) = uniTorus
1.2
probability
4
4 -6
6 -6
6 -6
6
Figure 5: Probability distribution p(y) of true models under various manifold latent state probability distributions p(z) and various distribution coordinates (z). For each model, we show three slices (XY, YZ, and XZ planes) that go through the 3D probability distribution for visualization. From the XY plane, it's clear that EuCOV makes p(y) more directional along the axis, and GeCOV is more cylindrically symmetric. Similarly, p(z) uniAng makes p(y) much denser in the inner ring compared to the outer ring, while p(z) uniTorus does not.
Mouse 12
(A) Manifoid: a fiteed loop is 2
IBI
PGPCA over cimonsions
10
Manisid fitted loop
(B)
PSPCA over dimensions
GeCOM Iramadd Rartic boon IP1
Ge00 boin (D1
PRCA
B9C%
5
drandes
Figure 6: PGPCA (GeCOV and EuCOV) better capture the distribution neural firing rates in mice head direction circuit compared with PPCA. (A) and (C) the same fitted loop manifolds in Figure 4A and 4C. (B) PGPCA models consistently has much higher log-likelihood than PPCA across all dimensions. (D) is the same (B) for second mouse, with conclusions being the same.
19


Reviews:
[{"review_id": "brSHKiR6Qc", "reviewer": "Reviewer_VSUw", "summary": "In this work, the author proposes a Probabilistic Geometric Principal Component Analysis (PGPCA) method which can be seen as an extension of PPCA with a better description of data distributions by incorporating a nonlinear manifold geometry. A data-driven EM algorithm is also proposed to solve the PGPCA problem. Experimental results verify the performance of the PGPCA is better than that of PPCA.", "questions": "1. As some works focus on nonlinear PPCA, such as \n\nLawrence, Neil, and Aapo Hyvärinen. \"Probabilistic non-linear principal component analysis with Gaussian process latent variable models.\" Journal of Machine Learning Research 6.11 (2005).\n\nZhang, Jingxin, et al. \"An improved mixture of probabilistic PCA for nonlinear data-driven process monitoring.\" IEEE transactions on cybernetics 49.1 (2017): 198-210.\n\nhow is the performance of the proposed PGPCA compared with these nonlinear PPCA methods?\n\n\n2. How does PGPCA compare with other methods in terms of computational efficiency?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The method and algorithms are technically sound, with reasonable insights and solutions.", "weaknesses": "Comparison with other nonlinear PPCA methods is not provided.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "qlNl39DfUb", "reviewer": "Reviewer_CyYb", "summary": "The paper's main idea is the PPCA extension containing nonlinear manifold terms. The reason for this assumption is not clear. \nBut the formulation and mathematical analysis are true.\nAlso, its performance should be seen in real data. However the experiments on real data are not clear and also they should apply it for more general data like well-known image data sets, or at least mention the class of real data that this method can work well for them in comparison with PPCA. By this version of experiments the importance and quality of the proposed method are not clear.\nAlso, they should report the complexity of the proposed method.", "questions": "_ experiments for real data should be done?\nwhat is the specific class of data that this idea works good?\nwhy they did not report the complexity?\nif they report one can compare them with methods like deep factor models", "limitations": "", "rating": 6, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "The writing is good and also the mathematics are true.", "weaknesses": "The paper's main idea is the PPCA extension containing nonlinear manifold terms. The reason for this assumption is not clear. \nBut the formulation and mathematical analysis are true.\nAlso, its performance should be seen in real data. However the experiments on real data are not clear and also they should apply it for more general data like well-known image data sets, or at least mention the class of real data that this method can work well for them in comparison with PPCA. By this version of experiments the importance and quality of the proposed method are not clear.\nAlso, they should report the complexity of the proposed method.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "qMQ0DJZKsK", "reviewer": "Reviewer_DfNu", "summary": "This paper proposes a generalization of probabilistic principal component analysis (PPCA) to data that lies on a (nonlinear) manifold rather than in Euclidean space around the mean. \n\nThe key contributions of the paper are (i) the introduction of the PGPCA model, (ii) an expectation-maximization algorithm to fit the proposed model and (iii) empirical experiments with real and simulated data to compare PGPCA with vanilla PPCA.", "questions": "- Please add literature on dimensionality reduction for data on a manifold.\n- How well does PGPCA work with other choices of $K(z)$ than geometric and Euclidean?\n- What theoretical guarantees for the hypothesis test on EuCOV vs. GeCOV are there? (e.g. in terms of statistical power/or asymptotic level)\n- How does PGPCA compare to PPCA for data in linear subspaces?\n- How does PGPCA compare against stronger competitors?\n- How well does PGPCA work on other real data?\n- Can the i.i.d. assumptions in model (1) be relaxed? E.g. by allowing (temporal) dependence or slightly varying distributions?\n- Can the assumption of normality be relaxed to allow for non-normal errors?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The problem of dimensionality reduction for data on a manifold is well motivated, and the proposed method closes the gap between PPCA and data on manifolds. \n- The proposed extension PPCA is innovative and well explained.\n- The mathematical derivation of the EM-algorithm is concise and understandable, and the empirical experiments show promising results.", "weaknesses": "- The simulated data seems to be favorable for PGPCA compared to PPCA as it lies on non-linear manifolds. It remains open how well PGPCA compares to PPCA for data in other settings, e.g. in linear subspaces.\n- PGPCA is only compared to PPCA and not other - potentially stronger - baseline methods. Especially, comparisons with other dimension reduction methods for data on manifolds is missing.\n- The evaluation on real data is rather limited. Only brain signals of two mice were used.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Spotlight)