A Fast and Accurate Estimator for Large Scale Linear Model via Data Averaging
Authors: Rui Wang
Yanyan Ouyang
Panpan Yu
Wangli Xu
Venue: NeurIPS 2023 poster (2023)
Keywords: Big data
Data averaging
Order statistic
Sampling method
Sketching method.
TL;DR: We propose a new sketching method for large scale linear model based on data averaging, which can achieve a faster convergence rate than the optimal convergence rate for sampling methods.

Abstract: This work is concerned with the estimation problem of linear model when the
sample size is extremely large and the data dimension can vary with the sample
size. In this setting, the least square estimator based on the full data is not feasible
with limited computational resources. Many existing methods for this problem are
based on the sketching technique which uses the sketched data to perform least
square estimation. We derive fine-grained lower bounds of the conditional mean
squared error for sketching methods. For sampling methods, our lower bound
provides an attainable optimal convergence rate. Our result implies that when the
dimension is large, there is hardly a sampling method can have a faster convergence
rate than the uniform sampling method. To achieve a better statistical performance,
we propose a new sketching method based on data averaging. The proposed
method reduces the original data to a few averaged observations. These averaged
observations still satisfy the linear model and are used to estimate the regression
coefficients. The asymptotic behavior of the proposed estimation procedure is
studied. Our theoretical results show that the proposed method can achieve a
faster convergence rate than the optimal convergence rate for sampling methods.
Theoretical and numerical results show that the proposed estimator has good
statistical performance as well as low computational cost.

A Fast and Accurate Estimator for Large Scale Linear Model via Data Averaging
Rui Wang
Yanyan Ouyang
Center for Applied Statistics
Center for Applied Statistics and School of Statistics
and School of Statisties
Renmin University of China
Renmin University of China Beijing 100872, China s staoyyygruc.edu. cn
Beijing 100872, China 446100240@qq.com com
Panpan Yu NavInfo
Wangli Xu
Center for Applied Statistics and School of Statistics Renmin University of China
Beijing 100094, China
yupanpan@navinfo.com
Beijing 100872, China wlxuêruc edu, cn
Abstract
This work is concerned with the estimation problem linear model when the sample size is extremely large and the data dimension can vary with the sample size. In this setting, the least square estimator based on the full data is not feasible with limited computational resources. Many existing methods for this problem are based on the sketching technique which uses the sketched data to perform least square estimation. We derive fine-grained lower bounds of the conditional mean squared error for sketching methods. For sampling methods, our lower bound provides attainable optimal convergence rate. Our result implies that when the dimension is large, there hardly a sampling method can have a faster convergence rate than the uniform sampling method. To achieve better statistical performance, we propose a new sketching method based on data averaging. The proposed method reduces the original data few averaged observations. These averaged observations still satisfy the linear model and are used to estimate the regression coefficients. The asymptotic behavior of the proposed estimation procedure is studied. Our theoretical results show that the proposed method can achieve a faster convergence rate than the optimal convergence rate for sampling methods. Theoretical and numerical results show that the proposed estimator has good statistical performance well low computational cost.
Introduction
Linear regression model is one of the simplest and the most fundamental models in statistics and machine learning. Suppose one collects independent and identically distributed (i.i.d)) observations {Z,, y where Z, IP' is the vector of the predictors and ys E R is the response. The linear model assumes
N
y: 80 + Z, B1 + E5. i i= 1= N,
(1)
where 3n E TR,33 T&P are the unknown coefficients and c EN random variables represent- .l.....yv) yN) and X (X.,... Xx)T The classical least square estimator equals arg min BERF+1 Ily - X . If X has full column rank, the least square estimator equals (ii", X.XT) MN' X.y% and direct computation costs O(Ny²)
ing noise. Let X; (1, ) T, (Ba, ,3 T. (yI.....
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

time. For large scale linear models with N > p. the computing time O(Np²) of the exact least square estimator is not negligible. Faster estimators of [3 can largely facilitate the practical data analysis pipelines.
Numerous research efforts have been devoted to the estimation problem for large scale linear model. Many existing work in this area can be understood as matrix sketching methods which explicitly or implicitly use matrix sketches as surrogates for the original observations reduce the data size. Specifically, sketching methods solve the sketched least square problem
min BERFIL
10Ty y O O'xB1²
where O E RRVxm is a sketching matrix with n < N. The solution to the problem (2) is the least square estimator based on the reduced data O X E [=*((++1) and O y E R". Since n < NN, the sketched least square problem can be solved much faster than the least square estimator for the full data. In this paper, we only consider the case that O is independent of 61, .1.....EN. That is, O may rely on X, but not rely on y. This guarantees that the solution to 2 is an unbiased estimator of f.. Note that sampling methods are special cases of the sketching framework 2 In fact, for sampling method, each column of O is a vector whose elements are all 0 except one that equals 1, Sketching
methods have been intensively researched in algorithmic aspect; see Mahoney (2010). Woodruff
[2114]. Drineas and Mahoney (016) for reviews. Recently, the statistical aspect of sketching methods
also draws much attention; see, e.g.. Ma et al. (2015), Raskutti and Mahoney [0166] Wang et al.
2017 Dobriban and Liu 2019). Ma et al. 2020). Ahfock et al. [2021].
Probably the simplest sketching method is the uniform sampling method which randomly selects n observations with equal probability to form the reduced data. Recently, Pilanci and Wainwright (2016) provides minimax lower bound for the mean squared prediction error of random sketching methods. Theorem 1 of PPllaci and Wainwright (2016) shows that for a large class of random sketching methods, including many existing data-oblivious sketching methods and sampling methods, the convergence rate of the mean squared prediction error can not be faster than the uniform sampling method. Hence it is nontrivial task construct a sketching method which has significantly better statistical performance than the uniform sampling method.
Recently, Wang et al. [019]] initiates the study of sampling methods based extreme values. Motivated by the D-optimal criterion. Waan et al. (2019) proposed the information-basee optimal subdata selection (IBOSS) algorithm which successively selects informative observations based on extreme values of variables. They showed that for fixed p. the estimator produced by the IBOSS algorithm can have faster convergence rate than the uniform sampling method. Meanwhile, the computation of the IBOSS algorithm can be completed within O(Np+ +np²) time which has the same order as the uniform sampling method if cN for some constant e >> O. Now the algorithm of Wang et al. [019 has become the building block of some recent methods for large scale problems. For example, Wang [019) proposed an algorithm which combines the algorithm of WWang al. [2019] and the divide and conquer strategy. [Cheng et al. [2020) extended the algorithm of Wang al. [2019] to the logistic regression model. Existing asymptotic results for the IBOSS algorithm are obtained in the setting of fixed and p. At present, there is still a lack of theoretical understanding of the behavior of the IBOSS algorithm in the setting of varying n and p.
The IBOSS algorithm is a sampling method and is therefore an instance of the sketching framework (2). Interestingly, the IBOSS algorithm can surpass the minimax lower bound of Pilanci and Wainwright 2016). In fact, a key condition for Theorem 1 Pilanci and Wainwright (2016) does not hold for the IBOSS algorithm. Thus, the IBOSS algorithm is not restricted by the minimax bound of Pilanci and Wainwright (2016). This fact is detailed in Section 2 Note that there are many potential sketching methods which are not restricted by the minimax bound of Pilanci and Wainwright (0016). To give a more comprehensive understanding of the behavior of these sketching methods, we derive fine- grained lower bounds for the conditional mean squared error of the sketched least square estimators produced by ( in the setting that Z; is a standard normal random vector. In particular, our result provides a lower bound for any sampling method which may possibly rely on X but does not rely on y. It turns out that if p < log (N/). then the optimal lower bound for sampling methods can have a faster convergence rate than the uniform sampling method. On the other hand, if log(N/n) < p. any sampling method can not largely surpass the uniform sampling method. Furthermore, we derive the asymptotic behavior of the IBOSS algorithm in the setting of varying n and p. It turns out that under certain conditions, the IBOSS algorithm can achieve the optimal rate for sampling methods.

---

Table 1: Theoretical performance of the ideal sampling method (abbreviated as ISM), the IBOSS
algorithm and the proposed method when Z1 [Op,, Ip). Assume that as N + 00, p -> oo,
(log(())) (p))* log(N) /NN - oo, n = O(N9) for some 0 < é < 1/2 and p = 0(n2/2-<') for some 0 <E E <1 1/2, log(N/n) = O See Theorems 2 3 The reported computing time is under the assumption that the multiplication of an m x n matrix and an n x matrix costs ((mpp) time, and
the inversion of a p x p matrix costs (pp") time.
Methods
Reduced sample size
E |8
-
z
Computing time
ISM
1
Op
"(p+log( 1)
IBOSS
Op
O(Np )
n(ptlog( 2
NEW
2p
(1 + op(()) 2 21og(2p)
O(Np+ +p")
For large scale linear models, is often the case that log(N/n) & p. In this case, any sampling method can not have a significantly better statistical performance than the uniform sampling method. Inspired by this phenomenon, we propose an alternative sketching method which can reduce the full data to just a few observations while the resulting estimator of ]3 may have smaller conditional mean squared error than sampling methods. The proposed method is based on data averaging. The main idea is to partition the observations into 2p groups such that the averages of Z; within groups are separated. The least square estimator based on 2p averaged observations is used to estimate B.. The computation of the proposed method can be completed within O(Np + p²) time. Our theoretical results show that the proposed method can have faster convergence rate than any sampling methods with comparative computing time. Also, the proposed method reduces the full data to merely 2p averaged observations. These averaged observations also satisfy the linear model I and have independent errors. Consequently, it is convenient to further compute other estimators or conduct statistical inferences using the reduced data. The good performance of the proposed estimator also verified by simulation results and a real data example. Table I summarizes the theoretical performance of the proposed method and compare with the ideal sampling method implied by Theorem 2 and the IBOSS algorithm.
The rest of the paper is organized as follows. Section investigates lower bounds for the conditional mean squared error of the sketched least square estimators produced by ( In Section 33, we propose a data averaging method to estimate 33 and investigate its asymptotic behavior. Section presents the simulation results briefly. Section 5 concludes the paper. The simulation results, a real data analysis and all proofs are deferred to the Supplement Material.
We close this section by introducing some notations and assumptions that will be used throughout the paper. For any real number w, let [w] denote the largest integer not larger than w. For any vector W. let ||W]| denote the Euclidean norm of W. For any matrix B, let |[B]] and [[]|| denote the operator norm and the Frobenious norm of B, respectively. Moreover, denote by B. the jth column of B. If B is symmetric, denote by (B) the ith largest eigenvalue of B. In this paper, the symmetric matrices are equipped with Loewner partial order. That is, for two symmetric matrices B and B.g B, B2 if and only if B B, is positive definite. For positive semidefinite matrix B, let B1/2 denote a positive semidefinite matrix such that (B1/2)2 B. For any set s, denote by cf its complement and (ard(of) ) its cardinality. Let (xx) and (px)) denote the distribution function and density function of the standard normal distribution, respectively. For random variables & R and > 0, g = 0p (π) means that e/n converges to 0 in probability, and = Op(n) means that e/n is bounded probability.
Let N denote the size of full sample, p denote the dimension of covariates. Let Z = (Z,,..., 2N) T be an N x p matrix of covariates. Denote by zı, the jth element of 7,, i 1,..., N, 1, ,,...,P. 2(N),, denote the order statistics of { V j 1,...,p. p. The following assumption for the data distribution is assumed throughout the paper.
Let 2(1),)
3

---

Assumption Suppose (Z, are i...d. and satisfy the linear model ① where E(E1) = 0, Var(EI) o and E is independent of Z, Suppose Z has a density function (1,,.1)...21..p))
with respect to the Lebesgue measure on RR, and E. = js, Cov(Z1) E are finite. Suppose Σ = (d%.3) is positive definite. Suppose r > O. As N - oc, the dimension p a function of N, while the distribution of (Z1,3) 91) only relies on Finally, assume o² is a constant which does not depend on N.
For simplicity, our notations suppress the dependence of p on N, and the dependence of the distribu- tion of (Zi,3) on p.
2 Risk bounds for sketched least square estimators
Theorem of Pilanci and Wainwright 2016] provides a minimax lower bound for the mean squared prediction error of random sketching methods. Their result implies that under certain conditions, there exists constant C' >> O such that for any estimator 3 which only relies (O OTy)
Cp de
sup B I 3)] Z BERr
The optimal convergence rate p/n can be achieved by the least square estimator based on n uniformly selected observations.A key condition for the above result is that
|| []B(O(OTO)
Z)|l S en/N
for some constant 0.
The result of Pillacii and Wainwright [2016] can be applied to general estimators based on (O O yTy.. In this paper, however, we focus on the least square estimator based on (0 (O X, O y). Let Bo denote the solution to the sketched least square problem (2). In this paper, use the condi- tional mean squared error El 31 ZZ to measure the performance of Bo. The following theorem gives a lower bound for the conditional mean squared error of Bo.
Ip), the sketching matrix O is an N x n matrix with full column rank. Assume that is independent of E..... SN and with probability 1, O X has full column rank. Suppose N + co, P/N - O. Then N -
Theorem Suppose Assumption D holds, Z ~ (O,,, Ip),
E
{ 133 - 3| z 2(1 - op))) |BO
|)
p+1 N --
Theorem I gives explicit characterization of the impact of B Z) on the lower bound of E1 - 2 Z). Pilanci and Wainwright (016) showed that the condition (3) is satisfied by many classical sketching methods. Under the conditions of Theorem 1 for sketching methods satisfying the condition (3) the convergence rate of E(ll 3 Z) is lower bounded by p/n, which is the convergence rate for the uniform sampling method. Thus, order achieve a faster convergence rate than the uniform sampling method, the condition 33 should be violated.
Many existing sketching methods are through sampling the observations. For sampling methods, O is a column orthogonal matrix and each column of O has single nonzero element with value Hence ( 0) 'OT is a diagonal matrix whose diagonal elements are zeros and ones. For the IBOSS algorithm of Wang et al. [2019), the selected observations are completely determined by X and does not rely on additional randomness. Consequently || EOOOOT
'OT Z) 1. In this case, the lower bound provided by Theorem 1 has rate p/ / which is too loose. The following theorem gives tighter lower bound of the mean squared error for sampling methods.
|OOO O) 10-1000-0
Theorem Suppose Assumption 7 hollds 7 ~ (0p: Ip), the sketching matrix O is an N x n matrix with full column rank. Assume that O is independent of E EN and with probability 1. O X has full column rank. Furthermore, suppose E(O(O 0) O Z) = diag(d),.. dy)..
Let dmaax=
max.e de Then
E |30 B1 2
6n (p + log Nd + Op(n)) 0

