Solving Tall Dense Linear Programs in Nearly Linear Time

Introduction

is one of the most fundamental and well-studied problems in computer science and optimization. Developing faster algorithms for (1.1) has been the subject of decades of extensive research and the pursuit of faster linear programming methods has lead to numerous algorithmic advances and the advent of fundamental optimization techniques, e.g. simplex methods [Dan51], ellipsoid methods [Kha80], and interior-point methods (IPMs) [Kar84].

With this IPM in hand, the problem of achieving our desired running time reduces to implementing this IPM efficiently. This problem is that of maintaining multiplicative approximations to vectors (the current iterates), leverage scores (a measure of importance of the rows under local rescaling), and the inverse of matrix (the system one needs to solve to take a step of the IPM) under small perturbations. While variants of vector maintenance have been considered recently [CLS19, LSZ19, Bra20] and inverse maintenance is well-studied historically [Kar84, NN89, Vai89a, NN91, LS14, LS15, CLS19, LSZ19, LS19, Bra20], none of these methods can be immediately applied in our setting where we cannot afford to pay too much in terms of $n$ each iteration.

Our second contribution is to show that these data-structure problems can be solved efficiently. A key technique we use to overcome these issues is heavy-hitters sketching. We show that it is possible to apply a heavy hitters sketch (in particular [KNPW11, Pag13]) to the iterates of the method such that we can efficiently find changes in the coordinates. This involves carefully sketching groups of updates and dynamically modifying the induced data-structure. These sketches only work against non-adaptive adversaries, and therefore care is need to ensure that the sketch is used only to save time and not affect the progression of the overall IPM in a way that break this non-adaptive assumption. To achieve this we use the sketches to propose short-lists of possible changes which we then filter to ensure that the output of the datastructure is deterministic (up to a low probability failure event). Coupling this technique with known Johnson-Lindestrauss sketches yields our leverage score maintenance data structure and adapting and simplifying previous inverse maintenance techniques yields our inverse maintenance data structure. We believe this technique of sketching the central path is powerful and may find further applications.

See [Ren88, LS14] on the discussion on converting such an approximate solution to an exact solution. For integral $\mathbf{A},b,c$, it suffices to pick $\delta=2^{-O(L)}$ to get an exact solution where $L=\log(1+d_{\max}+\|c\|_{\infty}+\|b\|_{\infty})$ is the bit complexity and $d_{\max}$ is the largest absolute value of the determinant of a square sub-matrix of $\mathbf{A}$. For many combinatorial problems $L=O(\log(n+\|b\|_{\infty}+\|c\|_{\infty}))$.

2 Previous Work

Linear programming has been the subject of extensive research for decades and it is impossible to completely cover to this impressive line of work in this short introduction. Here we cover results particularly relevant to our approach. For more detailed coverage of prior-work see, e.g. [NN94, Ye11].

IPMs: The first proof of a polynomial time IPM was due to Karmarkar in [Kar84]. After multiple running time improvements [Kar84, NN89, Vai89a, NN91] the current fastest IPMs are the aforementioned results of [LS19] and [CLS19]. Beyond these results, excitingly the work of [CLS19] was recently extended to obtain comparable running time improvements for solving arbitrary empirical risk minimization problems (ERM) [LSZ19] and was recently simplified and de-randomized by van den Brand [Bra20]. These works consider variants of the vector maintenance problem and our work is inspired in part by them. Our IPM leverages the barrier from [LS19] in a new way, that enables the application of robustness techniques from [LS19, CLS19, Bra20] and new techniques for handling approximately feasible points.

Leverage Scores and Lewis Weights: Leverage scores [SS11, Mah11, LMP13, CLM+15] (and more broadly Lewis weights [Lew78, BLM89, CP15, LS19]) are fundamental notions of importance of rows of a matrix with numerous applications. In this paper we introduce a natural online problem for maintaining multiplicative approximations to leverage scores and we show how to solve this problem efficiently. Though we are unaware of this problem being studied previously, we know that in the special case of leverage scores induced by graph problems, which are known as effective resistances, there are dynamic algorithms for maintaining them, e.g [DGGP19]. However, these algorithms seem tailored to graph structure and it is unclear how to apply them in our setting. Further, there are streaming algorithms for variants of this problem [AGM13, KLM+17, KNST20], however their running time is too large for our purposes.

Overview of Approach

We first give a quick introduction to primal-dual path IPMs, review the central path from [LS19] and from it derive the central path that our method is based on, and explain our new IPM.

The Lewis Weight Barrier and Beyond

where $\mathbf{W}=\mathbf{Diag}(w)$. In the special case when $p=2$,

is known as the leverage scores of the rows of $\mathbf{A}$ and is a fundamental object for dimension reduction and solving linear systems [SS11, Mah11, LMP13, NW14, Woo14, CLM+15, SWZ19].

Formally, in this paper we consider the following regularized-variant of this centrality measure:

Our Robust Primal-Dual Method:

The main guarantees of this new IPM are given by Theorem 5 (Section 4.1) and proven in Section 5 and Section B. This theorem formalizes the above discussion and quantifies how much the iterates, leverage scores, and Hessian can change in each iteration. These bounds are key to obtaining an efficient method that can maintain multiplicative approximations to these quantities.

2 Heavy Hitters, Congestion Detection, and Sketching the Central Path

We treat each of these problems as a self contained data-structure problem, the first we call the vector maintenance problem, the second we call the leverage score maintenance problem, and the third has been previously studied (albeit different variants) and is called the inverse maintenance problem. For each problem we build efficient solutions by combining techniques form the sketching literature (e.g. heavy hitters sketches and Johnson-Lindenstrauss sketches) and careful potential functions and tools for dealing with sparse and low rank approximations. (Further work is also needed to maintain the gradient of a potential used to measure proximity to the central path (Section D) and this is discussed in the appendix.)

In the remainder of this overview, we briefly survey how we solve each of these problems. A common issue that needs to addressed in solving each problem is that of hiding randomness and dealing with adversarial input. Each data-structure uses sampling and sketching techniques to improve running times. While these techniques are powerful and succeed with high probability, they only work against an oblivious adversary, i.e. one which provides input that does not depend on the randomness of the data structure. However, the output of our data-structures are used to take steps along the central path and provide the next input, so care needs to be taken to argue that the output of the data structure, as used by the method, doesn’t somehow leak information about the randomness of the sketches and samples into the next input. A key contribution of our work is showing how to overcome this issue in each case.

So far we only explained how we can detect large changes in $y^{(t)}$ that occur within a single iteration. However, it could be that some entry changes slowly over several iterations. It is easy to extend our heavy hitter technique to also detect these slower changes. The idea is that we not just detect changes within one iteration (i.e. $y_{i}^{(t)}\not\approx_{\epsilon/2}y_{i}^{(t-1)}$) but also changes within any power of two (i.e. $y_{i}^{(t)}\not\approx_{\epsilon/2}y_{i}^{(t-2^{i})}$for $i=1,...,\log t$). To make sure that this does not become too slow, we only check every $2^{i}$ iterations if there was a large change within the past $2^{i}$ iterations. One can prove that this is enough to also detect slowly changing entries.

Leverage Score Maintenance:

Inverse Maintenance:

Fortunately, we note that we do not need a sparsifier that satisfies both conditions at all time. Therefore, our final data structure has two ways to output the inverse of the sparsifier. The sparsifier that works only against oblivious adversary is used in implementing the Newton steps and the sparsifier that is slowly changing is used to compute the sketch $\mathbf{M}\mathbf{J}$ used in the leverage score maintenance problem mentioned above.

