Streaming Heteroscedastic Probabilistic PCA with Missing Data
Authors: Kyle Gilman
David Hong
Jeffrey A Fessler
Laura Balzano
Venue: Accepted by TMLR ()

Abstract: Streaming principal component analysis (PCA) is an integral tool in large-scale machine learning for rapidly estimating low-dimensional subspaces from very high-dimensional data arriving at a high rate. However, modern datasets increasingly combine data from a variety of sources, and thus may exhibit heterogeneous quality across samples. Standard streaming PCA algorithms do not account for non-uniform noise, so their subspace estimates can quickly degrade. While the recently proposed Heteroscedastic Probabilistic PCA Technique (HePPCAT) addresses this heterogeneity, it was not designed to handle streaming data that may exhibit non-stationary behavior. Moreover, HePPCAT does not allow for missing entries in the data, which can be common in streaming data. This paper proposes the Streaming HeteroscedASTic Algorithm for PCA (SHASTA-PCA) to bridge this divide. SHASTA-PCA employs a stochastic alternating expectation-maximization approach that jointly learns the low-rank latent factors and the unknown noise variances from streaming data that may have missing entries and heteroscedastic noise, all while maintaining a low memory and computational footprint. Numerical experiments demonstrate the superior subspace estimation of our method compared to state-of-the-art streaming PCA algorithms in the heteroscedastic setting. Finally, we illustrate SHASTA-PCA applied to highly heterogeneous real data from astronomy.

Published in Transactions on Machine Learning Research (08/2025)
Streaming Heteroscedastic Probabilistic PCA with Missing Data
Kyle Gilman
kgilman@umichhcoo edu
Deportment of Electrical Engineering and Computer Science University of Michigan
David Hong
hong & udel. edu
Department Electrical and Computer Engineering University of Delawore
Jeffrey A. Fessler
fesster@umich.. edu
Department of Electrical Engineering and Computer Science University of Michigan
Laura Balzano
girasole@umich.com edu
Department of Electrical Engineering and Computer Science Universily of Michigan
Reviewed on OpenReview: https: // optp//veenrre...e/ net/ forum? id= 1b2rPLuP9X
Abstract
Streaming principal component analysis (PCA) is an integral tool in large-scale machine learning for rapidly estimating low-dimensional subspaces from very high-dimensional data arriving at a high rate. However, modern datasets increasingly combine data from a variety of sources, and thus may exhibit heterogeneous quality across samples. Standard stream- ing PCA algorithms do not account for non-uniform noise, their subspace estimates can quickly degrade. While the recently proposed Heteroscedastic Probabilistic PCA Technique (HePPCAT) addresses this heterogeneity, it was not designed to handle streaming data that may cxhibit non-stationary behavior. Morcover, HcPPCAT does not allow for miss- ing entries in the data, which can common in streaming data. This paper proposes the Streaming HeteroscedASTic Algorithm for PCA (SHASTA-PCA) bridge this divide. SHASTA-PCA employs stochastic alternating expectation-maxxiizzaaoo approach that jointly learns the low-rank latent factors and the unknown noise variances from streaming data that may have missing entries and heteroscedastic noise, all while maintaining low memory and computational footprint. Numerical experiments demonstrate the superior subspace estimation of our method compared to statc-of-tho-art streaming PCA algorithms the heteroscedastie selting. Finally, we illustrate SHASTA-PCA applied highly hel- erogeneous real data from astronomy.
1 Introduction
Modern data are increasingly large in scale and formed by combining heterogeneous samples from diverse sources or conditions that exhibit heteroscodastic noise, or noises of different variances (Hong al., 2021),
Code for project be found K. Gilman & Balzano supported part NSF CAREER award CCF-1845076 and AFOSR YIP award FA9550- 1-0026. J.A. Fessler L. Balzano supported part by NSF award CCF-2331590. Hong supported by NSF
Mathematical Sciences Postdoctoral Research Fellowship DMS 2103353.

---

Published in Transactions on Machine Learning Research (08/2025)
such as in astronomy (Ahumada et al., 2020), medical imaging (Pruessmann et al., 1999; Anam et al., 2020); and cryo-electron microscopy imaging (Andén & Singer, 2017; Bendory et al., 2020). Principal component analysis (PCA) for visualization, exploratory data analysis, data compression, prediction, or other down- stream tasks is often a fundamental tool to process these high-dimensional data. However, several practical challenges arise when computing PCA on these types of data. In many applications, due to memory or phys- ical constraints, the full data cannot be observed in their entirety at computation time and are instead read partially into memory piece by piece, or observations may stream in continuously and indefinitely. More- over, the signal may evolve over time and require adaptive tracking algorithms for the low-rank component. Adding to these difficulties, it is also common for big data to contain missing entries, such as in magnetic resonance imaging (Mensch et al., 2017), collaborative filtering (Candes & Plan, 2010), and environmental sensing (Ni et al., 2009). Consequently, there is a need for scalable streaming PCA techniques that can handle heteroscedastic noise and missing data.
A tremendous body of work has studied streaming PCA techniques for learning a signal subspace from noisy incremental data observations that have missing entries. Streaming or online PCA algorithms oflen enjoy the advantages of computational efficiency, low memory overhead, and adaptive tracking abilities, making them very useful in real-world big-data applications. However, no existing streaming methods account for noise with differing variances across samples, i.e., sample-wise heteroscedastic noise, and their subspace estimates can be highly corrupted by the noisiest samples. Hong et al. (2021) developed a Heteroscedastic Probabilistic PCA technique (HePPCAT) for data with varying noise levels across the samples. НePPCAT estimation, but only
learns the low-rank factors and the unknown noise variances via maximum-likcllhoon in the batch setting with no missing entries. Other batch hederoscedlastic PCA algorithms, like weighted PCA studied in Jolliffe (2002); Young (1941); Hong et al. (2023) and HetcroPCA (Zhang et al., 2022) for data with heteroscedastic features, also lack streaming and adaptive tracking abilities. None of the existing methods handle all of the real-data complexities we study here: missing entries, heterosedastic noise, and streaming data. Tackling this non-trivial setting requires developing new algorithms.
To the best of our knowledge, this paper is the first work to develop a streaming PCA algorithm for data with missing entries and heteroscedastic noise. Our algorithm jointly estimates the factors and unknown noise variances in online fashion from strcaming incomplete data using alternating stochastic minorize- maximize (SMM) approach with small computational and memory overhead. We carefully design minorizers with a particular alternating schedule of stochastic updates that distinguishes our approach from existing SMM methods. Notably, handling missing entries and heteroscedastic noise involves more complex updates than the simpler algebraic formulations of algorithms like HºPPCAT or PETRELS (Chi et al., 2013). We demonstrate that our algorithm can estimate the signal subspace from subsampled data (even without know- ing the noise variances) better than state-of-the-art streaming PCA methods that assume homogeneous noise. Our algorithm is unique in that it not only tracks low-dimensional dynamic subspaces, but can also track dynamic noise variances that can occur, e.g., in sensor calibration (Jun-hua et al., 2003) and beamforming in nonstationary noise environments (Cohen, 2004). Finally, this work extends our understanding of streaming PCA to the setting of heteroscedastic noise and draws comnections to existing work in the literature. In par- ticular, the proposed method closely relates to the streaming PCA algorithm PETRELS (Chi et al., 2013), but differs by learning the noise variances of the data on-the-fly and downweighting noisier data samples in the factor updatos. Furthermore, we show that our proposed method implicitly optimizes a regularizod least-squares problem whose adaptive hyperparameter varies by the learned heterogeneous noise variances.
Section 2 discusses related works for both streaming PCA and heteteroscedastic PCA. Sections 3 to 6 describe the model we consider, define the resulting optimization problem, and derive the proposed algorithm. Section 7 presents synthetic and real data experiments that demonstrate the benefits of the proposed method over existing state-of-the-art streaming PCA algorithms.
1.1 Notation
We use bold upper-case letters A to denote matrices, bold lower-case letters v to denote vectors, and non- bold lower-case letters e for scalars. We denote the Hermitian transpose of a matrix as A' and the trace of a matrix as tr(.). The Euclidean norm is denoted by || ll2. The identity matrix of size d x d is denoted Ia- The notation i [k] means E (1,...,k).
2

---

Published in Transactions on Machine Learning Research (08/2025)
2 Related work
2.1 Streaming PCA
A rich body of work has investigated a variety of streaming PCA algorithms for learning signal subspace from incremental and possibly incomplete data observations. However, these methods assume that the data have homogeneous quality and do not model heteroscedastic properties like those considered in this paper. Since there are too many related works to detail here (scc, e.g., Balzano et al. (2018); Boutsidis et al. (2016), for related references), we highlight a few of the most related.
One prominent branch of algorithms uses stochastic gradients to update the learned subspace based on a new data obscrvation at each iteration; see, e.g., Bertsekas (2011) and Bottou (2010). Mardani et al. (2015) use stochastic gradient, descent to learn matrix and tensor factorization models in the presence of missing data and also include an exponentially weighted data term that trades off adapting to new data with fitting historical data. Stochastic gradient descent over Riemannian manifolds is also a popular approach; see, e.g., Bonnabel (2013); Balzano et al. (2010); He et al. (2012); Goes et al. (2014). Oja's method (Oja, 1982) takes a stochastic gradient step to update the subspace basis from the most recent data vector, obtaining a new orthonormal basis after orthogonalization. Balzano (2022) proved the equivalence of the GROUSE algorithm (Balzano et al., 2010) with Oja's method for certain step size. Other progress on improving and understanding Oja's method has recently been made, such rate in one pass over the data (Henriksen & Ward, 2019) and an analysis of the convergence rate for non-i.i.d. data sampled from a Markov chain (Kumar& & Sarkar, 2024).
algorithm to adaptively select the learning
Some streaming PCA methods share commonalities with quasi-second optimization methods. For example, the PETRELS algorithm proposed in Chi et al. (2013) fits a factor model to data with missing entries via a stochastic quasi-Newton method. PETRELS has computationally efficient. updates but can encounter numerical instability issues in practice after a large number of samples.
More recently, streaming PCA and its analogs have been extended to a variety of new problems. Giannakis et al. (2023) propose a streaming algorithm for forecasting dynamical systems that they show is a type of streaming PCA problem. Streaming algorithms have been proposed for robust PCA (Thanh et al., 2021; Diakonikolas et al., 2023), federated learning and differential privacy (Grammenos et al., 2020), a distributed Krasulina's method (Raja& & Bajwa, 2022), and probabilistic PCA to track nonstationary processes (Lu al., 2024). Although not strictly a streaming algorithm, Blocker al. (2023) also estimate dynamic subspaces but by using piecowisc- gicoosise-coodesic model on the Grassmann manifold. All of these approaches implicitly assume that the data quality is homogeneous across the dataset, in contrast to our proposed approach.
2.2 Stochastic MM methods
Another vein of work on streaming algorithms, which has the closest similarities to this paper, is stochastic pneriiatiionmmiimmzaate (SMM) algorithms for matrix and tensor factorization. MM methods construct surrogate functions that are more casily optimized than the original objective; algorithms that successively optimize these surrogates are guaranteed to converge to stationary points of the original objective function under suitable regularity conditions (Jacobson & Fessler, 2007; Lange, 2016). SMM algorithms such as Mairal (2013); Lyu et al. (2022); Lyu (2024) update and optimize a stochastic approximation of the surrogate function upon observing each new data sample. The work in Mairal (2013) proved almost sure convergence to stationary point for non-convex objectives with one block of variables for the SMM technique. The work in Mensch et al. (2017) proposes subsampled online matrix factorization (SOMF) for large-scale streaming subsampled data and gives convergence guarantees under mild assumptions. In Lyu et al. (2022); Lyu (2024), the authors extend SMM to functions that are multi-convex in blocks of variables for online tensor factorization. Their framework performs block-coordinate minimization of single majorizer at cach time point. They prove almost sure convergence of the iterates to a stationary point assuming certain regularily conditions of the loss function and learning rate and that the data tensors follow a Markovian process. and the data tensors form a Markov chain with a unique stationary distribution. More recently, Phan et al. (2024) propose and analyze several stochastic variance-reducee MM algorithms. Like these works, our paper also draws upon SMM techniques; we use an alternating SMM approach to optimize the log-likelihood function
3

