# DYNAMIC LEAST-SQUARES REGRESSION

**Anonymous authors**
Paper under double-blind review

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 R[t][⇥][d] _⇥R[t]_ evolve in an online fashion (t ≫ _d), and the goal is to ef-_
ficiently maintain an (approximate) solution of minx(t) kA[(][t][)]x[(][t][)] _−_ **b[(][t][)]k2 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(d[1+][o][(1)]), almost matching the running time of the static (sketching-based) so-_
lution. By contrast, for exact (or 1/ poly(n)-accuracy) solutions, we show a separation between the models, namely, that dynamic LSR requires ⌦(d[2][−][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.


1 INTRODUCTION

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 R[n][⇥][d], b 2 R[n],_
the goal is to find the solution vector x that minimizes the mean squared error (MSE)

min (1)

**x** R[n][ k][Ax][ −] **[b][k][2][.]**
_2_

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 high_accuracy 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)[−][1]A[>]b, which requires O(nd[2]) time to compute using naive matrixmultiplication, or O(nd[!][−][1]) ⇡ _O(nd[1][.][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

**Ax** **b** 2 (1 + ✏) min
_k_ _−_ _k_ __ **x2R[d][ k][Ax][ −]** **[b][k][2][.]**

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 _O(nnz(A)✏[−][1]_ + d[!]) [1], which is near optimal in the “low

1In this paper we use _O(·) to hide polylogarithmic terms, and we use[e]_ _O✏(·) to hide poly(log d, ✏[−][1]) terms._
e


-----

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 _O(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[e]

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 solu-_
tion 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(d[2]) 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(d[2]) 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 small-_
_approximate 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.

1.1 OVERVIEW OF OUR RESULTS

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 ⌦(d[2][−][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 ⌦(Td[2]) by Theorem 1.1. Indeed, the impli-_
cation 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-by-_

column (ABi), requires d[3][−][o][(1)] time, in sharp contrast to the batch setting where this can be done using FMM
in d[!] _⌧_ _d[3]_ 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 _O(nnz(A[(][T][ )])✏[−][1]_ +

_d[!]) =_ _O(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:[e]
**Theorem 1.2[e]** (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 ) + ✏[−][6]d[3] log[5](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(d[1+][o][(1)]).

1.2 RELATED WORK

**Sketching and sampling** The least squares regression as a fundamental problem has been exten
sively 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 _O(nnz(A) log(1/✏) + d[!]) time (Clarkson & Woodruff,_

2017). See Appendix A for detailed discussions.

**Regression in online, streaming, and sliding window models[e]** 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 ⌦(d[2]) (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 _O(d[2]✏[−][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

[e]

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.

2 PROBLEM FORMULATION


In a dynamic least-squares regression problem, initially, we are given a matrix A[(0)] _2 R[n][0][⇥][d]_

together with a vector b[(0)] _2 R[n][0]_ . At the t-th step, a new data of form ((a[(][t][)])[>], β[(][t][)]) 2 R[d] _⇥_ 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 parame-_
_ters. 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 R[d][+1]_ _in the_

_preprocessing phase._

-  For each iteration t 2 [T ], the algorithm receives updates a[(][t][)] _2 R[d]_ _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 R[d][+][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 R[d]:_

**A[(][t][)]x[(][t][)]** **b[(][t][)]** 2 (1 + ✏) min
_k_ _−_ _k_ __ **x2R[d][ k][A][(][t][)][x][ −]** **[b][(][t][)][k][2][.]**


-----

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[(]i,[T]⇤[ )][k][2][ ]_ _[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 σ_ (0, 1).
_≥_ _2_
**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(d[3]) time, and this is within our desired total_
_running time of O✏(nnz(A[(][T][ )]) + d[3]). 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._

3 DYNAMIC ✏-LSR DATA STRUCTURE

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 ex-**
_ists 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 )+ _✏[−][6]d[3]_ log[5](TD/σ) _. Our_

_data structure uses at most O_ _✏[−][2]d[2]_ log[2](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 R[n][⇥][d], i 2 [n] we define its leverage score
_⌧_ (A) 2 R[n] as ⌧i(A) := a[>]i [(][A][>][A][)][†][a][i][.][ We define the generalized leverage score (same as Cohen]

et al. (2015)) of A 2 R[n][⇥][d] with respect to another matrix B 2 R[n][0][⇥][d] as : ⌧i[B][(][A][) :=][ a]i[>][(][B][>][B][)][†][a][i][.]

For more properties of the leverage scores see Section B.1.

3.1 TECHNIQUE OVERVIEW

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

_⌧d[M]+[(]t[t]+1[−][1)]_ [(][M][(][t][)][) = (][m][(][t][)][)][>][((][M][(][t][−][1)][)][>][M][(][t][−][1)][)][−][1][m][(][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(d[2]) 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-embedding[3] 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(d[2]) time to update the maintained matrices using Woodburry identity.
Hence, the amortized update time of the sampled rows is O✏(d[3]/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 Wood-_
_burry identity (O(d[2]) time per iteration). A sketching-based solution maintains a sketched matrix_
**SM 2 R[O][✏][(][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 **B[(][t][−][1)]m[(][t][−][1)]** 2. To efficiently compute the leverage score of a new row,
_k_ _k_
we left-multiply by a JL matrix J[(][t][−][1)] and maintain a proxy **B[(][t][−][1)]** = J[(][t][−][1)]B[(][t][−][1)], with the

guarantee that kB[e] [(][t][−][1)]mk2 ⇡kB[(][t][−][1)]mk2 for any m 2 R[d][+1] with high probability.

[e]

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 Uminx2Rd kD[(][t][)]A[(][t][)]x − **DPDATE[(][t][)]b[(][t]M[)]kEMBERS2, which has the closed-form solution of).** To obtain the final solution x x[(][t][)][(][t]2[)] R= G[d], we solve[(][t][)] _· u[(][t][)]_

and can be efficiently maintained by taking G[(][t][)] = ((A[(][t][)])[>](D[(][t][)])[2]A[(][t][)])[−][1] _2 R[d][⇥][d]_ and
**u[(][t][)]** = (A[(][t][)])[>](D[(][t][)])[2]b[(][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 R[s][⇥][(][d][+1)]_

5: H ((N)[>]N)[−][1] _. H 2 R[(][d][+1)][⇥][(][d][+1)]_

6: B **N · H** _. B 2 R[s][⇥][(][d][+1)]_

# JL approximation

7: δ _O(1/T_ [2]), k _O(✏[−][2]_ log(T/δ))
8: J JL(s, ✏, δ, T ) . JL matrix J 2 R[k][⇥][s]

9: **B** **J · B** _._ **B 2 R[k][⇥][(][d][+1)]**

# Maintain solution

10: G (A[>]D[>]DA)[−][1] _. G_ R[d][⇥][d]

[e] [e] _2_

11: u **A[>]D[2]b** _. u 2 R[d]_

12: x **G · u** _. x 2 R[d]_

**Algorithm 2 UPDATE (a, β)**


**Algorithm 3 SAMPLE (m)**

1: ⌧ **B** **m** 2
_k_ [e] _·_ _k[2]_

2: p min{3(1 + ✏)[2]✏[−][2]⌧ log(1/δ), 1}
3: ⌫ 1/[p]p with probability p, and ⌫ 0

otherwise


**Algorithm 4 UPDATEMEMBERS (m)**

# Update spectral approximation

1: s _s + 1_
2: ∆H 1+m[>]Hm/p
_−_ **[Hmm][>][H][/p]**

3: H **H + ∆H**
4: B [(B + N · ∆H)[>], H · m/[p]p][>]

5: N [N[>], m/[p]p][>]

# Update JL approximation

6: J JL(s, ✏, δ, T )
7: **B** **J · B**
# Update solution

8: G[e] **G** 1+a[>]Ga/p
_−_ **[Gaa][>][G][/p]**

9: u **u + β · a/p**

10: x **G · u**


1: m [a[>], β][>] _. m 2 R[d][+1]_

2: ⌫ SAMPLE(m) _. ⌫_ _2 R_

**D** 0

3: D

0 _⌫_

 %

4: if ⌫ = 06 **then UPDATEMEMBERS(m)**
5: return x


We outline the proof of Theorem 3.1, and defer the detailed proof to Appendix C due to space limits.


-----

3.2 CORRECTNESS

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), R[s][(][t][)][⇥][(][d][+1)], where S[(][t][)] [d + t + 1] is defined as the set of
_⇤_ _2_ _⇢_

_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 R[s][(][t][)][⇥][(][d][+1)]._

_6._ **B[(][t][)]** = J[(][t][)] _· B[(][t][)]_ _2 R[k][⇥][(][d][+1)], where k = O(✏[−][2]_ log(T/δ)).

_7. G[(][t][)]_ = (A[(][t][)])[>](D[(][t][)])[2]A[(][t][)][#][−][1] R[d][⇥][d].

[e] _2_

_8. u[(][t][)]_ = ("A[(][t][)])[>](D[(][t][)])[2]b[(][t][)] _2 R[d]._

_9. x[(][t][)]_ = (A[(][t][)])[>](D[(][t][)])[2]A[(][t][)][#][−][1] _· (A[(][t][)])[>](D[(][t][)])[2]b[(][t][)]_ _2 R[d]._

The following lemma is key for our correctness analysis. It shows that we maintain a good approxi-"
mation 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]⌧d[M]+[(]t[t]+1[−][1)]_ [(][M][(][t][)][)][ ] _[⌧]_ [(][t][)][ ] [(1 +][ ✏][)][2][⌧]d[M]+[(]t[t]+1[−][1)] [(][M][(][t][)][)][,][ 8][t][ 2][ [][T] []][,] (2)

_and_

(M[(][t][)])[>](D[(][t][)])[2]M[(][t][)] _✏_ (M[(][t][)])[>]M[(][t][)], _t_ [0 : T ]. (3)
_⇡_ _8_ _2_

_Proof Sketch. We prove by induction and show that with probability 1 −_ 2tδ, Eq. (3) holds for all
_t[0]_ [0 : t] and Eq. (2) holds for all t[0] [t]. The base case t = 0 holds trivially, as D[(0)] = Id+1,
_2_ _2_
and therefore, (M[(0)])[>](D[(0)])[2]M[(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 − _✏)⌧d[M]+[(]t[t]+1[−][1)]_ [(][M][(][t][)][)][ k][B][(][t][−][1)][ ·][ m][(][t][)][k]2[2] _[]_ [(1 +][ ✏][)][⌧]d[M]+[(]t[t]+1[−][1)] [(][M][(][t][)][)][.]

-  We then use the JL lemma (Lemma B.4) to show that with probability 1 − _δ, the sketched_

covariance matrix kB[e] [(][t][−][1)]m[(][t][)]k returns good estimation kB[(][t][−][1)]m[(][t][)]k. That is

(1 − _✏)[2]_ _· ⌧d[M]+[(]t[t]+1[−][1)]_ [(][M][(][t][)][)][ k]B[e] [(][t][−][1)] _· m[(][t][)]k2[2]_ _[]_ [(1 +][ ✏][)][2][ ·][ ⌧]d[M]+[(]t[t]+1[−][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][)])[2]M[(][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 R[d]_ _such that_

**A[(][t][)]x[(][t][)]** **b[(][t][)]** 2 (1 + ✏) min
_k_ _−_ _k_ __ **x2R[d][ k][A][(][t][)][x][ −]** **[b][(][t][)][k][2][.]**


-----

3.3 TIME ANALYSIS

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][)]_ = 06, 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][)] = 06 _, then UPDATE takes O_ _✏[−][2]s[(][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_


_⌧d[M]+[(]t[t]+1[−][1)]_ [(][M][(][t][)][)][ ] _[O][(][d][ log(][TD/σ][))][.]_


_t=1_


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 ) + ✏[−][6]d[3] log[5](TD/σ) _._

" #

_Proof Sketch. In this proof sketch we simplify the second term as d[3]_ _· 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(


_p[(][t][)]) = O(✏[−][2]_ _· log(1/δ) ·_


_d · poly(✏[−][1]_ log(TD/σ))


_⌧_ [(][t][)])  _O_


_i=1_


_t=1_



-  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_


_d[3]_ _· poly(✏[−][1]_ log(TD/σ))


_s[(][T][ )]d · poly(✏[−][1]_ log(TD/σ))


= O


4 HARDNESS RESULT

We prove a ⌦(d[2][−][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 R[d]_ _at every time step. The conjectures states that there is no al-_
_gorithm that solves the OMv problem using poly(d) preprocessing time and O(d[2][−][γ]) 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** N, T =
1 _2_
poly(d), ✏ = _d[8]T_ [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 ⌦(d[2][−][γ]) 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 R[d][⇥][d]_ _be a symmetric matrix whose eigenvalues satisfy 1 _ _λd(H) · · · _
_λOMv conjecture is true, then there is no algorithm with1(H) _ 3. For any t 2 [T ], z[(][t][)] _2 R[d]_ _is revealed at the poly( t-th step, andd) preprocessing time and kz[(][t][)]k2 _ 1. Assuming the O(d[2][−][γ])
_amortized running time that can return an O(1/d[2])-approximate answer to Hz[(][t][)]_ _for all t, i.e., a_
_vector y[(][t][)]_ _2 R[d]_ _s.t. ky[(][t][)]_ _−_ **Hz[(][t][)]k2 ** _✏, and has an error probability _ 1/3.

_Proof Sketch. Given a Boolean matrix2Id_ _d1_ **[B]** **B 2 {0, 1}[d][⇥][d]** in the OMv conjecture, we construct a PSD

matrix H = 1 R[2][d][⇥][2][d]. We note that H is symmetric and 1 _λd(H)_ _λ1(H)_ 3.

_d_ **[B][>]** 2Id _2_ __ __ __

 %

Given a binary OMv query vector z[(][t][)], we construct z[(][t][)] = (0d, z[(][t][)]) 2 R[2][d]. Since H has a constant
condition number, we can prove that rounding an ✏ 1/d[2]-approximate answer **y** _✏_ **H** **z[(][t][)]** still
_⇠_ _⇡_ _·_

gives the correct binary answer to Bz[(][t][)].

b

**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 each1_ _t 2 [T_ ], the incoming row? abe the optimal solution at the[(][t][)] is a small scaled version of t-th step, z[(][t][)], i.e, H a[(][(][t][t][)][)] :==

_d[2][p]T_

_[·][ z][(][t][)][ 2][ R][d][, and the label is][ 1][. Let][ x][(][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][)] 1

_[±][ O][(]_ _d[4][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

_d[4][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

_d[6][p]T_ [)][.]


In particular, let y[(][t][)] = d[2][p]T (x[(][t][)] **x[(][t][−][1)]), Step 2 and 3 directly implies** **y** **Hz[(][t][)]** 2

_O(1/d[2]). This completes the proof._ _−_ _k_ _−_ _k_ __

5 EXPERIMENTS

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 R[d][⇥][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 R[d]_ 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 (= _pT_ ) 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 Repository[4].

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.

(a) Synthetic dataset (b) VirusShare

Figure 1: Experiment results. The x-axis shows the running time (unit: seconds), and the y-axis
shows the relative error (err/errstd 1), where err is the error of the particular approach, and errstd
is the error of the static Normal equation. The − _y-axis is on a log scale. For uniform sampling, we_
take sampling probability p = 0.05, 0.1, 0.2, 0.5. For row sampling and our algorithm, we take the
error parameter ✏ = 0.1, 0.2, 0.5, 1. Kalman’s approach has a relative error of 0.

6 CONCLUSION

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 ⌧ _d[2][−][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 d[3] 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](https://archive.ics.uci.edu/ml/datasets.php)


-----

REFERENCES

Thomas D Ahle, Michael Kapralov, Jakob BT Knudsen, Rasmus Pagh, Ameya Velingker, David P

Woodruff, and Amir Zandieh. Oblivious sketching of high-degree polynomial kernels. In Pro_ceedings of the Fourteenth Annual ACM-SIAM Symposium on Discrete Algorithms, pp. 141–160._
SIAM, 2020.

Nir Ailon and Bernard Chazelle. Approximate nearest neighbors and the fast johnson-lindenstrauss

transform. In Proceedings of the thirty-eighth annual ACM symposium on Theory of computing,
pp. 557–563, 2006.

Josh Alman and Virginia Vassilevska Williams. A refined laser method and faster matrix multipli
cation. In Proceedings of the 2021 ACM-SIAM Symposium on Discrete Algorithms (SODA), pp.
522–539. SIAM, 2021.

Haim Avron, Kenneth L Clarkson, and David P Woodruff. Faster kernel ridge regression using

sketching and preconditioning. SIAM Journal on Matrix Analysis and Applications, 38(4):1116–
1138, 2017.

Christos Boutsidis, David P Woodruff, and Peilin Zhong. Optimal principal component analysis in

distributed and streaming models. In Proceedings of the forty-eighth annual ACM symposium on
_Theory of Computing, pp. 236–249, 2016._

van den Jan Brand, Yin-Tat Lee, Danupon Nanongkai, Richard Peng, Thatchaphol Saranurak, Aaron

Sidford, Zhao Song, and Di Wang. Bipartite matching in nearly-linear time on moderately dense
graphs. In 2020 IEEE 61st Annual Symposium on Foundations of Computer Science (FOCS), pp.
919–930. IEEE, 2020a.

van den Jan Brand, Yin Tat Lee, Aaron Sidford, and Zhao Song. Solving tall dense linear programs

in nearly linear time. In Proceedings of the 52nd Annual ACM SIGACT Symposium on Theory of
_Computing (STOC), pp. 775–788, 2020b._

van den Jan Brand, Yin Tat Lee, Yang P Liu, Thatchaphol Saranurak, Aaron Sidford, Zhao Song, and

Di Wang. Minimum cost flows, mdps, and `1-regression in nearly linear time for dense instances.
In Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing (STOC),
pp. 859–869, 2021a.

van den Jan Brand, Binghui Peng, Zhao Song, and Omri Weinstein. Training (overparametrized)

neural networks in near-linear time. In 12th Innovations in Theoretical Computer Science Con_ference (ITCS 2021), 2021b._

Vladimir Braverman, Petros Drineas, Cameron Musco, Christopher Musco, Jalaj Upadhyay, David P

Woodruff, and Samson Zhou. Near optimal linear algebra in the online and sliding window

models. In 2020 IEEE 61st Annual Symposium on Foundations of Computer Science (FOCS), pp.
517–528. IEEE, 2020.

Sébastien Bubeck. Convex optimization: Algorithms and complexity. Foundations and Trends in

_Machine Learning, 8(3-4):231–357, 2015._

Nadiia Chepurko, Kenneth L Clarkson, Praneeth Kacham, and David P Woodruff. Near
optimal algorithms for linear algebra in the current matrix multiplication time. arXiv preprint
_arXiv:2107.08090, 2021._

Charles K. Chui. Estimation, control, and the discrete kalman filter (donald e. calin). SIAM Re
_[view, 32(3):493–494, 1990. doi: 10.1137/1032097. URL https://doi.org/10.1137/](https://doi.org/10.1137/1032097)_
[1032097.](https://doi.org/10.1137/1032097)

Kenneth L Clarkson and David P Woodruff. Numerical linear algebra in the streaming model. In

_Proceedings of the forty-first annual ACM symposium on Theory of computing, pp. 205–214,_
2009.

Kenneth L Clarkson and David P Woodruff. Low-rank approximation and regression in input spar
sity time. Journal of the ACM (JACM), 63(6):1–45, 2017.


-----

Michael B Cohen, Yin Tat Lee, Cameron Musco, Christopher Musco, Richard Peng, and Aaron

Sidford. Uniform sampling for matrix approximation. In Proceedings of the 2015 Conference on
_Innovations in Theoretical Computer Science (ITCS), pp. 181–190. ACM, 2015._

Michael B Cohen, Cameron Musco, and Jakub Pachocki. Online row sampling. Theory OF Com
_puting, 16(15):1–25, 2020._

Edgar Dobriban and Sifan Liu. Asymptotics for sketching in least squares. In Proceedings of the

_33rd International Conference on Neural Information Processing Systems, pp. 3675–3685, 2019._

Petros Drineas, Michael W. Mahoney, and S. Muthukrishnan. Sampling algorithms for l 2 regression

and applications. In Proceedings of the seventeenth annual ACM-SIAM symposium on Discrete
_algorithm (SODA), pp. 1127–1136, 2006a._

Petros Drineas, Michael W. Mahoney, and S. Muthukrishnan. Subspace sampling and relative
error matrix approximation: Column-based methods. In Approximation, Randomization, and
_Combinatorial Optimization. (APPROX-RANDOM), pp. 316–326, 2006b._

Petros Drineas, Michael W. Mahoney, and S. Muthukrishnan. Subspace sampling and relative
error matrix approximation: Column-row-based methods. In European Symposium on Algorithms
_(ESA), pp. 304–314, 2006c._

Trevor Hastie, Jerome H. Friedman, and Robert Tibshirani. The Elements of Statistical Learning:

_Data Mining, Inference, and Prediction. Springer Series in Statistics. Springer, 2001. ISBN 978-_
1-4899-0519-2. doi: 10.1007/978-0-387-21606-5. [URL https://doi.org/10.1007/](https://doi.org/10.1007/978-0-387-21606-5)

[978-0-387-21606-5.](https://doi.org/10.1007/978-0-387-21606-5)

Elad Hazan. Introduction to online convex optimization. arXiv preprint arXiv:1909.05207, 2019.


Monika Henzinger, Sebastian Krinninger, Danupon Nanongkai, and Thatchaphol Saranurak. Unify
ing and strengthening hardness for dynamic problems via the online matrix-vector multiplication
conjecture. In Proceedings of the forty-seventh annual ACM symposium on Theory of computing,
pp. 21–30, 2015.

Shunhua Jiang, Zhao Song, Omri Weinstein, and Hengjie Zhang. A faster algorithm for solving

general lps. In Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing
_(STOC), pp. 823–832, 2021._

William B Johnson and Joram Lindenstrauss. Extensions of lipschitz mappings into a hilbert space.

_Contemporary mathematics, 26(189-206):1, 1984._

Rudolph Emil Kalman. A new approach to linear filtering and prediction problems. Journal of Basic

_Engineering, 82(1):35–45, 1960._

François Le Gall. Powers of tensors and fast matrix multiplication. In Proceedings of the 39th

_international symposium on symbolic and algebraic computation, pp. 296–303, 2014._

Yin Tat Lee and Aaron Sidford. Path finding methods for linear programming: Solving linear

programs in _O(prank) iterations and faster algorithms for maximum flow. In 2014 IEEE 55th_

_Annual Symposium on Foundations of Computer Science, pp. 424–433. IEEE, 2014._

Yin Tat Lee, Zhao Song, and Qiuyi Zhang. Solving empirical risk minimization in the current matrix[e]

multiplication time. In Annual Conference on Learning Theory (COLT), 2019.

Edo Liberty, Zohar Karnin, Bing Xiang, Laurence Rouesnel, Baris Coskun, Ramesh Nallapati, Julio

Delgado, Amir Sadoughi, Yury Astashonok, Piali Das, et al. Elastic machine learning algorithms
in amazon sagemaker. In Proceedings of the 2020 ACM SIGMOD International Conference on
_Management of Data, pp. 731–737, 2020._

Aleksander Madry. Navigating central path with electrical flows: From flows to matchings, and

back. In 2013 IEEE 54th Annual Symposium on Foundations of Computer Science (FOCS), pp.
253–262. IEEE, 2013.

Nick Martin and Hermine Maes. Multivariate analysis. Academic press London, 1979.


-----

Xiangrui Meng and Michael W Mahoney. Low-distortion subspace embeddings in input-sparsity

time and applications to robust linear regression. In Proceedings of the forty-fifth annual ACM
_symposium on Theory of computing, pp. 91–100, 2013._

Jelani Nelson and Huy L Nguyên. Osnap: Faster numerical linear algebra algorithms via sparser

subspace embeddings. In 2013 ieee 54th annual symposium on foundations of computer science,
pp. 117–126. IEEE, 2013.

Rasmus Pagh. Compressed matrix multiplication. _ACM Transactions on Computation Theory_

_(TOCT), 5(3):1–17, 2013._

German I Parisi, Ronald Kemker, Jose L Part, Christopher Kanan, and Stefan Wermter. Continual

lifelong learning with neural networks: A review. Neural Networks, 113:54–71, 2019.

R. L. Plackett. Some theorems in least squares. _Biometrika, 37(1/2):149–157, 1950._ ISSN

[00063444. URL http://www.jstor.org/stable/2332158.](http://www.jstor.org/stable/2332158)

Lawrence R Rabiner and Bernard Gold. Theory and application of digital signal processing. Engle
_wood Cliffs: Prentice-Hall, 1975._

Ilya Razenshteyn, Zhao Song, and David P Woodruff. Weighted low rank approximations with

provable guarantees. In Proceedings of the forty-eighth annual ACM symposium on Theory of
_Computing, pp. 250–263, 2016._

Vladimir Rokhlin and Mark Tygert. A fast randomized algorithm for overdetermined linear least
squares regression. Proceedings of the National Academy of Sciences, 105(36):13212–13217,
2008.

Tamas Sarlos. Improved approximation algorithms for large matrices via random projections. In

_2006 47th Annual IEEE Symposium on Foundations of Computer Science (FOCS’06), pp. 143–_
152. IEEE, 2006.

Daniel A Spielman and Nikhil Srivastava. Graph sparsification by effective resistances. _SIAM_

_Journal on Computing, 40(6):1913–1926, 2011._

Daniel A Spielman and Shang-Hua Teng. Nearly-linear time algorithms for graph partitioning,

graph sparsification, and solving linear systems. In Proceedings of the thirty-sixth annual ACM
_symposium on Theory of computing, pp. 81–90, 2004._

Stephen M. Stigler. Gauss and the Invention of Least Squares. _The Annals of Statistics, 9(3):_

[465–474, 1981. doi: 10.1214/aos/1176345451. URL https://doi.org/10.1214/aos/](https://doi.org/10.1214/aos/1176345451)
[1176345451.](https://doi.org/10.1214/aos/1176345451)

Volker Strassen. Gaussian elimination is not optimal. Numerische mathematik, 13(4):354–356,

1969.

David Woodruff. A very sketchy talk (invited talk). In 48th International Colloquium on Au
_tomata, Languages, and Programming (ICALP 2021). Schloss Dagstuhl-Leibniz-Zentrum fuer_
Informatik, 2021.

David P. Woodruff. Sketching as a tool for numerical linear algebra. Foundations and Trends in

_Theoretical Computer Science, 10(1-2):1–157, 2014._


-----