Preliminaries

Here we discuss varied notation we use throughout the paper. We adopt similar notation to [LS14] and some of the explanations here are copied directly form this work.

Approximations: We use $x\approx_{\epsilon}y$ to denote that $\exp(-\epsilon)y\leq x\leq\exp(\epsilon)y$ and $\mathbf{A}\approx_{\epsilon}\mathbf{B}$ to denote that $\exp(-\epsilon)\mathbf{B}\preceq\mathbf{A}\preceq\exp(\epsilon)\mathbf{B}$.

Linear Programming Algorithm

The proof for Theorem 1 uses four intermediate results, which we formally state in the next Sections 4.1, to 4.4. Each of these intermediate results is self-contained and analyzed in its own section, Sections 5 to 8 correspondingly. In this section we show how these results can be combined to obtain Theorem 1. The first result is a new, improved IPM as outlined in Section 2.1. The exact statement is given in Section 4.1 and its analysis can be found in Section 5 and Section B. Here we give a rough summary to motivate the other three results used by our linear programming algorithm.

(ii) In Section 4.3 we present a data-structure that can main approximate leverage scores, required by the IPM. Its correctness is proven in Section 7.

With this we have all tools available for proving the main result Theorem 1. In Section 4.5 we show how to combine all these tools to obtain the fast linear programming algorithm.

Furthermore, throughout Algorithm 1, we have

$xs\approx_{4\epsilon}\mu\cdot\tau(x,s)$ for some $\mu$ where $(x,s)$ is immediate points in the algorithms

Note that the complexity of Theorem 5 depends on the cost of implementing Lines 1, 1, 1, and 1 of Algorithm 1, as well as the cost of Algorithm 10. Here Lines 1, 1 and 1 ask us to maintain an approximation of the primal dual solution pair $(x,s)$. A data-structure for this task is presented in Section 4.2. Additionally, to compute Lines 1 and 1, we must have access to an approximate inverse of $\mathbf{A}^{\top}\overline{\mathbf{S}}^{-1}\overline{\mathbf{X}}\mathbf{A}$ (see Line 1). The task of maintaining this inverse will be performed by the data-structure presented in Section 4.4. At last, consider Line 1. To implement this line, we must have an approximation of the leverage scores $\tau(x,s)$. In Section 4.3, we present a data-structure that can efficiently maintain such an approximation.

To help us analyze the cost of Algorithm 10, we prove the following in Section B.

2 Vector Data Structure

To maintain an approximation of $x^{(t)}$, we can then use the data-structure for maintaining an approximation of $y^{(t)}$ by choosing

Likewise, we can maintain an approximation of $s^{(t)}$ by a slightly different choice of parameters. The exact result, which we prove in Section 6, is the following Theorem 8:

There exists a Monte-Carlo data-structure (Algorithm 6), that works against an adaptive adversary, with the following procedures:

3 Leverage Score Maintenance

The IPM of Theorem 5, requires approximate leverage scores of some matrix of the form $\mathbf{G}\mathbf{A}$, where $\mathbf{G}$ is a diagonal matrix (see Line 1 of Algorithm 1, where $\mathbf{G}=(\overline{\mathbf{X}}/\overline{\mathbf{S}})^{1/2}$) . Here the matrix $\mathbf{G}$ changes slowly from one iteration of the IPM to the next one, which allows us to create a data-structure that can maintain the scores more efficiently than recomputing them from scratch every time $\mathbf{G}$ changes. In Section 7 we prove the following result for maintaining leverage scores:

There exists a Monte-Carlo data-structure (Algorithm 7), that works against an adaptive adversary, with the following procedures:

where $T_{\Psi}$ is the time required to multiply a vector with $\Psi^{(t)}$ (i.e. in case it is given implicitly via a data structure).

4 Inverse Maintenance

The amortized time per call of Update is $O(n+\epsilon^{-2}\cdot(d^{\omega-\frac{1}{2}}+d^{2})\cdot\log^{3/2}(n))$.

The time per call of Solve is $O(n+\delta^{-2}\cdot d^{2}\cdot\log^{2}(n/\delta))$.

Theorem 10 holds even if the input $w$ and $\widetilde{\tau}$ of the algorithm depends on the output of Solve. Furthermore, we have

5 Linear Programming Algorithm

Moving along these two paths of the first and second phase is performed via the IPM of Algorithm 1 (Section 4.1, Theorem 5). Note that Algorithm 1 does not specify in Line 1 how to obtain the approximate solution pair $(\overline{x},\overline{s})$, so we must implement this step on our own. Likewise, we must specify how to efficiently compute the steps in Lines 1 and 1. These implementations can be found in the second part of Algorithm 2. The high-level idea is to use the data-structures presented in Section 4.2 to 4.4.

To illustrate, consider Line 1 of Algorithm 1, which computes

for $\overline{\mathbf{W}}=\overline{\mathbf{X}}\overline{\mathbf{S}}$ and $\mathbf{Q}=\overline{\mathbf{S}}^{-1/2}\overline{\mathbf{X}}^{1/2}\mathbf{A}\mathbf{H}^{-1}\mathbf{A}^{\top}\overline{\mathbf{X}}^{1/2}\overline{\mathbf{S}}^{-1/2}$ for $\mathbf{H}\approx_{\epsilon}\mathbf{A}^{\top}\overline{\mathbf{S}}^{-1}\overline{\mathbf{X}}\mathbf{A}$. We split this task into three parts: (i) compute $r:=\mathbf{A}^{\top}\overline{\mathbf{X}}^{1/2}\overline{\mathbf{S}}^{-1/2}\overline{\mathbf{W}}^{1/2}h$, (ii) compute $v:=\mathbf{H}^{-1}r$, and (iii) compute $(1-2\alpha)\overline{\mathbf{S}}\overline{\mathbf{W}}^{-1/2}\overline{\mathbf{S}}^{-1/2}\overline{\mathbf{X}}^{1/2}\mathbf{A}v=(1-2\alpha)\mathbf{A}v$. Part (i), the vector $r$, can be maintained efficiently, because we maintain the approximate solutions $\overline{x}$, $\overline{s}$ (thus also $\overline{w}$) and vector $h$, (by Algorithm 12, Theorem 59, of Section D) in such a way, that per iteration only few entries change on average. Part (ii) is solved by the inverse maintenance data-structure of Section 4.4 (Theorem 10). The last part (iii) is solved implicitly by the data-structure of Section 4.2 (Theorem 8), which is also used to obtain the approximate solutions $\overline{x},\overline{s}$ in Line 1 of Algorithm 1. We additionally run the data-structure of Section 4.3 (Theorem 9) in parallel, to maintain an approximation of the leverage scores, which allows us to find the approximation $\overline{v}$ required in Line 1. These modifications to Algorithm 1 are given in the second part of Algorithm 2.

The following theorem shows how to reduce solving any bounded linear program to solving a linear program with a non-degenerate constrain matrix and an explicit initial primal and dual interior point. This theorem is proven in Appendix C.

Consider linear program $\min_{\mathbf{A}^{\top}x=b,x\geq 0}c^{\top}x$ with $n$ variables and $d$ constraints. Assume that 1. Diameter of the polytope : For any $x\geq 0$ with $\mathbf{A}^{\top}x=b$, we have that $\|x\|_{2}\leq R$. 2. Lipschitz constant of the linear program : $\|c\|_{2}\leq L$. 3. The constraint matrix $\mathbf{A}$ is non-degenerate.