---

Published in Transactions on Machine Learning Research (08/2025)
for our model. However, our setting and approach differ from these existing methods in key ways; we detail these differences and their impact on the related convergence theory in $6.3.1.
A close analogue to SMM is the Doubly Stochastic Successive Convex (DSSC) approximation algorithm (Mokhtari & Koppel, 2020) that optimizes convex surrogates to non-convex objective functions from stream- ing samples or minibatches. A key feature of their algorithm that decomposes the optimization variable into B blocks and operates on random subsets of blocks at each iteration. Specifically, the DSSC algorithm chooses block i [B], computes stochastic gradients with respect to the ith block of variables and then recursively updates the approximation to the ith surrogate function. From the optimizer to the approximate surrogate, their algorithm performs momentum updates of the iterates very similarly to SMM. Our own al- gorithm SHASTA-PCA in $6 can be interpreted as following a similar approach, but without using gradient methods since our problem does not have Lipschitz-continnov gradients.
2.3 Heterogeneous data
Several recently proposed PCA algorithms consider data contaminated by heteroscedastic noise across sam- ples, which is the setting we study in this paper. Weighted PCA is a natural approach in this context (Jolliffe, 2002), either weighting the samples by the inverse noise variances (Young, 1941) or by optimal weighting derived in Hong et al. (2023). In both instances, the variances must be known priorí esti- mated to compute the weights. Probabilistic PCA (PPCA) (Tipping & Bishop, 1999b) uses a probabilistic interpretation of PCA via a factor analysis model with isotropic Gaussian noise and latent variables. For a single unknown noise variance (i.e., homoscedastic noise), the learned factors and noise variance are solu- tions to a maximum-likclihond problem that can be optimized using an expectation maximization algorithm; these solutions correspond exactly to PCA. Hong et al. (2021) studied the heteroscedastic probabilistic PCA (HPPCA) problem that considers a factor model where groups of data may have different. (unknown) noise mexximum-likelihon
variances. Their method, HePPCAT, performs unknown noise variances (assuming knowledge of which samples belong to cach noise variance group); they consider various algorithms and recommend an alternating EM approach. Other batch heteroscedastic PCA methods have since followed HcPPCAT. ALPACAH (Salazar Cavazos et al., 2025) estimates the low-rank component and variances for data with sample-wise heteroscedasttcc noise, but no streaming counterpart currently exists. HeMPPCAT (Xu et al., 2023) extends mixtures of probabilistic PCA (Tipping & Bishop, 1999a) to the case of heterosccddastic noise across samples.
estimation of the latent factors and
More broadly, there is an increasing body of work that investigates PCA techniques for data contaminated by some sort of heterogeneous noise, including noise that is heteroscedastic across features. HeteroPCA (Zhang et al., 2022) iteratively imputes the diagonal entries of the sample covariance matrix to address bias that arises when the noise has featuro-wisc heteroscodasticityy and Zhou & Chen (2025) then extended HeteroPCA to the case of ill-conditioned low-rank data. Another line of work in Leeb & Romanov (2021) and Leeb (2021) has considered rescaling the data instead whiten the noise. Yan et al. (2024) develops inference and uncertainty quantification procedures for PCA with missing data and fcature-wise hetcroscodasticityy There has also been recent progress on developing methods to estimate the rank in heterogeneous noise contexts (Hong et al., 2020; Ke et al., 2021; Landa et al., 2022; Landa & Kluger, 2025) and on establishing fundamental limits for recovery in these settings (Behne & Reeves, 2022; Zhang & Mondelli, 2024).
A closely related problem in signal processing applications, such as heterogeneous clutter in radar, is data with heterogencous "textures", also called the "mixed effects" problem. Here, the signal modeled as a mixture of scaled Gaussians, each sharing a common low-rank covariance scaled by an unknown deterministic positive "texture" or power factor (Breloy et al., 2019). In fact, the heterogencous texture and HPPCA problems are related up to an unknown scaling (Hong et al., 2021). Ferrer et al. (2021) and Hippert-Ferrer et al. (2022) also studied variations of the heterogencous texture problem for robust covariance matrix estimation from batch data with missing entries. Collas et al. (2021) study the probabilistic PCA problem in the context of isotropic signals with unknown heterogeneous textures and a known noise floor. Their paper casts the maximum-likelioon estimation as an optimization problem over a Riemannian manifold, using gradient descent on the manifold to jointly optimize for the subspace and the textures. Their formulation also readily admits stochastic gradient algorithm for online optimization.
4

---

Published in Transactions on Machine Learning Research (08/2025)
Heteroscedastic data has also been investigated in the setting of supervised learning for fitting linear re- gression models with stochastic gradient descent (Song et al., 2015). The authors show that the model's performance given "clean" and "noisy" datasets depends on the learning rate and the order in which the datasots are processed. Further, they propose using separate learning rates that depend on the noise levels instead of using one learning rale as done in classical SGD. In the context generalized linear bandits, Zhao et al. (2023) propose an online algorithm for the heteroscodastic bandit problem using weighted lincar regression with weights selected as the inverse noise variance.
3 Probabilistic Model
Similar to Hong et al. (2019; 2021), we model data samples in R from L noise level groups as:
y: Fz; + Es.
for i = 1,2,
(1)
where F Bd is deterministic factor matrix estimate, z; NO Is) are independent and identically distributed (i..d.)) coefficient vectors, E; ~ NOa,, Ug. Ia) are independent noise vectors, 9: E fl,..., L} is the noise level group to which the ith sample belongs, and v1,-..,Uf, Ür. are deterministic noise variances to estimate. Typically, we assume & < to model data from a low-dimensional subspace. We also assume the group memberships 9i are known. This a valid assumption some real applications where we know the sources of the data, e.g., the sensor type in air quality monitoring (Hong et al., 2021), data from a particular coil in the MRI machine.
Let Ω C (1....,d) denote the set of entries observed for the ith sample, and let yo, R and Fo, E R |0, xk denote the restrictions of and F to the entries and rows defined by Ωg Then the observed entries of the data vectors are distributed as
yo, ~ NO Fo, Fo, + Fg. I..1))
We will express the joint log-likelihood over only the observed entries the data and maximize it for the unknown deterministic model parameters.
For batch of n vectors, the joint log- ligelikelhood over the observed batch data for Ω (S2,..., 2n) can be easily written in incremental form of log-likelihoods over the partially observed dataset Ya e (Ya)...., Vr.)):
c(Y;; F,
2 Σ Cy F, + C,
(2)
where v = =Iu,,..,, vnl', C a constant independent of and v, and
C, F, v) A det Fh, + " Ta.)) váa (Fo, F/, + rg yo.
(3)
is the loss for a single vector You- To jointly estimate the factor matrix and the variances v, we maximize this likelihood. Optimizing the log-likelihood (2) is a challenging non-concave optimization problem, so we propose an efficient alternating minorize-maximin (MM) approach.
4 Expectation Maximization Minorizer
This section derives minorizer for the log-likelihood (2) that will be used develop the proposed alternating MM algorithm the following sections. In particular, we derive a minorizer the point (F, 6) the style of expectation maximization methods. The minorizer follows from the work in Hong al. (2021); here we extend it to the case for data with missing entries.
5

---

