# An Improved Cutting Plane Method for Convex Optimization, Convex-Concave
Games and its Applications††thanks: A preliminary version of this paper
appears in Proceedings of the 52nd ACM Symposium on Theory of Computing (STOC
2020). The authors would like to express their sincere gratitude to matrix
multiplicationer Josh Alman and Jeroen Zuiddam for their patient and answers
of our exponential number of questions about fast matrix multiplication.

Haotian Jiang Yin Tat Lee Zhao Song Sam Chiu-wai Wong<EMAIL_ADDRESS>University of<EMAIL_ADDRESS>University of Washington and Microsoft
Research Redmond. Research supported in part by NSF Awards CCF-1740551,
CCF-1749609, and<EMAIL_ADDRESS>Princeton University and Institute
for Advanced Study. Work done while visiting University of
<EMAIL_ADDRESS>Microsoft Research Redmond.

Given a separation oracle for a convex set $K\subset\mathbb{R}^{n}$ that is
contained in a box of radius $R$, the goal is to either compute a point in $K$
or prove that $K$ does not contain a ball of radius $\epsilon$. We propose a
new cutting plane algorithm that uses an optimal $O(n\log(\kappa))$
evaluations of the oracle and an additional $O(n^{2})$ time per evaluation,
where $\kappa=nR/\epsilon$.

  * •

This improves upon Vaidya’s $O(\text{SO}\cdot
n\log(\kappa)+n^{\omega+1}\log(\kappa))$ time algorithm [Vaidya, FOCS 1989a]
in terms of polynomial dependence on $n$, where $\omega<2.373$ is the exponent
of matrix multiplication and SO is the time for oracle evaluation.

  * •

This improves upon Lee-Sidford-Wong’s $O(\text{SO}\cdot
n\log(\kappa)+n^{3}\log^{O(1)}(\kappa))$ time algorithm [Lee, Sidford and
Wong, FOCS 2015] in terms of dependence on $\kappa$.

For many important applications in economics, $\kappa=\Omega(\exp(n))$ and
this leads to a significant difference between $\log(\kappa)$ and
$\operatorname{poly}(\log(\kappa))$. We also provide evidence that the $n^{2}$
time per evaluation cannot be improved and thus our running time is optimal.

A bottleneck of previous cutting plane methods is to compute leverage scores,
a measure of the relative importance of past constraints. Our result is
achieved by a novel multi-layered data structure for leverage score
maintenance, which is a sophisticated combination of diverse techniques such
as random projection, batched low-rank update, inverse maintenance, polynomial
interpolation, and fast rectangular matrix multiplication. Interestingly, our
method requires a combination of different fast rectangular matrix
multiplication algorithms.

Our algorithm not only works for the classical convex optimization setting,
but also generalizes to convex-concave games. We apply our algorithm to
improve the runtimes of many interesting problems, e.g., Linear Arrow-Debreu
Markets, Fisher Markets, and Walrasian equilibrium.

###### Contents

  1. 1 Introduction
    1. 1.1 Applications
    2. 1.2 Previous works
    3. 1.3 Lee-Sidford-Wong method (LSW)
    4. 1.4 Overview of our approach
    5. 1.5 Discussion of optimality
    6. 1.6 Related works
  2. 2 Our Techniques
    1. 2.1 Efficient approximation of changes in leverage scores
    2. 2.2 Layered data structure
    3. 2.3 Batched low-rank update
    4. 2.4 Illustration of our analysis
    5. 2.5 Much faster rectangular matrix multiplication implies deterministic cutting plane method
  3. Acknowledgments
  4. 3 Preliminaries
    1. 3.1 Notations
    2. 3.2 Operators
    3. 3.3 Different types of running time
    4. 3.4 Basic results on matrices
    5. 3.5 Fast matrix multiplication
    6. 3.6 Multiple variable polynomial interpolation
  5. 4 Perturbed Volumetric Center Cutting Plane Method
  6. 5 Main Data Structure for Leverage Score Maintenance
  7. 6 Batched Low Rank Update
    1. 6.1 Simple low rank update
    2. 6.2 Batched low rank update
  8. 7 Simple Deterministic Leverage Score Maintenance
    1. 7.1 Main result
    2. 7.2 Approximate leverage score’s moving
    3. 7.3 $\ell_{2}$-error bound
  9. 8 Complicated Randomized Leverage Score Maintenance
    1. 8.1 Main result
    2. 8.2 Leverage score’s moving
    3. 8.3 Running time anlaysis
    4. 8.4 Upper bounding $\eta$, $\alpha$, $\beta$ and $\gamma$
      1. 8.4.1 Upper bounding $\sum_{i=1}^{m}\|\eta_{i}\|_{2}^{4}$
      2. 8.4.2 Upper bounding $\sum_{i=1}^{m}\|\alpha_{i}\|_{2}^{2}\|\beta_{i}\|_{2}^{2}$
      3. 8.4.3 Upper bounding $\sum_{i=1}^{m}\|\gamma_{i}\|_{2}^{4}$
    5. 8.5 Variance upper bound for random Gaussian matrices
      1. 8.5.1 A variance bound for random Gaussian matrices
      2. 8.5.2 Variance bound for $\eta^{\top}\eta$, $\alpha^{\top}\beta$ and $\gamma^{\top}\gamma$
    6. 8.6 Error upper bound for discrete sampling
  10. A Perturbed Volumetric Center Cutting Plane Method
    1. A.1 Our faster implementation via leverage score maintenance
  11. B Modified Projection Maintenance
  12. C Cutting Plane Method for Convex Minimization and Saddle Point Problems
    1. C.1 Notations and definitions
    2. C.2 Convex minimization
    3. C.3 Convex-concave games
    4. C.4 Applying cutting plane method to convex-concave games
    5. C.5 Cutting plane method for convex-concave games: generating solutions
  13. D Applications of Cutting Plane Method
    1. D.1 Linear Arrow-Debreu markets
    2. D.2 Fisher markets with spending constraint utilities
    3. D.3 Walrasian equilibrium for general buyer valuations and fixed supply

##  1 Introduction

The cutting plane methods are a class of optimization primitives for solving
Linear and Convex Programming, which is encapsulated by the _feasibility_
problem of finding a point in a given convex set $K$ equipped with a
separation oracle. In each iteration, cutting plane methods query the
_separation_ oracle which returns a hyperplane separating the query point from
$K$. Since Khachiyan’s breakthrough result [Kha80] on the _ellipsoid method_ ,
cutting plane methods have played a central role in theoretical computer
science, showing that a plethora of problems from diverse areas admit
polynomial time algorithms. Early prominent examples include linear
programming and submodular function minimization.

While more applications of the cutting plane methods are still being
discovered to this date, progress on faster cutting plane methods had
stagnated until the recent work of Lee, Sidford and Wong [LSW15]. Prior to
their work, cutting plane methods were considered too slow to be of relevance
in theory or in practice (other than establishing polynomiality). [LSW15]
debunked this by showing how their cutting plane method is faster than tailor-
made algorithms for a long list of well-studied problems.

The fastest algorithm before their work was Vaidya’s algorithm [Vai89a], which
runs in $O(n^{\omega+1})$ time. Here $\omega<2.373$ is the exponent of matrix
multiplication [CW87, Wil12, DS13, LG14]. Leveraging recent advances in
optimization and numerical linear algebra, LSW lowered the dependence on the
dimension $n$ at the expense of additional log factors in the accuracy
$1/\kappa$, namely $\log^{2}(\kappa)$. The following table compares the
runtimes of Vaidya’s and LSW methods, which both achieve the optimal number of
oracle calls of $n\log(\kappa)$.

Authors | #Oracle Calls | Runtimes
---|---|---
Vaidya | $n\log(\kappa)$ | $n^{\omega+1}\log(\kappa)$
Lee-Sidford-Wong | $n\log(\kappa)$ | $n^{3}\log^{O(1)}(n)\log^{3}(\kappa)$

This extra overhead of $\log^{2}(\kappa)$, as exemplified by various problems
in combinatorial optimization in their paper, translates into only a log-
squared factor in the maximum value $M$ of the input by taking
$\epsilon=O(1/M)$. This is because solutions to (polynomial-time solvable)
problems in combinatorial optimization are guaranteed to be integral. Despite
being a nuisance, this small overhead is relatively mild and can even be
absorbed completely for unweighted problems and strongly polynomial time
algorithms. Indeed, armed with this faster cutting plane method, they improved
the state-of-the-art running times for a host of problems such as semidefinite
programming, matroid intersection and submodular minimization.

##### Importance of $\log(\kappa)$ vs $\log^{3}(\kappa)$.

For many other problems including linear programming and market equilibrium
computation, $\epsilon$ must be taken as $1/M^{O(n)}$ resulting in an
additional factor $n^{2}$ between $\log(\kappa)$ and $\log^{3}(\kappa)$. For
these applications, LSW is slower than Vaidy’s by a factor of
$\widetilde{O}(n^{4-\omega})$. This raises a natural question:_is there a
cutting plane method that simultaneously runs in $O(n\log(\kappa))$ calls and
$O(n^{3}\log^{O(1)}(n)\log(\kappa))$ time? That is, can we achieve the best of
both worlds of Vaidya’s and LSW methods in terms of the dependence on $n$ and
$\kappa$?_

In this paper, we answer this question in the affirmative. Somewhat
surprisingly, we are able to remove $\log^{O(1)}(n)$ dependence in LSW as
well.

######  Theorem 1.1 (Main result).

There is a cutting plane method which runs in time
$O(n\cdot\operatorname{SO}\log(\kappa)+n^{3}\log(\kappa))$, where
$\operatorname{SO}$ is the time complexity of the separation oracle.

Reference | Year | Algorithm | Complexity
---|---|---|---
[Sho77, YN76, Kha80] | 1979 | Ellipsoid Method | $n^{2}\operatorname{SO}\log(\kappa)+n^{4}\log(\kappa)$
[KTE88, NN89] | 1988 | Inscribed Ellipsoid | $n\operatorname{SO}\log(\kappa)+(n\log(\kappa))^{4.5}$
[Vai89a] | 1989 | Volumetric Center | $n\operatorname{SO}\log(\kappa)+n^{\omega+1}\log(\kappa)$
[AV95] | 1995 | Analytic Center | $n\operatorname{SO}\log^{2}(\kappa)+n^{\omega+1}\log^{2}(\kappa)+(n\log(\kappa))^{2+\omega/2}$
[BV02] | 2004 | Random Walk | $n\operatorname{SO}\log(\kappa)+n^{7}\log(\kappa)$
[LSW15] | 2015 | Hybrid Center | $n\operatorname{SO}\log(\kappa)+n^{3}\log^{3}(\kappa)$
Theorem 1.1 | 2019 | Volumetric Center | $n\operatorname{SO}\log(\kappa)+n^{3}\log(\kappa)$
Table 1:  Algorithms for the Feasibility Problem. Let $\kappa=nR/\epsilon$. All methods can be used to solves a more general problem where only a membership oracle is given [LSV18]. Reference | Year | #Operations | Poly Type
---|---|---|---
[Eav75] | 1975 | Finite | Not poly
[Jai07] | 2007 | Polynomial | Weakly poly
[Ye08] | 2008 | $n^{6}\log(nU)$ | Weakly poly
[DPSV08] | 2008 | Polynomial | Weakly poly
[DM15] | 2015 | $n^{9}\log(nU)$ | Weakly poly
[DGM16] | 2016 | $n^{6}\log^{2}(nU)$ | Weakly poly
[GV19] | 2019 | $mn^{9}\log^{2}(n)$ | Strongly poly
Theorem 1.4 | 2019 | $mn^{2}\log(nU)$ | Weakly poly
Table 2: Linear Arrow-Debreu Markets. Let $n$ be the number of agents and $m$ the number of edges. Each operation of the given algorithms involves $O(n\log(nU))$-bit numbers.  Reference | Year | # Operations | Poly Type
---|---|---|---
[Vaz10] | 2010 | $n^{3}(n+m)^{2}\log(U)\cdot\mathcal{T}_{\text{max-flow}}$ | Weakly poly
[Vég16] | 2016 | $mn^{3}+m^{2}(m+n\log(n))\log(m)$ | Strongly poly
[Wan16] | 2016 | $m^{3}n+m^{2}\log(n)(n\log(n)+m)$ | Strongly poly
Theorem 1.5 | 2019 | $mn^{2}\log(nU)$ | Weakly poly
Table 3: Fisher Markets with Spending Constraint Utilities. Let $n$ denote the total number of buyers and sellers, and $m$ the total number of segments. $\mathcal{T}_{\text{max-flow}}$ denotes the number of operations needed for a max-flow computation. Each operation of the given algorithms involves $O(n\log(nU))$-bit numbers.  Reference | Year | # Operations | Poly Type
---|---|---|---
[KJC82] | 1982 | Finite | Not poly
[Par99] | 1999 | Finite | Pseudo poly
[PU02] | 2002 | Finite | Pseudo poly
[AM02] | 2002 | Finite | Pseudo poly
[dVSV07] | 2007 | Finite | Pseudo poly
[LW17] | 2017 | $n^{2}\mathcal{T}_{\mathrm{AD}}\log(SMn)+n^{6}\log^{O(1)}(SMn)$ | Weak poly
Theorem 1.6 | 2019 | $n^{2}\mathcal{T}_{\mathrm{AD}}\log(SMn)+n^{4}\log(SMn)$ | Weak poly
Table 4: Walrasian equilibrium for general buyer valuations and fixed supply.
Let $\mathcal{T}_{\mathrm{AD}}$ denotes the runtime of aggregate demand
oracle, $n$ denotes the number of goods.

As with previous methods, our result achieves the asymptotic optimal oracle
complexity of $O(n\cdot\operatorname{SO}\log(\kappa))$. Moreover, we
conjecture that the runtime is likely the best possible. Note that in each
iteration where a separation oracle call is made, even basic matrix operations
require $O(n^{2})$ time. Thus our runtime is essentially tight unless
properties like sparsity can be surprisingly exploited to update the feasible
region and compute the next query point for the separation oracle.

Our result is obtained by marrying Vaidya’s method with recent advances in
numerical linear algebra. Similar to LSW, we implement each iteration of
Vaidya’s via tools from fast numerical linear algebra. The main issue is that
such tools rely on approximation and would lead to errors which must be
carefully controlled.

Our first innovation is to run Vaidya’s method in _phases_ , each of which
consist of a certain number of iterations. Between phases we “recompute” to
eliminate the errors accumulated within a phase. Because of this recomputation
we can afford to tolerate higher errors in an iteration. Secondly, we present
a sophisticated data structure that enables us to implement each iteration of
Vaidya’s efficiently. Our data structure leverages recent advances on applying
numerical techniques to optimization. We hope that these numerical tools, as
well as our approach to applying them, would play a greater role in future
development of optimization.

###  1.1 Applications

We highlight some of the key applications of our faster cutting plane method.
Interestingly, even though our cutting plane method is a general purpose
algorithm, we are able to improve the runtimes of tailor-made algorithms for
various problems.

Using a standard reduction of convex minimization to the feasiblity problem
([Nem94] and Theorem 42 of [LSW15]), we can minimize a convex function with an
optimal $\widetilde{O}(n)$ subgradient oracle calls and an additional
$O(n^{2})$ time per oracle call.

######  Theorem 1.2 (Informal version of Theorem C.1).

Let $f$ be a convex function on $\mathbb{R}^{n}$ and $S$ be a convex set that
contains a minimizer of $f$. Suppose we have a subgradient oracle for $f$ with
cost $\mathcal{T}$ and $S\subset B(0,R)$. Using $B(0,R)$ as the initial
polytope for our Cutting Plane Method, for any $0<\alpha<1$, we can compute
$x\in S$ such that $f(x)-\min_{y\in S}f(y)\leq\alpha\left(\max_{y\in
S}f(y)-\min_{y\in S}f(y)\right)$, with high probability in $n$ and with a
running time of $O(\mathcal{T}\cdot n\log(\kappa)+n^{3}\log(\kappa))$, where
$\kappa=n\gamma/\alpha$ and $\gamma=R/(\mathrm{minwidth}(S))$.

Our convex minimization result can be further generalized to convex-concave
games.

######  Theorem 1.3 (Informal version of Theorem C.9).

Given convex sets $\mathcal{X}\subset B(0,R)\subset\mathbb{R}^{n}$ and
$\mathcal{Y}\subset B(0,R)\subset\mathbb{R}^{m}$ such that both $\mathcal{X}$
and $\mathcal{Y}$ contain a ball of radius $r$. Let
$f(x,y):\mathcal{X}\times\mathcal{Y}\rightarrow\mathbb{R}$ be an $L$-Lipschitz
function that is convex in $x$ and concave in $y$. Define
$\kappa=\frac{n+m}{\epsilon}\frac{R}{r}$. For any $\epsilon\in(0,1/2]$, we can
find $(\widehat{x},\widehat{y})$ such that
$\max_{y\in\mathcal{Y}}f(\widehat{x},y)-\min_{x\in\mathcal{X}}f(x,\widehat{y})\leq\epsilon
Lr$ in time
$O\left(\mathcal{T}\cdot(n+m)\log(\kappa)+(n+m)^{3}\log(\kappa)\right)$ with
high probability in $n+m$ where $\mathcal{T}$ is the cost of computing
subgradient $\nabla f$.

Leveraging this improved dependence on $\kappa$ (and hence $\epsilon$), our
cutting plane method can be used to improve the runtimes of a wide range of
problems, especially those on market equilibrium computation (see Table 2, 3
and 4 for a summary of previous runtimes). In all our applications,
$1/\epsilon$ needs to be exponentially large in $n$ which renders the
$\log(\kappa)$ factor polynomially large.

We show the following runtime improvement for the problem of computing a
market equilibrium in linear exchange markets:

######  Theorem 1.4 (Informal version of Theorem D.5).

There exists a weakly polynomial algorithm that computes a market equilibrium
in linear exchange markets in time $O(mn^{2}\log(nU))$.

The celebrated result of Arrow and Debreu [AD54] shows the existence of a
market equilibrium for a broad class of utility functions. Since then
researchers have attempted to design efficient algorithms to compute market
equilibria. One prominent special case of linear utilities has enjoyed
significant attention as demonstrated by the long line of work in Table 2.
Essentially, this problem can be captured by a convex program with linear
constraints [DGV16] (see (19)). While this convex program exhibits certain
advantageous features over previous ones [Cor89, Jai07, NP83], it had not led
to improved runtimes since the objective is not separable [GV19]. Moreover,
the number of variables in the convex program can be as large as $O(n^{2})$,
which prohibits a fast runtime for the cutting plane method if applied
directly. Our approach in Theorem 1.4 is to transform the convex program into
a convex-concave game with $O(n)$ variables and apply Theorem 1.3.

A similar technique can be applied to the problem of computing market
equilibrium for Fisher markets with spending utility constraints where a
convex program is given in [BDX10] (see (21)). However, directly transforming
it into a convex-concave game does not reduce the dimension of the variables,
as both the number of variables and constraints in the original convex program
is $\Theta(m)$, where $m$ is the total number of segments and it can be much
larger than $n^{2}$. In order to reduce the dimension of the convex-concave
game to $O(n)$, we express part of the variables as functions of $O(n)$
variables and show that this does not increase the time of the first-order
oracle. This leads to the following runtime improvement:

######  Theorem 1.5 (Informal version of Theorem D.7).

