Streaming PCA for Markovian Data
Authors: Syamantak Kumar
Purnamrita Sarkar
Venue: NeurIPS 2023 spotlight (2023)
Keywords: Streaming PCA
Markov Chain
Mixing
Oja's algorithm
TL;DR: We prove near-optimal rate of convergence of Oja's Streaming PCA algorithm for data sampled from a Markov chain and show that it is better than discarding data to reduce dependence at finding the leading eigenvector.

Abstract: Since its inception in 1982, Oja's algorithm has become an established method for streaming principle component analysis (PCA). We study the problem of streaming PCA, where the data-points are sampled from an irreducible, aperiodic, and reversible Markov chain starting in stationarity. Our goal is to estimate the top eigenvector of the unknown covariance matrix of the stationary distribution. This setting has implications in scenarios where data can solely be sampled from a Markov Chain Monte Carlo (MCMC) type algorithm, and the objective is to perform inference on parameters of the stationary distribution. Most convergence guarantees for Oja's algorithm in the literature assume that the data-points are sampled IID. For data streams with Markovian dependence, one typically downsamples the data to get a "nearly" independent data stream. In this paper, we obtain the first near-optimal rate for Oja's algorithm on the entire data, where we remove the logarithmic dependence on the sample size, $n$, resulting from throwing data away in downsampling strategies.