For any $\delta\in(0,1]$, the modified linear program $\min_{\overline{\mathbf{A}}^{\top}\overline{x}=\overline{b},\overline{x}\geq 0}\overline{c}^{\top}\overline{x}$ with

As outlined before, the initial points given in Theorem 12 do not satisfy

Similarly, we have $T(v)_{i}^{2/p}\geq e^{-|1-\frac{2}{p}|\alpha}T(w)_{i}^{2/p}$. Taking $p/2$ power of both sides, we have

Hence, $T(v)\approx_{|p/2-1|\alpha}T(w)$.

Consequently, let $w_{0}=\eta 1$ and consider the algorithm $w_{k+1}=T(w_{k})$. Since $\eta>0$ we have that

Now, we first prove the correctness of the Algorithm 2.

Algorithm 2 outputs $x$ such that w.h.p. in $n$

We define the primal dual point $(x,s)$ via the formula of lines 1 and 1. We start by showing that our implementation of Algorithm 1 in Algorithm 2 satisfies all required conditions, i.e. that we can apply Theorem 5. Throughout this proof, states hold only w.h.p. and therefore the restatement of this is often omitted for brevity.

We first show that throughout Centering, we have the invariant above, assuming the input parameter $\tau_{i}^{\text{(init)}}$ satisfied $\tau_{i}^{\text{(init)}}\approx_{\gamma/4}\sigma(\mathbf{S}^{-1/2-\alpha}\mathbf{X}^{1/2-\alpha}\mathbf{A})+\frac{d}{n}1$. Theorem 8 shows that $x^{\text{(tmp)}}\approx_{\gamma/8}x$ and $s^{\text{(tmp)}}\approx_{\gamma/8}s$ (Line 3 and 3). Then, by the update rule (Line 3 and 3), we have that $x^{\text{(tmp)}}\approx_{\gamma/8}\overline{x}$ and $s^{\text{(tmp)}}\approx_{\gamma/8}\overline{s}$. Hence, we have the desired approximation for $\overline{x}$ and $\overline{s}$.

$\Psi^{(\alpha)}\approx_{\gamma/(512\log n)}(\mathbf{A}^{\top}\mathbf{S}^{-1-2\alpha}\mathbf{X}^{1-2\alpha}\mathbf{A})^{-1}$ (Line 2) and

(Line 2 and Line 3). Again, by update rule (Line 3), we have the desired approximation for $\overline{\tau}$.

Invariant: x​s≈μ⋅τ⁡(x,s)xs\approx\mu\cdot\tau(x,s):

Initially, $x=1$ due to the reduction (Lemma 12), $\mu=1$ and $s\approx_{\epsilon}\sigma(\mathbf{S}^{-1/2-\alpha}\mathbf{A})+\frac{d}{n}1$ (Line 2). Hence, $xs\approx_{\epsilon}\mu\cdot(\sigma(\mathbf{S}^{-1/2-\alpha}\mathbf{X}^{1/2-\alpha}\mathbf{A})+\frac{d}{n}1)$ initially.

First, we bound the denominator. We have that $sx\approx_{1/2}\mu\tau$, so for all $i\in[n]$ we have

where we used $\|x\|_{\infty}\leq O(n)$ by Lemma 12 as we ensure $x$ is sufficiently close to feasible for the modified linear program by Theorem 5. For the numerator, we note that $\|c\|_{\infty}\leq 1$ for the modified linear program and $\|c^{\text{(tmp)}}\|_{\infty}=\|s\|_{\infty}\leq 3$ for the $s$ computed by Theorem 13 in Line 2.again by the definition of $s$ (Line 2) and the modified $\mathbf{A}$ and $y$. Hence, we have that

for any constant $c$ by choosing the constant in the $O(\cdot)$ in Line 2 appropriately. Thus $xs\approx_{2\epsilon}\mu\cdot\tau(x,s)$ and $\overline{\tau}\approx_{\gamma/4}\sigma(\mathbf{S}^{-1/2-\alpha}\mathbf{X}^{1/2-\alpha}\mathbf{A})+\frac{d}{n}1$, so when we call Centering in Line 2, we again obtain $xs\approx_{\epsilon}\mu\cdot\tau(x,s)$.

Conclusion:

Before the algorithm ends, we have $xs\approx_{1/4}\mu\tau$. In Line 2, we reduce $\Phi_{b}:=\|\mathbf{A}^{\top}x-b\|_{(\mathbf{A}^{\top}\mathbf{X}\mathbf{S}^{-1}\mathbf{A})^{-1}}^{2}$ to some small enough $\Phi_{b}=O(\delta/n^{2})$. By Corollary 56 this does not move the vector $x$ too much, i.e. we still have $x\cdot s\approx_{1/2}\mu\cdot\tau(x,s)$. Hence, by the choice of the new $\mu$ in Line 2, Lemma 12 shows that we can output a point $\widehat{x}$ with the desired properties. ∎ Finally, we analyze the cost of the Algorithm 2.

Algorithm 2 takes $O((nd+d^{3})\log^{O(1)}n\log(n/\delta))$ time with high probability in $n$.

Number of iterations:

Cost of DInverseD_{\textsc{Inverse}}:

Cost of DLeverageD_{\textsc{Leverage}}:

By Lemma 11, the total movement of the projection matrix is

where $K$ is the number of steps and $w=\mathbf{S}^{-1-2\alpha}x^{1-2\alpha}$. Then, Theorem 9 shows that the total cost is $O(\frac{n}{d}K^{2}\cdot d)=O(nK^{2})=O(nd\log^{2}(1/\delta))$.

Cost of DMatVec(x)D_{\textsc{MatVec}}^{(x)} and DMatVec(s)D_{\textsc{MatVec}}^{(s)}:

Cost of maintaining 𝐀⊤​𝐗¯​h\mathbf{A}^{\top}\overline{\mathbf{X}}h (DGradientD_{\text{Gradient}}):

Cost of implementing MaintainFeasibility (Algorithm 10):

Removing the extra log⁡(1/δ)\log(1/\delta) term:

We note that all the extra $\log(1/\delta$) terms are due to running the data structures for $\sqrt{d}\log(1/\delta)$ steps. However, we can we reinitialize the data structures every $\sqrt{d}$ iterations. This decreases the $K$ dependence from $K^{2}$ to $K\sqrt{d}$.

Independence and Adaptive Adversaries:

Randomized data-structures often can not handle inputs that depend on outputs of the previous iteration. For example Theorem 10 is such a case, where the input to Update is not allowed to depend on the output of any previous calls to Update. In our Algorithm 2 the input to the data-structures inherently depends on their previous output, so here we want to verify that this does not cause any issues.

Robust Primal Dual LS-Path Following

Similar to previous papers [CLS19, LSZ19, Bra20] our algorithms measure progress or the quality of this approximation by

To update our points and improve the potential we consider attempting to move $w$ in some direction $h$ suggested by $\nabla\Phi(v)$. To do so we solve for $\delta_{s},\delta_{x},\delta_{y}$ satisfying the following

where $\overline{\mathbf{W}}=\overline{\mathbf{X}}\overline{\mathbf{S}}$. Solving this system of equations, we have

Note that, since we use $\mathbf{H}$ and $\mathbf{Q}$ (rather than a true orthogonal projection) it is not necessarily the case that $\mathbf{A}^{\top}\delta_{x}=0$. Consequently, we will need to take further steps to control this error and this is analyzed in Section B.