Published in Transactions on Machine Learning Research (08/2025)
For the complete-data log-likelhoodd, the observed samples yn, and unknown coefficients 2;, leading to the following complete-data log-likelihood for the ith sample:
( v) _ In pynn, 2; F,
Inp |2 F. + In inp(zi F.
22
ynn. Fs, zi
1
In Vn.
(4)
20g:
where (4) drops the constants In(24)-1uu12² and
Next, we take the expectation of (4) with respect to the following conditional distribution x zly, derived in Appendix B:
ind N(E((F 8g (((, 6)), z; (yn, F F v=
where we define & and M, for use here and in following derivations as:
(FF & M,F, v) F) Yo. (. v) 2 (F Fo. + + Ix)-
(5)
(6)
Doing so yields the following minorizer for L, at (F, T):
1yn, 22
(FFFFFFFF8) F) 2
In yg
yn, Fs. 2 (F.
2
20g.
Ug.
(7)
Fa,
0)11
8g tr Fo (F..)))
20g:
where (7) drops terms that are constant with respect to F and v.
5 A Batch Algorithm
Before deriving the proposed streaming algorithm, SHASTA-PCA, we first derive batch method for com- parison purposes. Summing the sample-wise minorizer (7) across all the samples gives the following batch minorizer the point (F,0):
V U.
V.FFUu
(8)
12
11 11
In
Va. Fa, (F,
2
20
A= i:gimt
(F, & Fa, (, 5)),
2n
Similar to HePPCAT (Hong et al., 2021), which is a batch method for fully sampled data, in each iteration i, we first update v (with F fixed F(-1) then update F (with v fixed at v,), ie.,
v. arg max (F F2-1,0-111
(9)
F arg max (F ;y; F.11.
(10)
Here refers to index of the data sample, and t denotes only the algorithm iteration, in contrast to the strcaming algorithm $6, where t denotes both the time index (sample) and algorithm iteration. The following subsections derive efficient formulas for these updates and discuss the memory and computational costs.

---

Published in Transactions on Machine Learning Research (08/2025)
5.1 Optimizing v for fixed F
Here we derive an ellicient formula for the update in (9). While the update is similar to HePPCAT (Hong et al., 2021), the key difference lies in computing the minorizer parameters. Specifically, the missing data makes the samplo-wisc quantitics ())))) and Mc in (5) and (6) depend on the sampling patterns for the ith data vector and must be computed for every sample i E [n] per iteration compared to the single & and M used in HcPPCAT. This update separates into L univariate optimization problems, one in each variance Uy:
hatch
Pt,
Vt, = max 1/
In dE
2
2ve
where
Obatch
i:
batch PE.
|130. F--12, 02-1) - Ran (. p(-1))
::9:=t
F1-10 denotes the iterate F.-1 restricted to the defined by Q. The corresponding solutions
batch
Pt.
0t Chatch
We precompute @batch because it remains constant across iterations.
5.2 Optimizing F for fixed
Here derive efficient formula for the F update in (10). The update differs from the factor update in HcPPCAT again due to the missing entrics in the data. Specifically, this update separates into d quadratic optimization problems, one in each row fj of F:
batch 1 Bhatch fj+
St arg mag 81 2
f
where
Rbatch A
()) (t) (F
Ut.
8=1 iige
0,31
(t)
batch A
Σ
Mij
(F-.1,
Ft.t
8=1
0,31
capture the data-dependent terms in the minorizer in (8), and Yij the jth coordinate of the vector y. We compute the solutions for fy in parallel
fij (Rlatch)
-1 batch
St.
Vj E [d].
5.3 Memory and Computational Complexity
The batch algorithm above involves first accessing all Σl_ ny data vectors to compute the minorizer parameters M.(..) E R and 2,(::)) E R at a cost of Onk + E 24(k3) (n(122 + k)) memory elements. Computing R and batch incurs a cost of O En |0, R)) t.j
flops per iteration and
flops for all j = 1, 1....,d d, and finally solving for the rows of F costs (dd:3) flops per iteration. Updating v requires O Exgg,e (Q,122) computations.

---

