Bootstrapping the error of Oja’s algorithm
Year: 2021

Abstract:
We consider the problem of quantifying uncertainty for the estimation error of the leading eigenvector from Oja’s algorithm for streaming principal component analysis, where the data are generated IID from some unknown distribution. By combining classical tools from the U-statistics literature with recent results on high-dimensional central limit theorems for quadratic forms of random vectors and concentration of matrix products, we establish a weighted χ approximation result for the sin error between the population eigenvector and the output of Ojas algorithm. Since estimating the covariance matrix associated with the approximating distribution requires knowledge of unknown model parameters, we propose a multiplier bootstrap algorithm that may be updated in an online manner. We establish conditions under which the bootstrap distribution is close to the corresponding sampling distribution with high probability, thereby establishing the bootstrap as a consistent inferential method in an appropriate asymptotic regime.

heading: 1 Introduction
text: Since its discovery over a century ago [13], principal component analysis (PCA) has been a cornerstone of data analysis. In many applications, dimension reduction is paramount and PCA offers an optimal low-rank approximation of the original data. PCA is also highly interpretable as it projects the dataset onto the directions that capture the most variance known as principal components.
Important applications of PCA include image and document analysis, where the largest few principal components may be used to compress a large dimensional dataset to a manageable size without incurring much loss; for a discussion of some other applications of PCA, see for example, [28]. In these settings, the original dimensionality, which could be the number of pixels in an image or the vocabulary size after removing stop-words, is in the tens of thousands. An offline computation of the principal components would require the computation of eigenvectors of the sample covariance matrix. However, in high-dimensional settings, storing the covariance matrix and subsequent eigen-analysis can be challenging. Streaming PCA methods have gained significant traction owing to their ability to iteratively update the principal components by considering one data-point at a time.
One of the most widely used algorithms for streaming PCA is Oja’s algorithm, proposed in the seminal work of [41]. Oja’s algorithm involves the following update rule:
wt+1 − wt = η(wTt Xt)Xt; wTt+1wt+1 = 1, (1)
where Xt ∈ Rd is the tth data point and wt is the current estimate for the leading eigenvector of Σ = EXXT after t data-points have been seen. The parameter η can be thought of as a learning rate, which can either be fixed or varied as a function of t. In this paper we fix the learning rate, similar to [26].
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
Contribution: In the present work, we consider the problem of uncertainty quantification for the estimation error of the leading eigenvector from Oja’s algorithm, which is one of the most commonly used streaming PCA algorithms. Our contributions may be summarized as follows:
1. We derive a high-dimensional weighted χ2 approximation to the sin2 error for the leading eigenvector of Oja’s algorithm. We recover the optimal convergence rate O(1/n) while allowing d to grow at a sub-exponential rate under suitable structural assumptions on the covariance matrix, matching state-of-the-art theoretical results for consistency of Oja’s algorithm. Our result provides a distributional characterization of the sin2 error for Oja’s algorithm for the first time in the literature. The approximation holds for a wide range of step sizes.
2. Since the weighted χ2 approximation depends on unknown parameters, we propose an online bootstrap algorithm and establish conditions under which the bootstrap is consistent. Our bootstrap procedure allows the approximation of important quantities such as the quantiles of the error associated with Oja’s algorithm for the first time.
Prior analysis of Oja’s algorithm. While Oja’s algorithm was invented in 1982 it was not until recently that the theoretical workings of Oja’s algorithm have been understood. A number of papers in recent years have focused on proving guarantees of convergence of the iterative update in (1) toward the principal eigenvector of the (unknown) covariance matrix EXXT , which can be recast as stochastic gradient descent (SGD) on the quadratic objective function
min w
wTw=1
−trace(wTΣw), Σ = EXXT , (2)
projected onto the non-convex unit sphere. We assume that the data-points are mean zero. Despite being non-convex and thus falling outside the framework for which theory for stochastic gradient descent convergence is firmly established, the output of Oja’s algorithm be viewed as a product of random matrices and shares similar structure to other important classes of non-convex problems, such as matrix completion [27, 29], matrix sensing [27], and subspace tracking [4]. Thus, studying this optimization problem serves as a natural first step toward understanding the behavior of SGD in more general non-convex settings.
Let v1 denote the principal eigenvector of Σ, and let v̂1 = wn be the solution to the stochastic iterative method applying Eq 1. Finally, let λ1 > λ2 be the first and second principal eigenvalues of Σ. Sharp rates of convergence for Oja’s updates were established in [25]. Under boundedness assumptions on ‖XiXTi − Σ‖, they show that with constant probability, the square of the sine of the angle between v1 and w satisfies:
1− (vT1 v̂1)2 = O ( 1
n
) (3)
where the O hides a constant which depends in the optimal way on the eigengap between the top two eigenvalues, and independent of n or d, improving on previous error bounds for Oja’s algorithm [46, 18, 3, 47, 38, 2] which showed convergence rates that deteriorate with the ambient dimension d, and thus did not fully explain the efficiency of Oja’s update. This sharp rate is remarkable, as it matches the error of the principal eigenvector of the sample covariance matrix, which is the batch or offline version of PCA. Other notable work include [31, 33] for unbounded Xi, analysis of Oja’s algorithm for computing top k principal components [1, 24].
The bootstrap. The bootstrap, proposed by [9], is one of the most widely used methods for uncertainty quantification in machine learning and statistics and accordingly has a vast literature. We refer the reader to [17, 49] for expositions on the classical theory of the bootstrap for IID data. Recently, since the groundbreaking work of [7, 8], the bootstrap has seen a renewed surge of interest in the context of high-dimensional data where d can be potentially exponentially larger than n. Of particular relevance to the present work are high-dimensional central limit theorems (CLTs) for quadratic forms, which have been studied by [43, 51, 15]. In particular, our CLT for the estimation error of Oja’s algorithm invokes a modest adaptation of [51] to independent but non-identically distributed random variables. In machine learning, bootstrap methods have been used to estimate the uncertainty of randomized algorithms such as bagging and random forests [35], sketching for large scale singular value decomposition (SVD) [36], randomized matrix multiplication [37], and randomized least squares [34].
A standard notion of bootstrap consistency is that, conditioned on the data, the distribution of the suitably centered and scaled bootstrap functional approaches the true distribution with high probability in some norm on probability measures, typically the Kolmogorov distance, which is the supremum of the absolute pointwise difference between two CDFs. Bootstrap consistency is often established by deriving a Gaussian approximation for the sampling distribution and showing that the bootstrap distribution is close to the corresponding Gaussian approximation with high probability.
It may seem that if one knows that the approximating distribution of a statistic is Gaussian, this defeats the purpose of bootstrap. However, for most statistics, the parameters of the normal approximation depend on unknown model parameters, and have to be estimated if one intends to use the normal approximation. Furthermore, the CLT only gives a first-order correct approximation of the target distribution, i.e. with O(1/ √ n) error. In contrast, the bootstrap of a suitably centered and scaled statistic has been shown to be higher order correct for many functionals [16, 17, 19].
Quantifying uncertainty for SGD. Behind the recent success of neural networks in a wide range of sub-fields of machine learning, the workhorse algorithm has become Stochastic gradient descent (SGD) [42, 40, 44]. For establishing consistency of bootstrap, one requires to establish asymptotic normality [11, 42, 45, 39]. There has also been many works on uncertainty estimation of SGD [6, 32, 12, 48]. However, all these works are for convex, and predominantly strongly convex loss functions. Only recently, [52] has established asymptotic normality for nonconvex loss functions under dissipativity conditions and appropriate growth conditions on the gradient, which are weaker conditions than strong convexity but not significantly so.
Now, in Section 2 we present notation and do setup, present our main theoretical results in Section 3, followed by simulations in Section 4.

