Efficient Inverse Maintenance
and Faster Algorithms for Linear Programming

arXiv:1503.01752v3 [] 14 Oct 2015

Yin Tat Lee
MIT
yintat@mit.edu

Aaron Sidford
MIT
sidford@mit.edu

Abstract
In this paper, we consider the following inverse maintenance problem: given A ∈ Rn×d and a
number of rounds r, at round k, we receive a n×n diagonal matrix D(k) and we wish to maintain
an efficient linear system solver for AT D(k) A under the assumption D(k) does not change too
rapidly. This inverse maintenance problem is the computational bottleneck in solving multiple
optimization problems. We show how to solve this problem with Õ (nnz(A) + dω ) preprocessing
time and amortized Õ(nnz(A) + d2 ) time per round, improving upon previous running times.
Consequently, we obtain the fastest known running times for solving multiple problems
including, linear programming and computing a rounding of a polytope. In particular given
a feasible point in a linear program with n variables, d constraints, and constraint
matrix
√
A ∈ Rd×n , we show how to solve the linear program in time Õ((nnz(A) + d2 ) d log(−1 )).
We achieve our results through a novel combination of classic numerical techniques of low rank
update, preconditioning, and fast matrix multiplication as well as recent work on subspace
embeddings and spectral sparsification that we hope will be of independent interest.

1

Introduction

Solving a sequence of linear systems is the computational bottleneck in many state-of-the-art optimization algorithms, including interior point methods for linear programming [12, 29, 30, 19], the
Dikin walk for sampling a point in a polytope [9], and multiplicative weight algorithms for grouped
least squares problem [1], etc. In full generality, any particular iteration of these algorithms could
require solving an arbitrary positive definite (PD) linear system. However, the PD matrices involved
in these algorithms do not change too much between iterations and therefore the average cost per
iteration can possibly be improved by maintaining an approximate inverse for the matrices involved.
This insight has been leveraged extensively in the field of interior point methods for linear
programming. Combining this insight with classic numerical machinery including preconditioning,
low ranks update, and fast matrix multiplication has lead to multiple non-trivial improvements to
the state-of-the art running time for linear programming [28, 29, 33, 19, 12]. Prior to our previous
work [19], the fastest algorithm for solving a general linear program with d variables, n constraints,
and constraint matrix A ∈ Rn×d , i.e. solving minA~x≥~b ~cT ~x, depended intricately on the precise
ratio of d, n, and nnz(A) the number of non-zero entries in A (See In Figure 1.1). While, these
running times were recently
achieved was still a
√ improved by our work in [19], the running time we √
complicated bound of Õ( nβ(d2 + ndω−1 r−1 + β −ω r2ω + β −(ω−1) drω )) and Õ( d(nnz(A) + dω ))1
to compute an -approximate solution where β ∈ [ nd , 1] and r ≥ 1 are free parameters to be tuned
and ω < 2.3729 is the matrix multiplication constant [36, 7].
1

In this paper we use Õ to hide factors polylogarithmic in d, n and the width U . See Sec 7 for the definition of U .

1

def

z = nnz(A)
√

nz+

√

z = nd
nd2 + n1.34 d1.15
√

n(z + d2.38 )
d(z + d2.38 )

√

n z + nd1.38 + n0.69 d2

≈ d1.32

d2

n

Figure 1.1: The figure shows the fastest algorithms for solving a linear program, minA~x≥~b ~cT ~x,
before [19]. The horizontal axis is the number of constraints n written as a function of d, the
vertical access is the number of nonzero entries in A denoted by z. Each region of the diagram is
the previous best running time for obtaining one bit of precision in solving a linear program with
the given parameters. For more information on how these times were computed, see Appendix B.
In this paper we further improve upon these results. We cast this computational bottleneck
of solving a sequence of slowly changing linear systems as a self-contained problem we call the
inverse maintenance problem and improve upon the previous best√running times for solving this
problem. This yields an improved running time of Õ((nnz(A) + d2 ) d log(−1 )) for solving a linear
program and improves upon the running time of various problems including multicommodity flow
and computing rounding of an ellipse. We achieve our result by a novel application of recent
numerical machinery such as subspace embeddings and spectral sparsification to classic techniques
for solving this problem which we hope will be of independent interest.

1.1

The Inverse Maintenance Problem

The computational bottleneck we address in this paper is as follows. We are given a fixed matrix
A ∈ Rn×d and an number of rounds r. In each round k ∈ [r] we are given a non-negative vector
d~(k) ∈ Rn≥0 and we wish to maintain access to an approximate inverse for the matrix (AT D(k) A)−1
(k)

(k)

where D(k) ∈ Rn×n is a diagonal matrix with Dii = di . We call this problem the inverse
maintenance problem and wish to keep the cost of both maintaining the approximate inverse and
applying it as low as possible under certain stability assumptions on D(k) .
For our applications we do not require that the approximate inverse to be stored explicitly.
Rather, we simply need an algorithm to solve the linear system AT D(k) A~x = ~b efficiently for any
input ~b ∈ Rd . We formally define the requirements of this solver and define the inverse maintenance
problem below.2
Definition 1 (Linear System Solver). An algorithm S is a T -time solver of a PD matrix M ∈ Rd×d
if for all ~b ∈ Rd and  ∈ (0, 1/2], the algorithm outputs a vector S(~b, ) ∈ Rd in time O(T log(−1 ))
2
2
such that with high probability in n, S(~b, ) − M−1~b M ≤  M−1~b M . We call the algorithm S
linear if S(~b, ) = Q~b for some Q ∈ Rd×d that depends only on M and .
2

In Appendix 6 we show that the problem of maintaining a linear solver and maintaining an approximate inverse
are nearly equivalent.

2