There exists a weakly polynomial algorithm that computes a market equilibrium
in Fisher markets with spending constraint utilities in time
$O(mn^{2}\log(nU))$.

Yet another economics application of our cutting plane method is the problem
of computing a Walrasian equilibrium in a market with fixed supply. In this
economy, buyers may have arbitrary valuation functions and we would like to
compute prices so that the market clears, i.e. the aggregate demand of the
buyers matches the fixed supply. Recently Paes Leme and Wong [LW17] gave a
polynomial time algorithm for this problem _provided_ that an equilbrium
actually exists. They achieved this by showing that the cutting plane method
can be modified so that the convex program for the equilibrium can be solved
under the aggregate demand oracle.

By leveraging our faster cutting plane method we obtain an improved runtime:

######  Theorem 1.6 (Informal version of Theorem D.9).

There is an algorithm that runs in time
$O(n^{2}\mathcal{T}_{\mathrm{AD}}\log(SMn)+n^{4}\log(SMn))$ for computing a
market equilibrium in an economy with general buyer valuation in the aggregate
demand model.

###  1.2 Previous works

The cutting plane methods solve the following _feasbility_ problem that
conveniently abstracts the applications to specific scenarios.

  * Feasibility Problem: Given a separation oracle for a set $K$ contained in a box of radius $R$ either find a point $x\in K$ or prove that $K$ does not contain a ball of radius $\epsilon$.

All cutting plane methods maintain a candidate region $\Omega$ and solve the
feasibility problem by iteratively refining $\Omega$ based on the present
$\Omega$ and the new separating hyperplane. In each iteration:
1\. The separation oracle is queried at some point $x\in\Omega$.
2\. If $x\in K$ we have solved the feasibility problem.
3\. Otherwise, the separation oracle returns a separating hyperplane from
which $\Omega$ is further refined and the next query point $x$ is computed.

Previous works differ in how $x$ is selected and how $\Omega$ is refined. For
instance, the classic ellipsoid method maintains $\Omega$ as an ellipsoid and
$x$ as its center. Given $\Omega$ and the new separating hyperplane, the new
$\Omega$ is chosen to be the smallest ellipsoid containing their intersection.
Table 1 lists the running times for solving the feasibility problem in the
literature.

We focus our discussion on the trade-off between the oracle complexity and the
runtime per iteration. The oracle complexity has a lower bound
$\Omega(n\log(\kappa))$ [NY83]. While the ellipsoid method achieves a
suboptimal $O(n^{2}\log(\kappa))$ in oracle complexity, the runtime per
iteration is $O(n^{2})$ which is faster than all subsequent methods. The good
runtime follows from the simple calculations needed to update the ellipsoid,
whereas the suboptimal oracle complexity can be attributed to the “looseness”
of maintaining only an ellipsoid as a proxy to the intersection of past
separating half-spaces.

Indeed, one can attain the optimal oracle complexity $O(n^{2}\log(\kappa))$ by
maintaining _all_ previous separating half-spaces. This is the random walk
method [BV02] where the query point $x$ is chosen to be its (approximate)
center of gravity. Updating $x$ involves performing a random walk in this
polytope and is computationally expensive.

Other cutting methods improve oracle complexity by maintaining more fine-
grained information about past separating hyperplanes in a way that is
computationlly friendly. Of particular relevance is Vaidya’s volumetric center
method [Vai89a]. Vaidy’s $\Omega$, similar to the random walk method, is also
a polytope $\Omega=\\{x\in\mathbb{R}^{n}\,:\,Ax\geq b\\}$. As the name
suggests, the query point is chosen to be the _volumetric center_ , which is
the minimizer of the following convex function defined by the feasible region
$\Omega$:

$\displaystyle\frac{1}{2}\log\det\left(A^{\top}S_{x}^{-2}A\right)\enspace\text{
where }\enspace S_{x}:=\text{diag}(Ax-b)\text{ is the diagonal matrix of the
slacks}$

Nevertheless, by judiciously including only a representative subset of
previous half-spaces $ax\geq b$, Vaidya showed that the volumetric center and
$\Omega$ can be updated by basic matrix operations which run in
$O(n^{\omega})$ time.

We defer the discussion of LSW to the next subsection, which explains how
Vaidya’s method can be sped up using machineries from numerical linear
algebra.

###  1.3 Lee-Sidford-Wong method (LSW)

LSW’s key observation is that Vaidya’s method relies heavily on _leverage
scores_ , which measure the relative importance of each separating hyperplane.
In Vaidya’s method, naively updating leverage scores requires $O(n^{\omega})$
time and is a bottleneck. Inspired by the work of Spielman and Srivistava
[SS11], LSW attempted to address this by using random Johnson-Lindenstrauss
[JL84] (JL) projection to approximate _changes_ in leverage scores. Leverage
scores can then be updated by summing over the differences.

Approximating leverage scores changes via JL projection however still requires
solving a linear system which would still take $O(n^{\omega})$ time. LSW
further overcame this barrier by resorting to a recent work that efficiently
solves “slowly-changing” linear system in amortized $\widetilde{O}(n^{2})$
time. Thus after paying $O(n^{\omega})$ initially, they can solve such linear
systems in $\widetilde{O}(n^{2})$ time per iteration.

Error accumulation. Nevertheless, as an approximate method JL introduces
errors which would accumulate across iterations. While Vaidya’s method
tolerates small errors in leverage scores, the total errors incurred in JL
projection (as accumulated across iterations) would eventually become too big
and destroy the performance guarantee of Vaidya’s method.

LSW handled this by modifying Vaidya’s framework to take into account of the
error in the convex function to be minimized. This approach gives rise to a
“hybrid” center algorithm, which involves a complicated interplay of
optimization and linear algebra. In particular, to reduce the error
accumulated the error parameter $\epsilon_{0}$ in JL projection has to be as
small as $\epsilon_{0}=1/\log(\kappa)$. As the runtime of JL depends on
$1/\epsilon_{0}^{2}$, this unfortunately introduces the $\log^{2}\kappa$
overhead in LSW runtime when compared to Vaidya’s volumetric center method.

###  1.4 Overview of our approach

To achieve the desired $O(n\operatorname{SO}\log(\kappa)+n^{3}\log(\kappa))$
runtime, we build a data structure that approximates changes in leverage
scores to within $c$ in $\ell_{2}$ norm for small enough constant $c$. This
would suffice for Vaidya’s cutting plane method. The desired runtime then
directly follows from the performance guarantee of the data structure, which
handles each iteration in amortized $O(n^{2})$ time. Designing such a data
structure requires several new ideas to control the error accumulation issue.

Our data structure employs a layered approach where different layers are
associated with different error tolerances, and achieve different accuracy-
efficiency tradeoffs in approximating the changes in leverage scores. The more
inner a layer is, the more error it can tolerate and the faster is the
runtime. Whenever the error accumulated in a layer becomes too much, the layer
above would take over and produce a finer error estimate which would, of
course, be more time costly. But because of our layered approach, the higher
layer is called on less often and afford to spend more time.

Such a layered approach further leads to the following issue. In the middle
and outer layers, we batch the updates of multiple steps into one which allows
us to make use of fast rectangular matrix multiplication. However, our
algorithm needs to handle possibly exponential weight changes. We show there
are not too many such weights and can be handled separately in groups of
$\operatorname{poly}\log(n)$ size using low-rank update formula.

Our data structure also draws on various numerical tools such as fast
rectangular matrix multiplication, “tall” JL projection, preconditioning,
inverse maintenance, and polynomial interpolation for approximating integrals.

A more in-depth discussion of our techniques can be found in Section 2.

###  1.5 Discussion of optimality

Similar to previous methods our cutting plane method achieves the optimal
oracle complexity $n\log(\kappa)$ [NY83]. We present some evidence that our
running time of $O(n^{3}\log(\kappa))$ is also tight.

A bottleneck of Vaidya’s method is to solve the inverse maintenance problem.
Formally, given a sequence of positive vectors $w^{1},w^{2},\cdots w^{T}$, let
$P(w)$ be defined as

$P(w)=\sqrt{W}A(A^{\top}WA)^{-1}A^{\top}\sqrt{W},$

where $W$ is the diagonal matrix such that $W_{i,i}=w_{i}$. The goal is to
output a sequence of vectors $v^{1},v^{2},\cdots,v^{T}$ such that

$\displaystyle v^{t}\approx w^{t}\text{\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ and\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ }P(v^{t})\approx P(w^{t}),\forall
t\in[T].$

There is a long line of research on inverse maintenance and dynamic matrix
data-structure problems [Kha80, Vai89b, San04, LS15, HKNS15, CLS19, LSZ19,
Son19, BNS19, Bra20]. This task can be done naively by spending $n^{\omega}$
time so the goal is to achieve $o(n^{\omega})$ amortized cost per iteration.
For example, in the LP setting the number of iterations is $O(\sqrt{n})$ and
Vaidya [Vai89b] combined fast matrix multiplication with inverse maintenance
to achieve $O(n^{2})$ amortized cost per iteration, which gives an
$O(n^{2.5})$ time algorithm. This remained a barrier until recent works
[CLS19, LSZ19] combined sampling and sketching techniques with fast
rectangular matrix multiplication and inverse matrix maintenance to give an
$O(n^{\omega})$ time algorithm.

One of the major computation required in each step is matrix-vector
multiplication, e.g., $P(w)\cdot h$. Naively, this step takes $O(n^{2})$ time
per iteration. To achieve $o(n^{2})$ amortized cost per iteration, previous
works [CLS19, LSZ19] used an idea called “iterating and sketching” which was
formally described in [Son19]. This idea is very different from the classical
“sketch and solve” [CW13] and “guess a sketch” [RSW16]. The classical idea
usually applying sketching matrices only once without modifying the solver
itself. However, the “iterating and sketching” idea has to modify the solver
and applying sketching/sampling matrices over each iteration.

In [CLS19], they use the diagonal sampling matrix $D\in\mathbb{R}^{n\times n}$
which has roughly $\sqrt{n}$ nonzeros. They use that sampling matrix to sample
on the right hand side:

$\displaystyle\sqrt{W}A(A^{\top}WA)^{-1}A^{\top}\sqrt{W}\underbrace{D}_{\text{sample\leavevmode\nobreak\
right}}h.$

In [LSZ19], they use the subsampled randomized Hadamard/Fourier matrix
[LDFU13, PSW17] $R\in\mathbb{R}^{\sqrt{n}\times n}$. They use the sketching
matrix to sketch on the left hand side:

$\displaystyle\underbrace{R^{\top}R}_{\text{sketch\leavevmode\nobreak\
left}}\sqrt{W}A(A^{\top}WA)^{-1}A^{\top}\sqrt{W}h.$

Compared to the cutting plane method, LP is an easier maintenance task as the
matrix $A$ is fixed throughout. In the cutting plane method, however, rows get
inserted into or deleted from $A$ from continuously. One critical idea used in
all previous works on LP [Vai89b, CLS19, LSZ19] is to delay low-rank updates
on $(A^{\top}WA)^{-1}$. However, in the cutting plane method, the low rank
updates to $A$ cannot be delayed. Thus it appears that previous techniques are
inapplicable. Moreover, $n^{2}$ lower bounds have recently been established
for natural matrix maintenance tasks (e.g. determinant, inverse) with
row/column insertions/deletions under the Online Matrix-Vector conjecture
(e.g. [HKNS15, BNS19]). Therefore, we believe our algorithm is tight and
conjecture the following:

######  Conjecture 1.7.

Solving the feasibility problem requires $\Omega(n^{3}\log(\kappa))$ time.
Hence our cutting plane method achieves the optimal runtime.

###  1.6 Related works

##### Leverage scores.

Leverage scores are a fundamental concept in graph problems and numerical
linear algebra. There are many works about how to approximate leverage scores
[SS11, DMIMW12, CW13, NN13] or more general version of leverages, e.g. Lewis
weights [Lew78, BLM89, CP15] and ridge leverage scores [CMM17]. From graph
perspective, it was applied to solve max-flow [Mad13, Mad16], generate random
spanning trees [Sch18], and sparsify graphs [SS11]. From matrix perspective,
it was used to give matrix CUR decomposition [BW14, SWZ17, SWZ19] and tensor
CURT decomposition [SWZ19]. From optimization perspective, it was used for
approximating the John Ellipsoid [CCLY19], accelerating the kernel methods
[AKM+17, AKM+19], showing the convergence of the deep neural network [LSS+20],
cutting plane methods, e.g. [Vai89a, LSW15] and this paper.

##### Linear Program.

Linear Program is a fundamental problem in convex optimization and can be
treated as an special case where one can apply the cutting plane method. There
is a super long list of work focused on fast algorithms for linear program
[Dan47, Kha80, Kar84, Vai87, Vai89b, LS14, LS15, Sid15, Lee16, CLS19, LSZ19,
Son19, Bra20, BLSS20].

##### Membership oracle.

Besides the separation oracle considered in this paper, there is another line
of work on using the membership oracle to solve the feasibility problem
[Pro96, KV06, LV06, GLS12, LSV18]. For a query point $x$, this oracle outputs
$x\in K$ or $x\notin K$.

##  2 Our Techniques

Section 1.4 provided a quick overview. Here we take a deeper dive into our
techniques.

###  2.1 Efficient approximation of changes in leverage scores

The key to fast maintenance of leverage scores is an efficient way to
approximate their changes between consecutive steps. While a fine-grained
approximation leads to an accurate approximation, the time to compute such an
approximation might be unaffordable. On the other hand, a coarse-grained
approximation can be efficiently computed, but might lead to accumulating
errors that blow up after a small number of steps. This leads to a tradeoff
between accuracy and efficiency.

Central to our data structure are a coarse-grained formula (Lemma 7.4) and a
fine-grained formula for the change in leverage scores (Lemma 8.4). While the
coarse-grained formula approximates the leverage score’s change via a single
integral, the fine-grained formula is a cocktail involving integrals, matrix
inverse and matrix multiplication. For the coarse-grained formula, we simply
estimate the integral by a single point along the integral (Lemma 7.5). For
the fine-grained formula, however, the approximation (Lemma 8.10) is more
involved and requires appropreiate numerical tools.

Our course-grained and fine-grained formulas lead to two different data
structures for leverage score maintenance: a simple deterministic data
structure that achieves low running time but introduces a large error in each
step, and a more complicated randomized data structure that incurs small error
each step at the cost of efficiency.

Observe that we are allowed to recompute the leverage scores exactly after
every $O(n^{\omega-2+o(1)})$ steps as exact computation of leverage scores
takes $O(n^{\omega+o(1)})$ time. Therefore it suffices to control the error
accumulation in $O(n^{\omega-2+o(1)})$ steps. To achieve this, both the simple
and complicated data structures are crucial. While the simple data structure
achieves the desired $O(n^{2})$ time per step, the error accrued in the
$O(n^{\omega-2+o(1)})$ steps are too large for the cutting plane method. For
the complicated data structure, we are able to achieve $\ell_{2}$-error $o(1)$
in $O(n^{\omega-2+o(1)})$ steps, but the amortized running time would be
$O(n^{2+o(1)})$ per step.

###  2.2 Layered data structure

To achieve the best of both worlds of the simple and complicated data
structures, we propose a data structure that combines both data structures and
interpolates between accuracy and efficiency. Specifically, our data structure
approximates changes in leverage scores to $\ell_{2}$-error within
$1/\log^{O(1)}(n)$ which would suffice for Vaidya’s cutting plane method and
handles each iteration in amortized $O(n^{2})$ time.

Our data structure employs a layered approach (see Figure 1) where different
layers are associated with different error tolerances. The more inner a layer
is, the more error it can tolerate and the faster is the runtime. Whenever the
error accumulated in a layer becomes too high, the layer above would take over
and produce a finer error estimate which would, of course, be more time
costly. But because of our layered approach, the higher layer is called on
less often and can afford to spend more time. More specifically, our layered
data structure contains three layers. The inner and middle layers both employ
the simple data structure to achieve computational efficiency while the outer
layer uses the complicated data structure to ensure a low error.

Our layered approach builds on several fundamental results on matrix
multiplication which we summarize in Theorem 2.1. We remark that Vaidya
[Vai89b] used the first result, a recent LP solver [CLS19] used the first two,
while our cutting plane method crucially depends on all the results in the
table.

######  Theorem 2.1 (Fast rectangular matrix multiplication results).

For any $n,r>0$, denote $\mathcal{T}_{\mathrm{mat}}(n,n,r)$ the time to
compute the multiplication of an $n\times n$ matrix and an $n\times r$
matrix111Note that
$\mathcal{T}_{\mathrm{mat}}(n,n,r)=\mathcal{T}_{\mathrm{mat}}(n,r,n)$. Then we
have:

Reference | $r$ | $\mathcal{T}_{\mathrm{mat}}(n,n,r)$ | Layer
---|---|---|---
[LG14] | $n$ | $O(n^{\omega+o(1)})$ | Restart
[GU18] | $n^{0.31}$ | $O(n^{2+o(1)})$ | Outer layer
[Cop82] | $n^{0.17}$ | $O(n^{2}\log^{2}n)$ | Middle layer
[BD76] | $\log^{O(1)}n$ | $O(n^{2})$ | Inner layer

Our data structure also draws on various numerical tools such as fast
rectangular matrix multiplication, “tall” JL projection, preconditioning,
inverse maintenance, and polynomial interpolation for approximating integrals,
which we discuss in more details below.

Figure 1: Illustration of our three-level data structure with
$T_{\text{inn}}=2$, $T_{\text{mid}}=3$ and $T_{\text{out}}=4$ approximating
the leverage scores of the sequence $\\{w^{0},w^{1},\cdots\\}$. The three
levels maintain three approximate sequences
$\\{v_{\text{inn}}^{0},v_{\text{inn}}^{1},\cdots\\}$,
$\\{v_{\text{mid}}^{0},v_{\text{mid}}^{1},\cdots\\}$ and
$\\{v_{\text{out}}^{0},v_{\text{out}}^{1},\cdots\\}$, with errors
$\epsilon_{\text{inn}}=\left\lVert\log(w)-\log(v_{\text{inn}})\right\rVert_{\infty}$,
$\epsilon_{\text{mid}}=\left\lVert\log(w)-\log(v_{\text{mid}})\right\rVert_{\infty}$
and
$\epsilon_{\text{out}}=\left\lVert\log(w)-\log(v_{\text{out}})\right\rVert_{\infty}$
that satisfy
$1\gg\epsilon_{\text{inn}}\gg\epsilon_{\text{mid}}\gg\epsilon_{\text{out}}>0$.
The inner level takes a step for every $w$-update, the middle step takes a
step in every $T_{\operatorname{inn}}$ inner steps, and the outer step takes a
step in every $T_{\mathrm{mid}}$ middle steps. The entire data structure is
restarted after $T_{\operatorname{out}}$ outer steps. Each middle step refines
the inner approximation of leverage scores, and each outer step refines the
approximation of both the inner and middle approximations of the leverage
scores. For the actual choice of the parameters $T_{\text{inn}}$,
$T_{\text{mid}}$, $T_{\text{out}}$ and $\epsilon_{\text{inn}}$,
$\epsilon_{\text{mid}}$, $\epsilon_{\text{out}}$ in our data structure, see
Table 5.

###  2.3 Batched low-rank update