heading: 2 Preliminaries
text: We consider a row-wise IID triangular array, where the random vectors {Xi} in the nth row take values in Rdn , with E[Xi] = 0 and Var(Xi) = Σn. Note that the triangular array allows {X1, . . . , Xn} come from a different distribution for each n and the setting where d is fixed and n grows is a special case. For readability, we drop the subscript n from Σn. We use ‖ · ‖ to denote the Euclidean norm for vectors and the operator norm for matrices and ‖ · ‖F to denote the Frobenius norm. Expanding out the recursive definition in Eq 1, we see that Oja’s iteration can be expressed as wt+1 = (Id + ηXtX T t )wt. Thus, after n iterations the vector can be written as a matrix-vector product, where the matrix is a product of n independent matrices. Expanding out the recursive definition, we get:
Bn := n∏ i=1 (Id + ηXiX T i ) v̂1 = Bnu0 ‖Bnu0‖ , (4)
where Id is a d× d identity matrix. where u0 is a random unit vector in d dimensions. In the scalar case, when η = 1/n, for large n, the numerator of Eq 4 behaves like exp( ∑ iX 2 i /n), which in turn converges to exp(E[X21 ]). For matrices, one hopes that, by independence, a result of the same flavor will hold. And in fact if it does hold, then for η = lognn , the numerator in Eq 4 will concentrate around exp(log nΣ). The spectrum of this matrix is dominated by the principal eigenvector, i.e. the ratio of the first eigenvalue to the second one is exp(log n(λ1 − λ2)), where λi is the ith eigenvalue of the covariance matrix Σ. This makes it clear that Oja’s algorithm is essentially a matrix vector product of this matrix exponential (suitably scaled) and a random unit vector.
However, the intuition from the scalar case is nontrivial to generalize to matrices due to noncommutativity. Limits of products of random matrices have been studied in mathematics in the context of ergodic theory on Markov chains (see [14, 30, 5, 10] etc.). However, until recent results of [23], which extended and improved results in [21], there has not been much work on quantifying the exact rate of convergence, or finite-sample large deviation bounds for how a random matrix product deviates from its expectation.
We reparametrize η as ηn/n, where ηn is chosen carefully to obtain a suitable error rate. Note that this is not a scheme where we decrease η over time as in [20], but hold it as a constant which is a function of the total number of data-points.

