Year: 2021

Abstract:
In large-scale supervised learning, after a model is trained with an initial dataset, a common challenge is how to exploit new incremental data without re-training the model from scratch. Motivated by this problem, we revisit the canonical problem of dynamic least-squares regression (LSR), where the goal is to learn a linear model over incremental training data. In this setup, data and labels (A(t),b(t)) 2 Rt⇥d⇥Rt evolve in an online fashion (t d), and the goal is to efficiently maintain an (approximate) solution of minx(t) kA(t)x(t) bk2 for all t 2 [T ]. Our main result is a dynamic data structure which maintains an arbitrarily small constant approximate solution to dynamic LSR with amortized update time O(d1+o(1)), almost matching the running time of the static (sketching-based) solution. By contrast, for exact (or 1/ poly(n)-accuracy) solutions, we show a separation between the models, namely, that dynamic LSR requires ⌦(d2 o(1)) amortized update time under the OMv Conjecture (Henzinger et al., STOC’15). Our data structure is fast, conceptually simple, easy to implement, and our experiments demonstrate their practicality on both synthetic and real-world datasets.

heading: 1 INTRODUCTION
text: The problem of least-squares regression (LSR) dates back to Gauss in 1821 (Stigler, 1981), and is the backbone of statistical inference (Hastie et al., 2001), signal processing (Rabiner & Gold, 1975), convex optimization (Bubeck, 2015), control theory (Chui, 1990) and network routing (Lee & Sidford, 2014; Madry, 2013). Given an overdetermined (n d) linear system A 2 Rn⇥d,b 2 Rn, the goal is to find the solution vector x that minimizes the mean squared error (MSE)
min x2Rn kAx bk2. (1)
Among many other loss functions (e.g., `p) that have been studied for linear regression, `2-regression has been the most popular choice as it is at the same time robust to outliers, and admits a highaccuracy efficient solution.
The computational task of least-squares regression arises naturally in high-dimensional statistics and has been the central of focus. The exact closed-form solution is given by the well-known Normal equation x? = (A>A) 1A>b, which requires O(nd2) time to compute using naive matrixmultiplication, or O(nd! 1) ⇡ O(nd1.37) time using fast matrix-multiplication (FMM) (Strassen, 1969) for the current FMM exponent of ! ⇡ 2.37 (Le Gall, 2014; Alman & Williams, 2021). Despite the elegance and simplicity of this closed-form solution, in practice the latter runtime is often too slow, especially in modern data analysis applications where both the dimension of the feature space (d) and the size of datasets (n) are overwhelmingly large. A more modest objective in attempt to circumvent this computational overhead, is to seek an ✏-accurate solution that satisfies
kAx bk2  (1 + ✏) min x2Rd kAx bk2.
This was the primary motivation behined the development of the sketch-and-solve paradigm, where the idea is to first compress the matrix into one with fewer (⇠ d/✏2) rows and then to compute the standard LSR solution but over the smaller matrix. A long line of developments on this framework culminates in algorithms that run in close to input-sparsity time (Sarlos, 2006; Clarkson & Woodruff, 2017; Nelson & Nguyên, 2013; Chepurko et al., 2021). In particular, a direct application of sketchand-solve yields an algorithm runs in eO(nnz(A)✏ 1 + d!) 1, which is near optimal in the “low
1In this paper we use eO(·) to hide polylogarithmic terms, and we use O✏(·) to hide poly(log d, ✏ 1) terms.
precision” regime (✏ = 1/ poly(log d)). Interestingly, when combined with more sophisticated ideas of preconditioning and (conjugate) gradient descent, the runtime of this algorithm in terms of the error ✏ can be further improved to eO(nnz(A) log(1/✏) + d!), which yields a high precision algorithm, i.e., it can efficiently solve the problem to within polynomial accuracy ✏ = 1/ poly(d).
Dynamic least-squares In many real-world scenarios of the aforementioned applications, data is evolving in an online fashion either by nature or by design, and such applications require maintaining the solution (1) adaptively, where rows of the data matrix and their corresponding labels (A(t),b(t)) arrive one-by-one incrementally. This is known as the dynamic least-squares regression problem.
The origins of dynamic least-squares regression was in control theory of the 1950’s (Plackett, 1950), in the context of dynamical linear systems. In this setup, the data matrix [A(t),b(t)] corresponds to the set of measurement and it evolves in an online (incremental) fashion, and the goal is to efficiently maintain the (exact) solution to a noisy linear system b := A(t)x(t) + ⇠(t) without recomputing the LSR solution from scratch. The recursive least-squares (RLS) framework and the celebrated Kalman filter (Kalman, 1960) provide a rather simple update rule for maintaining an exact solution for this problem, by maintaining the sample covariance matrix and using Woodburry’s identity (which assert that an incremental update to [A(t),b(t)] translates into a rank-1 update to the sample covariance matrix), and hence each update can be computed in O(d2) time (Kalman, 1960).
Beyond this classic motivation for dynamic LSR, a more timely motivation comes from modern deep learning applications: Most neural networks need to be frequently re-trained upon arrival on new training data, in order to improve prediction accuracy, and it is desirable to avoid recomputing weights from scratch. This problem of efficient incremental training of DNNs has been studied before in elastic machine learning (Liberty et al., 2020) and in the context of continual learning (Parisi et al., 2019). Our work sheds light on this question by analyzing the minimal computational resources required for `2 loss-minimization.
Despite the rich and versatile literature on static LSR, the understanding of the dynamic counterpart was so far quite limited: The previous best known result requires O(d2) amortized update time (by a direct application of the Woodbury identity). The basic questions we address in this papers are:
Is it possible to achieve faster update time for maintaining an exact solution? How about a smallapproximate solution – Is it then possible to achieve amortized O(d) or even input-sparsity time?
In this paper, we settle both of these questions and present an essentially complete characterization of the dynamic complexity of LSR.

heading: 1.1 OVERVIEW OF OUR RESULTS
text: Our first result is a negative answer to the first question above of maintaining exact (or polynomialaccuracy) LSR solutions in the dynamic setting – We prove that Kalman’s approach is essentially optimal, assuming the popular Online Matrix-Vector (OMv) Conjecture (Henzinger et al., 2015) 2: Theorem 1.1 (Hardness of exact dynamic LSR, informal). Assuming the OMv Conjecture, any dynamic algorithm that maintains an ✏ = 1/ poly(d)-approximate solution for the dynamic LSR problem over T = poly(d) iterations, must have ⌦(d2 o(1)) amortized update time per iteration.
Theorem 1.1 separates the static and the dynamic complexities of the exact LSR problem: As mentioned above, the static problem can be solved by batching rows together using FMM in time O(Td! 1), whereas the dynamic problem requires ⌦(Td2) by Theorem 1.1. Indeed, the implication Theorem 1.1 is stronger, it also separates the static and dynamic complexity of approximate LSR problem under the high precision regime, it asserts that a polylogarithmic dependence on the precision (i.e. d poly(log(1/✏))) on update time is impossible (assuming OMv), in sharp contrast to the static case.
We next focus on an approximate version of this classic online problem, dynamic ✏-LSR, where the goal is to efficiently maintain, during all iterations t 2 [T ], an ✏-approximate solution under incremental row-updates to A(t) and labels b(t), where efficiency is measured by the amortized
2This conjecture postulates that multiplying a fixed d ⇥ d matrix A with an online matrix B, column-bycolumn (ABi), requires d3 o(1) time, in sharp contrast to the batch setting where this can be done using FMM in d! ⌧ d3 time. See Section 4.
update time for inserting a new row. A natural complexity benchmark for this dynamic problem is the aforementioned best static sketch-and-solve solution, which for n = T is eO(nnz(A(T ))✏ 1 + d!) = eO(nnz(A)✏ 1) for T d. Our main result is a provably efficient and practical dynamic data structure, whose total running time essentially matches the complexity of the offline problem: Theorem 1.2 (Main result, informal version of Theorem 3.1). For any accuracy parameter ✏ > 0, there is a randomized dynamic data structure which, with probability at least 0.9, maintains an ✏- approximate solution to the dynamic LSR problem simultaneously for all iterations t 2 [T ], with total update time
O(✏ 2 nnz(A(T )) log(T ) + ✏ 6d3 log5(T )).
Theorem 1.2 almost matches the fastest static sketching-based solution, up to polylogarithmic terms and the additive FMM term. When T d, this theorem shows that amortized update time of our algorithm is O(d1+o(1)).

heading: 1.2 RELATED WORK
text: Sketching and sampling The least squares regression as a fundamental problem has been extensively studied in the literature. A long line of work (Ailon & Chazelle, 2006; Clarkson & Woodruff, 2017; Nelson & Nguyên, 2013; Avron et al., 2017; Cohen et al., 2015; Woodruff, 2014; 2021) have focused on using dimension reduction technique (sketching or sampling) to speedup the computation task, culminates into algorithms that run in eO(nnz(A) log(1/✏) + d!) time (Clarkson & Woodruff, 2017). See Appendix A for detailed discussions.
Regression in online, streaming, and sliding window models Least-squares regressions have also been studied in various computational models, though the focus of these models are generally not the (amortized) running time. Our algorithm uses techniques developed by Cohen et al. (2020), where they study the regression problem in the online model, with the goal of maintaining a spectral approximation of data matrix in the online stream. Their algorithm only needs to store O✏(d) rows but the amortized running time is still ⌦(d2) (see Section 3.1 for detailed discussions). In the streaming model, the main focus is the space complexity, and a direct application of random Gaussian sketch or count sketch reduces the space complexity to eO(d2✏ 1) and it is shown to be tight (Clarkson & Woodruff, 2009). Recent work of (Braverman et al., 2020) studies regressions and other numerical linear algebra tasks in the sliding window model, where data come in an online stream and only the most recent updates form the underlying data set. The major focus of a sliding window model is still the space complexity, and there is no amortized running time guarantee.
Disparity from online learning Our work crucially differs from online learning literature (Hazan, 2019), in that the main bottleneck in online regret-minimization and bandit problems is informationtheoretic, whereas the challenge in our loss-minimization problem is purely computational. See Appendix A for detailed discussions.

heading: 2 PROBLEM FORMULATION
text: In a dynamic least-squares regression problem, initially, we are given a matrix A(0) 2 Rn0⇥d together with a vector b(0) 2 Rn0 . At the t-th step, a new data of form ((a(t))>, (t)) 2 Rd ⇥ R arrives, and the goal is to maintain an ✏-approximate solution. A formal description is provided below, where we assume n0 = d+ 1 for simplicity (see Remark 2.3). Definition 2.1 (Dynamic least-squares regression). Let d 2 N+ and ✏ 2 [0, 1) be two fixed parameters. We say an algorithm solves ✏-approximate dynamic least squares regression if
• The data structure is given a matrix A(0) 2 R(d+1)⇥d and a vector b(0) 2 Rd+1 in the preprocessing phase.
• For each iteration t 2 [T ], the algorithm receives updates a(t) 2 Rd and (t) 2 R. Define A
(t) := [(A(t 1))>,a(t)]> 2 R(d+t+1)⇥d to be A(t 1) appended with a new row (a(t))>, and b(t) := [(b(t 1))>, (t)]> 2 Rd+t+1 to be b(t 1) appended with a new entry (t). After this update, the algorithm outputs an ✏-approximate solution x(t) 2 Rd:
kA(t)x(t) b(t)k2  (1 + ✏) min x2Rd kA(t)x b(t)k2.
We write [0 : T ] = {0, 1, . . . , T}, and for any t 2 [0 : T ], we denote M(t) := [A(t),b(t)] 2 R(d+t+1)⇥(d+1). We make the following assumptions. Assumption 2.2. We assume 1. Each data have bounded `2 norm, i.e., 8i 2 [T +d+1], the i-th row of M(T ) satisfies kM(T )
i,⇤ k2  D. 2. The initial matrix M(0) has full rank, and its smallest singular value is bounded by d+1(M(0)) for some polynomially small 2 (0, 1). Remark 2.3. We remark that these assumptions are essentially w.l.o.g. for the following reasons: 1. Real world data inherently have bounded `2 norm, and in applications like machine learning, data are often normalized. 2. We can assume the initial matrix A(0) has d+1 rows because brute-forcely adding these d + 1 initial rows would only take O(d3) time, and this is within our desired total running time of O✏(nnz(A(T )) + d3). 3. To satisfy the assumption that d+1(M(0)) for some polynomially small , we could let the initial matrix M(0) = · Id+1. This is equivalent to adding a small regularization term of · kxk2 and this incurs only a polynomially small additive error.

heading: 3 DYNAMIC ✏-LSR DATA STRUCTURE
text: In this section, we provide an approximation algorithm for the dynamic least squares regression. Notably, our algorithm maintains an ✏-approximate solution in near input sparsity time. Theorem 3.1 (Data structure for dynamic least squares regression). Let ✏ > 0, d, T 2 N. There exists a randomized algorithm for dynamic least-squares regression (Algorithm 1–4). With probability at least 0.9, the algorithm maintains an ✏-approximation solution for all iterations t 2 [T ] and the total update time over T iterations is at most O ✏ 2 nnz(A(T )) log(T )+ ✏ 6d3 log5(TD/ ) . Our
data structure uses at most O ✏ 2d2 · log2(TD/ ) space.
Notations We use a superscript (t) to denote the matrix/vector/scaler maintained by the data structure at the end of the t-th iterations. In particular, the superscript (0) represents the variables after the preprocessing step. For any matrix A 2 Rn⇥d, i 2 [n] we define its leverage score ⌧(A) 2 Rn as ⌧i(A) := a>i (A>A)†ai. We define the generalized leverage score (same as Cohen et al. (2015)) of A 2 Rn⇥d with respect to another matrix B 2 Rn0⇥d as : ⌧B
i (A) := a> i (B>B)†ai.
For more properties of the leverage scores see Section B.1.

heading: 3.1 TECHNIQUE OVERVIEW
text: Our approach is formally described in Algorithm 1–4, we first overview the ideas behind it.
From a high level view, our approach follows the online row sampling framework (Cohen et al., 2015; 2020; Braverman et al., 2020): When a new row arrives, we sample and keep the new row with probability proportional to the (approximated version of) online leverage score
⌧M (t 1) d+t+1 (M (t)) = (m(t))>((M(t 1))>M(t 1)) 1m(t).
The sampled matrix is a spectral approximation to the true data matrix. We maintain an approximate least-squares regression solution using this sampled matrix.
Naively computing the online leverage score takes O(d2) time. In order to accelerate this computation, we use two approximations:
1. Similar to Cohen et al. (2020), we compute the online leverage scores with respect to the sampled matrix instead of the true data matrix. However, this idea alone is still not enough to achieve sub-quadratic time.
2. We use a JL-embedding3 trick (Spielman & Srivastava, 2011) to compress the size of the d⇥d matrix to⇡ ✏ 2⇥d. In this way, in each iteration it only takes O✏(d) time to compute the approximate online leverage score.
We further use an inductive analysis to bound the overall error (Lemma 3.4).
Finally, we adopt a similar strategy as Cohen et al. (2020) to prove that sampling according to the approximate online leverage score still keeps at most O✏(d) sampled rows (Lemma 3.7). Whenever a row is sampled, it takes O(d2) time to update the maintained matrices using Woodburry identity. Hence, the amortized update time of the sampled rows is O✏(d3/T ) = o(d) when T d.
3Johnson-Lindenstrauss (JL) Lemma shows a way to embed high-dimensional vectors to low-dimensional space while preserving the distances between the vectors. A rigorous statement is shown in Lemma B.4.
Remark 3.2 (Difference from sketching-based solutions). Our approach crucially differs from the sketching-based solutions, which do not provide any speedup over the direct application of Woodburry identity (O(d2) time per iteration). A sketching-based solution maintains a sketched matrix SM 2 RO✏(d)⇥d, where S is a sketching matrix (e.g. SRHT (Ailon & Chazelle, 2006) or Count Sketch (Clarkson & Woodruff, 2017)) that mixes the rows of M. When a new row of M arrives, at least one row of the sketched matrix SM needs to be updated, in contrast to our sampling-based approach where the sampled matrix is not updated in most of the iterations.
Implementation We explain the detailed implementation of our algorithm. At the beginning of the t-th iteration, a sampling matrix D(t 1) is derived based on the online leverage score, and the subsampled matrix N(t 1) = D(t 1)M(t 1) maintains a spectral approximation on the column space of M(t 1) = [A(t 1),b(t 1)]. Let s(t 1) denote the number of sampled rows. To obtain spectral approximation, we maintain the approximate covariance matrices H(t 1) = ((N(t 1))>N(t 1)) 1 and B(t 1) = N(t 1)H(t 1). The online leverage score ⌧ (t) of a new row m(t) = ((a(t))>, (t)) can be approximated as kB(t 1)m(t 1)k2. To efficiently compute the leverage score of a new row, we left-multiply by a JL matrix J(t 1) and maintain a proxy eB(t 1) = J(t 1)B(t 1), with the guarantee that keB(t 1)mk2 ⇡ kB(t 1)mk2 for any m 2 Rd+1 with high probability.
When a new row m(t) is inserted at the t-th iteration, we sample it via the approximate online leverage score (Line 3 of SAMPLE). We only perform update if the new row is sampled. In that case, we renew the JL matrix and perform a series of careful updates on all the variables that we maintain (See UPDATEMEMBERS). To obtain the final solution x(t) 2 Rd, we solve minx2Rd kD(t)A(t)x D(t)b(t)k2, which has the closed-form solution of x(t) = G(t) · u(t) and can be efficiently maintained by taking G(t) = ((A(t))>(D(t))2A(t)) 1 2 Rd⇥d and u (t) = (A(t))>(D(t))2b(t).
Algorithm 1 PREPROCESS (A, b, ✏, T )
1: M [A,b] . M 2 R(d+1)⇥(d+1) 2: D Id+1
# Spectral approximation 3: s d+ 1 4: N D ·M . N 2 Rs⇥(d+1) 5: H ((N)>N) 1 . H 2 R(d+1)⇥(d+1) 6: B N ·H . B 2 Rs⇥(d+1)
# JL approximation 7: O(1/T 2), k O(✏ 2 log(T/ )) 8: J JL(s, ✏, , T ) . JL matrix J 2 Rk⇥s 9: eB J ·B . eB 2 Rk⇥(d+1)
# Maintain solution 10: G (A>D>DA) 1 . G 2 Rd⇥d 11: u A>D2b . u 2 Rd 12: x G · u . x 2 Rd
Algorithm 2 UPDATE (a, )
1: m [a>, ]> . m 2 Rd+1 2: ⌫ SAMPLE(m) . ⌫ 2 R 3: D  D 0 0 ⌫
4: if ⌫ 6= 0 then UPDATEMEMBERS(m) 5: return x
Algorithm 3 SAMPLE (m)
1: ⌧ keB ·mk22 2: p min{3(1 + ✏)2✏ 2⌧ log(1/ ), 1} 3: ⌫ 1/pp with probability p, and ⌫ 0
otherwise
Algorithm 4 UPDATEMEMBERS (m) # Update spectral approximation
1: s s+ 1 2: H Hmm > H/p
1+m>Hm/p
3: H H+ H 4: B [(B+N · H)>, H ·m/pp]> 5: N [N>,m/pp]>
# Update JL approximation 6: J JL(s, ✏, , T ) 7: eB J ·B
# Update solution 8: G G Gaa > G/p
1+a>Ga/p
9: u u+ · a/p 10: x G · u
We outline the proof of Theorem 3.1, and defer the detailed proof to Appendix C due to space limits.

heading: 3.2 CORRECTNESS
text: We show the correctness of our algorithm and prove it maintains an ✏-approximate solution for all iterations with high probability. We start with closed-form formulas for all the variables we maintain.
Lemma 3.3 (Closed-form formulas). At the t-th iteration of Algorithm 1 – 4, we have
1. M(t) = [A(t),b(t)] 2 R(d+t+1)⇥(d+1).
2. D(t) 2 R(d+t+1)⇥(d+t+1) is a diagonal matrix with s(t) non-zero entries.
3. N(t) = (D(t)M(t))S(t),⇤ 2 Rs (t)⇥(d+1), where S(t) ⇢ [d + t + 1] is defined as the set of
non-zero entries of D(t).
4. H(t) = (N(t))>N(t) 1 2 R(d+1)⇥(d+1).
5. B(t) = N(t)H(t) 2 Rs(t)⇥(d+1).
6. eB(t) = J(t) ·B(t) 2 Rk⇥(d+1), where k = O(✏ 2 log(T/ )).
7. G(t) = (A(t))>(D(t))2A(t) 1 2 Rd⇥d.
8. u(t) = (A(t))>(D(t))2b(t) 2 Rd.
9. x(t) = (A(t))>(D(t))2A(t) 1 · (A(t))>(D(t))2b(t) 2 Rd.
The following lemma is key for our correctness analysis. It shows that we maintain a good approximation on online leverage scores and a spectral approximation of M(t) throughout all iterations. Lemma 3.4 (Spectral approximation via leverage score maintenance). With probability at least 1 2T ,
(1 ✏)2⌧M (t 1) d+t+1 (M (t))  ⌧ (t)  (1 + ✏)2⌧M (t 1) d+t+1 (M (t)), 8t 2 [T ], (2)
and (M(t))>(D(t))2M(t) ⇡✏ (M(t))>M(t), 8t 2 [0 : T ]. (3)
Proof Sketch. We prove by induction and show that with probability 1 2t , Eq. (3) holds for all t0 2 [0 : t] and Eq. (2) holds for all t0 2 [t]. The base case t = 0 holds trivially, as D(0) = Id+1, and therefore, (M(0))>(D(0))2M(0) = (M(0))>M(t). Given the induction hypothesis upon t 1, we proceed in the following three steps.
• We first use the induction hypothesis to prove that kB(t 1)m(t)k is a good estimate on the online leverage score, that is
(1 ✏)⌧M (t 1) d+t+1 (M (t))  kB(t 1) ·m(t)k22  (1 + ✏)⌧M (t 1) d+t+1 (M (t)).
• We then use the JL lemma (Lemma B.4) to show that with probability 1 , the sketched covariance matrix keB(t 1)m(t)k returns good estimation kB(t 1)m(t)k. That is
(1 ✏)2 · ⌧M (t 1) d+t+1 (M (t))  keB(t 1) ·m(t)k22  (1 + ✏)2 · ⌧M (t 1) d+t+1 (M (t)). (4)
• Finally, we wrap up the proof by proving the second part of induction. In particular, we show that conditioned on Eq. (4) holds, we have
(M(t))>(D(t))2M(t) ⇡✏ (M(t))>M(t).
The proof then follows by an union bound over failure events.
It is well known that spectral approximations of (M(t))>M(t) give approximate solutions to least squares regressions (Woodruff, 2014), so we have proved the correctness of our algorithm. Lemma 3.5 (Correctness of Algorithm 1–4). With probability at least 1 O(1/T ), in each iteration, UPDATE of Algorithm 2 outputs a vector x(t) 2 Rd such that
kA(t)x(t) b(t)k2  (1 + ✏) min x2Rd kA(t)x b(t)k2.

heading: 3.3 TIME ANALYSIS
text: Next, we bound the overall update time of our algorithm. We first compute the worst case update time of Algorithm 2. When ⌫(t) = 0, i.e., the t-th row is not sampled, the UPDATE procedure only needs to compute the approximate leverage score ⌧ (t) (SAMPLE, Algorithm 3), and it takes O(k · nnz(m(t))) time. When ⌫(t) 6= 0, i.e., the t-th row is sampled, the UPDATE procedure makes a call to UPDATEMEMBERS (Algorithm 4), and it takes O(k · s(t) · d) time. Plugging in the value of k, we have the following lemma. Lemma 3.6 (Worst case update time). At the t-th iteration of the UPDATE procedure (Algorithm 2),
• If ⌫(t) = 0, then UPDATE takes O ✏ 2 log(T/ ) · nnz(a(t)) time.
• If ⌫(t) 6= 0, then UPDATE takes O ✏ 2s(t)d log(T/ ) time.
To bound the amortized update time, we need to bound the total number of sampled rows, and this is closely related to the sum of online leverage scores. Such an upper bound was already established by Cohen et al. (2020), here we present a slightly generalized version of it. Lemma 3.7 (Sum of online leverage scores, generalization of Theorem 2.2 of (Cohen et al., 2020)). If the matrix M(T ) satisfy Assumption 2.2, then
TX
t=1
⌧M (t 1) d+t+1 (M (t))  O(d log(TD/ )).
Now we are ready to bound the amortized update time of our algorithm. Lemma 3.8 (Amortized update time). With probability at least 0.99, the total running time of UPDATE over T iterations is at most O ✏ 2 nnz(A(T )) log(T ) + ✏ 6d3 log5(TD/ ) .
Proof Sketch. In this proof sketch we simplify the second term as d3 · poly(✏ 1 log(TD/ )). The first term comes from the computation cost of querying leverage score, which takes O ✏ 2 log(T/ )·
nnz(a(t)) time in the t-th iteration even if the t-th row is not sampled. The second term bounds the total update time for the sampled rows:
• From Lemma 3.4 and Lemma 3.7, with high probability the sum of the approximate online leverage scores ⌧ (t) are bounded by O(d log(TD/ )).
• Using Markov inequality, the total number of sampled rows is bounded by
s(T ) = O( TX
i=1
p(t)) = O(✏ 2 · log(1/ ) · TX
t=1
⌧ (t))  O d · poly(✏ 1 log(TD/ )) .
• Since there are s(T ) sampled rows, and for each sampled row we update data structure members in  O s(T )d · poly(✏ 1 log(TD/ )) time, the total update time for sampled
rows is
s(T ) ·O s(T )d · poly(✏ 1 log(TD/ )) = O d3 · poly(✏ 1 log(TD/ )) .

heading: 4 HARDNESS RESULT
text: We prove a ⌦(d2 o(1)) amortized time lower bound for dynamic least squares regression with high precision, assuming the OMv conjecture. The OMv conjecture was originally proposed by Henzinger et al. (2015), and it is widely accepted in the theoretical computer science community. Conjecture 4.1 (OMv conjecture, (Henzinger et al., 2015)). Let d 2 N, T = poly(d). Let > 0 be any constant. B 2 {0, 1}d⇥d is a Boolean matrix. 8t 2 [T ], a Boolean vector z(t) 2 {0, 1}d is revealed at the t-th step. We say an algorithm solves the OMv problem if it returns the Boolean matrix-vector product Bz(t) 2 Rd at every time step. The conjectures states that there is no algorithm that solves the OMv problem using poly(d) preprocessing time and O(d2 ) amortized running time, and has an error probability  1/3.
The results in this section are all under the Word RAM model where the word size w = O(log d). Our main result is formally stated below. Theorem 4.2 (Hardness of dynamic-least squares regression with high precision). Let d 2 N, T = poly(d), ✏ = 1
d8T 2 = 1/ poly(d), and let > 0 be any constant. Assuming the OMv conjecture is
true, any dynamic algorithm that maintains an ✏-approximate solution of the least squares regression requires at least ⌦(d2 ) amortized time per update. Our lower bound is proved by first reducing the standard OMv conjecture for Boolean matrices to OMv-hardness for well-conditioned positive semidefinite (PSD) matrices over real numbers. Then we use this new OMv-hardness result to prove our lower bound for dynamic least squares regression. We only provide a proof sketch here and detailed proof are delayed to Section D.
OMv-hardness for well-conditioned PSD matrix The OMv conjecture asserts the hardness of solving online Boolean matrix-vector product exactly. We extend it to solving online real-valued matrix-vector product for well-conditioned PSD matrices, while allowing polynomially small error. Lemma 4.3 (Hardness of approximate real-valued OMv). Let d 2 N, T = poly(d). Let > 0 be any constant. Let H 2 Rd⇥d be a symmetric matrix whose eigenvalues satisfy 1  d(H)  · · ·  1(H)  3. For any t 2 [T ], z(t) 2 Rd is revealed at the t-th step, and kz(t)k2  1. Assuming the OMv conjecture is true, then there is no algorithm with poly(d) preprocessing time and O(d2 ) amortized running time that can return an O(1/d2)-approximate answer to Hz(t) for all t, i.e., a vector y(t) 2 Rd s.t. ky(t) Hz(t)k2  ✏, and has an error probability  1/3.
Proof Sketch. Given a Boolean matrix B 2 {0, 1}d⇥d in the OMv conjecture, we construct a PSD matrix H =  2Id 1 d B
1 d B > 2Id
2 R2d⇥2d. We note that H is symmetric and 1  d(H)  1(H)  3.
Given a binary OMv query vector z(t), we construct z(t) = (0d, z(t)) 2 R2d. Since H has a constant condition number, we can prove that rounding an ✏ ⇠ 1/d2-approximate answer by ⇡✏ H · z(t) still gives the correct binary answer to Bz(t).
Reducing OMv to dynamic least-squares regression We next wrap up the proof of Theorem 4.2 by reducing OMv to dynamic ✏-LSR.
Proof Sketch of Theorem 4.2. Given a PSD matrix H and a sequence of query {Hz(t)}T t=1 of the problem in Lemma 4.3, we reduce it to a dynamic ✏-LSR, where the initial A is such that A>A = H
1 (this preprocessing step of the reduction takes ⇠ d! time), and the label is 0 for the initial d data. For each t 2 [T ], the incoming row a(t) is a small scaled version of z(t), i.e, a(t) =
1 d2 p T · z(t) 2 Rd, and the label is 1. Let x(t)? be the optimal solution at the t-th step, H(t) := ((A(t))>A(t)) 2, the reduction is complete via the following three steps:
• Step 1. x(t) and x(t)? are close, i.e., x(t) = x (t) ? ±O( 1
d4 p T ).
• Step 2. x(t) x(t 1) recovers H(t 1)a(t), i.e., x(t) x(t 1) = H(t 1)a(t) ±O( 1 d4 p T ).
• Step 3. H(t 1)a(t) is close to Ha(t), i.e., H(t 1)a(t) = Ha(t) ±O( 1 d6 p T ).
In particular, let y(t) = d2 p T (x(t) x(t 1)), Step 2 and 3 directly implies ky Hz(t)k2  O(1/d2). This completes the proof.

heading: 5 EXPERIMENTS
text: Our method is most suitable for data distributions that are non-uniform. Indeed, if the data has low coherence (they are all similar to each other), then the naive uniform sampling is as good as leverage score sampling. We perform empirical evaluations on our algorithm over both synthetic and real-world datasets.
Synthetic dataset We follow the empirical study of (Dobriban & Liu, 2019) and generate data from the elliptical model. In this model a(t) = w(t)⌃z(t), where z(t) ⇠ N(0, Id) is a random Gaussian vector, ⌃ 2 Rd⇥d is a PSD matrix, and w(t) is a scaler. The label is generated as b(t) =
ha(t),x?i + w(t)⇠(t), where x? 2 Rd is a hidden vector and ⇠ ⇠ N(0, 1) is standard Gaussian noise. This model has a long history in multivariate statistics, see e.g. (Martin & Maes, 1979). In our experiments, we set ⌃ = Id for simplicity. In order to make the dataset non-trivial, we set w(t)
to be large (= p T ) for a few (= d/10) iterations, and small (= 1) for the rest of the iterations. We set T = 400000 and d = 500.
Real-world dataset We use the VirusShare dataset from the UCI Machine Learning Repository4. We select this dataset because it has a large number of features and data points, and has low errors when fitted by a linear model. The dataset is collected from Nov 2010 to Jul 2014 by VirusShare (an online platform for malware detection). It has T = 107888 data points and d = 482 features.
Baseline algorithms We compare with three baseline methods. 1. Kalman’s approach makes use of Woodburry identiy and gives an exact solution. 2. The uniform sampling approach samples new rows uniformly at random. 3. The row sampling approach samples new rows according to the exact online leverage scores. (Cohen et al., 2020)
Our experiments are executed on an Apple M1 CPU with codes written in MATLAB. We repeat all experiments for at least 5 times and take the mean. On both datasets, we initiate the model based on the first 10% of the data. The experiment results are formally presented in Figure 1 and more details can be found in Appendix E. Our algorithm consistently outperforms baseline methods: Our algorithm runs faster when achieving comparable error rates.

heading: 6 CONCLUSION
text: We provide the first practical and provably fast data structure for dynamic least-squares regression, obtaining nearly tight upper and lower bounds for this fundamental problem. On the algorithmic side, we design an ✏-approximation dynamic algorithm whose total update time almost matches the input sparsity of the (online) matrix. On the lower bound side, we prove that it is impossible to maintain an exact (or even high-accuracy) solution with⌧ d2 o(1) amortized update time under the OMv conjecture. As such, this result exhibits the first separation between the static and the dynamic LSR problems.
Our paper sets forth several interesting future directions. On the theoretical side, a very interesting question is whether it is possible to reduce the additive term d3 of our algorithm to matrixmultiplication time d!? A second open problem—of interest in both theory and practice—is whether it is possible to achieve input-sparsity amortized update time in the fully dynamic setting, i.e., when allowing both addition and deletion of data rows? Finally, it would be interesting to find connections between dynamic least-squares regression and incremental training of more complicated models, such as dynamic Kernel-ridge regression and to deep neural networks.
4https://archive.ics.uci.edu/ml/datasets.php

Review summaries:
This paper studies incremental least-squares regression, where the goal is to maintain an $(1+\epsilon)$-approximate solution to $\min_x \left\Vert Ax-b\right\Vert_2^2$ for some $A\in\mathbb{R}^{n\times d}$, under row insertions to $\begin{pmatrix}A & b\end{pmatrix}$, while keeping the total runtime as low as possible.   It builds on a work of Cohen et al. (2020), where it is shown that if, upon insertion of a row, we only keep (a multiple of) it with probability roughly equal to its current leverage score (otherwise we discard it), then the total number of rows that will be inserted is only $\widetilde{O}(d)$ and the matrix spectrally approximates the "true" matrix where all rows are kept. The main contribution of this paper is to make this algorithm faster, by improving the leverage score computation procedure. Instead of directly computing the leverage scores, the authors use the Johnson-Lindenstrauss lemma to be able to quickly approximate them. Then the main thing to be taken care of is maintaining certain dimension-reduced matrices that arise. In total the asymptotic amortized runtime goes down from $\widetilde{O}(d\cdot nnz(A))$ to $\widetilde{O}(nnz(A))$.   Additionally, the authors give a negative result, which shows that, under the OMv conjecture, the amortized runtime is $\Omega(d^2)$ if we require a high-precision solution ($\epsilon=1/\mathrm{poly}(d)$).  The theoretical result is accompanied by some synthetic and real experiments. In both, it seems that the runtime improvement is significant (e.g. up to 4x speedup in the real data when $\epsilon=1$).

The paper considers solving the regression problem min_x |Ax-b|_2 in the online setting, where A in R^{n x d} and b in R^n are given row by row, one at each time. The main task is to maintain a good approximate solution x (meaning that |Ax-b|_2 <= (1+eps)*OPT) throughout this process, with the update time as small as possible. The paper shows that if there are T updates to A (the initial A may not be empty), the total update time will be O(eps^{-2} nnz(A) log T + poly(eps^{-1} d log T)).   The first term O(eps^{-1}nnz(A)log T) matches the runtime in a typical sketching algorithm up to the log T factor, while the second term poly(eps^{-1} d log T) is much better in the dependence on T than the runtime of a simple sketching algorithm, which would be poly(d/eps)*T. This means that the amortized cost of each update is much smaller when T >> d.  The paper also proves a lower bound of Omega(d^{2-o(1)}) amortized update time for T = poly(d) updates, assuming the Online Matrix-vector conjecture (which is on the amortized runtime of multiplying a matrix with a vector in the online setting).

The paper discusses _Dynamic Least Squares_: the problem where the rows of a overdetermined least squares problem are revealed incrementally, and an algorithm has to maintain an accurate solution to the least squares system as these rows are revealed.  In prior work, this problem has been studied in the context of space complexity, where tools like leverage score sampling minimize the number of rows an algorithm must store. In this work, the focus shifts to time complexity, where the authors focus on making those leverage score sampling algorithm more computationally efficient.  The algorithm provided is matched with a compelling conditional lower bound.  Some experiments are provided.