An Improved Cutting Plane Method for Convex Optimization, Convex-Concave Games and its Applications

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$. $\bullet$ 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 $\text{SO}$ is the time for oracle evaluation. $\bullet$ 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 $\mathrm{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.


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 ω+1 ) time.
Here ω < 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/κ, namely log 2 (κ). The following table compares the runtimes of Vaidya's and LSW methods, which both achieve the optimal number of oracle calls of n log(κ).

Authors
#Oracle Calls Runtimes Vaidya n log(κ) n ω+1 log(κ) Lee-Sidford-Wong n log(κ) n 3 log O(1) (n) log 3 (κ) This extra overhead of log 2 (κ), 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 = 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(κ) vs log 3 (κ). For many other problems including linear programming and market equilibrium computation, must be taken as 1/M O(n) resulting in an additional factor n 2 between log(κ) and log 3 (κ). For these applications, LSW is slower than Vaidy's by a factor of O(n 4−ω ). This raises a natural question: is there a cutting plane method that simultaneously runs in O(n log(κ)) calls and O(n 3 log O(1) (n) log(κ)) 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 κ?
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 · SO log(κ) + n 3 log(κ)), where SO is the time complexity of the separation oracle.

Reference
Year # Operations Poly Type [Vaz10] 2010 n 3 (n + m) 2 log(U ) · T 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. T 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.
As with previous methods, our result achieves the asymptotic optimal oracle complexity of O(n · SO log(κ)). 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.

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 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 R n and S be a convex set that contains a minimizer of f . Suppose we have a subgradient oracle for f with cost T and S ⊂ B(0, R). Using B(0, R) as the initial polytope for our Cutting Plane Method, for any 0 < α < 1, we can compute x ∈ S such that f (x) − min y∈S f (y) ≤ α (max y∈S f (y) − min y∈S f (y)), with high probability in n and with a running time of O(T · n log(κ) + n 3 log(κ)), where κ = nγ/α and γ = R/(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 X ⊂ B(0, R) ⊂ R n and Y ⊂ B(0, R) ⊂ R m such that both X and Y contain a ball of radius r. Let f (x, y) : X × Y → R be an L-Lipschitz function that is convex in x and concave in y. Define κ = n+m R r . For any ∈ (0, 1/2], we can find ( x, y) such that max y∈Y f ( x, y) − min x∈X f (x, y) ≤ Lr in time O T · (n + m) log(κ) + (n + m) 3 log(κ) with high probability in n + m where T is the cost of computing subgradient ∇f .
Leveraging this improved dependence on κ (and hence ), 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/ needs to be exponentially large in n which renders the log(κ) 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 Θ(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 T AD log(SM n)+n 4 log(SM n)) for computing a market equilibrium in an economy with general buyer valuation in the aggregate demand model.

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 ∈ K or prove that K does not contain a ball of radius .
All cutting plane methods maintain a candidate region Ω and solve the feasibility problem by iteratively refining Ω based on the present Ω and the new separating hyperplane. In each iteration: 1. The separation oracle is queried at some point x ∈ Ω.
2. If x ∈ K we have solved the feasibility problem. 3. Otherwise, the separation oracle returns a separating hyperplane from which Ω is further refined and the next query point x is computed.
Previous works differ in how x is selected and how Ω is refined. For instance, the classic ellipsoid method maintains Ω as an ellipsoid and x as its center. Given Ω and the new separating hyperplane, the new Ω 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 Ω(n log(κ)) [NY83]. While the ellipsoid method achieves a suboptimal O(n 2 log(κ)) 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(κ)) 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 Ω, similar to the random walk method, is also a polytope Ω = {x ∈ R n : Ax ≥ 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 Ω: is the diagonal matrix of the slacks Nevertheless, by judiciously including only a representative subset of previous half-spaces ax ≥ b, Vaidya showed that the volumetric center and Ω can be updated by basic matrix operations which run in O(n ω ) 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.

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 ω ) 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 ω ) time. LSW further overcame this barrier by resorting to a recent work that efficiently solves "slowly-changing" linear system in amortized O(n 2 ) time. Thus after paying O(n ω ) initially, they can solve such linear systems in 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 0 in JL projection has to be as small as 0 = 1/ log(κ). As the runtime of JL depends on 1/ 2 0 , this unfortunately introduces the log 2 κ overhead in LSW runtime when compared to Vaidya's volumetric center method.

Overview of our approach
To achieve the desired O(n SO log(κ) + n 3 log(κ)) runtime, we build a data structure that approximates changes in leverage scores to within c in 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 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.