heading: 2.1 The Hoeffding decomposition
text: The Hoeffding decomposition, attributed to [22], is a key technical tool for studying the asymptotic properties of U-statistics. However, the idea generalizes far beyond U-statistics; see Supplement Section A for further discussion. In the present work, we use Hoeffding decompositions for matrix and vector-valued functions of independent random variables taking values in Rd to facilitate analysis for Bn.
A concept closely related to the Hoeffding decomposition is the more well-known Hájek projection, which gives the best approximation (in an L2 sense) of a general function of n independent random variables by a function of the form ∑ i gi(Xi), where gi are measurable functions satisfying a square integrability condition. The Hájek projection facilitates distributional approximations for complicated statistics since this linear projection is typically more amenable to analysis. However, establishing a central limit theorem requires showing the negligibility of a remainder term, which can be large if the projection is not accurate enough.
The Hájek projection may be viewed as the first-order term in the Hoeffding decomposition, a general way of representing functions of independent random variables. The Hoeffding decomposition consists of a sum of projections onto a linear space, quadratic space, cubic space, and so on. Each new space is chosen to be orthogonal to the previous space. Thus, the Hoeffding decomposition can be thought of as a sum of terms of increasing levels of complexity. Even if the remainder of the Hájek projection turns out to be small, the Hoeffding decomposition can be easier to work with due to the orthogonality of the projections.
The Hoeffding decomposition for the matrix product. Let Yi = XiXTi −Σ and let S ⊆ {1, . . . n}. By Corollary A.1 of the Supplement Section A, the Hoeffding Decomposition for Bn is given by:
Bn = n∑ k=0 Tk, Tk = ∑ |S|=k H(S). (5)
where H(S) = ∏n i=1A (S) i and A (S) i is given by: A (S) i = {ηn n Yi if i ∈ S I + ηnn Σ otherwise .
The above expansion has favorable properties that facilitate second-moment calculations. In fact, as a consequence of the orthogonality property of Hoeffding projections, we have that
E [ ‖Bn‖2F ] = n∑ k=0 ∑ |S|=k E [ ‖ n∏ i=1 H (S) i ‖ 2 F ]
E [ ‖Bnx‖2 ] = n∑ k=0 ∑ |S|=k E
[ ‖
n∏ i=1 H (S) i x‖ 2 ] where the second statement holds for any x ∈ Rd; see Proposition A.2 in Supplement Section A.

heading: 2.2 Online bootstrap for streaming PCA
text: To approximate the sampling distribution, we consider a Gaussian multiplier bootstrap procedure. As observed by [7], a Gaussian multiplier random variable eliminates the need to establish a Gaussian approximation for the bootstrap since conditional on the data, it is already Gaussian. It is not hard to see that this is a natural candidate for the online setting; the multiplier bootstrap has been used for bootstrapping the stochastic gradient descent estimator in [12].
We present our bootstrap in Algorithm 1. In our procedure, we update m+1 vectors at every iteration. The first one is v̂, which will result in the final Oja estimate of the first principal component. The other vectors {v∗(j), j = 1, . . .m} are obtained by perturbing the basic Oja update (Eq 1). The Wi’s are the multiplier random variables, which are scaled mean zero scaled Gaussians with variance 1/2. The update of the v∗(j) is novel because it preserves the mean and the variance of the original Oja estimator while not requiring access to the full sample covariance matrix. Consequently, we can make our updates online and attain both a point estimate and a confidence interval for the principal eigenvector, while increasing the computation and storage by only a factor of m.
Algorithm 1: Bootstrap for Oja’s algorithm Input: Datapoints X1, . . . , Xn, stepsize η, number of bootstrap replicates m Output: Oja’s solution v̂1 and m bootstrapped versions of it v ∗(1) 1 , . . . , v ∗(m) 1 Draw g ∼ N(0, Id) Create unit vector u0 ← g/‖g‖ Initialize v̂1, v ∗(1) 1 , . . . , v ∗(m) 1 ← u0 for t=2,. . . , n do Update v̂1 ← v̂1 + η(XTt v̂1)v̂1 Normalize v̂1 to have unit norm; for i=1:m do
Draw Wi ∼ N(0, 1/2); Let h(i) ← (XTt v ∗(i) 1 )Xt; Let g(i) ← (XTt−1v ∗(i) 1 )Xt−1; Update v∗(i)1 ← v ∗(i) 1 + η ( h(i) +Wi(h (i) − g(i)) ) ;
Normalize v∗(i)1 to have unit norm; end
end

heading: 3 Main results
text: In this section we present our main contributions: a CLT for the error of Oja’s algorithm and consistency of an online multiplier bootstrap for error.

heading: 3.1 Central limit theorem for the error of Oja’s algorithm
text: We start by stating a CLT for the error of Oja’s algorithm. To state this theorem, we will need to introduce some notation.
Let v̂1 denote the Oja vector and V⊥ the d×d−1 matrix with 2, . . . , d eigenvectors of Σ on its columns. Note that V⊥ is not uniquely defined, but V⊥V T⊥ = I − v1vT1 is if the leading eigenvalue is distinct and consequently, norms of the form ‖V T⊥ x‖ for x ∈ Rd are well-defined. Let λ1 ≥ · · · ≥ λd denote the eigenvalues of Σ and Λ⊥ be a diagonal matrix with Λ⊥(i, i) = (1 + ηnλi+1/n)/(1 + ηnλ1/n), i = 1, . . . , d− 1. Also let
M := E [ V T⊥ (X T 1 v1) 2X1X T 1 V⊥ ] (6)
Now we define
V̄n = ηn n ∑ i E[V⊥Λi−1⊥ V T ⊥ (XiX T i − Σ)v1vT1 (XiXTi − Σ)V⊥Λi−1⊥ V T ⊥ ]
= ηn n V⊥ (∑ i Λi−1⊥ MΛ i−1 ⊥ ) V T⊥ (7)
We have the following result: Theorem 1. Suppose that u0 is drawn from the uniform distribution on Sd−1, λ1 = O(1). Choose ηn → ∞ such that nd · exp(−ηn(λ1 − λ2)) → 0, (ηn∨ log d) η 2 n(M 2 d∨1)
n → 0, where Md = E[ ∥∥XiXTi − Σ∥∥2]. Further, let Z̃n be a mean 0 Gaussian matrix such that Var(Z̃n) = Var((X1XT1 − Σ)v1) and suppose that:
‖M‖F ≥ c > 0 (8)
E [∥∥V T⊥ Z̃n∥∥6] ∨ E [∥∥V T⊥ (X1XT1 − Σ)v1∥∥6]
‖M‖3F = o(n) (9)
Then, for a sequence of Gaussian distributions {Zn}n≥1 with mean 0 and covariance matrix V̄n (see Eq 7), the following holds:
sup t∈R ∣∣P (n/ηn sin2(v̂1, v1) ≤ t)− P (ZTn Zn ≤ t)∣∣→ 0 (10) Theorem 1 is very general. We allow the dimension to grow with the number of observations, which is typical in the high-dimensional bootstrap literature. Note that the case of fixed d and growing n is also a special case of this setup.
We want to point out that while previous literature obtained sharp bounds on the sin2 error 1−(vT1 v̂1)2, we go a step further. We establish an approximating distribution for n/ηn(1− (vT1 v̂1)2). Remark 1 (Condition on norm). For simplicity, we assume λ1 = O(1), which can be easily relaxed to grow slowly with n. We do not assume that the ‖XiXTi − Σ‖2 is bounded almost surely. However, the norm of XiXTi − Σ comes into play implicitly via the assumption in Eq 9. Consider the case where Xi are drawn from some multivariate Gaussian distribution. We use this to build intuition about the assumptions in Eq 8 and 9. In this case, XT1 V⊥ is a Gaussian of
independent entries and thus E [∥∥V⊥(X1XT1 − Σ)v1∥∥6] = E‖XT1 v1‖6E(∑j>1(XTj vj)2)3. Note
that ∑ j>1((X T j vj) 2 − λj) is a sub-exponential random variable with parameters (c1 ∑ j>1 λj , c2).
Furthermore, ‖M‖2F = λ1 ∑ i>1 λi. Thus Eq 9 reduces to checking if
λ 3/2 1 ( ∑ j>1 λj) 3
( ∑ i λ 2 i ) 3/2 = o(n)
Remark 2 (Coordinates with summable sub-Gaussian parameters). Eq 9 imposes a growth condition on the moments of both the data and a Gaussian analog. One setting for which both growth rates are in fact bounded is if the coordinates of X are sub-Gaussian and the sub-Gaussian parameters satisfy∑d i=1 νi < C <∞ following similar arguments to Proposition 1. Remark 3 (Constant vs Adaptive Learning Rate). Adaptive learning rates are also commonly studied in the literature on Oja’s algorithm and have the advantage that they require no prior knowledge of the sample size. It should be noted that our results hold for a wide range of learning rates, ranging from log(nd) ηn n1/3, so our results will still apply so long as in the initial guess of the sample size is not off by orders of magnitude. We leave a detailed study of the adaptive learning rate setting to future work.
As a corollary of our main theorem, we obtain the following error bound on the sin2 error. Corollary 1. Under the conditions in Theorem 1, we have
sin2(v̂1, v1) = OP
( ηnMd
n(λ1 − λ2) ) Remark 4 (Comparison with previous work). As a byproduct of our analysis, we recover the sharpest convergence rates for Oja’s algorithm in the literature. If we set ηn = c1 log nd/(λ1 − λ2), for large
enough c1, the dominating term in the error is OP
( Md log nd
n(λ1 − λ2)2
) under mild conditions on d. This
matches the bound in [25]. Remark 5 (Rate of convergence in Kolmogorov distance). To simplify the theorem statement, we have stated Theorem 1 without giving an explicit rate of convergence in the Kolmogorov distance. Convergence rates depend on the rate of decay of the remainder terms, which are worked out in Supplement Section B.3, and the magnitude of the quantity in Eq 9. The contribution of the latter quantity to the rate is worked out in the IID case in [51]. Remark 6 (Lower bound on norm). While our rate matches the sharp bounds in literature and our assumptions on norm upper bounds are similar or weaker than previous work, we do assume a lower bound on the Frobenius norm of the covariance matrix as in Eq 8. Note that if indeed all Xi’s were a scalar multiple of v1, then the V̄n matrix in Eq 7 will be zero. This will lead to a perfect point estimate, but there will not be any variability from the data and hence there will be no non-degenerate approximation. The lower bound on the norm is not resulting from loose analysis. Similar lower bounds on the variance are imposed in the high-dimensional CLT literature [7, 8].
Now we provide a proof sketch of Theorem 1 below.
Proof sketch for Theorem 1. We provide the main steps in our derivation. The detailed calculations can be found in Supplement Section B.
1. We start by expressing the sin2 error as a quadratic form:
sin2(v1, v̂1) = 1− uT0 B T n v1v T 1 Bnu0
uT0 B T nBnu0
= uT0 B T n (I − v1vT1 )Bnu0 uT0 B T nBnu0
= (V⊥V
T ⊥Bnu0) T (V⊥V T ⊥Bnu0)
‖Bnu0‖2
(11)
where in the last line we used the fact that V⊥V T⊥ is idempotent. Our proof strategy for the central limit theorem involves further approximating Eq 11 with an inner product of the Hájek projection (first-order) term in Eq 5.
2. Our second step is to show that ‖Bnu0‖ concentrates around its expectation (1+ηnλ1/n)n|vT1 u0|.
3. Next we establish that ‖V⊥V⊥BnV⊥V T ⊥ u0‖2
‖Bnu0‖ is OP
(√ d · exp{−ηn(λ1 − λ2)}+ √ η3nM 2 d log d
n2
) .
This is achieved by using a similar recursive argument as in [25], but with the crucial observation that the residual or common difference term is of a lower order because it can be replaced by a matrix product minus its expectation.
4. Now we go back to the expansion in Eq 5.
(vT1 u0)V⊥V T ⊥Bnv1 = (v T 1 u0) ∑ k V⊥V T ⊥ Tkv1
Since T0 = (I + ηn/nΣ)n, V⊥V T⊥ T0v1 is the zero vector. Now we examine the (vT1 u0)V⊥V T ⊥ (Bn − T1)v1 term. Here we use the structure of the higher order terms Tk. In particular, we use the fact that it is a matrix product interlaced with k XiXTi − Σ matrices. For example, for k = 2 we have
T2 = η2n n2 ∑ i<j ( I + ηn n Σ )i−1 Yi ( I + ηn n Σ )j−i−1 Yj ( I + ηn n Σ )n−j
We show that the norm of (vT1 u0)V⊥V T ⊥ (Bn − T1)v1, normalized by the denominator, is O(η2nM 2 d/n
2). The fact that the summands of Tk are uncorrelated and Tk and T` are uncorrelated for k 6= ` makes this possible.
5. Finally, we are left with V⊥V T⊥ T1v1(v T 1 u0). Note that this is of the following form:
ηn n
(vT1 u0)V⊥V T ⊥ T1v1 |vT1 u0|(1 + λ1ηn/n)n = ηnsgn(vT1 u0) n n∑ i=1 V⊥Λ i−1 ⊥ V T ⊥ (XiX T i − Σ)v1
It is not hard to see that this is a sum of independent random vectors with covariance matrix ηn/nV̄n (see Eq 7).
6. We adapt a result of distributional convergence of squared norm of sums of IID random vectors in [51] to squared norm of sums of independent random vectors. Under the assumptions 9 and 8, the conditions of distributional convergence are satisfied.
7. Finally, all the error terms are combined along with an anti-concentration argument for χ2 to establish the final result. The full proof and accompanying lemmas are in Section B of the Supplement.

heading: 3.2 Bootstrap consistency
text: Using the weighted χ2 approximation for inference requires estimating the eigenvalues of Σ and other population quantities; however, accurate estimates may not be available in a streaming setting. Instead, we propose a streaming bootstrap procedure that mimics the properties of the original Oja algorithm. While a similar structure leads to error terms that are similar to the CLT, the analysis of the bootstrap presents its own technical challenges. In what follows let P ∗ denote the bootstrap measure, which is conditioned on the data, and let E∗[·] denote the corresponding expectation operator. A common strategy for establishing consistency of the Gaussian multiplier bootstrap is to invoke a Gaussian comparison lemma. Since the multipliers are themselves Gaussian and the data is treated as fixed, the idea is that one can use specialized results for comparing the distributions of two Gaussians (bootstrapped Z∗n and approximating Zn from the CLT) that only depend on how close the covariance matrices E∗[Z∗nZ∗Tn ] and E[ZnZTn ] are in an appropriate metric. Using a Gaussian comparison lemma for quadratic forms (see Supplement Section C.3), we have the following result for the bootstrapped sin2 error:
Lemma 1. [Bounding the difference between the bootstrap covariance and true covariance] Let:
Z∗n = sgn(v T 1 u0) √ ηn n ∑ i WiV⊥Λ i−1 ⊥ V T ⊥ (XiX T i −Xi−1XTi−1)v1. (12)
Recall the definition of V̄n from Eq 7. We have,
|trace(E∗[Z∗nZ∗Tn ]− V̄n)|, ‖E∗[Z∗nZ∗Tn ]− V̄n‖F = OP √E‖X1XT1 − Σ‖4 n(λ1 − λ2)  With this lemma in hand, we are ready to state our bootstrap result.
Theorem 2 (Bootstrap Consistency). Suppose that the conditions of Theorem 1 are satisfied. Furthermore, let αn be a sequence such that P (Acn) → 0, where An is defined as An ={
maxi≤i≤n ‖Xi‖2 ≤ αn } . Further suppose that Md log 2 d η2n n → 0, (α3n∨Md log d) αnη 3 n
n → 0, αnMdη 2 n
n(λ1−λ2) → 0, and E[‖X1XT1 −Σ‖4] n(λ1−λ2) → 0. Then,
sup t∈R ∣∣P ∗(n/ηn sin2(v∗1 , v̂1) ≤ t)− P (n/ηn sin2(v̂1, v1) ≤ t)∣∣ P−→ 0 Proof sketch of Theorem 2. The proof follows a similar route to Theorem 2. We provide a detailed analysis in Supplementary Section. We use a bootstrap version of the Hoeffding decomposition conditioned on the data, stated in Supplement Section. In step one we have B∗n replace Bn, where B∗n is given by:
B∗n = n∏ i=1 ( I + ηn/n(XiX T i +Wi(XiX T i −Xi−1XTi−1) ) We work out Step 1 using concentration of matrix products [23]. For steps 2-3, we see that T ∗k has the same structure as Tk with the difference that (I + ηnΣ/n)i is replaced by its sample counterpart which is a product of i independent matrices of the form I + ηn/nXjXTj . Concentration of these terms in operator norm are established with results from [23]. Finally for step 4, we see that the main term that approximates the bootstrap residual V̂⊥V̂ T⊥B ∗ nu0 is given by √ ηn/nZ ∗ n, where Z ∗ n is given in Eq 12. Conditioned on the data, this is already Normally distributed since the multiplier random variables Wi are themselves Gaussian. We then invoke the Gaussian comparison result Lemma 1 to obtain convergence to the weighted χ2 approximation.
We now make a couple of points regarding our analysis. It should be noted that the terms in the product are weakly dependent, which is different from the CLT and would seem to complicate concentration arguments used to establish bootstrap consistency. However, the dependence is not strong and second-moment methods may be used. We also operate on a good set in which the norms
of the the updates are not too large, which is far less restrictive than assuming an almost sure bound on the norm.
In theorem above, we have stated the good set An in an abstract manner, but one may wonder how stringent the condition is in various problem settings. Below, we describe a general setup with sub-Gaussian entries of Xi in which αn grows as log n; under milder forms of various decay, all we need is for αn to grow slowly with n. Here ‖·‖ψ1 is the sub-Exponential Orlicz norm and ‖·‖ψ2 is the sub-Gaussian Orlicz norm (see, for example [50]).
Proposition 1 (The effect of variance decay on the norm). For each 1 ≤ j ≤ p, suppose that X1j satisfies ‖X1j‖ψ2 ≤ νj ∑p j=1 νj ≤ C1 < ∞. Then, for some universal constant C2 > 0,∥∥∥∑pj=1(X21j − E[X21j ])∥∥∥
ψ1 < C2, and for some c1, c2 > 0,
P ( max
1≤i≤n ‖Xi‖2 > c1 log n
) ≤ c2
n
We now present experimental validation of our bootstrap procedure below.

heading: 4 Experimental validation of the online multiplier bootstrap
text: We draw Zij IID∼ Uniform(− √ 3, √
3), for i = 1, . . . , n and j = 1, . . . d. Consider a PSD matrix Kij = exp(−|i− j|c) with c = 0.01. We create a covariance matrix such that Σij = K(i, j)σiσj . We consider σi = 5i−β for β = 0.2 and β = 1. Now we transform the data to introduce dependence by letting Xi = Σ1/2Zi. By construction, we have that E[XiXTi ] = Σ for all 1 ≤ i ≤ n. Our goal is to simply demonstrate that the bootstrap distribution of sin2 errors closely match that of the sampling distribution. To this effect, we fix u0 and draw 500 datasets and run streaming PCA on each and then construct an empirical CDF (F ) from the sin2 error with the true v1. This is the point of comparison for the bootstrap distribution (F ∗), for which we fix a datasetX . We then invoke algorithm 1 to obtain 500 bootstrap replicates v̂∗1 as well as the Oja vector for the dataset v̂1. The bootstrap distribution is the empirical CDF of 1 − (v̂T1 v̂∗1)2. We use ηn = log n. In Figure 1, we see that for β = 0.2 (see (A) and (B)), where the variance decay is slow and therefore the error bounds of the residual terms are expected to be large, the quality of approximation is poorer compared to (C) and (D), where β = 1. However, even for β = 0.2, increasing n improves performance. Also note that, for (A) and (B) the variance decay does not satisfy our theorems conditions and thus, the normalized error does not behave like a OP (1) random variable. However, for (C) and (D) the variance decay satisfies the conditions and in this case the normalized error is OP (1), which happens to be in the [0,1] range for this example.

heading: 5 Discussion
text: Modern tools in non-asymptotic random matrix theory have given rise to recent breakthroughs in establishing pointwise convergence rates for stochastic iterative methods in optimizing certain nonconvex objectives, including the classic Oja’s algorithm for online principal component analysis. By synthesizing modern random matrix theory tools with classic results from the U-statistics literature and recently developed high-dimensional central limit theorems, we extend the error analysis of Oja’s algorithm from pointwise convergence rates to distributional convergence and moreover establish an efficient online bootstrap method for Oja’s algorithm to quantify the error on the fly. Our results are a first step toward incorporating uncertainty estimation into the general framework of stochastic optimization algorithms, but we acknowledge the present limitations of our analysis: new tools will be needed to extend the current analysis to estimating higher-dimensional principal subspaces, and additional tools will be needed to account for non-independent matrix products which appear beyond the setting of online PCA.

heading: Acknowledgment
text: P.S. and R.L. are supported in part by NSF 2019844 and NSF HDR-1934932. R.W. is supported in part by AFOSR MURI FA9550-19-1-0005, NSF DMS 1952735, NSF HDR-1934932, and NSF 2019844.

Review summaries:
This paper develops an online bootstrap method to quantify the error of estimating the leading eigenvector of a matrix via Oja's algorithm. The authors a) find a closed from expression of the leading eigenvector from Oja's algorithm, b) approximate this closed from expression via Hoeffding decomposition for U-statistics, and c) devise a Gaussian multiplier bootstrap (algorithm) to approximate the distribution of the Hoeffding decomposition. The authors establish Gaussian approximation and bootstrap consistency results and include a small scale simulation study.

This paper study the problem of uncertainty quantification for the estimation error of the leading eigenvector from Oja's algorithm. The authors establish the distribution of $\sin^2$ error of the leading eigenvetor of Oja's algorithm. Also, the authors propose an online bootstrap algorithm for Oja's algorithm and provide the conditions under which the bootstrap is consistent.  

Oja’s principal subspace algorithm is a well-known and powerful technique for learning and tracking principal information in a time series setting. In this paper, the authors use U-statistics and a concentration of matrix products to obtain an error estimate between the leading eigenvector and the output of Oja’s algorithm. The other eigenvectors are probably beyond the theoretical understanding of random matrix theory at the moment. There are two main contributions: (1) Establishes the distribution of the sin^2 error estimate for the leading eigenvector computed via Oja’s algorithm and (2) A stream-friendly bootstrap algorithm to compute quantiles of the error estimate. The paper contains an interesting theoretical analysis involving modern matrix random theory, and since they have the distribution, a potential application to bootstrapping. The main novelties in this paper are theoretical, and the author could better demonstrate the advantages of their online multiplier bootstrapping technique.