---

If the matrix E(O(OTO) 0 Z d is diagonal, then
dmax
E(0 0) T Z) I B aER |||a|==1 a To =
aER [a
Under the conditions of Theorem the optimal convergence rate for sampling methods is lower bounded by
n log
Note that p < log(N/n), then the rate (4) is faster than the uniform sampling method. Theorem Section 3.2 will show that under certain conditions, the method of Wang et al. (2019) can achieve the optimal rate (
worth mentioning that Theorems D and are obtained under the condition Z; NN Ip). Perhaps, these results can be extended to the case that Z: has a general multivariate distribution. However, such results may not valid the distribution of Z₃ has heavier tail than normal distribution. In fact, our numerical results imply that faster convergence rate may be achieved when the distribution of Z; has heavy tail.
An estimator via data averaging
In this section, we would like to propose new sketching method which can hopefully have good statistical performance with low computational cost. To simple, when considering computation time, it is understood that the multiplication of an m X n matrix and n X matrix costs O (mmp) time, and the inversion of x p matrix costs (p²)) time. Note that the computation time ONN)) is essential if each observation is accessed at least once, e.g., to be loaded into the memory. The sketched least square problem (2) involves n reduced observations and the direct computation costs O(np² p²) computation time. The computation time O(p²) comes from the inversion of (p + 1) (p + 1) matrix which essential no matter how n is chosen. Hence the direct computation of any reasonable estimator which uses the information of the full data requires least O(Np p²) time. Thus, we restrict our attention to algorithms that can be completed within O(Np + p²) time.
To complete the computation within O(N++p') time, one needs to take n = O(N/p+p) Theorem impples that n c cap where 01, C2 > O are constants, then the optimal convergence rate for sampling methods reduces to p/n which is equal to the convergence rate for the uniform sampling method. Also, for large N. the reduced sample size n N/N/p + cep may still be large. To achieve faster convergence rate and better reduction of data, we would like to consider sketching methods other than sampling methods. This motivates us to propose a new data averaging method.
3.1 Methodology
Let F....., Ix C (I,...,NN N} be k mutually disjoint index sets, each containing indices, and .....,, N). To use the information of the full data, we assume N = kr. Let Z, = Eie ie M, Z; and D, = r-1 Eie ys be the averaged observation within the jthh index set, j =
U =1 I: (1,...,
1. k. It can be seen that
О 2 +
where E, = r-111 5, E.. Suppose that the choice of the index sets fI...., Ln is based on the covariates { and does not rely on the responses (ui), Then , Eg are mutually independent and are independent of (- Also, E, has mean 0 and variance o Thus, the averaged observations also satisfy the linear model and one can estimate B by the least square Xyūı), where X, r- "ie X., j 1,...,k. k. We would like to choose fI...., It such that f is a fast and accurate estimator of B. Let H (Z= () T)T and Z N Zi. The
estimator base on à reduced observations as P Xx,))
5

---