Discussion of optimality
Similar to previous methods our cutting plane method achieves the optimal oracle complexity n log(κ) [NY83]. We present some evidence that our running time of O(n 3 log(κ)) 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 , · · · w T , let P (w) be defined as 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 ω time so the goal is to achieve o(n ω ) amortized cost per iteration. For example, in the LP setting the number of iterations is O( √ n) and Vaidya [Vai89b] combined fast matrix multiplication with inverse maintenance to achieve O(n 2 ) amortized cost per iteration, 8 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 ω ) time algorithm. One of the major computation required in each step is matrix-vector multiplication, e.g., P (w)·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 ∈ R n×n which has roughly √ n nonzeros.
They use that sampling matrix to sample on the right hand side: In [LSZ19], they use the subsampled randomized Hadamard/Fourier matrix [LDFU13, PSW17] R ∈ R √ n×n . They use the sketching matrix to sketch on the left hand side: 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 W A) −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 Ω(n 3 log(κ)) time. Hence our cutting plane method achieves the optimal runtime. 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 ∈ K or x / ∈ K.

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

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 coarsegrained 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 ω−2+o(1) ) steps as exact computation of leverage scores takes O(n ω+o(1) ) time. Therefore it suffices to control the error accumulation in O(n ω−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 ω−2+o(1) ) steps are too large for the cutting plane method. For the complicated data structure, we are able to achieve 2 -error o(1) in O(n ω−2+o(1) ) steps, but the amortized running time would be O(n 2+o(1) ) per step.

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 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 T mat (n, n, r) the time to compute the multiplication of an n × n matrix and an n × r matrix 1 . Then we have: 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. .
out restart refine σmid restart restart restart Figure 1: Illustration of our three-level data structure with T inn = 2, T mid = 3 and T out = 4 approximating the leverage scores of the sequence {w 0 , w 1 , · · · }. The three levels maintain three approximate sequences {v 0 inn , v 1 inn , · · · }, {v 0 mid , v 1 mid , · · · } and {v 0 out , v 1 out , · · · }, with errors inn = log(w) − log(v inn ) ∞ , mid = log(w) − log(v mid ) ∞ and out = log(w) − log(v out ) ∞ that satisfy 1 inn mid out > 0. The inner level takes a step for every w-update, the middle step takes a step in every T inn inner steps, and the outer step takes a step in every T mid middle steps. The entire data structure is restarted after T 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 inn , T mid , T out and inn , mid , out in our data structure, see Table 5.

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 ∈ R m + satisfies an 2closeness property: log(w new ) − log(w) 2 = O(1). For a sequence of T updates w (0) , w (1) , · · · , w (T ) , the 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 mid that is ∞ -close to the vector w (T ) and an accurate estimate of the leverage score change σ(v mid ) − σ(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 poly log(n). We ensure the weight changes by only a quasi-polynomial factor and this decreases the cost of solving linear systems to poly log(n) steps. Since we batch the task of handling 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).

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 new , one of the integral terms is In order to approximate such an integral, we take a set T ⊆ [0, 1] of N = log O(1) (n) points along the integration together with weights {ω t } t∈T as in Theorem 3.10, and approximate the integral by Theorem 3.10 then shows that the 2 -error can be bounded as σ cts − σ dis 2 ≤ 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 , · · · w T ∈ R m + , let P (w) be defined as , and efficiently computes P (v t )u, for query vector u ∈ 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 W A) −1 A and M (w) −1 = (A W A) −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 γ i,s,t for all i, where Notice that computing γ i,s,t for all i is essentially computing the matrix products which would take O(n ω+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 where R γ,s,s ,t ∈ R n c ×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, σ i,jl should be a good estimate with a small variance. We note that σ i,jl is indeed an unbiased estimator of σ i,dis and its variance can be bounded as m i=1 Var[σ i,jl ] ≤ O( 2 /n c ), where 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 = n −0.1 . The variance would then be bounded by n −0.51 which is sufficiently small after n ω−2+o(1) steps before the data structure restarts.

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 T mat (n, n, r) = O(n 2 log O(1) (n)) for r = n β with β > 2/3, then there is a deterministic cutting plane method which runs in time where SO is the time complexity of the separation oracle.
Let α denote the dual exponent of matrix multiplication, which is the largest number α > 0 such that T mat (n, n, n α ) = n 2+o(1) . Let β denote the largest number such that T mat (n, n, n β ) = 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 α < 0.625 for certain tensors. We believe our work initiated two interesting open questions in the area of fast matrix multiplication:

13
(1) whether one can prove a better upper bound on β (compared to α) for certain tensor techniques, and (2) if there is a non-trivial inequality between β and α.
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 datastructure 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. [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.
[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.

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

Notations
For a positive integer n, let [n] denote the set {1, 2, · · · , n}. For any function f , we define O(f ) to be f · log O(1) (f ). In addition to O(·) notation, for two functions f , g, we use the shorthand f g (resp. ) to indicate that f ≤ C · g (resp. ≥) for some absolute constant C.
For vector x ∈ R n , we use x 2 to denote its 2 norm, we use x 1 to denote its 1 norm, we use x ∞ to denote its ∞ norm. For matrix A ∈ R m×n , we use A to denote the spectral norm of A, we use A F to denote its Frobenius norm ( For square matrix A, we use tr[A] to denote the trace of A. We say A is positive-semidefinite (psd) if A = A and x Ax ≥ 0, ∀x ∈ R n . We use , to denote the semidefinite ordering, e.g. A 0 means that A is psd.
For square full-rank matrix A, we use A −1 to denote the inverse of A.

Operators
Definition 3.1 (Operators τ and σ). Given a matrix A ∈ R m×n and a vector v ∈ R m , we define the leverage scores σ(v) i and (unnormalized) leverage scores τ (v) i as follows Definition 3.2 (Operators M , Q and P ). Given a matrix A ∈ R m×n , we define three operators :

Basic results on matrices
Fact 3.6 (Woodbury matrix identity, [Woo49, Woo50]). Given a square invertible n × n matrix A, an n × k matrix U and a k × n matrix V , let B be an n × n matrix such that B = A + U CV . Then,

Fast matrix multiplication
Definition 3.7 (Matrix multiplication time). For any n, r > 0, denote T mat (n, n, r) the time to compute the multiplication of an n × n matrix and an n × r matrix.

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 : R → R. There exits N points s 1 , · · · , s N and N weights ω 1 , · · · , ω N ≥ 0 such that i ω i = 1 and Now, we explain how to use one variable interpolation result (Theorem 3.9) to prove a multiple variable interpolation result (Theorem 3.10).
Then there exists N points s 1 , · · · , s N and N weights ω 1 , · · · , ω N such that Using one variable theorem 3.9, we have that : for all t 1 , t 2 , · · · , t d−1 ∈ [0, 1], 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 , · · · , t j .
For each j ∈ {d − 1, · · · , 1}, we define function g j : R j → R such that 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 ∈ [d], Proof. First, using one variable theorem 3.9, we can upper bound

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 ∈ K or prove that K does not contain a ball of radius .
Let K ⊂ 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 Ω and a query point z. In each iteration, z is the input to the separation oracle which then either certifies z ∈ K or produces a hyperplane a x = b separating z from K. The algorithm then refines Ω and computes the next query point z (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 ⊂ R n that is contained in a box of radius R and a parameter > 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 in O((n SO +n 3 ) log(κ)) time, where SO is the complexity of the separation oracle and κ = nR/ .
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.

26
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 SO log(κ) + n 3 log(κ)) 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 ∈ R m×n with m = O(n), initial weight w ∈ R m + . Then for any constant c > 0, there is a randomized data structure (Algorithm 1) that approximately maintains the leverage scores for positive diagonal matrices W = diag(w) through the following operations: (1) ) time to intialize the data structure.
2. Update(act): updates the data structure for the single update act.

Update w
such that m (k−1) = m (k) and that log(w (k) ) − log(w (k−1) ) 2 ≤ 0.01, then the function Update(act) takes an amortized O(n 2 ) time, and the vector σ (k) output by Query() at each step k ∈ [K] satisfies that 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 Init(A, w), and Query() are immediate corollaries of Theorem 7.2 and 8.1. We prove running time upper bound for the function Update(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 Update(act) makes one call of Update(act) with act having a single action act, it follows from Theorem 7.2 that the inner phase makes one call to pm.update and takes extra time at most O(n 2 ). Since we choose our parameter as 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 pm.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 inn = log 10 (n) ≤ n 0.08 calls to Update(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 inn calls to the function Update(act), the middle phase makes one call to pm.update and takes extra time at most O(n 2 T inn ). Since we choose mid = 1/n 0.08 , it follows from Theorem B.2 with C = T inn = log 10 (n) that the time for one call to pm.update is O(n 2 log 2 (n)). Amortized over the T inn = log 10 (n) calls to the function Update(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 inn · T mid = n 0.01 log 10 (n) calls to Update(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 inn · T mid calls to the function Update(act), the outer phase makes O(N ) = O(log 2 (n)) calls to pm.update and takes extra time at most O(n 2 · T inn · T mid + T (n, n, r out ) · N 3 · log(n)) = O(n 2 · T inn · T mid + n 2+o(1) log 7 (n)).
Since we choose out = 1/n 0.1 , it follows from Theorem B.2 with C = T inn · T mid that the time for O(N ) calls to pm.update is O(n 2+o (1) ). Amortized over the T inn · T mid = n 0.01 log 10 (n) calls to the function Update(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 ω+o(1) ) every T inn · T mid · T out = Ω(n ω−1.99 ) calls to the function Update(act). It follows that the amortized time to restart the data structure is O(n ω+o(1) )/(T inn · T mid · T 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 2 -error of the inner phase is O( inn ) = O(log −25 (n)) per call to Update(act) and the inner phase runs for at most T 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( inn ) · T inn = O(log −15 (n)). For the middle phase, it follows from Theorem 7.2 that the 2 -error introduced by each middle step is at most O(( mid + T inn · n −0.08 ) · T inn ) ≤ n −0.07 . Since the middle phase run for at most T 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 · 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 Since the outer phase run for at most T out = n ω−2 steps before the data structure gets restarted, the variance accumulated in the T out steps is at most O(n −0.1 ). Therefore, we have σ (k) − σ(w (k) ) 2 ≤ O(log −15 (n)) for each k ∈ [K]. This completes the proof of Theorem 5.1. ctr inn , ctr mid , ctr out ∈ Z + Step counters for different phases 5: acts mid , acts out ∈ Vector Action Delayed sequence of updates 6: Action includes insertion, deletion and w-update 7: MaintainLeverageScoreSimple simp inn , simp mid For inner and middle phases, Theorem 7.2, Algorithm 4 8:

MaintainLeverageScoreComplicated comp out
For outer phase, Theorem 8.1, Algorithm 5 9: end members 10: n ω−2 Number of outer step before re-computing r out n 0.31 Matrix multiplication dimension for outer step inn 1/ log 25 n Inner step projection maintenance error mid 1/n 0.08 Middle step projection maintenance error 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.
6 Batched Low Rank Update

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 × n PSD matrices A and M such that 0 M A κ · M . For any integer t ≥ 1 such that Furthermore, if M −1 and A are given explicitly, for any V ∈ R n×r , then we can compute f (M, t) · V in O(T mat (n, n, r) · t) time.
Proof. Note that Using Multiplying 1 κ M −1 on both sides of the above equation with a formula A −1 (Eq. (2)), we get Using the definition of f (M, t), we have that 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 ∈ R m×n , non-negative weights w, w new ∈ R m such that w new ≥ w with m = O(n) and w new − w 0 = k. Let β > 1 denote the parameter such that A W new A β · A W A. Given some explicit matrix U ∈ R n×n and U new ∈ R n×n such that Then, for any ∈ (0, 1/2), we can output a vector c ∈ R m such that O(T mat (n, n, k) · (1 + log log(n) (βk/ ))).
The same statement holds for decreasing w, namely Proof. We note that the decreasing case follows from the increasing case by swapping w and w new . Hence, we focus on the increasing case. Using Lemma 6.1, we can construct U such that U is low-degree polynomial of U and A W A and that We define U new similarly. We will show how to use U and U new to approximate σ(w new ) − σ(w). First, we note that We define c i , c 1,i , c 2,i as follows: We estimate the right hand side by where Cost of computing c For the first term in (5). We note that w new − w has only k non-zeros. Hence, we only need to compute AU new A on k of the diagonals. Let A S ∈ R k×n where each rows of A that w new = w. Then, we can compute the k of the diagonals via the whole matrix Lemma 6.1 shows that U new has O(1 + log log log(n) ) terms and that it takes O T mat (n, n, k) · 1 + log log log(n) to compute U new A S . Finally, it takes extra T mat (k, n, k) = O(T mat (n, n, k)) to do the left multiplication on A S .
For the second term in (5). Woodbury matrix identity shows that where ∆ W ∈ R k×k be W new −W restricted on non-zeros. By same argument above, we can compute A S U A S in O T mat (n, n, k) · 1 + log log log(n) time.
Then, we can compute and Now, we multiple the two terms above together in time T mat (n, k, n) = O(T mat (n, n, k)). Hence, the total time is O T mat (n, n, k) · 1 + log log log(n) = O(T mat (n, n, k)(1 + log log(n) (βk/ ))).
The accuracy of c For the first term in c, we have where we used that w new ≥ w in the second inequality and leverage score is upper bounded by 1 in the last equality. Hence the ∞ norm error of the first term is given by . Hence, the 2 norm error is given by For the second term in c, we have that the error is given by where M = A W A. To simplify the notation, we define M t,s = M +t·∆ 2 +s·∆ where ∆ 2 = U −1 −M . Then, the error term becomes Note that Combining (7) and (8), we have that the 2 norm error for the second term is bounded by To bound the Frobenius norm, we note that ∆ 0 by the assumption and that ∆ 2 = U −1 − M 0 (4). Hence, we have that M t,s M and

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 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 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 W A can change by a factor of e Θ(T ) after T iterations. So, Lemma 6.2 takes at least Ω(T mat (n, n, k) T log log(n) ) time. This forces us to do rank k changes where T mat (n, n, k) ≤ n 2 log log(n). Hence, for 2 update, we need to handle coordinates with notso-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 Θ(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 ∈ R m×n , initial weight w ∈ 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 ∈ R n to the matrix A ∈ R m×n and weight w a to the vector w • Delete: Delete a row a ∈ R n from the matrix A ∈ R m×n and weight w a from the vector w Let the sequence of matrices and vectors be A (0) , A (1) , · · · , A (T ) ∈ R m×n and w (0) , w (1) , · · · , w (T ) ∈ R m . We assume that T ≤ n 0.08 , the maximum number of rows is bounded by O(n), and that the changes satisfy Then, there is an algorithm ( BatchedUpdate in Algorithm 2 and 3) that finds a vector v ∈ R n and a vector c ∈ R n such that 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 η = 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 T n −0.08 . Therefore, we have for some v such that log v − log w (T ) ∞ ≤ T n −0.08 . Furthermore, we output the good enough approximation of σ A (T ) (v) − σ A (0) (w (0) ) = σ(M (3,End) ) − σ(M (3,0) ) where we abused the notation to use σ(M ) to denote the leverage score of the rows insides M . This is simply because we split the difference into σ(M (3,End) ) − σ(M (3,0) ) = T k=1 σ(M (3,k) ) − σ(M (3,k−1) ) and estimate each step with 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) M (0,k) 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 ∆ (1,k) = M (1,k+1) −M (1,k) . By the relation, we have that Note that step 2 simply permutes ∆ (1,k) , namely, for each k there is an unique k such that ∆ (1,k) = ∆ (2,k ) . Finally, we note that M (1,k) M (2,k ) for all k because M (2,k ) is simply M (2,k ) 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,k) be two matrices in the phase. Then, we have that 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 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 O(T mat (n, n, k)(1 + log log(n) (βk/ ))) where = n −1000 , β = 2 O(L) and k = O(L). By the choice of L = log 3 (n), we have T 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 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− 1 2 ) , we have that log w (k) − log w (k− 1 2 ) 2 = O(L) and that | log w i | 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− 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 O(T mat (n, n, k)(1 + log log(n) (βk/ ))) where = n −1000 , β = 2 O(L) and k = O(L 2 ). By the choice of L = log 3 (n), we have T mat (n, n, k) = n 2 and hence the cost is 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 2 changes in the weight (if we pick = 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) Easier to program via the n × n matrices A W A 4:

56:
Step 4: Add up the estimate change of leverage scores 57: c ← 0 58: for k = 1, · · · , T 3 do 59: if There is no insert/delete between M (3,k) and M (3,k−1) then Update 61: δ c ← Estimate of σ(w (k− 1 2 ) ) − σ(w (k−1) ) using Lemma 6.2 with error 1 n 1000 63:  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 subprocedure, 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 acts be A (0) , A (1) , · · · , A (T ) ∈ R m×n and w (0) , w (1) , · · · , w (T ) ∈ R m . We assume that T ≤ n 0.08 , the maximum number of rows is bounded by O(n), and that the changes satisfy

Main result
Theorem 7.2 (Simple leverage score maintenance). Given an initial matrix A ∈ R m×n with m = O(n), initial weight w ∈ R m + , and error parameter simp ≤ 0.01. There is a deterministic data structure (Algorithm 4) that approximately maintains the leverage scores for positive diagonal matrices W = diag(w) through the following operations: (1) ) time to initialize the data structure.
2. Update(acts): update the matrix A, vector w and the approximation of the leverage scores for the sequence of updates in acts. Moreover, if acts satisfies Assumption 7.1, then the function Update(acts) makes one call to pm.update 2 and takes extra time O(n 2 · |acts|) to compute a vector ∆ σ ∈ R m such that ∆ σ − ∆σ 2 ≤ 1 n 100 + O simp + |acts| · n −0.08 · |acts|.
3. RefineEstimate( σ new ): takes O(n) time to update the approximation of σ(w) to σ new . 4. Query(): takes O(n) time to output the approximation σ ∈ R m maintained by the data structure.
Proof of Theorem 7.2. We only need to prove the guarantee for the function Update(acts), as the guarantees for other functions are straightforward. Notice that each call to the function Update(acts) makes at most one call to pm.update, and the extra running time of O(n 2 · |acts|) follows directly from Theorem 6.3. To prove the error guarantee, we let the sequence of matrices and vectors in acts be A (0) , A (1) , · · · , A (T ) and w (0) , w (1) , · · · , w (T ) . notice that by Theorem 6.3, we have It follows from Lemma 7.6 that The error guarantee then follows by summing up the two error upper bounds above.

39
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/ )) 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 ω < 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 simp = n −α/2 / log O(1) (n). It follows from Theorem 7.2 that the error accumulated in the n 0.2 steps is O(n −α/2 n 1/3 log O(1) (n)) = n −Ω(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.

Approximate leverage score's moving
Lemma 7.4 (Leverage score's moving). Given any matrix A ∈ R m×n , vectors w mid , w new ∈ R m + . The change in leverage score from w mid to w new can be written as Then the second term is computed as where the last step follows by the definition of Lemma 7.5 (Computing ∆ σ new ). Given any matrix A ∈ R m×n , vectors w mid , w new ∈ R m + , where m = O(n). Define an approximation ∆ σ new i to the change in leverage score to be  w mid , ∆ σ mid ← BatchedUpdate(acts) Let w mid , ∆ σ mid be computed as in Theorem 6.3 (Algorithm 2 and 3) with the sequence of updates in acts Proof. Since we know the matrix Q(v mid ), the first term is essentially reading all the diagonal entries of matrix Q(v mid ), which can be computed for all i ∈ [m] in time O(n). For the second term, we denote Q = Q(x t ) and Q (2) the entry-wise square of Q. Then we have Therefore, in order to compute the second term for all i ∈ [m], it suffices to compute the matrixvector product Q (2) (w mid − w new ), which takes time O(n 2 ).

2 -error bound
Let A (T ) and A (0) be the final and initial matrix in the sequence of updates in acts, and w (0) be the initial vector in acts. 0) ). It follows from Theorem 6.3 that if Assumption 7.1 holds, then we have It therefore suffices to bound the 2 -error ∆ σ new − ∆σ new 2 . For simplicity, we use A to denote A (T ) . Recall from Lemma 7.4 that where x t = w mid + t(w new − w mid ) for each t ∈ [0, 1]. Also recall that in Lemma 7.5, our approximation ∆ σ new is computed as We use ∆ 1 and ∆ 2 to denote the two corresponding error terms We prove the following bound on the 2 -error.
Lemma 7.6 ( 2 -error bound). Assume Assumption 7.1 holds and that simp ≤ 0.01, where simp is the error parameter in Algorithm 4. Then the 2 -error is upper bounded as Proof. For simplicity, we assume here that w new − w mid ≥ 0 and that w mid − v mid ≥ 0. This is without loss of generality, since in the case where these assumptions don't hold, we can break the corresponding terms into positive and negative parts and bound each one of them. The lemma then follows immeditely from the following Lemma 7.7 and 7.8.
Lemma 7.7 ( ∆ 1 2 ). Assume Assumption 7.1 holds and that simp ≤ 0.01. Then we have the following upper bound on ∆ 1 2 Proof. Recall the definition of the error term ∆ 1 as We define y s = w mid + s(v mid − w mid ) for s ∈ [0, 1]. The term (∆ 1 ) i can be rewritten as for some ξ ∈ [0, 1]. It follows that where in the last step we define diagonal matrices B, C ∈ R m×m as and we use P to denote P (y ξ ). Then we have where the last step follows from Theorem 6.3.

43
Proof. Recall the definition of the error term ∆ 2 as where in the last line we define v t,s = x t + s(v mid − x t ) for s ∈ [0, 1]. Therefore, in order to bound the ∆ 2 i , it suffices to bound the 2 -norm of the derivative: It follows that for any s ∈ [0, 1], we have where in the last step we define diagonal matrices B, C ∈ R m×m as and we use P to denote P (v t,s ). By Assumption 7.1, we have This together with the assumption that simp ≤ 0.01 further implies that The lemma then follows from Theorem 6.3 which gives log(w new ) − log(w mid ) ∞ ≤ |acts| · n −0.08 and log(w new ) − log(w mid ) 2 ≤ 0.01|acts|.

Complicated Randomized Leverage Score Maintenance
In this section, we present a complicated randomized leverage score maintenance data structure which is used by the outer phase in our main data structure in Section 5. We again make crucial use of the batched low-rank update procedure in Section 6.

Main result
Theorem 8.1 (Complicated leverage score maintenance). Given an initial matrix A ∈ R m×n with m = O(n), initial weight w ∈ R m + , JL dimension r > 0, discrete sampling parameter N ∈ Z + , and error parameter comp < 0.01. There is a randomized data structure (Algorithm 5) that approximately maintains the leverage scores for positive diagonal matrices W = diag(w) through the following operations: 1. Init(A, w, r, N, comp ): takes O(n ω+o(1) ) time to initialize the data structure.
2. Update(acts): update the matrix A, vector w and the approximation of the leverage scores for the sequence of updates in acts. Moreover, if acts satisfies Assumption 7.1, then the function Update(acts) makes O(N ) calls to pm.update and takes extra time O(n 2 · |acts| + T mat (n, n, r) · N 3 · log(n)) to compute a random vector ∆ σ ∈ R m with error upper bounded by

and variance upper bounded by
3. RefineEstimate( σ new ): takes O(n) time to update the approximation of σ(w) to σ new .
4. Query(): takes O(n) time to output the approximation σ ∈ R m maintained by the data structure.
Proof of Theorem 8.1. We only need to prove the guarantee for the function Update(acts), as the guarantees for other functions are straightforward. The running time upper bound for the function Update(acts) is given in Lemma 8.7. The upper bound on the 2 -error ∆σ − E[∆ σ] 2 follows from Theorem 6.3 and the lemmas in Section 8.6. The upper bound on the variance follows from the lemmas in Section 8.4 and 8.5.

Leverage score's moving
Definition 8.2 (x t , z t and y s,t ). Given vectors w mid , w new ∈ R m + and v mid ∈ R m + . For all t ∈ [0, 1], we define x t to be w mid + t(w new − w mid ).
To define z t for t ∈ [0, 1], we first define z t for all t ∈ T as in Algorithm 5. We then extend the definition of z t to the entire interval [0, 1] by connecting consecutive points with segments.

15:
Set S, T with size N and corresponding weights ω as in Theorem 3.10 16: pm.init(A, w, comp ) 17: end procedure 18: 19: procedure Update(acts ∈ Vector Action ) Update estimate for a sequence of updates 20: Let A and w new be the final matrix and vector in acts 21: w mid , ∆ σ mid ← BatchedUpdate(acts) Let w mid , ∆ σ mid be computed as in Theorem 6.3 (Algorithm 2, Algorithm 3) with the sequence of updates in acts for t in T do 24:

26:
Compute θ i,t θ i,t for all i ∈ [m] as in Lemma 8.8

27:
for s in S do 28:

30:
Compute R αβ,s,t α i,s,t , R αβ,s,t β i,s,t , R η,s η i,s for all i ∈ [m] as in Lemma 8.9

31:
for s in S do

33:
Compute R γ,s,s ,t γ i,s,t and R γ,s,s ,t γ i,s ,t for all i ∈ [m] as in Lemma 8.9  Definition 8.3 (η, θ, α, β and γ). We define Here we make the simplifying assumption that (Z t − X t ) is PSD for every t ∈ [0, 1] and that (W mid − W new ) is PSD. In the case where these might not hold, we can simply break these matrices into positive and negative parts and bound both parts using essentially the same calculations.
Lemma 8.4 (Leverage score's moving). The change in leverage score can be written as Proof. We first rewrite σ(w new ) i − σ(w mid ) i into two terms The lemma then follows immediately from the following Lemma 8.5 and 8.6.
Lemma 8.5 (First term in leverage score's moving). For the first term in Lemma 8.4, we have Proof. We can compute τ (w new ) i − τ (w mid ) i as where the last step follows by the definition of in a similar way as follows: where recall the definition of y s,t = z t + s(x t − z t ). Using the above expressions, we can rewrite Lemma 8.6 (Second term in leverage score's moving). For the second term in Lemma 8.4, we have

49
The second term can be computed as

Running time anlaysis
Lemma 8.7 (Update time). Assume the sequence acts satisfies Assumption 7.1 and comp ≤ 0.01.
Then each call to the function Update(acts) makes O(N ) calls to pm.update and takes an extra O(n 2 · |acts| + T mat (n, n, r) · N 3 · log(n)) time.
Proof. The number of calls to pm.update directly follows from Algorithm 5. The time to compute ∆ σ mid is O(n 2 · |acts|) by applying Theorem 6.3. The remaining part of the running time is given by the following Lemma 8.8, 8.9 and 8.10.
Lemma 8.8 (Computing θ θ). Assume m = O(n). Then for any t ∈ T , the time to compute Proof. Recall from Definition 8.3 that Denote Q = Q(z t ) and Q (2) the entry-wise square of Q. We have Therefore, in order to compute θ i,t θ i,t for all i ∈ [m], it suffices to compute the matrix-vector product Q (2) (w mid − w new ), which takes time O(n 2 ).
Proof. We only describe how to compute R αβ,s,t α i,s,t for all i ∈ [m] in time O(T mat (n, n, r) · log(n)) as the rest of the calculations are similar. Notice that that computing R αβ,s,t α i,s,t for all i ∈ [m] is essentially computing the matrix Assume we know the matrix M (y s,t ) −1 , then we can perform the computation from left to right, and it follows that each matrix multiplication here can be done in time O(T mat (n, n, r)). To remove the assumption that we know M (y s,t ) −1 , we pre-condition on the matrix M (z t ) −1 and apply Lemma 6.1 to compute R αβ,s,t α i,s,t for all i ∈ [m] (up to negligible error) in time O(T mat (n, n, r) · log(n)).
Lemma 8.10 (Computing ∆ σ new ). Assume m = O(n) and comp ≤ 0.01. Define an approximation ∆ σ new i to the change in leverage score to be Then the time to compute ∆ σ new i for all i ∈ [m] (up to negligible error) is O(T mat (n, n, r)·N 3 ·log(n)).
, we can compute the first term for all i ∈ [m] in O(n) time. By Lemma 8.8, the time to compute the second term for all i ∈ [m] is O(n 2 · N ). It follows from Lemma 8.9 that the time to compute the rest of the terms for all i ∈ [m] is O(T mat (n, n, r) · N 3 · log(n)).

Upper bounding η, α, β and γ
From the computation in Section 8.5, the variance of the three terms where we applied JL matrices are bounded by i η i 4 2 , i α i 2 2 β i 2 2 and i γ i 4 2 . Our bounds for these terms are summarized in Table 6.

Quantity
Bound Lemma Lemma 8.13 Table 6: Upper bounds for i η i Proof. We recall the definition of vector η i,s ∈ R m from Definition 8.3 as follows: where in the last step we define diagonal matrices B, C ∈ R m×m as and for simplicity, we use P to denote the projection matrix P (y s,0 ). Since comp ≤ 0.01, it follows that log(y s,0 ) − log(x 0 ) ∞ ≤ 0.01. Proof. We recall the definition of vector α i,s,t ∈ R m and β i,s,t ∈ R m from Definition 8.3 as follows:

It follows that
We therefore can rewrite α i,s,t and β i,s,t as Our assumptions comp ≤ 0.01 and log(w new ) − log(w mid ) ∞ ≤ 0.01 imply that Therefore we can upper bound α i,s,t as where in the last step we define diagonal matrices Thus, we have

Upper bounding
Lemma 8.13. Assume comp ≤ 0.01 and log(w new ) − log(w mid ) ∞ ≤ 0.01. Then for any s ∈ S and t ∈ T , we have Proof. We recall the definition of vector γ i,s,t ∈ R m from Definition 8.3 as follows:

It follows that
where in the last step we define diagonal matrices B, C, D ∈ R m×m , and for simplicity, we use P to denote P (y s,t ). By our assumptions that simp ≤ 0.01 and we have log(y s,t ) − log(w new ) ∞ ≤ 0.1.

54
Therefore we have  Our variance bounds for the terms where we applied JL matrices are summarized in Table 7. These results are consequences of the following Lemma 8.14.

A variance bound for random Gaussian matrices
Lemma 8.14 (JL variance bound). Given vectors x, y ∈ R m . Let R ∈ R r×m be a random Gaussian matrix with R i,j ∼ N (0, 1/r). Then we have Next we prove the bound on the variance. We use R i for each i ∈ [r] to denote the column vector that corresponds to the ith row of R. We have 55 where the last inequality is because x R 1 is a Gaussian random variable with variance x 2 2 /r.

Variance bound for η η, α β and γ γ
The lemmas below all follow immediately from Lemma 8.14.
Lemma 8.15 (Variance bound for η η). Let σ i,dis and σ i,jl be defined as follows where R η,s ∈ R r×m for s ∈ S are independent random Gaussian matrices with (R η,s ) i,j ∼ N (0, 1/r). Then, we have Lemma 8.16 (Variance bound for α β). Let σ i,dis and σ i,jl be defined as follows where R αβ,s,t ∈ R r×m for s ∈ S, t ∈ T are independent random Gaussian matrices with (R αβ,s,t ) i,j ∼ N (0, 1/r). Then, we have Lemma 8.17 (Variance bound for γ γ). Let σ i,dis and σ i,jl be defined as follows where R γ,s,s ,t ∈ R r×m are independent random Gaussian matrices with (R γ,s,s ,t ) i,j ∼ N (0, 1/r).

Error upper bound for discrete sampling
Our error upper bounds for approximating integral by discrete sampling are summarized in Table 8. We only prove Lemma 8.18 in the following. The rest of the lemmas follows from similar arguments.
Lemma 8.18 (Error upper bound for η η). Assume comp ≤ 0.01. Let σ i,cts and σ i,dis be defined as follows where |S| = N . Then we have Proof. Define function f : [0, 1] → R to be f (s) = η i,s η i,s . Applying Theorem 3.9 with N = |S|, the error is bounded as In order to bound M 2N , we need the following Cauchy's estimates.
Theorem 8.19 (Cauchy's Estimates). Suppose f is holomorphic on a neighborhood of the ball B = {z ∈ C : |z − z 0 | ≤ r}, then we have that Since f (x) is a rational polynomial, we can extend the definition of f (s) to the complex plane and the resulting function, which we also denote as f (s). Since log(z 0 ) − log(x 0 ) ∞ ≤ comp ≤ 0.01, M (y s,0 ) is invertible for |s| ≤ 1. Hence f is holomorphic on the unit ball on the complex plane around 0. Applying Theorem 8.19 with r = 1, we have that Therefore, the error is bounded as |σ i,dis − σ i,cts | ≤ poly(n) 2 2N , which finishes the proof of the lemma.
Lemma 8.20 (Error upper bound for θ θ). Assume comp ≤ 0.01. Let σ i,cts and σ i,dis be defined as follows where |T | = N . Then we have Lemma 8.21 (Error upper bound for α β). Assume comp ≤ 0.01. Let σ i,cts and σ i,dis be defined as follows where |S| = |T | = N . Then we have where |S| = |T | = N . Then we have In this section we present an overview of Vaidya's cutting plane method [Vai89a] and illustrate how our leverage score maintenance data structure in Section 5 implies a faster implementation in O(n SO log(κ) + n 3 log(κ)) time. Formally, we prove the following Theorem 4.1 from Section 4.
Theorem 4.1. Given a separation oracle for a convex set K ⊂ R n that is contained in a box of radius R and a parameter > 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 in O((n SO +n 3 ) log(κ)) time, where SO is the complexity of the separation oracle and κ = nR/ .
Our cutting plane method essentially replaces the leverage scores σ in Vaidya's method by estimates σ from our leverage score maintenance data structure in Section 5, which satisfies that σ − σ 2 ≤ 1/ log O(1) (n). To justify the validility of such a replacement, we first give an overview of Vaidya's method.

Feasible region Vaidya's algorithm maintains a polytope P
Volmetric barrier function and volumetric center Vaidya's algorithm maintains an approximate minimizer z (k) ∈ R n , known as the volumetric center, of the volumetric barrier function: x . Here z (k) ∈ R n serves as the query point to the separation oracle.
Leverage score Each constraint i of P (k) is associated with a leverage score which measures its relative importance (see preliminary for the definition). It is well-known that 0 ≤ σ i (z) ≤ 1, ∀i ∈ [m (k) ] and m (k) i=1 σ i (z) = n. We denote by σ and Σ the vector and diagonal matrix of leverage scores respectively.
Updating P (k) according to leverage scores P (k) can be updated in two different ways. Parameters are chosen such that F (z) − min z∈R n F (z) does not increased more than a multiplicative constant.
• Whenever the leverage score σ i (z) is smaller than some universal constant c 1 , constraint i is dropped and z is updated by the Newton method below. As i σ i (z) = n, this implies that the number of constraints is ≤ i σ i c 1 = n c 1 = O(n).
• Otherwise, the separation oracle is queried at the volumetric center z (k) and returns a new separating hyperplane a k x ≥ b k . However, P (k+1) is not the intersection of P (k) and a k x ≥ b k . Instead, a k x ≥ b k is added for some b k ≤ b k so that the leverage score of a k x ≥ b k is 0.5(δc 1 ) 1/2 , where δ ≥ 10 3 c 1 is another small universal constant.
Updating z (k) via Newton-type method Upon adding or removing a constraint the volumetric center z (k) ∈ R n must be recomputed. As z (k) should minimize F (k) (x), we iteratively update z (k) via a Newton-type step: where Q(z) := A S −1 z ΣS −1 z A ∈ R n×n can be shown to be a constant spectral approximation to the Hessian ∇ 2 F (z) ∈ R n×n . One can calculate ∇F (z) = A S −1 z σ ∈ R n . We perform this Newton step iteratively until z new is a good approximate minimizer of F (k) (x), i.e. F (z) − min z∈R n F (z) ≤ c 2 for some small universal constant c 2 . Theorem 1 in full version of [Vai89a] states that Since the decrease is multiplicative, this can be accomplished in only O(1) many iterations.
Performance guarantee of Vaidya's method

Number of iterations required
Vaidya showed that after T iterations, the volume of P (k) decreases by a factor of c T −O(n log(n)) for some constant c, i.e. vol(P (k) ) ≤ c T −O(n log(n)) vol(P (0) ) = c T −O(n log(n)) R n .
Therefore in T = O(n log(nR/ )) iterations, we have vol(P (k) ) ≤ O(n) showing that P (k) does not contain a ball of radius and hence solving the feasibility problem.

Running time per iteration
As i σ i (z) = n and the leverage scores are maintained so that σ i ≥ c 1 always holds, the number of constraints is ≤ i σ i c 1 = n c 1 = O(n). Thus all vectors and matrices above have dimension O(n) and O(n) × O(n). Moreover, recall that only O(1) steps of Newton method are needed within a iteartion of cutting plane. Therefore in one iteration, the running time of Vaidya is O(n 2 ) plus the time to compute σ and to solve a linear system in Q(z) −1 (from the Newton step), which naively requires O(n ω ) time. In the rest of this section we explain speed up these two bottlenecks using our leverage score maintenance data structure.

A.1 Our faster implementation via leverage score maintenance
We provide a faster implementation of Vaidya's method via our leverage score maintenance data structure, which efficiently updates leverage scores. Specifically, We design a data structure which, upon updates to the volumetric center z, maintains an estimate σ of the leverage scores σ in amortized O(n 2 ) time (Theorem 5.1). Our error guarantee satisfies To apply our data structure, we plug in W = S −2 z as the weight in Theorem 5.1. To establish the validity of our method, we show that conditions (1) and (2) required for our data structure are satisfied for sufficiently small parameters c 1 , δ, c 2 . We then prove that Vaidya's performance guarantee is preserved in the presence of a small perturbation to the leverage score.

B Modified Projection Maintenance
Theorem B.1 ([CLS19, Appendix]). Let t k = T mat (m, m, k) denote the time to multiply a m × m and a m × k matrix. Given a matrix A ∈ R m×n with m ≥ n, and k * ∈ [m], there is a deterministic data structure that approximately maintains the projection matrices √ W A(A W A) −1 A √ W and the inverse matrices (A W A) −1 for positive diagonal matrices W through the following operations: • Initialize(A, w, ): Initialize the data structure of the matrix A, the weight w and the target accuracy ∈ (0, 1/4) in m ω+o(1) time.
• Update(w): Output a vector v ∈ R m along with (A V A) −1 and • Insert(a, w a ): Insert a column a into A, a weight w a into w in O(m 2 ) time.
• Delete(a, w a ): Delete a column a from A and its corresponding weight w a from w in O(m 2 ) time.
Suppose that the number of columns is O(m) during the whole algorithm and that for any call of where w is the input of call, w (old) is the weight before the call. Then, the amortized expected time We will use this theorem with difference algorithms for rectangular matrix multiplication. Recall that it takes O(m 2 log 2 m) time to multiply an m × m and an m × m 0.17 matrix. By splitting the matrix into blocks (See e.g. [CLS19, Lemma A.5]), one can check it takes time to multiply a m × m and a m × k matrix with α = 0.17. Using this and putting k * = m α , we have

C Cutting Plane Method for Convex Minimization and Saddle Point Problems
We show in this section that cutting plane methods can be applied to not only convex minimization, but also the more general problem of computing a saddle point in a convex-concave game with essentially the same guarantee.

C.1 Notations and definitions
For a convex set X ⊆ R n , the interior of X , denoted as int(X ), is the subset of points in X that has a small neighborhood fully contained inside X . We denote by vol(X ) the volume of X . For any vector x and r > 0, we use B(x, r) to denote the 2 ball of radius r centered at x, i.e. Theorem C.1. Let f be a convex function on R n and S be a convex set that contains a minimizer of f . Suppose we have a subgradient oracle for f with cost T and S ⊂ B ∞ (0, R). Using B ∞ (0, R) as the initial polytope for our Cutting Plane Method, for any 0 < α < 1, we can compute x ∈ S such that f (x) − min y∈S f (y) ≤ α (max y∈S f (y) − min y∈S f (y)) , with high probability in n and with a running time of O(T · n log(κ) + n 3 log(κ)), where κ = nγ/α and γ = R/minwidth(S). Here the minwidth of S is defined by minwidth(S) := min a∈R n : a 2 =1 max y∈S a y − min y∈S a y .

C.3 Convex-concave games
In this subsection we show that a similar guarantee holds for solving convex-concave games. Much of the materials in this section are modified from [Nem95, Lecture 5]. For completeness, we will explain both the standard theory and the various changes needed for our promised runtime. In the convex-concave game, we are asked to solve where X ∈ R n , Y ∈ R m are convex sets and f (x, y) : X × Y → R is a continuous function that is convex in x ∈ X and concave in y ∈ Y. Von Neumann's minimax theorem states that

C.4 Applying cutting plane method to convex-concave games
In the cutting plane framework, we maintain a polytope P (k) ⊆ R n+m that contains a saddle point (x * , y * ) of f . The algorithm initially starts with the polytope P (0) := B ∞ (0, R) ⊃ X × Y.
In each iteration, we call the first-order oracle at a certain point z (k) = (x (k) , y (k) ) ∈ P (k) depending on P (k) . If z (k) / ∈ X × Y is not feasible, then we obtain a supporting hyperplane to separate z (k) from the feasible region X × Y; and if z (k) ∈ X × Y is feasible, Lemma C.3 shows that any saddle point (x * , y * ) lies in the half-space H (k) := {z ∈ R n+m : (z−z (k) ) g(z (k) ) ≤ 0} and hence (x * , y * ) ∈ P (k) ∩ H (k) . The algorithm continues by choosing P (k+1) which contains P (k) ∩ H (k) .
For simplicity, we extend the gradient function g as follows: where n(x, y) is some vector of length β > L such that X × Y ⊂ {z ∈ R n+m : (z − (x, y)) n(x, y) ≤ 0}. The cutting plane framework is given in Algorithm 6.

10:
Construct the next z (k+1) , P (k+1) according to cutting plane method 11: end for 12: return a point according to Lemma C.7 13: end procedure Theorem C.4. Let c < 1 be a universal constant. For T ≥ n + m, we can find z (k) ∈ R n+m and P (k) ⊆ R n+m in total time O((n + m) 2 T ) with high probability in n + m such that vol(P (T ) ) ≤ c T −O(n+m) log(n+m) vol(P (0) ).
Proof. By Lemma C.3, any saddle point lies in H (k) showing that it is indeed a separating hyperplane. With this separation oracle, we apply our faster implementation of Vaidya's cutting plane method from section 4.
The second part of the lemma follows from the fact that this method always maintains a polytope with O(n) constraints. For the first part, note that this just paraphases Lemma A.4.

C.5 Cutting plane method for convex-concave games: generating solutions
Recall that for minimizing a convex function f (z), one can simply output the best z (k) found in all iterations. The argument here is that as long as the volume vol(P (T ) ) is small enough, then in some iteration k, a point close to the optimal solution z * gets removed from P (k) which indicates that f (z (k) ) is also close to optimal. For the convex-concave game, however, such a naive approach of outputting the "best" z (k) would fail as illustrated by [Nem95, Section 5.3]. The crucial difference here in the convex-concave game is that although we have an objective function f , but we are not interested in optimizing it. Instead, we are interested in finding an -saddle point (x, y) that satisfies (14).
The correct idea is to output some convex combination of all z (k) 's that does converge to the saddle point (x * , y * ) ∈ X × Y which we describe in the following.
Assume that we have performed T steps of the cutting plane method.
Proof. We run our cutting plane method for T = (n + m) log n+m R r iterations. By Theorem C.4, we obtain P (T ) with volume vol(P (T ) ) ≤ n + m r R
Next we give a slightly refined runtime in the case where the domain is a ball using a result of Nemirovski [Nem95].
Corollary C.10. Given an L-Lipschitz function f (x, y) : X × Y → R that is convex in x and concave in y with X = B(0, R) ⊂ R n and Y = B(0, R) ⊂ R m . For any 0 < ≤ 1 2 , we can find ( x, y) such that By using the first result when ≥ 1 n+m and using the second result when ≤ 1 n+m , we have the promised runtime.

D.1 Linear Arrow-Debreu markets
Previous work The celebrated result of Arrow and Debreu [AD54] shows the existence of a market equilibrium for a broad class of utility functions. From the computational aspects, computing the equilibrium of a market equilibrium in the case of linear utility functions admits a long line of research (see Table 9), leading to both weakly and strongly polynomial runtimes.
A convex formulation The problem of computing a market equilibrium in linear exchange markets can be formulated as the following convex program due to [DGV16], with variables p i representing the prices, β i the inverse best bang-per-bucks, and y ij the money paid by agent i to agent j. Given u ∈ R n×n + as an input matrix, the goal is to solve the following minimization problem, min p∈R n ,β∈R n ,y∈R n×n i∈ [n] p i log(p i /β i ) − Convex-concave game formulation The number of variables in the above convex program can be as large as n 2 and thus a naive application of the cutting plane method would lead to lead to a large runtime. To solve the convex program more efficiently, we transform it into a convex-concave game with a reduced number of variables. Let λ, η ∈ R n be the Lagrange multiplier of constraint 19 and 20 respectively. We have that the above convex program is equivalent to: min p≥1∈R n ,β∈R n + ,u ij β i ≤p j min y∈R n×n which is equivalent to min p≥1,β≥0,u ij β i ≤p j max λ,η∈R n min . Assume all utilities are integers ≤ U and we let ∆ = (nU ) n . Then there exists equilibrium prices p that are quotient of two integers ≤ ∆, along with allocations x that are quotients of two integers ≤ ∆ 2 .
Now using Theorem C.9, we need O(n 2 log(nU )) iterations and the first order oracle requires O(m) operations. Therefore, the total number of operations is O(mn 2 log(nU )). We formally state our result as follows.
Theorem D.5 (Formal version of Theorem 1.4). There exists a weakly polynomial algorithm that computes a market equilibrium in linear exchange markets in time O(mn 2 log(nU )).

D.2 Fisher markets with spending constraint utilities
Market equilibrium An allocation x ∈ R n B ×n G ×L + and price vector p ∈ R n G + are an equilibrium if two conditions are satisfied. The first condition, buyer optimality, is that given p, each player maximizes his utility subject to his budget constraint. In other words, for all i ∈ [n B ], the allocation x optimizes the following program (where p is a constant): The second equilibrium condition, market clearance, is that i,l x i,j,l = 1 for all j ∈ [n G ].
Previous work The Fisher market with spending constraint utilities problem was introduced in [Vaz10] where they gave a weakly poly time algorithm that takes O(n 3 (n + m) 2 log U ) max-flow computations, where U = max i∈[n B ],j∈[n G ],l∈[L] u i,j,l . A convex formulation of the problem was first given in the unpublished manuscript [BDX10]. A strongly poly time algorithm (which is also the best weakly polynomial running time) was given in [Vég16]which achieves O(mn 3 +m 2 (m+n log n) log m) number of operations.
A convex formulation The problem of computing a market equilibrium in Fisher markets with spending constraint utilities problem can be captured by the following convex program due to [BDX10]. Theorem D.6 ( [BDX10]). An optimum solution to the above convex program corresponds to an equilibrium for the Fisher market with spending constraint utilities with allocation given by x i,j,l = b i,j,l /p j .
Convex-concave game formulation We transform the convex program above to a convexconcave saddle point which allows us to apply the cutting plane method. Let η j , λ i be the Lagrange multipliers for the two equality constraints and µ i,j,l ≥ 0 be the Lagrange multipliers for the inequalities b i,j,l ≤ B i,j,l in the above convex formulation. We have the above convex program is equivalent to max p≥0 max b≥0 min µ≥0,η,λ − n G j=1 p j log p j + i,j,l b i,j,l log u i,j,l where the last step is because B i,j,l ≥ 0. Now applying Theorem C.9, we need O(n 2 log(nU )) iterations and the first order oracle takes O(m) operations. So the total number of operations is O(mn 2 log(nU )). We formally state our result as follows: Theorem D.7 (Formal version of Theorem 1.5). 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 )).

D.3 Walrasian equilibrium for general buyer valuations and fixed supply
They showed that such equilibria are characterized by the minima of the following convex program, which can be solved in polynomial time via a suitable extension of the cutting plane methods. Let [m] and [n] be the set of buyers and goods respectively. The setting of the problem is as follows.
Buyers' valuation and utility Each buyer i ∈ [m] may have an arbitrary valuation v i : N n → R over goods, i.e. for x ∈ N n , v i (x) is the valuation of i if she gets x j units of good j ∈ [n]. Given prices p ∈ R [n] on goods, the utility of i getting x ∈ N n is given by v i (x) − p · x. To maximize her utility, given p, i would buy a bundle Aggregate demand oracle To study this problem computationally, one must specify the information about the market available to the algorithm. Paes Leme and Wong showed that this problem can be solved in polynomial time assuming the aggregate demand oracle. Given prices p, the aggregate demand is given by i∈ [m] x (i) , where is a utility-maximizing bundle for buyer i.
Walrasian equilibrium In this economy the supply of goods is fixed at s ∈ N n . Prices p are said to be at equilibrium if there exist some aggregate demand that matches the supply, i.e. s = i∈[m] x (i) .
Convex programming One approach to computing Walrasian equilibrium is by formalating the problem as the following convex program. Previous works applied subgradient descent methods to this program to obtain pseudopolynomial time algorithms [Par99, PU02, AM02]. In contrast, Paes Leme and Wong showed that this program can in fact be solved in weakly polynomial time.
Here p denotes the prices and u i can be thought of as the utility of buyer i, ∀i ∈ [m]. This convex program can clearly be solved with access to the demand oracle of each buyer i, i.e. given p, return arg max x v i (x) − p · x, which serves as a separation oracle.
The main result of Paes Leme and Wong is that this can be done via the weaker aggregate demand oracle, i.e. given p, return i x (i) , where Their algorithm is based on a suitable extension of standard cutting plane methods: Theorem D.8 ([LW17]). There is an algorithm that runs in time O(n 2 T AD log(SM n) + n 6 log O(1) (SM n)) to compute Walrasian prices whenever it exists using only access to an aggregate demand oracle. Here T AD denotes the runtime of the aggregate demand oracle, n denotes the number of goods, M = max i∈[m],0≤x≤s |v i (x)| and S = max j∈[n] s j (s and v i are integer-valued).
As in previous applications, n 2 iterations of cutting plane are required as has to be taken to be exponentially small which results in an overall runtime overhead of O(n 6 ). By leveraging our faster cutting plane method, we obtain the following improved result: Theorem D.9 (Formal version of Theorem 1.6). There is an algorithm that runs in time O(n 2 T AD log(SM n) + n 4 log(SM n)) to compute Walrasian prices whenever it exists using only access to an aggregate demand oracle.
Proof. Same as Paes Leme-Wong except that we invoke our cutting plane method for convex minimization in Theorem C.1 instead of LSW's.