Much of the remaining analysis is bounding the effect of such a step and leveraging this analysis to tune $\phi$ and $h$. In Section 5.1 we bound the change in the point when we take a Newton step, in Section 5.2 we bound the effect of a Newton step on centrality for arbitrary $\Phi$, and in Section 5.3 we analyze the particular structure of $\Phi$ and use this to prove Theorem 1.

Next, we relate leverage scores of the projection matrix $\mathbf{P}(\mathbf{S}^{-1/2}\mathbf{X}^{1/2}\mathbf{A})$ considered in taking a Newton step to the leverage scores used to measure centrality, i.e. $\sigma(\mathbf{S}^{-1/2-\alpha}\mathbf{X}^{1/2-\alpha}\mathbf{A})$.

For $(\mu,\epsilon)$-centered point $(x,s)$ with $\epsilon\in[0,1/80)$ we have

By assumption $w=\mathbf{X}s$ satisfies $w\approx_{\epsilon}\mu[\sigma(\mathbf{S}^{-1/2-\alpha}\mathbf{X}^{1/2-\alpha}\mathbf{A})+\frac{d}{n}]$. Since leverage scores lie between $0$ and $1$, this implies $e^{-\epsilon}\frac{d}{n}\leq\mu^{-1}w\leq e^{\epsilon}2$. Consequently $e^{-\alpha\epsilon}2^{-\alpha}\mathbf{I}\preceq\mu^{\alpha}\mathbf{W}^{-\alpha}\preceq e^{\alpha\epsilon}(\frac{d}{n})^{-\alpha}\mathbf{I}$ and since $\alpha=1/(4\log(4n/d))$ we have that entrywise

The result then follows from the assumption $\epsilon\leq 1/80$. ∎

Leveraging Lemma 19 we obtain the following bounds on $\mathbf{Q}$ and $\overline{\mathbf{W}}^{-1/2}\mathbf{Q}\overline{\mathbf{W}}^{1/2}h$ for Newton directions.

Since $\mathbf{H}\approx_{\epsilon}\mathbf{A}^{\top}\overline{\mathbf{S}}^{-1}\overline{\mathbf{X}}\mathbf{A}$, we have $\mathbf{H}^{-1}\preceq e^{\epsilon}(\mathbf{A}^{\top}\overline{\mathbf{S}}^{-1}\overline{\mathbf{X}}\mathbf{A})^{-1}$ and

Since $\mathbf{Q}$ is PSD this implies $\|\mathbf{Q}\|_{2}\leq e^{\epsilon}$ and since $\overline{w}\approx_{2\epsilon}w$ and $w\approx_{\epsilon}\mu\tau$ this implies

Similarly, we have $0\preceq\mathbf{Q}\preceq e^{\epsilon}\mathbf{I}$ and hence $\|\mathbf{I}-\mathbf{Q}\|_{2}\leq 1$ (using that $\epsilon\leq 1/80$). Therefore, we have

Further, by Cauchy Schwarz and $\mathbf{Q}=\mathbf{Q}^{1/2}\mathbf{Q}^{1/2}$ we have

Since $\mathbf{Q}\preceq e^{\epsilon}\mathbf{I}$ we have

Further, by $\mathbf{Q}\preceq e^{\epsilon}\mathbf{P}(\overline{\mathbf{S}}^{-1/2}\overline{\mathbf{X}}^{1/2})\preceq e^{3\epsilon}\mathbf{P}(\mathbf{S}^{-1/2}\mathbf{X}^{1/2})$ and Lemma 19 we have that

Combining these and using $\epsilon\leq 1/80$ yields the desired bounds. ∎

Leveraging Lemma 19 we obtain the following bounds on the multiplicative stability of Newton directions.

Note that $\mathbf{S}^{-1}\delta_{s}=(1-2\alpha)\mathbf{S}^{-1}\overline{\mathbf{S}}\overline{\mathbf{W}}^{-1/2}\mathbf{Q}\overline{\mathbf{W}}^{1/2}h$. Since $s\approx_{\epsilon}\overline{s}$ and $\overline{w}\approx_{2\epsilon}w$ the bounds on $\|\mathbf{S}^{-1}\delta_{s}\|_{\tau}$ and $\|\mathbf{S}^{-1}\delta_{s}\|_{\infty}$ follow immediately from Lemma 18 and Lemma 20.

Similarly, $\mathbf{X}^{-1}\delta_{x}=(1+2\alpha)\mathbf{X}^{-1}\overline{\mathbf{X}}\overline{\mathbf{W}}^{-1/2}(\mathbf{I}-\mathbf{Q})\overline{\mathbf{W}}^{1/2}h$. Since $\overline{x}\approx_{\epsilon}x$ and $\overline{w}\approx_{2\epsilon}w$, we have $\|\mathbf{X}^{-1}\overline{\mathbf{X}}h\|_{\tau}\leq e^{\epsilon}\|h\|_{\tau}$ and $\|\mathbf{X}^{-1}\overline{\mathbf{X}}h\|_{\infty}\leq e^{\epsilon}\|h\|_{\infty}$. The bounds on $\|\mathbf{X}^{-1}\delta_{x}\|_{\tau}$ and $\|\mathbf{X}^{-1}\delta_{x}\|_{\infty}$ follow by triangle inequality and the same derivation as for $\|\mathbf{S}^{-1}\delta_{s}\|_{\tau}$ and $\|\mathbf{S}^{-1}\delta_{s}\|_{\infty}$. ∎

2 Newton Step Progress

Here we analyze the effect of a Newton step from a $(\mu,\epsilon)$ centered point on the potential $\Phi(x,s,\mu)$. We show that up to some additive error on $h$ in the appropriate mixed norm the potential decreases by $\nabla\Phi(v)^{\top}h$. The main result of this section is the following.

where $[\nabla\Phi(v^{\prime})]_{i}$ is the derivative of $\phi$ at $(v^{\prime})_{i}$ for all $i\in[n]$.

We prove this theorem in multiple steps. First, in Lemma 23 we directly compute the change in the potential function by chain rule and mean value theorem. Then, in Lemma 24 and Lemma 25 we bound the terms in this change of potential formula and use this to provide a formula for the approximate the change in the potential in Lemma 26. Finally, in Lemma 26 we bound this approximate change and use this to prove Theorem 22.

In the setting of Theorem 22 for all $t\in$ let $x_{t}=x+t\delta_{x}$ and $s_{t}=s+t\delta_{s}$. Then for some $t_{*}$, we have that

By the mean value theorem, there is $t\in$ such that

and Lemma 45 combined with chain rule implies that

The result then follows by the definition of $\delta_{s}$ and $\delta_{x}$. ∎

On the other hand, we know that $\boldsymbol{0}\preceq\boldsymbol{\Lambda}_{t}\preceq\boldsymbol{\Sigma}_{t}\preceq\mathbf{T}_{t}$ by Lemma 44 and that $\sigma_{t}\leq\tau_{t}$. Consequently, $\boldsymbol{0}\preceq\mathbf{T}_{t}^{-1/2}\boldsymbol{\Lambda}_{t}\mathbf{T}_{t}^{-1/2}\preceq\mathbf{I}$ and Lemma 18 yields

Since for all $c\in$ we have $\min\{\|y\|_{\infty},\|y\|_{\tau_{t}}\}\leq c\|y\|_{\infty}+(1-c)\|y\|_{\tau_{t}}$ and $\max\{1-\beta,\beta\}\leq 1-\beta$ as $\beta\leq 1/2$ combining (5.1) and (5.2) yields that for all $c\in$ it is the case that