The layered approach in the previous subsection leads to the following
technical hurdle. While the inner layer reacts to each update, the middle and
outer layers perform one step only after a certain number of updates. In the
cutting plane method, each update of the vector $w\in\mathbb{R}_{+}^{m}$
satisfies an $\ell_{2}$-closeness property:
$\left\lVert\log(w^{\operatorname{new}})-\log(w)\right\rVert_{2}=O(1)$. For a
sequence of $T$ updates $w^{(0)},w^{(1)},\cdots,w^{(T)}$, the
$\ell_{2}$-closeness property is no longer satisfied for $w^{(0)}$ and
$w^{(T)}$, and the weight changes between $w^{(0)}$ and $w^{(T)}$ can be
exponential in $T$.

To resolve this issue, we employ a batched low-rank update method that
computes a vector $v^{\mathrm{mid}}$ that is $\ell_{\infty}$-close to the
vector $w^{(T)}$ and an accurate estimate of the leverage score change
$\sigma(v^{\mathrm{mid}})-\sigma(w^{0})$. This accurate estimate is done by a
low-rank update rule for matrices (Woodbury matrix identity). Unfortunately,
the Woodbury formula involves the inverse of certain matrices but the inverse
maintenance data structure only maintains the inverse of a nearby matrix.
Therefore, we need to use the maintained inverse as a preconditioner to
approximately compute the Woodbury formula. Due to the exponential changes in
the weight, we need to solve certain linear systems to exponential accuracy
and it can be too expensive even with a very good preconditioner (Lemma 6.2).

To rescue this preconditioner idea, we transform and split the sequence into
pieces of size $\operatorname{poly}\log(n)$. We ensure the weight changes by
only a quasi-polynomial factor and this decreases the cost of solving linear
systems to $\operatorname{poly}\log(n)$ steps. Since we batch the task of
handling $\operatorname{poly}\log(n)$ weight changes into one rectangular
matrix multiplication which can be performed in $O(n^{2})$ time, we make sure
the cost per weight change is exactly $O(n^{2})$ time (Theorem 6.3).

###  2.4 Illustration of our analysis

We describe the numerical tools used in our analysis, and provide simple
illustrations of our applications of these tools. The actual way in which they
are used in Sections 7-8 are more involved.

##### Discrete sampling for multiple variable integrals

The most standard way to approximate an integral is by discretization, which
takes a weighted sum of the integrand over a set of points in the domain.
Unlike common discretization tools like the trapezoidal method, for our
purpose we need to interpolate multiple variable integral using a polynomial
for higher accuracy (Theorem 3.10). We give a simple example to illustrate our
application of polynomial interpolation for multiple variable integrals as
follows. In our fine-grained formula of the leverage score change from $w$ to
$w^{\operatorname{new}}$, one of the integral terms is

$\displaystyle\sigma_{i,\textrm{cts}}=\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}\gamma_{i,s,t}^{\top}\gamma_{i,s^{\prime},t}\mathrm{d}s\mathrm{d}s^{\prime}\mathrm{d}t,$

where

$\displaystyle\gamma_{i,s,t}=$ $\displaystyle\leavevmode\nobreak\
\sqrt{W^{\mathrm{mid}}-W^{\operatorname{new}}}\cdot
Q(y_{s,t})\cdot(Z_{t}-X_{t})\cdot
Q(y_{s,t})\cdot\sqrt{W^{\operatorname{new}}}\cdot e_{i}.$

In order to approximate such an integral, we take a set
$\mathcal{T}\subseteq[0,1]$ of $N=\log^{O(1)}(n)$ points along the integration
together with weights $\\{\omega_{t}\\}_{t\in\mathcal{T}}$ as in Theorem 3.10,
and approximate the integral by

$\displaystyle\sigma_{i,\textrm{dis}}=\sum_{t\in\mathcal{T}}\sum_{s\in\mathcal{S}}\sum_{s^{\prime}\in\mathcal{S}}\omega_{t}\omega_{s}\omega_{s^{\prime}}\gamma_{i,s,t}^{\top}\gamma_{i,s^{\prime},t}.$

Theorem 3.10 then shows that the $\ell_{2}$-error can be bounded as
$\left\lVert\sigma_{\textrm{cts}}-\sigma_{\textrm{dis}}\right\rVert_{2}\leq\operatorname{poly}(n)/2^{2N},$
which is negligible by our choice of $N=\log^{O(1)}(n)$.

##### Projection maintenance and preconditioning

Inverse maintenance was first proposed in [Kha80] as a method for solving
“slowly-changing” linear system.

Given a sequence of positive vectors $w^{1},w^{2},\cdots
w^{T}\in\mathbb{R}_{+}^{m}$, let $P(w)$ be defined as
$P(w)=\sqrt{W}A(A^{\top}WA)^{-1}A^{\top}\sqrt{W}$, where
$W\in\mathbb{R}^{m\times m}$ is the diagonal matrix with $W_{i,i}=w_{i}$. The
goal is to output a sequence of vectors $v^{1},v^{2},\cdots,v^{T}$ such that
$(1-\epsilon)v^{t}\leq w^{t}\leq(1+\epsilon)v^{t}$, $\forall t\in[T]$, and
efficiently computes $P(v^{t})u$, for query vector $u\in\mathbb{R}^{n}$. The
recent work of [CLS19] gave an efficient way to perform such a task. For our
purpose, however, rather than just computing $P(v^{t})u$, we also need
explicit approximations to the matrices $Q(w)=A(A^{\top}WA)^{-1}A^{\top}$ and
$M(w)^{-1}=(A^{\top}WA)^{-1}$. The matrices $Q(w)$ and $M(w)^{-1}$ appear
frequently in our formula for the changes in leverage scores, and we need
their approximations as pre-conditioners for accelerating the computation of
certain matrix rectangular multiplication involving $Q(w)$ and $M(w)^{-1}$
(Lemma 6.1).

##### “Tall” JL & fast rectangular matrix multiplication

From the previous paragraph on discrete sampling, it suffices to compute
$\gamma_{i,s,t}$ for all $i$, where

$\displaystyle\gamma_{i,s,t}=$ $\displaystyle\leavevmode\nobreak\
\sqrt{W^{\mathrm{mid}}-W^{\operatorname{new}}}\cdot
Q(y_{s,t})\cdot(Z_{t}-X_{t})\cdot
Q(y_{s,t})\cdot\sqrt{W^{\operatorname{new}}}\cdot e_{i}.$

Notice that computing $\gamma_{i,s,t}$ for all $i$ is essentially computing
the matrix products

$\displaystyle\sqrt{W^{\mathrm{mid}}-W^{\operatorname{new}}}\cdot
Q(y_{s,t})\cdot(Z_{t}-X_{t})\cdot
Q(y_{s,t})\cdot\sqrt{W^{\operatorname{new}}},$

which would take $O(n^{\omega+o(1)})$ time if computed exactly. To improve the
time while ensuring keeping the error small, we invoke JL with dimension
$n^{c}$ for small constant $c$ by computing

$\displaystyle\sigma_{i,\textrm{jl}}=\sum_{t\in\mathcal{T}}\sum_{s\in\mathcal{S}}\sum_{s^{\prime}\in\mathcal{S}}\omega_{t}\omega_{s}\omega_{s^{\prime}}\gamma_{i,s,t}^{\top}R_{\gamma,s,s^{\prime},t}^{\top}R_{\gamma,s,s^{\prime},t}\gamma_{i,s^{\prime},t},$

where $R_{\gamma,s,s^{\prime},t}\in\mathbb{R}^{n^{c}\times n}$ is a random
matrix. It is essential that $c$ is picked such that the rectangular matrix
multiplication can be done in time roughly $n^{2}$. To obtain a small error,
$\sigma_{i,\textrm{jl}}$ should be a good estimate with a small variance. We
note that $\sigma_{i,\textrm{jl}}$ is indeed an unbiased estimator of
$\sigma_{i,\textrm{dis}}$ and its variance can be bounded as
$\sum_{i=1}^{m}\operatorname*{{\bf{Var}}}[\sigma_{i,\textrm{jl}}]\leq
O(\epsilon^{2}/n^{c}),$ where $\epsilon$ is the error for the projection
maintenance used by the data structure (Sections 8.4 and 8.5). Leveraging fast
rectangular matrix multiplication, we pick $c=0.31$ and $\epsilon=n^{-0.1}$.
The variance would then be bounded by $n^{-0.51}$ which is sufficiently small
after $n^{\omega-2+o(1)}$ steps before the data structure restarts.

###  2.5 Much faster rectangular matrix multiplication implies deterministic
cutting plane method

Our algorithm crucially relies on different kinds of fast rectangular matrix
multiplication results. We also show that if these results are improved, then
we are able to get a deterministic cutting plane method immediately.

######  Corollary 2.2 (Informal version of Corollary 7.3).

If $\mathcal{T}_{\mathrm{mat}}(n,n,r)=O(n^{2}\log^{O(1)}(n))$ for
$r=n^{\beta}$ with $\beta>2/3$, then there is a deterministic cutting plane
method which runs in time

$\displaystyle O(n\cdot\operatorname{SO}\log(\kappa)+n^{3}\log(\kappa)),$

where $\operatorname{SO}$ is the time complexity of the separation oracle.

Let $\alpha$ denote the dual exponent of matrix multiplication, which is the
largest number $\alpha>0$ such that ${\cal
T}_{\mathrm{mat}}(n,n,n^{\alpha})=n^{2+o(1)}$. Let $\beta$ denote the largest
number such that ${\cal T}_{\mathrm{mat}}(n,n,n^{\beta})=n^{2}\log^{O(1)}n$. A
very recent result by Christandl, Le Gall, Lysikov and Zuiddam [CGLZ20] showed
the limitations of several tensor techniques: they proved that $\alpha<0.625$
for certain tensors. We believe our work initiated two interesting open
questions in the area of fast matrix multiplication: (1) whether one can prove
a better upper bound on $\beta$ (compared to $\alpha$) for certain tensor
techniques, and (2) if there is a non-trivial inequality between $\beta$ and
$\alpha$.

Organization. We introduce basic notations, backgrounds and tools in Section
3. In Section 4, we present the statement that Vaidya’s cutting plane method
tolerates perturbed leverage scores. We present our main data-structure for
maintaining leverage scores in Section 5. Our main data-structure uses two
different leverage score maintenance data-structures in Sections 7 and 8 with
three different settings of the error parameter. Both our data-structures in
Section 7 and 8 reply on the batched low rank algorithm in Section 6.

In Section A, we prove that Vaidya’s cutting plane method tolerates perturbed
leverage scores. We provide several modified versions of projection
maintenance data-structures in Section B. Finally, we explain how to handle
convex-concave game optimization in Section C, and present our applications to
market equilibrium computations in Section D.

## Acknowledgments

The authors would like to express their sincere gratitude to matrix
multiplicationer Josh Alman for his patient and answers of our exponential
number of questions about fast matrix multiplication.

The authors would like to thank Swati Padmanabhan for very useful discussions
at the early stage of this project. The authors would like to thank Lijie
Chen, Nikhil Devanur, Irit Dinur, Simon S. Du, Wei Hu, Jason Lee, Jerry Li,
Ruoqi Shen, Aaron Schild, Aaron Sidford, Santosh Vempala, Xin Yang, Peilin
Zhong, and Danyang Zhuo.

The authors would like to thank Josh Alman, Irit Dinur, Avi Wigderson, and
Jeroen Zuiddam for useful discussion about limitation of fast matrix
multiplication.

The authors would like to thank Sanjeev Arora, Alexandr Andoni, Ainesh Bakshi,
Yangsibo Huang, Rajesh Jayaram, Ravindran Kanna, Michael Kapralov, Adam
Klivans, Kai Li, Christos Papadimitriou, Eric Price, Daniel Roy, Clifford
Stein, Omri Weinstein, David P. Woodruff, and Hengjie Zhang for asking
interesting questions at the end of the talk of this paper.

This project was supported in part by NSF awards CCF-1749609, CCF-1740551,
DMS-1839116, and Microsoft Research Faculty Fellowship.

This project was supported in part by Special Year on Optimization,
Statistics, and Theoretical Machine Learning (being led by Sanjeev Arora) at
Institute for Advanced Study.

