Bootstrapping the Error of Oja's Algorithm 
Authors: Robert Lunde
Purnamrita Sarkar
Rachel Ward
Venue: NeurIPS 2021 Spotlight (2021)
Keywords: Gaussian approximation
bootstrap
Streaming PCA

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 $\chi^2$ approximation result for the $\sin^2$ error between the population eigenvector and the output of Oja’s 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. 

Bootstrapping the error of Oja's algorithm
Robert Lunde University of Michigan rlunde@umico cemniih.eee edu
Purnamrita Sarkar
University of Texas at Austin purna. sarkar@austin utexas.edi
Rachel Ward
University of Texas at Austin rward@math. utexas.
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-statisties literature with recent results on high-dimensional central limit theorems for quadratic forms of random vectors and concentration matrix products, we establish a weighted x² approximation result for the sin error between the population eigenvector and the ouiput of Ojas algorithm. Since estimating the covariance matrix associated with the approximat- propose
ing distribution requires knowledge of unknown model parameters, multiplier bootstrap algorithm that may be updated in an online manner. We estab- lish 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.
1
Introduction
Since its discovery over a century ago 13] principal component analysis (PCA) has been a corne-- stone data analysis. many applications, dimension reduction paramount and PCA offers an optimal low-rank approximation of the original data. PCA is also highly interpretable as it projects 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 compress a large dimensional dataset a manageable size without incurring much loss; for 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 siigh-iimnsionall settings, storing the covariance matrix and subsequent 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 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:
W++1 wt nn Xx)
T W1+1 = 1,
(I)
where X, E Rd is the th data point and we is the current estimate for the leading eigenvector of E = E after t dota-points have been seen. The parameter 1) can be thought of as learning rate, which can either be fixed or varied as a function of In this paper we fix the learning rate, similar to 266
35th Conference on Neural Information Processing Systems (NcurIPS 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 high-dimensional weighted x² approximation to the sin" error for the leading eigenvector of Oja's algorithm. We recover the optimal convergence rate O(1/n) while allowing d to grow at 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 distributional characterization of the sin' 2 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 x" 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 eтror associated with Oja's 's algorithm for the first time.
Prior analysis of Oja's algorithm. While Oja's algorithm was invented in 1982 was not until recently that the theoretical workings of Oja's 's algorithm have been understood. A number of papers in recent years have focused on proving guarantees of convergence of the iterative update (I) toward the principal eigenvector of the (unknown) covariance matrix EX which can be recast as stochastic gradient descent (SGD) on the quadratic objective function
min traco(w" w Ew), T w==1
Σ EXXX
projected onto the non-convex unit sphere. We assume that the data- mata-pooins 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 is algorithm be viewed
product of random matrices and shares similar structure to other important classes of non-convex problems, such as matrix completion [27 22], 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 d1 denote the principal eigenvector of E, and let 01 = w,, be the solution to the stochastic iterative method applying Eq[T] Finally, let X > X2 be the first and second principal eigenvalues of Σ. Sharp rates of convergence for Oja's 's updates were established in (25). Under boundedness assumptions on ||x;X?' Σ| they show that with constant probability, the square of the sine of the angle between 01 and w satisfies:
0 1
1 - (of 6,
where the O hides a constant which depends in the optimal way on the eigengap between the top two eigenvalues, and independent of 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 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 [3 33) for unbounded X,, analysis of Oja's
algorithm for computing top k principal components 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 [17 44 for expositions on the classical theory of the bootstrap for IID data. Recently, since the groundbreaking work of [7] [). 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 's algorithm invokes modest adaptation of [5I] to independent but non-identically distributed random variables. 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 [4].

---

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 appoximattio. Furthermore, the CLT only gives a first-order correct approximation of the target distribution, i.e. with OLL/ error. In contrast, the bootstrap of suitably centered and scaled statistic has been shown to be higher order correct for many functionals [6 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 44). For establishing consistency of bootstrap, one requires to establish asymptotic
(SGD) [4 40
45 39]. There has also been many works on uncertainty estimation of SGD (6] 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
normality
32
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
2 Preliminaries
We consider row-wise IID triangular array, where the random vectors {X,} in the nith row take val- ues in Rd,, with [[X,]= =x]]=00 0 and Var(X,) = Σn. Note that the triangular array allows ( X come from different distribution for each n and the setting where is fixed and grows is special case. For readability, we drop the subscript n from E,, We use || || to denote the Euclidean norm for vectors and the operator norm for matrices and -| || to denote the Frobenius norm.
Expanding out the recursive definition in Eq [ we see that Oja's iteration can be expressed as Wr+1 (la + xx7 )u,, 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:
B u0
B, = II ( (1++
01 =
(4)
|,,uvll
where la is d x d identity matrix. where uo is random unit vector in d dimensions. In the scalar case, when 17 = 1/n,, for large n. the numerator of Eq 4 behaves like exp(E; X²/n), which in turn converges to exp (xx))( For matrices, one hopes that, by independence, result of the same flavor will hold. And in fact if does hold, then for 7 Jug the numerator in Eq 4 will concentrate around cxp(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(Ai - ->2)). where X is the jth eigenvalue of the covariance matrix Σ. This makes it clear that Oja's 's algorithm is essentially a matrix vector product of this matrix exponential (suitably scaled) and random unit vector.
However, the intuition from the scalar case is nontrivial to generalize to matrices due to non- commutativity. Limits of products of random matrices have been studied in mathematics in the context of ergodic theory Markov chains (see 14 30 5 I0 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-saple large deviation bounds for how random matrix product deviates from its expectation.
We reparametrize 1) as nn /n, where Th is chosen carefully to obtain a suitable enror rate. Note that this is not a scheme where we decrease 1) over time as in 20 but hold as a constant which is a function of the total number of data-points.
3

---

2.1 The Hoeffding decomposition
The Hoeffding decomposition, attributed to [22] is 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-valuee functions of independent random variables taking values in R to facilitate analysis for B,,
concept closely related to the Hoeffding decomposition is the more well-known Hájek projection, which gives the best approximation (in L2 sense) of general function of independent random variables by a function of the form >; 9: (X;), where gi are measurable functions satisfying square integrability condition. The Há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 Hájek projection may be viewed as the first-order term in the Hoeffding decomposition, 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 Each new space is chosen to be orthogonal to the previous space. Thus, the Hoeffding decomposition can be thought of a sum of terms of increasing levels of complexity. Even if the remainder of the Há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 Y; X;X and let S (L....n.. By Corollary of the Supplement Section the Hoeffding Decomposition for B, is given by:
n...
A
Bn = Tk
Tk
H())
(5)
k=0
|
where (()) =
( and 1 is given by: 1
Um Y
i
Σ
otherwise
The above expansion has favorable properties that facilitate second- consequence of the orthogonality property of Hoeffding projections, we have that
second-moment moment calculations. In fact, as a
B
k=0
|S|
=
i=1
E
||
Σ
H (3) x||
where the second statement holds for any E Rdd see Proposition A.2 in Supplement Section A
2.2 Online bootstrap for streaming PCA
To approximate the sampling distribution, we consider a Gaussian multiplier bootstrap procedure. As observed by 77). Gaussian multiplier random variable eliminates the need to establish a Gaussian approximation for the bootstrap since conditional on the data, already Gaussian. It is not hard to see that this is natural candidate for the online setting; the multiplier bootstrap has been used for bootstrapping the stochastic gradient descent estimator in 112.
We present our bootstrap in Algorithm|1 [ In our procedure, we update m++ 1 vectors at every iteration. The first one is é, which will result the final Oja estimate of the first principal component. The other vectors {v* .j = ..... m) are obtained by perturbing the basic Oja update (Eq Ib.
The Wi's are the multiplier random variables, which are scaled mean zero scaled Gaussians with variance 1/2. The update of the -²( 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 factor of m.

---

Algorithm 1: Bootstrap for Oja's algorithm
Input:
Datapoints X1, Xi--.., Xns stepsize 7/, number of bootstrap replicates m
Output: Oja's solution 81 and m bootstrapped versions vi (11 Draw g (O,
(mm)
0j
Create unit vector uo + 9//llgl
Initialize 0,, 01 ((1) ,..., ((m) ug
for 1= R1 do
Update 81 < 01 (x? 61)
Normalize in to have unit norm;
for i= do
Draw W; N(0,
Let h(() +- (X7 vi (1) X,:
Let g(t) < (X7 U1 Xx-11
Update "i .=()) t U1 (1)
+ (h + Wi(h() g(o)):
(())
Normalize vi to have unit norm;
end end
Main results
In this section present our main contributions: CLT for the error of Oja's algorithm and consistency of an online multiplier bootstrap for error.
3.1 Central limit theorem for the error of Oja's algorithm
We start by stating CLT for the error of Oja's algorithm. To state this theorem, we will need to introduce some notation.
Let úJ denote the Oja vector and V1 the dl 1xd11- matrix with 2..... . eigenvectors of Σ on its columns. Note that VI is not uniquely defined, but Vi VT" = V is if the leading eigenvalue is distinct and consequently, norms of the form ||V]x =ll for IRd are dell-defined. Let À1 the cigenvalues of 2 and A be a diagonal matrix with (i, = (1 + 7h A4111/ /n////n//// (1 + The / nn).
X denote
i ,,...,,,,,,,,,,,,, Also let
M E [V xx
(6)
Now we define
V
Un Σ EVV V - E) - Σ)V -
Un V1 EaT VT
(7)
We have the following result:
Theorem 1. Suppose that uo is drawn from the uniform distribution on Sd-1 A1 O(1).
Choose On - oo such that nd exp(-1m (A1 X2)) + 0, (on' log nf (M331) + 0. where
Ma = XX - >. Further, let Zn, be a mean 0 Gaussian matrix such that Var(Zn)
XX - E) and suppose that:
|M] 0
(8)
E V 2m 110
VE - 11
o(n)
3
5

---

Then, for sequence of Gaussian distributions (Znbx2 with mean 0 and covariance matrix Vn (see 7 the following holds:
sup SR
P (n/nn, sin 11) s P(Z Zn S -
(10)
Theorem T is very general. We allow the dimension to grow with the number of observations, which typical in the high-dimensional bootstrap literature. Note that the case of fixed d and growing is also a special case of this setup.
We want to point out that while previous literature obtained sharp bounds on the sin" error 1- (1 01)². we go a step further. We establish an approximating distribution for 11/1n(1 (1 o1)).
Remark 1 (Condition on norm). For simplicity, we assume X1 O(), which can be easily relaxed to grow slowly with n. We do not assume that the ||X,X;" - Dill 2 is bounded almost surely. However, the norm of X;XT E comes into play implicitly via the assumption in Eq 9 Consider the case where X; are drawn from some multivariate Gaussian distribution. We use this to build intuition about the assumptions in Eq and 9
In this case, X is a Gaussian of 3
independent entries and thus TE (
11 [X) v)) NNtee
that E}>1((33 0z)² X;) is a sub-exponential random variable with parameters (c1 Es>1 Ag. 02). Furthermore, ||MI] = A1 Ei>ı X Thus Eq 9 reduces to checking if
2//2 (E}1 ((()
o(n)
E, E,230//2
Remark 2 (Coordinates with summable sub-Gaussian parameters). Eq 9 imposes a growth condition on the moments of both the data and Gaussian analog. One setting for which both growth rates are in fact bounded is if the coordinates of X are nab-Gaussian and the sub-Gaussian parameters satisfy
E 1/ < C < oo 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 's algorithm and have the advantage that they require 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) < The < n f3, 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 future work.
As a corollary of our main theorem, we obtain the following error bound on the sin error.
Corollary 1. Under the conditions in Theorem 7
we have
The Ma
sin (01, 01) Op
(1())-22)
Remark 4 (Comparison with previous work). As a byproduct of our analysis, we recover the sharpest convergence rates for Oja's 's algorithm in the literature. If we set 1hn c log rdd X2), for large
Ma lognd
enough C1, the dominating term in the error is Op ne\i 22)2
under mild conditions on d. This
matches the bound in 125).
Remark 5 (Rate of convergence in Kolmogorov distance). To simplify the theorem statement, have stated Theorem D 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 (5T].
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 lower bound on the Frobenius norm of the covariance matrix as in Eq g Note that if indeed all X, 's were a scalar multiple of Ui, then the V a mairix 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 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 digh-dimensional CLT literature 77 8).

---

Now we provide proof sketch of Theorem below.
Proof sketch for Theorem ! We provide the main steps in our derivation. The detailed calculations can be found in Supplement Section
1. We start by expressing the sin error quadratic form:
sin sin"(u1, 01) = " B 01 ofee Bo
w6 BT (1 Bn
m/ BTT u0
ub B Bn
(11)
(V, VV Bn un) )"() VF BauD)
|B
where in the last line we used the fact that V1 V is idempotent. Our proof strategy for the central limit theorem involves further approximating Eq 11 with an inner product of the Hájek projection (first-order) term in Eq 5
2. Our second step is to show that || uo concentrates around its expectation (1 mm>1/n) unl-
Next we establish that [ V. V.V, is Op (Vd expl expl-rm(A (A1 A2)}
2: M log d
n
This is achieved by using similar recursive argument as in 225] but with the crucial observation that the residual common difference term is of lower order because it can be replaced by matrix product minus its expectation.
4. Now go back to the expansion in Eq
V. V Bav = (vf vo) Σ
EViVis Vi
Since To (I + nn nn/nE)" V V is the vector. Now examine the (uf uo) VVLVV Ti term. Here use the structure of the higher order terms Tg. In particular, use the fact that it is matrix product interlaced with k X,X? - Σ matrices. For example, for we have
Nh n
Y, Da 5-1-1 1 Y, (1 " n-j
Σ KE
I+
[I+
We show that the norm of (vf VL Ti)e1, normalized by the denominator, is O M3 /n²). The fact that the summands of Tj uncorrelated and T₃ and Te uncor- related for + makes this possible.
5. Finally, we are left with Vi ViTie (oTu))) uo). Note that this is of the following form:
Th (of uo VT snsggn uo)
VIAA
-Emr
n lof unl + Amp./))n
not hard to see that this sum of independent random vectors with covariance matrix Th/ nn/nV, (see Eq 7
We adapt 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 D and the conditions of distributional convergence satisfied.
Finally, all the error terms are combined along with an cntieconcentration argument for x² to establish the final result. The full proof and accompanying lemmas are in Section B of the Supplement.

---

3.2 Bootstrap consistency
Using the weighted x² approximation for inference requires estimating the eigenvalues of Σ and other population quantities: however, accurate estimates may not be available in streaming setting. Instead, we propose a streaming bootstrap procedure that mimics the properties of the original Oja algorithm. While similar structure leads to error terms that are similar to the CLT, the analysis the bootstrap presents its own technical challenges, In what follows let P* denote the bootstrap measure, which is conditioned on the data, and let [.] 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 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 and approximating Zn, from the CLT) that only depend on how close the covariance matrices E* [Z, Ze"" and EL Z, Zn] 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 sin error:
Lemma 1. [Bounding the difference between the bootstrap covariance and true covariance] Let:
Zn sgn(u] V The Σ V X
(12)
Recall the definition of from Eq] We have,
B| - Σ11
[trace (TE* 17
V
m(Ai - x22)
With this lemma in hand, we are ready to state our bootstrap result.
Theorem 2 (Bootstrap Consistency). Suppose that the conditions of Theorem are satisfied. Furthermore, let On be sequence such that P(As) - O, where An is defined An (a VM log am + 0,
{maxigigo ||x:1|2 S an). Further suppose that Malog' o2 + 0, May? + and n(A1-2a A2) m(A1-)2)
x||xxxx-=|[1 + Then,
Pd sinnnniinn sin sin (ef, (1), st) - P(nn/
v1) 1 =() 0
sup tER
Proof sketch of Theorem analysis in Supplementary Section. We use a bootstrap version of the Hoeffding decomposition conditioned the data, stated in Supplement Section. In step have B,, replace B,, where B", is given by:
The proof follows similar route to Theorem We provide detailed
7 B = II (1+ 1=1
WXXXXX - X-IXZ(
We work out Step using concentration of matrix products (23 For steps 2-3, we see that Tg has the same structure as Ta with the difference that (I + 1)n E/n)* 2 is replaced by its sample counterpart which is a product of independent matrices of the form I + Tm, / 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 B", uo is given by Vmn/nZ; where Z is given in Eq 12 Conditioned the data, this already Normally distributed since the multiplier random variables W; are themselves Gaussian. We then invoke the Gaussian comparison result Lemma obtain convergence the weighted x² approximation.
We now make a couple of points regarding our analysis. 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-mont momontt methods may be used. We also operate 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 A. in an abstract manner, but one may wonder how stringent the condition is in various problem setrings. Below, we describe a general setup with grows as log n; under milder forms of various decay, all we
sub- Gau-sass entries of X; in which need for an to grow slowly with n. Here ||-1|a, is the sub-Exponential Orlicz norm and Illlaa is the sub-Gaussian Orlicz norm (see, for example (50))
Proposition (The effect variance decay on the norm). For each 1 S j S p, suppose that X1j satisfies [X1jll S Mj E)=1"} I Ci cc. Then, for some universal constant C2 > 0.
- E(Xi,D) < C2 and for some c1, C2 > 0,
max I8i5n
|x,1|
01 log
S
C2
n
We now present experimental validation of our bootstrap procedure below.
4 Experimental validation of the online multiplier bootstrap
We draw Zıj IID Uniform(-V3, 3). for i 1.... n and j = 1,... d. Consider a PSD matrix Kij exp( - 3 jlc) with 0.01. We create covariance matrix such that Enj k(i, (i, j)o We consider 0; 5i for 3 = 0.2 and 1. Now transform the data introduce dependence by letting X, = E1/2 74. By construction, we have that E[X,
Σ for all 1 < i S n. Our goal is to simply demonstrate that the bootstrap distribution of sin" errors closely match that of the sampling distribution. To this effect, we fix ug and draw 500 datasets and run streaming PCA on each and then construct empirical CDF (P)) from the sin? error with the true V1. This is the point of comparison for the bootstrap distribution (F"), for which we dataset X. We then invoke algorithm Too obtain 500 bootstrap replicates of well as the Oja vector for the dataset 0.. The bootstrap distribution is the empirical CDF of - [6] 0.) We use 77n = log n. In Figure ! 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 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 Pp(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.
Discussion
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-dimensiooal 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-dimension principal subspaces, and additional tools will be needed to account for non-independent matrix products which appear beyond the setting of online PCA.

---

Besteraped GBF
werping
08
0.
compiin GUP
900
400
600
E00
1000
(A)
(B)
True wwwgling CD
True sampling CDT
0.5
0.5
6.7
0.
(C)
Figure Bootstrapped and sampling CDF for 1000, d 500 in (A) and (C) and for 10, 000,00044 d = 500 in (B) and (D). (A) and (B) use 8 = 0.2 whereas (C) and (D) use 3 = 1.
Acknowledgment
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 1199550-19-100000 NSF DMS 1952735, NSF HDR-1934932, and NSF 2019844.
References
[1] Z. Allen-Zhu, Y. Li, R. Oliveira, and A. Wigderson. Much Faster Algorithms for Matrix Scaling. In Proceedings of the 58th Symposium on Foundations of Computer Science, FOCS '17, 2017.
[2] M.-F. Balcan, S. S. Du, Y. Wang, and A. W. Yu. An improved gap-dependency analysis of the noisy power method. In V. Feldman, A. Rakhlin, and O. Shamir, editors, 29th Annual Conference on Learning Theory, volume 49 of Proceedings of Machine Learning Research, pages 284-309, Columbia University, New York, New York, USA, 23-26 Jun 2016. PMLR.
[3] A. Balsubramani, S. Dasgupta, and Y. Freund. The fast convergence of incremental PCA. In C. J. C. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K. Q. Weinberger, editors, Advances in Neural Information Processing Systems 26, pages 17843182.. Curran Associates, Inc., 2013.
[4] L. Balzano, R. Nowak, and B. Recht. Online Identification and Tracking of Subspaces from Highly Incomplete Information. arXiv e-prints, page arXiv: 1006.4046, June 2010.
[5] Y. Benoist and J.-F. Quint. Random Walks on Reductive Groups. Springer International Publishing, June 2016.
[6] X. Chen, J. D. Lee, X. T. Tong, and Y. Zhang. Statistical inference for model parameters in stochastic gradient descent. Ann. Statist,, 4811)2251- 02 2020.
[ V. Chernozhukov, D. Chetverikov, and K. Kato. Gaussian approximations and multiplier bootstrap for maxima of sums of bigh-dimensional random vectors. The Annals of Statisties, 41(6):2786 2819, 2013.
10

---

[8] V. Chernozhukov, D. Chetverikov, and K. Kato. Central limit theorems and bootstrap in high dimensions. Ann. Probab., 45(4):2309-2333 07 2017.
[9] B. Efron and R. Tibshirani. Bootstrap methods for standard errors, confidence intervals, and other measures of statistical accuracy. Statist. Sci., 1(1):54-75, 02 1986.
[10] J. Emme and P. Hubert. Limit laws for random matrix products. arXiv e-prints, page arXiv:1712.03668, Dec. 2017.
[11] V. Fabian. On asymptotic normality in stochastic approximation. Ann. Math. Statist., 39(4):3327- 1332, 08 1968.
[12] Y. Fang, J. Xu, and Yang. Online bootstrap confidence intervals for the stochastic gradient descent estimator. J. Mach. Learn. Res., 19(1:30553377, Jan. 2018
[13] P.R.S. Liii on lines and planes of closest fit to systems of points in space. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 2(11):559-572, 1901.
[14] A. Furman. Chapter 12 random walks on groups and random transformations. volume I of Handbook of Dynamical Systems, pages 931 1014. Elsevier Science. 2002.
[15] A. Giessing and J. Fan. Bootstrapping e statistics high dimensions, 2020.
[16] Gtze and H. R. Knsch. Second-order correctness of the blockwise bootstrap for stationary observations. Ann. Statist., 24(5):1914-1933 10 1996.
[17] P. Hall. The Bootstrap and Edgeworth Expansion. Springer- Verlag, New York, 1992.
[18] M. Hardt and E. Price. The noisy power method: A meta algorithm with applications. In Z. Ghahramani, M. Welling, C. Cortes, N. D. Lawrence, and K. Q. Weinberger, editors, Advances Curran Associates, Inc., 2014.
Neural Information Processing Systems 27, pages
[19] R. Helmers. the edgeworth expansion and the bootstrap approximation for a studentized u-statistic. Ann. Statist., 19(1):470-444 03 1991.
[20] A. Henriksen and R. Ward. AdaOja: Adaptive Learning Rates for Streaming PCA. arXiv e-prints, page arxiv: 1919912111,,, 12115, May 2019.
[21] A. Henriksen and R. Ward. Concentration inequalities for random matrix products. arXiv e e-prints. page arXiv:1907, 00803,0533, July 2019.
[22] W. Hoeffding. A class of statistics with asymptotically normal distribution. Annals of Mathe- matical Statistics, 19(3):293-325, 09 1948.
[23] D. Huang, J. Weed, J. A. Tropp, and R. Ward. Matrix concentration for products, 2020.
[24] D. Huang, Niles-Weed, and R. Ward. Streaming k-PCA: Efficient guarantees for Oja's algorithm, beyond rank-one updates, 2021.
[25] P. Jain, C. Jin, S. Kakade, P. Netrapalli, and A. Sidford. Streaming PCA: Matching matrix bernstein and npar-optimal finite sample guarantees for Oja's algorithm. In Proceedings of The 29th Conference on Learning Theory (COLT), June 2016.
[26] P. Jain, D. Nagaraj, and P. Netrapalli. SGD without Replacement: Sharper Rates for General Smooth Convex Functions. arXiv e-prints, page 1eXiv:1903..m 01463,, Mar. 2019.
[27] P. Jain, P. Netrapalli, and S. Sanghavi. Low-rank matrix completion using alternating mini- mization. In Proceedings of the Forty-Fifth Annual ACM Symposium on Theory of Computing, STOC 13, page 665674, New York, NY, USA, 2013. Association for Computing Machinery.
[28] I. T. Jolliffe and J. Cadima. Principal component analysis: a review and recent developments. Phil. Trans. Soc. A., 2016.
[29] R. Keshavan, A. Montanari, and S. Oh. Matrix completion from a few entries. Information Theory, IEEE Transactions on, 56:2980 2998, 07 2010.

---

[30] F. Ledrappier. Some asymptotic properties of random walks on free groups, pages 117-152. 06 2001.
[31] C. J. Li, M. Wang. H. Liu, and T. Zhang. Near-optimal stochastic approximation for online principal component estimation. Math. Program., 167(1):75-97, 2018.
[32] Li, Liu, Kyrillidis, and Caramanis. Statistical inference using SGD. arXiv e-prints, page N.Kiv:1705.07444 17777,,a, May 2017.
[33] X. Liang. On the optimality of the Oja's algorithm for online PCA, 2021.
[34] M. Lopes, S. Wang, and M. Mahoney. Error estimation for randomized least-squares algorithms via the bootstrap. volume 80 of Proceedings of Machine Learning Research, pages 3217-3226, Stockholmsmssan, Stockholm Sweden, 10-15 Jul 2018, PMLR.
[35] M. E. Lopes. Estimating the algorithmic variance of randomized ensembles via the bootstrap. Ann. Statist., 4782):1088-1112 04 2019.
[36] M. E. Lopes, N. B. Erichson, and M. W. Mahoney. Error Estimation for Sketched SVD via the Bootstrap. In Proceedings of the International Conference of Machine Learning, 2020.
[37] M. E. Lopes, S. Wang, and M. W. Mahoney. A bootstrap method for error estimation in randomized matrix multiplication. Journal of Machine Learning Research, 20(39):1-40) 2019.
[38] I. Mitliagkas, C. Caramanis, and P. Jain. Memory limited, streaming PCA. In Proceedings of the 26th International Conference on Neural Information Processing Systems Volume 2, NIPS' 13, page 28862894, Red Hook, NY, USA, 2013. Curran Associates Inc.
[39] E. Moulines and R. Bach. Non-asymptotic analysis of stochastic approximation algorithms for machine learning. In J. Shawe-Taylor, R. S. Zemel, P. L. Bartlett, F. Pereira, and K. Q. Weinberger, editors, Advances in Neural Information Processing Systems 24, pages 451-459.
Curran Associates, Inc., 2011.
[40] A. Nemirovski, A. B. Juditsky, G. Lan, and A. Shapiro. Robust stochastic approximation approach to stochastic programming. SIAM J. Optimization, 19(4):1574-109, 2009.
[41] E. Oja. Simplified neuron model as principal component analyzer. Journal of Mathematical Biology, 15(3)2267-273, Nov. 1982.
[42] B. T. Polyak and A. B. Juditsky, Accelcration of stochastic approximation by averaging. SIAM
J. Control Optim., 30(4):838855, July 1992.
[43] D. Pouzo. Bootstrap consistency for quadratic forms of sample averages with increasing dimension. Electronic Journal of Statistics, 9(2)30046 - 3097, 2015.
[44] H. Robbins and Monro. A stochastic approximation method. Ann. Math. Statist., 22(3))400- 407, 09 1951.
[45] D. Ruppert. Efficient estimations from slowly convergent Robbins-Momet process. 1988.
[46] C. Sa, Olukotun, and C. Ré. Global convergence of stochastic gradient descent for some nonconvex matrix problems. CoRR, abss/4111 1111111144 2014.
[47] O. Shamir. Convergence of stochastic gradient descent for pea. In Proceedings of the 33rd International Conference on International Conference on Machine Learning Volume 48, ICML' 16, page 257265. JMLR.org. 2016.
[48] W. Su and Y. Zhu. Uncertainty quantification for online learning and stochastic approximation via hierarchical incremental gradient descent, 2018.
[49] A. der Vaart. Asymptotic statistics. Cambridge University Press, 2000.
[50] R. Vershynin. High-Dimensional Probability. Cambridge University Press, Cambridge, UK,
2018.
12

---

[51] M. Xu, D. Zhang, and W. B. Wu. Pearsons chi-squared statistics: approximation theory and beyond. Biometrika, 106(3):16-723, 04 2019.
[52] L. Yu, K. Balasubramanian, S. Volgushev, and M. A. Erdogdu. An analysis of constant step size sgd in the non-convex regime: Asymptotic normality and bias, 2020.
13


Reviews:
[{"review_id": "xIfTTWIem_Y", "reviewer": "Reviewer_Kj1d", "summary": "In this paper the bootstrapped Oja scheme is considered. For the squared sin distance between the Oja’s scheme output and the first eigenvector of the covariance matrix of observations, the authors provide a Gaussian approximation and demonstrate bootstrap consistency. ", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The authors consider a problem of quantification of uncertainty for Oja scheme in the framework for the streaming PCA. With this aim, the multiplier bootstrap is introduced in the update step of the Oja scheme. The proof of asymptotical Gaussian approximation for the squared sin distance and the asymptotical first term expansion of the squared sin error is based on the Hoeffding decomposition of the update matrices and high-dimensional CLT for the quadratic forms of random vectors and concentration and anti-concentration bounds. For the bootstrapped Oja the standard techniques lead to the bootstrap consistency result. In the end, the authors demonstrate the convergence of the bootstrapped CDF to the sampling CDF for the case of well-posed problems and further discuss the limitations.\n\nThe paper is overall quite well written. The following points could be corrected and clarified:  \n1) S and S_k in the Eq. 5 are not defined \n2) In Algorithm 1 in the update step a vector estimate is missing in the brackets\n3) In Algorithm 1 Z_i -> W_i \n4) In Remark 2, which is the referenced result? In [28] the error has log (d) in the nominator, and in the Remark 2 it is log(n)  \n5) I am not sure I understand the condition (10) in view of (11), does the left-hand side of (10) contain anything depending on n?  \n6) In supplementary material in Proposition A2 a number of the equation is missing\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The learning rate in the considered version of Oja scheme assumes the knowledge of the number of upcoming observations (e.g. see [28]). \nWould not it be a limitation for the streaming PCA problem?   ", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "n-A_t1WGFnz", "reviewer": "Reviewer_tCtM", "summary": "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. ", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper's main contribution is theoretical, and so to truly appreciate this paper, one must get into the technical details. The authors take the standard point-of-view and analyze Oja's algorithm as a limit of products of random variables. I am very impressed that they were able to establish the distribution of the sin^2 error between the leading eigenvector of X and the computed vector via Oja’s algorithm (see Theorem 1). It is implicit as the probability distribution is connected to the eigenvalues and eigenvectors of Sigma, limiting its practical application. To have some practice application of Theorem 1, the authors propose a streaming bootstrap procedure (see Theorem 2). \n\nThe paper explains an overview of their analysis. After n iterations, Oja's algorithm can be expressed as a normalized matrix-vector product with a matrix B_n. Using Hoeffding decomposition, the authors write B_n as a sum of matrices that are uncorrelated (see (5)). This allows them to compute second-moment information of norms of B_n and consequential the distribution of the sin^2 error. For the bootstrapping idea, the authors use a comparison lemma. \n\nWhile I don't fully appreciate all the 36 pages of technical details in the supplementary materials, I hope that the authors can better describe their final practical algorithm. In fact, I have several confusions involving Algorithm 1:  \n\nThe input is u_0, but at the same time u_0 is assigned g/||g||. I think the authors want u_0 to be uniform on the sphere, so it shouldn't be an input. Also, what is Z_i in the algorithm (how to estimate their distribution)? I encourage the author's to write \"Compute Y_t = X_tX_t^T - \\Sigma\" in the algorithm? How does one compute Sigma in Algorithm 1? What is the integer p in Algorithm 1? Theorem 1 needs the stepsize to go to infinite, but in Algorithm 1, it is constant. While Algorithm 1 presents the stepsize as constant, Section 4 uses log(n) (but what happens at n = 1?). Also, Theorem 1 presents an asymptotic result, but do the authors understand the optimal step size for the fastest convergence in probability? \n\nHere are a few grammar errors that I spotted while reading. I am sure the authors will want to correct them. \n \nl30: Define X. \n\nl65: \"convergence rates which deteriorate\" should read \"convergence rates that deteriorate\"\n\nl105: \"where the the random\" should read \"where the random\" \n\nl161: \"approximation to the distribution\" should read \"approximation of the distribution\"\n\nl163: \"wew update\" should read \"we update\"\n\nl189: \"growly wih n\" should be \"grow with n\"\n\nl191: \"multivariate Gaussian distribution\", do the authors mean \"multivariate standard Gaussian distribution\"? Otherwise, is the statement true?  \n\nl249: \"parameters of this distributions\" should read \"parameters of this distribution\"\n\nl260: \"we use the later\" should read \"we use the latter\"\n\nl263: \"that that of the\" should read \"that of the\"\n\nl276: \"For step 2-3\" should read \"For steps 2-3\"\n\nReferences: put captials for: oja, bernstein, pca, fpca, robbins-monro, etc.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The expected limitations were discussed such as the analysis is only for the first principal vector of Oja's algorithm.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6"}, "extra_text": {}}, {"review_id": "TDbAUzrmFWM", "reviewer": "Reviewer_H32z", "summary": "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. \n", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "\nStrength:\n1. Bootstrapping the error of Oja's algorithm is an interesting and important problem. The paper provides solid contributions to this problem.\n2. The distribution of $\\sin^2$ error of the leading eigenvetor of Oja's algorithm is novel.\n\n\nWeakness:\n1. The writing is not good. Some notations are confusing. For example, the variables $Y_t$ and $Z_i$ in Algorithm 1 are not defined. Also, some expectation notations lack brackets. For instance, in line 259, $\\mathbb{E}ZZ^T$ should be $\\mathbb{E}[ZZ^T]$.\n2. The second comment is still about the writing. In Figure 1, what does the x-axis represent? Why the range of the x-axis of figures (C) and (D) are [0,1]? I suggest the authors can provide more descriptions on the figures.\n3. I'm confused about the data generation of the experiment. The paper chooses $X=Z(\\Sigma^{1/2})^T$, which indicates that $\\mathbb{E}[X^T X]=\\Sigma^{1/2}\\mathbb{E}[Z^T Z]\\Sigma^{1/2}$. Since $\\mathbb{E}[Z^T Z]=nI$, the expectation of$X^T X$ should be $n\\Sigma$ rather than $\\Sigma$. Would you please tell me if I have missed something?\n\nOverall, I think the paper provides some solid contributions to an important problem, but the writing of the paper should be polished.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "See main review", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6"}, "extra_text": {}}, {"review_id": "Sf5gIw00dxY", "reviewer": "Reviewer_2S7D", "summary": "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.", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "This is an interesting and well-written paper which balances rigorous theory and intuition in just the right way. I particularly enjoyed how the authors draw from classical statistical theory (U-statistics, Hoeffding decomposiiton) and modern theory (Gaussian approximation, hiogh-dim. bootstrap) to solve a relevant problem in machine learning.\n\nI have the following comments:\n1) Unfortunately, I do not understand Algorithm 1: Where are the $W_i$'s used? What are $Y_t$ and $Z_i$? I would appreciate if the authors could explain and clarify the notation.\n2) Reg Theorem 1: To choose $\\eta_n$ one has to gauge $\\lambda_1 - \\lambda_2$ and $M_d$. How do the authors do this in practice? How sensitive is their algorithm to values of $\\eta_n$?\n3) Line 187: Technically, convergence of the Kolmogorov distance in eq. (11) does not imply that the authors \"establish the [limiting] distribution of $n/\\eta_n(1- (v_1^T\\hat{v}_1)^2)$\" because the limiting distribution of $Z_n^TZ_n$ need not to exist. It would be more accurate to state that in the limit the difference of the distributions of $n/\\eta_n(1- (v_1^T\\hat{v}_1)^2)$ and $Z_n^TZ_n$ in Kolmogorov distance vanishes. And that's good enough for most (all?) statistical purposes.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "I wished the authors had discussed how to choose the tuning parameters of their algorithm and included a more comprehensive simulation study. (The current simulation is based on a covariance matrix with very light correlation/ almost sparse. How does their method work on dense and/ or equicorrelated covariances)?", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "1.5"}, "extra_text": {}}]

Author rebuttal:

Meta-review: The reviewers were unanimous in their support for this paper, with a clear theoretical contribution to our understanding of an important algorithm. Most reviewers also felt the paper was well-written, and we thus feel it will be broadly interesting to the NeurIPS community.   
Decision: Accept (Spotlight)