Published in Transactions on Machine Learning Research (08/2025)
Since each complete update depends on all n samples, the batch algorithm must have access to the entire dataset at run-time, either by reading over all the data in multiple passes while accumulating the computed terms used to parameterize the minorizers, or by storing all the data at once, which requires OiEi= [22,] + dk²) memory. This requirement, combined with the ((n) inversions & x k matrices in each iteration, significantly limits the practicality of the batch algorithm for massive-soale or high-arrival-rate data as well as in infinite-streaming applications.
6 Proposed Algorithm: SHASTA-PCA
The structure of the log-likelihodd in (2) suggests natural way to perform incremental (in the finite data setting) or stochastic (in expectation) updates. If each data sample from the 9th group is drawn i.i.d. from the model in (1) then under uniform random sampling of the data entries, each L; is an unbiased estimator of f. Hence, we leverage the work in Mairal (2013), which proposed a stochastic MM (SMM) technique for optimizing empirical loss functions from large-scale or possibly infinite data sets. For the loss function we consider, these online algorithms have recursive updates with a light memory footprint that is independent of the number of samples.
In the streaming setting with only a single observation yn, at cach time index t, we do not have access to the full batch minorizer in (8), but rather only a single (F u; F. ö). Key to our approach, for each t, our pro- posed algorithm uses (F v; F, 8) to update two separate approximations to (FF, u; F, 8) parameterized by F and v, respectively, i..,, ( (F) and (v), in alternating way. While other optimization ap- proaches are possible for example, performing block coordinate maximization of single joint approximate majorizer our novel approach of alternating between the two separate approximate minorizers reduces mem- ory usage and computational overhead in our setting. Given a sequence of non-increasing positive weights fwe E (0,1) and positive scalars Co and CF, we first update v (with F fixed at Fi-1) with one SMM iteration, then update (with v fixed at vx) with another. Namely, we have the noise variance update:
(v),
T (w) -1- w.) w.)b(( (v) + w, (F(FF-1, F 1t-1);
(11)
(1 -1-ca (1-c,)or- cu arg max 700) (v),
(12)
followed by the factor update:
(F) = (1 w.) w.)) we V, ve: F.) vz)
(13)
F, = (1 CF) F F + OF arg max T (F).
(14)
F
The iterate averaging updates in (12) and (14) are important to control the distance between iterates and have both practical and theoretical significance in SMM algorithms (Lyu et al., 2022; Lyu, 2024). Empirically, we found that using constant er and o, worked well, but other iterate averaging techniques are also possible, such as those discussed in Mairal (2013). Other ways to control the iterates include optimizing over a trust region, as done in Lyu al. (2022).
Since the iterate and the time index the same in the streaming setting, i.e., t= i, we now denote both the sample and the SMM iteration by in the remainder of this section. We now derive efficient recursive updales and compare the memory and computational costs to the batch algorithm.
6.1 Optimizing e for fixed F
Now note that the minorizer in (7) for i = 1 and 0) = (F.-1,0--) can be re-written
V. F--1, Ct - 12p In In
Pt.
Uge
where C₁ does not depend on a and
Pt P yo. F1-1. -
01 tr( tr(Fi-,,0, F (F,1,,,
8

---

Published in Transactions on Machine Learning Research (08/2025)
Recall that 9r [L] is the group index of the tth data vector. Thus, it follows from Equation (11) that
(v) C2 De
PL.
where C is constant that does not depend on
or. - (1 - We 0c 1.00 + we (9tl
(15)
Pt.ge - (1 - w/ w, 1.9 + WpPt.
(16)
and for & # 9t
02, - (1 Wi ww)0:-11. Pe, = (1 we
(17)
The fth term in the sum is optimized by pt,e/0e,e
Pr.
Wr = (1 Cvp Ut-
1,e
Co
(18)
Here the vectors 0, E RL and P: E RL aggregate past information to parameterize the approximate minorizer in v. Since there is no past information at t = O, we chose to initialize them with zero vectors. However, in the initial iterations where no data vectors have been observed for the fth group, (18) is undefined, so a valid argument maximizing (12) is simply 00,6, ie., the initialized value.
6.2 Optimizing F for fixed v
Maximizing the approximate minorizer V (FF with respect to F reduces to maximizing quadratics in the rows of F for j E [d]:
d
1 (F) Σ f f,
(19)
where for j S
Re,j = (1 - wi R + we Re.j.
(20)
Btj = (1 - we 3 1, + wr
(21)
where
( 1 -) + (F
Re,
A 2 Ut,g
8tJ
2
Mto = ve),
(22)
and for j st
R (1 -w--ww)) W. Rt. (1 w 1.j-
(23)
The parameters (Re,, Bt,) for E [d] of the approximate minorizer / (F) aggregate past information from previously observed samples, permitting our algorithm to stream over an arbitrary amount of data while using a constant amount memory.
Maximizing the approximate minorizer w ((()(F)
(F) with respect to each row of F yields
f, = R i 1....d d.
(24)
Since the problem separates in each row of F, this form permits efficient parallel computations. Further, because f, R $t, = R 3t-1, for j & Ωe we solve the linear systems in (24) only for the rows indexed by j 9r. After obtaining the candidate iterate F above, the final step updates F, by averaging in (14).
9

---

Published in Transactions on Machine Learning Research (08/2025)
Algorithm 1: SHASTA-PCA
Input: Rank h, weights (wi) E (0, 1], parameters cp, Cu > O, initialization parameter 8 8> 0.
Data: Ui....,yr]] Yt E IRd, group memberships Ir E 21........ ... L. for all i, and sets of observed indices (A,,..., St), where Ω S (I,...,d)
Output: E Rd v _
Initialize Fo and vo via random initialization;
2 Initialize surrogate parameters R ôl for and SI,, 02 for j E [d]. 6g =Pu =0
for 1= ,,...,? do
Fixing F at Fi-1, Compute 0, and P. via (15)-(17).
2. Compute from (18).
Fixing at Us,
1. Update Reg and Bry via (20)-(23).
2. Compute F via (24) in parallel.
3. F. = (1 Fr-1 F.
6.3 Algorithm and Memory/Computational Complexity
Together, these alternating updates form the Streaming HeteroscodASTic Algorithm for PCA (SHASTA- PCA), detailed in Algorithm 1.
The primary memory requirement of SHASTA-PCA is storing d 1 many k x k matrices and k-length vectors for the surrogate parameters and two additional L-length vectors for the v parameters. Thus, the dominant memory requirement of SHASTA-PCA is (d + k)) memory elements throughout the runtime, which is independent of the number of data samples.
The primary sources of computational complexity arise from: i) forming Mc(..) at cost of 0(19/k2 + k²) flops, ii) computing (,,-) at a cost of 0(12, flops, iii) forming -_= ) (.))))) + M,(,,) at a cost of Ft.gs
O(k²) flops, iv) computing Pt. at cost of O k22 + k²) flops when using an efficient implementation with matrix-vector multiplications, and v) updating F; at a cost of 0(12, (R. + R²)) for the multiplications and inverses. total, each iteration of SHASTA-PCA incurs O1 + +²)) flops.
As discussed below, PETRELS (Chi et al., 2013) uses rank-one updates the pseudo-inverses of the matrices in (26) to avoid computing new pseudo-inverse each iteration, but that approach does not apply in our case since the updates to Ru in (20) are not rank-one. Still, updating F only requires inverting [St many k x k matrices each iteration, which remains relatively inexpensive since k < d and is often small in practice. Note that the complexity appears to be the worst for [2p] = d with the implementation described above, but in this setting, since all the surrogate parameters are the same, one can verify that only single surrogate parameter R, (and its inverse) is necessary. In realily, the worst-case complexity happens when [2,] 1, i.e., when single entry per column missing. Identifying an approach with improved computational complexity in the highly-sampled setting remains an interesting future research direction.
6.3.1 Convergence
Empirically, we observe that SHASTA-PCA converges to a stationary point as the number of samples grows. Several factors influence how fast the algorithm converges in practice, Similar to stochastic gradient descent, the choices of weights (w;) and iterate averaging parameters OF and e, affect both how fast the algorithm converges and what level of accuracy it achieves. Using larger weights tends lead faster convergence but only to within larger, suboptimal local region of an optimum. Conversely, using smaller weights tends
10

---

Published in Transactions on Machine Learning Research (08/2025)
to lead to slower progress but to tighter region around an optimum. The amount of missing data also plays a major role. A higher percentage of missing entries generally requires more samples more passes over the data to converge to an optimum..
Our setting differs in several key ways from the prior works discussed in $2 that establish convergence for SMM algorithms. First, our minorizers are neither Lipschitz smooth nor strongly concave (in fact, the minorizer for vg is nonconcave), so the theory in Mairal (2013) and Mensch et al. (2017) does not directly apply. Second, our algorithm maximizes the log-likelihood in two blocks of variables, and does in an alternating fashion with two separate approximate minorizers (F)) and (v), which distinct from the work Lyu et al. (2022) that alternates updates over the blocks of a single joint approximate minorizer. Notably, such as the case of v, we update an aggregation ( (v) of the restricted minorizers V. (F v; Fi v;-1):
(e) (v) (1 - (1-w)) (2) (v) + w. (F Fi-1, v1-1)
W (1 wy))
(25)
The dependence of 1 (v) on all past iterates {F in the first argument of each V, in (25) precludes using the analysis of Lyu et al. (2022) for BCD of single approximate minorizer each ileration. Likewise, the use of this restricted minorizer minorizes the objective with respect to v with F fixed at Fş-1, which breaks, e.g., the induction argument used in Mairal (2013, Proposition 3.3). To our knowledge, existing work establishes convergence for this setting.
(1)
()) precludes using the analysis of Mairal (2013) since only
However, we conjecture that similar convergence guarantees are possible for SHASTA-PCA since both (v) and (v) retain many of the same MM propertics used to analyze SMM algorithms. It is likely that theory of convergence to a stationary point would require adjustments to our algorithm, such as the addition of constraints on F and v to provide some control over the curvature of the objective. Since the convergence results possible for non-convex problems using stochastic optimization are typically weak, we leave the proof of convergence to future work.
6.4 Connection to recursive least squares and HePPCAT
The SHASTA-PCA update for the factors in (24) resembles recursive least squares (RLS) algorithms like PE- TRELS (Chi et al., 2013). The objective function considered in Chi et al. (2013) is, in fact, the maximization of our complete log-likelihood in the homoscodastic setting with respect to the factors F and latent variables 2, without the & penalty on z, in (4). PETRELS first estimates the minimizer to 2, via the pseudo-inverse solution and then updates each row f, by computing (24) using similar updates to the minorizer parameters:
for j jEB. Q Raj = AR-11 + 2,2//
J.j = \5 + De.j
for j & Ω: Ri, = AR-11.j
. = \m
where &, = FF 1,,, Yn, and A E (0, 1) is a forgetting factor that exponentially downweights the importance of past data. In the stochastic MM framework, wi plays an analogous role to À by exponentially down-weighting surrogates constructed from historical data.
However, there are some important differences. Here, the complete data log-likelihood effectively introduces Tikhonov regularization on že, where the Tikhonov regularization parameter is learned by estimating the the weights,
noise variances. The PETRELS objective function can similarly incorporate regularization but with uscr-specified hyperparameter. It is well known that the appropriate hyperparameter in Tikhonov regularization depends on the noise variance of the data (O'Leary, 2001; Cao et al., 2021). Here, SHASTA- PCA implicitly learns this hyperparameter as part of the maximun-likelihood estimation problem for the unknown heterogeneous noise variances.
PETRELS can also be thought of as a stochastic second-order method that quadratically majorizes the function in F at each time l using the pseudo-inverse solution of the weights given some estimate F-11-
11

---

Published in Transactions on Machine Learning Research (08/2025)
Our algorithm optimizes similar quadratic majorizer in F for each t. While the iseudo-inverse solution for maximizing the complete data log-likehood in (4) with respect to 2t. equivalently the conditional mean of 2;, appears in the update of F through žt, the update additionally leverages the covariance of the latent variable's conditional distribution and, perhaps most importantly, an inverse weighting according to the learned noise variances that downweights noisier data samples.
SHASTA-PCA also has connections to the HePPCAT algorithm (Hong al., 2021). Indeed, the SHASTA- PCA updates of F and closely resemble- can be interpreted as stochastic approximations to- HePPCATTs EM updates of F and in Hong et al. (2021, equ. (8)) and Hong et al. (2021, eqn. (15)). respectively. More preciscly, cach j, approximates each column of
Z
Vt,E
of Hong et al. (2021, eqn. (8)), where Zrd Mi,FF.Y, and M., A (F(F. + v.ll)) and each Ra. approximates the matrix
+ -nMMe..
However, each of these terms in SHASTA-PCA depends on the observed data coordinate in the update of the corresponding row of F. Since each row of depends on a different Reg for each j E [d] due to missing data, we cannot use the SVD factorization of F; to expedite the inverse computation as in HePPCAT (Hong et al., 2021, eqn. (9)). Other minorizers for the update that were considered in Hong et al. (2021), such as the difference of concave, quadratic solvable, and cubic solvable minorizers, may also have possible stochastic implementations. We leave these possible approaches to future work since they did not appear to result in more efficient updates.
7 Experimental Results
7.1 Incremental computation with static subspace
This section considers the task of estimating a static planted subspace from low-rank data corrupted by heterogeneous noise. We generate data according to the model in (1) with ambient dimension 100 from a rank-3 subspace with squared singular values 7 = [4, 2, 1]. drawing 500 samples with noisc variance 10-2 and 2,000 samples with noise variance 10-1 We draw an orthonormal subspace basis matrix U E R R uniformly at random from the Stiefel manifold, and set the planted factor matrix to be F U diag(v/))
After randomly permuting the order of the data vectors, we compare SHASTA-PCA to PETRELS and the streaming PCA algorithm GROUSE (Bakano et al., 2010) (which has recently been shown to be equivalent to Oja's method (Oja, 1982) in Balzano (2022)) that estimates a subspace from rank-one gradient steps on the Grassmann manifold, with a tuned step size of 0.01. SHASTA-PCA jointly learns both the factors F and noise variances e from each streaming observation. For SHASTA-PCA, we use w; = 1/t (where the time index), Cp = G. = 0.1 and initialize the parameters RG 5I with ò 0.1 for both SHASTA-PCA and PETRELS. We initialize each streaming algorithm with the same random Fo, and each entry of for SHASTA-PCA uniformly at random between 0 and 1. We set the forgetting parameter in PETRELS to A = 1, corresponding to the algorithm's 's batch mode. As a baseline, we compare to batch algorithms for fully-observed data: HcPPCAT (Hong et al., 2021) with 100 iterations, which we found to be sufficient for convergence, and homoscedastic probabilistic PCA (PPCA) (Tipping & Bishop, 1999b) on the full data. In addition, we compute PPCA on each data group individually, denoted by "G1" ("G2") in the legend of Fig. la corresponding group (2) with 500 (2,000) samples with noise variances 10-2 (10-1) respectively.
The first experiment in Fig. la compares each algorithm in the fully observed data setting, where the streaming algorithms compute F and v incrementally using single vector in each iteration. Given the
1 All experiments were performed in Julia on 2021 Machook Pro with the Apple MI Pro processor and 16 GB of memory. We reproduced and implemented all algorithms ourselves from their original source works.
12

---

Published in Transactions on Machine Learning Research (08/2025)
planted model parameters F' and vr and their log- log-likelihood value = LFFF ,0),, the left plot in Fig. la shows the normalized log-likelihood C(F,, - with respect to the full dataset in (2) for each iteration of SHASTA-PCA compared to the batch algorithm baselines, Because GROUSE and PETRELS do not estimate the noise variances, we omit them from this plot. The right plot Fig. la shows convergence of the iterates with respect to the normalized subspace error - for the estimate Ü, RRdk of the planted subspace U; for SHASTA-PCA and PETRELS, orthogonalizing F. Each figure plots the mean of 50 random initializations in bold dashed traces, where their standard deviations are displayed as ribbons. The experiment in Fig. 1b then subsamples 50% of the data entries uniformly at random, inserting zeros for the missing entries for methods that require fully sampled data, and compares the same stalistics across the algorithms.
compute the estimated subspace Ü, by
As expected, when the samples fully observed, PETRELS converges to the log-likelihood and subspace error for each set of training data the batch algorithms that assume homoscedastic noise, and SHASTA-PCA converges to the same log-likelihood value and subspace error as HePPCAT. Here see the advantage of using heteroscedastic data analysis. Instead of discarding samples from either data group or combining them in a single PPCA, the hetcroscedaste PPCA algorithms leverage both the "clean" samples, the additional "noisy" samples, and the noise variance estimates to produce better subspace estimates. With many missing entries (imputed with zeros), the batch algorithms' subspace estimates quickly deteriorate, as seen the right-hand side of Fig. 1b. Of the streaming PCA algorithms for missing data, SHASTA-PCA again attains the best subspace estimate compared to GROUSE and PETRELS.
Log like ihood: 100% observed
Subspace Error: 100% observed
1d
PETRELS
PETRELS-G1
2.50x10'
PETRELS-G2
GROUSE
SHASTA-PCA
-5.00x10'
HePPCAT
PPCA
10a
PPCA
-7.50x10*
-1.00x10"
Iteration
Iteration
(a) Fully observed data.
Log likelihoad: 50% observed
Subspace Error: 50% observed
00.0
PETRELS
PETRELS-G1
-2.50x10*
PETRELS-G2
GROUSE
10 0.5
SHASTA-PCA HEPPCAT
-5.00x10"
-7.50x10"
10.
PPCA-2
-1.00x10'
500 1000 1500 2000 2500
0 500 1000 1500 2000 2500
(b) Data with 50% entries observed uniformly at random.
Figure Incremental computation (with one pass) over batch data generated from a static subspace for 100, ni = 500, nı 2,000, ii = 10-2, and 02 = 10-1, Horizontal dashed lines show the terminal values for the batch algorithms, and the horizontal axis shows the iteration index for the online algorithms.
7.2 Dynamic subspace
This section studies how well SHASTA-PCA can track timo-varying subspace. We generate 20,000 stream- ing data samples according to the model in (1) for 2 groups with noise variances 01 10-4 and
U2
10-2
We use randomly drawn F U ddag((/)) where d 100 and 7 2, 1]. The data
13

---