## References

  * [AD54] Kenneth J Arrow and Gerard Debreu.  Existence of an equilibrium for a competitive economy.  Econometrica: Journal of the Econometric Society, pages 265–290, 1954.
  * [AKM+17] Haim Avron, Michael Kapralov, Cameron Musco, Christopher Musco, Ameya Velingker, and Amir Zandieh.  Random Fourier features for kernel ridge regression: Approximation bounds and statistical guarantees.  In International Conference on Machine Learning (ICML). https://arxiv.org/pdf/1804.09893.pdf, 2017.
  * [AKM+19] Haim Avron, Michael Kapralov, Cameron Musco, Christopher Musco, Ameya Velingker, and Amir Zandieh.  A universal sampling method for reconstructing signals with simple Fourier transforms.  In 51st Annual ACM SIGACT Symposium on Theory of Computing (STOC). https://arxiv.org/pdf/1812.08723.pdf, 2019.
  * [AM02] Lawrence M Ausubel and Paul R Milgrom.  Ascending auctions with package bidding.  Advances in Theoretical Economics, 1(1), 2002.
  * [AV95] David S Atkinson and Pravin M Vaidya.  A cutting plane algorithm for convex programming that uses analytic centers.  Mathematical Programming, 69(1-3):1–43, 1995.
  * [BD76] Roger W Brockett and David Dobkin.  On the number of multiplications required for matrix multiplication.  SIAM Journal on Computing, 5(4):624–628, 1976.
  * [BDX10] Benjamin Birnbaum, N Devanur, and Lin Xiao.  New convex programs and distributed algorithms for fisher markets with linear and spending constraint utilities.  Unpublished manuscript, 2010.
  * [BLM89] Jean Bourgain, Joram Lindenstrauss, and V Milman.  Approximation of zonoids by zonotopes.  Acta mathematica, 162(1):73–141, 1989.
  * [BLSS20] Jan van den Brand, Yin Tat Lee, Aaron Sidford, and Zhao Song.  Solving tall dense linear programs in nearly linear time.  In 52nd Annual ACM SIGACT Symposium on Theory of Computing (STOC). https://arxiv.org/pdf/2002.02304.pdf, 2020.
  * [BNS19] Jan van den Brand, Danupon Nanongkai, and Thatchaphol Saranurak.  Dynamic matrix inverse: Improved algorithms and matching conditional lower bounds.  In 60th Annual IEEE Symposium on Foundations of Computer Science (FOCS). https://arxiv.org/pdf/1905.05067.pdf, 2019.
  * [Bra20] Jan van den Brand.  A deterministic linear program solver in current matrix multiplication time.  In ACM-SIAM Symposium on Discrete Algorithms (SODA). https://arxiv.org/pdf/1910.11957.pdf, 2020.
  * [BV02] Dimitris Bertsimas and Santosh Vempala.  Solving convex programs by random walks.  In Proceedings of the thiry-fourth annual ACM symposium on Theory of computing (STOC), pages 109–115. ACM, 2002.
  * [BW14] Christos Boutsidis and David P Woodruff.  Optimal cur matrix decompositions.  In Proceedings of the 46th Annual ACM Symposium on Theory of Computing (STOC), pages 353–362. ACM, https://arxiv.org/pdf/1405.7910, 2014.
  * [CCLY19] Michael B Cohen, Ben Cousins, Yin Tat Lee, and Xin Yang.  A near-optimal algorithm for approximating the John ellipsoid.  In Annual Conference on Learning Theory (COLT). https://arxiv.org/pdf/1905.11580.pdf, 2019.
  * [CGLZ20] Matthias Christandl, François Le Gall, Vladimir Lysikov, and Jeroen Zuiddam.  Barriers for rectangular matrix multiplication.  In arXiv preprint. https://arxiv.org/2003.03019.pdf, 2020\.
  * [CLS19] Michael B Cohen, Yin Tat Lee, and Zhao Song.  Solving linear programs in the current matrix multiplication time.  In Proceedings of the 51st Annual ACM Symposium on Theory of Computing (STOC). https://arxiv.org/pdf/1810.07896.pdf, 2019.
  * [CMM17] Michael B Cohen, Cameron Musco, and Christopher Musco.  Input sparsity time low-rank approximation via ridge leverage score sampling.  In Proceedings of the Twenty-Eighth Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), pages 1758–1777. SIAM, https://arxiv.org/pdf/1511.07263.pdf, 2017.
  * [Cop82] Don Coppersmith.  Rapid multiplication of rectangular matrices.  SIAM Journal on Computing, 11(3):467–471, 1982.
  * [Cor89] Bernard Cornet.  Linear exchange economies.  Cahier Eco-Math, Université de Paris, 1, 1989.
  * [CP15] Michael B. Cohen and Richard Peng.  $\ell_{p}$ row sampling by lewis weights.  In Proceedings of the Forty-Seventh Annual ACM on Symposium on Theory of Computing (STOC), pages 183–192. https://arxiv.org/pdf/1412.0588, 2015.
  * [CW87] Don Coppersmith and Shmuel Winograd.  Matrix multiplication via arithmetic progressions.  In Proceedings of the nineteenth annual ACM symposium on Theory of computing (STOC), pages 1–6. ACM, 1987.
  * [CW13] Kenneth L. Clarkson and David P. Woodruff.  Low rank approximation and regression in input sparsity time.  In Symposium on Theory of Computing Conference (STOC), pages 81–90. https://arxiv.org/pdf/1207.6365, 2013.
  * [Dan47] George B Dantzig.  Maximization of a linear function of variables subject to linear inequalities.  Activity analysis of production and allocation, 13:339–347, 1947\.
  * [DGM16] Ran Duan, Jugal Garg, and Kurt Mehlhorn.  An improved combinatorial polynomial algorithm for the linear arrow-debreu market.  In Proceedings of the twenty-seventh annual ACM-SIAM symposium on Discrete algorithms (SODA), pages 90–106. SIAM, https://arxiv.org/pdf/1510.02694.pdf, 2016.
  * [DGV16] Nikhil R Devanur, Jugal Garg, and László A Végh.  A rational convex program for linear arrow-debreu markets.  volume 5(1), page 6. https://arxiv.org/pdf/1307.8037.pdf, 2016.
  * [DM15] Ran Duan and Kurt Mehlhorn.  A combinatorial polynomial algorithm for the linear arrow–debreu market.  Information and Computation, 243:112–132, 2015.
  * [DMIMW12] Petros Drineas, Malik Magdon-Ismail, Michael W Mahoney, and David P Woodruff.  Fast approximation of matrix coherence and statistical leverage.  Journal of Machine Learning Research, 13(Dec):3475–3506, 2012.
  * [DPSV08] Nikhil R Devanur, Christos H Papadimitriou, Amin Saberi, and Vijay V Vazirani.  Market equilibrium via a primal–dual algorithm for a convex program.  Journal of the ACM (JACM), 55(5):22, 2008.
  * [DS13] Alexander Munro Davie and Andrew James Stothers.  Improved bound for complexity of matrix multiplication.  Proceedings of the Royal Society of Edinburgh Section A: Mathematics, 143(2):351–369, 2013.
  * [dVSV07] Sven de Vries, James Schummer, and Rakesh V Vohra.  On ascending vickrey auctions for heterogeneous objects.  Journal of Economic Theory, 132(1):95–118, 2007.
  * [Eav75] B Curtis Eaves.  A finite algorithm for the linear exchange model.  Technical report, STANFORD UNIV CALIF SYSTEMS OPTIMIZATION LAB, 1975.
  * [GLS12] Martin Grötschel, László Lovász, and Alexander Schrijver.  Geometric algorithms and combinatorial optimization, volume 2.  Springer Science & Business Media, 2012.
  * [GU18] Francois Le Gall and Florent Urrutia.  Improved rectangular matrix multiplication using powers of the coppersmith-winograd tensor.  In Proceedings of the Twenty-Ninth Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), pages 1029–1046. SIAM, https://arxiv.org/pdf/1708.05622.pdf, 2018.
  * [GV19] Jugal Garg and László A Végh.  A strongly polynomial algorithm for linear exchange markets.  In Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing (STOC), pages 54–65. ACM, https://arxiv.org/pdf/1809.06266.pdf, 2019.
  * [HKNS15] Monika Henzinger, Sebastian Krinninger, Danupon Nanongkai, and Thatchaphol Saranurak.  Unifying 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 (STOC), pages 21–30. ACM, https://arxiv.org/pdf/1511.06773.pdf, 2015.
  * [IK94] Eugene Isaacson and Herbert Bishop Keller.  Analysis of numerical methods.  1994\.
  * [Jai07] Kamal Jain.  A polynomial time algorithm for computing an arrow–debreu market equilibrium for linear utilities.  SIAM Journal on Computing, 37(1):303–318, 2007.
  * [JL84] William B Johnson and Joram Lindenstrauss.  Extensions of lipschitz mappings into a hilbert space.  Contemporary mathematics, 26(189-206):1, 1984.
  * [Kar84] Narendra Karmarkar.  A new polynomial-time algorithm for linear programming.  In Proceedings of the sixteenth annual ACM symposium on Theory of computing (STOC), pages 302–311. ACM, 1984.
  * [Kha80] Leonid G Khachiyan.  Polynomial algorithms in linear programming.  USSR Computational Mathematics and Mathematical Physics, 20(1):53–72, 1980.
  * [KJC82] Alexander S Kelso Jr and Vincent P Crawford.  Job matching, coalition formation, and gross substitutes.  Econometrica: Journal of the Econometric Society, pages 1483–1504, 1982.
  * [KTE88] Leonid G Khachiyan, Sergei Pavlovich Tarasov, and I. I. Erlikh.  The method of inscribed ellipsoids.  In Soviet Math. Dokl, volume 37, pages 226–230, 1988.
  * [KV06] Adam Tauman Kalai and Santosh Vempala.  Simulated annealing for convex optimization.  Mathematics of Operations Research, 31(2):253–266, 2006.
  * [LDFU13] Yichao Lu, Paramveer Dhillon, Dean P Foster, and Lyle Ungar.  Faster ridge regression via the subsampled randomized hadamard transform.  In Advances in neural information processing systems (NIPS), pages 369–377, 2013.
  * [Lee16] Yin Tat Lee.  Faster algorithms for convex and combinatorial optimization.  PhD thesis, Massachusetts Institute of Technology, 2016.
  * [Lew78] D. Lewis.  Finite dimensional subspaces of $\ell_{p}$.  Studia Mathematica, 63(2):207–212, 1978.
  * [LG14] François Le Gall.  Powers of tensors and fast matrix multiplication.  In Proceedings of the 39th international symposium on symbolic and algebraic computation (ISSAC), pages 296–303. ACM, 2014.
  * [LS14] Yin Tat Lee and Aaron Sidford.  Path finding methods for linear programming: Solving linear programs in ${O}(\sqrt{rank})$ iterations and faster algorithms for maximum flow.  In 55th Annual IEEE Symposium on Foundations of Computer Science (FOCS), pages 424–433. https://arxiv.org/pdf/1312.6677.pdf, https://arxiv.org/pdf/1312.6713.pdf, 2014.
  * [LS15] Yin Tat Lee and Aaron Sidford.  Efficient inverse maintenance and faster algorithms for linear programming.  In 56th Annual IEEE Symposium on Foundations of Computer Science (FOCS), pages 230–249. https://arxiv.org/pdf/1503.01752.pdf, 2015.
  * [LSS+20] Jason Lee, Ruoqi Shen, Zhao Song, Mengdi Wang, and Zheng Yu.  Leverage score sampling, neural tangent kernel, and kernel ridge regression.  In Manuscript, 2020.
  * [LSV18] Yin Tat Lee, Aaron Sidford, and Santosh S Vempala.  Efficient convex optimization with membership oracles.  In Annual Conference on Learning Theory (COLT). https://arxiv.org/pdf/1706.07357.pdf, 2018.
  * [LSW15] Yin Tat Lee, Aaron Sidford, and Sam Chiu-wai Wong.  A faster cutting plane method and its implications for combinatorial and convex optimization.  In 56th Annual IEEE Symposium on Foundations of Computer Science (FOCS), pages 1049–1065. https://arxiv.org/pdf/1508.04874.pdf, 2015.
  * [LSZ19] Yin Tat Lee, Zhao Song, and Qiuyi Zhang.  Solving empirical risk minimization in the current matrix multiplication time.  In Annual Conference on Learning Theory (COLT). https://arxiv.org/pdf/1905.04447, 2019.
  * [LV06] László Lovász and Santosh Vempala.  Fast algorithms for logconcave functions: Sampling, rounding, integration and optimization.  In 47th Annual IEEE Symposium on Foundations of Computer Science (FOCS), pages 57–68. IEEE, 2006.
  * [LW17] Renato Paes Leme and Sam Chiu-wai Wong.  Computing walrasian equilibria: fast algorithms and structural properties.  In Proceedings of the Twenty-Eighth Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), pages 632–651. SIAM, https://arxiv.org/pdf/1511.04032.pdf, 2017.
  * [Mad13] Aleksander Madry.  Navigating central path with electrical flows: From flows to matchings, and back.  In 54th Annual IEEE Symposium on Foundations of Computer Science (FOCS), pages 253–262. https://arxiv.org/pdf/1307.2205.pdf, 2013.
  * [Mad16] Aleksander Madry.  Computing maximum flow with augmenting electrical flows.  In 2016 IEEE 57th Annual Symposium on Foundations of Computer Science (FOCS), pages 593–602. IEEE, https://arxiv.org/pdf/1608.06016.pdf, 2016.
  * [Nem94] Arkadi Nemirovski.  Efficient methods in convex programming.  Lecture notes, 1994.
  * [Nem95] Arkadi Nemirovski.  Information-based complexity of convex programming.  Lecture Notes, 1995.
  * [NN89] YE Nesterov and AS Nemirovskii.  Self-concordant functions and polynomial time methods in convex programming. preprint, central economic & mathematical institute, ussr acad.  Sci. Moscow, USSR, 1989.
  * [NN13] Jelani Nelson and Huy L Nguyên.  OSNAP: Faster numerical linear algebra algorithms via sparser subspace embeddings.  In 54th Annual IEEE Symposium on Foundations of Computer Science (FOCS), pages 117–126. IEEE, https://arxiv.org/pdf/1211.1002, 2013.
  * [NP83] EI Nenakov and ME Primak.  One algorithm for finding solutions of the arrow-debreu model.  Kibernetica, 3:127–128, 1983.
  * [NY83] Arkadii Semenovich Nemirovsky and David Borisovich Yudin.  Problem complexity and method efficiency in optimization.  1983\.
  * [Par99] David C Parkes.  ibundle: An efficient ascending price bundle auction.  1999\.
  * [Pro96] V Yu Protasov.  Algorithms for approximate calculation of the minimum of a convex function from its values.  Mathematical Notes, 59(1):69–74, 1996.
  * [PSW17] Eric Price, Zhao Song, and David P. Woodruff.  Fast regression with an ${\ell}_{\infty}$ guarantee.  In International Colloquium on Automata, Languages, and Programming (ICALP), 2017.
  * [PU02] David C Parkes and Lyle H Ungar.  An ascending-price generalized vickrey auction.  2002\.
  * [RSW16] Ilya Razenshteyn, Zhao Song, and David P Woodruff.  Weighted low rank approximations with provable guarantees.  In Proceedings of the 48th Annual Symposium on the Theory of Computing (STOC), 2016.
  * [San04] Piotr Sankowski.  Dynamic transitive closure via dynamic matrix inverse.  In 45th Annual IEEE Symposium on Foundations of Computer Science (FOCS), pages 509–517. IEEE, 2004.
  * [Sch18] Aaron Schild.  An almost-linear time algorithm for uniform random spanning tree generation.  In Proceedings of the 50th Annual ACM SIGACT Symposium on Theory of Computing (STOC), pages 214–227. ACM, https://arxiv.org/pdf/1711.06455.pdf, 2018.
  * [Sho77] Naum Z Shor.  Cut-off method with space extension in convex programming problems.  Cybernetics and systems analysis, 13(1):94–96, 1977.
  * [Sid15] Aaron Daniel Sidford.  Iterative methods, combinatorial optimization, and linear programming beyond the universal barrier.  PhD thesis, Massachusetts Institute of Technology, 2015.
  * [Son19] Zhao Song.  Matrix Theory : Optimization, Concentration and Algorithms.  PhD thesis, The University of Texas at Austin, 2019.
  * [SS11] Daniel A Spielman and Nikhil Srivastava.  Graph sparsification by effective resistances.  In SIAM Journal on Computing, volume 40(6), pages 1913–1926. https://arxiv.org/pdf/0803.0929.pdf, 2011.
  * [SWZ17] Zhao Song, David P Woodruff, and Peilin Zhong.  Low rank approximation with entrywise $\ell_{1}$-norm error.  In Proceedings of the 49th Annual Symposium on the Theory of Computing (STOC). ACM, https://arxiv.org/pdf/1611.00898, 2017.
  * [SWZ19] Zhao Song, David P Woodruff, and Peilin Zhong.  Relative error tensor low rank approximation.  In ACM-SIAM Symposium on Discrete Algorithms (SODA). https://arxiv.org/pdf/1704.08246, 2019.
  * [Vai87] Pravin M Vaidya.  An algorithm for linear programming which requires $o(((m+n)n^{2}+(m+n)^{1.5}n)l)$ arithmetic operations.  In 28th Annual IEEE Symposium on Foundations of Computer Science (FOCS), 1987.
  * [Vai89a] Pravin M Vaidya.  A new algorithm for minimizing convex functions over convex sets.  In 30th Annual IEEE Symposium on Foundations of Computer Science (FOCS), pages 338–343, 1989.
  * [Vai89b] Pravin M Vaidya.  Speeding-up linear programming using fast matrix multiplication.  In 30th Annual Symposium on Foundations of Computer Science (FOCS), pages 332–337. IEEE, 1989.
  * [Vaz10] Vijay V Vazirani.  Spending constraint utilities with applications to the adwords market.  Mathematics of Operations Research, 35(2):458–478, 2010.
  * [Vég16] László A Végh.  A strongly polynomial algorithm for a class of minimum-cost flow problems with separable convex objectives.  volume 45(5), pages 1729–1761. https://arxiv.org/pdf/1110.4882.pdf, 2016.
  * [Wan16] Zi Wang.  A New Strongly Polynomial Algorithm for Computing Fisher Market Equilibria with Spending Constraint Utilities.  PhD thesis, New York University, 2016.
  * [Wil12] Virginia Vassilevska Williams.  Multiplying matrices faster than coppersmith-winograd.  In Proceedings of the forty-fourth annual ACM symposium on Theory of computing (STOC), pages 887–898. ACM, 2012.
  * [Woo49] Max A Woodbury.  The stability of out-input matrices.  Chicago, IL, 9, 1949.
  * [Woo50] Max A Woodbury.  Inverting modified matrices.  Memorandum report, 42(106):336, 1950.
  * [Ye08] Yinyu Ye.  A path to the arrow–debreu competitive market equilibrium.  Mathematical Programming, 111(1-2):315–348, 2008.
  * [YN76] David B Yudin and Arkadii S Nemirovski.  Evaluation of the information complexity of mathematical programming problems.  Ekonomika i Matematicheskie Metody, 12:128–142, 1976.

##  3 Preliminaries

In this section we introduce the notions and tools used throughout this paper.

###  3.1 Notations

For a positive integer $n$, let $[n]$ denote the set $\\{1,2,\cdots,n\\}$. For
any function $f$, we define $\widetilde{O}(f)$ to be $f\cdot\log^{O(1)}(f)$.
In addition to $O(\cdot)$ notation, for two functions $f$, $g$, we use the
shorthand $f\lesssim g$ (resp. $\gtrsim$) to indicate that $f\leq C\cdot g$
(resp. $\geq$) for some absolute constant $C.$

For vector $x\in\mathbb{R}^{n}$, we use $\|x\|_{2}$ to denote its $\ell_{2}$
norm, we use $\|x\|_{1}$ to denote its $\ell_{1}$ norm, we use
$\|x\|_{\infty}$ to denote its $\ell_{\infty}$ norm. For matrix
$A\in\mathbb{R}^{m\times n}$, we use $\|A\|$ to denote the spectral norm of
$A$, we use $\|A\|_{F}$ to denote its Frobenius norm
($\|A\|_{F}=(\sum_{i=1}^{m}\sum_{j=1}^{n}A_{i,j}^{2})^{1/2}$), we use
$\|A\|_{1}$ to denotes its entry-wise $\ell_{1}$ norm
$\|A\|_{1}=\sum_{i=1}^{m}\sum_{j=1}^{n}|A_{i,j}|$. We use $A^{\dagger}$ to
denote its Moore–Penrose inverse. We use $A^{\top}$ denote the transpose of
$A$.

For square matrix $A$, we use $\operatorname{tr}[A]$ to denote the trace of
$A$. We say $A$ is positive-semidefinite (psd) if $A=A^{\top}$ and
$x^{\top}Ax\geq 0$, $\forall x\in\mathbb{R}^{n}$. We use $\succeq,\preceq$ to
denote the semidefinite ordering, e.g. $A\succeq 0$ means that $A$ is psd.

For square full-rank matrix $A$, we use $A^{-1}$ to denote the inverse of $A$.

For vectors $a,b\in\mathbb{R}^{n}$ and accuracy parameter $\epsilon\in(0,1)$,
we use $a\approx_{\epsilon}b$ to denote that $(1-\epsilon)b_{i}\leq
a_{i}\leq(1+\epsilon)b_{i}$, $\forall i\in[n]$. Similarly, for any scalar $t$,
we use $a\approx_{\epsilon}t$ to denote that $(1-\epsilon)t\leq
a_{i}\leq(1+\epsilon)t$, $\forall i\in[n]$.

###  3.2 Operators

######  Definition 3.1 (Operators $\tau$ and $\sigma$).

Given a matrix $A\in\mathbb{R}^{m\times n}$ and a vector $v\in\mathbb{R}^{m}$,
we define the leverage scores $\sigma(v)_{i}$ and (unnormalized) leverage
scores $\tau(v)_{i}$ as follows

$\displaystyle\tau(v)_{i}=(A(A^{\top}VA)^{-1}A^{\top})_{i,i}\mathrm{\leavevmode\nobreak\
and\leavevmode\nobreak\
}\sigma(v)_{i}=(\sqrt{V}A(A^{\top}VA)^{-1}A^{\top}\sqrt{V})_{i,i}.$

######  Definition 3.2 (Operators $M$, $Q$ and $P$).

Given a matrix $A\in\mathbb{R}^{m\times n}$, we define three operators :
$M:\mathbb{R}^{m}\rightarrow\mathbb{R}^{n\times n}$,
$Q:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m\times m}$ and
$P:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m\times m}$ such that for any vector
$v\in\mathbb{R}^{m}$

$\displaystyle M(v)=$ $\displaystyle\leavevmode\nobreak\ A^{\top}VA$
$\displaystyle Q(v)=$ $\displaystyle\leavevmode\nobreak\
A(A^{\top}VA)^{-1}A^{\top}$ $\displaystyle P(v)=$
$\displaystyle\leavevmode\nobreak\
\sqrt{V}A(A^{\top}VA)^{-1}A^{\top}\sqrt{V}.$

###  3.3 Different types of running time

######  Definition 3.3 (Pseudo-polynomial time).

An algorithm runs in pseudo-polynomial time if its running time is a
polynomial in the length of the input (the number of bits required to
represent it) and the numeric value of the input (the largest integer present
in the input).

######  Definition 3.4 (Strongly-polynomial time).

An algorithm runs in strongly polynomial time if (1) the number of operations
in the arithmetic model of computation is bounded by a polynomial in the
number of rational numbers in the input instance, and (2) the space used by
the algorithm is bounded by a polynomial in the size of the input.

######  Definition 3.5 (Weakly-polynomial time).

An algorithm runs in weakly-polynomial time if its running time is upper
bounded by a polynomial in the size of the input, and it doesn’t run in
strongly polynomial time.

###  3.4 Basic results on matrices

######  Fact 3.6 (Woodbury matrix identity, [Woo49, Woo50]).

Given a square invertible $n\times n$ matrix $A$, an $n\times k$ matrix $U$
and a $k\times n$ matrix $V$, let $B$ be an $n\times n$ matrix such that
$B=A+UCV$. Then, assuming $(I_{k}+VA^{-1}U)$ is invertible, we have

$\displaystyle B^{-1}=A^{-1}-A^{-1}U(C^{-1}+VA^{-1}U)^{-1}VA^{-1}.$

###  3.5 Fast matrix multiplication

######  Definition 3.7 (Matrix multiplication time).

For any $n,r>0$, denote $\mathcal{T}_{\operatorname{mat}}(n,n,r)$ the time to
compute the multiplication of an $n\times n$ matrix and an $n\times r$ matrix.

######  Theorem 3.8 (Fast matrix multiplication).

We have the following upper bound on
$\mathcal{T}_{\operatorname{mat}}(n,n,r)$:

  1. 1.

[GU18]. For $r=n^{0.31}$, we have
$\mathcal{T}_{\operatorname{mat}}(n,n,r)=O(n^{2+o(1)})$.

  2. 2.

[Cop82]. For $r=n^{0.17}$, we have
$\mathcal{T}_{\operatorname{mat}}(n,n,r)=O(n^{2}\log^{2}n)$.

  3. 3.

[BD76]. For $r=\log^{c}n$ for any constant $c>0$, then we have
$\mathcal{T}_{\operatorname{mat}}(n,n,r)=O(n^{2})$.

###  3.6 Multiple variable polynomial interpolation

We first state a one-variable version interpolation theorem.

######  Theorem 3.9 (One variable [IK94, Eq 10 in Page 331]).

Given any $2N$ times differentiable function
$f:\mathbb{R}\rightarrow\mathbb{R}$. There exits $N$ points
$s_{1},\cdots,s_{N}$ and $N$ weights $\omega_{1},\cdots,\omega_{N}\geq 0$ such
that $\sum_{i}\omega_{i}=1$ and

$\displaystyle\left|\int_{0}^{1}f(t)\mathrm{d}t-\sum_{i=1}^{N}\omega_{i}f(s_{i})\right|\leq
O(1)\cdot\frac{M_{2N}}{(2N)!4^{N}}$

where

$\displaystyle M_{2N}=\max_{t\in[0,1]}\left|\frac{\partial^{(2N)}f}{\partial
t^{(2N)}}(t)\right|.$

Now, we explain how to use one variable interpolation result (Theorem 3.9) to
prove a multiple variable interpolation result (Theorem 3.10).

######  Theorem 3.10 ($d$-variable case).

