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

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)$.

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.

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.

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.

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

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

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:

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:

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:

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.

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.

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.

4 Overview of our approach

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.

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

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

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.

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:

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

6 Related works

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 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].

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$.

Our Techniques

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

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.

2 Layered data structure

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.

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.

3 Batched low-rank update

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).

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.

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

In order to approximate such an integral, we take a set $\mathcal{T}\subseteq$ 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

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

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

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

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

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.

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

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

Preliminaries

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

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 square full-rank matrix $A$, we use $A^{-1}$ to denote the inverse of $A$.

2 Operators

3 Different types of running 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).

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.

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.

4 Basic results on matrices

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

5 Fast matrix multiplication

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.

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

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

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

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

6 Multiple variable polynomial interpolation

We first state a one-variable version interpolation theorem.

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

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

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}$.

where the last step follows by Claim 3.11. ∎

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

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

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$.

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.

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.

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.

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

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

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

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

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

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

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})$.

Batched 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.

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

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

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

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

The result follows from some rearranging. ∎

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

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$.

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

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

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

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

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

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

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

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

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

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

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

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

(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

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

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)

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

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)}$

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

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

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

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

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).

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

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)$.

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)$.

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

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)$.

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:

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)}$ .

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

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

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

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

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.

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.

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. ∎

2 Approximate leverage score’s moving

We use $\Delta_{1}$ and $\Delta_{2}$ to denote the two corresponding error terms

Assume Assumption 7.1 holds and that $\epsilon_{\operatorname{simp}}\leq 0.01$. Then we have the following upper bound on $\left\lVert\Delta_{1}\right\rVert_{2}$

Recall the definition of the error term $\Delta_{1}$ as

and we use $P$ to denote $P(y_{\xi})$. Then we have

where the last step follows from Theorem 6.3. ∎

Assume Assumption 7.1 holds and that $\epsilon_{\operatorname{simp}}\leq 0.01$. Then we have the following upper bound on $\left\lVert\Delta_{2}\right\rVert_{2}$

Recall the definition of the error term $\Delta_{2}$ as

and we use $P$ to denote $P(v_{t,s})$. By Assumption 7.1, we have

This together with the assumption that $\epsilon_{\operatorname{simp}}\leq 0.01$ further implies that

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.

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

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

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

2 Leverage score’s moving

To define $z_{t}$ for $t\in$, we first define $z_{t}$ for all $t\in{\cal T}$ as in Algorithm 5. We then extend the definition of $z_{t}$ to the entire interval $$ by connecting consecutive points with segments.

For any $s\in$, we define $y_{s,t}=z_{t}+s(x_{t}-z_{t})$.

The change in leverage score can be written as

The lemma then follows immediately from the following Lemma 8.5 and 8.6. ∎

For the second term in Lemma 8.4, we have

3 Running time anlaysis

Assume the sequence $\mathsf{acts}$ satisfies Assumption 7.1 and $\epsilon_{\operatorname{comp}}\leq 0.01$. Then each call to the function $\textsc{Update}(\mathsf{acts})$ makes $O(N)$ calls to $\textsc{pm}.\textsc{update}$ and takes an extra $O(n^{2}\cdot|\mathsf{acts}|+\mathcal{T}_{\operatorname{mat}}(n,n,r)\cdot N^{3}\cdot\log(n))$ time.

Assume $m=O(n)$. Then for any $t\in{\cal T}$, the time to compute $\theta_{i,t}^{\top}\theta_{i,t}$ for all $i\in[m]$ is $O(n^{2})$.

Denote $Q=Q(z_{t})$ and $Q^{(2)}$ the entry-wise square of $Q$. We have