Published in Transactions on Machine Learning Research (08/2025)
Subspace Error: 50% observed
10
PETRELS GROUSE SHASTA-PCA
10"
5,00x10²
1.00x10*
1.50x10*
2.00x10
Iteration
Figure 2: Dynamic tracking of rapidly shifting subspace with 50% of the entries observed uniformily at random using SHASTA-PCA versus streaming PCA algorithms that assume homoscedastic noise. Here, d 100 and 20% of the data has noise variance 10-4 and 80% of the data. has noise variance 10-2. Iterations refers to the number of streamed data vectors.
samples are drawn from the two groups with 20% and 80% probability, respectively. We then observe 50% of the entries selected uniformly at random. To simulate dynamic jumps of the model, we set the planted subspace U to a new random draw every 5000 samples and compare the subspace errors of the various methods with respect to the current U over time. Here, empirically selected the constant paramcters = cF = 0.01, and C, = 0.1 for SHASTA-PCA that do not decay with time to adaptively track the dynamics of the subspace. After hyperparameter tuning, we set the step size of GROUSE to be 0.02 and set the forgetting factor for PETRELS to be X = 0.998. Each algorithm is initialized with the same random factors Fu, and SHASTA-PCA's 's noise variances are initialized uniformly at random between 0 and 1. Fig. 2 shows SHASTA-PCA outperforms the reamingg PCA algorithms that assumee homoscedastic noise by half an order of magnitude. The results highlight how the largest noise variance dominates the streaming PCA algorithms" subspace tracking performance while SHASTA-PCA obtains more faithful estimates by accounting for the heterogeneity.
7.3 Dynamic noise variances
In some applications, due to temperature, age, change in calibration, the quality of the sensor measure- ments may also change with time (Jun-hua et al., 2003), thereby affecting the levels of noise in the data. To study the performance of SHASTA-PCA in these settings, we generate samples from the planted model described above where we change the noise variances over time while keeping the subspace stationary. As before, SHASTA-PCA initialized from a random (Fo, vg). Figs. 3a and 3c show the estimated noise vari- ances and the subspace error as we double the noise variance of the first group every 5,000 samples. Figs. 3b and 3d repeat the experiment but double the noise variance of the second group instead. As Vi increases and the cleaner group becomes noisier, the data becomes noisier overall and also closer to homoscedastic. SHASTA-PCA's 's estimate of the subspace degrades and approaches the estimates obtained by PETRELS and GROUSE. On the other hand, as the noisier group gels even noisier, the quality of the PETRELS subspace estimate deteriorates over time, whereas SHASTA-PA
remains robust to the added noise by leveraging the cleaner data group. In both instances, GROUSE appears to oscillate about an optimum in a region whose size depends on the two noise variances. The variance estimates demonstrate how SHASTA-PCA can quickly adapt changes in the noise variances; SHASTA-PCA adapted here less than 1,000 samples.
7.4 Comparison to streaming robust PCA
This section experimentally compares heteroscedasticc PCA and robust PCA methods. A variety of robust subspace tracking techniques exist in the literature, including for problems with missing data; see He et al. (2012); Bouwmans & Zahzah (2014); Zhan et al. (2016); Vaswani et al. (2018); Liu et al. (2020); Dung et al. (2021); Thanh et al. (2021). Many of these works model the data a true signal lying in low-dimensional
14

---

Published in Transactions on Machine Learning Research (08/2025)
SHASTA Variances: 50% observed
SHASTA Variances: 50% observed
10
V2
v1
10
10
10
10001100
1.00x10'
1.50x10
2.00x10
5. 5.00x10
1.00x10'
1.50x10
2.00x
lteration
(a) Estimated variances for varying #i.
(b) Estimated variances for varying 22
Subspace Error: 50% observed
Subspace Error: 50% observed
10
PETRELS GROUSE SHASTA-PCA
10a
PETRELS GROUSE SHASTA-PCA
10
10"
5,00x10'
1.00x10
1.50x10'
2.00x10
5. 3000x10'
1.00x10'
1.50x10
2.00x10
Iteration
Iteration
(c) Subspace crror with varying v1
(d) Subspace error with varying 02.
Figure Experiments with changing variances across time (iterations) for a single F R 100x3 starting with planted noise variances e [10-1, 10-2) with 50% of the entries observed uniformly random. The vertical dashed lines indicate points which we double one of the planted variances. The top plots show the estimated variances from SHASTA-PCA, and the bottom plots show the subspace error.
subspace, plus Gaussian noise, plus a sparse outlier vector. Given mild conditions on the incoherence and dimension of the subspace well as the sparsity and distribution of outliers' support, robust PCA provably decomposes low-rank plus sparse data (Chandrasckaran et al., 2011).
Here, we compare SHASTA-PCA with PETRELS-ADMM (Thanh et al., 2021), a streaming robust PCA algorithm². We consider subsets of data with increasingly larger differences in noise variances and varying ratios of "clean" " to "noisy" samples. one scenario, a minority of samples with a very large noise variance can be interpreted varying choices of and missing data rates. For each experiment, we vary the sparsity hyperparameter of PETRELS-ANE ADMM in increments of 0.1 between 0.1 and 2, choosing the hyperparameter that gave the lowest subspace error. We use a forgetting factor of À = 1 for both PETRELS and PETRELS-ADMM, and we use the same SHASTA-PCA hyperparameters as the experiments in Section 7.1.
outliers. We repeat the experiments in Section 7.1 for one random data draw with
Although PETRELS-ADMI models for sparse entry-wise outliers, our experiments Fig. show that it achieved lower subspace estimation error than PETRELS when the two noise variances differ significantly. This suggests that the low-rank-pp www-rank-luus-ppaee PCA model improves subspace estimation even under model mismatch and when entire groups of samples heavily contaminated with noise. Notably, Fig. shows that in low-SNR scenarios and consistently achieved lower subspace error when the data contained missing entries. However, SHASTA-PCA assumes and requires knowledge the data group memberships. When this assumption is met, streaming heteroscedastic PCA offers a robust solution for data with outlier groups.
SHASTA-PCA matched or outperformed
7.5 Computational timing experiments
Computational and/or storage considerations can inhibit the use of batch algorithms for large datasets, especially on resource constrained devices. To demonstrate the benefit SHASTA-PCA in such settings,
We used ChatGPT the authors' original MATLAB code (attps://github. cotps//ttt//tetttb//ssr to Julia, and validated its outputs.
15

---