Given a function $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ such that for any
$j\in[d]$, for any $t_{1},\cdots,t_{j-1},t_{j+1},\cdots t_{d}\in[0,1]$,
$f(t_{1},\cdots,t_{j-1},t_{j},t_{j+1},\cdots,t_{d})$ is $N+1$ times
differentiable respect to $t_{j}$. Then there exists $N$ points
$s_{1},\cdots,s_{N}$ and $N$ weights $\omega_{1},\cdots,\omega_{N}$ such that

$\displaystyle\left|\int_{[0,1]^{d}}f(t)\mathrm{d}t-\sum_{i_{1}=1}^{N}\cdots\sum_{i_{d}=1}^{N}\omega_{i_{1}}\cdots\omega_{i_{d}}f(s_{i_{1}},\cdots,s_{i_{d}})\right|\leq\frac{1}{(2N)!2^{2N}}\cdot\sum_{j=1}^{d}\max_{t\in[0,1]^{d}}\left|\frac{\partial^{(2N)}f}{\partial
t_{j}^{(2N)}}(t)\right|.$

###### Proof.

Let $\tau$ denote

$\displaystyle\tau=\frac{1}{(2N)!2^{2N}}$

Using one variable theorem 3.9, we have that : for all
$t_{1},t_{2},\cdots,t_{d-1}\in[0,1]$,

$\displaystyle\left|\int_{0}^{1}f(t_{[d-1]},t_{d})\mathrm{d}t_{d}-\sum_{i_{d}=1}^{N}\omega_{i_{d}}f(t_{[d-1]},s_{i_{d}})\right|\leq\tau\cdot\max_{t_{[d-1]}\in[0,1]^{d-1}}\max_{t_{d}\in[0,1]}|f_{t_{d}}^{(2N)}(t)|$

To write recursive thing in an easy way, we let $g_{d}$ denote function $f$.
We use $t_{[j]}$ to denote $t_{1},t_{2},\cdots,t_{j}$.

For each $j\in\\{d-1,\cdots,1\\}$, we define function
$g_{j}:\mathbb{R}^{j}\rightarrow\mathbb{R}$ such that

$\displaystyle
g_{j}(t_{[j]})=\sum_{i_{j+1}=1}^{N}\omega_{i_{j+1}}g_{j+1}(t_{[j]},s_{i_{j+1}})$

Let $g_{0}=\sum_{i_{1}=1}^{N}\omega_{i_{1}}g_{1}(s_{i_{1}})$, we can also
rewrite $g_{0}$ as $\int_{[0,1]^{0}}g_{0}(t_{[0]})\mathrm{d}t_{[0]}$.

Finally, we want to bound

$\displaystyle\leavevmode\nobreak\
\left|\int_{0}^{1}\cdots\int_{0}^{1}f(t_{1},\cdots,t_{d})\mathrm{d}t_{1}\cdots\mathrm{d}t_{d}-\sum_{i_{1}=0}^{N}\cdots\sum_{i_{d}=1}^{N}\omega_{s_{i}}f(s_{i_{1}},\cdots,s_{i_{d}})\right|$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
\sum_{j=1}^{d}\left|\int_{[0,1]^{j}}g_{j}(t_{[j]})\mathrm{d}t_{[j]}-\int_{[0,1]^{j-1}}g_{j-1}(t_{[j-1]})\mathrm{d}t_{[j-1]}\right|$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
\tau\sum_{j=1}^{d}\max_{t\in[0,1]^{d}}\left|\frac{\partial^{(2N)}f}{\partial
t_{j}^{(2N)}(t)}\right|,$

where the last step follows by Claim 3.11. ∎

######  Claim 3.11 (Bounding the difference between $j$-th term and $j-1$-th
term).

For each $j\in[d]$,

$\displaystyle\left|\int_{[0,1]^{j}}g_{j}(t_{[j]})\mathrm{d}t_{[j]}-\int_{[0,1]^{j-1}}g_{j-1}(t_{[j-1]})\mathrm{d}t_{[j-1]}\right|\leq\tau\cdot\max_{t\in[0,1]^{d}}\left|\frac{\partial^{(2N)}f}{\partial
t_{j}^{(2N)}(t)}\right|.$

###### Proof.

First, using one variable theorem 3.9, we can upper bound

$\displaystyle\leavevmode\nobreak\
\left|g_{j}(t_{[j-1]},t_{j})\mathrm{d}t_{j}-\sum_{i_{j}=1}^{N}\omega_{i_{j}}g_{j}(t_{[j-1]},s_{i_{j}})\right|$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
\tau\cdot\max_{t_{[j-1]}\in[0,1]^{j-1}}\max_{t_{j}\in[0,1]}|g_{j,t_{j}}^{(N+1)}(t_{[j]})|$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\ \cdots$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
\tau\cdot\max_{t_{[j-1]}\in[0,1]^{j-1}}\max_{t_{j}\in[0,1]}\max_{t_{[d]\backslash[j]}\in[0,1]^{d-j}}|f^{(N+1)}_{t_{j}}(t_{[j-1]},t_{j},t_{[d]\backslash[j]})|$

where the last step follows by re-using $\omega_{i}>0$ and
$\sum_{i=1}^{N}\omega_{i}=1$.

We have

$\displaystyle\leavevmode\nobreak\
\left|\int_{[0,1]^{j}}g_{j}(t_{[j]})\mathrm{d}t_{[j]}-\int_{[0,1]^{j-1}}g_{j-1}(t_{[j-1]})\mathrm{d}t_{[j-1]}\right|$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\left|\int_{[0,1]^{j}}g_{j}(t_{[j]})\mathrm{d}t_{[j]}-\int_{[0,1]^{j-1}}\sum_{i_{j}=1}^{N}\omega_{i_{j}}g_{j}(t_{[j-1]},s_{i_{j}})\mathrm{d}t_{[j-1]}\right|$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\left|\int_{[0,1]^{j-1}}\left(g_{j}(t_{[j-1]},t_{j})\mathrm{d}t_{j}-\sum_{i_{j}=1}^{N}\omega_{i_{j}}g_{j}(t_{[j-1]},s_{i_{j}})\right)\mathrm{d}t_{[j-1]}\right|$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
\int_{[0,1]^{j-1}}\left|g_{j}(t_{[j-1]},t_{j})\mathrm{d}t_{j}-\sum_{i_{j}=1}^{N}\omega_{i_{j}}g_{j}(t_{[j-1]},s_{i_{j}})\right|\mathrm{d}t_{[j-1]}$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
\tau\cdot\max_{t_{[j-1]}\in[0,1]^{j-1}}\max_{t_{j}\in[0,1]}\max_{t_{[d]\backslash[j]}\in[0,1]^{d-j}}|f^{(2N)}_{t_{j}}(t_{[j-1]},t_{j},t_{[d]\backslash[j]})|$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\tau\cdot\max_{t\in[0,1]^{d}}\left|\frac{\partial^{(2N)}f}{\partial
t_{j}^{(2N)}(t)}\right|.$

∎

##  4 Perturbed Volumetric Center Cutting Plane Method

Recall that the feasibility problem is defined as follows.

  * Feasibility Problem: Given a separation oracle for a set $K$ contained in a box of radius $R$ either find a point $x\in K$ or prove that $K$ does not contain a ball of radius $\epsilon$.

Let $K\subset\mathbb{R}^{n}$ be a convex set that is contained in a box of
radius $R$ and equipped with a separation oracle. To solve the feasibility
problem, a cutting plane method iteratively refines a feasible region $\Omega$
and a query point $z$. In each iteration, $z$ is the input to the separation
oracle which then either certifies $z\in K$ or produces a hyperplane
$a^{\top}x=b$ separating $z$ from $K$. The algorithm then refines $\Omega$ and
computes the next query point $z^{\text{(new)}}$ from the new separating
hyperplane and any past information.

At a high level, Vaidya’s method heavily employs _leverage scores_ to decide
which hyperplane to add or to drop from the feasible region. The main
innovation behind our faster implementation is a data structure that maintains
an estimate of the leverage scores in amortized $O(n^{2})$. This removes the
bottleneck in Vaidya’s algorithm and yields the following result.

######  Theorem 4.1.

Given a separation oracle for a convex set $K\subset\mathbb{R}^{n}$ that is
contained in a box of radius $R$ and a parameter $\epsilon>0$, there is a
cutting plane method that either computes a point in $K$ or proves that $K$
does not contain a ball of radius $\epsilon$ in
$O((n\operatorname{SO}+n^{3})\log(\kappa))$ time, where $\operatorname{SO}$ is
the complexity of the separation oracle and $\kappa=nR/\epsilon$.

Since the proof of the validity of Vaidya’s method for our result is mostly a
perturbed version of his analysis, we defer the proof of Theorem 4.1 to
Section A.

##  5 Main Data Structure for Leverage Score Maintenance

In this section, we present our main data structure for leverage score
maintenance that achieves an amortized $O(n^{2})$ time per update. Combined
with Theorem 4.1 in Section 4, our main data structure implies a faster
$O(n\operatorname{SO}\log(\kappa)+n^{3}\log(\kappa))$ time cutting plane
method. Our main data structure further uses the simple leverage score
maintenance data structure in Section 7 and the complicated leverage score
maintenance data structure in Section 8.

######  Theorem 5.1 (Leverage score maintenance).

Given an initial matrix $A\in\mathbb{R}^{m\times n}$ with $m=O(n)$, initial
weight $w\in\mathbb{R}_{+}^{m}$. Then for any constant $c>0$, there is a
randomized data structure (Algorithm 1) that approximately maintains the
leverage scores

$\displaystyle\sigma_{i}(w)=(\sqrt{W}A(A^{\top}WA)^{-1}A^{\top}\sqrt{W})_{i,i},$

for positive diagonal matrices $W=\operatorname{diag}(w)$ through the
following operations:

  1. 1.

$\textsc{Init}(A,w)$: takes $O(n^{\omega+o(1)})$ time to intialize the data
structure.

  2. 2.

$\textsc{Update}(\mathsf{act})$: updates the data structure for the single
update $\mathsf{act}$.

  3. 3.

$\textsc{Query}()$: takes $O(n)$ time to output a vector
$\widetilde{\sigma}\in\mathbb{R}^{m}$.

Moreover, if the sequence of $K$ updates $(A^{(k)}\in\mathbb{R}^{m^{(k)}\times
n},w^{(k)}\in\mathbb{R}_{+}^{m^{(k)}})$ with $m^{(k)}=O(n)$ satisfies that
each update $k$ is one of the following:

  1. 1.

Insertion (resp. deletion) of row $a$ with weight $w_{a}$ into (resp. from)
$(A^{(k-1)},w^{(k-1)})$ that satisfies

$\displaystyle w_{a}aa^{\top}\preceq
0.01(A^{(k-1)})^{\top}W^{(k-1)}A^{(k-1)}.$

  2. 2.

Update $w^{(k-1)}\in\mathbb{R}_{+}^{m^{(k-1)}}$ to
$w^{(k)}\in\mathbb{R}_{+}^{m^{(k)}}$ such that $m^{(k-1)}=m^{(k)}$ and that

$\displaystyle\left\lVert\log(w^{(k)})-\log(w^{(k-1)})\right\rVert_{2}\leq
0.01,$

then the function $\textsc{Update}(\mathsf{act})$ takes an amortized
$O(n^{2})$ time, and the vector $\widetilde{\sigma}^{(k)}$ output by
$\textsc{Query}()$ at each step $k\in[K]$ satisfies that

$\displaystyle\left\lVert\widetilde{\sigma}^{(k)}-\sigma(w^{(k)})\right\rVert_{2}\leq
O(1/\log^{c}(n)).$ (1)

The proof of Theorem 5.1 relies on the data structures in Sections 6, 7 and 8.

###### Proof of Theorem 5.1.

We first prove the running time upper bound of the different functions. Notice
that the running time upper bound for $\textsc{Init}(A,w)$, and
$\textsc{Query}()$ are immediate corollaries of Theorem 7.2 and 8.1. We prove
running time upper bound for the function $\textsc{Update}(\mathsf{act})$ in
the following. In particular, we show that the inner, middle and outer phases
all run in amortized $O(n^{2})$ time, and the amortized time to restart the
data structure in Step 29 is $O(n^{2})$.

We start by analyzing the running time of the inner phase. Since each call to
$\textsc{Update}(\mathsf{act})$ makes one call of
$\textsc{Update}(\mathsf{act})$ with $\mathsf{act}$ having a single action
$\mathsf{act}$, it follows from Theorem 7.2 that the inner phase makes one
call to $\textsc{pm}.\textsc{update}$ and takes extra time at most $O(n^{2})$.
Since we choose our parameter as
$\epsilon_{\operatorname{inn}}=1/\log^{25}(n)$ (see Table 5), it then follows
from Theorem B.4 with $C=O(1)$ that the amortized time per call to
$\textsc{pm}.\textsc{update}$ is $O(n^{2})$. Therefore, the inner phase runs
in amortized time $O(n^{2})$.

Next we analyze the running time of the middle phase. Notice that each middle
update happens once every $T_{\operatorname{inn}}=\log^{10}(n)\leq n^{0.08}$
calls to $\textsc{Update}(\mathsf{act})$, so Assumption 7.1 holds for the
middle update in Step 23 of Algorithm 1. It then follows from Theorem 7.2 that
for every $T_{\operatorname{inn}}$ calls to the function
$\textsc{Update}(\mathsf{act})$, the middle phase makes one call to
$\textsc{pm}.\textsc{update}$ and takes extra time at most
$O(n^{2}T_{\operatorname{inn}})$. Since we choose
$\epsilon_{\mathrm{mid}}=1/n^{0.08}$, it follows from Theorem B.2 with
$C=T_{\operatorname{inn}}=\log^{10}(n)$ that the time for one call to
$\textsc{pm}.\textsc{update}$ is $O(n^{2}\log^{2}(n))$. Amortized over the
$T_{\operatorname{inn}}=\log^{10}(n)$ calls to the function
$\textsc{Update}(\mathsf{act})$, the amortized running time for the middle
phase is thus $O(n^{2})$.

Finally we analyze the running time of the outer phase. Each outer update
happens once every $T_{\operatorname{inn}}\cdot
T_{\mathrm{mid}}=n^{0.01}\log^{10}(n)$ calls to
$\textsc{Update}(\mathsf{act})$, so Assumption 7.1 holds for the outer update
in Step 26 of Algorithm 1. It then follows from Theorem 8.1 that for every
$T_{\operatorname{inn}}\cdot T_{\mathrm{mid}}$ calls to the function
$\textsc{Update}(\mathsf{act})$, the outer phase makes $O(N)=O(\log^{2}(n))$
calls to $\textsc{pm}.\textsc{update}$ and takes extra time at most

$\displaystyle O(n^{2}\cdot T_{\operatorname{inn}}\cdot
T_{\mathrm{mid}}+T(n,n,r_{\operatorname{out}})\cdot
N^{3}\cdot\log(n))=O(n^{2}\cdot T_{\operatorname{inn}}\cdot
T_{\mathrm{mid}}+n^{2+o(1)}\log^{7}(n)).$

Since we choose $\epsilon_{\operatorname{out}}=1/n^{0.1}$, it follows from
Theorem B.2 with $C=T_{\operatorname{inn}}\cdot T_{\mathrm{mid}}$ that the
time for $O(N)$ calls to $\textsc{pm}.\textsc{update}$ is $O(n^{2+o(1)})$.
Amortized over the $T_{\operatorname{inn}}\cdot
T_{\mathrm{mid}}=n^{0.01}\log^{10}(n)$ calls to the function
$\textsc{Update}(\mathsf{act})$, the amortized running time for the outer
phase is again $O(n^{2})$.

It remains to upper bound the amortized time to restart the data structure due
to Step 29. We note that restarting takes time $O(n^{\omega+o(1)})$ every
$T_{\operatorname{inn}}\cdot T_{\mathrm{mid}}\cdot
T_{\operatorname{out}}=\Omega(n^{\omega-1.99})$ calls to the function
$\textsc{Update}(\mathsf{act})$. It follows that the amortized time to restart
the data structure is $O(n^{\omega+o(1)})/(T_{\operatorname{inn}}\cdot
T_{\mathrm{mid}}\cdot T_{\operatorname{out}})=O(n^{2})$. This finishes the
running time analysis.

Now we proceed to prove the error guarantee given in the second part of the
theorem. We only prove (1) for $c=15$ as it’s easy to adjust the parameter
settings in Table 5 to achieve the guarantee for any constant $c>0$. By
Theorem 7.2, the $\ell_{2}$-error of the inner phase is
$O(\epsilon_{\operatorname{inn}})=O(\log^{-25}(n))$ per call to
$\textsc{Update}(\mathsf{act})$ and the inner phase runs for at most
$T_{\operatorname{inn}}=\log^{10}(n)$ steps before its estimate gets refined
by the estimate of the middle or outer phases (Step 37 and 33). Therefore, the
error introduced by the inner phase is at most
$O(\epsilon_{\operatorname{inn}})\cdot
T_{\operatorname{inn}}=O(\log^{-15}(n))$. For the middle phase, it follows
from Theorem 7.2 that the $\ell_{2}$-error introduced by each middle step is
at most $O((\epsilon_{\mathrm{mid}}+T_{\operatorname{inn}}\cdot
n^{-0.08})\cdot T_{\operatorname{inn}})\leq n^{-0.07}$. Since the middle phase
run for at most $T_{\mathrm{mid}}=n^{0.01}$ steps before its estimate gets
refined by the estimate of the outer phase (Step 34), the error of the middle
phase is at most $n^{-0.07}\cdot n^{0.01}=n^{-0.06}$. For the outer phase, the
error in Theorem 8.1 is negligible so it suffices to bound the variance. From
Theorem 8.1, the variance in each outer step is
$O(N^{3}/r_{\operatorname{out}})\cdot\epsilon_{\operatorname{out}}^{2}\cdot(T_{\operatorname{inn}}\cdot
T_{\mathrm{mid}})^{2}\leq O(n^{-0.48})$. Since the outer phase run for at most
$T_{\operatorname{out}}=n^{\omega-2}$ steps before the data structure gets
restarted, the variance accumulated in the $T_{\operatorname{out}}$ steps is
at most $O(n^{-0.1})$. Therefore, we have
$\left\lVert\widetilde{\sigma}^{(k)}-\sigma(w^{(k)})\right\rVert_{2}\leq
O(\log^{-15}(n))$ for each $k\in[K]$. This completes the proof of Theorem 5.1.

∎

Parameter | Value | Description
---|---|---
$T_{\operatorname{inn}}$ | $\log^{10}n$ | Number of inner steps per middle step
$T_{\mathrm{mid}}$ | $n^{0.01}$ | Number of middle step per outer step
$T_{\operatorname{out}}$ | $n^{\omega-2}$ | Number of outer step before re-computing
$r_{\operatorname{out}}$ | $n^{0.31}$ | Matrix multiplication dimension for outer step
$\epsilon_{\operatorname{inn}}$ | $1/\log^{25}n$ | Inner step projection maintenance error
$\epsilon_{\mathrm{mid}}$ | $1/n^{0.08}$ | Middle step projection maintenance error
$\epsilon_{\operatorname{out}}$ | $1/n^{0.1}$ | Outer step projection maintenance error
$N$ | $100\log^{2}n$ | Number of discrete sampling points in an outer step
Table 5: Values of parameters used in the main data structure for leverage
score maintenance in Algorithm 1. Algorithm 1

1:data structure MaintainLeverageScoreMain $\triangleright$ Theorem 5.1

2:

3:members