Streaming PCA for Markovian Data
Syamantak Kumar!
Purnamrita Sarkar2
Department of Computer Science, UT Austin Department of Statistics and Data Sciences, UT Austin
syamantak@utexas.c edu, purna sarkareasstiin utexas edu
Abstract
Since its inception in 1982, Oja's algorithm has become an established method for streaming principle component analysis (PCA). We study the problem of streaming PCA, where the data-points are sampled from an irreducible, aperiodic, and reversible Markov chain starting in stationarity. Our goal is to estimate the top eigenvector of the unknown covariance matrix the stationary distribution. This setting has implications in scenarios where data solely be sampled from a Markov Chain Monte Carlo (MCMC) type algorithm, and the objective is to perform inference on parameters of the stationary distribution. Most convergence guarantees for Oja's algorithm in the literature assume that the data-points are sampled IID. For data streams with Markovian dependence, one typically downsamples the data to get a "nearly" independent data stream. In this paper, we obtain the first oetr-oltimal rate for Oja's 's algorithm on the entire data, where we remove the logarithmic dependence on the sample size, resulting from throwing data away in downsampling strategies.
Introduction
Streaming Principal Component Analysis (PCA) is an important and well studied problem where the principal eigenvector of the sample covariance matrix of a dataset is computed one data-point at time, One of the most popular algorithms for streaming PCA was introduced by Erkki Oja in 1982 Most existing analyses of Oja's algorithm are done when the data is sampled IID.
29 30
However, in many practical applications, the data-points are dependent and are sampled from an MCMC process converging to a target stationary distribution. This naturally arises in the context of token algorithms for Federated PCA settings [10
with multiple machines communicating via fixed and connected graph topology. Each machine contains an arbitrary fraction of data-points and the goal is design streaming algorithm that respects this topology and returns the principal component of the whole dataset. This is typically achieved using Metropolis- Metropolis-Hasttings scheme that uses local information to design the transition matrix of Markov chain with any desired stationary distribution. The stationary distribution, of the random walk is chosen so that the distribution of the samples under π matches the uniform distribution over data-points. Governed by this Markov chain, random walker then travels the network of machines and samples one data-por point at time from the current machine and computes the update. However, even under the stationary distribution, the data- kata-points are dependent, which deviates from the IID setup. Our goal is to obtain near-optimal analysis of the sin' error of the estimated vector with respect to the true top eigenvector of the unknown covariance matrix in the Markovian setting.
12 13
Estimating the first principal component with streaming PCA: Let X, be mean zero d dimensional vector with covariance matrix E, and 7) be a decaying learning rate. The update rule of Oja's 's algorithm is given as
Wr
W: f (I+ (1++)XXX?+) w wi |
(I)
where w, is the estimate of U1, the top eigenvector of Σ, and ne is the step-size at timestep t. We aim analyse the sin? error of Oja's iterate timestep 1, defined
(w,, v2
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

log-free main-term
Paper
Markovian?
Online?
sin error rate
Sample complexity
Y
Y
V gap Ing(d) gap"
Jain et al.
N
gap? log(d)
16
N
N
gap
Chen al.
Y
Y
N
O
G gap² e log?
@gapp r
3
Vlog d
Vlog d==
Neeman al.
Y
N
N
O
O
[I (1-122 gap
(1-1A2(P)) (P)))ppp>
28
Theorem
Y
Y
[1 (P2(P)])
(I (2( (P) p
Table I: Comparison of sin error rates and sample complexities. Here gap := (A1 - A2), where \1, are the top 2 eigenvalues of D and V represents suitably defined variance parameter (see assumption
2). The sample complexity represents the number of samples required to achieve sin? error most We note that []] and [I5] also match the online sample complexity bound provided in [16]. Further, for the offline algorithm with IID data, [17] removes the log (d) factor in exchange for constant probability of success for large enough
Streaming PCA in the IID setting: For an IID data stream with E [X,] = 0 and E [ Σ. there has been a lot of work on determining the non-asymptotic convergence rates for Oja's algorithm
and its various adaptations [16] 1 3 38 14, 15 26 21 25 Amongst these, [16], 1 and [15] match
the optimal offline sample complexity bound suggested by the independent and identically distributed (IID) version of Theorem I (See Theorem 1.1 in [16).
We consider Oja's algorithm in the setting where the data is generated from reversible, irreducible, and aperiodic Markov chain with stationary distribution æ. We denote by [-[-] the expectation under the stationary distribution. In this setting our goal
estimate the principal eigenvector of E,, [ As in the IID setting, E,, [X,] = O. The challenge is that the data, even when it reaches stationarity, is dependent. Here the degree of dependence is captured by the second eigenvalue in the magnitude of the transition matrix P' (denoted as (P)( of the Markov chain. This is closely related to the mixing time of Markov chain [9]] denoted as Tmix which is the time after which the conditional distribution of state is close in total variational distance to its stationary distribution, (See Section 2.I).
Using a series of approximations, we obtain an optimal error rate for the sin? error, which worse by factor of 1/(1 - [A2(P)]) from the corresponding error rate of the IID case. Previous work 3 has established rates worse by a poly-logarithmic factor by using downsampling, i.c. applying the update on every 4 datapoint. In Figure [ we compare Oja's 's algorithm with its downsampled and offline variants (see Section for more details on setup). We see that Oja's algorithm performs significantly better than the downsampled variant, and similarly to the offline variant where for the 1 data point we compute the eigenvector of the sample covariance matrix of all data-points up-to Our work provides a concrete and novel result that explains these observations. In Table we compare our bounds with related analyses of Oja's 's algorithm. The last row shows that we are the first to obtain an error whose main term is free of logarithmic dependence on n or d for streaming PCA in the Markovian case.
Our contribution:
We break the logarithmic barrier in previous work by considering a series of approximations of finer granularity which uses reversibility of the Markov chain and standard mixing conditions of irreducible and aperiodic Markov chains. Our rates are comparable to the recent work of (28) (Proposition I that establishes an offline eтror analysis for estimating the principal component of the empirical covariance matrix of Markovian data by using Matrix Bernstein inequality. Therefore, our results are nearly optimal in terms of the dependence on sample size, the dimension d, and the model parameters V. (A1 >x2) and D (P) (See Section 2 for definitions). Our results also imply a linearly convergent decentralized algorithm for streaming PCA in distributed setting. As simple byproduct of our theoretical result, we also obtain a rate for Oja's algorithm applied on downsampled data, which is worse by a factor of log n. as shown in Figure To our knowledge, this is the first work that analyzes the Markovian streaming PCA problem without any downsampling that matches the error of the offline algorithm.

---

The crux of our analysis uses the mixing properties of the Markov chain. Strong mixing intuitively says that the conditional distribution of state in timestep given the starting state exponentially close to the stationary distribution of 8, the closeness being measured using the total variation distance. All previous work on Markovian data exploits this property by conditioning on states many time steps before. However, it is crucial to a) adaptively find how far to look back and b) bound the error of the sequence matrices we ignore between the current state and the state we are conditioning on. Observe that these two components are related. Looking back too fär makes the dependence very small but increases the error resulting from approximating a larger matrix product of intermediate matrices. We present fine analysis that balances these two parts and then uses spectral theory to bound the second part within a factor of a variance parameter that characterizes the variability of the matrices and shows up in the analysis of [6 28
Related work on streaming PCA and online matrix decomposition on Markovian data: Amongst recent work, 3 is very relevant to our setting, since it analyzes Oja's algorithm with Markovian Data samples. Inspired by the ideas of [ the authors propose downsampled version of Oja's algorithm to reduce dependence amongst samples and provide Stochastic Differential Equation (SDE) based analysis to achieve a sample complexity of O (A1->23² log²
(oi- for sin? error smaller than e, where G a variance parameter. We obtain similar rate Corollary I [trouuhh our techniques. However, comparing with Theorem I we observe that downsampling leads to an extra O (log (n)) factor. is important to point out that [3] provides an analysis for estimating top k principal components, whereas this paper focuses on obtaining log-free error rate for the first principal component.
G
[22] consider the harder problem of online non-negative matrix factorization for Markovian data. Their analysis establishes asymptotic convergence of error, but does not provide rate.
Downsampled
Oflline
Oja's 's algorithm
Stochastic Optimization with Markovian Data : Markovian models are often considered in Reinforcement Learning and Linear Dynamic Systems 2 5
0.6
9 31 35
18
24]. There
have been many notable nonasymptotic bounds for stochastic gradient descent (SGD) methods for general convex and nonconvex functions with Markovian data 88 32 6. 7 TO 39] 34]
0.2
The convergence rates (sample complexities) obtained in these works apply to more general problems but do exploit the matrix product structure inherent to Oja's algorithm. In this work, we develop novel techniques to show that a sharper analysis is possible for the PCA objective. The paper is organized as follows. Section 2 contains the problem setup and preliminaries about Markov Chains. Section 3 contains Theorem I We present sketch of the main technical tools in Section 4 intermediate theorems needed for the main theorem in Section and conclude with simulations in Section 6
0.0
1000
2000 3000
4000
5000
Comparison of Oja's algorithm with downsampled and offline variants. The X-axis represents the sample size and the Y axis represents the sin't error of each algorithm's 's estimate the leading eigenvector. The experimental setup is available in Section
Figure 1:
2 Problem Setup and Preliminaries
This section presents the problem setup and outlines important properties of the Markov chain that will be utilized subsequently. We assume that:
Assumption The Markov chain irreducible, aperiodic, reversible, and starts in stationarity,
with state distribution !
In practice, one can let the Markov chain mix until a burn-in period of O (7mix steps to achieve state distributions close to the stationary one. A simulation experiment with different burn-in in periods can be found in the Supplement Section S.6
3

---

Such a Markov chain can arise in various situations, for e.g., while performing random walks on expander graphs which are used extensively in fields such computer networks, error-correcting codes, and pseudorandom generators. Each state 8 of the Markov chain is associated with distribution D(s) over d-dimensional vectors with mean μts IR and covariance matrix Eg Rdxdd
For a random walk 51, 82, S5 St on the Markov chain, we define the sequence of random variables X2, ....... Xt, where conditioned the state 8;, X; ~ (s;)). We represent the mean as fl := Bown [ Es TT and the covariance matrix as E E Raxd, which, for E [ra] can be expressed as:
E: ES (X: - (x; - 11)²] Egwn [S)]] Big [1s
In this work, we assume μ = 0, which is a common assumption in the IID setting (see [T) D. it may be possible to extend our analysis the non-zero mean case, is out of the scope of this paper. Therefore. Σ = En Ep(s) [X,X;]] for i [n]
While
Let the eigenvalues of Σ be denoted λ A2 2 .. Aa- Let denote the leading eigenvector of Σ and VL denote the [ddx(d-1) matrix with the remaining eigenvectors columns. In this work, Il-llz denotes the Euclidean L2 norm for vectors and the operator norm for matrices unless otherwise specified and I denotes the identity matrix of the appropriate dimensionality. We proceed under the following standard assumptions for E [n]. (see for eg. (5)..
Assumption 2.
E) S |X,x7 - E|12 < M with probability 1.
Assumption 3.
Assumption 3 also implies |[ + Mapu"" - Σ1le I M with probability 1. WLOG, we assume M + À 2 I. We use E == Eg Es(s)) [-] to denote the expectation over state ~ т and over the state-specific distributions D (-), unless otherwise specified. We observe that Assumption 3 can be relaxed to accommodate sub-gaussian data distributions with appropriate variance decay (see Proposition in (2I)) using standard truncation-based arguments. We do not pursue this here. Define the matrix product
B (T+nx))??)) (T xx-IIx( ... ... xix))
Unrolling the recursion in Eq I the output of Oja's algorithm at timestep t is given as w, [o
B,wn
2.1 Markoy chain mixing times
Now we will discuss some well-known properties of an irreducible, aperiodic, and reversible Markov chain (see [19 for details and derivations). Let \>2 (P) denote the second largest absolute eigenvalue of the Markov chain; let the state-distri of the Markov chain at timestep with 81 = x be P((x. .). For any two probability distributions and recall that the total variational distance is TV (2) = |||| - V2 lTv Exealvi (x) - (z)). The distance from т at the tth timestep defined as dmix (1) == supzen TV (xx .). T). For irreducible and aperiodic Markov chains, by
Theorem 4.9 in [19], we have (l) < C' C'ep(-ee)) for some C, >> 0. The mixing time of the
Markov chain is defined as:
Tinix() : inf duni (t) (c)
As defined in Section 4.5 of [19] will denote Tmix Tmix (t). Then, have Tiix (c) [1gc(1 Tiix. It is worth mentioning the useful relationship between dmix (.) and E No. These results about mixing time are valid for general irreducible and aperiodic Markov chains. A reversible Markov chain satisfies V v, E Ω. (z) P (z. y) x (y) P (). For a reversible, irreducible, and aperiodic Markov chain, the gap (2 P) is inversely proportional to ( (Theorem 12.4, Tmix [9).
Tmix- given as dmix (ITmix) < 2-
3 Main Results
In this section, we present our main result, near-optimal convergence rate for Oja's algorithm on Markovian data. As corollary, we also establish rate of convergence for Oja's 's algorithm applied

---

on downsampled data, where every hthh data-point is considered. Supplement Section S.5 contains
comprehensive proofs Theorem I and Corollary [ while the proof of Proposition T can found in
Supplement Section S.2
Theorem
Fix a ò E (0, 1) and let the step-sizes be qj = (A with 110
S a. >> 2. Under
assumptions
/ and for sufficiently large number of samples v. such that
n
log
log
1000a max Tiiix log (1) 700
(M+
(M + +11)²
11g177 +4i
100
B
(4)
(A1 x2) log (1 + 200
the output wn of Oja's is algorithm satisfies
1- (-(w 01) 2 82
log (3)
d
2B
2x
CIV
+ C2M +112
Timix
1m
(A1- 1-- A2 (P))
(A1 -22)²
n1a
with probability atleast (1 - 8). Here C is an absolute constant and
C = a² (3 + 7/22 (P)D.
C2 : 35a
2a
Remark 1. (Interpreting the sin error) Theorem simply establishes that, under suitable choices of a. B, as long assumptions 7 andd 3 hold, for decaying step sizes ni (21-)2)(8+d))
with
probability at least 3 4. we have,
V ( 122 (P)
1- 1=
1
(>1 ( -2)² (1
for large enough This be seen by using 6. setting 3 as in Eq E and taking Ω (B)1.)) where Si hides logarithmic in These choices make the first and third terms of the upper bound the sin' error in Theorem smaller order than the second term, which is
( For the first term, we have d
28
20
=0 (
The third term involves
Tmie min
which is
n
O log" (ig((m) n
using the definition of Tmix (e) in Eq 3
Remark (Knowledge of mixing times) note that setting the step-sizes Y6s Theorem ① requires the knowledge of the mixing time of the Markov chain, Timx- An offline algorithm for estimating the mixing time is provided in [37]. They show that at least a constant Emix samples are needed for estimating Tiis within an absolute error e. Therefore, it is possible to use the first O (+) samples to estimate the mixing time and then set the learning rate accordingly.
Remark 3. (Regarding non-reversible Markov chains) We assume the reversibility of Markov chains simplify our analyses. Our techniques and results extend the non-reversible case using reversibilization tools (see (20), (II]). Specifically, one can obtain bounds in terms of the second-largest eigenvalue,
magnitude, of the modified transition matrix P.P. the multiplicative reversibilization of P. Here, P. represents the probability transition matrix of the time reversed chain, given Vz, Ω P. (x,v)
z(v) P (y. xz). We leave this extension for future work.
x(z)
Next, we compare the rate of convergence proposed in Theorem T with the offline algorithm having access to the entire dataset {X,}* 1 using recent result from 128]. Here, the authors extend the
Matrix Bernstein inequality [6 33 to Markovian random matrices. Their setup is much like ours except that the matrix at any state is fixed, i.e., there is no data distribution D(s) as in our setup. However, it easy to extend their result to our setting by observing that conditioned on the state sequence, the matrices X, X;XF, E [r]] are independent under our model, and we can push in the expectation over the state-specific distributions, (s)), whenever required. Therefore, we have the following result-
5

---

Proposition 1 (Theorem 2.2 of [28 theorem). Fix 8 E (0. 1). Consider an irreducible and
aperiodic Markov chain. Under assumptions 3 with probability 1 8. the leading eigenvector è of Σ To
== X;X? /n satisfies,
Vlog L- 7
Mlog d +
A2 (PI)) 1 X2 (P)
1 =
1 - vi" S sc (A1 A2)
C2
n
(A1 22) (1 12 (P)))
'n2
for absolute constants Cj and C2
Observe that Theorem matches the leading term
(A1 (11-Aa)P11-Paa (1 in Theorem except the log (d)
v
term. We believe, much like the IID case (also see footnote in [16]. this logarithmic term in 28 result is removable for large n and a constant probability of success.
Remark 4. (Comparison with IID algorithm) Fix å E (0, the data-points {Xi)-= are sampled IID from the stationary distribution πT, then under assumptions 2and using Theorem 4.1
from 166, we have that the output Wp, of Oja' 's algorithm! 7 satisfies.
Coo (2) 82
8
2a
a 2²
1
1 (=
d
(2a' - 1) (AI x2)
The leading term of Theorem worse by a factor of 1- (P)) Further, it has an additive lower
order term O (nn n²
due to the covariance between data-points in the Markovian case.
Corollary 1. (Downsampled Oja's 's algorithm) Fix a 6 E (0. If Oja's algorithm is applied on the downsampled data-stream with every 1shh data-point, where == Tmix (1/) then under the conditions of Theorem 7 with appropriately modified a and 3, the output wn satisfies,
1 (w 01)
C'log (+) 82
(2
2mB log (n)
C1
log (n)
(MM(++AI log? (n) Tmix (7hm) 2
d
+
n
(A1 A2)
n
(A1 -22)
112
with probability atleast (1- Here C is an absolute constant and C1 := 300
Ia-I C2
35a
Remark 5. Data downsampling to reduce dependence amongst samples has been suggested in recent work 127 23 3y. In Corollary [ establish that the bound the rate obtained is sub sub-optimal compared to Theorem]] [ by a log (n) factor. We prove this by a simple yet elegant observation: the downsampled data stream can be considered to be drawn from Markov chain with transition kernel pk (..-) since each data-point is l steps away from the previous one. For sufficiently large ki this implies that the mixing time of this chain is θ (1). These new parameters are used select the modified values of a, according to Lemma S. in the Supplement.
The proof of Theorem follows the same general recipe as in [16] for obtaining a bound on the sin? error. However, the original proof techniques heavily rely on the IID setting. We carry out a refined analysis for each step under the Markovian data model by careful control of error terms arising out of dependence. The first step involves obtaining high-probability bound on the sin' 2 error, by noting that Oja's 's algorithm on n data-points can be viewed as single iteration of the power method on B, Therefore, fixing 8 (0, 1) using Lemma 3.1 from [16). we have with probability at least (1 8).
C Clog (+) Tr (V? BY, B of Bi BT
sin? (Wn, w1)
5
where C' is an absolute constant. The numerator is bounded by first bounding its expectation (see Theorem 3 and then using Markov's inequality. To bound the denominator, similar to [16]. we will Chebyshev's inequality. Theorem provides lower bound for the expectation E. [u] B,, B vn]. Chebyshev's inequality also requires upper- bounding the variance [o] B, BT v1]. which requires (uf B, v1)
us to bound
(see Theorem

---

4 Main Technical Tools
In this section, we provide a sketch of the main arguments used in our proof.
Warm-up with downsampled Oja's algorithm: We start with the simple downsampled Oja's algorithm to build intuition. Here, one applies Oja's update rule (Eq to every k*h data-point, for a suitably chosen
For &- Tmix log n. the total variation distance between any consecutive data-points in the downsampled data stream is (nn-),,, As we show in Corollary | the error of this algorithm is similar to the error of Oja's s algorithm applied to n/k data-points in the IID setting, i.e., m logni n/n).
We will take E [u] B. B. U1| as an example. Let us introduce some notation.
By. (T My Xxxz) (1 + n
(I x:X;)?))
(7)
We peel this quantity one matrix a time from the inside. Note that for a reversible Markov chain, standard results imply (see Lemma D that the mixing conditions apply to the conditional distribution of state given another state steps the "future" (see Supplement Section S.3 for a proof). Recall dmix (k) from Section 2.1
1 sup Σ (Z, s Zit t) π (s) dmni (k). 2 TEn
Lemma 1. Under Assumption
It will be helpful explain our analysis by comparing it with the IID setting. For this reason, will use [m][-] to denote the expectation under the IID data model. Define On,, := |u] Bn,; BT, ; [n)]. Then we have.
0n,1 Bm. (I +
Σ)) (I + m m (Xxxx Σ)) B BT,201
ui Bn (I + E)² 2i + TT + T2,
(8)
where the first term is smaller than (1 + ))>1) Qn,2- We define Ti and T2 as follows. [u] B, (I + m (X,XT' -)) 21 and T2 E of Ba, 2 (X1x7 2)² BT
For the IID setting, the second term is zero, and the third term can be bounded as follows:
EuD of Bn,2 (x,x7 Σ)² BTT 2T
EaD
uf Br. E (X1x7->7 =))" 2T
I VEED [u] B, Br 21
Let us denote the IID version of an,, by a == End Lou BT, v1). The final recursion for the IID .nD
.ID
case becomes: 22 (1 An 'n (Af + aa So, for our Markovian data model, the hope that the cross term Tj (which has multiplicative factor of m) is (()) and T2 O(n)). We will start with the 7i term, which is zero in the IID setting.
E
I+n
x,x,'x
Figure 2: If we could replace the intermediate products (white matrices) by 1, the conditional expectation of the noise matrix X Σ conditioned on the gray matrices would be nearly zero.
We will show that Ti. while not zero like the IID case, is still sufficiently small. Intuitively, if we could replace the matrices between X, x - Σ and B,,2+ for some suitably large integer by identity (see Figure 2), then using (reverse) mixing properties of the Markov chain, could argue using Lemmal that E [x xf - _ 8,] is very close to zero (see Figure The following lemma formally bounds the deviation of the length-k matrix product from identity.
Lemma 2. Let Assumption 3 hold. If Vi [m], mike (M + x1) S E,E (0, and ny forms non-increasing sequence then V m S n km
Bm+ 1,m I122 (1 e) mm (M + A1) and
(9)
m m -1
Bm44m=1.nm
MxX,X
< S (M 21)²
(10)
2
7

---

Lemma 2 bounds the norm of the matrix product B at two levels. The first result provides coarse bound, approximating linear and higher-order terms. The second result provides finer bound, preserving the linear term and approximating quadratic and higher-order terms. The proofs involve straightforward combinatorial expansion of B and are deferred to the Supplement Section (S.3
Approximating 11 (2 + xxx; requires 7) be small. Since this is recursive argument, we would need 1j.k to be small for i 1.... which is satisfied by the strong condition 77 is small. To obtain a tight analysis, we choose adaptively, We set k; Tmix (1/² (see definition in Eq 3
As we will show in detail in the Supplement, Lemma Eq 10 along with the adaptive choice k; gives sharp error bound. Using it, can bound Ti (see Eq as:
k+1
Ti S
Σ 7 j=2
T Br,k+2 [xxx]]) (I+ E) (X1 - Σ) Xx+2.... ................................
,,,,,,,,] Xn) B
Ti
(vp
Naively bounding the Ti, term by (1) leads to the same rate as downsampled Oja's 's algorithm.
In the following lemma, we will establish that, indeed, Ti.; has a much smaller norm. The novelty of our bound is not just in using the mixing properties of the Markov chain but also in teasing out the variance parameter V. We will state the lemma, slightly more general form as
Lemma 3. Under Assumptions
3 for < j S i k, LE [(x,X? - 2) X S 5n 5n1]11 - aAe
7
2
and
Ar)) ||S||2
where k; as defined Lemma 8.77 and S is constant symmetric positive semi-definite matrix.
Lemma bounds the norm of the covariance between matrices (x,X7 Σ) S and x,XT. In particular, this implies that the norm of Ti,j decays AAz (P) The proof uses spectral argument that replaces a coarse approximation by a sum of k2 ((1) terms to sum exponentially decaying terms, thereby removing the dependence on k;, which can be large as log(n). The proof is deferred to the Supplement Section S.4
Let {c1, C22 C3: C1} be positive constants for ease of notation. Coming back to Eq g we can bound T1
ci A2 (P))) + onikkt A similar argument can be applied to bound
as follows: Ti & an,k+2
[P()P]
T2 as: T2 S (2,kkt2 (V k3)) Putting everything together in
we have
(P)
(1+md1)
(1- (P))))
Vn an,k+2 +
CAll k
3
On,1
Onn,2
(P))
Error due approximation of matrix product
Recursion for ITD setting
Eruer due to Markavian dependence
Recursing on this inequality gives us our bound on E |u B, BT to present all our accompanying theorems.
v1l
(Theorem We are now ready
5 Intermediate Theorems for Convergence Analysis
In this section, we present our accompanying theorems which are used obtain the main result in Theorem 1 But before doing so, we will need to establish some notation. Let k; == Tmix (n²), and the step-sizes be set 17
(B with a, 8 defined in Theorem Let == 100 As
[A
shown in Lemma S.
Supplement Section S.3 our choice of step-sizes satisty, Vi E [n].
C.1 mi (M + A1) S
C.2 (Slow decay) S ni-k; S (1 + 2c) m S 2mh
Further, we define scalar variables
2 (1 + e) KnTn (M + A1),
:= 40k+41 (M 21)
V' = 1 1+ (3 + 4c) Ao (P) v
UK, = 6M [+332 (M 21)
(II)
1 Da (P)
8

---

and recall the definitions of B, and By, in Eqs and 7 respectively. We are now ready to present the theoretical results needed to prove our main result. For simplicity of notation, we present versions of the results by using 1/2 :- (A1-)a)(0++)) with B as defined in Theorem T However, these theorems in fact valid under more general step-size schedules. We state and prove the more general versions in the Supplement Section SS.4
Theorem 2. Under Assumptions
and for all n >> kn and nι satisfying C.I and C.2, we have,
B B S (1+ exp Σ
(2 + n/ (V* + Xi) + n (&ve.)
The three primary differences with the IID case are a) the (1 + r) term, which arises since the recursion sketched in Section llavees out the last k, terms which are bounded by (1 + r)² (b) the new factor of [P]] with D' due to the Markovian dependence between terms; and c) the extra lower order term T/ th., arising from the use of Lemmas and
Theorem 3. Let u min (t E [n].t K 2 0}. Under Assumptions 7 2 and for all n > u. and nn satisfying C.I and C.2,
n
E [Tr (VF Bn B Vi ) (1 5c) exp
Σ
21 + ni (V' x7) + ,ie Wh.d
x d Σ
(V' + mmbke) C2 exp
Σ iiutt1
2n (11 X2))
tu+11
where C₁2, = (1+ 2 exp (2A1 Σ 77))
Here, the difference is mainly in the new variable arising since the recursion stops u, not 1. (1 + 5c) represents the approximation of the first u terms.
Theorem Under Assumptions 7
3 for all n kn Ti satisfying C.I and C.2. and
: 2r
5
we have,
1000.
E
BRB
(1 5) exp
2m
4m
This differs from its IID counterpart by a multiplicative factor of (1 - s) for the same reason as before, which also makes the sums go up to (n k,) instead of Note that for sufficiently large
(Lemma.S., S.12 O log(n) is very small and 5 E (0, 1). Therefore, (1 (1s) s) Rs 1 as large
Theorem 5. Under Assumptions and 3 for all n >> kn and r7 satisfying C.1 and C.2, we have,
k.
n-k
(ff By Br un)
< + exp
Σ
4m +
Σ
The Sa,t
The differences are similar to the last theorems involving V1- Surprisingly, for this, the coarse approximation suffices, leading to an absence the term the bound. Having established these results, the final step is to substitute them into E follow the proof recipe described earlier. This requires significant calculations and deferred to the Supplement Section S.5
6 Experimental Validation
In this section, we present some simple experiments to validate our theoretical results. For more detailed experiments, see the Supplement Section S..6] We design Markov chain with ||/|= 10 states, where the transition matrix entries P.y equal p/[[5//-1] for of and 1 p for = j. Smaller values lead to larger mixing times. It can be verified that the stationary distribution
24 (92 uniform over the state-space and D2 (P) j (I p). We set 0.2 for Figures and 3a and vary in Figure 3b) Each point in the plot averaged over 20 random runs over different Markov chains, datasets, and initialization.

---

1.0
1-10
140
= =1-Se2
0.0
PuPt -1-
0.6
(uPt) =1-1e3
0.4
PuPhi Se
0.
(a) Variation of sin" error with 3
(b) Variation of sin" emor with [Az (P)
Figure 3: X axis represents the sample size, and Y axis represents the sin" error.
Each state E Ω is associated with (s) == Bernouli(pp.) distribution. We set d 1000 and select Ps ~ U (0, 0.05) at the start each random run. The covariance matrix, Mg, for each state set as E, (i,5) exp (- jlc,) 0,a, where c, ==
[2| os == 5 We start with the stationary distribution æ. and for each state st, draw IID samples Z; D (st). We standardize Z; such that all components have zero mean and unit variance under the state distribution, D (s:). We then generate the sample data-point for PCA as X, = E Z4. By construction, ED(,.) and E[X,] = 0 The step sizes for Oja's ' algorithm are set
[ = E, = (8+1)(A1-Aa] for a 5, B 1- (P) For the downsampled variant, every 10th Stat-pooia is considered, and 33 is accordingly divided by 10. For the offline algorithm, we recompute the leading cigenvector of the sample covariance matrix of data-points seen so far.
Figure T compares the performance of different algorithms for the Bernoulli distribution. Here, we checking if the results obtained in Theorem T Proposition [1] and Corollary T are reflected in the experiments. The experimental results demonstrate that Oja s algorithm performs significantly better than the downsampled version, consistent with the theoretical results. It also shows that Oja j algorithm performs similarly to the offline algorithm, which is also confirmed by our theoretical results and that of (28). Figure 3a compares the performance of Oja's algorithm for different covariance matrices. Smaller values 3 decrease the eigengap a - A2, and hence lead to a slower convergence. Figure 3b} confirms that smaller values of p (larger values of (A (P) |) also worsen the rate, which matches with our thcoretical results.
7 Conclusion
We have considered the problem of streaming PCA for Markovian data, which has implications in various settings like decentralized optimization, reinforcement learning, etc. The analysis of streaming algorithms in such settings has seen renewed surge of interest in recent years. However, the dependence between data-points makes it difficult to obtain sharp bounds. We provide. to our knowledge, the first near-optimal bound for obtaining the first principal component from a Markovian data stream that breaks the logarithmic barrier present the analysis done for downsampled data. We believe that the theoretical tools that we have developed in this paper would enable one to obtain sharp bounds for other dependent data settings, learning top principal components, and online inference algorithms with updates involving products of matrices.
8 Acknowledgements
We gratefully acknowledge NSF grants 2217069, 2019844, and DMS 2109155. We are also grateful to Rachel Ward, Bobby Shi, and Soumendu Sundar Mukherjee for useful discussions and to the anonymous reviewers for their valuable feedback.
10

---

References
[1] Zeyuan Allen-Zhu and Yuanzhi Li. First efficient convergence for streaming k-pca: a global. gap-free, and near-optimal rate. In 2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS), pages 487-492. IEEE, 2017.
[2] Jalaj Bhandari, Daniel Russo, and Raghav Singal. A finite time analysis of temporal difference learning with linear function approximation. CoRR, foss1806.02450. 2018.
[3] Minshuo Chen, Lin Yang, Mengdi Wang, and Tuo Zhao. Dimensionality reduction for stationary time series via stochastic nonconvex optimization. Advances in Neural Information Processing Systems, 31, 2018.
[4] Shuhang Chen, Adithya Devraj, Ana Busic, and Sean Meyn. Explicit mean -quare error bounds for conte-carlo and linear stochastic approximation. International Conference on Artificial Intelligence and Statisties, pages 4173 4183. PMLR, 2020.
[5] Thinh T. Doan, Lam M. Nguyen, Nhan H. Pham, and Justin Romberg. Convergence rates of accelerated markov gradient descent with applications in reinforcement learning, 2020.
[6] Thinh T. Doan, Lam M. Nguyen, Nhan H. Pham, and Justin Romberg. Finite-time analysis of stochastic gradient descent under markov randomness, 2020.
[7] Ron Dorfman and Kfir Yehuda Levy. Adapting mixing time in stochastic optimization with markovian data. In International Conference on Machine Learning, pages 5429-5446. PMLR, 2022.
[8] John C. Duchi, Alekh Agarwal, Mikael Johansson, and Michael I. Jordan. Ergodic mirror descent. SIAM J. Optim., 22(4):1549-1578, 2012.
(9] Alain Durmus, Eric Moulines, Alexey Naumov, Sergey Samsonov, and Hoi-To Wai. On the stability of random matrix product with markovian noise: Application linear stochastic approximation and td learning. In Conference on Learning Theory, pages 1711-1752. PMLR, 2021.
[10] Mathieu Even. Stochastic gradient descent under markovian sampling schemes. 2023.
[11] James Allen Fill. Eigenvalue Bounds on Convergence to Stationarity for Nonreversible Markov Chains, with an Application to the Exclusion Process. The Annuls of Applied Probubility, 1(1):62 - 87, 1991.
[12] Andreas Grammenos, Rodrigo Mendoza-Smith, Cecilia Mascolo, and Jon Crowcroft. Federated PCA with adaptive rank estimation. CoRR, hetp/9070080555 2019.
[13] Anne Hartebrodt, Reza Nasirigerdeh, David B. Blumenthal, and Richard Röttger. Federated principal component analysis for genome wide association studies. In 2021 IEEE International Conference on Data Mining (ICDM), pages 1090-1095, 2021.
[14] Amelia Henriksen and Rachel Ward. AdaOja: Adaptive Learning Rates for Streaming PCA. arXiv e-prints, page arXiv:1905.12111 May 2019.
[15] De Huang, Jonathan Niles-Weed. Weed, and Rachel Ward. Streaming k-pea: Efficient guarantees for oja's 's algorithm, beyond rank-one updates. CoRR, abs/2102.03644 2021.
[16] Prateek Jain, Chi Jin, Sham Kakade, Praneeth Netrapalli, and Aaron Sidford. Streaming pea: Matching matrix bernstein and opar-optimal finite sample guarantees for oja's algorithm. In Proceedings of The 29th Conference on Learning Theory (COLT), June 2016.
[17) Chi Jin, Sham M Kakade, Cameron Musco, Praneeth Netrapalli, and Aaron Sidford. Robust shift-and-invert preconditioning: Faster and more sample efficient algorithms for eigenvector computation. arXiv preprint arXiv: 15 10.08896, 2015.
[18] Harold J Kushner and G George Yin. Applications in signal processing, communications, and adaptive control. Stochastic Approximation and Recursive Algorithms and Applications, pages 63-93, 2003.

---

[19] David A Levin and Yuval Peres. Markov chains and mixing times, volume 107. American Mathematical Soc., 2017.
[20] Pascal Lezaud. Chernoff-type bound for finite Markov chains. The Annals of Applied Probability, 8(3):849 867, 1998.
[21] Robert Lunde, Pumamrita Sarkar, and Rachel Ward. Bootstrapping the error of oja's algorithm. Advances in Neural Information Processing Systems, 34:6240-6233 2021.
[22] Hanbaek Lyu, Deanna Needell, and Laura Balzano. Online matrix factorization for Markovian data and applications to Network Dictionary Learning. arXiv e-prints, page 11Xiv:1911.01931, November 2019.
[23] Shaocong Ma, Ziyi Chen, Yi Zhou, Kaiyi Ji, and Yingbin Liang. Data sampling affects the complexity of online sgd over dependent data. In Uncertainty Artificial Intelligence, pages 1296-1305. PMLR, 2022.
[24] Abdelkader Mokkadem. Mixing properties of arma processes. Stochastic Processes and their Applications, 29(2):309-315, 1988.
[25] Jean-Marie Monnez. Stochastic approximation of eigenvectors and eigenvalues of the q- symmetric expectation of a random matrix. Communications in Statistics-Thorry and Methods, pages 1-15, 2022.
[26] Nikos Mouzakis and Eric Price. Spectral guarantees for adversarial streaming pea, 2022.
[27] Dheeraj Nagaraj, Xian Wu, Guy Bresler, Prateek Jain, and Praneeth Netrapalli. Least squares regression with markovian data: Fundamental limits and algorithms. Advances in neural information processing systems, 33: 13:16666-1667, 2020.
[28] Joe Neeman, Bobby Shi, and Rachel Ward. Concentration inequalities for sums of markov dependent random matrices, 2023.
[29] Erkki Oja. Simplified neuron model a principal component analyzer. Journal of Mathematical Biology, 15(3):267-273, November 1982.
[30] Erkki Oja and Juha Karhunen. On stochastic approximation of the eigenvectors and eigenvalues of the expectation of a random matrix. Journal of mathematical analysis and applications, 106(1):69-84, 1985.
[31] Rayadurgam Srikant and Lei Ying. Finite-time error bounds for linear stochastic approximation andtd learning. In Conference Learning Theory, pages
PMLR, 2019.
[32] Tao Sun, Yuejiao Sun, and Wotao Yin. On markov chain gradient descent. Advances in neural information processing systems, 31, 2018.
[33] Joel A Tropp. User-friendly tail bounds for sums of random matrices. Foundations of computational mathematics, 12:389 2012.
[34] Lan V Truong. Generalization enror bounds on deep learning with markov datasets. Advances in Neural Information Processing Systems, 35:23452-23462 2022.
[35] John Tsitsiklis and Benjamin Van Roy. Analysis of temporal lemporal-difffereee learning with function approximation. In M.C. Mozer, M. Jordan, and T. Persche, editors, Advances in Neural Information Processing Systems, volume 9. MIT Press, 1996.
[36] Roman Vershynin. Introduction to the non- non-asymptotic analysis of random matrices. arXiv preprint arXiv:1011.3022) 2010.
[37] Geoffrey Wolfer and Aryeh Kontorovich. Estimating the mixing time ergodic markov chains. In Alina Beygelzimer and Daniel Hsu, editors, Proceedings of the Thirty-Second Conference on Learning Theory, volume 99 of Proceedings of Machine Learning Research, pages 31 3120-3159. PMLR, 25-28 Jun 2019.
12

---

[38] Puyudi Yang, Cho-Jui Hsich, and Jane-Ling Wang. History pea: A new algorithm for streaming pea. arXiv preprint arXiv:1802005447, 2018.
[39] Ingvar Ziemann and Stephen Tu. Learning with little mixing. S. Koyejo, S. Mohamed,
A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, Advances in Neural Information Processing Systems, volume 33, pages 4626 4637. Curran Associates, Inc., 2022.
13


Reviews:
[{"review_id": "QdxPcznlsW", "reviewer": "Reviewer_M3wf", "summary": "This paper analyzes the performance of Oja'a algorithm for streaming PCA when the samples come from an irreducible, aperiodic, and reversible Markov chain. In this work the authors avoid using downsampling to reduce dependency, and their result improves on that proposed in [3] by a poly-logarithmic multiplicative factor. ", "questions": "1. What is C in line 128? C also appear as a constant later. \n2. $\\pi$ is the stationary distribution of $X$ or $s$? In line 46 $\\pi$ seems to be the stationary distribution for $X$, while in section 2 $\\pi$ seems to be the stationary distribution for $s$. \n3. Assumption is not satisfied in several basic settings for example $X_i$ are Gaussian, which is a little disappointing. Can this assumption be relaxed using some additional arguments, for example, truncation? ", "limitations": "Yes. Societal impact not applicable. ", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "This paper gives a near optimal upper bound on the estimation error of Oja's algorithm for streaming PCA with Markovian data. Most past results consider iid data, while in many practical applications the data are dependent. Their result improve upon previous error rate bounds, and nearly matches an offline algorithm. Careful comparisons of the error rates are also provided in the paper. Leveraging the mixing properties of Markov chains, they present novel analysis that might find applications in other online problems. ", "weaknesses": "It seems that knowledge of the transition matrix is required to produce step sizes. In practice this must be estimated, it will be nice to include this as a part of the algorithm. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your kind words regarding the novelty of our analyses and wider applications for online problems.\n\n1. *Regarding the estimation of the mixing time*: [A] provides an offline algorithm for estimating the mixing time of Markov chains. They show that at least a constant times  $\\frac{\\tau_{\\text{mix}}|\\Omega|}{\\epsilon^{2}}$ samples are needed for estimating $\\tau_{\\text{mix}}$ within an absolute error $\\epsilon$, where $|\\Omega|$ denotes the size of the state space. Therefore, it is possible to use $O\\left(\\frac{1}{\\epsilon^{2}}\\right)$ samples to estimate the mixing time within accuracy $\\epsilon$. We will add this as a remark in the manuscript.\n\n\n\n\n2. Questions-\n\n    Q1. *Regarding $C$ on line 128*: $C$ here denotes the Markov Chain. We will remove this to avoid confusion.\t\n\n    Q2. *Regarding $\\pi$ being the stationary distribution of $s$ or $X$*: With $E_{\\pi}$, we meant the expectation over the stationary distribution of the states. We will replace this notation with $E_{s \\sim \\pi}E_{D(s)}\\left[.\\right]$ (the notation on line 140).\n\n    Q3. *Regarding truncation*: Thank you for asking this question. Indeed, we can extend our results to sub-gaussian data distributions with appropriate variance parameter decay. For space constraints and ease of exposition, we will provide an easy argument that changes the algorithm by considering \"truncated\" datapoints $\\\\{Y_i:= X_i1(\\lVert X_i \\rVert ^2\\leq \\alpha_n), i=1,\\dots, n\\\\}$, where the truncation parameter $\\alpha_n$ is set to be $c\\log n$. Thus we replace any datapoint whose squared norm exceeds the truncation parameter by zero. We provide an argument to show that the principal eigenvector of the covariance matrix of this truncated distribution is sufficiently close to that of the original one. Standard arguments can then show that the sin-squared error of the output of this algorithm w.r.t the original principal eigenvector has the same theoretical guarantee.\n\n    To be concrete, for each state, consider a sub-gaussian distribution, where $\\lVert X_{ij} \\rVert_{\\psi_2} \\leq \\nu_j$, and $\\sum_j\\nu_j<C$. Here $X_{ij}$ is the $j^{th}$ entry of random vector $X_i$, and $\\lVert.\\rVert_{\\psi_2}$ is the sub-Gaussian Orlicz norm (see [D]). A similar setup for unbounded distributions has been used in the IID setting by [B] and [C]. \n\n    We can show that this constrained distribution's covariance matrix $\\Sigma' := \\mathbb{E}\\left[Y_{i}Y_{i}^{T}\\right]$ is close in Frobenius norm to the original covariance matrix $\\Sigma = \\mathbb{E}\\left[X_{i}X_{i}^{T}\\right]$.\n\n    Next,  the Davis-Kahan theorem ([E]) can be applied to show that the population eigenvectors of $\\Sigma$ and $\\Sigma'$ are within a small sin-squared error, which can be made to be smaller than $1/n^b$ for a large $b$ (depending on the choice of $\\alpha_n$). Thus the sin-squared error of this algorithm is still $O\\left(\\frac{\\mathcal{V}}{n(1-|\\lambda_2\\left(P\\right)|)(\\lambda_1-\\lambda_2)^2}\\right)$ with high probability.  This is because the norm of the vectors (which is now $\\Theta(\\sqrt{\\log n})$) only appears in the higher-order error terms. We also believe that it is possible to analyze the original algorithm directly by looking at the intersection of any relevant event with the set $\\mathcal{A}_n=\\\\{\\max_i \\|X_i\\|^2\\leq \\alpha_n\\\\}$. Proposition 1 of [B] shows that $P(\\mathcal{A}_n)=1-o(1)$. We defer this argument for clarity.\n\nReferences\n\n[A] Wolfer, Geoffrey, and Aryeh Kontorovich. Estimating the mixing time of ergodic Markov chains. Conference on Learning Theory. PMLR, 2019.\n\n[B] Lunde, R., Sarkar, P. and Ward, R., 2021. Bootstrapping the error of Oja's algorithm. Advances in Neural Information Processing Systems, 34, pp.6240-6252.\n\n[C] Liang, X., 2023. On the optimality of the Oja's algorithm for online PCA. Statistics and Computing, 33(3), p.62.\n\n[D] R. Vershynin. High-Dimensional Probability. Cambridge University Press, Cambridge, UK,\n2018.\n\n[E]. Yu, Yi, Tengyao Wang, and Richard J. Samworth. A useful variant of the Davis–Kahan theorem for statisticians. Biometrika 102.2 (2015): 315-323."}, {"review_id": "kONVQDJr4s", "reviewer": "Reviewer_1vdj", "summary": "This paper focuses on a particular algorithm for PCA known as the Oja's algorithm. \nThis is an iterative algorithm that operates in a streaming manner. \nThe goal is to, when the data stream is Markovian, estimate the principal component of the covariance matrix of the unknown covariate vectors. \nThe main question that the authors addressed is to obtain sharp (up to multiplicative constant) rate guarantee for such an algorithm. \nThere are miscellaneous results in other parts of the paper that compare favorably to the main result of the paper. \nI found the most impressive aspect of this paper to be the fact that authors managed to shave the log factor in the rate (compared to the down-sampled version commonly used in Markovian setting). ", "questions": "1. The sin^2 guarantees in various theorem statements are somewhat heavy. Does it make sense to the authors to state a more compact bound (say using big-O notation) in the main text and defer the exact bound to appendices? On the other hand, I do understand that the authors would like to present completely nonasymptotic bounds...\n\n2. I did not check the proof carefully, but my impression upon scanning the appendices is that the proof seems to rely heavily on the simplicity of the form of the Oja's iteration. \nFor instance, there're intimate relation between this work and a line of work on concentration for random matrix products. \nGoing beyond Oja's iteration, it possible to extrapolate the proof techniques to other iterations where there are nonlinearities?\nE.g., in generalized power iteration, approximate message passing and so on, it's important to pick judicious nontrivial nonlinearities. \nOf course, not all of these algorithms are applicable in the Markovian setting, but I'm sure the authors are aware of applicable examples. \n\n3. The authors have already discussed the down-sampling trick in details. I would like to add a few more examples where there's loss of log factors using down-sampling. Bresler and his students have a line of work on inference in the presence of temporal correlation, notably this one https://arxiv.org/abs/2006.08916 in the context of linear regression. The idea there is precisely down-sampling and they showed down-sampling is optimal *up to log factors*. Clearly the authors of this submission is not satisfied with log factors. But it's still interesting to see similar phenomenon being recurrently observed. \nAnother \"example\" is in the setting of parameter estimation from Gaussian mixtures where the mixing coefficients form a binary Markov chain; see https://arxiv.org/abs/2206.02455. \nOver there, curiously, one cannot waste samples. Instead, one takes coherent combination within blocks of size being roughly the mixing time of the chain. ", "limitations": "N/A", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The result of this paper seems rather satisfactory and tells a relatively complete story. \nThe main text is somewhat dense but it's never getting obscure. \nContributions, techniques to achieve them and difficulties for doing so are clearly explained. ", "weaknesses": "I don't see obvious weaknesses. \nSo I'll just propose a few minor suggestions. \n\nI personally found it useful to keep the following binary Markov chain in mind when reading the results and proof sketches: the states are -1, +1 and the transition kernel is given by a binary symmetric channel with flip probability delta. \nThen the spectral gap is delta (perhaps up to constant) and the mixing time is essentially 1/delta. \nI encourage the authors to use this as a running example in the proof sketch, though I don't strongly insist on this. \n\nIf I understand correctly, the Markov chain does not have to be finite-state right?\nIf so, this is perhaps worth commenting on, if it has not been done yet. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your encouraging words about removing the logarithmic dependence in our rate. \n\n1. *Regarding the example Markov chain*: Thank you for this suggestion. Due to the interpretable relationship between the spectral gap and the mixing time, we do use a similar Markov chain for our experimental setup (See lines 302-308). For added clarity, as per your suggestion, we will also try to include this in the proof sketch as a running example.\n\n2. *Regarding finite state-space Markov chains*: We believe that our analysis should extend to the general state-space Markov chains, similar to the extensions in Section 4 in [A]. However, for ease of exposition, we have focused on the finite-state-space case in this paper.\n\n3. Questions \n\n    Q1. *Regarding the exposition of theoretical guarantees*: We will add another version of the statement of Theorem 1 with reasonable values substituted for the learning rate. We do believe this will help the exposition.\n\n    Q2. *Regarding extrapolation to nonlinear updates*: We are currently working on the concentration of matrix products in the Markovian setting, which is an open problem to the best of our knowledge. Even in the IID case, non-linearities in the update may make the analysis deviate significantly from the Oja-type update, as in the case of streaming sparse PCA (see [B]).\n\n    Q3. *Regarding loss of log factors using down-sampling*: We are grateful to the reviewer for pointing out these references. We tried implementing the setting of Bressler et al.’s work [25] without discarding samples along with a decaying step size. Indeed, the algorithm that looked at the entire dataset performed better. Therefore, we believe this is a common phenomenon across multiple settings. We were unaware of the second paper and will cite it in our manuscript.\n\nReferences\n\n[A]. Lezaud, Pascal. Chernoff-type bound for finite Markov chains. Annals of Applied Probability (1998): 849-867.\n\n[B]. Yang, W. and Xu, H., 2015, June. Streaming sparse principal component analysis. In International Conference on Machine Learning (pp. 494-503). PMLR."}, {"review_id": "YE8OmxpVEC", "reviewer": "Reviewer_iipj", "summary": "This paper studies the problem of streaming PCA, where the data points are revealed to the algorithm one at a time and the goal is to estimate the top eigenvector of the underlying covariance matrix.\nCrucially, the data points are sampled from a Markov chain and thus do not satisfy the typical assumption of being independent (they are still identically distributed as the Markov chain is stationary).\nThe paper analyzes the Oja's algorithm, which has shown a great empirical (and theoretical) success in streaming PCA with i.i.d. data. \nCrucially, the paper does not do any downsampling of the data, which is a common practice in the literature to curb the influence of dependence in the data.\n\n----\n# Update\nI thank the authors for their detailed responses. Based on their promise to include a nuanced discussion on the results in the final version of the paper (along with other presentation improvements suggested by the reviewers), I have now increased my score to 7 from 5.", "questions": "# Questions\n\n1. Why are the error rates in the paper claimed to be sharp? Although it is reasonable to expect that a dependence on the mixing time appears, but does it have to be a multiplicative dependence? Are there lower bounds?\n\n2. Does the algorithm need to know the mixing time of the chain to declare victory? That is, Theorem 1 requires n to be larger than the mixing time. If yes, is it inherent?\n\n3. Isn't Corollary 1 immediate from the the i.i.d. result by using downsampling argument (without using Theorem 1)?\n\n# Other Comments\n\n1. Theorem 1 is too dense and difficult to parse. \n\n   i. Line 158 says that the gap $1- \\lambda_2(P)$ is inversely proportional to $\\tau_\\mathrm{mix}$. Then, why does Theorem 1 use both of them in the statement? It would be simple to use $\\tau_\\mathrm{mix}$.\n\n   ii. The notation with $\\alpha$ and $\\beta$ is convoluted. For the version in the main theorem, I suggest choosing a good enough value of $\\alpha$ and present a simplified choice of $\\beta$.\n\n2.  Line 194: Only the upper bounds obtained are sub-optimal. No lower bound has been shown.\n\n3. Section 5: What is the connection between the series of theorem statements.\n\n\n4. There are many typos and grammatical errors in the paper and many articles are missing. Several mathematical quantities are described after they are used. For example,\n    i. Line 31:  w.r.t.\n\n    ii. Line 31:  the true\n\n    iii. Line 36: v_1 is not defined.\n\n    iv. Line 75: What is meaning of 'stationary distribution of s\n\n    v. Assumption 1: What exactly is the last assumption is ambiguous.\n\n    vi. What norms are used in Assumptions 2 and 3?\n\n    vii. Line 139: WLOG\n\n    viii. Line 140: Why is it without loss of generality?\n\n    vix. Line 142: Braces around the equation are missing.\n\n    x. Section 2.1: Give specific references to the theorem statements in the book [18]. \n\n    xi. Line 181: Eq 1 should be Theorem 1.\n\n    xii.Line 182: Which remark in [15]?\n\n    xiii. Line 227 and 230: $\\alpha_{n,2}$ is not defined.\n\n    xiv. Line 228: $\\beta_{n,2}$ is not defined.\n", "limitations": "What are some instances where $\\mathcal{V}$ is much smaller than $\\mathcal{M}$. Using some simple calculations, I am getting similar expressions for them in the case of Gaussian data.", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 1, "contribution": 3, "strengths": "The paper studies an important problem and the paper's results are novel. In particular, they seem to improve over the de-facto approach of downsampling. I like the qualitative aspects of the paper's results.\n", "weaknesses": "In terms of weaknesses, I have the following comments:\n\n1. (Writing) The paper's writing is sadly not up to the mark. In addition to the typos (some of which are listed below), theorem statements need to be simplified and the paper (especially Sections 4 and 5) need major re-writing. Even after spending a lot of time reading through the technical proof sketches in Sections 4 and 5, I was a bit clueless. \n\n2. (Gap-free bounds) There are no results for the gap-free setting, which is needed in several settings. I recommend adding them to strengthen the paper.\n\n3. (Dependence on failure probability) At several places, the paper claims to match the existing bounds in the literature for the i.i.d. setting. However, the dependence on the failure probability in the paper is much worse, polynomial as opposed to logarithmic in [15]. \n\n4. (Non-zero mean setting) Can the results be generalized to the zero mean setting? In the iid setting, one can simply take a pair of data and take their difference as the sample. Can one do the same trick here in the streaming setting (without compromising the mixing time)? ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your kind words regarding our results' novelty and the importance of the problem.\n\n1. *Regarding writing*: We will correct the typographical errors and incorporate the suggestions. We will also add a simplified version of Theorem 1 with reasonable values substituted for the parameters.\n\n2. *Regarding Gap-free bounds*: To our knowledge, this is the first sharp rate for Oja's algorithm for the setting with an eigengap. The gap-free bounds, although important in various scenarios, are beyond the scope of this work.\n\n3. *Regarding dependence on failure probability*: We are comparing against Theorem 4.1 from Jain et al. [15], which has an **identical dependence** on the failure probability as us, i.e., $\\frac{\\log(\\frac{1}{\\delta})}{\\delta^2}$. Hence our dependence on the failure probability *exactly matches* the existing results in the IID case.\n\n4. *Regarding the non-zero mean case*: Thank you for raising this insightful point. Indeed, in the Markovian setting, this construction leads to a modified Markov chain with the same mixing time. Let the modified mean-zero dataset be $\\\\{Y_i , i = 1,\\dots,\\lfloor \\frac{n}{2} \\rfloor\\\\}$ where $Y_{i} := X_{2i-1} - X_{2i}$. This data comes from another Markov chain with states $(a,b)$ such that $P(a,b)>0$ in the original Markov chain. The new transition matrix $P'$ is $P'((a,b),(c,d))=P(b,c)P(c,d).$ When the original Markov chain is irreducible and aperiodic so is this one. The stationary distribution of $P'$ is $\\pi'((a,b))=\\pi(a)P(a,b)$. Finally, to see the relationship between the mixing times, observe that, $P(s'_t=(c,d)|s'_0=(a,b)) = P(s_t=c|s_0=b)P(c,d)$.  \n\n   Thus, \n\n    $\\sup_{(a,b)}\\sum_{(c,d)}|P(s'_t=(c,d)|s'_0=(a,b))-\\pi(c)P(c,d)|$  \n\n    $=\\sup_{b}\\sum_c|P(s_t=c|s_0=b)-\\pi(c)|$.\n\n    Therefore, the mixing time stays unchanged. Thus our analysis applies to non-zero mean distributions similar to the IID case.\n\n5. Questions \n\n    Q1: *Regarding the sharpness of our results*: Neeman et al. [26] obtain an upper bound for offline estimation of the principal eigenvector of the empirical covariance matrix of a Markov chain. Their bound is claimed to be optimal up to logarithmic factors. Proposition 1 (line 177) provides this bound in our setting. Our error bound for streaming PCA improves this $\\log(d)$ factor. Thus, in the same sense that Jain et al.’s result in [15] is sharp, we have claimed ours to be sharp. The multiplicative factor of $\\tau_{\\text{mix}}$ also arises in lower bounds for various convex optimization problems with Markovian data (see [8, 25]).\n\n    Q2:  *Regarding knowledge of the mixing time*:  $\\tau_{\\text{mix}}$ is needed to set the learning rate optimally. However, [A] provides an offline algorithm for estimating the mixing time of Markov chains. They show that at least a constant times  $\\frac{\\tau_{\\text{mix}}|\\Omega|}{\\epsilon^{2}}$ samples are need for estimating $\\tau_{\\text{mix}}$ within an absolute error $\\epsilon$ where $|\\Omega|$ denotes the state space size. Therefore, one may use $O\\left(\\frac{1}{\\epsilon^{2}}\\right)$ samples to estimate the mixing time within accuracy $\\epsilon$. We will add this as a remark.\n\n    Q3: *Regarding Corollary 1 being immediate from the i.i.d. result*: Intuitively, this should be the case. However, even after downsampling, the data points are not completely independent, so IID results do not apply immediately. For example, the data drop algorithm in [25] for online linear regression for Markovian data also requires a significant amount of analysis to handle dependence. \n\n6. Other Comments\n\n   (a) *Regarding Theorem 1 presentation* : \n   Theorem 1 is presented in the most general form for easy comparison with the IID version in [25]. We will add a version including the suggestions. \n\n    *Regarding $\\tau_{\\text{mix}}$ and $1 - |\\lambda_2\\left(P\\right)|$ in Theorem 1*: The relationship between $\\tau_{\\text{mix}}$ and $1 - |\\lambda_2\\left(P\\right)|$ for reversible Markov chains (Eq~S.58, supplement line390) gives rise to an additional factor logarithmic in the size of the state-space, which is why we used both.\n\n   (b) *Regarding sub-optimality shown only for upper bounds (Line 194)*: We will clarify this in the final manuscript.\n\n   (c) *Regarding the connection between theorem statements (Section 5)*: Lines 200-210 show how the results in Section 5 help in proving the main result.\n\n   (d) Clarifications:\n\n   (iii) *(Line 36) Regarding $v_{1}$ not being defined*: $v_1$ is defined in line 37, which will be moved earlier. \n\n   (iv)  *(Line 75) Regarding the meaning of stationary distribution of $s$*: Stationary distribution of $s$ means $\\pi(s)$.\n\n   (v) *Regarding the last assumption being ambiguous in Assumption 1*:  We mean the requirement for starting at stationarity in Assumption 1.\n\n   (vi) *Regarding norms used in Assumptions 2 and 3*: A in lines 143 and 144, $\\lVert.\\rVert_2$ denotes the Euclidean norm for vectors and the operator norm for matrices. \n\n   (viii) *(Line 140) Why is it without loss of generality*: We assumed $M+\\lambda_1\\geq 1$ because when it is less than 1, it only makes our bounds tighter.\n\n   (xii) *(Line 182) Which remark in [15]?*: We made an error in calling this a remark while it is footnote 1 on page 2.\n\n   (xiii, xiv) *(Lines 227, 228, and 230)*: $B_{n,i}$ is defined on Line 217 (Eq 7). $\\alpha_{n,i}$ is defined as $\\alpha_{n,i} := \\mathbb{E}\\left[v_1^TB_{n,i}B_{n,i}^Tv_1\\right]$. PDF document in the global author rebuttal also contains a pictorial representation of $\\alpha_{n,2}$.  We could not find a $\\beta_{n,2}$ notation anywhere. \n\n7. Limitations: \n\n   For Bernoulli random variables with small p, $\\mathcal{V}$ is much smaller than $M^2$. We show this for a one-dimensional $Y\\sim Bernoulli(p)$ as a proof of concept. Let $X=Y-p$. So $|X^2-p(1-p)|\\leq 1$, i.e. $M=1$ and $\\mathbb{E}(X^2-p(1-p))^2=O(p)$ for small $p$.\n \nReferences\n\n[A]. Wolfer, Geoffrey, and Aryeh Kontorovich. Estimating the mixing time of ergodic Markov chains. COLT 2019."}, {"review_id": "X17iojpD2n", "reviewer": "Reviewer_eXbt", "summary": "The submission considers the task of *streaming* PCA, estimating the largest eigenvector of a covariance matrix in an online manner as new data points arrive. The focus is on theoretical analysis of a classical method, Oja's Algorithm, in the hitherto unconsidered non-IID setting. \n\nA seemingly strong restriction is made that the data come from a chain that is not only ergodic, but also reversible. The latter seems a potentially strong condition, but the justification appears as this is exactly the setting arising in algorithms in Federated PCA. \n\nThe primary contribution appears to be a theoretical justification for using Oja's algorithm as per normal in the considered setting, as opposed to \"downsampling\" points to decrease dependence between them. A multiplicative factor of how optimal error rates differ in the setting considered from the IID setting is also established. ", "questions": "\n* The title of the submission in my opinion would benefit from being changed to being clearer about the specific setting considered. At present it masks the fact that the data are from a reversible chain, which could be quite niche.\n\n* Suggestion: It would be clearer to say from the outset (in the abstract even) that the points are assumed to come from a reversible process in stationary. \n\n* Why was the $\\sin^2$ error used in the analysis, as opposed to say, $L_p$ norms? Is it standard in analysis of such algorithms, if so, why? The paper could benefit from discussion of this aspect.\n\n* Issue: What are $\\mathcal{V}$ and $\\mathcal{M}$? These are not defined before their introduction (or even after, it seems), which is a serious problem.\n\n* The motivation for the work is indeed Metropolis-Hastings (for example) output, which justifies the reversible setting, but how much more difficult is the problem without it? Would one expect the result to change? A small discussion of this could possibly be beneficial.\n\n* Abstract, l2: \"principle\" should be \"principal\".\n\n* p2, l47: \"*the second eigenvalue in the magnitude of the transition matrix $P$*\" reads awkwardly and incorrectly (it sounds like the magnitude is taken wrt the matrix, not the eigenvalues).\n\n* Suggestion: Figure 1 would be easier to parse if the downsampling rate is mentioned in the legend.  \n\n* p5, l192: \"atleast\" should be \"at least\".\n\n* p7, equation after line 250 is significantly out of the page margins. ", "limitations": "For the most part yes, but additional clarity (as per other comments) regarding certain aspects (including the numerical experiments) would be beneficial for the paper to be up to NeurIPS standard in my opinion. ", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The theoretical problem appears to be tackled in a manner involving some level of sophistication, and the results of potential interest. As the authors themselves state, many intermediate results (as contained in the supplement) may be of wider interest. \n\nThat said, whilst I am familiar with many of the individual constituent elements of this work, I am not familiar with the dependent streaming PCA literature so much, so it is difficult for me to assess the contribution. \n", "weaknesses": "* The submission is, on the whole, not particularly well written in my opinion. It could benefit from another round of two of editing and polish. The implications of *why* the particular problem are being studied and clearer discussions of its potential use would be beneficial. There are also some obvious typos and undefined notations.\n\n* The problem considered appears quite niche and it is not made clear in the paper the extent to which the results are of wider interest (a cursory mention is made of reinforcement learning in the conclusion, but there is no discussion of this).\n\n* As the motivation for the paper is for settings where the output comes from MCMC algorithms, I would expect an experiment exploring the effect of different burn in periods and how well Oja's method performs in practical settings and when the chain has converged to what is generally considered close enough to its stationary distribution for iterates to be used in practice (see comment below about burn in). \n\n* As far as I can tell, simply downsampling as done in the experiment is not equivalent to the cited Chen et al. (2018) approach, who propose further refinements to the algorithm (see their Algorithm 1). As the papers primary contribution is supposedly justification for using Oja's algorithm as is over other approaches for dependent data, I would expect at a minumum a comparison is made to the actual Chen algorithm (this includes for the experiments in the supplement), as well as any other recent competitors. \n\n* The assumption of requiring the data come from a reversible process appears somewhat strong to be of general interest. I understand the specific application, but the authors should discuss how/if this holds in a wider setting. That is, under what conditions are multivariate time series reversible? Does this hold for certain wide cases of interest? \n\n* Issue: Page 1, footnote: It is not a fair assumption to assume the assumption holds after a sufficient period of burn in. Indeed, a burn in period may suffice to ensure that a Markov chain has burned in long enough to ensure future iterates have distribution that are sufficiently *close* to the stationary, but in general the chain never actually reaches it. Unless an analysis is provided that accounts for some initial distance to stationarity and mixing rate of the chain, I think it is a tad disingenuous to claim that the requirement for stationarity is a non-issue. I don't think not starting in stationarity is a condition that must be considered theoretically, but I **do** think at least one experiment exploring the effect of different burn ins in practical settings is required for the paper to tell a cohesive story. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for commenting on the novelty of our analyses and the fact that they can find applications in other online problems. We provide answers to specific queries below.\n\n1. *Regarding wider interest*:\nLines 19-32 mention the connection to distributed PCA where the Markov chain arising from a random walk on a connected undirected graph is naturally reversible. Furthermore, as Reviewer 1vdj pointed out, our analysis is inherently connected to the concentration of random matrix products, which has only been studied in the IID setting. We believe that the tools and results in this paper would help establish results in the Markovian case.\n\n2. *Regarding burn-in*: Bressler et al. [25], Chen et al. [3], etc., also make the stationarity assumption. We will remove the footnote to avoid further confusion. Please see the plot in the PDF document attached with the global author rebuttal for the detailed experiment observing the effect of different burn-in periods. In general, we observed that the higher burn-in periods lead to slightly faster error decay. \n\n3. *Regarding further refinements to Chen et al. 2018's algorithm*: We have rechecked this. The only difference we see is using the $\\Pi_{\\text{orth}}$ operator, which (as pointed out by Chen et al.) is the same as the orthogonalization step of Oja’s algorithm. We are also unaware of any recent competitors for streaming PCA in the Markovian setting. We would be very grateful if the reviewer points us to the refinement and other competitors they have in mind.\n\n4. *Regarding reversibility*: Thank you for raising this important point. We assumed reversibility to simplify our analyses. However, one can extend our results to the non-reversible case following reversibilization tools (also see [A] and [B]). We use reversibility only in Lemma 3, which is then used by all other theorems. Here we look at the time-reversed chain, which is identical to the original Markov chain in the reversible case. The probability transition matrix $P_{\\*}$ of the time-reversed Markov chain is given by $P_{\\*}(x,y) = \\frac{\\pi(y)}{\\pi(x)}P(y,x)$. It can also be expressed as $P_{\\*}:= \\Pi^{-1}P^{T}\\Pi$, where $\\Pi$ is a diagonal matrix with $\\Pi(x,x)=\\pi(x)$.  In the non-reversible case, we consider the matrix $R_{\\*}:=P_{\\*}P$, which is also known as the multiplicative reversibilization of $P$ (see [B]). It corresponds to the transition matrix of a reversible Markov chain and has real and nonnegative eigenvalues. In a brief sketch, we show how this allows us to replace the eigengap of the transition matrix $P$, i.e., $1-|\\lambda_2(P)|$ in our current bound by $1-\\sqrt{|\\lambda_2(P_{\\*}P)|}$. \n\n   The proof of Lemma 3 essentially reduces to computing the norm of the matrix $Q$ (Supplement line 97): $  Q:= \\Pi^{\\frac{1}{2}}\\left(P^{t}-11^{T}\\Pi\\right)\\Pi^{-\\frac{1}{2}}. $\n   For a non-reversible Markov chain, the original transition matrix $P$ is replaced by the transition matrix $P_{\\*}$ of the time-reversed Markov chain: $Q_{\\*} := \\Pi^{\\frac{1}{2}}\\left(P_{\\*}^{t}-11^T\\Pi\\right)\\Pi^{-\\frac{1}{2}}$.\n\n   We will show how the norm of $Q_{\\*}$ is computed in the non-reversible case. Note that $P_{\\*}^t=\\Pi^{-1}(P^t)^T\\Pi$, and $\\Pi^{\\frac{1}{2}}\\left(P_{\\*}^{t}P^{t}\\right)\\Pi^{-\\frac{1}{2}}$ is a symmetric matrix, with first eigenvector $\\Pi^{1/2}1$.\n\n   Define: $ Q_{\\*}Q_{\\*}^T = \n     \\Pi^{\\frac{1}{2}}\\left(P_{\\*}^{t}P^{t} -11^{T}\\Pi\\right)\\Pi^{-\\frac{1}{2}}$. Then, we finally have,\n    $ \\lVert Q_{\\*}\\rVert = \\sqrt{\\lVert Q_{\\*}Q_{\\*}^{T}\\rVert } = \\sqrt{| \\lambda_{2}\\left(\\Pi^{1/2}P_{\\*}^{t}P^{t}\\Pi^{-1/2}\\right)|} = \\sqrt{|\\lambda_{2}\\left(P_{\\*}^{t}P^{t}\\right)|} \\leq \\left(\\sqrt{|\\lambda_{2}\\left(P_{\\*}P\\right)|}\\right)^{t},$ \n   where the last inequality follows from Section 3.2 of [A].\n\n5. *Regarding reversibility conditions for multivariate time series*: Stationary Gaussian autoregressive models are time-reversible (see [C]).\n\n6. Questions: \n\n    Q1: *Regarding the title of the paper*:\n\t We will add a remark stating how our argument extends to the non-reversible case.\n\n    Q2: *Regarding making the reversibility assumption clearer*: The submitted manuscript does mention the reversibility assumption in the abstract (see line 4), the introduction (see line 43), and all the theoretical results.\n\n    Q3: *Regarding the use of the $\\sin^{2}$  error as opposed to $L_{p}$  norms*: The $\\sin^{2}$ error is always used in the related literature since that is a standard metric to measure the distance between two unit vectors. We also want to point out that $\\sin^2(v_1,\\hat{v})= \\frac{1}{2}\\lVert v_1v_1^T-\\hat{v}\\hat{v}^T\\rVert_F^2$.\n\n    Q4: *What are $M$ and $\\mathcal{V}$?* \n\t $\\mathcal{V}$ and $M$ are defined in Assumptions 2 (line 137) and 3 (line 138), respectively.\n\n    Q5: *Regarding reversibility*: Please see our discussion to *point 4* about the extension to the non-reversible case.\n\n    Q6-10: Thank you for pointing out these typographical issues. We will update our manuscript to reflect these changes.\n\nReferences \n\n[A]. Lezaud, Pascal. Chernoff-type bound for finite Markov chains. Annals of Applied Probability (1998): 849-867.\n\n[B]. J. A. Fill. Eigenvalue bounds on convergence to stationarity for non-reversible Markov chains,\nwith an application to the exclusion process. The Annals of applied probability pages 62–87,\n1991. \n\n[C].  A. J. Lawrance, Directionality and Reversibility in Time Series, International Statistical Review / Revue Internationale de Statistique, Vol. 59, No. 1 (Apr., 1991), pp. 67-79"}]

Author rebuttal:
We want to first thank all the reviewers for their valuable suggestions and insightful feedback. We believe we have addressed nearly all of their main technical questions. We give brief sketches to show how to extend our analysis to **non-reversible** Markov chains, the **non-zero mean** data settings, and **sub-gaussian data distributions**.  In what follows, we will address some important points each reviewer has raised. We will correct all the typographical issues pointed out and will not address them here.
 
1. *Regarding reversibility (Reviewer eXbt)*:
We assumed reversibility to simplify our analyses. However, one can extend our results to the non-reversible case following reversibilization tools (also see [I] and [II]).  We use reversibility only in Lemma 3, which is then used by all other theorems. Here we look at the time-reversed chain, which is identical to the original Markov chain in the reversible case. The probability transition matrix $P_{\*}$ of the time-reversed Markov chain is given by $P_{\*}(x,y) = \frac{\pi(y)}{\pi(x)}P(y,x)$. In the non-reversible case, we consider the matrix $R_{\*}:=P_{\*}P$, which is also known as the multiplicative reversibilization of $P$ (see [II]). It corresponds to the transition matrix of a reversible Markov chain and has real and nonnegative eigenvalues. This allows us to replace the eigengap of the transition matrix $P$, i.e., $1-|\lambda_2(P)|$ in our current bound by $1-\sqrt{|\lambda_2(P_{\*}P)|}$. In our response to Reviewer eXbt, we provide more details.

2. *Regarding non-zero mean data (Reviewer iipj)*: Like the IID case, we can also deal with non-zero mean data by taking the difference of disjoint pairs of datapoints as our new data. This works for the streaming Markovian setting without increasing the mixing time. In our response to Reviewer iipj, we show that this essentially reduces to a modified Markov chain, which shares many important properties of the original Markov chain. Define the modified mean zero dataset as $\\{Y_i , i = 1,\dots,\lfloor \frac{n}{2} \rfloor \\}$ where $Y_{i} := X_{2i-1} - X_{2i}$. It is not hard to argue that this dataset is coming from another Markov chain with states $(a,b)$ such that $P(a,b)>0$ in the original Markov chain. We will denote the transition matrix of this Markov chain by $P'$, which is defined as follows: $P'((a,b),(c,d))=P(b,c)P(c,d) $. This is a valid transition matrix, and when the original Markov chain is irreducible and aperiodic, so is this one. For an argument on why the mixing time stays unchanged, please see our response to Reviewer iipj.
 
3. *Regarding sub-gaussian data distributions (Reviewer M3wf)*: 
 One can extend our results to sub-gaussian data distributions for each state using a truncation-based argument as the reviewer points out. We consider sub-gaussian distributions with appropriate variance parameter decay similar to the setup used in the IID setting by [IV] and [V]. In our response to Reviewer M3wf, we show a roadmap to an argument, which makes a small modification to the algorithm by considering $\\{X_i1(\lVert X_i \rVert ^2\leq \alpha_n), i=1,\dots, n\\}$, where the truncation parameter $\alpha_n$ is set to be $c\log n$. Proposition 1 of [IV] shows that the probability of the largest squared norm of $n$ sub-gaussian vectors exceeding $\alpha_n$ is $o(1)$. Using this, we can extend our method to incorporate sub-gaussian distributions because the truncation parameter with logarithmic dependence on $n$ only appears in the higher-order error terms.

4. *Regarding estimation of mixing time (Reviewers M3wf, iipj)*:  [III] provides an offline algorithm for estimating the mixing time of Markov chains. They show that at least a constant times  $\frac{\tau_{\text{mix}}|\Omega|}{\epsilon^{2}}$ samples are need for estimation of the $\tau_{\text{mix}}$ within an absolute error $\epsilon$, where $|\Omega|$ denotes the size of the state space. Therefore, it is possible to use $O\left(\frac{1}{\epsilon^{2}}\right)$ samples to estimate the mixing time within accuracy $\epsilon$. We will add this as a remark in the manuscript.
 
5. *Regarding burn-in (Reviewer eXbt)*:  Bressler et al. [25] and Chen et al. [3], etc. also make the stationarity assumption when analyzing Markovian data. Please see the attached pdf for experiments observing the effect of different burn-in periods when starting from non-stationary distributions. We observe that the longer burn-in periods lead to a slightly faster error decay.

6. *Regarding presentation of Theorem 1 (Reviewers 1vdj, iipj, eXbt)*: We will add a simpler version of Theorem 1 with reasonable values substituted for the parameters for ease of exposition.
 
References
 
[I] Lezaud, Pascal. Chernoff-type bound for finite Markov chains. Annals of Applied Probability (1998): 849-867.
 
[II] J. A. Fill. Eigenvalue bounds on convergence to stationarity for non-reversible Markov chains,
with an application to the exclusion process. The Annals of applied probability pages 62–87,
1991. 
 
[III] Wolfer, Geoffrey, and Aryeh Kontorovich. Estimating the mixing time of ergodic Markov chains. Conference on Learning Theory. PMLR, 2019.
 
[IV] Lunde, R., Sarkar, P. and Ward, R., 2021. Bootstrapping the error of Oja's algorithm. Advances in Neural Information Processing Systems, 34, pp.6240-6252.
  
[V] Liang, X., 2023. On the optimality of the Oja’s algorithm for online PCA. Statistics and Computing, 33(3), p.62.


Meta-review: All reviewers appreciate the important problem of analyzing algorithms in the Markovian setting as opposed to iid data. The paper gives a nearly tight bound without losing a log factor for Oja's algorithm. This log factor is common in bounds for many algorithms in the Markovian setting due to a simpler analysis using down-sampling. On the downside, the paper needs to improve the presentation, which the authors promised in the rebuttals.
Decision: Accept (spotlight)