conditional mean squared error of 3 is
k
Σラ 2j
-1
ko?
F(1133 - 3||2 Z)
k
N
Σ 2j
Σ 2,2,
j=1
j=1
ko?
1 #
- -2 2HH
N
-H
H
0²
= (b (tr((t 2 H
N
In order to achieve a good statistical accuracy, we would like to choose the index sets such that tr(H--) + 2 H Z is minimized.
First we consider the simplest case of p = In this case, the matrix H reduces to a real number. Since / E IR², one needs at least two observations to estimate 3. To achieve maximum reduction of data, we take
2. Then takes its maximum when f1 {i E (I,.... N} = S =(N/2),1)
and f2 = {i E {I. [I,...,AN) .N. : 21,1 2 2(1/2+1)a)- The least square estimator of based on the averaged observations is 3 ((22) 2102)/(2, 21). (i2 31 (22 2)) The above estimator (ý2- 5/1 / Zn) of 3, is considered in Barton and Casley 1958) as a quick estimate of B,. which only considered the case of p = 1. To the best of our knowledge, no previous study generalized this estimator of Barton and Casley 1958) to the case p > 1.
For the general case of p 2 1, the exact minimizer of E(II@ - B|1 Z) may not be easy obtain. A simpler criterion to choose the index sets is to maximize the trace tr(H)= = E Zz 212.. This problem is equivalent to minimizing Σx -==1 CieA I - 1 and is instance of the balanced k- means clustering problem; see Lin et al. 2019) and the references therein. Unfortunately, algorithms for the k-means clustering problem are computationally intensive. In fact, for the vanilla k-means algorithm, each iteration takes PkN time which even exceeds the computing time of the least square estimator based on the full data. To achieve a balance between the statistical accuracy and the computing time, we deal with each variable in turn. We take k = 2p and for j ......... P, we determine two index sets, namely Lr,j and R based on the jth variable. Hence the set (1. N}
partitioned into 2p index sets 1r.1,..... Lr.p and R R. Rr.pr each containing r (2p) indices. The choice of these index sets is based on the following lower bound of tr (H).
N
N
tr
Σ 2.
1
Σ
25
+
Σ iE
2
ZJ
N
T
N
i=1
f=1
N
2
Σ
Σ
1
Σ
1
Zu
Zi,
+
Zi.
2h,
N
N
1=1
iE
1=1
2
2
$
max(2z
1
Σ
25.g
+
max Σ
PL. - 0):
r
V.
Zi.j. For j = I,...,,, Ps we choose Lys and Rr. to maximize the jth term of the above lower bound. Specifically, the first term of the above lower bound is {max((] - Eier., 2,1/r,0)) 0)}² + (max(X 2.1/r ž1, 0)}2, which takes its maximum when Lt.1 = EiE f1. [1....,N] N}: w} 2,1 I } and LT (ie {1....,N) 2,1 2 72,1 where 71.1 = 2(r),1 and 72.1 = 11-- After obtaining the index sets L, L....... L and E 2r,1,-..,- -1, we choose Lrj and Rr to maximize the jth term of the above lower bound, which is equivalent to maximizing {max(ž; - - Eies, zg//r, 0))² + {max(( R/a/r 2y, 0)}2, Hence we choose Lrg and R,, to be the indices of the remaining observations whose jtth variable is no larger than Mi.j and no less than 72.j respectively, where 71. and y2,5 are the rth smallest and the rth largest element of {=., : (1,..., NN
where 2 = -2r(p- - j 1)) "
R,,e))) respectivel. We average the observations within the groups Lr,,,.., Ltp and R. R Pr.p- Finally, we use the least square estimator based on the 2p averaged observations to estimate B. The proposed estimation procedure is summarized in Algorithm [

---

Algorithm 1: Data averaging algorithm
Input: Observations {Z,, 94 J.m.. covariate dimension p
Output: Estimator of /3
is assumed to be an integer
for E {1, .) do
n. + the rth smallest element of [z., : E (1, NJ
72.j f the rth largest element of [Z, : (1....,
Lrj - {i E AAL,..,,NN(U.. (Zn U..))) 1 124 5 ns}
Rr, + [i (I,....,N][(U
U.
:24 r2,)
7/2 +- r- r Eiee, 7, X = (1,7,
8 f- Liezy V/
2 rr EieR; Zi X = (1,23T)T,
8) (+-rr EiE V1
3A + Ej-1 X + E)
(En 83 + Ei X,
return BA
In Algorithm our strategy to select the index sets Lt. j and Lr. is closely related to the IBOSS algorithm of Wang et al. (2019). In fact, the index sets in Algorithm [ is exactly the index sets selected by IBOSS algorithm with subdata size == N. Of course, for IBOSS algorithm, taking n N unreasonable since the sample size not reduced. In fact, for IBOSS algorithm, one needs to take O + p) to complete the computation within O(Np + p²) time. Thus, the selection procedures of the proposed method and the IBOSS algorithm have different behavior. Theorem 4 will show that under certain conditions, IBOSS can achieve the optimal convergence rate (4) among all sampling methods. We shall see that the statistical performance of Algorithm is even better than the IBOSS algorithm.
Now we give an analysis of the computing time of Algorithm I Note that V1,j and Y2. are order statistics of no more than N elements. It is known that the selection of an order statistic among m elements can be completed within (m) time even the worst case; see Paterson (1996). Hence the computation of 7..,,.... ",,, and ,. 72.p can be completed within O(Np) time in total. It takes only one scan of the full data to compute the averaged observations, which takes O(N)) time.
Finally, the computation of B, based on 2p averaged observations can be completed within (pp)) time. In summary, Algorithm can be completed within O(Np + +p²) time and reduces the full data to merely 2p observations.
3.2 Asymptotic results
Now we investigate the asymptotic behavior of the conditional mean squared error of PA- In our asymptotic results, we treat p a function of N, and N tends to infinity. Let Z (21, (21....: 2p) ) be a random vector which is independent of Z and y and has the same distribution as Zz. The following theorem gives the exact limit of B 31 Z} when Z is a standard normal random vector.
holds, = N/(2p) is an integer, N > 2p², and Z Theorem Suppose that Assumption
(O Ip)- Also suppose that as N xc, + xo and p' (log(p)) log(N)) + 0. Then as N- 00,
p²o²
E (11@,, - 3| z = (1 + Op (1))
2 2log(2p)
Theorem iimlies that when is a standard normal random vector, the conditional mean squared error of BA has convergence rate p² /(log(2p) N). On the other hand, for sampling methods with N/ c2p for constants c1, C2 such that the computing time may be comparable. Theorem 2 implies that the optimal convergence rate of the conditional mean squared error is p² /N. In this view, the convergence rate of the proposed estimator is faster than sampling methods for p -> oc.
The proposed algorithm is closely related to the IBOSS algorithm. We would like to derive the asymptotic behavior of the conditional mean squared error of 3 Theorem 6(i) of Wang et al. 2019 gives asymptotic expression of the conditional covariance of (3, under the assumption that Z is
7

---

normally distributed. It is implied that if n and are fixed, then as N + 20,
( 113 - 3313 z = (1+ or))))
(l)g(N) tr( -1 diagi Σ)Σー + o
(5)
Now we derive the fine-grained limiting behavior of E: |1 - 3|| z for varying n and Let
Piá ozs (o1/(0 denote the correlation coefficient between and 2j. i = ,,..., Define
WN -aw a (1 AN aN diag( ー) E. -.
= p 2 log (N/r)'
We have the following theorem.
Theorem Suppose AssumptionPholle holds and Z ~ (M, r= is an integer, there exist
constants C1, C22 C3 O such that C₁ < > > < C2 and ||l| < Cz, there exists constant < 1/v2 such that maxigicisp Pig S p. there exist 61,62 (0, 1) such that for sufficiently large N, 4r I N S and
(1 + 262 1/ + [(q + 202) 1211112 < 1/2
(6)
Furthermore, suppose as N + oo,
0 and p² (log(n))
(7)
max (n, r log (N/r))
Then as N + 00.
E 11111- 2 -(1 (1+op(1)) (aN tr(W + aND W p + +1 o2
Remark If and p are fixed, then the condition (6) is sutisfied for sufficiently large N. On the other hund, if Pi.j = 0 for all S < j p. that is, the variables are independent, then the condition (6) becomes E1 Eg 1/2. In this case, the condition (6) holds for = O(N)) for some 0<E E 1/2.
Remark The condition (7) is satisfied if + and p = O(n2/2-c) for some > Also, the condition ( can be satisfied for arbitrary n, P provided N is sufficientll large.
Compared with Theorem 6()) Wang al. 2019]. our Theorem 4 gives more comprehensive characterization of the asymptotics of Var(B X). If μ = Op and ax - 0, then Theorem implics that
of
B 1111 B112 2 01
E
(r diag( (Σ)Σ
+
8
2 (N/r)
If we further assume that log(r)/ (og()))og())))))) + 0. then the expressions 5 and 8 are equivalent. However, theorem@iimllee that the expression S is not valid if aN does not converge to
Now we consider the special case that Z ~ (O,. I,). In this case, Theorem implies that if r/N - 0, n O(N*) for some 0 < < 1/2 and O(n-/2-e") for some 0 < < < 1/2, then
B B| Z} has convergence rate (a.p + 1)/n. We have
aN p 2ogglo222)) (p) + 2 (og(N/n)
p log(N/))
Hence if log(N/n) O then E (1a 1 z
Orlp" (n(p + log(N/n))) which
matches In this case, 3 achieves the optimal rate of sampling methods given by Theorem 2 and hence the optimal rate given by Theorem 2 tight.

---

4 Simulation results
In this section, we conduct simulations to examine the performance of the proposed estimator 3 A For comparison, the simulations also include the vanilla data averaging algorithm (abbreviated as VDA) where the full data is uniformly divided into A := 2p groups in random and the observations averaged within groups, the least square estimator based the uniform sampling method (abbreviated as UND, the leverage score sampling estimator (abbreviated as LEV) as described in Ma et al. 2015]] the sketched least square estimator based on the subsampled randomized Hadamard transform (abbreviated as SRHT), the estimator 31 produced by the IBOSS algorithm, and the least square estimator based on the full data (abbreviated as FULL.. The methods VDA, UNI, LEV, SRHT and IBOSS instances of the sketching framework (2) For these three methods, take n N/p. For SRHT, if N is power of 2, then the sketching matrix O (PHD) where P is an N matrix whose rows are uniformly sampled from the standard bases of IR H is an x N Walsh- Hadamard matrix (see, e.g., Dobriban and Liu (2119) and D is an N x N diagonal matrix whose diagonal elements are i.idd. Rademacher random variables; and N is not a power of 2, we pad zeros the original data make N reach power of 2. The computation of the proposed method and the IBOSS estimator rely on certain order statistics. For these two methods, the algorithm SELECT of Floyd and Rivest 1975) is used select the order statistics. The algorithms are implemented by C++. To be fair, for all algorithms, the estimators of B are solved by Gaussian elimination. The simulations are performed on a CPU with 3.30 GHz.
The statistical performance of an estimator 33 of / is evaluated by the empirical mean squared error based on 100 independent replications. Specifically, the empirical mean squared error is defined as ==
100-1>E3 100 |1@ |1@( 311 where 330) is the estimator in the ith replication. In all simulations, the
ground truth of B is vector with all elements equal to 1. We consider distributions of E1: the normal distribution E1 (0, 1) and the normalized chi-squared distribution E (x) We consider the following distributions of Z.
Case 1: {2y}5=1 are i..d. with uniform distribution Uniform 1).
Case 2: {=y} are i.i.d. with normal distribution (0,
Case { are i.i.d. with lognormal distribution, that is, log(z)) ~ N (0,
Case 4: { are i.i.d. with student distribution with 3 degrees of freedom 13.
Case 5: ~ (Op: Σ), where the diagonal elements of Σ all equals and the off diagonal elements all equals 0.5.
Case 6: Z is distributed as a mixture of (w, Σ) and (- Σ is defined as in Case 5, and the mixing proportions of the two component distributions are both 0.5.
p4 Σ) where y has all 1 entries,
Table 2 and Tables A.1-A.3 in Supplementary Material list the empirical mean squared errors of various estimators, where the proposed method referred as NEW. Among the implemented methods, VDA, UNI and SRHT are data-oblivious sketching methods while NEW, LEV and IBOSS are data-aware sketching methods. can be seen that VDA has the worst performance. This implies that the selection procedure is necessary for the proposed method. The simulation results show that UNI, SRHT. LEV have similar statistical performance. can be seen that the proposed estimator can achieve substantial improvement over the competing sketching methods. Especially, the proposed method shows superiority when p is large.
We also evaluate the computing time for various algorithms. Table 3 lists the computing time for Case with c ~ A (0. 1). Results for other settings are similar. It can be seen that the proposed method is slower than VDA and UNI. Compared with VDA and UNI, however, the proposed method has significantly better statistical performance and can achieve better data reduction. Compared with IBOSS, the proposed method has comparable computing time but a much better statistical performance. In summary, the new estimator has good performance in both speed and statistical performance.

---

Table 2: Empirical mean squared errors (multiplied by 10*) of various algorithms with N 8 10 and ε ~ N (0, 1).
p
NEW
VDA
UNI
SRHT
LEV
IBOSS
FULL
Case
50
182.438
16134
492.376
497.665
488.99
507.642
9.38737
100
651.596
14751.4
2149.19
2163
2114.55
2186.94
19.2125
200
2469.56
15014
13398.1
15788.1
14584.2
15159.2
38.0159
Case
50
6,95906
1039.08
33.6247
33,5017
32,9702
27.03
0.624944
100
21.9735
1039.08
144,342
144.579
141.629
125.218
1.28266
200
69.4913
1034.62
1033.37
1021.28
1035.9
926.762
2.54829
Case 3
50
5.99607
1008.38
28.2232
32.817
34.4454
13.635
0,465752
100
11,6466
998.573
105.39
126.361
168.775
40,3553
1,00904
200
23.0075
866.214
920.989
1003.37
982.843
235.444
1,95179
Case
50
2.00812
345.559
12.0477
11,2138
11.6838
2.36205
0.213812
100
4. 4.66736
346,259
53.6057
49.8285
50.3076
7.96878
0.439798
200
10.658
345.624
372.099
346,143
344.283
37.7542
0.856995
Case
50
19.8421
2127.66
65, 1897
63.4604
61.4904
59.3276
1.21801
100
62.2782
1922.4
284,473
281,714
275.586
264.896
2,4893
200
192.085
2015.87
2022.8
2031.46
1979.68
1973.69
5.00354
Case 6
50
19,5601
1987.22
62.558
63,5263
62.6909
55.7036
1,22561
100
60.9593
2022.88
282.951
286.522
280.046
263.008
2.47379
200
194.566
2007.68
1988.91
2022.94
2000.11
1938.25
5,09764
Table 3: Computing time (in seconds) of various algorithms.
N
P
NEW
VDA
UNI
SRHT
LEV
IBOSS
FULL
X 10*
50
0.05127
0,00302
0.00204
0.06567
0,12722
0,03279
0,05261
10a
100
0.10580
0.00635
0.00332
0.15607
0.52488
0.07127
0.22732
104
200
0.21261
0.01618
0.00653
0.32740
2.14122
0.16388
0.90082
6.4 X 105
100
1.32952
0.05715
0,03292
1.75414
4.51072
1.16810
1.95327
6.4 105
200
2.47827
0.10503
0.05132
3.31281
17.8339
2.31180
7.56397
6.4 105
400
4.88818
0.25092
0.12785
6.86445
85.0793
4.71242
36.3755
Conclusion
In this paper, presented new sketching method which is based on data averaging. The com- putation the proposed method can be completed within O(Np + y²) time. We proved that the proposed method can achieve faster convergence rate than sampling methods.
In the proposed algorithm, we need select certain order statistics of variables. This selection procedure is adapted from the IBOSS algorithm and thus allows us compare the performance of these two methods in fair manner. In theory, the selection order statistics can be completed within O(Np) time. However, this procedure may cost a lot of time in practice. It is interesting to investigate other selection procedures for data averaging. Also, this work focuses on the data averaging method for the linear model. It is interesting apply the data averaging method to regularized linear models and generalized linear models. We leave these topics for possible future research.
Supplementary Material
The Supplementary Material includes additional numerical results, all proofs and codes.
10

---

Acknowledgments and Disclosure of Funding
Wangli Xu is the corresponding author of this paper. The majority of this work was done when the first author was a postdoc in Renmin University of China. At the time the paper was accepted, the first author was with Inspur (Beijing) Electronic Information Industry Co., Ltd. This work supported by Beijing Natural Science Foundation (No Z200001), National Natural Science Foundation of China (No 11971478), and by Public Health & Disease Control and Prevention, Major Innovation & Planning Interdisciplinary Platform for the "Double-First Class" Initiative, Renmin University of China (No. 2023PDPC), and the MOE Project of Key Research Institute of Humanities and Social Sciences(No. 22JJD910001).
References
D. C. Ahfock, W. Astle, and S. Richardson. Statistical properties of sketching algorithms. Biometrika, 108(2):283-297, jul 2021.
D. E. Barton and D. J. Casley. A quick estimate of the regression coefficient. Biometrika, 45(3-4): 431-435, 1958.
Q. Cheng, H. Wang, and M. Yang. Information based optimal subdata selection for big data logistic regression. Journal of Statistical Planning and Inference, 209: 1 12-122, 2020.
E. Dobriban and S. Liu. Asymptotics for sketching in least squares regression. In H. Wallach, H. Larochelle, A. Beygelzimer, d' Alché-Buse E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
Drineas and M. W. Mahoney, RandNLA. Communications of the ACM, 59(6):80 90, may 2016.
R. W. Floyd and R. L. Rivest. Algorithm 489; The algorithm select-for finding the ith smallest of elements [ml]. Communications of the ACM, 18(3):173, Mar. 1975.
W. Lin, Z. He, and M. Xiao. Balanced clustering: A uniform model and fast algorithm. In Proceedings of the Twenty-Eigtth International Joiní Conference on Artificial Intelligence, IJCAI-19 pages 2987-2993. International Joint Conferences on Artificial Intelligence Organization, 2019.
P. Ma, M. W. Mahoney, and B. Yu. A statistical perspective algorithmic leveraging. Journal of Machine Learning Research, 16(27):861-91 2015.
Ma, X. Zhang, Xing. J. Ma. and M. Mahoney. Asymptotic analysis of sampling estimators for randomized numerical linear algebra algorithms. In Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, pages 1026-1035, Online, 2020.
M. W. Mahoney. Randomized algorithms for matrices and data. Foundations and Trends@ in Machine Learning, 3(2):123-22) 224,, 2010.
M. Paterson. Progress in selection. In Algorithm Theory SWAT'96, pages 368-379, 1996.
M. Pilanci and M. J. Wainwright. Iterative Hessian sketch: fast and accurate solution approximation for constrained least-squares. Journal of Machine Learning Research, 17:1-38, 2016.
G. Raskutti and M. W. Mahoney. A statistical perspective on randomized sketching for ordinary least-squares. Journal of Machine Learning Research, 17:1-31, 2016.
H. Wang. divide-and-conquu information- optimal subdata selection algorithm. Journal of Statistical Theory and Practice, 13(3):1-19, 2019.
H. Wang, M. Yang, and J. Stufken. Information- based optimal subdata selection for big data linear regression. Journal of the American Statistical Association, 114(525):393 405, 2019.
J. Wang, J. D. Lee, M. Mahdavi, M. Kolar, and N. Srebro. Sketching meets random projection in the dual: A provable recovery algorithm for big and high-dimensional data. Electronic Journal of Statistics, 1(2):489649944, 2017.
D. P. Woodruff. Sketching as a tool for numerical linear algebra. Foundations and Trends in Theoretical Computer Science, (0(1-))iv+157 2014.


Reviews:
[{"review_id": "qMIZn95m1u", "reviewer": "Reviewer_neFN", "summary": "This paper studies linear regression where the number of samples N is much larger than the number of predictors p (N>>p), which is computationally costly due to large N. The paper investigates lower bounds for existing sampling based methods, and proposes a novel sketching method based on data averaging which reduces the original data to a few averaged observations. Theoretical results show that the method has a faster convergence rate than the optimal convergence rate for sampling methods, and experiments show that the method reduces mean squared error over previous methods.", "questions": "Questions:\n\n- How does the memory usage compare for competing methods? Please add the CPU memory usage either in Table 5, or in a new table. \n\nSuggestions:\n\n- Minor grammatical and typographical errors should be fixed with proofreading. The citation style also does not match the Neurips style which should be fixed.\n\n- One example for the poor presentation is that the matrix $A$ is defined after it is referenced in the text as $\\beta_A$. Fixing such issues would help the paper.\n\n- The presentation of the paper can be improved by deferring some of the related work in page 2, and details of the derivations for the theorems to the supplemental material. The generated space can be used to add numerical results, which could include some more, real data examples.", "limitations": "There are no clear potential negative societal impacts of the work.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "- The proposed method is a new sketching method for large-scale linear regression with large N, and it is original.\n\n- The method is well-motivated from a theoretical perspective, and the lower-bounds for the proposed method under certain assumptions are better than recent sampling based methods.\n\n- The method has $\\mathcal{O}(Np + p^3)$ which is favorable when $N>>p$.", "weaknesses": "- There are no results in the main paper, which are quite important for a sketching paper where claims are improved computing time with reduced mean squared error. The results which are in the supplemental are not extensive, and there is only one real data example, which makes it difficult to quantify the benefits of the proposed method.\n- The setup of the method is constrained to the setting with no regularization, and only the setting where $N>>p$. It might be good to add a discussion on how applicable this method is to other settings for linear models with different N and p, and with regularization.\n- The presentation of the paper is quite difficult to follow.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely appreciate your valuable comments and suggestions, particularly the positive feedback regarding the proposed method. \n\nBelow, we will provide a response that centers around three main aspects: **weaknesses, questions, and suggestions** related to our paper. These relevant sections have been indicated with italicized text in the following. We have noticed that your feedback primarily focuses on the presentation of the paper. We believe and hope that our revisions will ultimately meet your satisfaction, and that the paper will have the opportunity to be published in NeurIPS.\n\n\n**Weaknesses:**  \n\n  * *There are no results in the main paper...*  \n  \n     **Reply:**  Thank you for your suggestions. In the original paper, experimental results were deferred to the appendix due to the space constraint. We will consider including some experimental results in the main text.\n \n     The present paper is mainly a theory/exploratory paper. The merit of the newly proposed method is that it break the ice and proves that there indeed exists a method that can break the lower bound in Theorem 2. This proves: better is possible. This fact may be much important than the practical performance of the proposed method. While the proposed method may be far from optimal, it may shed light on the future development.\n\n  * *The setup of the method is constrained to the setting with no regularization...*\n\n    **Reply:**\n    Our results carries particular significance in the case of $N>>p$, where the proposed estimation method in the paper exhibits a faster convergence rate.\n    You give two interesting direction: regression with regularization and the case beyond $N>>p$.\n    These directions may be highly nontrivial.\n    We will consider add brief discussions on these topics in the paper.\n    \n    For regularization, a straightforward method is to use the proposed method to reduce the observations to $2p$, and use the $2p$ observations to perform regularized regression.\n    One may immediately obtain some theoretical results since the reduced $2p$ observations also satisfy the linear regression model.\n    But much effort may be paid to seriously explore the fine-grained behavior of such approach.\n    \n    For the case beyond $N>>p$, the methodology may be entirely different since in this case, reduction of $N$ may not be a good strategy.\n    In fact, perhaps most existing work on sketching methods for least squared problem can only work well for the case $N>>p$ since the core idea is reducing $N$.\n\n    Overall, these two directions may be good topics to consider in the future.\n  \n  * *The presentation of the paper is quite difficult to follow.*\n\n     **Reply:**  Thank you for your comments. We will improve the presentation of our paper based on your following suggestions. Specifically, the presentation of the paper can be improved by deferring some of the related work in page 2, some theoretical results, and some details of the derivations for the theorems to the supplemental material. The saved space can be used to show numerical results.\n\n\n\n\n**Questions:** \n\n\n  * *How does the memory usage compare for competing methods? Please add the CPU memory usage either in Table 5, or in a new table.*\n\n    **Reply:** \n    Now we report the memory usage results.\n     The setting is as follows: the CPU frequency is 3.1 GHz; the memory is 16 GiB; the operation system is Ubuntu 22.04.2; the compiler is  gcc 11.4.0.\n     The data type is double.\n     We record the maximum resident set size (maxrss) used (in kilobytes) which is obtained via the function getrusage in sys/resource.h.\n     Since the result is slightly different in each run, we report the average result of 10 runs.\n\n    Note that the data itself takes $N \\times (p+1) * 8$ bytes of memory.\n    For the case of $N = 8 \\times 10^4$ and $p = 50$, the data takes 31875 KiB memory.\n    For the case of $N = 6.4 \\times 10^5$ and $p = 400$, the data takes 2005000 KiB memory.\n    It can be seen that except for SRHT, the memory cost is mainly for storing the data.\n    For SRHT, we need to transform the data, so additional memory is used.\n    In our experiments, the full data is in the main memory.\n    The present work is mainly a theoretical work, and is confined to this case.\n    One interesting scenario is that the full data is too large to fit in the main memory. \n    This is a topic worth exploring in the future.\n\n\n\n    $N$ | $p$ | NEW | VDA | UNI | SRHT | LEV | IBOSS | FULL |\n     -------- | -------- | -------- | -------- | -------- | -------- | -------- | -------- | -------- |\n    $8 \\times 10^4$ | $50$ | 36242 |   36516 |  35989.6 | 89050.8 | 36839.2 | 36260  | 35456 |\n    $6.4 \\times 10^5$ | $400$ | 2017482 | 2019960.4 |  2012389.6 |  5300062.8  | 2025016 | 2015057.6  | 2009964 |  \n\n\n\n\n**Suggestions:**\n\n\n  * *Minor grammatical and typographical errors should be fixed with proofreading. The citation style also does not match the Neurpis style which should be fixed.*\n\n      **Reply:** \n      Sorry for that!\n      We will address the presentation issue carefully.\n      Thank you for pointing out the issue of citation style.\n      We will fix it in the camera-ready version.\n    \n      \n   * *One example for the poor presentation is that the matrix $A$ is defined after it is referenced in the text as $\\beta\\_A$. Fixing such issues would help the paper.*\n\n       **Reply:** \n       Thank you for pointing out the confusion about notation.\n       The subscript of $\\beta\\_A$ was intended to be the initial letter of \"average\".\n       We will consider improving the notations in the camera-ready version.\n      \n  * *The presentation of the paper can be improved...*\n\n      **Reply:** Thank you for providing detailed suggestions on the paper's presentation. We will improve the presentation in the camera-ready version."}, {"review_id": "K2rx7XUCy8", "reviewer": "Reviewer_EUfX", "summary": "This submission studies the asymptotic estimation risk of linear regression \nwith various data sketching strategies. \nThe authors start by refining an existing lower bound to show that sketching\nby uniform sub-sampling is only minimax optimal when the feature dimension\nis very large. otherwise, improvement is possible and the authors prove the\npreviously proposed IBOSS algorithm is minimax optimal. \nIn order to circumvent their own lower-bound, the authors then exit the \nrandom sketching framework to define an new method based on data averaging.\nThis method is shown to outperform uniform random sketching and IBOSS\nboth theoretically and empirically.\n", "questions": "- Assumption 1: I don't see where $r > 0$ is used in the assumption.\n\n- Theorem 1/2: How restrictive is the requirement that $Z$ be normally distributed? \n    Does the result of Pilanci and Wainright (2016) require similar\n    distributional assumptions on $Z$ in order to derive minimax lower bounds?\n\n- Theorem 2: Why does $\\\\mathbf{E}[O (O^\\\\top O)^{-1} O^\\\\top | Z]$ being diagonal\n    allow for a tighter lower bound to be derived? \n    In particular, is there a chance that sketching matrices not satisfying \n    this condition can out-perform the lower bound in Theorem 2?\n\n- Displays after lines 187 and 208: Where do these expressions come from? \n    The lower bound on the trace after line 208 is used to guide the development\n    of Algorithm 1. Since it plays a critical role, some intuition into where\n    this expression comes from is is important for the reader to follow\n    the argument.\n\n- Theorem 3: It is somewhat strange that $r = N / 2p$ is said to be an integer when\n    the conditions on $p$ and $N$ require that $r \\\\rightarrow \\\\infty$. \n    I suppose what you really mean is that the divergent sequences $N_i$, $p_i$ \n    are chosen so that $r_i = N_i / 2p_i$ is always integer?\n    Is this condition simply for convenience when working with Algorithm, i.e.\n    to avoid working with floor operators? \n\nRemark 1: I think some additional comments on the tail bounds in (5)-(7)\n    must be provided. Firstly, what is $\\\\mathcal{A}$? I only see the definition\n    of $\\\\mathcal{A}_j$ in lines 236/237. Is this any measurable set? \n    Secondly, what to do the tail bounds imply about the original dataset $Z$? \n    It seems that they say something like \"tail examples concentrate \n    faster as dimension increases,\" but it's hard to tell since \n    the $\\\\mathcal{A}$ sets also depend on $p$. \n    Do these conditions make sense for general data? Can you give an example\n    of a data-generating distribution for which these conditions hold, e.g.\n    sub-Gaussian distributions? \n\nTheorem 4: So, essentially $z$ follows a symmetric, isotropic distribution \n    with strict conditions on the tails?\n    Equation 9 implies that the tails need to be sufficiently large, but \n    not too large (Equation 10); how do these conditions compare to (5)-(7)? \n    Again, how they are satisfied by real data?\n    I cannot understand the impact of Theorems 3/4 without understanding\n    how and when these tail conditions are satisfied.\n\nTheorem 5: Logarithmic improvement in the convergence when the number of \n    features tends to infinity seems quite weak. Asymptotic regimes where \n    $N \\\\rightarrow \\\\infty$ are somewhat justified since more data can be \n    collected, but rarely can more features be collected.\n    Without that said, I do like that the $\\\\log(N / n) \\\\ll p$ condition \n    is not required, since this implies even more unrealistic asymptotics\n    for $p$.\n    Additionally, why not just assume $Z$ is standard normal from the start,\n    rather than developing Theorems 3, 4 under awkward tail assumptions?\n    Theorem 5 is much cleaner and Theorems 1, 2 already require the \n    standard normal assumption.\n\nTheorem 6: It would be nice to remind the reader the difference between\n    $\\\\beta_I$ and $\\\\beta_A$ here. I think $\\\\beta_I$ is the IBOSS algorithm,\n    but this should probably be stated in the Theorem.\n\nTables 1/2: It seems like every method but the one proposed experiences catastrophic failure \n    when $p = 200$, performing about as well as the naive VDA method.\n    IBOSS performs slightly better than competitors in Cases 3,4, but otherwise\n    it also fails.\n    I have two questions here: (i) Why do all other methods approach the same\n    MSE for $p = 200$? \n    Is this the error of an zeroth-order estimator, e.g. the mean? \n    (ii) The baseline methods are all sketching-based, meaning Theorem 2 would\n    lower-bound their risk, correct? \n    If so, does their poor performance reflect this lower \n    bound and the good performance of the proposed method relies on circumventing\n    it by exiting the sketching framework?\n\nLine 83 (Appendices): Using commas to separate methods and commas to denote\n    the thousands place makes this list difficult to read. I suggest putting the \n    results in a standard table with one row for improved clarity.\n", "limitations": "As mentioned in the \"Weaknesses\" section, the major limitation of this \npaper is the presentation and quality of the writing. Specifically,\n\n- The assumptions in Theorems 3 and 4 need to be clarified before their \n    utility is clear to me. Since Theorem 5 makes a standard normal assumption \n    just as Theorem 2 does, it's not clear to me why it is interesting to present\n    Theorems 3/4 in the main paper when they require specific and unjustified tail bounds\n    on the data distribution.\n\n- The explanation of Algorithm should be improved. I suggest including a figure\n    illustrating Algorithm 1, if possible.\n\n- The rates in Theorems 5 and 6 should be compared and some commentary given.\n    In general, more explanation of the theorems ought to be given. \n\n- I strongly suggest reducing the technicality of the main paper and introducing\n    some of the empirical results from the appendix to clarify the utility\n    of the proposed averaging method. The method works well in practice,\n    which should be highlighted, rather than hidden.\n", "rating": 7, "confidence": 2, "soundness": 4, "presentation": 2, "contribution": 3, "strengths": "This submission has two main strengths:\n\n- It completes the landscape of minimax lower bounds for random sketching\n    applied to linear regression, establishing optimality of IBOSS\n    in a new regime where the number of features is not too large. \n\n- It improves on the performance of sketching methods by developing a new framework\n    based on data averaging. This framework is novel, achieves faster\n    (theoretical) convergence of the estimation risk, and outperforms sketching\n    methods in practice. It is also comparable to IBOSS in computation time. \n\nI cannot comment on the correctness of the theoretical developments or the \nnovelty of the proofs since this paper is out of my research area. Related\nwork seems appropriately referenced.\n", "weaknesses": "The main weakness of this paper is the presentation, which confuses the novelty\nof the results and makes much of the theoretical development hard to follow.\nIn particular,\n\n- The discussion is dense, mathematical, and does not give any intuition into the\n    bounds proved or and the assumptions required. This makes it difficult to\n    evaluate the novelty and utility of some results, in particular Theorems 3 and 4.\n    No effort is made to compare the risk bounds given throughout the paper,\n    although it would be very useful for Theorems 5 and 6.\n\n- The experimental results, which make a strong argument for the utility of the \n    proposed data averaging method, are all deferred to the appendix.\n\n- The description of the data averaging method in Section 3.1a is difficult\n    follow so that Algorithm 1 is more useful than the text itself. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely appreciate your valuable comments and suggestions. Due to space constraints, we display the discussion about non-normal cases, presentation and tables in the global rebuttal. Below, we will provide a response for main **questions and limitations**. \n\n**Questions:**  \n* *Assumption 1...*  \n**Reply:** In the latter part of the text, we define $r = \\frac N{2p}$, and the condition $r>0$ can mathematically avoid the case $N<2p$. We would like to delete this condition and make clear where necessary that the case $N<2p$ is not considered.\n\n* *Theorem 1/2...*  \n**Reply:**  The response to the first question can be seen in the global rebuttal.\nAs for the second question, the minimax result of Pilanci and Wainright (2016) does not assume $Z$ is random, and hence is more general indeed. However, as we pointed out in the main text, their restriction on $\\mathbf{E}_o = \\mathbf{E}\\left[O(O^{\\top} O)^{-1}O^T|Z\\right]$ is restricted, which has a different focus. As a result, the corresponding conditions for establishing the lower bound differ.\n    \n* *Theorem 2...*  \n**Reply:**  The diagonal condition of $\\mathbf{E}_o$ restricts the class of legal sketching methods. And it is possible to obtain fine-grained result for a restricted class of methods. Technically, in the proof of Theorem 2, the diagonal condition is used to obtain \nthat $\\text{tr}[\\textbf{X}^T\\mathbf{E}_o\\textbf{X}]=\\sum\\_{i=1}^Nd_i||X\\_i||^2$.  \nThis equality does not generally hold if $\\mathbf{E}_o$ is not diagonal. So there is a chance that sketching matrices not satisfying this condition can out-perform the lower bound in Theorem 2. Indeed, the proposed data averaging method is such a method, where corresponding $\\mathbf{E}_o$ is not diagonal in general. According to Theorem 5, the proposed method can outperform the lower bound established in Theorem 2. So this is the merit of the proposed method, its existence proves it is possible to go beyond the lower bound of Theorem 2.\n\n* *Displays after lines 187 and 208...*  \n**Reply:** The details for formula derivation can be seen in the PDF file attached to the global rebuttal.\n\n* *Theorem 3...*  \n**Reply:**  Yes, what we really mean is that the divergent sequences $N_i,p_i$ are chosen so that $r_i=N_i/2p_i$ is always integer for convenience. We will clarify this in the main text.\n\n* *Remark 1...*  \n**Reply:** Yes, $\\mathcal{A}$ is a any measurable set in $p$-dimensional Euclidean space, and $\\{\\mathcal{A}: \\Pr(Z\\in\\mathcal A)\\leq\\frac 1p\\}$ represents the class of measurable set $\\mathcal{A}$ satisfying $\\Pr(Z\\in\\mathcal A)\\leq\\frac 1p$.  \nRoughly speaking, tail bounds (5)-(7) require that the dimension $p$ is not too large and/or the tail of $Z$ is sufficiently light. They are are satisfied by standard Gaussian distribution under some conditions on dimension, see lines 429-437 in Appendix. These bounds may be satisfied by more general distributions, sub-Gaussian distributions, but to establish such results may require a lot of mathematical work. We would like to follow the suggestion and confine ourselves to the Gaussian setting. \n\n* *Theorem 4...*  \n**Reply:**  Equation 9 is a weak condition. Roughly speaking, Equation 9 merely say that no $z_j$ will shrink to $\\mu_j$. Equation 10 is another bound for the tail. So the conditions (5)-(7) and (10) all require the tail to be not too large. However, there seems no simple condition to unify these three conditions. Under certain conditions on dimension, these conditions are satisfied by standard normal distribution. We would like to confine ourselves to the normal distribution setting.\n\n\n* *Theorem 5...*  \n**Reply:**  We agree that the logarithmic improvement is quit weak. But the merit is that it proves that the improvement is possible. Now we know there is indeed some methods with good computing time, other than samplings, may break the lower bound of Theorem 2. So if the theoretical contribution is taken into account, the present work may be meaningful.  \nPreviously, one of our goal is to try to make the results general and go beyond normal distribution. It turns out that the cost is greatly reduced readability. \n* *Theorem 6...*  \n**Reply:** Yes, $\\beta_I$ is the IBOSS algorithm.  \n* *Tables 1/2...*  \n**Reply:** For question (1): For competing sketching methods, we take $n=N/p$ so that the typical computing time is $O(Np+p^3)$. In this case, the error rate is of order $p^2/N$. In Tables 1 and 2, $N=8\\times10^4$, $p=200$ and $n=400=2p$. For VDA method, all data are reduced to $2p$ averaged data. In this view, the result in Tables 1 and 2 is reasonable.  \nFor question (2): The baseline methods would be either bounded by Theorem 1 or Theorem 2. Yes, the poor performance reflect the lower bounds of Theorems 1 and 2. The newly proposed method can also be treated as a sketching method, but it violates the key conditions of Theorems 1 and 2. So it is possible to break the lower bounds of Theorem 1 and 2, resulting better performance.\n  \n**Limitations:**\n\n* *The rates...*  \n**Reply:**   Theorem 6 implies that, under some technical conditions, IBOSS achieves the lower bound given by Theorem 2, and hence the bound of Theorem 2 is tight. One may not expect a sampling method significantly better than IBOSS. For the typical implementation of sampling methods, one needs to take $n\\asymp N/p$ to make the computing time within $O(Np+p^3)$. Theorem 5 implies that, under some technical conditions, the proposed method converges faster than sampling methods with comparable computing time, including IBOSS. Roughly speaking, the proposed method breaks the lower bound in Theorem 2. Hence better is possible.\n"}, {"review_id": "FFO2thB3Y7", "reviewer": "Reviewer_DaqC", "summary": "This paper gives lower bounds of the conditional mean squared error for sketching methods. They focus on least square estimator, and show that when the problem dimension is sufficiently large, the optimal error rate among all sampling reductions is achieved by uniform sampling. They also propose a sketching method based on data averaging that achieves better performance under some scenarios. Fast implementation of their proposal is also given in the paper.  ", "questions": "1. Line 80 \"ruduce -> reduce\"\n2. In the equation between line 123 and 124, X stands for the training data or the testing data? \n3. Why focus on mean squared error instead of prediction error? The latter objective seems more interesting to me. If the authors have specific reason for that, it would be helpful to state it in the paper. \n4. The Gaussian assumption is a bit restrictive. Does the theory work for other light-tail distributions, for example sub-Gaussian? \n5. I might have missed this, but where is $\\hat\\beta_I$ in Theorem 6 defined? \n", "limitations": "Yes. Societal impact not applicable. ", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "This paper gives a careful analysis on the lower bound of estimation error caused by sampling methods, which is tight in many cases. They also proposes a simple yet useful sketching algorithm that has better statistical guarantees comparing to sampling based algorithms in many interesting regimes. They make the comment that the associated optimization problem could be hard to solve, and they come up with a reasonable alternative approach. This work gives nice theoretical insights on data averaging algorithms, which could help with the design of data reduction methods in practice. ", "weaknesses": "This paper considers solely least square estimators after sketching, which could restrict its applicability in comparison to results that involve more general estimators. In the paper, they proposed a sketching method based on data averaging that is supposed to be more efficient and achieves better accuracy. However, they didn't compare their proposal with other sketching algorithms other than sampling based methods.  ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "\n**Comment:**  We sincerely appreciate your valuable comments and suggestions, particularly the positive feedback regarding the proposed method. \n\nBelow, we will provide a response that centers around two main aspects: **weaknesses, questions** related to our paper. These relevant sections have been indicated with italicized text in the following.  We believe and hope that our revisions will ultimately meet your satisfaction, and that the paper will have the opportunity to be published in NeurIPS.\n\n\n**Weaknesses:**  \n\n* *This paper considers solely least square estimators after sketching, which could restrict its applicability in comparison to results that involve more general estimators. In the paper, they proposed a sketching method based on data averaging that is supposed to be more efficient and achieves better accuracy. However, they didn't compare their proposal with other sketching algorithms other than sampling based methods.*\n\n    **Reply:** As you pointed out, in this work, we confine ourselves to the least square estimators after sketching, and do not consider some more general estimators such as iteration methods.\n    There are two benefits to consider such a class of estimators.\n    \n    First, the considered class, although restricted, is important.\n    Indeed, within the considered class of estimators, there are not only some recent sampling and sketching methods, but also some new methods such as the proposed one.\n    Even for iteration methods, the initial point is often chosen as a sketched least square estimator which is in the considered class.\n    \n    Second, the considered class allows for theoretical analysis.\n    The considered class of estimators is not too large, so it is possible prove some theorems to understand certain properties of the class, just like we did in Section 2.\n    Such results may be very difficult to obtain if one considers a larger class including more general estimators such as iterations methods.\n    Our theoretical analysis does provide insights: it shows that sampling methods may be sub-optimal.\n    \n    Indeed, in the supplementary material, we provided comparisons with the sketched least square estimator based on subsampled randomized Hadamard transform (SRHT).\n    For sketched least square estimators with other sketchings, \n    the convergence property may have a similar behavior in view of the lower bound of Pilanci and Wainwright (2016).\n\n\n**Questions:**  \n\n\n\n\n  1. *Line 80 ''ruduce $\\rightarrow$ reduce\"*\n  \n      **Reply:** Thank you for pointing out the typo, and we  correct it in the updated version.\n\n  2. *In the equation between line 123 and 124, $X$ stands for the training data or the testing data?*\n\n      **Reply:** Thank you for raising this question. In lines 123-124, the symbol $X$ refers to the training data set.\n\n\n  3. *Why focus on mean squared error instead of prediction error? The latter objective seems more interesting to me. If the authors have specific reason for that, it would be helpful to state it in the paper.*\n\n      **Reply:** Mean squared error (MSE) and prediction error (within training data) can both measure the performance of an estimator.\n      MSE may be more convenient in our setting.\n      In fact, prediction error involves the data matrix $\\textbf{X}$ which is random.\n      There is also randomness in $\\hat \\beta\\_{A}$.\n      These two randomness are not independent.\n      The two dependent sources of randomness make the analysis of prediction error complicated.\n      For some other methods, either $\\textbf{X}$ is assumed to be constant or $\\textbf{X} (\\hat \\beta\\_A - \\beta)$ can be simplified by explicit formula of $\\hat \\beta\\_A$.\n      For out method,  $\\textbf{X}$ is random, and $\\hat \\beta\\_A$ does not have a simple formula. So we choose MSE for convenience.\n\n\n  \n  4. *The Gaussian assumption is a bit restrictive. Does the theory work for other light-tail distributions, for example sub-Gaussian?*\n\n      **Reply:**   Theorem 1 may be generalized to non-Gaussian case.\n      We do not think Theorem 2 can be easily generalized to non-Gaussian distributions since in that case, one may need to deal with the concentration inequality of more general order statistics. Theorems 3 and 4 do not assumes $\\textbf{Z}$ Gaussian. They imposes some conditions on the tail probability of $\\textbf{Z}$.\n      But the imposed conditions may not be easily verified for sub-Gaussian data.\n      Theorems 5 and 6 rely on Gaussian assumptions. There are some $\\log$ in the statement of Theorem 2 and 3. These $\\log$ comes from Gaussian data.\n      For sub-Gaussian data, the $\\log$ term may not be correct.\n  \n  5. *I might have missed this, but where is $\\hat{\\beta}\\_I$ in Theorem 6 defined?*\n\n      **Reply:** $\\hat{\\beta}_I$ denotes the estimator of IBOSS.\n      And we add a diagram which may demonstrate IBOSS more clearly.\n\n**Reference**\n\nMert Pilanci and Martin J. Wainwright. Iterative Hessian sketch: fast and accurate solution approximation for constrained least-squares. Journal of Machine Learning Research, 17: 1–38, 2016."}, {"review_id": "365enWGbqY", "reviewer": "Reviewer_sefb", "summary": "This paper considers a new estimation method for a large scale linear regression model. Specifically, the regression coefficients are estimated by least squares estimation of averaged observations for which data are partitioned via a method similar to the information-based optimal subdata selection (IBOSS) algorithm proposed in Wang et al. [2019]. The paper develops both lower and upper bounds on the mean squared error of the estimator and compares it with the existing methods theoretically as well as numerically.  ", "questions": "- The expression of \"fine-grained lower bounds\" in the abstract is a bit misleading because regressors are assumed to be jointly normal with mean zero vector and an identity covariance matrix. It might be better to say this limitation more explicitly in the abstract.\n- Line 115: $r > 0$ is introduced but it seems that it is not defined before.\n- The paper focuses on the normal regressor case. Would it be possible to extend to sub-Gaussian or other more general cases?\n- In view of the remarks after Theorem 5, the proposed method has the convergence rate $p^2/(\\log(2p)N)$ for the mean squared error; whereas, the convergence rate for the sampling methods is $p^2/N$. Thus, the difference is only by $1/\\log(2p)$, which may not matter much when $p$ is relatively small. The abstract states that \"our theoretical results show that the proposed method can achieve a faster convergence rate than the optimal convergence rate for sampling methods.\" However, the difference is only at the log rate and so it might be better to explicitly state that the faster rate is only up to the log factor. \n- On the other hand, if we look at the numerical results in the supplement, there are huge differences between the proposed method and existing methods in terms of mean squared errors. It would be good to know where this large difference comes from. For example, is this because the subsample size $n = N/p$ in the existing methods is chosen to be too small? Some clarification would be helpful.\n     ", "limitations": "- One limitation is that the paper does not provide any guidance regarding how to conduct statistical inference, e.g., construction of confidence intervals. It would be helpful to provide a method for inference if possible and difficulties if not readily available. \n- The supplementary material does not include replication files. It would be desirable to provide them if the paper is accepted. ", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "- The paper proposes a new estimation method for linear regression, which is arguably one of the most important estimation problems in statistics and machine learning, although the new method shares lots of similarities with the information-based optimal subdata selection (IBOSS) algorithm proposed in Wang et al. [2019].\n- The paper contains a number of interesting theoretical results.   ", "weaknesses": "- The main text of the paper does not contain any experimental results, although numerical results in the supplement are quite impressive. It would benefit the readers if the paper contains a concise summary of key numerical results in the main text.\n- It is very important to compare the proposed method with the information-based optimal subdata selection (IBOSS) algorithm proposed in Wang et al. [2019] because data selection is very similar between the two methods. The paper includes lots of remarks on the IBOSS algorithm; nonetheless, it is still unclear exactly what sense they are different. It seems to me that the averaging method in the paper induces different type of weights across included observations relative to the IBOSS algorithm, but the exact link between the two methods is illusive.\n- On page 25 in the supplement, a different algorithm is introduced. There are some possible typos: (i) a variant of Algorithm ??, which I presume refers to the algorithm in the main text; (ii) $r \\leftarrow\\left\\lfloor\\frac{n}{2 r}\\right\\rfloor$: perhaps $2p$ in the denominator. Please check them. This variant of Algorithm ?? looks much simpler but it is only introduced in the supplement. It might be useful to include this in the main text to improve the understanding of the proposed method.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely appreciate your valuable comments and suggestions, particularly the positive feedback regarding the proposed method. \n\nBelow, we will provide a response that centers around three main aspects:  **weaknesses, questions, and limitations** related to our paper. We have taken note that your feedback primarily emphasizes the presentation and explanation of the paper. We believe and hope that our revisions will ultimately meet your satisfaction, and that the paper will have the opportunity to be published in NeurIPS.\n\n**Weaknesses:** \n\n  * *The main text...*\n  \n    **Reply:**  Thank you for your suggestions. Due to space constraints, all experimental results were deferred to the appendix in the original paper. In response to your feedback, we will address this by moving some content from the main paper to the appendix. If space permits, we will include Table 1 and Table 5 in the main text and include a concise summary of key experimental results.\n\n  * *It is...*\n\n    **Reply:**  Thank you for your suggestion. As you noted, the proposed method has a similar selection method as IBOSS. However, the two algorithms has essential differences which results different convergence rate. Below we elaborate the key differences. The discussion in our paper will be improved in camera-ready version.\n    \n     1. The selected observations of IBOSS is only a subset of the full data. This is an inherent limitation method of sampling methods: sampling methods only select a (typically small) subsample of the full data, other samples are discarded. In comparison, the proposed averaging method selects all samples: it clusters the full data into $2p$ clusters (via the IBOSS way), and averaging within each group. Hence the averaged samples are computed from the full data, no observation is discarded.\n    \n    2. The proposed method and IBOSS has different convergence rate. Theorem 6 implies that under certain technical conditions, the convergence rate of IBOSS can reach the lower bound in Theorem 2. On the other hand, Theorem 5 implies that under certain technical conditions, the convergence rate of the proposed method can break the ice and be lower than the lower bound in Theorem 2.\n    \n    3. The proposed method and IBOSS have different degree of data reduction.\n    IBOSS reduces $N$ samples to $n$ samples where $n\\asymp N / p$ to achieve $O(Np + p^3)$ computing time.\n    In comparison, the proposed method reduces $N$ samples to merely $2p$ observations.\n\n  * *On page 25...*\n\n    **Reply:** We deeply apologize for these typos, and will improve the presentation of our paper. The ''Algorithm ??\" should be  IBOSS algorithms. The  denominator of $\\frac{n}{2r}$: is  $2p$ instead of $2r$. This algorithm is coupled with IBOSS algorithm to prove Theorem 6. Since it is not related to the proposed algorithm, we decide to not put it in the main text.\n\n**Questions:**  \n\n  * *The expression...*\n  \n    **Reply:** Thank you for your reminder. We have  addressed this limitation in the revised paper by clarifying it in the abstract\n\n  * *Line 115...*\n\n    **Reply:** Sorry, $r = \\frac{N}{2p}$ is defined latter. We will fix this mistake.\n\n  * *The paper...*\n\n    **Reply:** Theorem 1 may be generalized to non-Gaussian case.\n    We do not think Theorem 2 can be easily generalized to non-Gaussian distributions since in that case, one may need to deal with the concentration inequality of more general order statistics. Theorems 3 and 4 do not assumes $\\textbf{Z}$ Gaussian. They imposes some conditions on the tail probability of $\\textbf{Z}$. But the imposed conditions may not be easily verified for sub-Gaussian data. Theorems 5 and 6 rely on Gaussian assumptions. There are some $\\log$ in the statement of Theorem 2 and 3. These $\\log$ comes from Gaussian data. For sub-Gaussian data, the $\\log$ term may not be correct.\n\n  * *In view...*\n\n    **Reply:** Following your advice, we will make clear that our theoretical results demonstrate that the proposed method can achieve a faster convergence rate, up to a $\\log(p)$ factor, compared with the optimal convergence rate of sampling methods.\n    \n    We would like argue a bit that the $\\log(p)$ term may not be as weak as it appear.\n    In fact, from our theoretical results, the $\\log(p)$ improvement can already break the lower bound of sampling methods stated in Theorem 2.\n\n  * *On...*\n\n    **Reply:**  In fact, for large $p$, the term $2\\log(2p)$ in Theorem 5 may not be a small number. For example, if $p = 200$, then $2\\log(2p)$ is approximately 12, which may look \"huge\". In this view, the good performance of the proposed method is reasonable.\n    \n    Our experimental results show that compared with UNI, the proposed method has particularly good performance in Case 3 and Case 4.\n    Note that in Case 3 and Case 4, data distribution has heavy tail. Theorem 5 says nothing about the case of heavy tail.\n    Heuristically, since the proposed method relies on order statistics, it may be expected that it has excellent performance when data distribution has heavy tail. In this view, the results in Case 3 and 4 are reasonable.\n\n**Limitations:** \n\n  * *One limitation...*\n\n    **Reply:** It is convenient to use the reduced $2p$ observations to conduct statistical inference since these observations also satisfy the linear model; see the formula between line 179 and 180 of the main text. From taht formula, it can be seen that the reduced observations satisfy a regression whose error term $\\bar {\\varepsilon}\\_j$ is a mean of independent stochastic errors. Hence from central limit theorem, it can be expected that $\\bar {\\varepsilon}\\_j$ behaves just like a normal random error. In this view, classical statistical method for Gaussian data can be used which may produce asymptotically correct inference result.\n\n  * *The supplementary...*\n\n    **Reply:** Thank you for your suggestion. We promise that if the paper is accepted, we will open-source our code.\n\n"}, {"review_id": "JTPCaIfmAL", "reviewer": "Reviewer_A4mU", "summary": "This paper studies the linear regression problem and proposes a new sketching method based on data averaging.", "questions": "This topic falls outside my current expertise. My review below is pretty limited.\n\n- Comparison against some standard sketching methods such as Gaussian sketch could be valuable. Gaussian sketch is computationally expensive to apply but easier to analyze. I'm curious how the convergence rates would compare.\n\n- I wonder if some plots would help improve the presentation. Perhaps something similar to Figure 7 of Pilanci and Wainwright [2016]?\n\nMinor:\n- There is a latex-related issue with the reference numbering in the appendix.\n- typo in line 166: \"no mater\"", "limitations": "Please see the \"questions\" section.", "rating": 6, "confidence": 1, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Please see the \"questions\" section.", "weaknesses": "Please see the \"questions\" section.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "\n**Comment:**  We sincerely appreciate your valuable comments and suggestions, particularly the positive feedback regarding the proposed method. \n\nBelow, we will provide a response that centers around  **questions** related to our paper. These relevant sections have been indicated with italicized text in the following.  We believe and hope that our revisions will ultimately meet your satisfaction, and that the paper will have the opportunity to be published in NeurIPS.\n\n\n**Questions:** \n\n\n\n  * *Comparison against some standard sketching methods such as Gaussian sketch could be valuable. Gaussian sketch is computationally expensive to apply but easier to analyze. I'm curious how the convergence rates would compare.*\n  \n    **Reply:** Thanks a lot for your suggestion. \n    Specifically, sketching methods solve the sketched least square problem\n    \\begin{eqnarray}\n    \\min\\_{\\boldsymbol{\\beta} \\in \\mathbb{R}^{p+1}}||\\mathbf{O}^{\\top} \\mathbf{y}-\\mathbf{O}^{\\top} \\mathbf{X} \\boldsymbol{\\beta}||^2 \n    \\end{eqnarray}\n    where $\\mathbf{O} \\in \\mathbb{R}^{N \\times n}$ is a sketching matrix with i.i.d.\\ standard Gaussian entries. \n    Suppose $p \\ll n \\ll N$.\n    From Theorem 1 of Ahfock et. al. (2021),\n    under the conditions of our Theorem 5,\n    the conditional mean squared error for Gaussian sketching least square estimator has convergence rate $\\sigma_{\\varepsilon}^2 p / n$, which is the same as the uniform sampling method.\n    To make the computation of the least square estimator based on sketched data computed within $O(n p + p^3)$ time, we let $n \\asymp N / p$.\n    In this case, the convergence rate of Gaussian sketching is $O(\\sigma_{\\varepsilon}^2 p^2 / N)$.\n    In this case,  the convergence rate of the proposed method, which is\n    $(1 + o_P (1)) \n    \\frac{p^2\\sigma\\_{\\varepsilon}^2}{\n    2 log(2p)N\n    }\n    $,\n    is faster than Gaussian sketchign by a factor of order $\\log(p)$.\n\n  * *I wonder if some plots would help improve the presentation. Perhaps something similar to Figure 7 of Pilanci and Wainwright [2016]?*\n  \n    **Reply:** Thank you for your suggestion. \n    We will consider improving the presentation of our paper.\n\n  \n  * *There is a latex-related issue with the reference numbering in the appendix.    \n  Typo in line 166: \"no mater\"*\n  \n    **Reply**:  Thank you for pointing out the typos. \n    We will correct typos we found, \n    and fix the reference numbering problem.\n\n  \n\n**Reference**\n\nAhfock D , Astle W J , Richardson S .Statistical properties of sketching algorithms[J].  2021."}]

Author rebuttal:
Thank all reviewers warmly for the time you took to review and understand our paper.

Most reviewers pointed out that the presentation of the paper should be improved. The discussion is dense, mathematical. Due to space constraints, the experimental results are all deferred to the appendix. Following the reviewers suggestion for the presentation of our paper, we would like move some details of the derivations for the theorems to the supplemental material and introduce some of the empirical results from the appendix to clarify the utility of the proposed averaging method.

Specifically, the primary objective of presenting Theorems 3 and 4 in the main paper is to provide a comprehensive understanding of the theoretical framework. It is important to note that Theorem 5 assumes a standard normal distribution, which can be seen as a specific case of the more general results provided by Theorems 3 and 4. In the camera-ready version, we would like to move Theorems 3 and 4 to the appendix. 
Then we can include a concise summary of key experimental results. If space permits, we will move Table 1 and Table 5 to the main text.

To make the presentation of the results clear, we plan to add Table 1 in the PDF file which summarizes the theoretical performance of the proposed method and compare it with the ideal sampling method implied by Theorem 2 and the IBOSS algorithm.
In short, the merit of Theorem 6 is that it shows that IBOSS can match the lower bound provided by Theorem 2; 
the merit of Theorem 5 is that it shows that there exists a method (i.e., the proposed method) which can break the lower bound of Theorem 2.

Reviewer EUfX, Reviewer DaqC and Reviewer sefb showed that the Gaussian assumption is a bit restrictive.
Firstly, we clarify the assumptions in Theorems. 
Theorem 1 may be generalized to non-Gaussian case. We do not think Theorem 2 can be easily generalized to non-Gaussian distributions since in that case, one may need to deal with the concentration inequality of more general order statistics.
Theorems 3 and 4 are more general than Theorem 5, and are used for the proof of Theorem 5. These two theorems do not assume Gaussian distribution, the cost is some involved characterizations of the tail behavior. While these two theorems may shed some lights on the behavior for non-Gaussian distribution, we actually did not rigorously give results for concrete non-Gaussian distributions.
They imposes some conditions on the tail probability of $\textbf{Z}$.
But the imposed conditions may not be easily verified for sub-Gaussian data.
Theorems 5 and 6 rely on Gaussian assumptions. 
There are some $\log$ in the statement of Theorem 2 and 3.
These $\log$ comes from Gaussian data.
For sub-Gaussian data, the $\log$ term may not be correct.
Previously, one of our goal is to try to make the results general and go beyond normal distribution.
It turns out that the cost is the greatly reduced readability, which may not be a good trade-off.


Meta-review: Dear Authors,

Thank you for your valuable contribution to Neurips and the ML community. Your submitted paper has undergone a rigorous review process, and I have carefully read feedback provided by the reviewers and considered the author rebuttal in detail.

This paper is focused on the asymptotic estimation risk of large-scale linear regression with various data sketching methods. The authors refine an existing estimation lower-bound and develop a new method based on data averaging. The method is shown to outperform uniform subsampling and the IBOSS method in theory and in practice. The reviewers all agree that the contribution is solid and interesting.

Given this positive assessment, I am willing to recommend the acceptance of your paper for publication.

I would like to remind you to carefully review the reviewer feedback and the resulting discussion. While most reviews were positive, the reviewers have offered valuable suggestions that can further strengthen the quality of the paper. Please take another careful look a the 'weaknesses' section of each reviewer comment. Please also review the discussions with Reviewer sefb, whose review was the most critical. I encourage you to use this feedback to make any necessary improvements and refinements before submitting the final version of your paper.

Once again, thank you for submitting your work to Neurips.

Best,
Area Chair
Decision: Accept (poster)