4:
$\operatorname{ctr}_{\operatorname{inn}},\operatorname{ctr}_{\mathrm{mid}},\operatorname{ctr}_{\operatorname{out}}\in\mathbb{Z}_{+}$
$\triangleright$ Step counters for different phases

5:
$\mathsf{acts}_{\mathrm{mid}},\mathsf{acts}_{\operatorname{out}}\in\textsc{Vector}\langle\textsc{Action}\rangle$
$\triangleright$ Delayed sequence of updates

6: $\triangleright$ Action includes insertion, deletion and $w$-update

7: MaintainLeverageScoreSimple $\textsc{simp}_{\operatorname{inn}}$,
$\textsc{simp}_{\mathrm{mid}}$ $\triangleright$ For inner and middle phases,
Theorem 7.2, Algorithm 4

8: MaintainLeverageScoreComplicated $\textsc{comp}_{\operatorname{out}}$
$\triangleright$ For outer phase, Theorem 8.1, Algorithm 5

9:end members

10:

11:procedure Init($A\in\mathbb{R}^{m\times n},w\in\mathbb{R}_{+}^{m}$)
$\triangleright$ Initialization

12:
$\operatorname{ctr}_{\operatorname{inn}},\operatorname{ctr}_{\mathrm{mid}},\operatorname{ctr}_{\operatorname{out}}\leftarrow
0$ $\triangleright$ Reset counters

13:
$\mathsf{acts}_{\mathrm{mid}},\mathsf{acts}_{\operatorname{out}}\leftarrow\textsc{NULL}$

14:
$\textsc{simp}_{\operatorname{inn}}$.Init($A,w,\epsilon_{\operatorname{inn}}$)

15: $\textsc{simp}_{\mathrm{mid}}$.Init($A,w,\epsilon_{\mathrm{mid}}$)

16:
$\textsc{comp}_{\operatorname{out}}$.Init($A,w,r_{\operatorname{out}},N,\epsilon_{\operatorname{out}}$)
$\triangleright$ Initialize member data structures

17:end procedure

18:

19:procedure Update($\mathsf{act}\in\textsc{Action}$) $\triangleright$ Update
estimate for a single update

20: Add $\mathsf{act}$ to $\mathsf{acts}_{\mathrm{mid}}$ and
$\mathsf{acts}_{\operatorname{out}}$

21: $\textsc{simp}_{\operatorname{inn}}$.Update($\mathsf{act}$),
$\operatorname{ctr}_{\operatorname{inn}}\leftarrow\operatorname{ctr}_{\operatorname{inn}}+1$
$\triangleright$ Inner update

22: if $\operatorname{ctr}_{\operatorname{inn}}=T_{\operatorname{inn}}$ then
$\triangleright$ Number of inner steps reaches limit

23: $\textsc{simp}_{\mathrm{mid}}$.Update($\mathsf{acts}_{\mathrm{mid}}$),
$\operatorname{ctr}_{\mathrm{mid}}\leftarrow\operatorname{ctr}_{\mathrm{mid}}+1$
$\triangleright$ Middle update

24: $\mathsf{acts}_{\mathrm{mid}}\leftarrow\textsc{NULL}$

25: if $\operatorname{ctr}_{\mathrm{mid}}=T_{\mathrm{mid}}$ then
$\triangleright$ Number of middle steps reaches limit

26:
$\textsc{comp}_{\operatorname{out}}$.Update($\mathsf{acts}_{\operatorname{out}}$),
$\operatorname{ctr}_{\operatorname{out}}\leftarrow\operatorname{ctr}_{\operatorname{out}}+1$
$\triangleright$ Outer update

27: $\mathsf{acts}_{\operatorname{out}}\leftarrow\textsc{NULL}$

28: if $\operatorname{ctr}_{\operatorname{out}}=T_{\operatorname{out}}$ then
$\triangleright$ Number of outer steps reaches limit

29: Self.Init($\textsc{simp}_{\operatorname{inn}}.A,w^{\operatorname{new}}$)
$\triangleright$ Restart the data structure

30: end procedure

31: end if

32: $\operatorname{ctr}_{\operatorname{inn}}\leftarrow
0,\operatorname{ctr}_{\mathrm{mid}}\leftarrow 0$ $\triangleright$ Restart
counters for inner/middle phases

33:
$\textsc{simp}_{\operatorname{inn}}.\textsc{RefineEstimate}(\textsc{comp}_{\operatorname{out}}.\textsc{Query}())$

34:
$\textsc{simp}_{\mathrm{mid}}.\textsc{RefineEstimate}(\textsc{comp}_{\operatorname{out}}.\textsc{Query}())$
$\triangleright$ Refine $\widetilde{\sigma}_{\operatorname{inn}}$ and
$\widetilde{\sigma}_{\mathrm{mid}}$

35: else

36: $\operatorname{ctr}_{\operatorname{inn}}\leftarrow 0$ $\triangleright$
Restart counter for inner phase

37:
$\textsc{simp}_{\operatorname{inn}}.\textsc{RefineEstimate}(\textsc{simp}_{\mathrm{mid}}.\textsc{Query}())$
$\triangleright$ Refine $\widetilde{\sigma}_{\operatorname{inn}}$

38: end if

39: end if

40:end procedure

41:

42:procedure Query()

43: return $\textsc{simp}_{\operatorname{inn}}$.Query() $\triangleright$
Return the estimate maintained by $\textsc{simp}_{\operatorname{inn}}$

44:end procedure

45:

46:end data structure

##  6 Batched Low Rank Update

###  6.1 Simple low rank update

In this section, we show how to maintain leverage score under low rank update.
First, we start with a preconditioning lemma showing that given two spectrally
similar matrix, we can invert one matrix faster by knowing the inverse of the
other matrix.

######  Lemma 6.1.

Given $n\times n$ PSD matrices $A$ and $M$ such that $0\preceq M\preceq
A\preceq\kappa\cdot M$. For any integer $t\geq 1$ such that
$\kappa(1-1/\kappa)^{t+1}<1$, we have

$\displaystyle f(M,t)\preceq
A^{-1}\preceq\frac{1}{1-\kappa(1-1/\kappa)^{t+1}}\cdot f(M,t)$

where $f:\mathbb{R}^{n\times n}\times\mathbb{N}\rightarrow\mathbb{R}^{n\times
n}$ is defined as

$\displaystyle
f(M,t)=\frac{1}{\kappa}M^{-1}\sum_{i=0}^{t}(I-\frac{1}{\kappa}AM^{-1})^{i}.$

Furthermore, if $M^{-1}$ and $A$ are given explicitly, for any
$V\in\mathbb{R}^{n\times r}$, then we can compute $f(M,t)\cdot V$ in
$O(\mathcal{T}_{\operatorname{mat}}(n,n,r)\cdot t)$ time.

###### Proof.

Note that

$\displaystyle A^{-1}=$ $\displaystyle\leavevmode\nobreak\
M^{-\frac{1}{2}}(M^{-\frac{1}{2}}AM^{-\frac{1}{2}})^{-1}M^{-\frac{1}{2}}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\frac{1}{\kappa}M^{-\frac{1}{2}}(I-(I-\frac{1}{\kappa}M^{-\frac{1}{2}}AM^{-\frac{1}{2}}))^{-1}M^{-\frac{1}{2}}.$
(2)

Using

$\displaystyle 0\preceq I-(1/\kappa)\cdot
M^{-\frac{1}{2}}AM^{-\frac{1}{2}}\preceq(1-1/\kappa)\cdot I$

and

$\displaystyle\sum_{i=0}^{t}x^{i}\leq(1-x)^{-1}\leq\sum_{i=0}^{t}x^{i}+\kappa(1-1/\kappa)^{t+1}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \forall
x\in[0,1-1/\kappa],$

we have

$\displaystyle\sum_{i=0}^{t}\Big{(}I-\frac{1}{\kappa}M^{-\frac{1}{2}}AM^{-\frac{1}{2}}\Big{)}^{i}$
$\displaystyle\preceq\Big{(}I-(I-\frac{1}{\kappa}M^{-\frac{1}{2}}AM^{-\frac{1}{2}})\Big{)}^{-1}$
$\displaystyle\preceq\sum_{i=0}^{N}\Big{(}I-\frac{1}{\kappa}M^{-\frac{1}{2}}AM^{-\frac{1}{2}}\Big{)}^{i}+\kappa\cdot(1-1/\kappa)^{t+1}I$
(3)

Multiplying $\frac{1}{\kappa}M^{-1}$ on both sides of the above equation with
a formula $A^{-1}$ (Eq. (6.1)), we get

$\displaystyle\frac{1}{\kappa}M^{-1}\sum_{i=0}^{t}\Big{(}I-\frac{1}{\kappa}M^{-\frac{1}{2}}AM^{-\frac{1}{2}}\Big{)}^{i}\preceq
A^{-1}\preceq\frac{1}{\kappa}M^{-1}\sum_{i=0}^{t}\Big{(}I-\frac{1}{\kappa}M^{-\frac{1}{2}}AM^{-\frac{1}{2}}\Big{)}^{i}+(1-1/\kappa)^{t+1}M^{-1}$

Using the definition of $f(M,t)$, we have that

$\displaystyle f(M,t)\preceq A^{-1}\preceq f(M,t)+(1-1/\kappa)^{t+1}M^{-1}.$

Using $M^{-1}\preceq\kappa A^{-1}$, we have

$\displaystyle f(M,t)\preceq A^{-1}\preceq
f(M,t)+\kappa\cdot(1-1/\kappa)^{t+1}A^{-1}.$

The result follows from some rearranging. ∎

Now, we show how to do leverage score update under monotone updates.

######  Lemma 6.2.

Given a matrix $A\in\mathbb{R}^{m\times n}$, non-negative weights
$w,w^{\operatorname{new}}\in\mathbb{R}^{m}$ such that
$w^{\operatorname{new}}\geq w$ with $m=O(n)$ and
$\|w^{\operatorname{new}}-w\|_{0}=k$. Let $\beta>1$ denote the parameter such
that $A^{\top}W^{\operatorname{new}}A\preceq\beta\cdot A^{\top}WA$. Given some
explicit matrix $U\in\mathbb{R}^{n\times n}$ and
$U^{\operatorname{new}}\in\mathbb{R}^{n\times n}$ such that

$\displaystyle U^{-1}\preceq$ $\displaystyle\leavevmode\nobreak\
A^{\top}WA\preceq(1+1/\log n)U^{-1},$
$\displaystyle(U^{\operatorname{new}})^{-1}\preceq$
$\displaystyle\leavevmode\nobreak\
A^{\top}W^{\operatorname{new}}A\preceq(1+1/\log
n)(U^{\operatorname{new}})^{-1}.$

Then, for any $\epsilon\in(0,1/2)$, we can output a vector
$c\in\mathbb{R}^{m}$ such that

$\displaystyle\|c-(\sigma(w^{\operatorname{new}})-\sigma(w))\|_{2}\leq\epsilon$

in time

$\displaystyle
O(\mathcal{T}_{\operatorname{mat}}(n,n,k)\cdot(1+\log_{\log(n)}(\beta
k/\epsilon))).$

The same statement holds for decreasing $w$, namely
$w^{\operatorname{new}}\leq w$, with
$A^{\top}W^{\operatorname{new}}A\succeq\beta^{-1}A^{\top}WA$.

###### Proof.

We note that the decreasing case follows from the increasing case by swapping
$w$ and $w^{\operatorname{new}}$. Hence, we focus on the increasing case.
Using Lemma 6.1, we can construct $U_{\epsilon}$ such that $U_{\epsilon}$ is
low-degree polynomial of $U$ and $A^{\top}WA$ and that

$U_{\epsilon}\preceq(A^{\top}WA)^{-1}\preceq(1+\epsilon)U_{\epsilon}.$ (4)

We define $U_{\epsilon}^{\operatorname{new}}$ similarly. We will show how to
use $U_{\epsilon}$ and $U_{\epsilon}^{\operatorname{new}}$ to approximate
$\sigma(w^{\operatorname{new}})-\sigma(w)$. First, we note that

$\displaystyle\sigma(w^{\operatorname{new}})-\sigma(w)=(w^{\operatorname{new}}-w)\cdot\tau(w^{\operatorname{new}})+w\cdot(\tau(w^{\operatorname{new}})-\tau(w)).$

We define $c_{i}$, $c_{1,i}$, $c_{2,i}$ as follows:

$\displaystyle c_{i}=$ $\displaystyle\leavevmode\nobreak\
\sigma(w^{\operatorname{new}})_{i}-\sigma(w)_{i}$ $\displaystyle=$
$\displaystyle\leavevmode\nobreak\
\underbrace{(w^{\operatorname{new}}-w)_{i}\cdot\tau(w^{\operatorname{new}})_{i}}_{c_{1,i}}+\underbrace{w_{i}\cdot(\tau(w^{\operatorname{new}})-\tau(w))_{i}}_{c_{2,i}}$

We estimate the right hand side by

$\widetilde{c}_{i}:=\underbrace{(w^{\operatorname{new}}-w)_{i}\cdot(AU_{\widetilde{\epsilon}}^{\operatorname{new}}A^{\top})_{i,i}}_{\widetilde{c}_{1,i}}+\underbrace{w_{i}\cdot(A((U_{\widetilde{\epsilon}}^{-1}+\Delta)^{-1}-U_{\widetilde{\epsilon}})A^{\top})_{i,i}}_{\widetilde{c}_{2,i}}$
(5)

where $\Delta=A^{\top}W^{\operatorname{new}}A-A^{\top}WA$ and
$\widetilde{\epsilon}=\frac{\epsilon}{3\beta\sqrt{k}}$.

##### Cost of computing $\widetilde{c}$

For the first term in (5). We note that $w^{\operatorname{new}}-w$ has only
$k$ non-zeros. Hence, we only need to compute
$AU_{\widetilde{\epsilon}}^{\operatorname{new}}A^{\top}$ on $k$ of the
diagonals. Let $A_{S}\in\mathbb{R}^{k\times n}$ where each rows of $A$ that
$w^{\operatorname{new}}\neq w$. Then, we can compute the $k$ of the diagonals
via the whole matrix

$\displaystyle
A_{S}U_{\widetilde{\epsilon}}^{\operatorname{new}}A_{S}^{\top}.$

Lemma 6.1 shows that $U_{\widetilde{\epsilon}}^{\operatorname{new}}$ has
$O(1+\frac{\log\widetilde{\epsilon}}{\log\log(n)})$ terms and that it takes

$\displaystyle
O\left(\mathcal{T}_{\operatorname{mat}}(n,n,k)\cdot\Big{(}1+\frac{\log\widetilde{\epsilon}}{\log\log(n)}\Big{)}\right)$

to compute $U_{\widetilde{\epsilon}}^{\operatorname{new}}A_{S}^{\top}$.
Finally, it takes extra
$\mathcal{T}_{\operatorname{mat}}(k,n,k)=O(\mathcal{T}_{\operatorname{mat}}(n,n,k))$
to do the left multiplication on $A_{S}$.

For the second term in (5). Woodbury matrix identity shows that

$\displaystyle(U_{\widetilde{\epsilon}}^{-1}+\Delta)^{-1}-U_{\widetilde{\epsilon}}=U_{\widetilde{\epsilon}}A_{S}^{\top}(\Delta_{W}^{-1}+A_{S}U_{\widetilde{\epsilon}}A_{S}^{\top})^{-1}A_{S}U_{\widetilde{\epsilon}}$

where $\Delta_{W}\in\mathbb{R}^{k\times k}$ be $W^{\operatorname{new}}-W$
restricted on non-zeros. By same argument above, we can compute
$A_{S}U_{\widetilde{\epsilon}}A_{S}^{\top}$ in

$\displaystyle
O\left(\mathcal{T}_{\operatorname{mat}}(n,n,k)\cdot\Big{(}1+\frac{\log\widetilde{\epsilon}}{\log\log(n)}\Big{)}\right)$

time.

Then, we can compute
$(\Delta_{W}^{-1}+A_{S}U_{\widetilde{\epsilon}}A_{S}^{\top})^{-1}$ in
$\mathcal{T}_{\operatorname{mat}}(k,k,k)=O(\mathcal{T}_{\operatorname{mat}}(n,n,k))$
time. Then, we compute

$\displaystyle
AU_{\widetilde{\epsilon}}A_{S}^{\top}(\Delta_{W}^{-1}+A_{S}U_{\widetilde{\epsilon}}A_{S}^{\top})^{-1}\qquad\text{in
time
}\mathcal{T}_{\operatorname{mat}}(m,n,k)=O(\mathcal{T}_{\operatorname{mat}}(n,n,k))$

and

$\displaystyle
A_{S}U_{\widetilde{\epsilon}}=(U_{\widetilde{\epsilon}}A_{S}^{\top})^{\top}\qquad\text{in
time
}O\left(\mathcal{T}_{\operatorname{mat}}(n,n,k)\cdot\Big{(}1+\frac{\log\widetilde{\epsilon}}{\log\log(n)}\Big{)}\right).$

Now, we multiple the two terms above together in time
$\mathcal{T}_{\operatorname{mat}}(n,k,n)=O(\mathcal{T}_{\operatorname{mat}}(n,n,k))$.
Hence, the total time is

$\displaystyle
O\left(\mathcal{T}_{\operatorname{mat}}(n,n,k)\cdot\Big{(}1+\frac{\log\widetilde{\epsilon}}{\log\log(n)}\Big{)}\right)=O(\mathcal{T}_{\operatorname{mat}}(n,n,k)(1+\log_{\log(n)}(\beta
k/\epsilon))).$

##### The accuracy of $\widetilde{c}$

For the first term in $\widetilde{c}$, we have

$\displaystyle A(A^{\top}W^{\operatorname{new}}A)^{-1}A^{\top}\preceq
AU_{\widetilde{\epsilon}}^{\operatorname{new}}A^{\top}\preceq(1+\widetilde{\epsilon})\cdot
A(A^{\top}W^{\operatorname{new}}A)^{-1}A^{\top}.$

Hence, we have

$\displaystyle|\widetilde{c}_{1,i}-c_{1,i}|=$
$\displaystyle\leavevmode\nobreak\
\left|(w^{\operatorname{new}}-w)_{i}\left(AU_{\widetilde{\epsilon}}^{\operatorname{new}}A^{\top}-A(A^{\top}W^{\operatorname{new}}A)^{-1}A^{\top}\right)_{i,i}\right|$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
\widetilde{\epsilon}\cdot|w^{\operatorname{new}}_{i}-w_{i}|\cdot(A(A^{\top}W^{\operatorname{new}}A)^{-1}A^{\top})_{i,i}$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
\widetilde{\epsilon}\cdot
w_{i}^{\operatorname{new}}\cdot(A(A^{\top}W^{\operatorname{new}}A)^{-1}A^{\top})_{i,i}$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\ \widetilde{\epsilon}$

where we used that $w^{\operatorname{new}}\geq w$ in the second inequality and
leverage score is upper bounded by $1$ in the last equality. Hence the
$\ell_{\infty}$ norm error of the first term is given by
$\widetilde{\epsilon}$. Hence, the $\ell_{2}$ norm error is given by

$\|\widetilde{c}_{1}-c_{1}\|_{2}\leq\sqrt{k}\cdot\widetilde{\epsilon}.$ (6)

For the second term in $\widetilde{c}$, we have that the error is given by