Published in Transactions on Machine Learning Research (08/2025)
Subspace Error: 100% observed
Subspace Erron: 100% observed
Suospace Error: 100% observed
10
10
10
PETRELS
13
PETRELS-ADMI SHASTA-PCA
PPCA HoPPCAT
10 9.6
10-05
10
10-09
500
1000
1500
2000
2500
500
1000
1500
2000 2500
0 500 1000
1500
2000 2500
Iteration
Iteration
teration
(a)
12
[500,
2000].
(b) n 500, 2000],
[1, 10-1
(c) [000, 2000 110-..1
[10
10
Subspace Error: observed
Subspace
Subspace
10
10
1600
PETRELS
PETRELS-ADMM
16
SHASTA
105
HePPCAT
PPCA
10.0.
10-93
PPCA G1
PPCA-G2
10
167
0
500
1000
1500
2000
2500
0
500
1000
1500
2000
2500
0
500
1000
1500
2000
2500
(d)
[500,
2000],
[500, 2000],
[1, 10
( 5000, 2000]
[10-1.1]
[10
11-11
Figure 4: Experiments comparing SHASTA-PCA (this work), PETRELS (Chi et al., 2013), and PETRELS- ADMM (Thanh al., 2021) for synthetic data.
we generate a 2GB dataset according to our model, where 1,000, - 14, 24 n |50,000, 200,000). 1), and we observed only 20% of the entries uniformly random. For this experiment, we set W: 0.01/Vt, OF 0.01, and C2 0.1 for SHASTA's 's hyperparameters, and pass over the entire data once. Fig. 5 compares the convergence in ligeliiel values and subspace errors by elapsed wallclock time for SHASTA and the batch algorithm in $5, where both algorithms are randomly initialized from the same random starting itcrate (Fo, vo). SHASTA-PCA rapidly obtains good estimate of the model, using only roughly 60% of the time that it took the batch method, all while using only 0.0048% of the memory per iteration.
7.6 Real data from astronomy
We illustrate SHASTA-PCA on real astronomy data from the Sloan Digital Sky Survey (SDSS) Data Release 16 (Ahumada et al., 2020) using the associated DR16Q quasar catalog (Lyke et al., 2020). In particular, we consider the subset that was considered Hong et al. (2023, Section 8); see Hong et al. (2023, Supplementary Material SM5) for details about the subset selected and the preprocessing performed. The dataset contains n = 10,459 quasar spectra, where each spectrum is vector of 281 flux measurements across wavclengths and the data come with associated noise variances.
Ordering the samples from the smallest to the largest noise variance estimates, obtained a "ground-truth" signal subspace by taking the left & = 5 singular vectors of the data matrix for the first 2,000 samples with the smallest noise variance estimates. We then formed a training dataset with two groups: first, we collected samples starting from sample index 6,500 to the last index where the noise variance estimate less than or oqual to (7,347); second, we collected training data beginning estimate is greater than or equal 2 (8,839) up to the sample index 10,449, excluding the last 10 samples that grossly corrupted. The resulting training dataset. had n 848 and na 1,611 samples for the
the first index where the noise variance
16

---

Published in Transactions on Machine Learning Research (08/2025)
Log liilliiioodd:: 20% observed
Subspace Error: 20% observed
10
SHASTA PCA Batch method
-
5.00x10
0.5
10
-7.00x10
-1.00x10
100
200
300
100 150 200 250 300 350
Time (s)
Time (s)
(a)
Figure 5: Log- ligeliiel values and subspace errors versus elapsed wall clock time for one run of SHASTA-
PCA versus the batch method in $5 on 2GB of synthetic data: d = 1,000, k = 3, n = |50,000, 200,000), and == [0.1,1]. Both algorithms used the same random initialization. Markers for SHASTA-PCA are for every 10,000 vector samples, and markers for the batch method are for each algorithm iteration.
Truth
Variance
10a
SHASTA (L-n) Estimated variance
Subspace Error: 100%
10
PETRELS
500
1000
2500
SHASTA (L-n)
BBtch method (L-n)
HePPCAT (L=2)
250
HePPCAT (L-n)
Wavelength
200
1.0
PPCA Group
PPCA
150
PPCA PPCA-Group 2
100
0.5
50
500
1000 1500 2000 2500
500
1000
1500
2000
2500
teration
Quasar
(a)
(b)
Figure 6: (a) Visualization of quasar data with true associated variances of each sample plotted by quasar. The estimated variances from SHASTA-PCA L, n model streaming over the samples (fully sampled, in randomized order) closely matched the true variances. (b) For fully sampled data, the subspace error of SHASTA-PCA converged to the error of HePPCAT's estimated subspace for L = model. We repeated the experiment 50 times, each with different random initialization and order of the samples.
two groups, respectively, and had strong noise heterogeneity across the samples, where the second group was much noisier than the first. Fig. 6a shows the training datasct and the associated noise variance estimates for each sample.
Although we formed the data by combining two groups of consecutive samples, the data actually contain L - n groups since each spectrum has its own noise variance. This setting allows for heuristic computational simplification in SAAASAA-CAA Namely, we adapt the L = 1 model for a single variance u and use separate weights for the SMM updates of F and where w, 0.001, update is then equivalent to maximizing V, (F(Fi. v; Fi-1, 1,-1), i.e., the minorizer centered at the previous variance estimate with no memory of previous minorizers. The number of variance EM updates per data vector may be increased beyond just a single update, but in practice we observed little additional benefit.. Fig. 6 shows how SHASTA-PCA adaptively learns the unknown variances for each new sample and converges to the same level of subspace error the batch HePPCAT L model.
0.1, and w = = 1 Vt. The variance
17

---

Published in Transactions on Machine Learning Research (08/2025)
Subspace Error: 40% observed
200
PETRELS
SHASTA (L=n)
Batch method (L=n)
HePPCAT (L=2)
HePPCAT (L=n)
10
PPCA Group
PPCA Group
5.00x10' 1.00x10* 1.50x10' 2.00x10* 2.50x10 Iteration
Figure 7: Subspace error for quasar data with 40% of the entries observed uniformly at random. SHASTA- PCA streams over the full dataset 10 times, with the sample entries missing, where the order of the samples was randomized on each pass. Each experiment was initialized with factors chosen uniformly at random.
In many modern large datasets, entries may be missing in significant quantities due to sensor failure time and memory constraints that preclude acquiring complete measurements. Indeed, the experiment designer may only wish to measure a "sketch" of the full data to save time and resources and learn the underlying signal subspace from limited observations using an algorithm like SHASTA-PCA. To study this case, we randomly obscure 60% of the entries uniformly at random and perform 10 passes over the data, randomizing the order of the samples each time. We use the same choice weights described above to estimate a single variance for every new sample. We initialize with random F, and vo using the zero-padded data. As Fig. shows, SHASTA-PCA has better subspace estimates in this limited sampling setting than the saato-of-thh-al bascline methods with zero-filled missing entries and/or homoscedastic noise assumptions. Interostingly, the SHASTA-PCA subspace estimate is even better than the batch method in $5 for the L. = model.
8 Conclusion & Future Work
This paper proposes streaming PCA algorithm (SHASTA-PCA) that is robust to both missing data and heteroscedasstic noise across samples. SHASTA-PCA requires only a modest amount of memory that is independent of the number of samples and has efficient updates that can scale to large datasets. The results showed significant improvements over state-of-the-art streaming PCA algorithms in tracking nonstationary subspaces under heteroscedosticc noise and significant improvement over batch algorithm in speed.
There are many future directions building on this work. An interesting line of future work to establish convergence guarantees for SHASTA-PCA, particularly since our optimization approach unique among other works using stochastic MM. Second, while each update of inverting [2, -many & x & matrices per data vector, which is particularly wasteful when |2, = d since this is equivalent to the fully-sampled ligllike,,hood, which only requires inverting single k & matrix to update F. It may be possible to find other surrogate functions that would avoid this large number of small inverses in each iteration. Finally, although SHASTA-PCA enjoys lightweight computations in each iteration, achieving rapid convergence can depend on carcfully tuning the weights w; and the parameters op and c,- Adaptively selecting these parameters with stochastic MM techniques and developing theory to guide the selection of these parameters remain open problems.
in F is relatively cheap, it still requires
9 Acknowledgements
The authors acknowledge and thank Hanback Lyu for his helpful discussions on convergence analysis for stochastic MM algorithms.
18

---

Published in Transactions on Machine Learning Research (08/2025)
References
Romina Ahumada, Carlos Allende Prieto, Andrés Almeida, Friedrich Anders, Scott F Anderson, Brett H. Andrews, Borja Anguiano, Riccardo Arcodia, Eric Armengaud, Maric Aubert, et al. The 16th data release of the sloan digital sky surveys: first release from the APOGEE-2 southern survey and full release of eBOSS spectra. The Astrophysical Journal Supplement Series, (4) 2020.
Choirul Anam, Kusworo Adi, Heri Sutanto, Zaenal Arifin, Wahyu Budi, Toshioh Fujibuchi, and Geoff Dougherty. Noise reduction in CT images with selective mcan filter. J Biomed Phys Eng., 10:623- 634, 09 2020. doi:
Joakim Andén and Amit Singer. Factor analysis for spectral estimation. In 2017 International Conference Sampling Theory and Applications (SampTA), pp. 169 173. IEEE, 2017.
Laura Balzano. On the equivalence of Oja's algorithm and GROUSE. International Conference on Artificial Intelligence and Statislies, pp. 7014-7030. PMLR, 2022.
Laura Balzano, Robert D. Nowak, and Benjamin Recht. Online identification and tracking of subspaces from 2010 48th Annual Allerton Conference on Communication, Control, and Computing (Allerton), pp. 704-711, 2010.
highly incomplete information.
Laura Balzano, Yucjic Chi, and Yue M Lu. Streaming PCA and subspace tracking: The missing data case. Proceedings of the IEEE, 106(8):1293 1310, 2018.
Joshua K Behne and Galen Reeves. Fundamental limits for rank-one matrix estimation with groupwise heteroskedasticity. In International Conference on Artificial Intelligence and Statistics, pp. 8650-8672. PMLR, 2022.
Tamir Bendory, Alberto Bartesaghi, and Amit Singer. Single-particle cryo-electron microscopy: Mathemat- ical theory, computational challenges, and opportunities. JEEE Signal Processing Magazine, 37(2):58 76, 2020.
Dimitri P Bertsckas. Ineremental gradient, subgradient, and proximal methods for convex optimization: A survey. Oplimization for Machine Learning, 2010(1-38):3, 2011.
Cameron J Blocker, Haroon Raja, Jeffrey A Fessler, and Laura Balznoo. Dynamic subspace estimation with Grassmannian geodesics. arXiv preprint Xix:::333314855 2023.
Silvere Bonnabel. Stochastic gradient descent on Riemannian manifolds. IEEE Tronsactions on Automatic Control, (8))):2117-2229, 2013. ISSN 00189286. doi: 10.109//AAC2211325551199
Léon Bottou. Large-scale machine learning with stochastic gradient descent. In Yves Lechevallier and Gilbert Saporta (eds.), Proceedings of COMPSTAT'2010 pp. 177 186, Heidelberg, 2010. Physica- HD.
Christos Boutsidis, David P Woodruff, and Peilin Zhong, Optimal principal component analysis in distributed and streaming models. In Proccedings of the forty-eightt annual ACM symposium on Theory of Computing. pp. 236-249 2016.
Thierry Bouwmans and El Hadi Zahzah. Robust. PCA via principal component pursuit: A review for a comparative evaluation in video surveillance. Computer Vision and Image Understanding, 122:22 34, 2014.
Stephen P Boyd and Lieven Vandenberghe. Conver optimization. Cambridge University Press, 2004.
Arnaud Breloy, Guillaume Ginolhac, Alexandre Renaux, and Florent Bouchard. Intrinsic Cramér Rao bounds for scatter and shape matrices estimation in CES distributions. IEEE Signal Processing Letters, 26(2):262-25) 266, February 2019. doi: 10.1109/1sp2201..8887707
Emmanuel J. Candes and Yaniv Plan. Matrix completion with noise. Proceedings of the IEEE, (8(6))922599)) 936, 2010. doi: J0.1000..RR...2000.0333...
19

---

Published in Transactions on Machine Learning Research (08/2025)
Kaidi Cao, Yining Chen, Junwei Lu, Nikos Arechiga, Adrien Gaidon, and Tengyu Ma. Heteroskodastic and imbalanced deep learning with adaptive regularization. In International Conference on Learning Representations, 2021. URL https://openrevier. http://www....oooo.m/om//aaaaaaa
Venkat Chandrasckaran, Sujay Sanghavi, Pablo A. Parrilo, and Alan S. Willsky. Rank Dank-Sparsitty Incoherence for Matrix Decomposition. STAM Journal on Oplimization, (2):5572 596, April 2011. doi: 10.1137/
090761793.
Y. Chi, Y. C. Eldar, and R. Calderbank. PETRELS: Parallel subspace cstimation and tracking by recursive least squares from partial observations. IEEE Transactions on Signal Processing, (23))59944 5959, 2013. doi: TTw.:S..P...22222.
I. Cohen.
Multichannel post-filtering nonstationary noise environments. IEEF Transactions on Signal Processing, 52(5):1149 1160, 2004. doi: N0.1109/TSP:2004.55
Antoine Collas, Florent Bouchard, Arnaud Breloy, Guillaume Ginolhac, Chengfang Ren, and Jean-Philippe Ovarlez. Probabilistic PCA from heteroscedastie signals: Geometric framework and application to clus- tering. IEEE Trunsactions Signal Processing, 69:6546 6560, 2021. doi: 2TLL0//TSS..021.11100
Ilias Diakonikolas, Daniel Kane, Ankit Pensia, and Thanasis Pittas. Nearly-linear time and streaming algorithms for outlier-robust PCA. In International Conference Machine Learning, pp. 1886-7921.
PMLR, 2023.
Nguyen Viet Dung, Nguyen Linh Trung, Karim Abed-Meraim, et al. Robust subspace tracking algorithms in signal processing: A brief survey. REV Journal on Electronics and Communicattinns, 11(1-)), 2021.
A Hippert Ferrer, Mohamed Nabil El Korso, Arnaud Breloy, and Guillaume Ginolhac. Robust mean and covariance matrix estimation under heterogeneous mixed-effects model with missing values. Signal Pro- cessing, 188:108195, 2021.
Dimitrios Giannakis, Amelia Henriksen, Joel A Tropp, and Rachel Ward. Learning to forecast dynamical systems from streaming data. SIAM Journal on Applied Dynamical Systems, 22(():527755) 2023.
John Goes, Teng Zhang, Raman Arora, and Gilad Lerman. Robust stochastic principal component analysis. In Artificial Intelligence and Statistics, pp. 266 274. PMLR, 2014.
Andreas Grammenos, Rodrigo Mendoza Smith, Jon Crowcroft, and Cecilia Mascolo. Federated principal component analysis. Aduances in neural information processing systems, 33:6453 6464, 2020.
Jun He, Laura Balzano, and Arthur Szlam. Incremental gradient on the Grassmannian for online foreground and background separation in subsampled video. In 2012 IEEE Conference on Computer Vision and Pallern Recognilion, pp. 1568 1568-1575. 2012. doi:
Amelia Henriksen and Rachel Ward. AdaOja: Adaptive learning rates for streaming PCA. arXiv preprint 1XXi:19955.1115 2019.
A. Hippert-Ferrer, M.N. El Korso, A. Breloy, and G. Ginolhac. Robust low-rank covariance matrix estima- tion with general pattern of missing values. Signal Processing, 195:10 2022. ISSN 0165- doi: http://d....rg/....1//.//.g///// 108400 stticle/pp//0166664420000
URL httpp://wwww.sceenedd...e. com/science/
David Hong, Laura Balzano, and Jeffrey Fessler. Probabilistic PCA for heteroscedastic data. 8th IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing, pp. 26 30, 2019.
doi: 10.1109 10.1109/CAMSAPP /AMSSAA5577666201111 902294
David Hong, Yue Sheng, and Edgar Dobriban. Selecting the number of components in PCA via random signflips, 2020. URL attt:://axxer http::///ariw..o//aa//////ss /202.. 02985v1.
20