Next, since by assumptions $w_{0}\approx_{\epsilon}\mu\cdot\tau_{0}$ and $\overline{w}\approx_{2\epsilon}w_{0}$ we have that $\overline{w}\approx_{6\epsilon}\mu\cdot\tau_{t}$ and therefore the above implies

In the setting of Lemma 23 there is a matrix $\mathbf{E}_{t}$ with $\|\mathbf{E}_{t}\|_{\tau_{t}+\infty}\leq 1000\epsilon$ such that

By Lemma 25, $x_{t}\approx_{\frac{3}{2}\epsilon}x_{0}$, $s_{t}\approx_{\frac{3}{2}\epsilon}s_{0}$, and $\tau_{t}\approx_{3\epsilon}\tau_{0}$. By definition of an $\epsilon$-Newton direction this further implies that $x_{t}\approx_{\frac{5}{2}\epsilon}\overline{x}$ and $s_{t}\approx_{\frac{5}{2}\epsilon}\overline{s}$. Further, this implies that $w_{t}\approx_{3\epsilon}w_{0}$ and since $w_{0}\approx_{\epsilon}\mu\cdot\tau_{0}$ by definition of $(\mu,\epsilon)$-centered this implies that $w_{t}\approx_{7\epsilon}\mu\cdot\tau_{t}$. Consequently, Fact 18, Lemma 46, and $\epsilon<1/80$ implies

Further, $\|\overline{\mathbf{W}}^{-1/2}\mathbf{Q}\overline{\mathbf{W}}^{1/2}\|_{\tau_{t}+\infty}\leq 2$ by Lemma 25 with $\beta=0$ and $\|\mathbf{T}_{t}^{-1}\boldsymbol{\Lambda}_{t}\|_{\tau_{t}+\infty}\leq 2$ by Lemma 24 with $\beta=0$. Combining these facts and applying Lemma 18 we have that there are $\mathbf{E}_{t}^{s}$ and $\mathbf{E}_{t}^{x}$ with $\|\mathbf{E}_{t}^{s}\|_{\tau_{t}+\infty}\leq 500\epsilon$ and $\|\mathbf{E}_{t}^{x}\|_{\tau_{t}+\infty}\leq 500\epsilon$ with

In the setting of Lemma 26 we have $\|\mathbf{J}_{t}-\mathbf{E}_{t}-\mathbf{I}\|_{\tau_{t}+\infty}\leq 1-\alpha$.

Lemma 26 implies that for $\beta=2\alpha(1-4\alpha^{2})^{-1}$

Further, Lemma 24 and $\alpha\in(0,1/5)$ implies that $\beta\in[0,1/2]$ and therefore

We now have everything we need to prove the main theorem.

where by Lemma 26 and Lemma 26 we know that for some matrix $\mathbf{E}_{t_{*}}$ with $\|\mathbf{E}_{t_{*}}\|_{\tau_{t_{*}}+\infty}\leq 1000\epsilon$ it holds that $\|\mathbf{J}_{t_{*}}-\mathbf{E}_{t_{*}}-\mathbf{I}\|_{\tau_{t_{*}}+\infty}\leq 1-\alpha$. Hence, we have $\mathbf{J}_{t_{*}}h=h+e$ where

Lemma 21 implies that $x_{t}\approx_{\frac{3}{2}\epsilon}x_{0}$, $s_{t}\approx_{\frac{3}{2}\epsilon}s_{0}$, and $\tau_{t}\approx_{3\epsilon}\tau_{0}$. Consequently, $w_{t}\approx_{3\epsilon}w$ and $v_{t_{*}}\approx_{6\epsilon}\mu\mathbf{W}^{-1}\tau$. The definition of $(\mu,\epsilon)$-centered implies that $\|\mu\mathbf{W}^{-1}\tau-1\|_{\infty}\leq\epsilon$ and therefore $\|v_{t_{*}}-\mu\mathbf{W}^{-1}\tau\|_{\infty}\leq 10\epsilon$. The bound on $\|e\|_{\tau_{t_{*}}+\infty}$ follows from (5.5) and $\|h\|_{\tau_{*}+\infty}\leq(1+5\epsilon)\|h\|_{\tau+\infty}$. ∎

3 Following the Central Path

Here we use the particular structure of $\Phi$ to analyze the effect of Newton steps and changing $\mu$ on centrality. First we state the following lemma about the potential function from [LS19]. Then we use it to analyze the effect of centering for one step, Lemma 31, and we conclude the section by proving a simplified variant of Theorem 5.

When $\tau$ is clear from context we also write $v^{\flat}$ for $v^{\flat(\tau)}$.

Here we provide a general lemma bounding how much $\Phi$ can increase for a step of bounded size in the mixed norm.

where in the second line we applied (5.8) of Lemma 28 as $t\delta\leq 1/(5\lambda)$. ∎

Further, since $\|\mathbf{X}^{-1}\delta_{x}\|_{\infty}\leq\epsilon\leq 1/(100)$ we have that

Now, the proof of Lemma 24 shows that $\|\left(1-2\alpha\right)\mathbf{T}_{t}^{-1}\boldsymbol{\Lambda}_{t}-\mathbf{I}\|_{\tau_{t}+\infty}\leq 2$ and that $\|\mathbf{X}^{-1}\delta_{x}\|_{\infty}\leq\epsilon\leq 1/100$ and $(x_{0},s)$ is $(\mu,\epsilon)$-centered shows that $\|\mu\mathbf{W}_{t}^{-1}\mathbf{T}_{t}\|_{\tau_{t}+\infty}\leq 2$ and $\|\mathbf{X}_{t}^{-1}\mathbf{X}_{0}\|_{\tau_{t}+\infty}\leq 2$. Combining yields the result. ∎

where we used (5.7), $\|\overline{v}-v^{\prime}\|_{\infty}\leq\|\overline{v}-v\|_{\infty}+\|v-v^{\prime}\|_{\infty}\leq 11\gamma$ and Theorem 22. Since $\gamma\leq\frac{\alpha}{50\lambda}$ and $\epsilon\leq\frac{\alpha}{4000}$, (5.8) implies that

Further $\|\tau_{1}\|_{1}=2d$ and $\|v_{2}\|_{\infty}\leq 2$ imply that

Using the bounds on $\delta_{\mu}$ and $\|\mathbf{X}^{-1}e_{x}\|_{\tau_{1}+\infty}\leq\gamma\alpha 2^{-15}$ this implies that $\|v_{3}-v_{1}\|_{\tau_{1}+\infty}\leq\gamma\alpha 2^{-10}$. Consequently, applying Lemma 29 yields

where in the second line we applied Lemma 28 and $\|v-v_{1}\|_{\infty}\leq 10\lambda$. Combining with (5.9) yields

Finally, for any vector $u$, we have \|u\|_{*}\geq\left\langle u,\frac{\alpha\text{\cdotsign}(u)}{40\sqrt{d}}\right\rangle=\frac{\alpha}{40\sqrt{d}}\|u\|_{1} as

and combining yields the desired result. ∎

Furthermore, during the algorithm 1, we have

$xs\approx_{4\epsilon}\mu\cdot\tau(x,s)$ for some $\mu$ where $(x,s)$ is immediate points in the algorithms