$\displaystyle\widetilde{c}_{2,i}-c_{2,i}=w_{i}\cdot(A((U_{\widetilde{\epsilon}}^{-1}+\Delta)^{-1}-U_{\widetilde{\epsilon}})A^{\top})_{i,i}-w_{i}\cdot(A((M+\Delta)^{-1}-M^{-1})A^{\top})_{i,i}$

where $M=A^{\top}WA$. To simplify the notation, we define
$M_{t,s}=M+t\cdot\Delta_{2}+s\cdot\Delta$ where
$\Delta_{2}=U_{\widetilde{\epsilon}}^{-1}-M$. Then, the error term becomes

$(\sqrt{W}A(M_{1,1}^{-1}-M_{1,0}^{-1}-M_{0,1}^{-1}+M_{0,0}^{-1})A^{\top}\sqrt{W})_{i,i}$
(7)

Note that

$\displaystyle M_{1,1}^{-1}-M_{1,0}^{-1}-M_{0,1}^{-1}+M_{0,0}^{-1}$
$\displaystyle=$
$\displaystyle\int_{0}^{1}\int_{0}^{1}\frac{\mathrm{d}^{2}}{\mathrm{d}t\mathrm{d}s}M_{t,s}^{-1}\mathrm{d}t\mathrm{d}s$
$\displaystyle=$
$\displaystyle\int_{0}^{1}\int_{0}^{1}M_{t,s}^{-1}\Delta_{2}M_{t,s}^{-1}\Delta
M_{t,s}^{-1}+M_{t,s}^{-1}\Delta
M_{t,s}^{-1}\Delta_{2}M_{t,s}^{-1}\mathrm{d}t\mathrm{d}s$ (8)

Combining (7) and (8), we have that the $\ell_{2}$ norm error for the second
term is bounded by

$\displaystyle\|\widetilde{c}_{2}-c_{2}\|_{2}\leq$
$\displaystyle\leavevmode\nobreak\
\left\|\int_{0}^{1}\int_{0}^{1}\sqrt{W}A(M_{t,s}^{-1}\Delta_{2}M_{t,s}^{-1}\Delta
M_{t,s}^{-1}+M_{t,s}^{-1}\Delta
M_{t,s}^{-1}\Delta_{2}M_{t,s}^{-1})A^{\top}\sqrt{W}\mathrm{d}t\mathrm{d}s\right\|_{F}$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
2\max_{t\in[0,1],s\in[0,1]}\left\|\sqrt{W}AM_{t,s}^{-1}\Delta_{2}M_{t,s}^{-1}\Delta
M_{t,s}^{-1}A^{\top}\sqrt{W}\right\|_{F}.$

To bound the Frobenius norm, we note that $\Delta\succeq 0$ by the assumption
and that $\Delta_{2}=U_{\widetilde{\epsilon}}^{-1}-M\succeq 0$ (4). Hence, we
have that $M_{t,s}\succeq M$ and

$\displaystyle\|\sqrt{W}AM_{t,s}^{-1}\Delta_{2}M_{t,s}^{-1}\Delta
M_{t,s}^{-1}A^{\top}\sqrt{W}\|_{F}^{2}$ $\displaystyle=$
$\displaystyle\operatorname{tr}\Big{[}A^{\top}WAM_{t,s}^{-1}\Delta_{2}M_{t,s}^{-1}\Delta
M_{t,s}^{-1}A^{\top}WAM_{t,s}^{-1}\Delta
M_{t,s}^{-1}\Delta_{2}M_{t,s}^{-1}\Big{]}$ $\displaystyle=$
$\displaystyle\operatorname{tr}\Big{[}M^{\frac{1}{2}}M_{t,s}^{-1}\Delta_{2}M_{t,s}^{-1}\Delta
M_{t,s}^{-1}MM_{t,s}^{-1}\Delta
M_{t,s}^{-1}\Delta_{2}M_{t,s}^{-1}M^{\frac{1}{2}}\Big{]}$ $\displaystyle\leq$
$\displaystyle\operatorname{tr}\Big{[}M^{\frac{1}{2}}M_{t,s}^{-1}\Delta_{2}M_{t,s}^{-1}\Delta
M_{t,s}^{-1}\Delta M_{t,s}^{-1}\Delta_{2}M_{t,s}^{-1}M^{\frac{1}{2}}\Big{]}$
$\displaystyle=$
$\displaystyle\operatorname{tr}\Big{[}M_{t,s}^{-\frac{1}{2}}\Delta
M_{t,s}^{-1}\Delta_{2}M_{t,s}^{-1}M^{\frac{1}{2}}M^{\frac{1}{2}}M_{t,s}^{-1}\Delta_{2}M_{t,s}^{-1}\Delta
M_{t,s}^{-\frac{1}{2}}\Big{]}$ $\displaystyle\leq$
$\displaystyle\operatorname{tr}\Big{[}M_{t,s}^{-\frac{1}{2}}\Delta
M_{t,s}^{-1}\Delta_{2}M^{-1}\Delta_{2}M_{t,s}^{-1}\Delta
M_{t,s}^{-\frac{1}{2}}\Big{]}.$

(4) shows that $\Delta_{2}\preceq\widetilde{\epsilon}\cdot M$ and hence
$(M^{-\frac{1}{2}}\Delta_{2}M^{-\frac{1}{2}})^{2}\preceq\widetilde{\epsilon}^{2}\cdot
I$. Continuing the last equation, we have

$\displaystyle\|\sqrt{W}AM_{t,s}^{-1}\Delta_{2}M_{t,s}^{-1}\Delta
M_{t,s}^{-1}A^{\top}\sqrt{W}\|_{F}^{2}$
$\displaystyle\leq\widetilde{\epsilon}^{2}\cdot\operatorname{tr}\Big{[}M_{t,s}^{-\frac{1}{2}}\Delta
M_{t,s}^{-1}MM_{t,s}^{-1}\Delta M_{t,s}^{-\frac{1}{2}}\Big{]}$
$\displaystyle\leq\widetilde{\epsilon}^{2}\cdot\operatorname{tr}\Big{[}M_{t,s}^{-\frac{1}{2}}\Delta
M_{t,s}^{-1}\Delta M_{t,s}^{-\frac{1}{2}}\Big{]}$
$\displaystyle\leq\widetilde{\epsilon}^{2}\cdot\|M^{-\frac{1}{2}}\Delta
M^{-\frac{1}{2}}\|_{F}^{2}.$

Finally, using $\Delta$ has rank $k$, we have

$\displaystyle\|\widetilde{c}_{2}-c_{2}\|_{2}\leq$
$\displaystyle\leavevmode\nobreak\
2\widetilde{\epsilon}\cdot\|M^{-\frac{1}{2}}\Delta M^{-\frac{1}{2}}\|_{F}$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
2\widetilde{\epsilon}\cdot\sqrt{k}\|M^{-\frac{1}{2}}\Delta M^{-\frac{1}{2}}\|$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
2\widetilde{\epsilon}\cdot\beta\sqrt{k}.$ (9)

Combining (6) and (9), the total $\ell_{2}$ error is bounded by
$2\widetilde{\epsilon}\beta\sqrt{k}+\sqrt{k}\cdot\widetilde{\epsilon}\leq
3\widetilde{\epsilon}\beta\sqrt{k}$. This explains the choice of
$\widetilde{\epsilon}$. ∎

###  6.2 Batched low rank update

When we are given a sequence of update and we need to compute the change of
leverage score, we cannot apply Lemma 6.2 directly for the following reasons:

  * •

The update in general involves both increasing some weight and decreasing some
weight. This can be fixed by splitting the update into positive update and
negative update. (See Step 1 in Algorithm 2)

  * •

To get error $\frac{1}{n^{O(1)}}$, Lemma 6.2 takes at least $n^{2}\log(n)$
time. Hence, it is too costly to do a rank 1 update, which can happen if we
have alternating positive and negative updates (since we can only apply Lemma
6.2 for a monotone change). This can be fixed by moving all positive/update
and insert in the front of the update sequence. (See Step 2 in Algorithm 2)

  * •

Our $\ell_{2}$ update are generally dense. This can be fixed by ignoring all
small multiplicative update. (See Line 44 and 49 of Step 3 in Algorithm 2)

  * •

In general, the matrices $A^{\top}WA$ can change by a factor of
$e^{\Theta(T)}$ after $T$ iterations. So, Lemma 6.2 takes at least
$\Omega(\mathcal{T}_{\operatorname{mat}}(n,n,k)\frac{T}{\log\log(n)})$ time.
This forces us to do rank $k$ changes where
$\mathcal{T}_{\operatorname{mat}}(n,n,k)\leq n^{2}\log\log(n)$. Hence, for
$\ell_{2}$ update, we need to handle coordinates with not-so-small
multiplicative updates separately. This can be fixed by cutting the update
sequence into phases of length $L$ and noticing there are not too many
coordinate that is changed by $e^{\Theta(L)}$ factor. (See the rest of Step 3
in Algorithm 2)

Given the above intuitions, Theorem 6.3 follows from multiple use of Lemma
6.2.

######  Theorem 6.3 (Batched Update of Leverage Score).

Given an initial matrix $A\in\mathbb{R}^{m\times n}$, initial weight
$w\in\mathbb{R}_{+}^{m}$ and a sequence of $T$ updates that belong to one of
the following forms:

  * •

Update: Update the weight to $w$

  * •

Insert: Insert a row $a\in\mathbb{R}^{n}$ to the matrix
$A\in\mathbb{R}^{m\times n}$ and weight $w_{a}$ to the vector $w$

  * •

Delete: Delete a row $a\in\mathbb{R}^{n}$ from the matrix
$A\in\mathbb{R}^{m\times n}$ and weight $w_{a}$ from the vector $w$

Let the sequence of matrices and vectors be
$A^{(0)},A^{(1)},\cdots,A^{(T)}\in\mathbb{R}^{m\times n}$ and
$w^{(0)},w^{(1)},\cdots,w^{(T)}\in\mathbb{R}^{m}$. We assume that $T\leq
n^{0.08}$, the maximum number of rows is bounded by $O(n)$, and that the
changes satisfy

  * •

Update: $\|\log w^{(k)}-\log w^{(k-1)}\|_{2}\leq 0.01$

  * •

Insert/Delete: $w_{a}aa^{\top}\preceq 0.01\cdot A^{(k-1)}W^{(k-1)}A^{(k-1)}$

Then, there is an algorithm (BatchedUpdate in Algorithm 2 and 3) that finds a
vector $v\in\mathbb{R}^{n}$ and a vector $c\in\mathbb{R}^{n}$ such that

$\displaystyle\|\log v-\log w^{(T)}\|_{\infty}\leq T\cdot n^{-0.08},\|\log
v-\log w^{(T)}\|_{2}\leq 0.01\cdot T$

and that

$\displaystyle\|c-(\sigma_{A^{(T)}}(v)-\sigma_{A^{(0)}}(w^{(0)}))\|_{2}\leq
1/n^{100}.$

in time $O(n^{2}T)$.

###### Proof.

##### Correctness

Note that the algorithm 2 involves reducing the sequence $\\{M^{(0,k)}\\}$ to
$\\{M^{(1,k)}\\}$ to $\\{M^{(2,k)}\\}$ and finally to $\\{M^{(3,k)}\\}$. All
steps maintains the matrix at $k=0$. The first two steps maintains the matrix
at the last step. For the last step, we only ignore the rows with
multiplicative changes less than $\eta=n^{-0.08}$. Since there are at most
$T/L=n^{0.08}/\log^{3}(n)$ phases, the total accumulated multiplicative
changes for one row is $Tn^{-0.08}$. Therefore, we have

$\displaystyle M^{(3,\textsc{End})}=(A^{(T)})^{\top}VA^{(T)}$

for some $v$ such that $\|\log v-\log w^{(T)}\|_{\infty}\leq Tn^{-0.08}$.
Furthermore, we output the good enough approximation of
$\sigma_{A^{(T)}}(v)-\sigma_{A^{(0)}}(w^{(0)})=\sigma(M^{(3,\textsc{End})})-\sigma(M^{(3,0)})$
where we abused the notation to use $\sigma(M)$ to denote the leverage score
of the rows insides $M$. This is simply because we split the difference into

$\displaystyle\sigma(M^{(3,\textsc{End})})-\sigma(M^{(3,0)})=\sum_{k=1}^{T}\sigma(M^{(3,k)})-\sigma(M^{(3,k-1)})$

and estimate each step with $\ell_{2}$ error $n^{-1000}$. There are at most
$8T$ many terms, so the total error is less than $n^{-100}$. (The number of
terms can increases by at most a factor of 8 due to step 1, step 2 and step
4.)

##### Ratio on $M$ within each phase

By the assumptions on insert/delete and update, we see that
$0.8M^{(0,k-1)}\preceq M^{(0,k)}\preceq 1.2M^{(0,k-1)}$ for all $k$. The first
step clearly maintains this relation. We claim that the second step also
maintains this relation. To see this, we let
$\Delta^{(1,k)}=M^{(1,k+1)}-M^{(1,k)}$. By the relation, we have that

$-0.2M^{(1,k)}\preceq\Delta^{(k)}\preceq 0.2M^{(1,k)}.$ (10)

Note that step 2 simply permutes $\Delta^{(1,k)}$, namely, for each $k$ there
is an unique $k^{\prime}$ such that $\Delta^{(1,k)}=\Delta^{(2,k^{\prime})}$.
Finally, we note that $M^{(1,k)}\preceq M^{(2,k^{\prime})}$ for all $k$
because $M^{(2,k^{\prime})}$ is simply $M^{(2,k^{\prime})}$ plus some PSD
matrix (which due to some insert/positive update move into the sum or some
delete/negative update removed from the sum). Hence, we still have the
relation (10).

For the step 3, we should not expect the relation (10). However, we claim that
within each phase, let $M^{(3,k)}$ and $M^{(3,\overline{k})}$ be two matrices
in the phase. Then, we have that

$2^{-O(L)}M^{(3,\overline{k})}\preceq M^{(3,k)}\preceq
2^{O(L)}M^{(3,\overline{k})}.$ (11)

The proof for this is same for the phase with increasing $M$ and the phase
with decreasing phase. Hence, we only discuss the first case. Say $k_{1}$ is
the first matrix in the phase and $k_{2}$ is the last matrix in the phase.
Note that the step $3$ does not change the matrix $M^{(3,k_{1})}$ and
$M^{(3,k_{2})}$ because we only swapping operation order within a phase. Hence
by (10), we have

$\displaystyle M^{(3,k_{2})}=M^{(2,k_{2})}\preceq
2^{O(L)}M^{(2,k_{1})}=2^{O(L)}M^{(3,k_{1})}.$

Since all matrix in the phase is sandwiched between $M^{(3,k_{1})}$ and
$M^{(3,k_{2})}$ (due to the monotonicity of the sequence of $M$). Hence, we
have (11).

##### Cost of Insert/Delete (Line 65)

Now, we bound the runtime. We will show the each phase takes $O(n^{2}L)$ time.
Since there are $O(T/L)$ phases, this shows that the total cost is
$O(n^{2}T)$.

Each phase contains at most $L$ insert or delete. Hence, Line 65 involves
estimating leverage score up to rank $L$ update. By (11), we know the matrix
is changed by a $2^{O(L)}$ factor. Theorem 6.3 shows the cost is

$\displaystyle
O(\mathcal{T}_{\operatorname{mat}}(n,n,k)(1+\log_{\log(n)}(\beta
k/\epsilon)))$

where $\epsilon=n^{-1000}$, $\beta=2^{O(L)}$ and $k=O(L)$. By the choice of
$L=\log^{3}(n)$, we have $\mathcal{T}_{\operatorname{mat}}(n,n,k)=n^{2}$ and
hence the cost is $O(n^{2}L)$.

##### Cost of Large Update (Line 63)

Each phase contains at most $L$ updates. By the assumption, each update change
the $\ell_{2}$ norm of $\log w$ by at most $O(1)$. Note that this is maintains
during all steps. Hence, after merging all updates into $1$, we have $\|\log
w^{(k)}-\log w^{(k-1)}\|_{2}=O(L)$. By the construction of
$w^{(k-\frac{1}{2})}$, we have that $\|\log w^{(k)}-\log
w^{(k-\frac{1}{2})}\|_{2}=O(L)$ and that $|\log w_{i}^{(k-\frac{1}{2})}-\log
w_{i}^{(k)}|$ is either $0$ or at least $\log 2$ for all $i$. Hence the number
of coordinates that is non-zero in $\log w^{(k)}-\log w^{(k-\frac{1}{2})}$ is
at most $O(L^{2})$. Again, by (11), we know the matrix is changed by a
$2^{O(L)}$ factor. Hence, Theorem 6.3 shows the cost is

$\displaystyle
O(\mathcal{T}_{\operatorname{mat}}(n,n,k)(1+\log_{\log(n)}(\beta
k/\epsilon)))$

where $\epsilon=n^{-1000}$, $\beta=2^{O(L)}$ and $k=O(L^{2})$. By the choice
of $L=\log^{3}(n)$, we have $\mathcal{T}_{\operatorname{mat}}(n,n,k)=n^{2}$
and hence the cost is $O(n^{2}L)$.

##### Cost of Small Update (Line 62)

Similar to above, we have that $\|\log w^{(k-\frac{1}{2})}-\log
w^{(k-1)}\|_{2}=O(L)$. Furthermore, we have that $\|\log
w^{(k-\frac{1}{2})}-\log w^{(k-1)}\|_{\infty}=O(1)$. Due to Line 44 and 49, we
removed all rows with multiplicative changes less than $\eta$. Hence, the
number of coordinates that is non-zero in $\log w^{(k-\frac{1}{2})}-\log
w^{(k-1)}$ is at most $O(L^{2}/\eta^{2})$. Hence, Theorem 6.3 shows the cost
is

$\displaystyle
O(\mathcal{T}_{\operatorname{mat}}(n,n,k)(1+\log_{\log(n)}(\beta
k/\epsilon)))$

where $\epsilon=n^{-1000}$, $\beta=O(1)$ and $k=O(L^{2}/\eta^{2})$. By the
choice of $L=\log^{3}(n)$ and $\eta=n^{-0.08}$, we have
$\mathcal{T}_{\operatorname{mat}}(n,n,k)=n^{2}\log^{2}(n)$ and hence the cost
is $O(n^{2}\log^{3}(n))=O(n^{2}L)$.

##### Summary

For the cost of the rest of the algorithm, we note that step 1, 2 and 3 simply
involves rearrangements and it takes $O(nT)$ for all of them, much faster than
$O(n^{2}T)$. Finally, we note that Theorem 6.3 requires having an approximate
explicit matrix inverse. This can be done by matrix maintenance (Theorem B.4)
which takes $n^{2}L$ per phase. To see this, we note that Theorem B.4 takes
$O(n^{2})$ time per rank 1 update and it takes $O(n^{2})$ time per unit
multiplicative $\ell_{2}$ changes in the weight (if we pick $\epsilon=0.01$ in
Theorem B.4). Hence, each phase, the cost is just $O(n^{2}L)$. ∎

Algorithm 2 Batched Update of Leverage Score (Theorem 6.3)

1:Procedure BatchUpdate($\\{A^{(k)}\in\mathbb{R}^{m^{(k)}\times
n},w^{(k)}\in\mathbb{R}^{m^{(k)}}\\}_{k=0}^{T}$)

2:

3:$T_{0}\leftarrow T$ and $M^{(0,k)}\leftarrow
A^{(k)\top}W^{(k)}A^{(k)}$$\triangleright$ Easier to program via the $n\times
n$ matrices $A^{\top}WA$