---

Published in Transactions on Machine Learning Research (08/2025)
David Hong, Kyle Gilman, Laura Balzano, and Jeffrey A. Fessler. HePPCAT: Probabilistic PCA for data with heteroscedastic noise. JEEE Transactions on Signal Processing, 69:4819-448 48844, 2021. doi: 10.1109/ TSP.2021.3104979.
David Hong. Fan Yang, Jeffrey A. Fessler, and Laura Balzano. Optimally weighted PCA for high-dimensional heteroscodastic data. SIAM Journal on Mathematics of Data Science, 5(1):222-250, March 2023. doi:
10.1137/22m14402444
M. W. Jacobson and J. A. Fossler. An expanded thcoretical treatment of iteration-dependent majorize- mininize algorithms. IEEE Trans. Im. Proc., 16(10):2411-2 22, October 2007. doi: 10.1109/TIP.2000 904387.
I. T. Jolliffe. Principal Component Anolysis. Springer-Verlee 2002. doi: 10.1007 /b98355
Liu Jun-hua, Shen Zhong-ru, et al. Drift reduction gas sensor by wavelet and principal component analysis. Sensors and Actuators B: Chemical, (6(1-2):354 363, 2003.
Zheng Tracy Ke, Yucong Ma, and Xihong Lin. Estimation the number of spiked eigenvalues in a covariance matrix by bulk eigenvalue matching analysis. Journal of the American Statistical Association, 118(541): 374 392, jul 2021. doi: 10.1080/01660 20t.://1221549.21.1111199
Syamantak Kumar and Purnamrita Sarkar. Streaming PCA for Markovian data. Advances in Neural Information Processing Systems, 36, 2024.
Boris Landa and Yuval Kluger. The Dyson equalizer: Adaptive noise stabilization for low-rank signal detection and recovery. Information and Inference: A Journal of the IMA, 14(1):iaac036, 2025.
Boris Landa, Thomas T. C. K. Zhang, and Yuval Kluger. Biwhitening reveals the rank of a count matrix.
SIAM Journal on Mathematics of Data Science, 4(4):1420-1446, dec 2022. doi: 10.1137 (0.1137/1mm14568)) 1211155607
Kenneth Lange. MM optimization algorithms. SIAM, 2016.
William Leeb and Flad Romanov. Optimal spectral shrinkage and PCA with heteroscedastic noise. Transactions on Information Theory, 67(5):3009 3037, 2021.
IEFE
William E Leeb. Matrix denoising for weighted loss functions and heterogeneous signals. SIAM Journal on Mathematics of Data Science, (3):9887 1012, 2021.
Ying Liu, Konstantinos Tountas, Dimitris A Pados, Stella N Batalama, and Michael J Medley. Ll-subspace tracking for streaming data. Pattern Recognition, 97:106992, 2020.
Cheng Lu, Jiusun Zeng Yuxuan Dong, and Xiaobin Xu. Strcaming variational probabilistic principal com- ponent analysis for monitoring of nonstationary process. Journal of Process Control, 133:103134, 2024.
Brad W Lyke, Alexandra N Higley, JN McLane, Danielle P Schurhammer, Adam D Myers, Ashley Ross, Kyle Dawson, Solene Chabanier, Paul Martini, Helion Du Mas Des Bourboux, el al. The Sloan Digital Sky Survey Quasar Catalog: Sixteenth data release. The Astrophysical Journal Supplement Series, 250 (1):8, 2020.
Hanback Lyu. Stochastic regularized majorizattion-miiiiiza with weakly convex and multi-convex sur- rogates. Journal of Machine Learning Research, 25(306):1-83, 2024.
Hanback Lyu, Christopher Strohmeier, and Deanna Needell. Online nonnegative CP-dictionary learning for Markovian data. Journal of Machine Learning Research, 23(148):1 50, 2022.
Julien Mairal. Stochastic majoriaation-minnnnzaaioo algorithms for large-scale optimization. Advonces in
Neural Information Processing Systems, 26, 2013.
M. Mardani, G. Matcos, and G. B. Giannakis. Subspace learning and imputation for streaming big data matrices and tensors. IEEE Transactions on Signal Processing. (6(10)226632677) 2015. doi: 10.1109/
21

---

Published in Transactions on Machine Learning Research (08/2025)
Arthur Mensch, Julien Mairal, Bertrand Thirion, and Gaël Varoquaux. Stochastic subsampling for factorizing huge matrices. IEEF Transactions on Signal Processing, (6(1):113-12)) 2017.
Aryan Mokhtari and Alec Koppel. High-dimensional nonconvex stochastic optimization by doubly stochastic successive convex approximation. IEEE Transactions on Signal Processing, 68:6287 6302, 2020. doi: N0.1109/TSP:2000.30002
Kevin Ni, Nithya Ramanathan, Mohamed Nabil Haji Chehade, Laura. Balzno,, Sheela Nair, Sadaí Zahedi, Eddie Kohler, Greg Pottie, Mark Hansen, and Mani Srivastava. Sensor network data fault types. ACM Trans. Sen. Netw., 5(3), jun 2009. ISSN 1550-4859. doi: 00.1145/122556.1525566 URL https://doi. org/10 1/5//0011140 11111111111111 1525863.
E. Oja. Simplified neuron model as a principal component analyzer. Journal of Mathematicaal Biology, 15:
267 273, 1982.
Dianne P O'Leary. Near-optimal parameters for Tikhonov and other regularization methods. SIAM Journal on scientific computing, 23(4):116111171, 2001.
Duy Nhat Phan, Sedi Bartz, Nilabja Guha, and Hung M Phan. Stochastic varianco-reduced majorization- minimization algorithms. SIAM Journal on Malhemalics of Dala Science, 6(4):926 952, 2024.
Klaas Pruessmann, Markus Weiger, Markus Scheidegger, and Peter Boesiger. SENSE: sensitivity
encoding for fast MRI. Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine, 42(5):952 962, 1999.
Haroon Raja and Waheed Bajwa. Distributed stochastic algorithms for high-rate streaming principal component analysis. Transactions on Machine Learning Research, 2022. ISSN 2835-8856. URL
Javier Salazar Cavazos, Jeffrey A. Fessler, and Laura Balzano. ALPCAH: Subspace learning for sample-wise wise heteroscedastic data. IEEE Transactions on Signal Processing, pp. 11, 2025. doi: 10.1109/TSP.2023 3537867.
Shuang Song, Kamalika Chaudhuri, and Anand D. Sarwate. Learning from data with heterogeneous noise using SGD. Journal of Machine Learning Research, 38:894-902, 2015. ISSN 15337928.
Le Trung Thanh, Nguyen Viet Dung, Nguyen Linh Trung, and Karim Abed-Meraie Robust subspace tracking with missing data and outlicrs: Novel algorithm with convergence guarantee. IEEE Transactions on Siynal Processing, 69:2070- 2085, 2021. doi: 10.1109/TSS TTP09///TP2221111 SSS.2001.3306799)
Michael E Tipping and Christopher M Bishop. Mixtures of probabilistic principal component analyzers.
Newral computation, 11(2):443-485 1999a.
Michael E. Tipping and Christopher M. Bishop. Probabilistic Principal Component Analysis. Journal of the Royal Statistical Society: Series B (Statistical Methodology), (1()):611 622, August 1999b. doi: 10711111146--986.8
Namrata Vaswani, Thierry Bouwmans, Sajid Javed, and Praneeth Narayanamurthy. Robust subspace learn- ing: Robust PCA, robust subspace tracking, and robust subspace recovery. IEEE signal processing mag- azine, 35 35(4):32256 2018.
Alec S. Xu, Laura Balzano, and Jeffrey A. Fessler. HeMPPCAT: Mixtures of probabilistic principal compo- nent analysers for data with heteroscedastic noise. In ICASSP 2023 2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 1-5, 2023. doi: T0.1109/TCAASS449337.. 10094719.
Yuling Yan, Yuxin Chen, and Jianqing Fan. Inference for heteroskodastic PCA with missing data. The Annals of Stalistics, 52(():729-775 -56,, 2024.
22

---

Published in Transactions on Machine Learning Research (08/2025)
Gale Young. Maximum likelihood estimation and factor analysis. Psychometrika, 6(1):49-53, February 1941. doi: 10. 10.1007/bM0228857
Jinchun Zhan, Brian Lois, Han Guo, and Namrala Vaswani. Online (and offline) robust PCA: Novel algorithms and performance guarantees. In Arthur Gretton and Christian C. Robert (eds.), Proceed- ings of the 19th International Conference on Artificial Intelligence and Statistics, volume 51 of Pro- ceedings of Machine Learning Research, pp. 1488 1496, Cadiz, Spain, 09 11 May 2016. PMLR. URL http://proceeeingss wttp:///g.......ees///5/s// Saat...ree/v/////haaa
html.
Anru R Zhang, T Tony Cai, and Yihong Wu. Heteroskedastic PCA: Algorithm, optimality, and applications.
The Annals of Stotistics, 50(1):53-8)) 80, 2022.
Yihan Zhang and Marco Mondelli. Matrix demoising with doubly heteroscedastic noise: Fundamental limits and optimal spectral methods. Advances in Neural Information Processing Systems, 77:93060-93117, 2024.
Heyang Zhao, Dongruo Zhou, Jiafan He, and Quanquan Gu. Optimal online generalized linear regression with stochastic noise and its application to heteroscedastic bandits. In International Conference Machine Learning, pp. 42259-42279 PMLR, 2023.
Yuchen Zhou and Yuxin Chen. Deflated HeteroPCA: Overcoming the curse of ill-conditioning in het- eroskedastic PCA. The Annals of Stalislics, 53(1):91 116, 2025. doi: 10.1214/44A0052450 URL https://doi.org org. 1211/11-1111111 200.///4-AA2256...
A Notation table
Table summarizes notation in the paper.
B
Minorizer derivation
We give a short proof for the expectation stop to derive the ith minorizer in (7). For notational simplicity, we drop the indexing on i and give the proof for the case with fully-sampled data, without loss of generality. Let
I
0
+
(26)
y
F
n
where I) and (V0, vf). Then ~ N
(..
I
vI vI+FF .FFFF)))
F
1)
Using conditional expectations
F
of multivariate Gaussian random variables,
ally=a~~Nii E),
(27)
where
== FF aa== I- FF') F (M
& M(F,
(28)
(29)
The two righthand equalities follow from applying the Sherman- Seaaan-Woobbrry-oooo
formula, ie., the matrix inversion lemma (Boyd & Vandenberghe, 2004, Section C.4.3). Finally, expanding the square quadratic terms of the complete likeliiidoood (4), taking the expectation with respect to zly, plugging the expressions for i and Σ above, and dropping terms constant with respect to F and yields the minorizer in (7).
23

---

Published in Transactions on Machine Learning Research (08/2025)
Table 1: Summary of notation.
Notation
Description
Yi Rd
ith data vector
F E R dd 2, E R G 1
factor matrix (k < d) ith coefficients distributed N ith noise vector distributed (0,, Ug I) noise variances of the L data groups vector of noise variances
e E RR_
9i
1, ...., }}
group index of the ith data vector set of entries observed for the ith sample observed indices over the batch data restriction of y, to the rows defined by Q restriction of F to the rows defined by Ω partially observed dataset.
Ω S
Ω
yn, E Rlai
Fa, E R xk
Yn
C(Yn; F, v)
log-likelihood over Yn with respect to F and v log-likelihood of the ith sample complete log-likelihood the ith sample function defined in (5) for intermediate calculations function defined in (6) for intermediate calculations minorizer of the ith log-likelihood at point (F,
Lyyn, F.
())
(F)
(F 0; F. 5)
( v; F.
minorizer over the batch data
tth approximate minorizer with respect to in the streaming algorithm tth approximate minorizer with respect to F in the streaming algorithm total count of observed entries in the fth group
(F)
Obatch
0, batch
term computed from the count of observed entries in the streaming algorithm data-dependent term computed in the batch algorithm's v update data-dependent term computed in the streaming algorithm's v update data-dependent & x k matrix computed in the batch algorithm's F update data-dependent h k matrix computed in the streaming algorithm's 's F update data-dependent k-length vector computed in the batch algorithm's F update data-dependent k-length vector computed in the streaming algorithm's F update sequence of weights in the streaming algorithm positive scalars the streaming algorithm for updating and v.
Pt.
Pt
R
Rij 8 batch Lg
SrJ
fweheo E 1)
cF and Co
C Hyperparameter sensitivity tests
Figures 8 and 9 summarize the hyperparameter sweep experiments. For the experiments in Section 7.1, we varied the constant cu in the numerator for decaying weights w, cw/t, where t is the algorithm iteration. We plotted the log- log-likelshoods and subspace errors in Fig. 8, showing the algorithm's performance is largely stable across many choices of hyperparameters. Similarly, for the dynamic subspace tracking experiments in Section 7.2, we plotted the subspace error of the final estimate for across varying values of a constant weight w = we for all and the constants Cp = Cy-
24