Initially, we have $\Phi(v)\leq 2n\cdot e^{3\epsilon\lambda}$ because of $x^{\textrm{(init)}}s^{\textrm{(init)}}\approx_{2\epsilon}\mu^{\textrm{(init)}}\cdot\tau(x^{\textrm{(init)}},s^{\textrm{(init)}})$ and (5.6) of Lemma 28.

We proceed to show that in each iteration $\Phi(v)\leq\min\{2n\cdot e^{3\epsilon\lambda},2^{16}n\sqrt{d}\alpha^{-2}\}$. Note that, by Lemma 28 in any iteration this holds, this implies that $\|v-1\|_{\infty}\leq\frac{1}{\lambda}\log\left(\Phi(v)\right)$. Since $\lambda=\frac{2}{\epsilon}\log(2^{16}n\sqrt{d}/\alpha^{2})\geq\frac{2}{\epsilon}\log(2n)$ this implies that $\|v-1\|_{\infty}\leq 3.5\epsilon$ which by the choice of $\epsilon$ implies that $x,s$ is $(\mu,4\epsilon)$ centered. Further, Lemma 31 and the design of the algorithm then imply that if $\Phi(v)\geq 2^{16}n\sqrt{d}\alpha^{-2}$ in this iteration, then

Now recall that $\alpha=1/(4\log(4n/d))\geq 1/(2^{8}n)$ and consequently as $d\leq n$ we have

Consequently, the reasoning in the preceding paragraph implies that in each step $\mu$ is moved closer to $\mu^{\textrm{(target)}}$ by a $1-\frac{\gamma\alpha}{2^{15}\sqrt{d}}=1-\Omega(\frac{\epsilon\alpha}{\log(n)\sqrt{d}})$ multiplicative factor. Hence, it takes $O(\frac{\sqrt{d}\log(n)}{\epsilon\alpha}\log(\frac{\mu^{\textrm{(target)}}}{\mu^{\textrm{(init)}}}))$ iterations to arrive $\mu^{\textrm{(target)}}$. Further, (5.10) shows that it takes $O(\sqrt{d}\log(n)/\alpha^{3})$ iterations to decrease $\Phi$ to $2^{16}n\sqrt{d}/\alpha^{2}$. Hence, in total, it takes

iterations for the algorithm to terminate. Further, when the algorithm terminates, we have that $x^{(\textrm{final})}s^{(\textrm{final})}\approx_{\epsilon}\mu^{\textrm{(target)}}\cdot\tau(x^{(\textrm{final})},s^{(\textrm{final})})$ because $\Phi\leq 2^{16}n\sqrt{d}\alpha^{-2}$, $\lambda=2\epsilon^{-1}\log(2^{16}n\sqrt{d}\alpha^{-2})$ and (5.6).

Finally, the bounds on the movement of $x$, $s$, follow from Lemma 21, the condition $\|\mathbf{X}^{-1}e_{x}\|_{\tau+\infty}\leq\frac{\gamma\alpha}{2^{20}}$, and the choice of $\gamma$. The movement of $\tau$ follows from Lemma 14 in [LS15]. ∎

Vector Maintenance

Here we show how to preprocess any matrix $\mathbf{A}$, such that we can build a data structure which supports quickly computing the large entries of the product $\mathbf{G}\mathbf{A}h$ for changing diagonal $\mathbf{G}$.

There exists a Monte-Carlo data-structure (Algorithm 4), that works against an adaptive adversary, with the following procedures:

for $\mathbf{G}=\mathbf{Diag}(g)$ in time $O(\|\mathbf{G}\mathbf{A}h\|^{2}\cdot\varepsilon^{-2}\cdot d\log^{3}n)$.

Scale($i,u$): Sets $g_{i}\leftarrow u$ in $O(d\log^{4}n)$ time.

Johnson-Lindenstrauss lemma is a famous result about low-distortion embeddings of points from high-dimensional into low-dimension Eculidean space. It was named after William B. Johnson and Joram Lindenstrauss [JL84]. The most classical matrix that gives the property is random Gaussian matrix, it is known that several other matrices also suffice to show the guarantees, e.g. FastJL[AC06], SparseJL[KN14], Count-Sketch matrix [CCFC02, TZ12], subsampled randomized Hadamard/Fourier transform [LDFU13].

It is worth to mention that JL lemma is tight up to some constant factor, i.e., there exists a set of points of size $n$ that needs dimension $\Omega(\epsilon^{-2}\log n)$ in order to preserve the distances between all pairs of points [LN17].

Our proof of Lemma 33 is through the analysis of Algorithm 4 which works as follows. The data-structure creates $\log_{2}n$ sketching matrices $\Phi_{1},...,\Phi_{\log_{2}n}$ according to Lemma 34, where each $\Phi_{j}$ is constructed for $\varepsilon=2^{-j}.$The data-structure then maintains the matrices $\mathbf{M}_{j}=\Phi_{j}\mathbf{G}\mathbf{A}$, whenever $\mathbf{G}$ is changed by calling Scale. Whenever $\textsc{Query}(h,\varepsilon)$ is called, the data-structure estimates $\|\mathbf{G}\mathbf{A}h\|_{2}$ via a Johnson-Lindenstrass matrix, and then pick the smallest $j$, for which $\mathbf{M}_{j}h=\Phi_{j}\mathbf{G}\mathbf{A}h$ is guaranteed to find all $i$ with $|(\mathbf{G}\mathbf{A}h)_{i}|>\varepsilon$, i.e. where $\Phi_{j}$ was initialized for $\varepsilon^{\prime}=2^{-j}\approx\|\mathbf{G}\mathbf{A}h\|_{2}$. To make sure that the algorithm works against an adaptive adversary, we compute these entries $(\mathbf{G}\mathbf{A}h)_{i}$ exactly and output only those entries whose absolute value is at least $\varepsilon$ (as we might also detect a few $i$ for which the entry is smaller).

We now give a formal proof, that this data-structure is correct.

Queries:

Now, If $r/\epsilon\geq\sqrt{n}$ then $\|\mathbf{G}\mathbf{A}h\|_{2}/\epsilon\geq\sqrt{n}/2$ and $O(nd)$ time is within the $O(\|\mathbf{G}\mathbf{A}h\|^{2}\cdot\varepsilon^{-2}\cdot d\log^{3}n)$ time-budget for the query. Hence, the algorithm simply outputs the vector by direct calculations, which takes $O(nd)$ time.

Otherwise, the algorithm uses Lemma 34 to find a list of indices $L$ which contains all $i$ such that $|(\mathbf{G}\mathbf{A}h)_{i}|\geq 2^{-j}\|\mathbf{G}\mathbf{A}h\|_{2}$ with probability at least $9/10$. Since $j=1+\left\lceil\log_{2}(r/\epsilon)\right\rceil$ and $r\leq 2\|\mathbf{G}\mathbf{A}h\|_{2}$, the heavy hitters algorithm outputs all $i$ such that $|(\mathbf{G}\mathbf{A}h)_{i}|\geq\epsilon$. Hence, the vector $v$ that we output is correct.

To bound the runtime, note that we find the list $L$ by first computing $y:=\mathbf{\Phi}_{j}\mathbf{G}\mathbf{A}h=\mathbf{M}_{j}h$ in $O(m_{j}d)=O((4^{j}\log^{2}n)d)$ time. Then, we use Lemma 34 to decode the list in time $O(4^{j}\log^{2}n)$ time. Hence, in total, it takes $O(\|\mathbf{G}\mathbf{A}h\|_{2}^{2}\varepsilon^{-2}\cdot d\log^{2}n)$.