We only describe how to compute $R_{\alpha\beta,s,t}\alpha_{i,s,t}$ for all $i\in[m]$ in time $O(\mathcal{T}_{\operatorname{mat}}(n,n,r)\cdot\log(n))$ as the rest of the calculations are similar. Notice that that computing $R_{\alpha\beta,s,t}\alpha_{i,s,t}$ for all $i\in[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(\mathcal{T}_{\operatorname{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_{\alpha\beta,s,t}\alpha_{i,s,t}$ for all $i\in[m]$ (up to negligible error) in time $O(\mathcal{T}_{\operatorname{mat}}(n,n,r)\cdot\log(n))$.

Assume $m=O(n)$ and $\epsilon_{\operatorname{comp}}\leq 0.01$. Define an approximation $\Delta\widetilde{\sigma}^{\operatorname{new}}_{i}$ to the change in leverage score to be

Then the time to compute $\Delta\widetilde{\sigma}^{\operatorname{new}}_{i}$ for all $i\in[m]$ (up to negligible error) is $O(\mathcal{T}_{\operatorname{mat}}(n,n,r)\cdot N^{3}\cdot\log(n))$.

4 Upper bounding η𝜂\eta, α𝛼\alpha, β𝛽\beta and γ𝛾\gamma

From the computation in Section 8.5, the variance of the three terms where we applied JL matrices are bounded by $\sum_{i}\left\lVert\eta_{i}\right\rVert_{2}^{4}$, $\sum_{i}\left\lVert\alpha_{i}\right\rVert_{2}^{2}\left\lVert\beta_{i}\right\rVert_{2}^{2}$ and $\sum_{i}\left\lVert\gamma_{i}\right\rVert_{2}^{4}$. Our bounds for these terms are summarized in Table 6.

Assume $\epsilon_{\operatorname{comp}}\leq 0.01$, where $\epsilon_{\operatorname{comp}}$ is the error parameter in Algorithm 5. Then for any $s\in{\cal S}$, we have

where $y_{s,t}=z_{t}+s(x_{t}-z_{t})$ and $Y_{s,t}=Z_{t}+s(X_{t}-Z_{t})$. Recall from Definition 3.2 that $P(v)=V^{1/2}A(A^{\top}VA)^{-1}A^{\top}V^{1/2}=V^{1/2}\cdot Q(v)\cdot V^{1/2}$. We can rewrite $\eta_{i,s}$ as

and for simplicity, we use $P$ to denote the projection matrix $P(y_{s,0})$. Since $\epsilon_{\operatorname{comp}}\leq 0.01$, it follows that

where recall that $y_{s,t}=z_{t}+s(x_{t}-z_{t})$ and $Y_{s,t}=Z_{t}+s(X_{t}-Z_{t})$. Recall from Definition 3.2 that $P(v)=V^{1/2}A(A^{\top}VA)^{-1}A^{\top}V^{1/2}=V^{1/2}\cdot Q(v)\cdot V^{1/2}$. We therefore can rewrite $\alpha_{i,s,t}$ and $\beta_{i,s,t}$ as

Therefore we can upper bound $\alpha_{i,s,t}$ as

where in the last step we define diagonal matrices

where recall that $y_{s,t}=(z_{t}+s(x_{t}-z_{t}))$ and $Y_{s,t}=(Z_{t}+s(X_{t}-Z_{t}))$. Recall from Definition 3.2 that $P(v)=V^{1/2}A(A^{\top}VA)^{-1}A^{\top}V^{1/2}=V^{1/2}\cdot Q(v)\cdot V^{1/2}$. We therefore can rewrite $\gamma_{i,s,t}$ as

and for simplicity, we use $P$ to denote $P(y_{s,t})$. By our assumptions that $\epsilon_{\operatorname{simp}}\leq 0.01$ and

5 Variance upper bound for random Gaussian matrices

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.

Since $\operatorname*{{\bf{E}}}[R^{\top}R]=I$, we have

Next we prove the bound on the variance. We use $R_{i}$ for each $i\in[r]$ to denote the column vector that corresponds to the $i$th row of $R$. We have

where the last inequality is because $x^{\top}R_{1}$ is a Gaussian random variable with variance $\left\lVert x\right\rVert_{2}^{2}/r$.

The lemmas below all follow immediately from Lemma 8.14.

Let $\sigma_{i,\operatorname{dis}}$ and $\sigma_{i,\operatorname{jl}}$ be defined as follows

Let $\sigma_{i,\operatorname{dis}}$ and $\sigma_{i,\operatorname{jl}}$ be defined as follows

Let $\sigma_{i,\operatorname{dis}}$ and $\sigma_{i,\operatorname{jl}}$ be defined as follows

6 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.

Assume $\epsilon_{\operatorname{comp}}\leq 0.01$. Let $\sigma_{i,\operatorname{cts}}$ and $\sigma_{i,\operatorname{dis}}$ be defined as follows

In order to bound $M_{2N}$, we need the following Cauchy’s estimates.

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 $\left\lVert\log(z_{0})-\log(x_{0})\right\rVert_{\infty}\leq\epsilon_{\operatorname{comp}}\leq 0.01$, $M(y_{s,0})$ is invertible for $|s|\leq 1$. Hence $f$ is holomorphic on the unit ball on the complex plane around . Applying Theorem 8.19 with $r=1$, we have that

Assume $\epsilon_{\operatorname{comp}}\leq 0.01$. Let $\sigma_{i,\operatorname{cts}}$ and $\sigma_{i,\operatorname{dis}}$ be defined as follows

Assume $\epsilon_{\operatorname{comp}}\leq 0.01$. Let $\sigma_{i,\operatorname{cts}}$ and $\sigma_{i,\operatorname{dis}}$ be defined as follows

where $|{\cal S}|=|{\cal T}|=N$. Then we have

Assume $\epsilon_{\operatorname{comp}}\leq 0.01$. Let $\sigma_{i,\operatorname{cts}}$ and $\sigma_{i,\operatorname{dis}}$ be defined as follows

where $|{\cal S}|=|{\cal T}|=N$. Then we have

Appendix

Appendix A Perturbed Volumetric Center Cutting Plane Method

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\operatorname{SO}\log(\kappa)+n^{3}\log(\kappa))$ time. Formally, we prove the following Theorem 4.1 from Section 4.

Our cutting plane method essentially replaces the leverage scores $\sigma$ in Vaidya’s method by estimates $\widetilde{\sigma}$ from our leverage score maintenance data structure in Section 5, which satisfies that $\|\widetilde{\sigma}-\sigma\|_{2}\leq 1/\log^{O(1)}(n)$. To justify the validility of such a replacement, we first give an overview of Vaidya’s method.

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\leq\sigma_{i}(z)\leq 1$, $\forall i\in[m^{(k)}]$ and $\sum_{i=1}^{m^{(k)}}\sigma_{i}(z)=n$. We denote by $\sigma$ and $\Sigma$ the vector and diagonal matrix of leverage scores respectively.

Whenever the leverage score $\sigma_{i}(z)$ is smaller than some universal constant $c_{1}$, constraint $i$ is dropped and $z$ is updated by the Newton method below. As $\sum_{i}\sigma_{i}(z)=n$, this implies that the number of constraints is $\leq\frac{\sum_{i}\sigma_{i}}{c_{1}}=\frac{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}^{\top}x\geq b_{k}$. However, $P^{(k+1)}$ is not the intersection of $P^{(k)}$ and $a_{k}^{\top}x\geq b_{k}$. Instead, $a_{k}^{\top}x\geq b_{k}^{\prime}$ is added for some $b_{k}^{\prime}\leq b_{k}$ so that the leverage score of $a_{k}^{\top}x\geq b_{k}^{\prime}$ is $0.5(\delta c_{1})^{1/2}$, where $\delta\geq 10^{3}c_{1}$ is another small universal constant.

Since the decrease is multiplicative, this can be accomplished in only $O(1)$ many iterations.

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.

Therefore in $T=O(n\log(nR/\epsilon))$ iterations, we have $\mathsf{vol}(P^{(k)})\leq\epsilon^{O(n)}$ showing that $P^{(k)}$ does not contain a ball of radius $\epsilon$ and hence solving the feasibility problem.

As $\sum_{i}\sigma_{i}(z)=n$ and the leverage scores are maintained so that $\sigma_{i}\geq c_{1}$ always holds, the number of constraints is $\leq\frac{\sum_{i}\sigma_{i}}{c_{1}}=\frac{n}{c_{1}}=O(n)$. Thus all vectors and matrices above have dimension $O(n)$ and $O(n)\times 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 $\sigma$ and to solve a linear system in $Q(z)^{-1}$ (from the Newton step), which naively requires $O(n^{\omega})$ 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 $\widetilde{\sigma}$ of the leverage scores $\sigma$ in amortized $O(n^{2})$ time (Theorem 5.1). Our error guarantee satisfies

To apply our data structure, we plug in $W=S_{z}^{-2}$ 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},\delta,c_{2}$. We then prove that Vaidya’s performance guarantee is preserved in the presence of a small perturbation to the leverage score.

For any constraint $a^{\top}x\geq b$ added or removed, let $s=a^{\top}z-b$ be its slack. We have

Let $H(z)=A^{\top}S_{z}^{-2}A$. Our goal is to show

Recall that a constraint is removed when its leverage score is smaller than $c_{1}$ and added so that its leverage score is $(\delta c_{1})^{1/2}$. In Vaidya’s analysis, the only requirement on $c_{1}$ and $\delta$ is that $c_{1},\delta$ are sufficiently small constants and $\delta\geq 10^{3}c_{1}$. Hence in either case, we can make the leverage score of $a^{\top}x\geq b$ smaller than 0.01 by choosing $c_{1}$ and $\delta$ small enough, i.e. the leverage score of $a^{\top}x\geq b$ satisfies

Since $H(z)$ is PSD and the square root of a PSD matrix exists,

Note that the spectral norm of $(H(z)^{-1/2}a)(H(z)^{-1/2}a)^{\top}$ is $(H(z)^{-1/2}a)^{\top}(H(z)^{-1/2}a)$. Thus