---

Published in Transactions on Machine Learning Research (08/2025)
log-likelihood error: 100% observed
Log-likelihood error: 50% observed
1.0
1.0
0.9
10000
0.9
10000
20000
20000
0.7
-30000
0.7
30000
0.5
40000
0.5
0.6
40000
0.4
50000
0.4
50000
0.3
60000
0.3
60000
0.2
70000
0.2
70000
0.01
80000
0.01
80000
0.001
90000
0.001
900000
(a)
(b)
Logarithmic subspace erro: 100% observed
Logarithmic subspace 50% observed
1.0
1.0
0
0.5
0.9
0.8
0.8
0.5
0.7
0.7
0.6
0.6
110
0.5
1.5
0.5
0.4
0.4
1.5
0.3
2.0
8"
0.3
2.0
0.1
2.5
0.1
0.001
0.01
-300
0.001
0.01
2.5
0000
0,00
Cp
(c)
(d)
Figure Hyperparameter sweep for experiments in Section 7.1.
Logarithmic subspace error
0.1
0.09
-.0
0.08
0.07
2.5
0.06
0.05
3.0
0.04
0.03
0.02
3.5
0.01
0.001
4.0
001
0.00
02
0.03
0.04
05
D6
200
188
100%
0
0.0
Figure 9: Hyperparameter sweep for experiments in Section 7.2.
25


Reviews:
[{"review_id": "mbZZbw7JZn", "reviewer": "Reviewer_QkZF", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The paper tackles an interesting problem but unfortunately I find the writing of the paper to be quite dense. Below, I point towards specific instances:\n\nIn Section 3, it will be better if is clearly mentioned that k is usually much smaller than d. The goal of modeling via F is to ensure that all the observed vectors lie in a low dimensional subspace.\n\nWhy do the authors assume that group memberships are known? Why is this a valid assumption in the first place?\n\nWhere is the proof of equation 3? It seems non-trivial to me at first glance\n\nI was completely lost in Section 4 - the authors need to write sentences like the following in a significantly better fashion\n\nNext, we take the expectation of (4) with respect to the following conditionally independent distributions obtained from Bayes’ rule and the matrix inversion lemma:\n\nWhat is the matrix inversion lemma you are referring to? Cite or put the lemma in the appendix.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "In this paper, the authors look at the problem of modeling high dimensional data that lies in a lower dimensional subspace using gaussian noise with different noise levels. In particular, the authors look at the setting of streaming data with missing entries and heteroscedastic noise - in this setting, the authors consider varying levels of noise variance depending on the sample. In such an instance, the authors introduces an algorithm SHASTA-PCA that has low memory and computational overhead.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "See above", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "X1WjuBLZ6a", "reviewer": "Reviewer_mM8x", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "### Strengths\n- Demonstrates clear improvements in subspace recovery accuracy and runtime compared to GROUSE, PETRELS, and batch heteroscedastic PCA (HePPCAT) on both synthetic and SDSS DR16 datasets.  \n- Successfully recovers the underlying low-rank subspace even with up to 50 % randomly missing entries per sample.  \n- Requires only $O(dk^2)$ memory and $O(|\\Omega_t|k^3)$ compute per sample, independent of the total stream length.  \n- Integrates heteroscedastic noise modeling and online updates in a single end-to-end algorithm, generalizing PETRELS (streaming + missing data) and HePPCAT (batch + heteroscedastic noise).\n\n### Weaknesses\n-  Currently, the paper does not include a formal proof or sketch of conditions under which the stochastic majorize–minimize updates converge to a stationary point, even under a warm start as is typically the case for EM algorithms. Even a concise convergence guarantee or proven stability conditions would enhance confidence in the algorithm’s reliability.\n\n- The submission does not compare against in detail against state-of-the-art adversarial/outlier-robust PCA methods such as:  \n  - Ilias Diakonikolas *et al.*, “Nearly-Linear Time and Streaming Algorithms for Outlier-Robust PCA” (ICML 2023): develops a nearly-linear time algorithm with near-optimal error guarantees for robust PCA and a single-pass streaming variant using nearly-linear space [\\[1\\]](https://arxiv.org/abs/2305.02544).  \n  - Eric Price and Zhiyang Xun, “Spectral Guarantees for Adversarial Streaming PCA” (FOCS 2024): provides the first adversarial streaming PCA algorithm with spectral error bounds [\\[2\\]](https://arxiv.org/abs/2408.10332).  \n  Including either an empirical comparison or a discussion of how SHASTA-PCA complements or outperforms these approaches would clarify its place in the robust PCA landscape.\n\n- The roles of the step-size weights $w_t$ and regularization constants $(c_F, c_v)$ are explained, but practical tuning advice or an ablation study is not provided. A brief sensitivity analysis would help practitioners reproduce and adapt the method to different datasets.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "In this submission, the authors introduce **SHASTA-PCA**, a novel streaming algorithm for probabilistic principal component analysis that simultaneously handles heteroscedastic noise, missing data, and non-stationary observations. At each time step $t$, one observes a partially‐observed vector $x_t\\in\\mathbb{R}^d$ with index set $\\Omega_t\\subseteq[d]$ and noise $x_t = F v_t + \\varepsilon_t,\\quad \\varepsilon_t\\sim\\mathcal{N}\\bigl(0,\\mathrm{diag}(\\sigma_t^2)\\bigr),$ where $F\\in\\mathbb{R}^{d\\times k}$ is the low‐rank factor and $\\sigma_t^2\\in\\mathbb{R}^d$ encodes sample‐wise variances. The core technical innovation is a stochastic majorize–minimize (SMM) scheme that yields closed‐form, per‐sample updates for both $U$ and $\\sigma_t^2$, avoiding expensive batch EM steps. These updates incur only $O(dk^2)$ memory and $O(|\\Omega_t|k^3)$ computation per iteration. \n\nThe paper provides a a theoretical complexity analysis of per-step updates, discusses empirical convergence under mild conditions, and validates performance on synthetic data as well as real SDSS DR16 quasar spectra, demonstrating superior accuracy and substantial runtime gains over GROUSE, PETRELS, and batch heteroscedastic PCA (HePPCAT).", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "I believe the paper makes a good contribution towards combining the settings of heteroscedastic noise, missing data and streaming PCA. Please see the weaknesses section above for proposed adjustments. However, they are not critical to securing my recommendation for acceptance but would strengthen the paper in my opinion.", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "HS84OBMNqZ", "reviewer": "Reviewer_yPyN", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "### **Strengths**\n- The authors convincingly argue that real-world streaming data often contain both missing values and heteroscedastic noise, and no existing PCA method addresses both.\n- The proposed method builds on stochastic MM techniques and adapts them to heteroscedastic noise and missing data.\n- Experiments demonstrate good performance on both synthetic and real-world data.\n\n### **Weaknesses**\n- The paper lacks theoretical convergence results; but the authors discuss this and highlight the challenges.\n-  SHASTA-PCA depends on the tuning of hyperparameters. Coming with a more adaptive or self-tuning strategy could be beneficial, in particular, in the setting of streaming data.\n- The algorithm, while efficient, is relatively intricate and may be hard to reimplement without detailed code or pseudo-code beyond what is provided in the paper.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposes SHASTA-PCA, a streaming algorithm for probabilistic PCA that can handle both heteroscedastic noise across samples and missing data. Existing PCA methods either assume homoscedastic noise or require access to all data (full batch). Thus, a streaming variant of PCA for streaming data (with heteroscedastic noise and missing entries) is a natural extension.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "### **Major**\n- The paper should be heavily revised. It could be shortened (by adding some appendices), and remove repeating sentences. Also, the author use the same notation for different quantities which is confusing to read. E.g., the authors could clearly state that $t$ is time, $i$ is the iteration, and k,l,s is the index for their sums. The authors have eq. numbers on all eqs.! but rarely refer to them...\n- Add theoretical convergence analysis - similar results as in the literature - would improve the paper significantly.\n- Code release or pseudo-code: Include a more detailed algorithmic description or commit to releasing code for reproducibility.\n- Add notation table which summarize key symbols used throughout would improve readability.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)