Scaling:

When we set $g_{i}=u$, i.e by adding $u-g_{i}$ to entry $g_{i}$, then the product $\mathbf{\Phi}_{j}\mathbf{G}\mathbf{A}$ changes by adding the outer product $((u-g_{i})\cdot\mathbf{\Phi}_{j}1_{i})(1_{i}^{\top}\mathbf{G}\mathbf{A})$. This outer product can be computed in $O(d\log^{2}n)$ time, because each column of $\mathbf{\Phi}_{j}$ has only $O(\log^{2}n)$ non-zero entries. Since there are $O(\log n)$ many $j$, the total time is $O(d\log^{3}n)$. ∎

2 Maintaining Accumulated Sum

There exists a Monte-Carlo data-structure (Algorithm 5), that works against an adaptive adversary, with the following procedures:

Scale($i,u$): Sets $g_{i}=u$ in $O(d\log^{4}n)$ amortized time.

MarkExact($i$): Marks the $i$-th row to be calculated exactly in the next call to Query in $O(d\log^{4}n)$ amortized time.

$\textsc{ComputeExact}(i)$: Let $w=\sum_{k\in[t]}\mathbf{G}^{(k)}\mathbf{A}h^{(k)}$ be the vector approximated by the last call to Query.Then $\textsc{ComputeExact}(i)$ returns $w_{i}$ in $O(d)$ amortized time.

Finally, the cost of Query mainly consists of the following: (1) computation of sums in Line 5, (2) the calculation on indices in $F$ (i.e. calling ComputeExact), (3) the call of $D.\textsc{Query},$ and finally (4) the calls to $D.\textsc{Scale}$.

The cost of (1) is bounded by $O(td)$, which we charge as amortized cost to Update. The cost of (2) is bounded by $O(|F|d)$, which we charge to the amortized cost of Scale and MarkExact. The cost of (3) is exactly the runtime we claimed due to Lemma 33 since we zero out some of the rows using scale (Line 5) which causes the $(\mathbf{I}-\mathbf{F})$ term in the complexity.. The cost of (4) is bounded by $O(|F|)$ times the cost of $D.\textsc{Scale}$. We charge this cost to the Scale and MarkExact procedure of the current data structure.

We now prove that the output of Query is correct. Let $v$ be the output of the algorithm. The algorithm Query consists of two parts. From Line 5 to 5, we calculate $v$ on indices in $F$. From Line 5 to 5, we calculate $v$ on indices in $F^{c}$.

For any $i\notin F$, we note that $g_{i}$ is never changed. Hence, we have

and therefore we can re-write the $i$-th entry of $\sum_{k\in[t]}\left(\mathbf{G}^{(k)}-\mathbf{G}^{(0)}\right)\mathbf{A}h^{(k)}$ as

By the definition of $\bar{h}^{(k)}$ this is precisely what is computed in Line 5.

3 Vector Maintenance Algorithm and Analysis

We note that this assumption is satisfied for our interior point method. If this assumption is violated, we can detect it easily by Lemma 33 and propagate this change through the data-structure with small amortized cost. However, this would make the data structure have more special cases and make the code more difficult to read.

To analyze this data structure, recall that

and let $y^{(t)}$ be the output of the algorithm at iteration $t$. Our goal is to prove that the $x^{(t)}\approx_{\epsilon}y^{(t)}$. We prove this by induction.We assume $x^{(s)}\approx_{\epsilon}y^{(s)}$ for all $s<t$ and proceed to prove this is also true for $s=t$.

for all $i$. (Note that the sum on the left is using $t^{\prime}$ for its indices while the error bound on the right uses $x^{(t)}$.)

Correctness of xx:

ComputeExact:

Complexity:

The runtime of initialization and ComputeExact are that of Lemma 36, but increased by an $O(\log n)$ factor, as we run that many instances of these data-structures.

Adaptive Adversaries:

The algorithm works against adaptive adversaries, because the data structures of Lemma 36 work against adaptive adversaries. ∎

Leverage Score Maintenance

The IPM presented in this work requires approximate leverage scores of a matrix of the form $\mathbf{G}\mathbf{A}$, where $\mathbf{G}$ is a non-negative diagonal matrix that changes slowly from one iteration to the next. Here we provide a data-structure which exploits this stability to maintain the leverage scores of this matrix more efficiently than recomputing them from scratch every time $\mathbf{G}$ changes. Formally, we provide Algorithm 7 and show it proves the following theorem on maintaining leverage scores:

The high-level idea for obtaining this data-structure is that when $\mathbf{G}$ changes slowly, the leverage scores of $\mathbf{G}\mathbf{A}$ change slowly as well. Thus not too many leverage scores must be recomputed per iteration, when we are interested in some $(1+\varepsilon)$-approximation, as the previously computed values are still a valid approximation. The task of maintaining a valid approximation can thus be split into two parts: (i) for some small set $J$, quickly compute the $j$-th leverage score for $j\in J$, and (ii) detect which leverage scores should be recomputed, i.e. which leverage scores have changed enough so their previously computed value is no longer a valid approximation.

This section is structured according to the two tasks, (i) and (ii). In Section 7.1 we show how to solve (i), i.e. quickly compute a small number of leverage scores. In Section 7.2 we then show how to solve (ii) by showing how to detect which leverage scores have changed significantly and must be recomputed. The last Section 7.3 then combines these two results to prove Theorem 9.

Before we proceed to prove (i) and (ii), we first give a more accurate outline of how the algorithm of Theorem 9 works. An important observation is, that by standard dimension reduction tricks, a large part of Theorem 9 can be obtained by reduction to the vector data-structures of Section 6. Note that the $i$-th leverage score of $\mathbf{G}^{(t)}\mathbf{A}$ is

There exists a Monte-Carlo data-structure, that holds against an adaptive adversary, with the following procedures:

where $\mathbf{F}\in\{0,1\}^{n\times n}$ is a $0$-$1$ diagonal matrix with $\mathbf{F}_{ii}=1$ if and only if $i\in F$ in time

Scale($i,u$): Sets $g_{i}\leftarrow u$ in $O(d\log^{4}n)$ time.

The algorithm is directly obtained from Lemma 33. Let $\mathbf{H},F,\epsilon$ be the input to Query and let $D$ be an instance of Lemma 33. We perform $D.\textsc{Scale}(i,0)$ for all $i\in F$, so the data-structure $D$ represents $(\mathbf{I}-\mathbf{F})\mathbf{G}\mathbf{A}$. We then obtain $\mathbf{V}$ by setting the $i$-th column of $\mathbf{V}$ to the result of $D.\textsc{Query}(\mathbf{H}e_{i},\epsilon)$ for $i=1,...,b$. Afterward, we call $D.\textsc{Scale}$ again, to revert back to the previous scaling. The correctness and the time complexity of this new data-structure follow from Lemma 33. ∎