4:

5:$\triangleright$ Step 1: Split the update into positive and negative update

6:$M^{(1,0)}\leftarrow M^{(0,0)}$

7:for $k=1,\cdots,T_{0}$ do

8: if $\operatorname{rank}(M^{(0,k)}-M^{(0,k-1)})=1$ then$\triangleright$
Insert/Delete Step

9: $M^{(1,\textsc{End}+1)}\leftarrow M^{(0,k)}$$\triangleright$ Append
$M^{(0,k)}$ to the end of the list of matrix $\\{M^{(0,:)}\\}$

10: else$\triangleright$ Update Step

11: Write $M^{(0,k)}=A^{\top}W^{(k)}A$ and $M^{(0,k-1)}=A^{\top}W^{(k-1)}A$

12: $M^{(1,\textsc{End}+1)}\leftarrow
A^{\top}\max(W^{(k)},W^{(k-1)})A$$\triangleright$ positive update:
$M^{(1,\textsc{End}+1)}\succeq M^{(1,\textsc{End})}$

13: $M^{(1,\textsc{End}+1)}\leftarrow M^{(0,k)}$$\triangleright$ negative
update: $M^{(1,\textsc{End}+1)}\preceq M^{(1,\textsc{End})}$

14: end if

15: $\textsc{End}\leftarrow\textsc{End}+1$

16:end for

17:

18:$\triangleright$ Step 2: Move insert / positive update into front

19:$M^{(2,0)}\leftarrow M^{(1,0)}$, $L\leftarrow\log^{3}(n)$$\triangleright$
$L$ is the size of the phase

20:for $k=1,\cdots,T_{1}$ where $M^{(1,T_{1})}$ is the last element in
$\\{M^{(1,:)}\\}$ do

21: if $M^{(1,k)}\succeq M^{(1,k-1)}$ then$\triangleright$ Insert / Positive
Update Step

22: $M^{(2,\textsc{End}+1)}\leftarrow
M^{(2,\textsc{End})}+M^{(1,k)}-M^{(1,k-1)}$$\triangleright$ Put it into front

23: $\textsc{End}\leftarrow\textsc{End}+1$

24: end if

25:end for

26:

27:Append enough $M^{(2,\textsc{End})}$ into $\\{M^{(2,:)}\\}$ such that the
last index is a multiple of $L$

28:for $k=1,\cdots,T_{1}$ do

29: if $M^{(1,k)}\prec M^{(1,k-1)}$ then$\triangleright$ Delete / Negative
Update Step

30: $M^{(2,\textsc{End}+1)}\leftarrow
M^{(2,\textsc{End})}+M^{(1,k)}-M^{(1,k-1)}$$\triangleright$ Put it into end

31: $\textsc{End}\leftarrow\textsc{End}+1$

32: end if

33:end for

34:Append enough $M^{(2,\textsc{End})}$ into $\\{M^{(2,:)}\\}$ such that the
last index is a multiple of $L$

35:

Algorithm 3 Batched Update of Leverage Score (Theorem 6.3, Continuation of
Algorithm 2)

36:$\triangleright$ Step 3: Within each phase, merge rank-1 update together
and sparsify non rank-1 update

37:$M^{(3,0)}\leftarrow M^{(2,0)}$, $\eta\leftarrow n^{-0.08}$$\triangleright$
Ignore all updates with multiplicative change less than $\eta$

38:for $\ell=1,\cdots,T_{2}/L$ do$\triangleright$ Due to line 27 and 34, each
phase is monotone

39: if Phase $\ell$ involves only insert / positive update then

40: $\Delta_{\text{Insert}}\leftarrow\sum_{k}(M^{(2,k)}-M^{(2,k-1)})$ where it
sums rank $1$ update in phase $\ell$

41: $\Delta_{\text{PosUpdate}}\leftarrow\sum_{k}(M^{(2,k)}-M^{(2,k-1)})$ where
it sums rank $>1$ update in phase $\ell$

42: $M^{(3,\textsc{End}+1)}\leftarrow
M^{(3,\textsc{End})}+\Delta_{\text{Insert}}$

43: Write $M^{(3,\textsc{End})}=A^{\top}WA$ and
$\Delta_{\text{PosUpdate}}=A^{\top}\operatorname{diag}(\delta)A$

44: $M^{(3,\textsc{End}+1)}\leftarrow
M^{(3,\textsc{End})}+A^{\top}\operatorname{diag}(\delta_{\text{sparse}})A$
where $\delta_{\text{sparse}}=\delta 1_{|\delta|\geq\eta|w|}$.

45: else

46: $I\leftarrow$ the list of rows deleted in phase $\ell$

47: $\Delta_{\text{NegUpdate}}\leftarrow\sum_{k}(M^{(2,k)}-M^{(2,k-1)})$ where
it sums rank $>1$ update in phase $\ell$

48: Write $M^{(3,\textsc{End})}=A^{\top}WA$ and
$\Delta_{\text{NegUpdate}}=A^{\top}\operatorname{diag}(\delta)A$

49: $M^{(3,\textsc{End}+1)}\leftarrow
M^{(3,\textsc{End})}+A^{\top}\operatorname{diag}(\delta_{\text{sparse}})A$
where $\delta_{\text{sparse}}=\delta 1_{|\delta|\geq\eta|w|}$.

50: $M^{(3,\textsc{End}+1)}\leftarrow$ $M^{(3,\textsc{End})}$ with rows in $I$
removed

51: end if

52: $v\leftarrow w+\delta_{\text{sparse}}$

53: $\textsc{End}\leftarrow\textsc{End}+1$

54:end for

55:

56:$\triangleright$ Step 4: Add up the estimate change of leverage scores

57:$c\leftarrow 0$

58:for $k=1,\cdots,T_{3}$ do

59: Write $M^{(3,k)}=A^{\top}W^{(k)}A$ and $M^{(3,k-1)}=A^{\top}W^{(k-1)}A$

60: if There is no insert/delete between $M^{(3,k)}$ and $M^{(3,k-1)}$
then$\triangleright$ Update

61:
$w^{(k-\frac{1}{2})}\leftarrow\text{median}(\frac{1}{2}w^{(k-1)},w^{(k)},2w^{(k-1)})$$\triangleright$
Entry-wise median

62: $\delta_{c}\leftarrow$ Estimate of
$\sigma(w^{(k-\frac{1}{2})})-\sigma(w^{(k-1)})$ using Lemma 6.2 with error
$\frac{1}{n^{1000}}$

63: $\delta_{c}\leftarrow\delta_{c}+$ Estimate of
$\sigma(w^{(k)})-\sigma(w^{(k-\frac{1}{2})})$ using Lemma 6.2 with error
$\frac{1}{n^{1000}}$

64: else$\triangleright$ Insert/Delete

65: $\delta_{c}\leftarrow$ Estimate of $\sigma(w^{(k)})-\sigma(w^{(k-1)})$
using Lemma 6.2 with error $\frac{1}{n^{1000}}$

66: end if

67: $c\leftarrow c+\delta_{c}$

68:end for

69:return $v,c$

70:EndProcedure

##  7 Simple Deterministic Leverage Score Maintenance

In this section, we give a simple deterministic leverage score maintenance
data-structure which is used by both the inner and middle phases in our main
data structure in Section 5. As a sub-procedure, we make use of the batched
low-rank update procedure in Section 6 which requires the following assumption
for a sequence of updates acts:

######  Assumption 7.1 (Short sequence and small updates).

Let the sequence of matrices and vectors in $\mathsf{acts}$ be
$A^{(0)},A^{(1)},\cdots,A^{(T)}\in\mathbb{R}^{m\times n}$ and
$w^{(0)},w^{(1)},\cdots,w^{(T)}\in\mathbb{R}^{m}$. We assume that $T\leq
n^{0.08}$, the maximum number of rows is bounded by $O(n)$, and that the
changes satisfy

  * •

Update: $\|\log w^{(k)}-\log w^{(k-1)}\|_{2}\leq 0.01$.

  * •

Insert/Delete: $w_{a}aa^{\top}\preceq
0.01(A^{(k-1)})^{\top}W^{(k-1)}A^{(k-1)}$ .

###  7.1 Main result

######  Theorem 7.2 (Simple leverage score maintenance).

Given an initial matrix $A\in\mathbb{R}^{m\times n}$ with $m=O(n)$, initial
weight $w\in\mathbb{R}_{+}^{m}$, and error parameter
$\epsilon_{\operatorname{simp}}\leq 0.01$. There is a deterministic data
structure (Algorithm 4) that approximately maintains the leverage scores

$\displaystyle\sigma_{i}(w)=(\sqrt{W}A(A^{\top}WA)^{-1}A^{\top}\sqrt{W})_{i,i},$

for positive diagonal matrices $W=\operatorname{diag}(w)$ through the
following operations:

  1. 1.

$\textsc{Init}(A,w,\epsilon_{\operatorname{simp}})$: takes
$O(n^{\omega+o(1)})$ time to initialize the data structure.

  2. 2.

$\textsc{Update}(\mathsf{acts})$: update the matrix $A$, vector $w$ and the
approximation of the leverage scores for the sequence of updates in
$\mathsf{acts}$. Moreover, if $\mathsf{acts}$ satisfies Assumption 7.1, then
the function $\textsc{Update}(\mathsf{acts})$ makes one call to
$\textsc{pm}.\textsc{update}$ 222The running time of procedure Update depends
on $\epsilon_{\mathrm{simp}}$, see details in Theorem B.1 and takes extra time
$O(n^{2}\cdot|\mathsf{acts}|)$ to compute a vector
$\Delta\widetilde{\sigma}\in\mathbb{R}^{m}$ such that

$\displaystyle\left\lVert\Delta\widetilde{\sigma}-\Delta\sigma\right\rVert_{2}\leq$
$\displaystyle\leavevmode\nobreak\
\frac{1}{n^{100}}+O\left(\epsilon_{\operatorname{simp}}+|\mathsf{acts}|\cdot
n^{-0.08}\right)\cdot|\mathsf{acts}|.$

  3. 3.

$\textsc{RefineEstimate}(\widetilde{\sigma}^{\operatorname{new}})$: takes
$O(n)$ time to update the approximation of $\sigma(w)$ to
$\widetilde{\sigma}^{\operatorname{new}}$.

  4. 4.

$\textsc{Query}()$: takes $O(n)$ time to output the approximation
$\widetilde{\sigma}\in\mathbb{R}^{m}$ maintained by the data structure.

###### Proof of Theorem 7.2.

We only need to prove the guarantee for the function
$\textsc{Update}(\mathsf{acts})$, as the guarantees for other functions are
straightforward. Notice that each call to the function
$\textsc{Update}(\mathsf{acts})$ makes at most one call to
$\textsc{pm}.\textsc{update}$, and the extra running time of
$O(n^{2}\cdot|\mathsf{acts}|)$ follows directly from Theorem 6.3. To prove the
error guarantee, we let the sequence of matrices and vectors in
$\mathsf{acts}$ be $A^{(0)},A^{(1)},\cdots,A^{(T)}$ and
$w^{(0)},w^{(1)},\cdots,w^{(T)}$. notice that by Theorem 6.3, we have

$\displaystyle\left\lVert\sigma_{A^{(T)}}(w^{\mathrm{mid}})-\sigma_{A^{(0)}}(w^{(0)})-\Delta\widetilde{\sigma}^{\mathrm{mid}}\right\rVert_{2}\leq\frac{1}{n^{100}}.$

It follows from Lemma 7.6 that

$\displaystyle\left\lVert\sigma_{A^{(T)}}(w^{(T)})-\sigma_{A^{(T)}}(w^{\mathrm{mid}})-\Delta\widetilde{\sigma}^{\operatorname{new}}\right\rVert_{2}\leq$
$\displaystyle\leavevmode\nobreak\
O\left(\epsilon_{\operatorname{simp}}+|\mathsf{acts}|\cdot
n^{-0.08}\right)\cdot|\mathsf{acts}|.$

The error guarantee then follows by summing up the two error upper bounds
above. ∎

The following corollary is an immediate consequence of Theorem 7.2 and B.1. It
states that if matrix multiplication can be performed fast enough, then the
simple leverage score maintenance would imply a deterministic
$O(n^{3}\log(n/\epsilon))$ algorithm for the cutting plane method.

######  Corollary 7.3.

If $\mathcal{T}_{\operatorname{mat}}(n,n,r)=O(n^{2}\log^{O(1)}(n))$ for
$r=n^{\alpha}$ with $\alpha>2/3$, then there’s a deterministic
$O(n^{3}\log(n/\epsilon))$ time algorithm for the cutting plane method.

###### Proof.

This algorithm only uses the inner and middle phases. The inner phase is run
as is our main algorithm. For the middle phase, we notice that, in this case,
the exponent of matrix multiplication time can be bounded as
$\omega<3-2/3=7/3$. Therefore, the data structure is restarted after every
$n^{1/3}$ calls to the middle phase. For the middle phase, we use the simple
leverage score maintenance data structure in Theorem 7.2 with parameter
$\epsilon_{\operatorname{simp}}=n^{-\alpha/2}/\log^{O(1)}(n)$. It follows from
Theorem 7.2 that the error accumulated in the $n^{0.2}$ steps is
$O(n^{-\alpha/2}n^{1/3}\log^{O(1)}(n))=n^{-\Omega(1)}$. By Theorem 7.2 and
5.1, the running time is $O(n^{2}\log^{O(1)}(n))$ per step for middle phase
which is fast enough. This finishes the proof of the corollary. ∎

Algorithm 4

1:data structure MaintainLeverageScoreSimple $\triangleright$ Theorem 7.2

2:

3:members

4: $m\in\mathbb{Z}_{+}$ $\triangleright$ $m=\Theta(n)$

5: $A\in\mathbb{R}^{m\times n}$

6: $w\in\mathbb{R}^{m}$ $\triangleright$ Target vector

7: $\widetilde{\sigma}\in\mathbb{R}^{m}$ $\triangleright$ The approximate
leverage scores

8: ProjectionMaintain pm $\triangleright$ Theorem B.2 or B.4

9:end members

10:

11:procedure Init($A\in\mathbb{R}^{m\times
n},w\in\mathbb{R}^{m}_{+},\epsilon_{\operatorname{simp}}\in\mathbb{R}_{+}$)
$\triangleright$ Initialization

12: $A\leftarrow A$, $m\leftarrow m$, $w\leftarrow w$

13: $\widetilde{\sigma}\leftarrow\sigma(w)$

14: $\textsc{pm}.\textsc{init}(A,w,\epsilon_{\operatorname{simp}})$

15:end procedure

16:

17:procedure
Update($\mathsf{acts}\in\textsc{Vector}\langle\textsc{Action}\rangle$)
$\triangleright$ Update estimate for a sequence of updates

18: Let $A$ and $w^{\operatorname{new}}$ be the final matrix and vector in
$\mathsf{acts}$

19:
$w^{\mathrm{mid}},\Delta\widetilde{\sigma}^{\mathrm{mid}}\leftarrow\textsc{BatchedUpdate}$($\mathsf{acts}$)
$\triangleright$ Let
$w^{\mathrm{mid}},\Delta\widetilde{\sigma}^{\mathrm{mid}}$ be computed as in
Theorem 6.3 (Algorithm 2 and 3) with the sequence of updates in
$\mathsf{acts}$

20: $v^{\mathrm{mid}},M(v^{\mathrm{mid}})^{-1},Q(v^{\mathrm{mid}})\leftarrow$
pm.update$(w^{\mathrm{mid}})$

21: Compute $\Delta\widetilde{\sigma}^{\operatorname{new}}$ as in Lemma 7.5

22:
$\widetilde{\sigma}\leftarrow\widetilde{\sigma}+\Delta\widetilde{\sigma}^{\mathrm{mid}}+\Delta\widetilde{\sigma}^{\operatorname{new}}$

23: $w\leftarrow w^{\operatorname{new}}$

24:end procedure

25:

26:procedure
RefineEstimate($\widetilde{\sigma}^{\operatorname{new}}\in\mathbb{R}_{+}^{m}$)
$\triangleright$ Refine the estimation $\widetilde{\sigma}$ to
$\widetilde{\sigma}^{\operatorname{new}}$

27: $\widetilde{\sigma}\leftarrow\widetilde{\sigma}^{\operatorname{new}}$

28:end procedure

29:

30:procedure Query()

31: return $\widetilde{\sigma}$

32:end procedure

33:

34:end data structure

###  7.2 Approximate leverage score’s moving

######  Lemma 7.4 (Leverage score’s moving).

Given any matrix $A\in\mathbb{R}^{m\times n}$, vectors
$w^{\mathrm{mid}},w^{\operatorname{new}}\in\mathbb{R}_{+}^{m}$. The change in
leverage score from $w^{\mathrm{mid}}$ to $w^{\operatorname{new}}$ can be
written as

$\displaystyle\Delta\sigma^{\operatorname{new}}_{i}:=$
$\displaystyle\leavevmode\nobreak\
\sigma(w^{\operatorname{new}})_{i}-\sigma(w^{\mathrm{mid}})_{i}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
(w^{\operatorname{new}}_{i}-w^{\mathrm{mid}}_{i})\cdot\tau(w^{\mathrm{mid}})_{i}+w^{\operatorname{new}}_{i}\cdot\int_{0}^{1}e_{i}^{\top}Q(x_{t})(W^{\mathrm{mid}}-W^{\operatorname{new}})Q(x_{t})e_{i}\mathrm{d}t,$

where we define
$x_{t}=w^{\mathrm{mid}}+t(w^{\operatorname{new}}-w^{\mathrm{mid}})$ for each
$t\in[0,1]$.

###### Proof.

We can rewrite
$\sigma(w^{\operatorname{new}})_{i}-\sigma(w^{\mathrm{mid}})_{i}$ as

$\displaystyle\sigma(w^{\operatorname{new}})_{i}-\sigma(w^{\mathrm{mid}})_{i}=(w^{\operatorname{new}}_{i}-w^{\mathrm{mid}}_{i})\cdot\tau(w^{\mathrm{mid}})_{i}+w^{\operatorname{new}}_{i}(\tau(w^{\operatorname{new}})_{i}-\tau(w^{\mathrm{mid}})_{i}).$

Then the second term is computed as

$\displaystyle\leavevmode\nobreak\
w^{\operatorname{new}}_{i}(\tau(w^{\operatorname{new}})_{i}-\tau(w^{\mathrm{mid}})_{i})$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
w^{\operatorname{new}}_{i}\cdot((A(A^{\top}W^{\operatorname{new}}A)^{-1}A^{\top})_{i,i}-(A(A^{\top}W^{\mathrm{mid}}A)^{-1}A^{\top})_{i,i})$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
w^{\operatorname{new}}_{i}\cdot\int_{0}^{1}e_{i}^{\top}Q(x_{t})(W^{\mathrm{mid}}-W^{\operatorname{new}})Q(x_{t})e_{i}\mathrm{d}t.$

where the last step follows by the definition of
$Q(v)=A(A^{\top}VA)^{-1}A^{\top}\in\mathbb{R}^{m\times m}$. ∎

######  Lemma 7.5 (Computing $\Delta\widetilde{\sigma}^{\operatorname{new}}$).

Given any matrix $A\in\mathbb{R}^{m\times n}$, vectors
$w^{\mathrm{mid}},w^{\operatorname{new}}\in\mathbb{R}_{+}^{m}$, where