Definition 2 (Inverse Maintenance Problem). We are given a matrix A ∈ Rn×d and a number of
rounds r > 0. In each round k ∈ [r] we receive d~(k) ∈ Rn>0 and wish to find a T -time solver S(k) for
PD AT D(k) A such that both T and the cost of constructing S(k) is small. (See Algorithm 1.)
Algorithm 1: Inverse Maintenance Problem
Input: Full rank matrix A ∈ Rn×d , initial scaling d~(0) ∈ Rn≥0 , the number of rounds r > 0.
for Each round k ∈ [r] do
Input: Current point d~(k) ∈ Rn>0
(k)
(k)
Output: Solver S(k) for AT D(k) A where D(k) ∈ Rn×n is diagonal with Dii = di
end
The inverse maintenance problem directly encompasses a computational bottleneck in many
interior point methods for linear programming [17, 18, 30] as well as similar methods for sampling
a point in a polytope and computing approximate John Ellipsoids. (See Section 7 and 8.) Without
further assumptions, we would be forced to pay the cost of solving a general system, i.e. Õ(nnz(A)+
dω ) [26, 21, 3] where ω < 2.3729 is the matrix multiplication constant [36, 7]. However, in these
algorithms the D(k) cannot change arbitrarily and we develop algorithms to exploit this property.
We consider two different assumptions on how D(k) can change. The first assumption, we call
the `2 stability assumption, and is met in most short-step interior point methods and many of our
applications. Under this assumption, the multiplicative change from d~(k) to d~(k+1) is bounded in `2
norm, indicating that multiplicative changes to coordinates happen infrequently on average. The
second assumption, we call the σ stability assumption, is weaker and holds both for these algorithms
as well as a new interior point method we provided recently in [19] to improve the running time
for linear programming. Under this assumption, the multiplicative change from d~(k) to d~(k+1) is
bounded in a norm induced by the leverage scores, ~σA (d~(k) ), a common quantity used to measure
the importance of rows of a matrix (See Section 2.3). Here, multiplicative changes to these important
rows happen infrequently on average (See Section 4).
Definition 3 (`2 Stability Assumption). The inverse maintenance problem satisfies the `2 stability assumption if for all k ∈ [r] we have log(d~(k) ) − log(d~(k−1) ) 2 ≤ 0.1 and β −1 AT D(0) A 
AT D(k) A  βAT D(0) A for β = poly(n).
Definition 4 (σ Stability Assumption). We say that the inverse maintenance problem satisfies
the σ stability assumption if for each k ∈ [r] we have log(d~(k) ) − log(d~(k−1) ) ~σ (d~(k) ) ≤ 0.1,
A
log(d~(k) ) − log(d~(k−1) )
≤ 0.1, and β −1 AT D(0) A  AT D(k) A  βAT D(0) A for β = poly(n).
∞

Note that many inverse maintenance algorithms do not require the condition β −1 AT D(0) A 
AT D(k) A  βAT D(0) A; However, this is a mild technical assumption which holds for all applications we are interested in and it allows our algorithms to achieve further running time improvements.
Also note that the constant 0.1 in these definitions is arbitrary. If d~(k) changes more quickly
than this, then in many cases we can simply add intermediate d~(k) to make the conditions hold and
only changes the number of rounds by a constant.
Finally, note that the σ stability assumption is strictly weaker the `2 stability assumption as
σi ≤ 1 (See Section 2.3). We use `2 assumption mainly for comparison to previous works and as
a warm up case for our paper. In this paper, our central goal is to provide efficient algorithms for
solving the inverse maintenance under the weaker σ stability assumption.
3

Year

Author

Amortized Cost Per Iteration

Best known linear system solver [26, 21, 3]
1984 Karmarkar [12]
1989 Nesterov and Nemirovskii [27, 28, Thm 8.4.1]
1989 Vaidya [33]
2014 Lee and Sidford [17]
2015 This paper

Õ(nnz(A) + d2.3729 )
Õ(nnz(A) + n0.5 d1.3729 + n0.1865 d2 )
Õ(n0.5 d1.3729 + n0.1718 d1.9216 + n1.8987 )
O(nnz(A) + d2 + n0.8334 d1.1441 )
Õ(nnz(A) + d2 + n0.8260 d1.1340 )
Õ(nnz(A) + d2 )

Figure 1.2: Algorithms for solving the inverse maintenance problem under `2 stability guarantee.
To simplify comparison we let the amortized cost per iteration denote the worst of the average per
iteration cost in maintaining the T -time solver and T .

1.2

Previous Work

Work on the inverse maintenance problem dates back to the dawn of interior point methods, a broad
class of iterative methods for optimization. In 1984, in the first proof of an interior point method
solving a linear program in polynomial time, Karmarkar [12] observed that his algorithm needed to
solve the inverse maintenance problem under the `2 stability guarantee. Under this assumption if
one entry of d~(k) changes by a constant then the rest barely change at all. Therefore the updates
to AT D(k) A are essentially low rank, i.e. on average they are well approximated by the addition
of a rank 1 matrix. He noted that it sufficed to maintain an approximation (AT D(k) A)−1 and
by using explicit formulas for rank 1 updates he achieved an average cost of Õ(n(ω−2)/2 d2 ) for the
inverse maintenance problem. This improved upon Õ(ndω−1 ), the best running time for solving the
necessary linear system known at that time.
Nesterov and Nemirovskii [27] built on the ideas of Karmarkar. They showed that how to maintain an even cruder approximation to (AT DA)−1 and use it as a preconditioner for the conjugate
gradient method to solve the necessary linear systems. By using the same rank 1 update formulas
as Karmarkar and analyzing the quality of these preconditioned systems they were able to improve
the running time for solving the inverse maintenance problem.
Vaidya [33] noticed that instead of maintaining an approximation to (AT DA)−1 explicitly one
can maintain the inverse implicitly and group the updates together into blocks. Using more general
low-rank update formulas and fast matrix multiplication he achieved an improved cost of O(nnz(A)+
d2 + (nd(ω−1) )5/6 ). His analysis was tailored to the case where the matrix A is dense; his running
time is O(nd) which is essentially optimal when nnz(A) = nd. Focusing on the sparse regime, we
showed that these terms that were “lower order” in his analysis can be improved by a more careful
application of matrix multiplication [19].
Note that, for a broad setting of parameters the previous fastest algorithm for solving the inverse
maintenance problem was achieved by solving each linear system from scratch using fast regression
algorithms [26, 21, 3] These algorithms all provide an Õ(nnz(A) + dω )-time solver for AT D(k) A
directly by using various forms of dimension reduction. The algorithm in [26] directly projected
the matrix to a smaller matrix and the algorithms in [21, 3] each provide iterative algorithms for
sampling rows from the matrix to produce a smaller matrix.
Also note that we are unaware of a previous algorithm working directly with the σ stability
assumption. Under this assumption it is possible that many D(k) change by a multiplicative constant
in a single round and therefore updates to AT D(k) A are no longer “low rank” on average, making
low rank update techniques difficult to apply cheaply. This is not just an artifact of an overly weak
assumption; our linear programming algorithm in [19] had this same complication. In [19] rather

4

than reasoning about the σ stability assumption we simply showed how to trade-off the `2 stability
of the linear systems with the convergence of our algorithm to achieve the previous best running
times.

1.3

Our Approach

In this paper we show how to solve the inverse maintenance under both the `2 stability assumption
and the weaker σ stability assumption such that the amortized maintenance cost per iteration
Õ(nnz(A) + d2 ) and the solver runs in time Õ((nnz(A) + d2 ) log(−1 )).
To achieve our results, we show how to use low rank updates to maintain a spectral sparsifier for
AT D(k) A, that is a weighted subset of Õ(d) rows of A that well approximate AT D(k) A. We use the
well known fact that sampling the rows by leverage scores (see Section 2.3) provides precisely such
a guarantee and show that we can decrease both the number of updates that need to be performed
as well as the cost of those updates.
There are two difficulties that need to be overcome in this approach. The first difficulty is
achieving any running time that improves upon both the low rank update techniques and the
dimension reduction techniques; even obtaining a running time of Õ(nnz(A) + dω−c ) for the current
ω and c > 0 for the inverse maintenance problem under the `2 stability assumption was not known.
Simply performing rank 1 updates in each iteration is prohibitively expensive as the cost of such an
update is O(d2 ) and there would be Ω(dc ) such updates that need to be performed on average in
each iteration for some c > 0. Furthermore, while Vaidya [33] overcame this problem by grouping
the updates together and using fast matrix multiplication, his approach needed to preprocess the
rows of A by exactly computing (AT D(0) A)−1 A. This takes time O(ndω−1 ) and is prohibitively
expensive for our purposes if n is much larger than d.
To overcome this difficulty, we first compute an initial sparsifier of AT D(0) A using only Õ(d)
rows of A in time Õ(nnz(A)+dω ) [26, 21, 3]. Having performed this preprocessing, we treat low rank
updates to these original rows in the sparsifier and new rows separately. For the original rows we
can perform the preprocessing as in Vaidya [33] and the techniques in [33, 17] to obtain the desired
time bounds. For low rank updates to the new rows, we show how to use subspace embeddings
[26] to approximate their contribution to sufficient accuracy without hurting the asymptotic running
times. In Section 3 we show that this technique alone (even without use of sparsification), suffices to
mildly improve the running time of solving the inverse maintenance assumption with the `2 stability
assumption (but not to obtain the fastest running times in this paper).
The second difficulty is how to use the σ stability assumption to bound the number of updates to
our sparsifier; under this weak assumption where there can be many low rank changes to AT D(k) A.
Here we simply show that the σ stability assumption implies that the leverage scores do not change
too quickly in the norm induced by leverage scores (See Section 4). Consequently, if we sample rows
(k)
by leverage scores and only re-sample when either di or the leverage score changes sufficiently
then we can show that the expected number of low rank updates to our sparsifier is small. This
nearly yields the result, except that our algorithm will use its own output to compute the approximate leverage scores and the user of our algorithm might use the solvers to adversarially decide
the next d~(k) and this would break our analysis. In Section 5 we show how to make our solvers
indistinguishable from a true solver plus fixed noise with high probability, alleviating the issue.
Ultimately this yields amortized cost per iteration of Õ(nnz(A) + d2 ) for the inverse maintenance problem. We remark that our algorithm is one of few instances we can think of where an
algorithm needs to use both subspace embeddings [26] as well as iterative linear system solving and
sampling techniques [3, 21] for different purposes to achieve a goal; for many applications they are
interchangeable.
5

In Section 6 we also show that our approach can be generalized to accommodate for multiple
new rows to be added or removed from
the matrix. We also show that our inverse maintenance
√
problem yields an Õ((nnz(A) + d2 ) d log(−1 )) algorithm for solving a linear program improving
upon our previous work in [19] (See Section 7.1). In Section 7 we provide multiple applications of
our results to more specific problems including minimum cost flow, multicommodity flow, non-linear
regression, computing a rounding of a polytope, and convex optimization and conclude with an open
problem regarding the sampling of a random point in a polytope (Section 8).

2
2.1

Preliminaries
Notation

Throughout this paper, we use vector notation, e.g ~x = (x1 , . . . , xn ), to denote a vector and bold,
e.g. A, to denote a matrix. We use nnz(~x) or nnz(A) to denote the number of nonzero entries in a
vector or a matrix respectively. Frequently, for ~x ∈ Rd we let X ∈ Rd×d denote diag(~x),√the diagonal
def
matrix such that Xii = xi . Forq
a symmetric matrix, M and vector ~x we let k~xkM = ~xT M~x. For
P
def
2
vectors ~x and ~y we let k~y k~x =
i x i yi .
In most applications, we use d to denote the smaller dimension of the problem, which can be
either the number of variables or the number of constraints. For example, the linear programs
minA~x=~b,~l≤~x≤~u ~cT ~x we solve in Theorem 21 has n variables and d constraints (and use this formulation because it is more general than minA~x≥~b ~cT ~x).

2.2

Spectral Approximation

For symmetric matrices N, M ∈ Rn×n , we write N  M to denote that ~xT N~x ≤ ~xT M~x for all
~x ∈ Rn and we define N  M, N ≺ M and N  M analogously. We call a symmetric matrix M
positive definite (PD) if M  0. We use M ≈ N to denote the condition that e− N  M  e N.

2.3

Leverage Scores

Our algorithms make extensive use of leverage scores, a common measure of the importance of rows
of a matrix. We denote the leverage scores
Rn×d by ~σ ∈ Rn and say the leverage
−1 Tof a matrix A ∈
def
def
T
n×d
~ we use
score of row i ∈ [n] is σi = [A A A
A ]ii . For A ∈ R
, d~ ∈ Rn>0 , and D = diag(d)
1/2
~ to denote the leverage scores of the matrix D A. We frequently use well
the shorthand ~σA (d)
known facts regarding leverage scores, such as σi ∈ [0, 1] and ~σ 1 ≤ d. (See [32, 24, 21, 3] for a
more in-depth discussion of leverage scores, their properties, and their many applications.)
We use the following two key results regarding leverage scores. The first result states that one
can sample Õ(d) rows of A according to their leverage score and obtain a spectral approximation
of AT A [5]. In Lemma 5, we use a variant of a random sampling lemma stated in [3]. The second,
Lemma 6, is a generalization of a result in [32] that has been proved in various settings (see [17]
for example) that states that given a solver for a AT A one can efficiently compute approximate
leverage scores for A.
Lemma 5 (Leverage Score Sampling). Let A ∈ Rn×d and let ~u ∈ Rn be an overestimate of ~σ ,

def
the leverage scores of A; i.e. ui ≥ σi for all i ∈ [n]. Set pi = min 1, cs −2 ui log n for some
absolute constant cs > 0 and  ∈ (0, 0.5] and let H ∈ Rn×n be a random diagonal matrix where
independently Hii = p1i with probability pi and is 0 otherwise. With high probability in n, we have
nnz(H) = O( ~u 1 −2 log n) and AT HA ≈ AT A.
6

Lemma 6 (Computing Leverage Scores). Let A ∈ Rn×d , let ~σ denote the leverage scores of A, and
let  > 0. If we have a T -time solver for AT A then in time Õ((nnz(A) + T )−2 log(−1 )) we can
compute ~τ ∈ Rn such that with high probability in n, (1 − )σi ≤ τi ≤ (1 + )σi for all i ∈ [n].

2.4

Matrix Results

Our algorithms combine the sampling techniques mentioned above with the following two results.
Lemma 7 (Woodbury Matrix Identity). For any matrices A, U, C, V with compatible size, if A
−1
and C are invertible, then, (A + UCV)−1 = A−1 − A−1 U C−1 + VA−1 U
VA−1 .
Theorem 8 (Theorem 9 in [26] ). There is a distribution D over Õ(d−2 ) × n matrices such that for
any A ∈ Rn×d , with high probability in n over the choice of Π ∼ D, we have AT ΠT ΠA ≈ AT A.
Sampling from D and computing ΠA for any Π ∈ supp(D) can be done in Õ(nnz(A)) time.

3

Solving the Inverse Maintenance Problem Using `2 Stability

In this section we provide an efficient algorithm to solve the inverse maintenance problem under
the `2 stability assumption (See Section 1.1). The `2 stability assumption is stronger than the σ
stability assumption and the result we prove in this section is weaker than the one we prove under
the σ stability assumption in the next section (although still a mild improvement over many previous
results). This section serves as a “warm-up” to the more complicated results in Section 4.
The main goal of this section is to prove the following theorem regarding exactly maintaining
an inverse under a sequence of low-rank updates. We use this result for our strongest results on
solving the inverse maintenance problem under the σ stability assumption in Section 4.
Theorem 9 (Low Rank Inverse Maintenance). Suppose we are given a matrix A ∈ Rn×d , a vector
~b(0) ∈ Rd , a number of round r > 0, and in each round k ∈ [r] we receive d~(k) ∈ Rn such that
>0
>0
def
(k)
(k−1)
B(k) = AT D(k) A is PD. Further, suppose that the number of pairs (i, k) such that di 6= di
−1
(0)
(k)
(0)
is bounded by α ≤ d and suppose that there is β > 2 such that β B  B  βB . Then, in
−1
round k, we can implicit construct a symmetric matrix K such that K ≈O(1) B(k)
such that we
2
can apply K to an arbitrary vector in time Õ(d log(β)). Furthermore, the algorithm in total takes
ω−2
def
time Õ(sdω−1 + αrd αr
+ rαω ) where s = max{nnz(d~0 ), d}.
This improves upon the previous best expected running times in [33, 17] which had an additive
Õ(ndω−1 ) term which would be prohibitively expensive for our purposes.
To motivate this theorem and our approach, in Section 3.1, we prove that Theorem 9 suffices
to yield improved algorithms for the inverse maintenance problem under `2 stability. Then in
Section 3.2 we prove Theorem 9 using a combination of classic machinery involving low rank update
formulas and new machinery involving subspace embeddings [26].

3.1

Inverse Maintenance under `2 Stability

Here we show how Theorem 9 can be used to solve the inverse maintenance problem under the `2
stability assumption. Note this algorithm is primarily intended to illustrate Theorem 9 and is a
warm-up to the stronger result in Section 4.
Theorem 10. Suppose that the inverse maintenance problem satisfies the `2 stability assumption. Then Algorithm 2 maintains a Õ(nnz(A) + d2 )-time solver with high probability in total time
 2ω
Õ(sdω−1 + r ndω−1 2ω+1 ) where s = max{maxk∈[r] nnz(d(k) ), d} and r is the number of rounds.
7

Algorithm 2: Algorithm for the `2 Stability Assumption
Input: Initial point d~(0) ∈ Rn>0 .
Set d~(apr) := d~(0) .
def

Q(0) = AT D(apr) A.
Let K(0) be an approximate inverse of Q(0) computed using Theorem 9.
Output: A Õ(d2 + nnz(A))-time linear solver for AT D(0) A (using Lemma 31 on K(0) ).
for each round k ∈ [r] do
Input: Current point d~(k) ∈ Rn>0 .
for each coordinate i ∈ [n] do
(apr)
(k)
(apr)
if 0.9di
≤ di ≤ 1.1di
is false then
(apr)

di
end

(k)

:= di .

def

Q(k) = AT D(apr) A.
Let K(k) be an approximate inverse of Q(k) computed using Theorem 9.
Output: A Õ(d2 + nnz(A))-time linear solver for AT D(k) A (using Lemma 31 on K(k) ).
end

Proof. Recall that by the `2 stability assumption log(d~(k) )−log(d~(k−1) ) 2 ≤ 0.1 for all k. Therefore,
in the r rounds of the inverse maintenance problem at most O(r2 ) coordinates of d~(k) change by
(apr)
(k)
(apr)
any fixed multiplicative constant. Consequently, the condition, 0.9di
≤ di ≤ 1.1di
in
2
Algorithm 2 is false at most O(r ) times during the course of the algorithm and at most O(r2 )
coordinates of the vector d~(apr) change during the course of the algorithm.
Therefore, we can use Theorem 9 on AT D(apr) A with α = O(r2 ) and s as defined in the theorem
−1
statement. Consequently, the total cost of maintaining an implicit approximation of AT D(apr) A

for r rounds is Õ sdω−1 + drω+1 + r2ω+1 .
1

To further improve the running time we restart the maintenance procedure every (sdω−1 ) 2ω+1
iterations. This yields a total cost of maintenance that is bounded by
&
'
!

ω+1 
2ω+1 
1
1
r
Õ
sdω−1 + d (sdω−1 ) 2ω+1
+ (sdω−1 ) 2ω+1
.
1
ω−1
2ω+1
(sd
)
which is the same as

&
Õ

'

r



1

ω−1

sd

ω−1

+ d(sd

)

ω+1
2ω+1



!
.

(sdω−1 ) 2ω+1
√
ω+1
ω+1
Since ω ≤ 1 + 2, we have d(dω ) 2ω+1 ≤ dω and thus sdω−1 dominates d(sdω−1 ) 2ω+1 when s = d.
ω+1
Furthermore, as s ≥ d the term sdω−1 grows faster than d(sdω−1 ) 2ω+1 . Consequently, sdω−1 always
dominates and we have the desired result.

3.2

Low Rank Inverse Maintenance

Here we prove Theorem 9 and provide an efficient algorithm for maintaining the inverse of a matrix
under a bounded number of low rank modifications. The algorithm we present is heavily motivated
by the work in [33] and the slight simplifications in [17]. However, our algorithm improves upon

8

the previous best cost of Õ(ndω−1 + drω+1 + r2ω+1 ) achieved in [17] by a novel use of subspace
embeddings [26] that we hope will be of independent interest.3
We break our proof of Theorem 9 into several parts. First, we provide a simple technical lemma
about maintaining the weighted product of two matrices using fast matrix multiplication.
Lemma 11. Let A, B ∈ Rn×d and suppose that in each round k of r we receive ~x(k) , ~y (k) ∈ Rn .
(k)
(k−1)
(k)
(k−1)
Suppose that the number of pairs (i, k) such that xi 6= xi
or yi 6= yi
is upper bounded by
α. Suppose that nnz(~x(1) ) ≤ α, nnz(~y (1) ) ≤ α and α ≤ d. With O(dαω−1 ) time precomputation, we
ω−2
can compute X(k) ABT Y(k) explicitly in an average cost of O(αd αr
) per iteration.
Proof. For the initial round, we compute X(0) ABT Y(0) explicitly by multiplying an α × d and a
d × α matrix. Since α ≤ d, we can compute X(0) ABT Y(0) by multiplying O αd matrices of size
α × α. Using fast matrix multiplication this takes time O(dαω−1 ).
(k)
(k) def
For all k ∈ [r] let ∆X = X(k) − X(k−1) and ∆Y = Y(k) − Y(k−1) . Consequently,




(k)
(k)
X(k) ABT Y(k) = X(k−1) + ∆X ABT Y(k−1) + ∆Y
.
Note that nnz(~x(k) ) and nnz(~y (k) ) are less than 2α. Thus, if we let uk denote the number of
(k)
(k−1)
(k)
(k−1)
coordinates i such that xi 6= xi
or yi 6= yi
, then to compute X(k) ABT Y(k) we need only
multiply a uk × d with a d × uk matrix and multiply two 2α × d matrices with d × uk matrices. Since
uk ≤ α, the running time is dominated by the
 compute the 2α × d and a d × uk . Since
 time to
d
α
uk ≤ α ≤ d, we can do this by multiplying O uk · uk matrices of size uk × uk . Using fast matrix
).
multiplication, this takes time O(αduω−2
k
Summing over all uk we see that the total cost of computing the X(k) ABT Y(k) is

!
!ω−2 

r
r
 α ω−2 
X
X
1
ω−2
 ≤ O αdr
O
αduk
≤ O αdr
uk
,
r
r
k=1

k=1

where in the second inequality we used the concavity of xω−2 . Since this is at least the O(dαω−1 )
cost of computing X(0) ABT Y(0) we obtain the desired result.
Next, for completeness we prove a slightly more explicit variant of a Lemma in [17].
def
Lemma 12. Let A ∈ Rn×d , d~(0) , . . . , d~(r) ∈ Rn>0 , and B(k) = AT D(k) A for all k. Suppose that the
ω−2
(k)
(k−1)
number of pairs (i, k) such that di 6= di
is bounded by α and α ≤ d. In O(ndω−1 + αdr αr
)
total time, we can explicitly output C ∈ Rn×d and a V(k) ∈ Rn×n in each round k such that



B(k)

−1


−1
= B(0)
− CT ∆(k) V(k) ∆(k) C
(k)

(k)

(0)

where ∆(k) is a n × n diagonal matrix with ∆ii = d~i − d~i . Furthermore, V(k) is an Õ(α) × Õ(α)
matrix if we discard the zero rows and columns.
def

Proof. Let B = B(0) for notational simplicity. Note that we can compute B directly in O(ndω−1 )
time using fast matrix multiplication. Then, using fast matrix multiplication we can then compute
def
B−1 explicitly in O(dω ) time. Furthermore, we can compute C = AB−1 in O(ndω−1 ) time.
3

We require a slightly stronger assumption than [17] to achieve our result; [17] did not require any bound on β.

9

Let uk be the number of non-zero entries in ∆(k) and P(k) be a uk × n matrix such that
(k)
the Pij = 1 if j is the index of the ith non-zero diagonal entries in ∆(k) and 0 otherwise. Let

† (k) T
S(k) = P(k) ∆(k)
P
. By definition of S(k) and P(k) , we have

T
∆(k) = ∆(k) P(k) S(k) P(k) ∆(k) .
By the Woodbury matrix identity (Lemma 7), we know that
−1

T

−1 
(k) (k) (k)
T (k)
(k)
(k)
S P ∆ A
= B+A ∆
P
B

T
= B−1 − CT ∆(k) P(k) T(k) P(k) ∆(k) C
where T(k) = ((S(k) )−1 + P(k) ∆(k) CAT ∆(k) (P(k) )T )−1 . Now by Lemma 11, we know that we
ω−2
) time, using that C = AB−1 was
can maintain ∆(k) CAT ∆(k) in an additional O(αdr αr
precomputed in O(ndω−1 ) time. Having the matrix ∆(k) CAT ∆(k) explicitly, one can compute
T
P(k) ∆(k) CAT ∆(k) P(k) in O(u2k ) time and hence compute T(k) in O(uωk ) time using fast matrix

P ω
ω−1 + αdr α ω−2 +
multiplication.
Hence,
the
total
running
time
is
O(nd
uk ). The convexity of
r
P
xω yields
uωk ≤ αω . Hence, the total time is bounded by
O(ndω−1 + αdr

 α ω−2
r

+ αω ) = O(ndω−1 + αdr

 α ω−2
r

).

T
By setting V(k) = P(k) T(k) P(k) , we have the desired formula. Note that T(k) is essentially
V(k) with the zero columns and rows and hence we have the desired running time.
We now everything we need to prove Theorem 9.
Proof of Theorem 9. Let S ⊂ [n] denote the indices for which d~(0) is nonzero. Furthermore, let
us split each vector d~(k) into a vector ~e(k) for the coordinates in S and f~(k) for the coordinates
not in S, i.e. ~e(k) , f~(k) ∈ Rn>0 such that d~(k) = ~e(k) + f~(k) . By the stability guarantee, we know
1 ~(0)
β −1 B(0)  B(k)  βB(0) for some β. We make AT E(k) A invertible for all k by adding 10β
d to
all of d~(k) and ~e(k) ; we will show this only increases the error slightly.
−1
Now, we can compute B(0) and B(0)
in O(sdω−1 ) time using fast matrix multiplication
where s = |S|. Furthermore, using Lemma 12 we can compute C and maintain V(k) such that

−1 
−1
AT E(k) A
= B(1)
− CT ∆(k) V(k) ∆(k) C
in total time O(sdω−1 + αdr


α ω−2
). All that remains
r
−1
of AT E(k) A
and the

is to maintain the contribution from f~(k) .

Using our representation
Woodbury matrix identity (Lemma 7), we
have

−1 
−1
B(k)
= AT E(k) A + AT F(k) A

−1

†
T (k)
T (k)
(k)
(k)
= A E A+A F
F
F A

−1 
−1
−1


−1
= AT E(k) A
− AT E(k) A
AT F(k) M(k)
F(k) A AT E(k) A

10

(3.1)

where
(k)

M

=



F

(k)

† −1


−1
AT F(k)
+ F(k) A AT E(k) A

(3.2)

Now note that

−1 

−1

−1
AT F(k)
AT F(k) = F(k) A AT E(k) A
AT E(k) A AT E(k) A
F(k) A AT E(k) A

T
= N(k) N(k)
where
−1

1/2 
def
AT F(k)
A AT E(k) A
N(k) = E(k)


1/2 
1/2 
1/2  
−1
(0)
(k)
(0)
(0)
T (k) (k) (k)
=
D
+
E
− D
A
B
− C ∆ V ∆ C AT F(k) .
Note that computing the D(0)

1/2

−1

AT F(k) term directly would be prohibitively exT
pensive. Instead, we compute a spectral approximation to N(k) N(k) and show that suffices.
T
Since N(k) N(k) is a rank α matrix, we use Theorem 8 to Π ∈ RÕ(α)×d such that


N(k)

A B(0)

T

T

ΠT ΠN(k) ≈1 N(k) N(k)

(3.3)

for all k with high probability. Now, we instead consider the cost of maintaining

T
N(k) ΠT ΠN(k) .
To see the cost of maintaining ΠN(k) , we separate the matrix as follows:
1/2 
−1

ΠN(k) = Π D(0)
A B(0)
AT F(k)

−1
1/2  
1/2 
(0)
(k)
AT F(k)
A B(0)
− D
+Π
E

(3.4)


1/2
−Π D(0)
ACT ∆(k) V(k) ∆(k) CAT F(k)

1/2 
1/2 
(k)
(0)
−Π
E
− D
ACT ∆(k) V(k) ∆(k) CAT F(k) .
1/2
For the first term in (3.4), note that D(0)
A is a s × d matrix and hence we can precompute
1/2
−1
1/2
−1
(0)
(0)
ω−1
D
A B
in Õ(sd
) time and therefore precompute Π D(0)
A B(0)
in the
same amount of time. Note that Π is a Õ(α) × d matrix, so, we can write

1/2 
−1
Π D(0)
A B(0)
AT F(k) = ΛKAT F(k)
where K is some explicitly computed n × d matrix and Λ is a diagonal matrix with only Õ(α)
non-zeros. Consequently, we can use Lemma 11 to maintain Π(D(0) )1/2 A(B(0) )−1 AT F(k) in total
ω−2
time Õ(sdω−1 + αdr αr
).
11

−1
For the second term in (3.4), we can precompute A B(0)
in Õ(sdω−1 ) time. Therefore, using



1/2
1/2
−1 T (k)
ω−2 
Lemma 11 we can maintain E(k)
− D(0)
A B(0)
A F in total time Õ αdr αr

1/2
1/2 
−1 T (k)
and by Theorem 8 we can maintain Π E(k)
− D(0)
A B(0)
A F
in the same
amount of time.
1/2
For the last two terms in (3.4), we use Lemma 11 on ∆(k) CAT F(k) , Π D(0)
ACT ∆(k) and




ω−2
1/2
1/2
. Note that all of those matrices
− D(0)
ACT ∆(k) all in total time Õ αdr αr
E(k)
including V(k) are essentially Õ(α) × Õ(α) matrices if we discard the zero rows and columns. So,
having those matrices explicitly computed, we can compute the last two terms in an additional of
Õ(αω ) time per iteration.
T
Finally computing N(k) ΠT ΠN(k) only requires an additional Õ(αω ) time per-iteration.
T
Hence, the total cost of maintaining N(k) ΠT ΠN(k) is


 α ω−2
ω−1
ω
+ sd
+ rα .
Õ αdr
r
Using (3.1) and (3.2), we have shown how to approximate B(k)
compute a better approximation. Recall from (3.1) that


B(k)

−1

−1

. Now, we show how to


−1 
−1

−1

−1
= AT E(k) A
+ AT E(k) A
AT F(k) M(k)
F(k) A AT E(k) A
.

−1
Using Lemma 12 as we have argued, we can apply the first term AT E(k) A
exactly to a vector
2
in Õ(d ) time. The only difficulty in applying the second term to a vector comes from the term
−1
M(k)
. Using (3.3), we have
(k)

M



(k)

≈1 F

†



(k)

+ N

T

ΠT ΠN(k)

and hence


Since computing



F(k)

(k)

M

†

−1

+ N(k)

≈1
T



F

(k)

†

ΠT ΠN(k)



(k)

+ N

−1

T

T

(k)

Π ΠN

−1
.

only requires Õ(αω ) time and (3.2) shows that

we can apply M(k) to a vector exactly in Õ(d2 ) time, by Lemma 31 we have a symmetric matrix


(k)
(k) −1
M ≈ M(k) such that we can apply M
to a vector in Õ(d2 log(−1 )) for any  > 0. Hence,
we obtain

−1

−1

−1
T (k)
AT F(k) M(k)
Λ(k)
F(k) A AT E(k) A
 ≈ A E A
(k)

such that we can apply Λ to a vector in Õ(d2 log(−1 )) time. All that remains is to compute
−1
what value of  is needed for this to be a spectral approximation to the B(k)
.
−1
(0)
(k)
(0)
Recall that by the assumption β B  B  βB . As we mentioned in the beginning, we
1 ~(0)
replaced d~(k) with d~(k) + 10β
d and compute a constant spectral approximation to the new B(k)
that will suffice for the theorem statement. Consequently

−1 
−1

−1
def
Λ(k) = AT E(k) A
− B(k)
 10β B(0)

12

Furthermore, since f~(k) ≥ 0 we have 0  Λ(k) and therefore

−1

−1
Λ(k) − 10β B(0)
 Λ(k)  Λ(k) + 10β B(0)
.
Using the assumption β −1 B(0)  B(k)  βB(0) again, we have

−1

−1
 Λ(k)  Λ(k) + 10β 2 B(k)
.
Λ(k) − 10β 2 B(k)
1
T (k) A
Picking  = O( 20β
2 ), we have a matrix A E

−1

β −1 B(0)

Õ(d2 log(β))

(k)

≈O(1) B(k)

B(k)

βB(0)

− Λ

−1

which we can apply in

time. Furthermore, again since


we see that our replacement
1 ~(0)
(k)
(k)
~
~
of d d with d + 10β d only affected our approximation quality by a constant.

An algorithm for the σ Stable Case

4

In this section, we provide our algorithm for solving the inverse maintenance problem under the σ
stability assumption. The central result of this section is the following:
Theorem 13. Suppose that the inverse maintenance problem satisfies the σ stability assumption.
Then Algorithm 3 maintains a Õ(nnz(A)+d2 )-time solver with high probability in total time Õ(dω +
r(nnz(A) + d2 )) where r is the number of rounds.
To prove this we first provide a technical lemma showing that leverage scores are stable in
leverage score number assuming σ stability (See Section 4.1). Using this lemma we then prove the
Theorem 13 (See Section 4.2). This proof will assume that the randomness we use to maintain
our solvers is independent from the output of our solvers. In Section 5 we show how to make this
assumption hold.

4.1

Leverage Scores are Stable under σ Stability

Here we show that leverage scores are stable in the leverage score norm assuming σ stability. This
technical lemma, Lemma 14, is crucial to showing that we do not need to perform too many low-rank
updates on our sparsifier in our solution to the inverse maintenance problem under the σ stability
assumption. (See Section 1.3 for more intuition.)
Lemma 14. For all A ∈ Rn×d and any vectors ~x, ~y ∈ Rn>0 such that ln(~x) − ln(~y )
have
ln ~σA (~x) − ln ~σA (~y ) ~σ (~x) ≤ e ln ~x − ln ~y ~σ (~x) .
A

∞

≤ , we

A

Proof. For t ∈ [0, 1], let ln θ~t denote a straight line from ln ~x to ln ~y with θ~0 = ~x and θ~1 = ~y or
def
equivalently let θ~t = exp(ln x + t(ln y − ln x))). Since ln(~x) − ln(~y ) ∞ ≤  we have for all i ∈ [n]
√
−1 T √
~σA (~x)i = ~1Ti XA AT XA
A X~1i
√

−1 T √
≈2 ~1Ti ΘA AT ΘA
A Θ~1i
= ~σA (θ~i )
(4.1)
Consequently, for all ~z, we have ~z
ˆ
ln ~σA (~x) − ln ~σA (~y )

σA (~
x)

≤
0

~
σA (θ~t )
1

≤ e ~z

~
σA

and by Jensen’s inequality we have


d
~
ln ~σA (θt ) dt
dt
13

ˆ


≤e ·
σA (~
x)

0

1

d
ln ~σA (θ~t )
dt

dt (4.2)
σA (θ~t )

Now, for all z ∈ Rn>0 let Jln ~z (ln ~σA (~z)) denote the Jacobian matrix

h

i

∂ ln ~
σA (~
z )i
∂ ln ~
zj
ij

=

h

i

∂ ln ~
σA (~
z )i
zj
∂~
zj
ij

and suppose that for all ~z and ~u we have
Jln ~z (ln ~σA (~z))~u
Then
ˆ 1
0

d
ln ~σA (θ~t )
dt

ˆ

1

dt =
σA (θ~t )

0

~
σA

Jln θ~t (ln ~σA (θ~t ))



≤ ~u

~
σA (~
z)

d ~
ln θt
dt

.


~
σA (θ~t )

(4.3)

≤ e ln ~x − ln ~y

~
σA (~
x)

(4.4)

where again used 4.1 as well as the definition of θ~t . All that remains is to prove (4.3). For this,
we first note that in [17, Ver 1, Lemma 36] we showed that J~z (~σA (~z)) = (ΣA (~z) − M) Z−1 where
√
√ 2
def
Mij (~z) =
ZA(AT ZA)−1 AT Z . Consequently Jln ~z (ln ~σA (~z)) = ΣA (~z)−1 (ΣA (~z) − M(~z)).
Now note that
X
X √
√ 2
Mij =
ZA(AT ZA)−1 AT Z
i

ij

i

E
D√
√
√
√
ZA(AT ZA)−1 AT Z~1j , ZA(AT ZA)−1 AT Z~1j
=
√
√ 
ZA(AT ZA)−1 AT Z
=
jj

= (~σA (~z))j .
Consequently, ΣA (~z) − M is a symmetric diagonally dominant matrix and therefore ΣA (~z) 
ΣA (~z) − M  0 and 0  ΣA (~z)1/2 Jln ~z (ln ~σA (~z))ΣA (~z)−1/2  I. Using this fact, we have that for
all ~z ∈ Rn>0 and ~u ∈ Rn
Jln ~z (ln ~σA (~z))~u

~
σA (~
z)

= ΣA (~z)1/2 ~u

(ΣA (~z)−1/2 (ΣA (~z)−M)ΣA (~z)−1/2 )

2

≤ ~u

~
σA (~
z)

and this proves (4.3). Combining (4.2) and (4.4) yields the result.

4.2

Algorithm for σ Stability

Here we prove Theorem 13. The full pseudocode for our algorithm, with the exception of how we
compute leverage scores and maintain the inverse of AH(k) A can be seen in Algorithm 3.
First, we bound the number of coordinates H that will change during the algorithm.
Lemma 15. Suppose that the changes of d~ and the error occurred in computing ~σ is independent
of our sampled matrix. Under the σ stability guarantee, during first r iterations of Algorithm 3, the
expected number of coordinates changes in H(k) over all k is O(r2 log d).
Proof. Since our error in computing σ is smaller than the re-sampling threshold on how much change
~ i or di changes by more
of σ, the re-sampling process for the ith row happens only when either σA (d)
T
than a multiplicative constant. In order for re-sampling to affect A HA, it must be the case that
it is either currently in AT HA or about to be put in AT HA. However, since whenever re-sampling
occurs the resampling probability has changed by at most a multiplicative constant, we have that
~ i ) using the independence between d~ and the
both these events happen with probability O(γ · σA (d)
approximate ~σ . By union bound we have that the probability of sampling row i changing the matrix
~ i log(d)).
AT HA is O(σA (d)
14

Algorithm 3: Algorithm for σ Stability Assumption
Input: Initial point d~(0) ∈ Rn>0 .
def
Set d~(old) := d~(0) and γ = 1000cs log d where cs defined in Lemma 5.
(apr)
(apr)
Use Lemma 6 to find σ (apr) such that 0.99σ
≤ σ(d~(0) )i ≤ 1.01σ
.
i
(apr)

(0)

i

(apr)

For each i ∈ [n] : let hi := di / min{1, γ · σi
} with probability min{1, γ · σi
}
and is set to 0 otherwise.
def
Q(0) = AT H(0) A.
Let K(0) be an approximate inverse of Q(0) computed using Theorem 9.
Output: A Õ(d2 + nnz(A))-time linear solver for AT D(0) A (using Theorem 10 on K(0) ).
for each round k ∈ [r] do
Input: Current point d~(k) ∈ Rn>0 .
Use Lemma 6 and the solver from the previous round to find σ (apr) such that
(apr)
(apr)
0.99σi
≤ σ(d~(k) )i ≤ 1.01σi
.
for each coordinate i ∈ [n] do
(old)
(apr)
(old)
(old)
(k)
(old)
if either 0.9σi
≤ σi
≤ 1.1σi
or 0.9di
≤ di ≤ 1.1di
is violated then
(old)
(k)
di
:= di .
(old)
(apr)
σi
:= σi
.
(k)
(k)
(apr)
(apr)
hi := di / min{1, γ · σi
} with probability min{1, γ · σi
}
and is set to 0 otherwise.
else
(k)
(k−1)
hi := hi
.
end
end
def

Q(k) = AT H(k) A.
Let K(k) be an approximate inverse of Q(k) computed using Theorem 9.
Output: A linear Õ(d2 + nnz(A))-time solver for AT D(k) A (using Theorem 10 on K(k) ).
end
~ i or di has changed by more
Observe that whenever we re-sampled the ith row, either σA (d)
than a multiplicative constant. Let k1 be the last iteration we re-sampled the ith row and let
P 2 −1
k2 be the current iteration. Then, we know that kk=k
ln σA (d~(k+1) )i − ln σA (d~(k) )i = Ω(1) or
1
Pk2 −1
~(k+1) − ln d~(k) ≥ Ω(1). Since there are only r steps, we have |k2 − k1 | ≤ r and hence
i
k=k1 ln di
kX
2 −1 

ln σA (d~(k+1) )i − ln σA (d~(k) )i

2

 
 
kX
2 −1 

1
1
(k+1)
(k) 2
~
~
or
ln di
− ln di
=Ω
.
=Ω
r
r
k=k1

k=k1

~ i does not change more by a constant between re-sample (by σ-stability assumption),
Since σA (d)
we have
!
kX
2 −1

2
(k2 ) )
~
σ
(
d
i
A
or
σA (d~(k) )i ln σA (d~(k+1) )i − ln σA (d~(k) )i
= Ω
r
k=k1
!
kX
2 −1

2
(k2 ) )
~
σ
(
d
(k+1)
(k)
i
A
σA (d~(k) )i ln d~i
− ln d~i
= Ω
.
r
k=k1

15

In summary, re-sampling the ith row indicates that the sum of the σ norm square of the changes
~ i /r and with O(σA (d)
~ i log d) probability, the sampled matrix
of either ln σA or ln d is at least σA (d)
P
P
2
(k+1)
(k)
is changed. Since
ln d~
− ln d~
ln σA (d~(k+1) ) −
~(k) ≤ O(r) and by Lemma 14
k

ln σA (d~(k) )

2
σA (d~(k) )

σA (d

)

k

≤ O(r) we have the desired result.

We now everything we need to prove our main theorem assuming that the changes of d~ and the
error occurred in computing σ is independent of our sampled matrix. (In Section 5 we show how to
drop this assumption.)
Proof of Theorem 13. Note that by design in each iteration k ∈ [r] we have that D(old) ≈0.2 D(k) ,
def
Σ(old) ≈0.2 Σ(d~(k) ) where Σ = diag(σ). Thus, we see that the sample probability was chosen
precisely so that we can apply Lemma 5. Hence, we have AT H(k) A ≈0.1 AT D(k) A. Thus, the
algorithm is correct, it simply remains to bound the running time.
To maintain inverse of AT H(k) A, we can simply apply Theorem 9. By Lemma 5, we know that
with high probability nnz(H(0) ) ≤ Õ(d). By Lemma 15 we know there are only Õ(r2 ) coordinate
changes during the algorithm in expectation. Therefore, we can use Theorem 9 on AT H(k) A with
α = Õ(r2 ) and s = Õ(d). Hence, the average cost of maintain a linear Õ(nnz(A) + d2 )-time solver
−1

ω
of AT H(k) A
is Õ dr + drω + r2ω . Similar to Theorem 10 , we restart the algorithm every
√
1
2ω 2
ω + rd 2ω+1 ). Since ω ≤ 1 + 2,
r = (ndω−1 ) 2ω+1 and see that the total cost of maintenance
is
Õ(d

we have the total maintenance cost is Õ dω + rd2 .
Thus, all the remains is to bound the cost of computing ~σ (k) . However, since by the σstability
assumption log(d~k+1 ) − log(d~k ) ∞ ≤ 0.1 we know that D(k) ≈0.1 D(k+1) and thus Q(k) ≈0.2
AT D(k+1) A. Thus, using Lemma 6 we can compute ~σ (k) for all k ≥ 2 in Õ(nnz(A) + d2 log β) time
within 0.01 multiplicative factor. Furthermore, using this same Lemma and fast matrix multiplication, we can compute ~σ (1) in Õ(nnz(A) + dω ) time; therefore, we have our result.
Note that Lemma 15 assume that the changes of d~ and the error occurred in computing σ is
independent of our sampled matrix. Given any linear solver, Theorem 17 in Section 5 shows how to
construct a solver that has same running time up to Õ(1) factor and is statistically indistinguishable
−1
and thus circumvent this issue; thereby completing the proof.
from the true solver AT D(k) A

5

Hiding Randomness in Linear System Solvers

In many applications (see Section 7), the input to a particular round of the inverse maintenance
problem depends on our output in the previous round. However, our solution to the inverse maintenance problem is randomized and if the input to the inverse maintenance problem was chosen
adversarially based on this randomness, this could break the analysis of our algorithm. Moreover,
even within our solution to the inverse maintenance problem we needed to solve linear systems and
if the output of these linear systems was adversarially correlated with our randomized computation
our analysis would break (see Section 4)
In this section, we show how to fix both these problems and hide the randomness we use to approximately solve linear system. We provide a general transformation, NoisySolver in Algorithm 4,
which turns a linear solver for AT A into a nonlinear solver for AT A that with high probability
is indistinguishable from an exact solver for AT A plus a suitable Gaussian noise. The algorithm
simply solves the desired linear system using the input solver and then add a suitable Gaussian
noise.

16

Algorithm 4: NoisySolver(~b, )
Input: a linear T -time solver S of AT A, vector ~b ∈ Rn , and accuracy  ∈ (0, 1/2).
−2
~y1 := S(~b, 1 ) where 1 =  32n7
.
n
~
Let ~η ∈ R be sampled randomly from
−2 the normal distribution with mean 0 and covariance I.
T
6
~y2 := S(A ~η , 2 ) where 2 = 12dn
.
p
def
1
Output: ~y = ~y1 + α~y2 where α = 8 n k~y1 kAT A .
We break the proof that this works into two parts. First, in Lemma 16 we show that NoisySolver,
is in fact a solver and then in Theorem 17 we show that with high probability it is indistinguishable
from an exact solver plus a Gaussian noise.
Lemma 16. Let A ∈ Rn×d and let S be a linear T -time solver for AT A. For all ~b ∈ Rn and
 ∈ (0, 1/2), the algorithm NoisySolver(~b, ) is a (nnz(A) + T )-time solver for AT A.
Proof. By the definition of ~y and the inequality (a + b)2 ≤ 2a2 + 2b2 , we have
−1
~b
~y − AT A

2
AT A

≤ 2 ~y1 − AT A

−1

~b

2
AT A

To bound the first term, recall that by the definition of S
−1 2

~b T ≤ 1 AT A −1 ~b
~y1 − AT A
A A

+ 2α2 ~y2

2
.
AT A

2
.
AT A
2

To bound the second term, we note that by the definition of ~η , ~η 2 follows χ2 -distribution with n
degrees of freedom. It is known that [16, Lem 1] for all t > 0,


√
2
P ~η 2 > n + 2 nt + 2t ≤ exp(−t) .
Hence, with high probability in n, ~η
~y2

2
2

< 2n. Using this and the definition of S yields
−1 T 2
−1 T 2
≤ 2 ~y2 − AT A
A ~η AT A + 2 AT A
A ~η AT A

−1 T
2
≤ 2 (1 + 2 ) AT A
A ~η AT A

2
AT A

≤ 4 ~η

2
2

≤ 8n.

−1 T
where we used that 2 ≤ 1 and A AT A
A 2 ≤ 1. Using that 1 ≤ 1 and applying a similar
proof yields that
−1 2


2
~b
~y1 AT A ≤
AT A
α2 =
64n
32n
AT A
Consequently, with high probability in n,
~y − AT A

−1

~b

2
AT A

≤ 21

−1
~b
AT A

2
AT A

+ 16α2 n ≤ 

−1
~b
AT A

2
AT A

.

−1
~b + ~c where ~c follows the normal distribution with
Theorem 17. Let IdealSolver(~b) = AT A

−1
p
−1
~b
mean 0 and covariance β 2 AT A
where β = 81 n AT A
. Then, the total variaT
A A

tion distance between the outcome of NoisySolver and IdealSolver is less than 1/n3 . Therefore,
any algorithm calls IdealSolver less than n2 times cannot distinguish between NoisySolver and
IdealSolver.
17

Proof. Since S is linear, we have ~y2 = Q2 AT ~η for some matrix Q2 . Since ~η follows the normal
distribution with mean 0 and covariance matrix I, we have




E ~y2 ~y2T = Q2 AT E ~η ~η T AQT2 = Q2 AT AQT2
Lemma 32 shows that QT2 AT AQ2 ≈4√2 (AT A)−1 and hence


−1
E ~y2 ~y2T ≈4√2 AT A
.

(5.1)

Since ~y2 is a linear transformation of ~η , ~y2 follows is given by the normal distribution with mean 0
and covariance Q2 AT AQT2 , i.e. ~y2 ∈ N (0, Q2 AT AQT2 ).
Now let, ~y , ~z be the output of NoisySolver(~b) and IdealSolver(~b) respectively; we know that

~y ∈ N ~y1 , α2 Q2 AT AQT2 and

−1
−1 
T
2
T
~
~z ∈ N A A
b, β A A

p
1
T A −1 ~
A
b T . To bound ~y − ~z tv , Pinsker’s inequalk~
y
k
and
β
=
T
1
A
A
n
8
n
A A
q
1
ity shows that ~y − ~z tv ≤ 2 DKL (~y ||~z) and using an explicit formula for the KL divergence for
where α =

1
8

p

the normal distribution we in turn have that ~y − ~z tv is bounded by
v
!
u
 2
1u
α
2
ttr
AT AQ2 AT AQT2 + ~y1 − (AT A)−1 ~b β −2 AT A − d + ln
2
β
Hence, we need to prove α2 Q2 AT AQT2 ≈ β 2 AT A
−1
~b:
First, we first prove ~y1 ≈ AT A
−1
~b
~y1 − AT A

2
β −2 (AT A)

≤ β −2 1

−1

AT A

and ~y1 ≈ AT A

−1

~b

det β 2 (AT A)−1
det α2 Q2 AT AQT2

−1

~b

2

(Definition of Sand ~y1 )

AT A

64n1
.
(Definition of β)

−1
AT A , we note that by triangle inequality and the definition
≤

Next, to prove α2 Q2 AT AQT2 ≈ β 2
of S and ~y1 we have
−1
√
~b
(1 − 1 ) AT A

!

AT A

≤ k~y1 kAT A ≤ (1 +

√

1 )

AT A

−1

~b

AT A

.

Therefore by the definition of α and β we have
√
√
α2
(1 − 3 1 ) ≤ 2 ≤ (1 + 3 1 ) .
β
Using (5.1) then yields that
−1
α2 Q2 AT AQT2 ≈4√2 +4√1 β 2 AT A
.
Therefore, we have




1/2
1/2
α 2 T
α 2
T
T
T
T
T
T
tr ( ) A AQ2 A AQ2
= tr ( ) A A
Q2 A AQ2 A A
β
β
≤ de4

√

√
2 +4 1

√
√
≤ d + 8 2 d + 8 1 d
18

and
ln

−1 !


1/2 −1
1/2
det β 2 AT A
α 2
T
T
T
T
= ln det ( ) A A
Q2 A AQ2 A A
det α2 Q2 AT AQT2
β
 √
√ 
≤ d ln e4 2 +4 1
√
√
= 4 2 d + 4 1 d.

Combining inequalities above and using our choice of 1 and 2 we have
r
√
√
1
1
~y − ~z tv ≤
12 2 d + 12 1 d + 64n ≤ 1/n3 .
2


6

Row Insertion and Removal

For some applications we need to solve the inverse maintenance problem when the matrix A might
change between rounds. In particular some rows of A might be entirely removed or some new rows
might be added.
For example, in each iteration of cutting plane methods, a new constraint is added to the current
polytope {A~x ≥ ~b}; each iteration of quasi newton methods, the current approximate Hessian is
updated by a rank 1 matrix. If the matrix is updated only by a rank 1 matrix each iteration, the
inverse can be updated efficiently using the Sherman-Morrison formula. However, for the general
case, it is less obvious when the matrix A and the diagonal D can be changed by a high rank matrix.
Here we formally define the more general set of assumptions under which we would like to solve
the inverse maintenance problem and show how to solve this problem efficiently.
Definition 18 (K Stability Assumption). The inverse maintenance problem satisfies the K stability
assumption if
1. A row of A is revealed to the algorithm at round k if and only if the corresponding entry in
d~(k) is non-zero.
2. For each k ∈ [r], we have either
(a)

log(d~(k) ) − log(d~(k−1) )

~
σA (d~(k) )

≤ 0.1 and

log(d~(k) ) − log(d~(k−1) )

∞

≤ 0.1; or

(b) The vectors d~(k) and d~(k−1) differ in at most K coordinates.
3. β −1 AT D(0) A  AT D(k) A  βAT D(0) A for β = poly(n).
Some of the previous algorithms for the inverse maintenance, such as [33], rely on the assumption
that the entire matrix A is given explicitly. In particular these algorithms perform precomputation
on the entirety of A that they use to decrease the amortize cost of later steps. Here we show that
the Algorithm 3 we proposed does not have this drawback and can be easily modified to solve this
version of inverse maintenance problem under the K stability assumption for a fairly large K.
Theorem 19. Suppose that the inverse maintenance problem satisfies the K stability assumption
for K ≤ d(3−ω)/(ω−1) , where ω is the matrix multiplication constant. Then there is a variant
of Algorithm 3 that maintains a Õ(nnz(A) + d2 )-time solver with high probability in total time
Õ(dω + r(nnz(A) + d2 )) where r is the number of rounds.
19

Proof. From our solution to the problem under the `2 stability assumption, i.e. the proof of Theorem 9, we see that we do not need to know the entire matrix A to maintain an approximate inverse
provided that the number of low rank updates we need to perform is small. Therefore, it suffices
to show how to reduce the σ stability assumption with row addition and removal to the low rank
update problem mentioned in Theorem 9.
Our previous reduction under the σ stability assumption, i.e. Theorem 13, relied on two facts,
(1) that AT D(k−1) A ≈O(1) AT D(k) A for all k and that (2) the expected number of row changes is
O(r2 ). The first condition, (1) was used to ensure that we could efficiently compute approximate
leverage scores for AT D(k) A and the second condition (2) was used to ensure our invocation of
Theorem 9 was efficient. In the remainder of this proof we show how both condition (1) and (2)
can be achieved.
First we show how to achieve condition (1) under changes to d~(k) . First note that if the change
follows case (a) of the K stability of assumption then since log(d~(k) ) − log(d~(k−1) ) ∞ ≤ 0.1 clearly
AT D(k−1) A ≈O(1) AT D(k) A as before. On the other hand, for case (b) we claim that (1) can be
achieved by adding intermediate steps into the problem as follows. First, we can assume d~(k) ≥ d~(k−1)
or d~(k) ≤ d~(k−1) by splitting changes that corresponds to inserting or removing rows into two steps.
If we do the insertion first, it is easy to prove that γ −1 AT D(0) A  AT D(k) A  γAT D(0) A still
holds for some γ = O(β) after the splitting. For the insertion case, we consider the matrix



H(α) = AT D(k−1) + α D(k) − D(k−1) A.
Since D(k) − D(k−1)  0, it is easy to see that H(α) ≈O(1) H(2α). Since H(1) ≈poly(n) H(0), as

we have just argued, we have AT D(k) − D(k−1) A ≈poly(n) H(0) and hence H(poly(1/n)) ≈O(1)
H(0). Therefore, one can split the insertion case into O(log(n)) steps to ensure AT D(k−1) A ≈O(1)
AT D(k) A. The argument for removal is analogous.
Next we show how to achieve condition (2). Again we consider the case where a round is updated
by (a) and by (b) separately. For (a), the changes of d~ is small in σ norm, we already know from
our previous analysis of the σ stability assumption that it can only cause O(r2 ) many rows change.
For case (b) we define α be the total number of row changes due to (b). From the definition of
Algorithm 3, we see that there are two causes of resampling, namely, the change of d~ and the change
of ~σ . By the assumption (b) we know that there are at most K rows changes in d~ and therefore,
this contributes at most rK to α. For the change of ~σ , from first half of the proof of Lemma 15, we
know that the number of resampling that causes changes in ~h is bounded by
n
o
n
o
X
max σi (d~(k) ), σi (d~(k+1) ) min ln σi (d~(k) ) − ln σi (d~(k+1) ) , 1 .
i

Therefore, we have that
o
n
o
n
X
α ≤ rK +
max σi (d~(k) ), σi (d~(k+1) ) min ln σi (d~(k) ) − ln σi (d~(k+1) ) , 1 .
i,k

For any k, we define
Sk =



1
i such that ≤ σi (d~(k+1) )/σi (d~(k) ) ≤ 2 .
2

20

(6.1)

Then, we have that
σi (d~(k) ) − σi (d~(k+1) ) + O (1)

XX

k i∈Sk

k i∈S
/ k

≤ rK + O (1)

XX

σi (d~(k) ) − σi (d~(k+1) ) + O (1)

XX

≤ rK + O (1)

X

α ≤ rK + O (1)

XX

k i∈Sk

n
o
max σi (d~(k) ), σi (d~(k+1) )
σi (d~(k+1) ) − σi (d~(k) )

k i∈S
/ k

σi (d~(k) ) − σi (d~(k+1) ) .

i,k

Again, we can split the second case into insertion only and removal only. For the insertion case, we
P
P
(k+1)
(k)
note that σi (d~(k+1) ) ≤ σi (d~(k) ) for all i with d~i
= d~i . Since
σi (d~(k+1) ) = σi (d~(k) ) = d, we
have that
X
X
σi (d~(k) ) − σi (d~(k+1) ) = 2
σi (d~(k) ) − σi (d~(k+1) )
i

(k+1) ~(k)
d~i
6=di

≤ 4

n
o
(k+1)
(k)
i such that d~i
6= d~i
.

We have the same bound for the removal case. Putting it into (6.1), we have
α ≤ rK + O (1)

X

(k+1)
(k)
{i such that d~i
6= d~i } ≤ O(rK) .

k

This proves that the total expected number of row changes is bounded by O(r2 + rK). Note that
we use the fact that if there is an update triggered by case (a) after the update is triggered by case
(b) then we account for it by case (b), i.e. the interactions between case (a) and (b) only increase
the number of row changes by a constant.
Now, we can apply Theorem 9 and restart every r = dω−2 steps and get that the average cost
of maintenance is
!
 2
ω−2
r
+
rK
Õ d2 + (r2 + rK)d
+ (r2 + rK)ω
r


= Õ d2 + dr2(ω−1)−(ω−2) + rω−1−(ω−2) K ω−1 + r2ω + rω K ω


= Õ d2 + d1+(ω−2)ω + dω−2 K ω−1 + d2ω(ω−2) + dω(ω−2) K ω
√
Now using that ω ≤ 1 + 2 we have that 1 + (ω − 2)ω ≤ 2 and 2ω(ω − 2) ≤ 2 and therefore the
average cost of maintenance is


2
ω−2 ω−1
ω(ω−2) ω
Õ d + d
K
+d
K
Hence, as long as K ≤ d(3−ω)/(ω−1) , we have that the average cost is Õ(d2 ).
Remark 20. Using ω < 2.37287 [7], the above theorem shows how to solve inverse maintenance
problem with d0.4568 rows addition and removal in amortized Õ(d2 ) time.

21

7

Applications

In this section, we provide multiple applications of our algorithm for solving the inverse maintenance
problem under the `2 and σ stability assumptions. First in Section 7.1 we show how to use our results
to solve a linear program. Then in Section 7.2, Section 7.4, and Section 7.5 we state the consequences
of this linear programming algorithm for solving the non-linear regression, multicommodity flow,
and minimum cost flow respectively. In Section 7.3 we show how our results lead to new efficient
algorithms for computing a rounding of a polytope.

7.1

Linear Programming

Here we show how to use our solution to the inverse maintenance problem under the σ stability
assumption (See Section 4) to improve the running time of solving a linear program. Our main
result is the following:
Theorem 21. Let A ∈ Rd×n , ~c, ~l, ~u ∈ Rn , and ~b ∈ Rd for d ≤ n and suppose ~x ∈ Rn is an interior
point of the polytope
n
o
S = ~x ∈ Rn : A~x = ~b and li ≤ xi ≤ ui for all i ∈ [n]
.


def
−~l
~
u−~l
~l , ~c
. Then, consider the linear program
Let U = max ~u~u−~
~
u
−
,
,
~
x
∞
∞
∞

~
x−l ∞

def

OPT = min ~cT ~x
~
x∈S

In time Õ

√

.

(7.1)



d nnz(A) + d2 log (U/) , we can compute ~y such that A~y − ~b

AT S−2 A

≤ , li ≤

xi ≤ ui and ~cT ~y ≤ OP T +  where S is a diagonal matrix Sii = min (ui − yi , yi − li ).
Proof. In [18, ArXiv v2, Thm 28], we showed how to solve linear program of the form (7.1) by
solving a sequence of slowly changing linear systems AT D(k) A~x = ~q(k) where D(k) is a diagonal
matrix corresponding to a weighted distance of ~xk to the boundary of the polytope. In [18, ArXiv
v2, Lem 32] we showed that this sequence of linear systems satisfied the σ stability assumption.
Furthermore, we showed √
that it suffices to solve these linear systems to 1/poly(n) accuracy. Since,
the algorithm consists of d log(U/) rounds of this algorithm plus additional Õ(nnz(A) + d2 ) time
per round we have the desired result.
We remark that we can only output an almost feasible point but it is difficult to avoid because
finding any point ~x such that A~x = ~b takes O(ndω−1 ) which is slower than our algorithm when
n  d. Similarly, we have an algorithm for the dual of (7.1) as follows:
Theorem 22. Let A ∈ Rn×d where d ≤ n. Suppose we have an initial point ~x ∈ Rn such that
AT ~x = ~b and −1 < xi < 1. Then, we can find ~y ∈ Rd such that


~bT ~y + A~y + ~c ≤ min ~bT ~y + A~y + ~c
+ .
(7.2)
1
1
~
y

in Õ

√




d nnz(A) + d2 log (U/) time where U = max

2
1−~
x ∞,

2
~
x+1 ∞ ,

~c


∞

.

Proof. It is same as Theorem 21 except we invoke [18, ArXiv v2, Thm 29].
Remark 23. The existence of the interior point in Theorem 22 certifies the linear program has
bounded optimum value. Standard tricks can be used to avoid requiring such interior point but
may yield a more complicated looking running time (see Appendix E of [17] for instance).
22

`1 and `∞ Regression

7.2

The `p regression problem involves finding a vector ~x that minimize A~x − ~c p for some n × d
matrix A and some vector ~c. Recently, there has been much research [25, 26, 21, 2, 3] on solving
overdetermined problems (i.e. n  d) as these arises naturally in applications involving large
datasets. While there has been recent success on achieving algorithms whose running time is nearly
linear input plus something polynomial in d, in the case that p 6= 2 these algorithms achieve a
polynomial dependence on the desired
accuracy  [4]. Here we show how to improve the dependence
√
on  by paying a multiplicative d factor and improve upon previous algorithms in this regime.
Corollary 24. Let A ∈ Rn×d , ~c ∈ Rn , and p = 1 or p = ∞. There is an algorithm to find ~x such
that
A~x − ~c p ≤ min A~y − ~c p +  ~c p
~
y

in Õ

√


d nnz(A) + d2 log −1



time.

Proof. The `1 case is the special case of Theorem 22 with ~b = ~0 and an explicit initial point ~0.
For the `∞ case, we consider the following linear program


1
(7.3)
t + A~x − ~c − t~1 1 + A~x − ~c + t~1 1 .
min −2n +
2
~
x,t
where ~1 ∈ Rn is the all ones vector. For all t > 0, we have


if a ≥ t
2t + (a − t)
|a − t| + |a + t| = 2t
if − t ≤ a ≤ t .


2t + (−t − a) if a ≤ −t
Letting dist(a, [−t, t]) denote the distance from a to the interval [−t, t] (and |a| + |t| if t ≤ 0) we
then have that |a − t| + |a + t| = dist(a, [−t, t]) + 2t and consequently the linear program (7.3) is
equivalent to
n
X
def t
min f (~x, t) = +
dist ([A~x − c]i , [−t, t]) .
2
~
x,t
i=1

Note that when t ≥ A~x − ~c ∞ we have f (~x, t) = 2t
1
x − ~c ∞ . Consequently, the linear program (7.3)
2 A~
we rewrite it as follows



1
A
t+
min −2n +
A
2
~
x,t
Since


we can use

2n− 21
y
2n ~

A −~1
A ~1

T


~y =

~0d
−2n

and when t ≤ A~x − ~c ∞ we have f (~x, t) ≥
is equivalent to `∞ regression. To solve (7.3),
−~1
~1




for

~x
t

def

~y =


− ~c

.

(7.4)

1



~1
−~1


,

as the initial point for Theorem 22 and apply it to (7.3) to find ~x, t as desired.

Remark 25. We wonder if it is possible to obtain further running time improvements for solving `p
regression when p ∈
/ {1, 2, ∞}.

23

7.3

Õ(d) Rounding Ellipsoid for Polytopes

For any convex set K, we call an ellipsoid E is an α-rounding if E ⊂ K ⊂ αE. It is known that every
convex set in d dimension has a d-rounding ellipsoid and that such rounding have many applications.
(See [14, 35]) For polytopes K = {~x ∈ Rd : A~x ≥ ~b}, the previous best algorithm for finding a
3.5 −1
(1 + )d rounding takes time
and Õ(nd2 −1 ) [15]. Here we show how to compute an
√ Õ(n  ) [14]
2
O(d) rounding in time Õ(√ d(nnz(A) + d ) log(U )), which is enough for many applications. Note
that this is an at least O( d) improvement over previous results and for the case A is sparse, this
is an O(d1.5 ) improvement.
We split our proof into two parts. First we provide a technical lemma, Lemma 26, which shows
conditions under which we a point in a polytope is the center of a suitable ellipse and in Theorem 27
we show how to us this Lemma to compute the ellipse.
def
Lemma 26. Let P = {~x ∈ Rd : A~x ≥h ~b} be ai polytope for A ∈ Rn×d and d ≤ n. Furthermore, let
P
def
w
~ ∈ Rn>0 and let p(~x) = − ni=1 wi ln A~x − ~b for all ~x ∈ P . Now let ~x∗ = arg min~y∈Rd p(~x) and

i

suppose that we have γ ≥ 1 such that (A~h)i /(A~x∗ − ~b)i ≤ γ ~h
ellipsoid E = {~h ∈ Rd : ~h ∇2 p(~x∗ ) ≤ 1} satisfies
~x +

1
E ⊂ P ⊂ ~x + γ w
~
γ

1

for all ~h ∈ Rd . Then the

∇2 p(~
x∗ )

E.

(7.5)

Proof. For any ~z ∈ γ1 E, by the assumptions, we have that
(A~z)i
≤ γ ~z
(A~x∗ − ~b)i

max
i∈[n]

∇2 p(~
x∗ )

≤ 1.

Consequently, for all ~z ∈ γ1 E we have A(~x∗ + ~z) ≥ ~b and therefore ~x∗ + γ1 B∇2 p(~x∗ ) ⊂ P .
def
def
To prove the other side of (7.5) let ~x ∈ P be arbitrary and let ~s = A~x − ~b, ~s∗ = A~x∗ − ~b,
def
def
S = diag(~s), S∗ = diag(~x∗ ), and W = diag(w).
~ By the optimality conditions on ~x∗ we have
~
AWS−1
x∗ ) = ~0 .
∗ 1 = ∇p(~
Consequently,
0 = ~1T WS−1
x∗ − ~x) = ~1T WS−1
s∗ − ~s) = 0
∗ A (~
∗ (~
and therefore as ~1T WS−1
s∗ = w
~
∗ ~
n
X
i=1

wi

1

we have ~1T WS−1
s~x = w
~
∗ ~

[~s∗ − ~s]2i
= S−1
s∗
∗ ~
[~s∗ ]2i

2
W

−1
− 2~sT∗ S−1
s + S−1
s
∗ WS∗ ~
∗ ~

= S−1
s~x
∗ ~

2
W

− w
~

≤ WS−1
s~x
∗ ~

1

S−1
s
∗ (~

− S∗ )

∞

≤ w
~

· γ ~x − ~x∗
1

2
W

(Expanding the quadratic)
(Using ~1T WS−1
s= w
~
∗ ~

∞

1

and therefore

1

S−1
s~x
∗ ~

= w
~

1

(7.6)

− w
~

)

(Cauchy Schwarz)
(Using ~1T WS−1
s= w
~ )
∗ ~

1

1

(Assumption)

∇2 p(~
x∗ )

−1
Now, ∇2 p(~x∗ ) = AT S−1
∗ WS∗ A and therefore

1

[~s∗ −~s]2i
i=1 wi [~s∗ ]2i

Pn

= ~x − ~x∗

sides of the above equation by ~x − ~x∗ ∇2 p(~x∗ ) yields that ~x − ~x∗
right hand side of (7.5) and completing the proof.
24

∇2 p(~
x∗ )

2
.
∇2 p(~
x∗ )

≤γ w
~

Dividing both
1

yielding the

Theorem 27. Let A ∈ Rn×d for d ≤ n and suppose we have an initial point ~x0 ∈ Rd such that
d
~ ⊂ 100d · E in time
A~x0 ≥ ~b. Then, we can find an ellipsoid E such
n that E ⊂ {~x ∈ R : A~x ≥ b} o
√

1
Õ( d nnz(A) + d2 log(U )) where U = max maxA~x≥~b A~x − ~b ∞ ,
.
~
A~
x0 −b ∞

Proof. Given any interior point ~x0 such that A~x0 > ~b in [17] we showed how to find a weight w
~
such that w
~ 1 ≤ 3d and for any ~h ∈ Rd we have
max
i∈[n]

where p(~x) = −

[A~h]i
≤ 3 ~h
[A~x∗ − ~b]i

(7.7)

∇2 p(~
x∗ )

h
i
~
and ~x∗ = arg minA~y≥~b p(~y ). Furthermore, we showed how to
w
ln
A~
x
−
b
i=1 i

Pn

i

find ~x ∈ Rd such that

n
X

1
[A(~x − ~x∗ )]2i
≤
.
(7.8)
2
100
(A~x∗ − ~b)i
i=1
√

As discussed in Theorem 21, this can be done in time Õ( d nnz(A) + d2 log(U )). Consequently,
~ 1 = 3d and this gives us an ellipsoid E such that
we can use Lemma 26 with γ = 3 and w
wi

1
~x∗ + E ⊂ P ⊂ ~x∗ + 9dE
3
where E = {~h ∈ Rd : ~h

∇2 p(~
x∗ )

≤ 1}.

Now, we show how to approximate E using ~x. (7.8) shows that that ~x − ~x∗

∇2 p(~
x∗ )

≤

1
10 ,

1
E and hence
therefore, we have ~x − ~x∗ ∈ 10




1
1
1
1
~x +
−
E ⊂ ~x∗ + E ⊂ P ⊂ ~x∗ + 9dE ⊂ ~x + 9d +
E.
3 10
3
10
A(~
x−~
x∗ )
A~
x∗ −~b i

Now, using (7.8) and (7.7), we have

≤

3
10 ,

we have

1
1
2
−1
2
x)  ∇2 p(~x∗ ) = AT S−1
x) .
∗ WS∗ A 
3 2 ∇ p(~
3 2 ∇ p(~
(1 + 10 )
(1 − 10
)
Therefore, we have
~x +
where E 0 = {~h ∈ Rd : ~h

7.4

∇2 p(~
x)

1
3

1
1
− 10
9d + 10
0
0
E
⊂
P
⊂
~
x
+
3
3 E
1 + 10
1 − 10

≤ 1}. After rescaling the ellipsoid, we have the result.

Multicommodity Flow

Here we show how our algorithm can be used to improve the running time for solving multicommodity flow. Note that the result presented here is meant primarily to illustrate our approach, we
believe it can be further improved using the techniques in [10, 18].
For simplicity, we focus on the maximum concurrent flow problem. In this problem, we are
given a graph G = (V, E) k source sink pairs (si , ti ) ∈ RV ×V , and capacities ~c ∈ RE and wish
to compute the maximum α ∈ R such that we can simultaneously for all
Pik ∈ [k] route α unit of
E
flow fi ∈ R between si and ti while maintaining the capacity constraint i=1 |fi (e)| ≤ c(e) for all
25

e ∈ E. There are no combinatorial algorithms known and there are multiple algorithms to compute
a (1 − ) optimal flow in time polynomial in |E|, |V | and −1 [6, 11, 8, 23]. In this regime, the
fastest algorithm for directed graphs takes O((|E| + k)|V |−2 log U ) [23] and the fastest algorithm
for undirected graphs takes Õ(|E|1+o(1) k 2 −2 logO(1) U ) [13] where U is the maximum capacity. For
linear convergence,
the previous best algorithm is a specialized interior point method that takes
p
time O( |E|k|V |2 k 2 log(−1 )) [10].
To solve the concurrent multicommodity flow problem we use the following linear program

X

max

α

g(e)

=

X

fi (u, v)

=

0 for all vertices u ∈
/ {si , ti },

fi (si , v)

=

α for all i,

fi (e) for all edges e,

v∈V

X
v∈V

c(e) ≥ fi (e) ≥ 0 for all edges e,
c(e) ≥

g(e)

≥ 0 for all edges e.

Note that there are O(k|E|) variables, O(k|V | + |E|) equality constraints, and O(k|E|) non-zeroes
in the constraint matrix. Furthermore, it is easy to find an initial point by computing a shortest path for each si and ti and sending a small amount of flow along that path. Also, given
any almost feasible flow, one can make it feasible by scaling and send excess flow at every vertexpback to si along some spanning tree. Therefore, Theorem 21 gives an algorithm
that takes

2
2.5
Õ( |E| + k|V | k|E| + (|E| + k|V |) log(U/)) = Õ (|E| + k|V |) log(U/) time. Note that
this is faster than the previous best algorithm when k ≥ (|E|/|V |)0.8 .

7.5

Minimum Cost Flow

The minimum cost flow problem and the more specific, maximum flow problem, are two of the most
well studied problems in combinatorial optimization [31]. Many techniques have been developed for
solving this problem. Yet, our result matches the fastest algorithm for solving these problems on
dense graphs [18] without using any combinatorial structure of this problem, in particular, Laplacian
system solvers. We emphasize that this result is not a running time improvement, rather just a
demonstration of the power of our result and an interesting statement on efficient maximum flow
algorithms.
Corollary 28. There is an Õ(|V |2.5 logO(1) (U )) time algorithm to compute an exact minimum cost
maximum flow for weighted directed graphs with |V | vertices, |E| edges and integer capacities and
integer cost at most U .
Proof. The proof is same as [18, ArXiv v2, Thm 34,35] except that we use Theorem 21 to solve the
linear program. The proof essentially writes the minimum cost flow problem into a linear program
with an explicit interior point and shows how to round an approximately optimal solution to a
vertex of the polytope. To perform this rounding, we need to a fractional solution with error less
1
than O( poly(|V
|U ) ) and which yields the log(U ) term in the running time.

7.6

Convex Problems

Many problems in convex optimization can be efficiently reduced to the problem of finding a point
in a convex set K given a separation oracle. Recall that given a point ~x, the separation oracle
26

either outputs that ~x is in K or outputs a separating hyperplane that separates the input point ~x
and the convex set K. In [20], they showed that how to make use of our fast algorithms for inverse
maintenance problem under the K stability assumption to obtain the following improved running
time for this fundamental problem:
Theorem 29 ([20]). Given a non-empty convex set K ⊆ Rd that is contained in a box of radius R,
i.e. maxx∈K x ∞ ≤ R. We are also given a separation oracle for K that takes O(T ) time for each
call. For any 0 <  < R, we can either finds ~x ∈ K or proves that K does not contains a ball with
radius  in time O(dT log(dR/) + d3 logO(1) (dR/)).
Remark 30. [20] uses Theorem 19 to solve the linear systems involved in their cutting plane method.
However, their inverse maintenance problem satisfies the `2 stability assumption with 1 rows addition
and removal. Furthermore, the A involved has only O(d) many rows. Therefore, we believe a simple
variant of [33] may also suffice for that paper.

8

Open Problem: Sampling from a Polytope

Sampling a random point in convex sets is a fundamental problem in convex geometry with numerous
applications in optimization, counting, learning and rounding [35]. Here consider a typical case
where the convex set is a polytope that is explicitly given as {~x ∈ Rd : A~x ≥ ~b} for A ∈ Rn×d .
The current fastest algorithm for this setting is Hit-and-Run [22] and Dikin walk [9].
Given an initial random point in the set, Hit-and-Run takes O∗ (d3 ) iterations and each iteration
takes time O∗ (nnz(A)) while Dikin walk takes O∗ (nd) iterations and each iteration takes time
O∗ (ndω−1 ) where the O∗ notation omits the dependence on error parameters and logarithmic terms.
Each iteration of Dikin walk is expensive because it solves a linear system to obtain the next point
and computes determinants to implement an importance sampling scheme. The linear systems can
be solved in amortized cost O∗ (nnz(A) + d2 ) by the inverse maintenance machinery presented in
this paper. Unfortunately it is not known how to use this machinery to efficiently compute the
determinant to sufficient accuracy to suffice for this method.
We leave it as an open problem whether it is possible to circumvent this issue and improve the
running time of a method like the Dikin walk. In an older version of this paper, we mistakenly
claimed an improved running time for Dikin walk by noting solely the improved running time for
linear system solving and ignoring the determinant computation. We thank Hariharan Narayanan
for pointing out this mistake.

9

Acknowledgments

We thank Yan Kit Chim, Andreea Gane, and Jonathan A. Kelner for many helpful conversations.
This work was partially supported by NSF awards 0843915 and 1111109, NSF Graduate Research
Fellowship (grant no. 1122374). Part of this work was done while both authors were visiting the
Simons Institute for the Theory of Computing, UC Berkeley.

References
[1] Hui Han Chin, Aleksander Madry, Gary L Miller, and Richard Peng. Runtime guarantees
for regression problems. In Proceedings of the 4th conference on Innovations in Theoretical
Computer Science, pages 269–282. ACM, 2013.

27

[2] Kenneth L Clarkson and David P Woodruff. Low rank approximation and regression in input
sparsity time. In Proceedings of the 45th annual ACM symposium on Symposium on theory of
computing, pages 81–90. ACM, 2013.
[3] Michael B. Cohen, Yin Tat Lee, Cameron Musco, Christopher Musco, Richard Peng, and Aaron
Sidford. Uniform sampling for matrix approximation. CoRR, abs/1408.5099, 2014.
[4] Michael B Cohen and Richard Peng. Lp row sampling by lewis weights. arXiv preprint
arXiv:1412.0588, 2014.
[5] Petros Drineas, Michael W Mahoney, and S Muthukrishnan. Subspace sampling and relativeerror matrix approximation: Column-based methods. In Approximation, Randomization, and
Combinatorial Optimization. Algorithms and Techniques, pages 316–326. Springer, 2006.
[6] Lisa K Fleischer. Approximating fractional multicommodity flow independent of the number
of commodities. SIAM Journal on Discrete Mathematics, 13(4):505–520, 2000.
[7] François Le Gall. Powers of tensors and fast matrix multiplication.
arXiv:1401.7714, 2014.

arXiv preprint

[8] Naveen Garg and Jochen Koenemann. Faster and simpler algorithms for multicommodity flow
and other fractional packing problems. SIAM Journal on Computing, 37(2):630–652, 2007.
[9] Ravindran Kannan and Hariharan Narayanan. Random walks on polytopes and an affine
interior point method for linear programming. Mathematics of Operations Research, 37(1):1–
20, 2012.
[10] Sanjiv Kapoor and Pravin M Vaidya. Speeding up karmarkar’s algorithm for multicommodity
flows. Mathematical programming, 73(1):111–127, 1996.
[11] George Karakostas. Faster approximation schemes for fractional multicommodity flow problems. ACM Transactions on Algorithms (TALG), 4(1):13, 2008.
[12] Narendra Karmarkar. A new polynomial-time algorithm for linear programming. In Proceedings
of the sixteenth annual ACM symposium on Theory of computing, pages 302–311. ACM, 1984.
[13] Jonathan A Kelner, Yin Tat Lee, Lorenzo Orecchia, and Aaron Sidford. An almost-linear-time
algorithm for approximate max flow in undirected graphs, and its multicommodity generalizations. In SODA, pages 217–226. SIAM, 2014.
[14] Leonid G Khachiyan. Rounding of polytopes in the real number model of computation. Mathematics of Operations Research, 21(2):307–320, 1996.
[15] Piyush Kumar and E Alper Yildirim. Minimum-volume enclosing ellipsoids and core sets.
Journal of Optimization Theory and Applications, 126(1):1–21, 2005.
[16] Béatrice Laurent and Pascal Massart. Adaptive estimation of a quadratic functional by model
selection. Annals of Statistics, pages 1302–1338, 2000.
[17] Yin Tat Lee and Aaron Sidford. Path finding i: Solving linear programs with \˜ o (sqrt(rank))
linear system solves. arXiv preprint arXiv:1312.6677, 2013.
[18] Yin Tat Lee and Aaron Sidford. Path finding ii: An\˜ o (m sqrt (n)) algorithm for the minimum
cost flow problem. arXiv preprint arXiv:1312.6713, 2013.
28

[19] Yin Tat Lee and Aaron Sidford. Path-finding methods for linear programming : Solving linear
programs in õ(sqrt(rank)) iterations and faster algorithms for maximum flow. In 55th Annual
IEEE Symposium on Foundations of Computer Science, FOCS 2014, 18-21 October, 2014,
Philadelphia, PA, USA, pages 424–433, 2014.
[20] Yin Tat Lee, Aaron Sidford, and Sam Chiu-wai Wong. A faster cutting plane method and its
implications for combinatorial and convex optimization. In Foundations of Computer Science
(FOCS), 2015 IEEE 56th Annual Symposium on. IEEE, 2015.
[21] Mu Li, Gary L Miller, and Richard Peng. Iterative row sampling. 2012.
[22] László Lovász and Santosh Vempala. Fast algorithms for logconcave functions: Sampling,
rounding, integration and optimization. In Foundations of Computer Science, 2006. FOCS’06.
47th Annual IEEE Symposium on, pages 57–68. IEEE, 2006.
[23] Aleksander Madry. Faster approximation schemes for fractional multicommodity flow problems
via dynamic graph algorithms. In Proceedings of the forty-second ACM symposium on Theory
of computing, pages 121–130. ACM, 2010.
[24] Michael W. Mahoney. Randomized algorithms for matrices and data. Foundations and Trends
in Machine Learning, 3(2):123–224, 2011.
[25] Michael W Mahoney, Petros Drineas, Malik Magdon-Ismail, and David P Woodruff. Fast
approximation of matrix coherence and statistical leverage. In ICML, 2012.
[26] Jelani Nelson and Huy L Nguyên. Osnap: Faster numerical linear algebra algorithms via sparser
subspace embeddings. arXiv preprint arXiv:1211.1002, 2012.
[27] Yu Nesterov and Arkadi Nemirovskiy. Self-concordant functions and polynomial-time methods
in convex programming. USSR Academy of Sciences, Central Economic & Mathematic Institute,
1989.
[28] Yu Nesterov and A Nemirovsky. Acceleration and parallelization of the path-following interior point method for a linearly constrained convex quadratic problem. SIAM Journal on
Optimization, 1(4):548–564, 1991.
[29] Yurii Nesterov and Arkadii Semenovich Nemirovskii. Interior-point polynomial algorithms in
convex programming, volume 13. Society for Industrial and Applied Mathematics, 1994.
[30] James Renegar. A polynomial-time algorithm, based on newton’s method, for linear programming. Mathematical Programming, 40(1-3):59–93, 1988.
[31] Alexander Schrijver.
Springer, 2003.

Combinatorial optimization:

polyhedra and efficiency, volume 24.

[32] Daniel A Spielman and Nikhil Srivastava. Graph sparsification by effective resistances. SIAM
Journal on Computing, 40(6):1913–1926, 2011.
[33] Pravin M Vaidya. Speeding-up linear programming using fast matrix multiplication. In Foundations of Computer Science, 1989., 30th Annual Symposium on, pages 332–337. IEEE, 1989.
[34] Pravin M Vaidya. A new algorithm for minimizing convex functions over convex sets. Mathematical Programming, 73(3):291–341, 1996.
29

[35] Santosh S Vempala. Recent progress and open problems in algorithmic convex geometry. In
LIPIcs-Leibniz International Proceedings in Informatics, volume 8. Schloss Dagstuhl-LeibnizZentrum fuer Informatik, 2010.
[36] Virginia Vassilevska Williams. Multiplying matrices faster than coppersmith-winograd. In
Proceedings of the forty-fourth annual ACM symposium on Theory of computing, pages 887–
898. ACM, 2012.

A

Relationships Between T -time Linear Solver and Inverse Matrix

In this section, we prove relationships between linear T -time solvers and inverse matrices. In
Lemma 31 we show that obtaining a spectral approximation to the inverse of a matrix suffices to
obtain a linear solver. In Lemma 32, we show how a linear solver of a matrix yields a spectral
approximation to the inverse of that matrix.
Lemma 31. Given a PD matrix M and a symmetric matrix N ≈O(1) M−1 that can be applied to
a vector in O(T ) time we have a linear (nnz(M) + T )-time solver S of M.
Proof. Since N ≈O(1) M−1 there is a constant L such that L1 N−1  M  LN. Consider the
def
algorithm S(~b, ) = Q~b where

z 
X
MN k
def 1
I−
Q = N
L
L
k=0

1
2

log(L−4 )/ log(1

L−2 ).

where z =
−
This is the linear operator corresponding to performing z
steps of preconditioned gradient descent to solve the linear system in M. Clearly, applying Q can
be done in time O((nnz(A) + T )ze ) = O((nnz(A) + T ) log(−1 )). Therefore, all that remains is to
show that S(~b, ) is a solver for M.
Note that we can rewrite Q equivalently as
k
z 
1 1/2 X
1 1/2
1/2
Q = N
I − N MN
N1/2
L
L
k=0

and hence it is symmetric. Furthermore, since clearly
−1

M

1
I
L2



1 1/2
MN1/2
LN

 I we have that



−1
k
z 
1 1/2
1 1/2 X
1 1/2
1 1/2
1/2
1/2
1/2
− Q = N
I − I − N MN
N − N
I − N MN
N1/2
L
L
L
L
k=0

k
∞
X
1
1
= N1/2
I − N1/2 MN1/2 N1/2 .
L
L
k=z +1

Using the above two inequalities, we have that
S(~b, ) − M−1~b

2
M

=


k
∞
1 1/2 X
1 1/2
1/2
N
I − N MN
N1/2~b
L
L
k=z +1

≤

1
L

k

∞
X
1 1/2
1/2
I − N MN
N1/2~b
L

k=z +1

M
2

2

30

2


and since using that 0  I − L1 N1/2 MN1/2  1 − L12 I, we have
k



∞
X
1 z
1 1/2
1/2
2
0
I − N MN
 L 1− 2
I .
L
L
k=z +1

Consequently
S(~b, ) − M−1~b

2
M



1 2z
≤ L3 1 − 2
N1/2~b
L

2
2



1 2z
≤ L4 1 − 2
M−1~b
L

2
M

.

Thus, we see that z was chosen precisely to complete the proof.
Lemma 32. Let S be a linear solver of M such that S(~b, ) = Q~b for  ∈ [0, 0.1]. Then, we have
QT MQ ≈4√ M−1 and Q MQT ≈4√ M−1
Proof. By the definition of S and standard inequalities, we have that for all ~b,


√ 
1
2
2
2
~
S(~b, ) − M−1~b M + 1 +  M−1~b M
S(b, ) M ≤
1+ √

√
√
√ 
2
2
≤  +  + 1 +  M−1~b
≤ (1 + 3 ) M−1~b
M

M

On the another hand, we have that for all ~b


√ 
1
2
2
M−1~b M ≤
1+ √
S(~b, ) − M−1~b M + 1 +  S(~b, )

√ 
√
2
2
+ 1 +  S(~b, ) .
≤ 2  M−1~b
M

2
M

M

Combining these two inequalities and using the definition of S, we see that for all ~b we have

√ 
√ 
1 − 2  ~ T −1~ ~ T T
√
b M b ≤ b Q MQ~b ≤ 1 + 3  ~bT M−1~b
1+ 
Using a Taylor expansion of e and the definition of ≈ then yields QT MQ ≈4√ M−1 . In other
√
√
T

words, all eigenvalues of M1/2 Q M1/2
M1/2 Q M1/2 lies between e−4  and e4  . In general,
for any square matrix B, the set of eigenvalues of BT B equals to the set of eigenvalues of BBT

T
because of the SVD decomposition. Hence, all eigenvalues of M1/2 Q M1/2 M1/2 Q M1/2 also
lies in the same region. This proves Q MQT ≈4√ M−1 .

B

Remarks for Figure 1.1

Figure 1.1 shows the previous fastest algorithm for linear programming minA~x≥~b ~cT ~x where A ∈
Rn×d . The running time described in the figure comes from the following algorithms:

√
√
1. Õ nz + nd2 + n1.34 d1.15 is achieved by the interior point method of Vaidya [33].

√
2. Õ nz + nd1.38 + n0.69 d2 is achieved by using the Karmarkar acceleration scheme [12] on
the short step path following method. See [28] for the details.
√
3. Õ( n(z + d2.38 )) is achieved by using the currently best linear system solver [26, 21, 7, 3] on
the short step path following method [30].
4. Õ(d(z + d2.38 )) is achieved by the cutting plane method of Vaidya [34].

31