Now that we have Corollary 37, we can outline how we obtain Theorem 9 from it. Given a matrix $\Psi\approx(\mathbf{A}^{\top}(\mathbf{G}^{(t)})^{2}\mathbf{A})^{-1}$, it is easy to compute individual leverage scores (see Section 7.1, Lemma 38). However, such computation requires nearly linear work, which is quite slow for our purposes, and therefore we do not want to recompute every leverage score every time $\mathbf{G}^{(t)}$ changes a bit. Instead, we use Corollary 37 to detect which scores have changed a lot. This is done in Section 7.2, Lemma 39 and Lemma s40. The idea is that, if $\|1_{i}^{\top}\mathbf{G}^{(t)}\mathbf{A}(\mathbf{A}^{\top}(\mathbf{G}^{(t)})^{2}\mathbf{A})^{-1}\mathbf{A}^{\top}\mathbf{G}^{(t)}\|_{2}$ and $\|1_{i}^{\top}\mathbf{G}^{(k)}\mathbf{A}(\mathbf{A}^{\top}(\mathbf{G}^{(k)})^{2}\mathbf{A})^{-1}\mathbf{A}^{\top}\mathbf{G}^{(k)}\|_{2}$ differ substantially, for some $k<t$, then

should be quite large as well for some JL-matrix $\mathbf{R}$, i.e. something we can detect via Corollary 37.

We start our proof of Theorem 9 by showing how to quickly compute a small set of leverage scores. This corresponds to the function EstimateScore in Algorithm 7.

For the other case, at Line 7, we sample row $i$ with probability at least $\Theta(1)\frac{\log n}{\delta^{2}}\widetilde{\tau}_{i}^{(t)}\geq\Theta(1)\frac{\log n}{\delta^{2}}\sigma(\mathbf{G}^{(t)}\mathbf{A})_{i}$ (using the assumption on $\widetilde{\tau}_{i}^{(t)}$). Hence, the leverage score sampling guarantee shows that $\mathbf{A}^{\top}\widetilde{\mathbf{G}}^{2}\mathbf{A}\approx_{\delta/4}\mathbf{A}^{\top}(\mathbf{G}^{(t)})^{2}\mathbf{A}$ by choosing large enough constant in the $\Theta(1)$ factor in the probability. Hence, by the same argument as above, we have the result.

2 Detecting Leverage Scores Change

Now that we know how to quickly compute a small set of leverage scores, we must figure out which scores to compute. We only want to recompute the scores, that have changed substantially within some time-span. More accurately, we want to show that FindIndices (Algorithm 7) does indeed find all indices that changed a lot. The proof is split into two parts: We first show that we can detect large changes that happen within a sequence of updates of length $2^{i}$ in Lemma 39. We then extend the result to all large leverage score changes in Lemma 40.

Consider an execution of Line 7 (Algorithm 7) at iteration $t$. Suppose $\Psi^{(t)}\approx_{\epsilon/12\log n}(\mathbf{A}^{\top}(\mathbf{G}^{(t)})^{2}\mathbf{A})^{-1}$ and $\Psi^{(t-2^{i})}\approx_{\epsilon/36\log n}(\mathbf{A}^{\top}(\mathbf{G}^{(t-2^{i})})^{2}\mathbf{A})^{-1}$. Then $j\in J$ w.h.p. in $n$ if

Further, the set $J$ remains the same w.h.p. in $n$ if we replace Line 7 by $J_{i}\leftarrow[n]$.

First, we prove that for any $j\in[n]$ with

we have $j\in J_{i}$ (see Line 7). To do this, consider any $j\notin J_{i}$, i.e. $j$ is not in $F_{i}$ and not picked by the function $D.\textsc{Query}$. By Line 7 and Corollary 37, $j$ not picked by the function $D.\textsc{Query}$ implies that

where the second inequality follows from JL (Lemma 35)and the choice of $\mathbf{R}$. By the definition of $\approx$ and that $\epsilon_{\Psi}$ is small enough, we have

and similarly the fact that $\|e_{j}^{\top}\mathbf{G}^{(t-2^{i})}\mathbf{A}\Psi^{(t-2^{i})}\mathbf{A}\mathbf{G}^{(t-2^{i})})\|_{2}^{2}\approx_{\epsilon/(18\log n)}\sigma(\mathbf{G}^{(t-2^{i})}\mathbf{A})_{j}$, we have

Using $\epsilon_{\Psi}=\frac{\epsilon}{144\log n}\sqrt{\frac{d}{nb}}$, we showed that $j\notin J_{i}$ implies (7.2) is false. This proves the claim that for any $j$ such that (7.2) holds, we have that $j\in J_{i}$.

The previous lemma only showed how to detect leverage scores that changed within some time-span of length $2^{i}$, we now show that this is enough to detect all changes, no matter if they happen slowly (over some long time-span) or radically (within one iteration), or something in-between.

Under the assumption of Lemma 38 and Lemma 39 let $J^{(t)}$ be the output of FindIndices (Algorithm 7) at $t$-th iteration and

Then, we have that $J^{*(t)}\subset J^{(t)}$ w.h.p. in $n$.

for $z=0,1,\cdots,u-1$. Since there are only $2\log n$ steps, (7.3) shows that

because $\widetilde{\tau}$ has not been updated and that Lemma 38 shows

for the last update. Combining (7.4), (7.5) and (7.6), we see that $\widetilde{\tau}^{(t)}_{j}\approx_{\epsilon}\sigma(\mathbf{G}^{(t)}\mathbf{A})_{j}+\frac{d}{n}$ and hence $j\notin J^{*(t)}$. This shows that $j\notin J^{(t)}$ implies $j\notin J^{*(t)}$. Hence, $J^{*(t)}\subset J^{(t)}$. ∎

3 Leverage Score Maintenance

We can now combine the previous results Lemma 38, Lemma�39 and Lemma 40 to show that Algorithm 7 yields a proof of Theorem

We prove this statement by induction. At the start of the algorithm the vector $\widetilde{\tau}{}^{(0)}$ is a good enough approximation of the shifted leverage scores because of Lemma 38. We are left with proving that with each call to Query, $\widetilde{\tau}{}^{(t)}$ is updated appropriately.

By the induction, we have that $\widetilde{\tau}{}^{(t-1)}\approx_{\varepsilon}\sigma(\mathbf{G}^{(t-1)}\mathbf{A})+\frac{d}{n}$. Furthermore, we know that $g^{(t)}\approx_{1/16}g^{(t-1)}$. Hence, we have $\sigma(\mathbf{G}^{(t)}\mathbf{A})\approx_{1/4}\sigma(\mathbf{G}^{(t-1)}\mathbf{A})$. Using $\epsilon\in[0,1/4]$, we have that $\widetilde{\tau}{}^{(t-1)}\approx_{\varepsilon}\sigma(\mathbf{G}^{(t)}\mathbf{A})+\frac{d}{n}$. Hence, the $\widetilde{\tau}{}^{(t-1)}$ is a good enough approximation and satisfies the assumption of Lemma 38 and hence the assumption of Lemma 40. Lemma 40 shows that $\textsc{FindIndices}()$ finds all indices that is $\epsilon$ far away from the shifted leverage score. Lemma 38 then shows that those $\widetilde{\tau}$ will be correctly updated. Hence, all $\widetilde{\tau}$ is at most $\epsilon$ far away from the target. This finishes the induction. ∎

The correctness for the approximation ratio of Theorem 9 was already proven in Lemma 41. Here we bound the time complexity and prove that the data-structure works against an adaptive adversary.

To bound the complexity, we first bound how many leverage scores the data structure estimate, namely the size of $J$.

Size of all JJ:

where we used Johnson Lindenstrauss lemma 35 at the end. Now, summing up over $\log T$ scales and using that we restart every $O(n)$ iterations, we have

For the second term in (7.7), we bound it simply by $O(\log n)\cdot\sum_{t=1}^{T}\|g^{(t)}-g^{(t-1)}\|_{0}$. Hence, the sum of the size of all $J$ is bounded by