Hardness Results for Structured Linear Systems
Rasmus Kyng∗
Yale University
rasmus.kyng@yale.edu

Peng Zhang†
Georgia Institute of Technology
pzhang60@gatech.edu

arXiv:1705.02944v2 [cs.CC] 24 Sep 2017

September 26, 2017

Abstract
We show that if the nearly-linear time solvers for Laplacian matrices and their generalizations
can be extended to solve just slightly larger families of linear systems, then they can be used
to quickly solve all systems of linear equations over the reals. This result can be viewed either
positively or negatively: either we will develop nearly-linear time algorithms for solving all
systems of linear equations over the reals, or progress on the families we can solve in nearlylinear time will soon halt.

∗
†

Supported by ONR Award N00014-16-1-2374.
Partially supported by the NSF under Grant No. 1637566.

Contents
1 Introduction
1.1 Our Results . . . . . . . . . . . . . . . . . . . . .
1.2 Approximately Solving Linear Systems and Normal
1.3 Graph-Structured Block Matrices . . . . . . . . . .
1.4 Our Reduction: Discussion and an Example . . . .

. . . . . . .
Equations .
. . . . . . .
. . . . . . .

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

2
3
4
5
7

2 Preliminaries
9
2.1 Approximately Solving A Linear System . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2 Measuring the Difficulty of Solving a Linear System . . . . . . . . . . . . . . . . . . 11
2.3 Matrix Classes and Reductions Between Them . . . . . . . . . . . . . . . . . . . . . 11
3 Main Results
3.1 Outline of Remaining Sections

12
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

4 Reducing Zero-Sum Power Two Linear Systems to Two-Commodity Linear Systems
4.1 Reduction Between Exact Solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Relationship Between Schur Complements . . . . . . . . . . . . . . . . . . . . . . . .
4.3 Approximate solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4 Bounding Condition Number of the New Matrix . . . . . . . . . . . . . . . . . . . .
4.5 Putting it All Together . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13
17
21
22
24
30

5 MC2 Efficiently Reducible to MC2>0
5.1 Construction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 Bounding Condition Number of the New System in MC2>0 . . . . . . . . . . . . . . .
5.3 Putting it All Together . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31
31
36
37

6 Rounding and Scaling Weights to Integers

37

7 G Efficiently Reducible to Gz,2
41
7.1 G ≤f Gz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
7.2 Gz ≤f Gz,2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
8 2D Trusses

52

9 Connections with Interior Point Methods
9.1 Brief Overview of Interior Point Methods . . . . . . . . . . . . . . . . . . . . . . . .
9.2 2-Commodity Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9.3 Isotropic Total Variation Minimization . . . . . . . . . . . . . . . . . . . . . . . . . .

56
56
56
63

10 A Natural Decision Problem: Vector in the Image of a Matrix
65
10.1 Why Not Solve the Exact Decision Problem? . . . . . . . . . . . . . . . . . . . . . . 65
A Linear Algebra Background

67

B Using Complexity Parameters

67

1

1

Introduction

We establish a dichotomy result for the families of linear equations that can be solved in nearlylinear time. If nearly-linear time solvers exist for a slight generalization of the families for which
they are currently known, then nearly-linear time solvers exist for all linear systems over the reals.
This type of reduction is related to the successful research program of fine-grained complexity,
such as the result [WW10] which showed that the existence of a “truly subcubic” time algorithm for
All-Pairs Shortest Paths Problem is equivalent to the existence of “truly subcubic” time algorithm
for a wide range of other problems. For any constant a ≥ 1, our result establishes for 2-commodity
matrices, and several other classes of graph structured linear systems, that we can solve a linear
e a ) if and only if we can solve linear
system in a matrix of this type with s nonzeros in time O(s
e a ).
systems in all matrices with polynomially bounded integer entries in time O(s
n×n
n
In the RealRAM model, given a matrix A ∈ R
and a vector c ∈ R , we can solve the linear
ω
system Ax = c in O(n ) time, where ω is the matrix multiplication constant, for which the best
currently known bound is ω < 2.3727 [Str69, Wil12]. Such a running time bound is cost prohibitive
for the large sparse matrices often encountered in practice. Iterative methods [Saa03], first order
methods [BV04], and matrix sketches [Woo14] can all be viewed as ways of obtaining significantly
better performance in cases where the matrices have additional structure.
In contrast, when A is an n × n Laplacian matrix with m non-zeros, and polynomially
bounded entries, the linear system Ax = c can be solved approximately to -accuracy in
O((m + n) log1/2+o(1) n log(1/)) time [ST14, CKM+ 14]. This result spurred a series of major developments in fast graph algorithms, sometimes referred to as “the Laplacian Paradigm” of designing
graph algorithms [Ten10]. The asymptotically fastest known algorithms for Maximum Flow in directed unweighted graphs [Mad13, Mad16], Negative Weight Shortest Paths and Maximum Weight
Matchings [CMSV17], Minimum Cost Flows and Lossy Generalized Flows [LS14, DS08] all rely on
fast Laplacian linear system solvers.
The core idea of the Laplacian paradigm can be viewed as showing that the linear systems that
arise from interior point algorithms, or second-order optimization methods, have graph structure,
and can be preconditioned and solved using graph theoretic techniques. These techniques could
potentially be extended to a range of other problems, provided fast solvers can be found for the
corresponding linear systems. Here a natural generalization is in terms of the number of labels per
vertex: graph Laplacians correspond to graph labeling problems where each vertex has one label,
and these labels interact pairwise via edges. Multi-label variants of these exist in Markov random
fields [SZS+ 08], image processing [WYYZ08], Euclidean embedding of graphs [CLM+ 16], data
processing for cryo-electron microscopy [SS11, SS12, ZS14], phase retrieval [ABFM14, MTW14],
and many image processing problems (e.g. [OSB15, ANKKS+ 12]). Furthermore, linear systems
with multiple labels per vertex arise when solving multi-commodity flow problems using primaldual methods. Linear systems related to multi-variate labelings of graphs have been formulated
as the quadratically-coupled flow problem [KMP12] and Graph-Structured Block Matrices [Spi16].
They also occur naturally in linear elasticity problems for simulating the effect of forces on truss
systems [DS07].
Due to these connections, a central question in the Laplacian paradigm of designing graph algorithms is whether all Graph-Structured Block Matrices can be solved in (approximately) nearlylinear time. Even obtaining subquadratic running time would be constitute significant progress.
There has been some optimism in this direction due to the existence of faster solvers for special cases:
nearly-linear time solvers for Connection Laplacians [KLP+ 16], 1-Laplacians of collapsible 3-D simplicial complexes [CFM+ 14], and algorithms with runtime about n5/4 for 2D planar truss stiffness
matrices [DS07]. Furthermore, there exists a variety of faster algorithms for approximating multi2

commodity flows to (1+) accuracy in time that scales as poly(−1 ) [GK98, Mad10, Fle00, LMP+ 91],
even obtaining nearly-linear running times when the graph is undirected [She13, KLOS14, Pen16].
The subquadratic variants of these routines also interact naturally with tools that in turn utilize
Laplacian solvers [KMP12, LMP+ 91]. These existing tight algorithmic connections, as well as the
solver for Connection Laplacians, and the fact that combinatorial preconditioners partly originated
from speeding up interior point methods through preconditiong Hessians [Vai89], together provide
ample reason to hope that one could develop nearly-linear time solvers for linear systems related
to multicommodity flows. Any algorithm that solves such systems to high accuracy in m1+α time
would in turn imply multicommodity flow algorithms that run in about n1/2 m1+α time [LS14],
while the current best running times are about n2.5 [LS15].
Unfortunately, we show that if linear systems in general 2D truss stiffness matrices or 2commodity Laplacians can be solved approximately in nearly-linear time, then all linear systems
in matrices with polynomially bounded integer entries can be solved in nearly-linear time. In fact,
we show in a strong sense that any progress made in developing solvers for these classes of matrices
will translate directly into similarly fast solvers for all matrices with polynomially bounded integer
entries. Thus developing faster algorithms for these systems will be as difficult as solving all linear
systems faster.
Since linear system solvers used inside Interior Point Methods play a central role in the Laplacian
paradigm for designing high-accuracy algorithms, this may suggest that in the high-accuracy regime
the paradigm will not extend to most problems that require multiple labels/variables per edge or
vertex. Alternatively, an algorithmic optimist might view our result as a road-map for solving all
linear systems via reductions to fast linear system solvers for Graph-Structured Block Matrices.

1.1

Our Results

Fast linear system solvers for Laplacians, Connection Laplacians, Directed Laplacians, and 2D
Planar Truss Stiffness matrices are all based on iterative methods and only produce approximate
solutions. The running time for these solvers scales logarithmically with the error parameter , i.e.
as log(1/). Similarly, the running time for iterative methods usually depends on the condition
number of the matrix, but for state-of-the-art solvers for Laplacians, Connection Laplacians, and
Directed Laplacians, the dependence is logarithmic. Consequently, a central open question is
whether fast approximate solvers exist for other structured linear systems, with running times
that depend logarithmically on the condition number and the accuracy parameter.
Integral linear systems are reducible to Graph-Structured Block Matrices. Our
reductions show that if fast approximate linear system solvers exist for multi-commodity Laplacians,
2D Truss Stiffness, or Total Variation (TV) Matrices, then fast approximate linear system solvers
exist for any matrix, in the very general sense of minimizing minx kAx − ck22 . Thus our result
also applies to singular matrices, where we solve the pseudo-inverse problem to high accuracy.
Theorem 1.1 gives an informal statement of our main result. The result is stated formally in
Section 3 as Theorem 3.2.
Theorem 1.1 (Hardness for Graph-Structured Linear Systems (Informal)). We consider three
types of Graph-Structured Block Matrices: Multi-commodity Laplacians, Truss Stiffness Matrices,
and Total Variation Matrices. Suppose that for one or more of these classes, the linear system
e a ), for some constant a ≥ 1,
Ax = c in a matrix A with s non-zeros can be solved in time O(s
with the running time having logarithmic dependence on condition number and accuracy1 . Then
1

This is the kind of running time guarantee established for Laplacians, Directed Laplacians, Connection Laplacians,
and bounded-weight planar 2D Truss Stiffness matrices.

3

linear systems in all matrices with polynomially bounded integer entries and condition number can
e a ), where again s is the number of non-zero entries of the
be solved to high accuracy in time O(s
matrix.
Our results can easily be adapted to show that if fast exact linear system solvers exist for multicommodity Laplacians, then exact solvers exist for all non-singular integer matrices. However, this
is of less interest since there is less evidence that would suggest we should expect fast exact solvers
to exist.
The notion of approximation used throughout this paper is the same as that used in the Laplacian solver literature (see Section 2.1). To further justify the notion of approximate solutions to
linear systems, we show that it lets us solve a natural decision problem for linear systems:
We show that deciding if a vector is approximately in the image of a matrix can be
reduced to approximately solving linear systems. We show this in Section 10. We also show
that the exact image decision problem requires working with exponentially small numbers, even
when the input has polynomially bounded integral entries and condition number. This means that
in fixed-point arithmetic, we can only hope to solve an approximate version of the problem. The
problem of approximately solving general linear systems can be reduced the problem of approximately solving Graph-Structured Block Matrix linear systems. Together, these statements imply
that we can also reduce the problem of deciding whether a vector is approximately in the image
of a general matrix to the problem of approximately solving Graph-Structured Block Matrix linear
systems.
We establish surprising separations between many problems known to have fast
solutions and problems that are as hard solving general linear systems. Our results
trace out several interesting dichotomies: restricted cases of 2D truss stiffness matrices have fast
solvers, but fast solvers for all 2D truss stiffness matrices would imply equally fast solvers for all
linear systems. TV matrices can be solved quickly in the anisotropic case, but in the isotropic case
imply solvers for all linear systems. Fast algorithms exist for multi-commodity problems in the low
accuracy regime, but existing approaches for the high accuracy regime seem to require fast solvers
for multi-commodity linear systems, which again would imply fast solvers for all linear systems.
Our reductions only require the simplest cases of the classes we consider: 2-Commodity
Laplacians are sufficient, as are (non-planar) 2D Truss Stiffness matrices, and Total Variation Matrices with 2-by-2 interactions. Linear systems of these three classes have many applications, and faster solvers for these would be useful in all applications. Trusses have been
studied as the canonical multi-variate problem, involving definitions such as Fretsaw extensions [ST08] and factor widths [BCPT05], and fast linear system solvers exist for the planar
2D case with bounded weights [DS07]. Total Variation Matrices are widely used in image denoising [CS05]. The anisotropic version can be solved using nearly-linear time linear system
solvers [KMT11], while the isotropic version has often been studied using linear systems for
which fast solvers are not known [GY04, WR07, CMMP13a]. Multi-commodity flow problems
have been the subject of extensive study, with significant progress on algorithms with low accuracy [GK98, Mad10, Fle00, LMP+ 91, She13, KLOS14, Pen16], while high accuracy approaches use
slow general linear system solvers.

1.2

Approximately Solving Linear Systems and Normal Equations

The simplest notion of solving a linear system Ax = c, is to seek an x s.t. the equations are
exactly satisfied. More generally, if the system is not guaranteed to have a solution, we can ask for
an x which minimizes kAx − ck22 . An x which minimizes this always exists. In general, it may
4

not be unique. Finding an x which minimizes kAx − ck22 is equivalent to solving the linear system
A> Ax = A> c, which is referred to as the normal equation for the linear system Ax = c (see
[TBI97]). The problem of solving the normal equations (or equivalently, minimizing kAx − ck22 ), is
a generalization of the problem of linear system solving, since the approach works when A is nonsingular, while also giving meaningful results when A is singular. The normal equation problem
can also be understood in terms of the Moore-Penrose pseudo-inverse of a matrix M , which is
denoted M † as x = (A> A)† A> c is a solution to the normal equations. Taking the view of linear
system solving as minimizing kAx − ck22 also gives sensible ways to define an approximate solution
to a linear system: It is an x that ensures kAx − ck22 is close to minx kAx − ck22 . In Section 2,
we formally define several notions of approximate solutions to linear systems that we will use
throughout the paper.
An important special case of linear systems is when the matrix of coefficients of the system is
positive semi-definite. Since A> A is always positive semi-definite, solving the normal equations
for a linear system falls into this case. Linear systems over positive semi-definite matrices can
be solved (approximately) by approaches known as iterative methods, which often lead to much
faster algoritms than the approaches used for general linear systems. Iterative methods inherently
produce approximate solutions2 .

1.3

Graph-Structured Block Matrices

Graph-Structured Block Matrices are a type of linear system that arise in many applications.
Laplacian matrices and Connection Laplacians both fall in this category.
Suppose we have a collection of n disjoint sets of variables X1 , . . . , Xn , with each set having
the same size, |Xi | = d. Let x i denote the vector3 of variables in Xi , and consider an equation of
the form S x i − T x j = 0, where S and T are both r × d matrices. Now we form a linear system
Bx = 0 by stacking m equations of the form given above as the rows of the system. Note that,
very importantly, we allow a different choice of S and T for every pair of i and j. This matrix
B ∈ Rmr×nd we refer to as a Incidence-Structured Block Matrix (ISBM), while we refer to B > B
as a Graph-Structured Block Matrix (GSBM). Note that B is not usually PSD, but B > B is. The
number of non-zeros in B > B is O(md2 ). GSBMs come up in many applications, where we typically
want to solve a linear system in the normal equations of B.
Laplacian matrices are GSBMs where d = 1 and S = T = w, where w is a real number, and we
allow different w for each pair of i and j. The corresponding ISBM for Laplacians is called an edgevertex incidence matrix. Connection Laplacians are GSBMs where d = O(1) and S = T > = wQ,
for some rotation matrix Q and a real number w. Again, we allow a different rotation matrix
and scaling for every edge. For both Laplacians and Connection Laplacians, there exist linear
system solvers that run in time O(m polylog(n, −1 )) and produce  approximate solutions to the
corresponding normal equations.
We now introduce several classes of ISBMs and their associated GSBMs. Our Theorem 3.2
shows that fast linear system solvers for any of these classes would imply fast linear system solvers
for all matrices with polynomially bounded entries and condition number.
Definition 1.2 (2-commodity incidence matrix). A 2-commodity incidence matrix
 is an ISBM
where d = 2 and r = 1, and S = T , and we allow three types of S : S = w 1 0 , S = w 0 1
2
A seeming counterexample to this is Conjugate Gradient which is an iterative method that produces exact
solutions in the RealRAM model. But it requires extremely high precision calculations to exhibit this behaviour in
finite precision arithmetic, and so Conjugate Gradient is also best understood as an approximate method.
3
We use superscripts to index a sequence of vectors or matrices, and we use subscripts to denote entries of a vector
or matrix, see Section 2.

5


and S = w 1 −1 , where in each case w is a real number which may depend on the pair i and
j. We denote the set of all 2-commodity incidence matrices by MC2 . The corresponding GSBM
is called a 2-commodity Laplacian. The ISBM definition is equivalent to requiring the GSBM to
have the form






1 0
0 0
1 −1
1
2
1+2
L ⊗
+L ⊗
+L
⊗
0 0
0 1
−1 1
where ⊗ is the tensor product and L1 , L2 , and L1+2 are all Laplacian matrices.
We adopt a convention that the first variable in a set Xi is labelled u i and the second is
labelled v i . Using this convention, given a 2-commodity incidence matrix B, the equation Bx = 0
must consist of scalings of the following three types of equations: u i − u j = 0, v i − v j = 0, and
u i − v i − (u j − v j ) = 0.
Definition 1.3 (Strict 2-commodity incidence matrix). A strict 2-commodity incidence matrix is
a 2-commodity incidence matrix where the corresponding 2-commodity Laplacian has the property
that L1 , L2 , and L1+2 all have the same non-zero pattern. We denote the set of all strict 2commodity incidence matrices by MC2>0 . We denote the set of all strict 2-commodity incidence
>0
matrices with integer entries by MC2,Z
.
Linear systems in MC2>0 are exactly the systems that one has to solve to when solving
2-commodity problems using Interior Point Methods (IPMs). For readers unfamiliar with 2commodity problems or IPMs, we provide a brief explanation of why this is the case in Section 9.
>0
The MC2>0 is more restrictive than MC2 , and MC2,Z
in turn is even more restrictive. One could
>0
>0
, even if they do not exist for MC2 .
hope that fast linear system solvers exist for MC2 or MC2,Z
>0
with polyHowever, our reductions show that even getting a fast approximate solver for MC2,Z
nomially bounded entries and condition number will lead to a fast solver for all matrices with
polynomially bounded entries and condition number.
The next class we consider is 2D Truss Stiffness Matrices. They have been studied extensively
in the numerical linear algebra community [ST08, BCPT05]. For Planar 2D Trusses with some
bounds on ranges of edges, Daitch and Spielman obtained linear system solvers that run in time
e 5/4 log(1/)).
O(n
Definition 1.4 (2D Truss Incidence Matrices). Let G = (V, E) be a graph whose vertices are n
points in 2-dimension: s 1 , . . . , s n ∈ R2 . Consider X1 , . . . , Xn where d = 2. A 2D Truss Incidence
Matrix is an ISBM where d = 2 and r = 1, and for each i and j, we have S = T and S = w(s i −s j )> ,
and w is a real number that may depend on the pair i and j, but s i depends only on i and vice
versa for s j . We denote the class of all 2D Truss Incidence Matrices by T2 .
Another important class of matrices is Total Variation Matrices (TV matrices). TV matrices
come from Interior Point Methods for solving total variation minimization problem in image, see
for example [GY05] and [CMMP13b]. Not all TV matrices are GSBMs, but many GSBMs can be
expressed as TV matrices.
Definition 1.5 (TV matrix and 2-TV Incidence Matrices). Let E1 ∪ . . . ∪ Es be a partition of the
edge set of a graph. For each 1 ≤ i ≤ s, let B i be the edge-vertex incidence matrix of Ei , W i be a
diagonal matrix of edge weights, and r i be a vector satisfying W i < r i (r i )> . Given these objects,
the associated total variation matrix (TV matrix) is a matrix M defined as


X
M =
(B i )> W i − r i (r i )> B i .
1≤i≤s

A 2-TV Incidence Matrix is defined as any ISBM whose corresponding GSBM is a TV matrix with
W i ∈ R2×2 and r i ∈ R2 . We denote the class of all 2-TV incidence matrices by V2 .
6

1.4

Our Reduction: Discussion and an Example

In this section we give a brief sketch of the ideas behind our reduction from general linear systems,
over matrices in G, to multi-commodity linear systems, over matrices in MC2 , and we demonstrate
the most important transformation through an example.
The starting point for our approach is the folklore idea that any linear system can be written
as a factor-width 3 system by introducing a small number of extra variables. Using a set of multicommodity constraints, we are able to express one particular factor-width 3 equation, namely
2x00 = x + x0 . After a sequence of preprocessing steps, we are then able to efficiently express
arbitrary linear systems over integer matrices using constraints of this form. A number of further
issues arise when the initial matrix does not have full column rank, requiring careful weighting of
the constraints we introduce.
Given a matrix A with polynomially bounded integer entries and condition number, we reduce
the linear system Ax = c to a linear system By = d , where B is a strict multi-commodity edge>0
vertex incidence matrix with integer entries (i.e. in MC2,Z
), with polynomially bounded entries and
>
condition number. More precisely, we reduce A Ax = A> c to B > By = B > d . These systems
always have a solution. We show that we can find an -approximate solution to the linear system
A> Ax = A> c by a simple mapping on any y that 0 -approximately solves the linear system
B > By = B > d , where 0 is only polynomially smaller than . If A has s non-zero entries and the
maximum absolute value of an entry in A is U , then B will have O(s log(sU )) non-zero entries and
our algorithm computes the reduction in time O(s log(sU )). Note that B > B has r = O(s log(sU ))
non-zeros, because every row of B has O(1) entries. All together, this means that getting a solver
e a log(1/)) will give a solver for A with O(s
e a log(1/))
for B > Bx = B > d with running time O(r
running time.
We achieve this through a chain of reductions. Each reduction produces a new matrix and
vector, as well as a new error parameter giving the accuracy required in the new system to achieve
the accuracy desired in the original system.
1. We get a new linear system AZ,2 x Z,2 = c Z,2 where AZ,2 has integer entries, and the entries of
each row of AZ,2 sum to zero, i.e. AZ,2 1 = 0, and finally in every row the sum of the positive
coefficients is a power of two.
2. AZ,2 x Z,2 = c Z,2 is then transformed to By = d , where B is a 2-commodity edge-vertex
incidence matrix.
3. By = d is then transformed to B >0 y = d >0 , where B >0 is a strict 2-commodity edge-vertex
incidence matrix.
4. B >0 y = d >0 is then transformed to B >0,Z y = d >0,Z , where B >0,Z is a 2-commodity edgevertex incidence matrix with integer entries.
We will demonstrate step 2, the main transformation, by example. When the context is clear, we
drop the superscripts of matrices for simplicity. The reduction handles each row (i.e. equation) of
the linear system independently, so we focus on the reduction for a single row.
Consider a linear system Ax = c, and let us pick a single row (i.e. equation) Ai x = c i 4 . We
will repeatedly pick pairs of existing variables of x , say x and x0 , based on their current coefficients
in Ai x = c i , and modify the row by adding C(2x00 − (x + x0 )) to the left hand side where x00 is a
new variable and C is a real number we pick. As we will see in a moment, we can use this pairand-replace operation to simplify the row until it eventually becomes a 2-commodity equation. At
4

We use Ai to denote the ith row of A, and c i to denote the ith entry of c, see Section 2.

7

the same time as we modify Ai , we also store an auxiliary equation C(x + x0 − 2x00 ) = 0. Suppose
initially that Ai x = c i is satisfied. After this modification of Ai x = c i , if the auxiliary equation
is satisfied, Ai x = c i is still satisfied by the same values of x and x0 . Crucially, we can express the
auxiliary equation C(x + x0 − 2x00 ) = 0 by a set of ten 2-commodity equations, i.e. a “2-commodity
gadget” for this equation. Our final output matrix will not contain the equation C(x+x0 −2x00 ) = 0
as a row, but will instead contain 10 rows of 2-commodity equations from our gadget construction.
Eventually, our pair-and-replace scheme will also transform the row Ai x = c i into a 2-commodity
equation on just two variables.
Next, we need to understand how the pair-and-replace scheme makes progress. The pairing
handles the positive and the negative coefficients of Ai separately, and eventually ensures that
Ai x = c i has only a single positive and a single negative coefficient in the modified row Ai x = c i ,
in particular it is of the form ax − ax0 = c i for two variables x and x0 that appear in the modified
vector of variables x , i.e. it is a 2-commodity equation.
To understand the pairing scheme, it is helpful to think about the entries of A written using
binary (ignoring the sign of the entry). The pairing scheme proceeds in a sequence of rounds: In
the first round we pair variables whose 1st (smallest) bit is 1. There must be an even number of
variables with smallest bit 1, as the sum of the positive (and respectively negative) coefficients is
a power of 2. We then replace the terms corresponding to the 1st bit of the pair with a new single
variable with a coefficient of 2. After the first round, every coefficient has zero in the 1st bit. In
the next round, we pair variables whose 2nd bit is 1, and replace the terms corresponding to the
the 2nd bit of the pair with a new single variable with a coefficient of 4, and so on. Because the
positive coefficients sum to a power of two, we are able to guarantee that pairing is always possible.
It is not too hard to show that we do not create a large number of new variables or equations using
this scheme.
For example, let us consider an equation
3x 1 + 5x 2 + 4x 3 + 4x 4 − 16x 5 = 1
↓
Replace x 1 + x 2 by 2x 6 .
Add auxiliary equation x 1 + x 2 − 2x 6 = 0.
2(x 1 + 2x 2 + x 6 + 2x 3 + 2x 4 ) − 16x 5 = 1
↓
Replace 2(x 1 + x 6 ) by 4x 7 .
Add auxiliary equation 2(x 1 + x 6 − 2x 7 ) = 0.
4(x 2 + x 7 + x 3 + x 4 ) − 16x 5 = 1
↓
Replace 4(x 2 + x 7 ) by 8x 8 ,
and 4(x 3 + x 4 ) by 8x 9 .
Add auxiliary equations 4(x 2 + x 7 − 2x 8 ) = 0,
and 4(x 3 + x 4 − 2x 9 ) = 0.
8(x 8 + x 9 ) − 16x 5 = 1
↓
Replace 8(x 8 + x 9 ) by 16x 10 .
Add auxiliary equation 8(x 8 + x 9 − 2x 10 ) = 0.
16x 10 − 16x 5 = 1
In this way, we process Ax = c to produce a new set of equations By = d where B is a
8

2-commodity matrix. If Ax = c has an exact solution, this solution can be obtained directly from
an exact solution to By = d . We also show that an approximate solution to By = d leads to
an approximate solution for Ax = c, and we show that B does not have much larger entries or
condition number than A.
The situation is more difficult when Ax = c does not have a solution and we want
to obtain an approximate minimizer arg minx ∈Rn kAx − ck22 from an approximate solution to
arg miny ∈Rn0 kBy − d k22 . This corresponds to approximately applying the Moore-Penrose pseudoinverse of A to c. We deal with the issues that arise here using a carefully chosen scaling of each
auxiliary constraint to ensure a strong relationship between different solutions.
In order to switch from a linear system in a general 2-commodity matrix to a linear system
in a strict 2-commodity matrix, we need to reason very carefully about the changes to the null
space that this transformation inherently produces. By choosing sufficiently small weights, we are
nonetheless able to establish a strong relationship between the normal equation solutions despite
the change to the null space.

2

Preliminaries

We use subscripts to denote entries of a matrix or a vector: let Ai denote the ith row of matrix
A and Aij denote the (i, j)th entry of A; let x i denote the ith entry of vector x and x i:j (i < j)
denote the vector of entries x i , x i+1 , . . . , x j . We use superscripts to index a sequence of matrices
or vectors, e.g., A1 , A2 , . . . , and x 1 , x 2 , . . ., except when some other meaning is clearly stated.
We use A† to denote the Moore-Penrose pseudo-inverse of a matrix A. We use im(A) to denote
the image of a matrix A. We use k·k2 to denote the Euclidean norm on vectors and the spectral
norm on matrices. When M is an n × n positive semidefinite matrix, we define a norm on vectors
def √
x ∈ Rn by kx kM = x > M x . We let nnz(A) denote
Pthe number of non-zero entries
P in a matrix
A. We define kAkmax = maxi,j |Aij |, kAk1 = maxj i |Aij | and kAk∞ = maxi j |Aij |. We let
min+ (A) = mini,j s.t. Aij 6=0 |Aij |. Given a matrix A ∈ Rm×n and a vector c ∈ Rm for some m, n,
def

we call the tuple (A, c) a linear system. Given matrix A ∈ Rm×n , let Π A = A(AA> )† A> , i.e.
2
the orthogonal projection onto im(A). Note that Π A = Π >
A and Π A = Π A .

2.1

Approximately Solving A Linear System

In this section we formally define the notions of approximate solutions to linear systems that we
work with throughout this paper.
Definition 2.1 (Linear System Approximation Problem, lsa). Given linear system (A, c), where
A ∈ Rm×n , and c ∈ Rm , and given a scalar 0 ≤  ≤ 1, we refer to the lsa problem for the triple
(A, c, ) as the problem of finding x ∈ Rn s.t.
kAx − Π A ck2 ≤  kΠ A ck2 ,
and we say that such an x is a solution to the lsa instance (A, c, ).
This definition of a lsa instance and solution has several advantages: when im(A) = Rm , we get
Π A = I , and it reduces to the natural condition kAx − ck2 ≤  kck2 , which because im(A) = Rm ,
can be satisfied for any , and for  = 0 tells us that Ax = c.
When im(A) does not include all of Rm , the vector Π A c is exactly the projection of c onto
im(A), and so a solution can still be obtained for any . Further, as (I − Π A )c is orthogonal to
9

Π A c and Ax , it follows that
kAx − ck22 = k(I − Π A )ck22 + kAx − Π A ck22 .
Thus, when x is a solution to the lsa instance (A, c, ), then x also gives an 2 kΠ A ck22 additive
approximation to
(1)
minn kAx − ck22 = k(I − Π A )ck22 .
x ∈R

Similarly, an x which gives an additive 2 kΠ A ck22 approximation to Problem (1) is always a
solution to the lsa instance (A, c, ). These observations prove the following (well-known) fact:
Fact 2.2. Let x ∗ ∈ arg minx ∈Rm kAx − ck22 , then for every x ,
kAx − ck22 ≤ kAx ∗ − ck22 + 2 kΠ A ck22
if and only if x is a solution to the lsa instance (A, c, ).
When the linear system Ax = c does not have a solution, a natural notion of solution is any
minimizer of Problem (1). A simple calculation shows that this is equivalent to requiring that x
is a solution to the linear system A> Ax = A> c, which always has a solution even when Ax = c
does not. The system A> Ax = A> c is referred to as the normal equation associated with Ax = c
(see [TBI97]).
Fact 2.3. x ∗ ∈ arg minx ∈Rn kAx − ck22 , if and only if A> Ax ∗ = A> c, and this linear system
always has a solution.
This leads to a natural question: Suppose we want to approximately solve the linear system
A> Ax = A> c. Can we choose our notion of approximation to be equivalent to that of a solution
to the lsa instance (A, c, )?
A second natural question is whether we can choose a notion of distance between a proposed
solution x and an optimal solution x ∗ ∈ arg minx ∈Rn kAx − ck22 s.t. this distance being small is
equivalent to x being a solution to the lsa instance (A, c, )? The answer to both questions is yes,
as demonstrated by the following facts:
Fact 2.4. Suppose x ∗ ∈ arg minx ∈Rn kAx − ck22 then
1.

A> Ax − A> c

(A> A)†

= kAx − Π A ck2 = kx − x ∗ kA> A .

2. The following statements are each equivalent to x being a solution to the lsa instance
(A, c, ):
A> Ax − A> c (A> A)† ≤  A> c (A> A)† if and only if x is a solution to the lsa
instance (A, c, ).
(b) kx − x ∗ kA> A ≤  kx ∗ kA> A if and only if x is a solution to the lsa instance (A, c, ).

(a)

For completeness, we prove Fact 2.4 in Appendix A. Fact 2.4 explains connection between our
Definition 2.1, and the usual convention for measuring error in the Laplacian solver literature [ST14].
In this setting, we consider a Laplacian matrix L, which can be written as L = A> A ∈ Rn×n , and
a vector b s.t. ΠA> A b = b. This condition on b is easy to verify in the case of Laplacians, since for
the Laplacian of a connected graph, ΠA> A = I − n1 11> . Additionally, it is also equivalent to the
condition that there exists c s.t. b = A> c. For Laplacians it is possible to compute both A and
a vector c s.t. b = A> c in time linear in nnz(L). For Laplacian solvers, the approximation error
of an approximate solution x is measured by the  s.t. A> Ax − b (A> A)† ≤  kbk(A> A)† . By
Fact 2.4, we see that this is exactly equivalent to x being a solution to the lsa instance (A, c, ).
10

2.2

Measuring the Difficulty of Solving a Linear System

Running times for iterative linear system solvers generally depend on the number of non-zeros in
the input matrix, the condition number of the input matrix, the accuracy, and the bit complexity.
In this section, we formally define several measures of complexity of the linear systems we use.
This is crucial, because we want to make sure that our reductions do not rely on mapping into
extremely ill-conditioned matrices, and so we use these measures to show that this is in fact not
the case.
Definition 2.5.
1. Given a matrix A ∈ Rm×n , q
we define the maximum singular value σmax (A) in the usual way
as σmax (A) = maxx ∈Rn ,x 6=0

x > A> Ax
.
x >x

2. Given a matrix A ∈ Rm×n which is not all zeros,
q we define the minimum non-zero singular
> >
value σmin (A) as σmin (A) = minx ∈Rn ,x ⊥ null(A) x xA> xAx .
3. Given a matrix A ∈ Rm×n which is not all zeros, we define the non-zero condition number of
(A)
A as κ(A) = σσmax
.
min (A)
Definition 2.6. The sparse parameter complexity of an lsa instance (A, c, ) where A ∈ Zm×n
and nnz(A) ≥ max(m, n), and  > 0, is




1
1
def
−1
S(A, c, ) = nnz(A), max kAkmax , kckmax ,
,
, κ(A), 
.
min+ (A) min+ (c)
Note in the definition above that when A 6= 0 and c 6= 0 have only integer entries, we trivially
have min+ (A) ≥ 1 and min+ (c) ≥ 1. However, including min+1 (A) , and min1+ (c) in the definition
stated above is useful when working with intermediate matrices whose entries are not integer valued.

2.3

Matrix Classes and Reductions Between Them

We use the term matrix class to refer to an infinite set of matrices M. In this section, we formally
define a notion of efficient reduction between linear systems in different classes of matrices.
Definition 2.7 (Efficient f -reducibility). Suppose we have two matrix classes M1 and M2 , and
there exist two algorithms A1→2 and A1←2 s.t. given an lsa instance (M 1 , c 1 , ), where M 1 ∈ M1 ,
the call A1→2 (M 1 , c 1 , 1 ) returns an lsa instance (M 2 , c 2 , 2 ) s.t. if x 2 is a solution to lsa instance
(M 2 , c 2 , 2 ) then x 1 = A1←2 (M 1 , M 2 , x 2 ) is a solution to lsa instance (M 1 , c 1 , 1 ).
Consider a function of f : R4+ → R4+ s.t. every output coordinate is an increasing function of
every input coordinate. Suppose that we always have
S(M 2 , c 2 , 2 ) ≤ f (S(M 1 , c 1 , 1 )),
and the running times of A1→2 (M 1 , c 1 , 1 ) and A1←2 (M 1 , M 2 , x 2 ) are both bounded by
O(nnz(M 1 )).
Then we say that M1 is efficiently f -reducible to M2 , which we also write as
M1 ≤f M2 .
Lemma 2.8. Suppose M1 ≤f M2 and M2 ≤g M3 . Then M1 ≤g◦f M3 .
11

Proof. The proof is simply by the trivial composition of the two reductions.
Definition 2.9. We let G denote the class of all matrices with integer valued entries s.t. there is
at least one non-zero entry in every row and column5 .

3

Main Results

In this section, we use the notions of sparse parameter complexity and matrix class reductions
to prove our main technical result, Theorem 3.1, which shows that linear systems in general matrices with integer entries can be efficiently reduced to linear systems in several different classes
of Incidence Structured Block Matrices. From this result, we derive as corollary our main result,
Theorem 3.2, which states that fast high accuracy solvers for several types of ISBMs imply fast
high accuracy solvers for all linear systems in general matrices with integer entries.
Theorem 3.1. Let f (s, U, K, ) = (O(s log(sU )), poly(U K−1 s), poly(U K−1 s), poly(U K−1 s)),
then
>0
1. G ≤f MC2,Z
.

2. G ≤f T2 .
3. G ≤f V2 .
Theorem 3.2. Suppose we have an algorithm which solves every Linear System Approximation Problem (A, c, ) with sparse parameter complexity S(A, c, ) ≤ (s, U, K, −1 ) in time
O(sa polylog(s,oU, K, −1 )) for some a ≥ 1, whenever A ∈ R for at least one of R ∈
n
>0
>0
, T2 , or V2 .
, T2 , V2 . I.e. we have a “fast” solver6 for one of the matrix classes MC2,Z
MC2,Z
Then every Linear System Approximation Problem (A, c, ) where A ∈ G with sparse parameter
complexity S(A, c, ) ≤ (s, U, K, −1 ) can be solved in time O(sa polylog(s, U, K, −1 )).
Proof. The theorem is a immediate corollary of Theorem 3.1.
Definition 3.3. We let Gz,2 denote the class of all matrices with integer valued entries s.t. there
is at least one non-zero entry in every row and column, and every row has zero row sum, and for
each row, the sum of the positive coefficients is a power of 2.



Lemma 3.4. Let f (s, U, K, ) = O(s), O −1 s9/2 U 3 , O −1 s8 U 3 K , O s5/2 U 2 −1 , then
G ≤f Gz,2 .





Lemma 3.5. Let f (s, U, K, ) = O(s log(sU )), O s3/2 U log1/2 (sU ) , O Ks4 U 2 log2 (sU ) , O sU 2 −1 ,
then
Gz,2 ≤f MC2 .



Lemma 3.6. Let f (s, U, K, ) = O(s), O −1 U 2 K , O −1 s2 U 2 K , O(−1 ) , then
MC2 ≤f MC2>0 .
5
If there is a row or column with only zeros, then it can always be handled trivially in the context of solving linear
systems
6
The reduction requires only a single linear system solve, and uses the solution in a black-box way. So the reduction
also applies if the solver for the class R only works with high probability or only has running time guarantees in
expectation.

12

Lemma 3.7. Let f (s, U, K, ) = (s, −1 sU, 2K, O(−1 )), then
>0
MC2>0 ≤f MC2,Z
.

Lemma 3.8. Let f (s, U, K, ) be as defined in Lemma 3.5 then
Gz,2 ≤f T2 .
Lemma 3.9. Let f (s, U, K, ) = (s, U, K, −1 ), then
MC2 ≤f V2 .
Proof of Theorem 3.1. Follows by appropriate composition (Lemma 2.8) applied to the the Lemmas
above, i.e. 3.4, 3.5, 3.6, 3.7, 3.8 and 3.9.

3.1

Outline of Remaining Sections

In Section 4 is presents the proof of Lemma 3.5, i.e. Gz,2 ≤f MC2 . This statement is our most
important reduction. In Section 5, we prove Lemma 3.5, the surprising reduction MC2 ≤h MC2>0 ,
which shows that we can solve normal equations even while changing null-spaces involved substantially. In Section 5, we show how to round weights in a 2-commodity problem to integers, proving
Lemma 3.7. In Section 7 is presents the proof of Lemma 3.4. This is a simpler reduction that we
use to establish the properties that our Lemma 3.5 relies on. In Section 8, we present the proof
of Lemma 3.8. Section 9 describes how Interior Point Methods give rise to multi-commodity and
Total Variation linear systems. Section 9 also contains a proof of Lemma 3.9.

4

Reducing Zero-Sum Power Two Linear Systems to TwoCommodity Linear Systems

To prove Lemma 3.5, we need to provide mapping algorithms AGz,2 →MC2 for mapping linear system
approximation (lsa) instances over matrices in Gz,2 to lsa problem instances over matrices in MC2 ,
as well as AMC2 ←Gz,2 for mapping the resulting solutions back. These leads to the following main
components:
• Algorithm 1 states the pseudo-code for the algorithm Reduce Gz,2 toMC2 , which implements
the desired mapping of problem instances. Given lsa problem instance (A, c A , A ) where
A ∈ Gz,2 , the call Reduce Gz,2 toMC2 (A, c A , A ) returns an lsa problem instance (B, c B , B )
where B ∈ MC2 . (Strictly speaking, Reduce Gz,2 toMC2 also has a parameter α which we
will set before using the algorithm.)
• Algorithm 2 provides the pseudo-code for MC2 Gadget, a short subroutine used in
Reduce Gz,2 toMC2 to represent the equation a + b = 2c using two-commodity constraints.
• Algorithm 3 provides the (trivial) pseudo-code for MapSoln MC2 toGz,2 used to map solutions to lsa problems over MC2 back to solutions over Gz,2 by restricting onto the original
variables.

13

Pseducode of the key reduction routine that creates the new linear system, Reduce Gz,2 toMC2 ,
is shown
a single
 > in Algorithm1. >Note that given two singleton multi-sets eachcontaining
 > equation,
>x = c
e.g.
a
x
=
c
and
a
x
=
c
where
a
,
a
are
vectors,
we
define
a
+
a 2 x = c2 =
1
1
2
1
2
2
1
 > 1 >

 >
 >
>x = c .
a 1 x + a 2 x = c1 + c2 and we define a >
x
=
c
∪
a
x
=
c
=
a
x
=
c
,
a
1 2
2
1
2
1
2
1
Algorithm 1 Reduce Gz,2 toMC2
Input: (A, c A , 1 ) where A ∈ Gz,2 is an m × n matrix, c A ∈ Rn , 0 < 1 < 1, and α > 0.
Output: (B, c B , 2 ) where B ∈ MC2 is an m × n matrix, c B ∈ Rn , and 0 < 2 < 1.

−1/2
2 2

(A)
kc A k2 σmax
1 −1/2
1: 2 ← 1 1 + α
1+
α+1
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:

X ← {u 1 , . . . , u n , v 1 , . . . , v n }
{MC2 variables and index of new variables}
Let x be the vector of variables corresponding to the set of variables X
t←n+1
Initialize A ← ∅, B ← ∅
{Multisets of main and MC2 auxiliary equations respectively}
for each equation 1 ≤ i ≤ m in A do
Let I +1 ← {j : Aij > 0}, I −1 ← {j : Aij < 0}
if I +1 = 1 and I −1 = 1 then
Let the only elements in I +1 and I −1 be j+ and j− respectively. {Note Aij+ = −Aij− }
A ← A ∪ {Aij+ u t − Aij+ u j− = c i }
wi ← α
1/2
B ← B ∪ wi · {Aij+ u j+ − Aij+ u t = 0}
X ← X ∪ {u t , v t }, update x accordingly
t←t+1
else
Ai ← {Ai u = c i }
Bi ← ∅
for s = −1, +1 do
r←0
while Ai has strictly more than 1 coefficient with sign s do
b ij be the coefficient of u j in Ai .
For each j, let A
s
s
b ij |/2r c is odd}
Iodd ← {j ∈ I : b|A
s
into k disjoint pairs (jk , lk )
Pair the indices of Iodd
for each pair of indices (jk , lk ) do
Ai ← Ai + s · 2r {(2u t − (u jk + u lk )) = 0}
Bi ← Bi ∪ −s · 2r · MC2 Gadget(u jk , u lk , t)
X ← X ∪ {u t , . . . , u t+6 , v t , . . . , v t+6 }, update x accordingly
t←t+7
r ←r+1
end for
end while
end for
wi ← α |Bi |
1/2
B ← B ∪ wi · B i
A ← A ∪ Ai .
end if
end for
return 2 and B, c s.t. Bx = c corresponds to the equations in A ∪ B, on the variable set X .
14

The central object created by Algorithm 1 is the matrix B, which contains both new equations
and new variables. We will superscript the variables with A to distinguish variables appear in th
original equation Ax A = c A from new variables. We will term the new variables as x aux , and write
a vector solution for the new problem, x B , as:
 A 
x
B
x =
.
(2)
x aux
Let nA be the dimension of x A , and nB be the dimension of x B , respectively. We order the
variables so that for an appropriately chosen index h we have
1. x aux
1:h corresponds to the u-coordinates of the auxiliary variables created in MC2 -gadgets.
2. x aux
h+1:nB −nA corresponds to the v -coordinates of the auxiliary variables created in MC2 gadgets.
With this ordering x B
t = u t for t ≤ nA + h.
Furthermore, we will distinguish the equations in B into ones formed from manipulating A,
i.e. the equations added to the set A, from the auxiliary equations, i.e. the equations added to the
1/2
set B. We use W 1/2 = diagi (wi ) to refer to the diagonal matrix of weights wi applied to the
auxiliary equations B in Algorithm 1. In Algorithm 1, a real value α > 0 is set initially and used
1/2
when computing the weights wi . For convenience, thoughout most of this section, we will treat
α as an arbitrary constant, and only eventually substitute in its value to complete our main proof.
This leads to the following representation of B and c B which we will use throughout our analysis
of the algorithm:
!
b
A
B=
.
(3)
b
W 1/2 B
b corresponds to A in Reduce Gz,2 toMC2 , and B
b corresponds to the
Here the equations of A
auxiliary constraints, i.e. equations of B in Reduce Gz,2 toMC2 . Also, the vector c B created is
simply an extension of c A :
 A 
c
B
.
(4)
c =
0
Finally, as Algorithm 1 creates new equations for each row of A independently, we will use Si
b that’s created from Ai,: aka. the the auxiliary
to denote the subset of indices of the rows of B
constraints generated from the call to MC2 Gadget upon processing Ai,: . We will also denote the
size of these using
def
mi = |Si | ,
b i to denote these part of B
b that corresponds to these rows. The Gadget routine used in
and use B
the reduction, and the (trivial) solution mapper are stated below.

15

Algorithm 2 MC2 Gadget
Input: Scalar variables u jk , u lk , integer t.
Output: Set of equations for a multi-commodity gadget that computes the average of u jk and
u lk .
1: return
{u t+3 − v t+3 − (u t+4 − v t+4 ) = 0,
u t − u t+3 = 0,
u t+4 − u jk = 0,
v t+3 − v t+1 = 0,
v t+2 − v t+4 = 0,
u t+5 − v t+5 − (u t+6 − v t+6 ) = 0,
u t − u t+5 = 0,
u t+6 − u lk = 0,
v t+5 − v t+2 = 0,
v t+1 − v t+6 = 0}

Algorithm 3 MapSoln MC2 toGz,2
0

Input: m × n matrix A ∈ Gz,2 , m0 × n0 matrix B ∈ MC2 , vector c A ∈ Rm , vector x B ∈ Rn .
Output: Vector x A ∈ Rn .
1: if A> c A = 0 then
2:
return x A ← 0
3: else
4:
return x A ← x B
1:n
5: end if
We upper bound the number of nonzero entries of B by the number of nonzero entries of A.
Lemma 4.1.
• We have nnz(B) = O (nnz(A) log kAk∞ ).
• Both dimensions of B are O (nnz(A) log kAk∞ ).
• The number of iterations required by Algorithm 1 to process row Ai is at most log kAi k1 .
• When processing row Ai , every new variable that appears in the Ai equation is paired in the
following iteration of the while-loop in Line 20 (unless it is the final iteration), and then
disappears from Ai .
Proof. Since Algorithm 1 constructs new equations for each row of A independently, we!bound the
bi
A
number of new variables and new equations (that is, the size of the submatrix
b i ) for each
B
row i of A separately.
Let ni be the number of nonzero entries of the ith row of A. We count the number of variables
in the equation Ai , at each iteration of the while-loop in line 20 of Algorithm 1. Let X (r) be the
16

(r)

subset of variables with nonzero coefficients in Ai , at the end of iteration r. Let Xaux be the subset
of X (r) containing all auxiliary variables created in iteration r.
Note in each iteration, Algorithm 1 replaces two variables by a new auxiliary variable. It gives
that
ni
(1)
(1)
Xaux
≤ , and X (1) \ Xaux
≤ ni .
2
(1)

Since each auxiliary variable in Xaux has coefficient 2, in the 2nd iteration, together with another
variable of coefficient 2, it must be replaced by a new auxiliary variable. Thus, at the end of the
(1)
2nd iteration, all auxiliary variables Xaux will not appear in the equation Ai . This implies that
(2)
Xaux
≤

ni ni
(2)
+ , and X (2) \ Xaux
≤ ni .
2
4

Similarly, at the end of the tth iteration, we have
X ni
(r)
Xaux
≤
≤ ni ,
2s
1≤s≤r

and
(r)
X (r) \ Xaux
≤ ni .

Since Algorithm 1 pulls out a factor 2 in each iteration, the total number of iterations is at
most log kAi k1 . Since each auxiliary variable in Ai during the construction corresponds to O(1)
b i , the total number of auxiliary variables and equations
auxiliary variables and O(1) equations in B
for row i of A is
O (ni log kAi k1 ) .
Therefore, the number of variables and the number of equations in B (that is, the both dimensions
of B) are
O (nnz(A) log kAk∞ ) .
Since each row of B has O(1) nonzero coefficients, we have
nnz(B) = O (nnz(A) log kAk∞ ) .
This completes the proof.

4.1

Reduction Between Exact Solvers

The most important relation between A and B is given by the following claim.
Claim 4.2 (Reduction between exact solvers). Fix any x A ∈ Rn . Then
 A 
2
α+1
x
B
A
A 2
Ax − c 2 =
min
B
−
c
.
x aux
α x aux
2
As a Corollary of Claim 4.2, we observe the following:
Lemma 4.3. Given lsa problem instance (A, c A , A ) where A ∈ Gz,2 , suppose
(B, c B , B ) = Reduce Gz,2 toMC2 (A, c A , A ).
 A 
x
B
Then B ∈ MC2 and if x =
is a solution to the exact lsa problem (B, c B , 0), then x A
x aux
is a solution to the exact lsa problem (A, c A , 0).
17

Proof. This follows immediately from minimizing over x A on both sides of the equation established
by Claim 4.2 and then applying and Fact 2.2.
Before proving Claim 4.2 we first note a basic guarantee obtained by Reduce Gz,2 toMC2 :
!
 A !
b
A
ci
i
1>
xB −
= Ai x A − c A
(5)
i .
b
0
Bi
To verify this guarantee, we consider two cases separately. The first case is when the condition
I +1 = 1 AND I −1 = 1 is true (see Algorithm 1, Line 7). In this case, the main constraint in
B
A
the output corresponding to row i is {Aij+ x B
t − Aij+ x j− = c i }, while the auxiliary constraints
B
contain only a single row {Aij+ x B
j+ − Aij+ x t = 0}, and adding these proves the guarantee for this
case. The second case is when the condition in Algorithm 1, Line 7) is false. We consider the case
s = +1. The case s = −1 is proved similarly. Note that we will refer to variables jk and lk only
in the context of a fixed value of t, which always ensures
defined. In
 that they
 are unambiguosly

B
B
B
r
= 0, we we also use
the case s = +1, each time we modify Ai by adding 2 2x t+1 − x jk + x lk



B + xB
MC2 Gadget to create auxiliary constraints that sum up to exactly −2r 2x B
−
x
, so
t+1
jk
lk
adding these together will cancel out the changes.
Proof of Claim 4.2. For convenience, we also write
def b
B
1. δj = B
jx ,
def b
B
A
2. δ̂i = A
ix − ci ,
def

3. i = Ai x A − c A
i .
Thus
!

bi
A
1/2 b
Wi B
i

cA
i
0



xB −



2

= δ̂i2 + wi

X

δj2 .

j∈Si

2

Summing over all rows, we get
xA


B

x aux





2

− cB

=

X

2



δ̂i2 + wi

i

X

δj2  .

j∈Si

Similarly, the square of row i of Ax A − c A is 2i and summing over all rows we get
X
2
2i .
Ax A − c A 2 =
i

Recalling Equation (5), we have
>

1

bi
A
bi
B

!
B



x −

cA
i
0

!

Thus
δ̂i +

X

δj = i .

j∈Si

18

= Ai x A − c A
i .

By the Cauchy-Schwarz inequality (applied to two vectors a and b given by a 1 = δ̂i , a j =
√
√
wi δj , b 1 = 1, b j = 1/ wi ),
2


2i

= δ̂i +

X





δj  ≤ δ̂i2 + wi

j∈Si

X

δj2 

j∈Si



mi
1+
wi


.

the equality holds if and only if
∀j ∈ Si .

δ̂i = wi δj

(6)

Note that we have ensured that for all i, wi = αmi .
By summing over rows we conclude that for every x A and every x aux , we have
 A 
α+1
x
A
A 2
Ax − c 2 ≤
B
− cB .
x aux
α
2

(7)

The inequality above will be an equality if Equations (6) are satisfied. We now show that for every
fixed x A , minimizing over x aux ensures that (4.1) holds with equality.
In particular, we will momentarily prove the following Claim.
Claim 4.4. For any fixed x A and its associated i values, for each row i of A, the linear system
 A 
α
x
b
i ,
(8)
= cA
Ai
i +
aux
x
α+1
 A 
1
x
b
Bj
=
i , ∀j ∈ Si .
(9)
x aux
(α + 1)mi
has a solution (which may not be unique).
Since every auxiliary variable is associated with only one row i of A, Claim 4.4 implies that we
can choose x aux s.t. all these linear systems are satisfied at once.
Given such a choice of x aux , we get that Equations (6) are satisfied so
Ax A − c A

2
2

=

α+1
B
α



xA
x aux



2

− cB

.
2

Given Claim 4.4, this completes the proof of Claim 4.2.
Proof of Claim 4.4. We focus on the case of the condition I +1 = 1 AND I −1 = 1 in Algorithm 1, Line 7 being false. The case when then condition is true is very similar, but easier as it
deals with a set of just two equations.
We will construct an assignment to all the variables of x aux s.t. Equations (8) and (9) are
satisfied. We start with an assignment x A to the main variables, and we then assign values to
auxiliary variables in the order they are created by the algorithm Reduce Gz,2 toMC2 . Note that
we will refer to variables jk and lk only in the context of a fixed value of t, which always ensures
B
that they are unambiguosly defined. When the algorithm processes pair x B
jk , x lk = u jk , u lk , the
value of these variables will have been set already, while x B
t = u t and the other newly created
auxiliary variables have not. Every auxiliary variable is associated with only one row, so we never
get multiple assignments to a variable using this procedure.

19

Recall the constraints created by the MC2 Gadget call are
{u t+3 − v t+3 − (u t+4 − v t+4 ) = 0,
u t − u t+3 = 0,
u t+4 − u jk = 0,
v t+3 − v t+1 = 0,
v t+2 − v t+4 = 0,
u t+5 − v t+5 − (u t+6 − v t+6 ) = 0,
u t − u t+5 = 0,
u t+6 − u lk = 0,
v t+5 − v t+2 = 0,
v t+1 − v t+6 = 0}
Let z = (u jk , u lk , u t , u t+1 , u t+2 , u t+3 , . . . , u t+6 , v t , v t+1 , . . . , v t+6 ) ∈ R16 . Let G ∈ R10×16 be
the matrix s.t. Gz = 0 corresponds to the constraints listed above. Note that all coefficients of
u t+1 , u t+2 and v t are zero. We set these three variables to zero.
For some , which we will fix later, we choose u t such that
ut =

1
(u jk + u lk ) + 5.
2

Again, for the same , we fix the following values for u t+i , v t+j , 3 ≤ i ≤ 6, 1 ≤ j ≤ 6
u t+3 = u t − ,
u t+4 = u jk + ,
u t+5 = u t − ,
u t+6 = u lk + ,
v t+1 = 0,
v t+2 = 5 − (u t − u jk ),
v t+3 = ,
v t+4 = 4 − (u t − u jk ),
v t+5 = 6 − (u t − u jk ),
v t+6 = −.
b i x B = 2r Gz = 2r 1, so by choosing
This ensures Gz = 1. Note that for some r, B
1
−r
 = 2 (α+1)mi i , we can ensure Equations (9) are satisfied.
Also
!
bi
A
A
>
b i x B + 1> B
b ix B = A
b i x B + 1 i .
c i + i = 1
xB = A
b
α+1
Bi
Thus, we have
b ix B = cA +
A
i
which is Equation 8.
This completes the proof of the claim.

20

α
i ,
α+1

Remark. The optimal solutions for minx B Bx B − c B 2 and minx A Ax A − c A
one-to-one map, however, the optimal values are different:
From the proof of Claim 4.2 and Equation (1)


1
2
A
A 2
c − ΠA c 2 = 1 +
c B − ΠB c B 2 ,
α

2

have a

where we set the weight wi = αmi , ∀i. Note when α → ∞, the two optimal values approach the
same value.

4.2

Relationship Between Schur Complements

n×n be a symmetric PSD matrix. Write C =
Definition 4.5
 (Schur Complement). Let C ∈ R

C 11 C 12
where C 11 , C 12 , C 12 are block matrices, the Schur complement of C is
C>
12 C 22
def

>
C S = C 11 − C 12 C −1
22 C 12 .
def

If C 22 is not invertible, then we use pseudo-inverse, that is, C S = C 11 − C 12 C †22 C >
12 .
The following fact is important for Schur complement.
Fact 4.6 (Schur complement). For any fixed vector x ,



x
>
>
y
= x >C S x .
min x
C
y
y
Proof. We expand the left hand side,



x
>
>
x
y
= x > C 11 x + 2x > C 12 y + y > C 22 y .
C
y

(10)

Taking derivative w.r.t. y and setting it to be 0 give that
2C 22 y + 2C >
12 x = 0.
Plugging y = −C †22 C >
12 x into (10),



x
>
>
>
y
min x
C
= x > C 11 x − x > C 12 C †22 C >
12 x = x C S x .
y
y
This completes the proof.
Recall that B is the output coefficient matrix of
B 1 is the submatrix corresponding to x A and B 2 is
 >
B1 B1
>
B B=
B>
2 B1
Claim 4.7.

>
α
α+1 A A


Algorithm 1. Write B = B 1 B 2 , where
the submatrix corresponding to x aux . Then,

B>
1 B2
.
B>
2 B2

is the Schur complement of B > B.

Proof. By (4.2), for any fixed x A ,
min
x > B > Bx =
aux
x

By Claim 4.6,

>
α
α+1 A A

>
α
x A A> Ax A .
α+1

is the Schur complement of B > B.
21

4.3

Approximate solvers

We now show that approximate solvers for B also translate to approximate solvers for A. The
following Lemma about the length of projections involving integral matrices is crucial for our
bounds.
Lemma 4.8. Let A and c a matrix and a vector such that kAck22 is integral. If A> c 6= 0, then
kΠA ck22 ≥

1
2
σmax
(A)

=

λmax

1
.
A> A

Proof. As kAck22 is integral, the condition of A> c 6= 0 implies A> c
A:
X
A=
σi r i (s i )> .

2

≥ 1. Consider the SVD of

i
>

Substituting it in for A gives
A> c =

X


σi (r i )> c s i ,

i

(r i )> c

and we will use αi = σi
to denote the coefficients of c against thePsingularPvectors of A.
Note that σi = 0 implies αi = 0. The norm condition on A> c also means i αi2 = i:σi 6=0 αi2 ≥ 1,
which then gives

!
!

†
X
X
X
X α2
1
1
2
i i >
i
i
kΠA ck2 = c > A A> A A> c =
αi (s i )> 
.
s
(s
)
≥ 2
α
s
=
i
2
2
σmax (A)
σi
σi
i
i
i:σi 6=0

i:σi 6=0

This completes the proof.
Lemma 4.9. Let Bx B
=
c B be the linear system returned by a call to
A
A
Reduce Gz,2 toMC2 (A, c ,  ), and let B be the the error parameter returned by this call.
Let x B be a vector such that
Bx B − ΠB c B

≤ B ΠB c B

2

2

.

(11)

Let xe A be the vector returned by a call to MapSoln MC2 toGz,2 (A, B, c A , x B ). Then,
Ae
x A − ΠA c A

2

≤ A ΠA c A

2

.

Proof. We first consider the case A> c A = 0, then Algorithm 3 returns xe A = 0, which gives
Ae
x A − ΠA c A

2

= 0,

A
> A
A
and completes
 A  the proof for this case. When A c 6= 0, the Algorithm 3 returns xe = x , where
x
xB =
. Let
x aux
 A∗ 
x
B∗ def
x
=
∈ arg min Bx − c B 2 .
x aux∗
x

By Fact 2.2, the condition (11) is equivalent to
x B − x B∗

B >B

≤ B ΠB c B

22

2

.

By Claim (4.2), x A∗ ∈ arg minx Ax − c A 2 .
We now upper bound the error of x A − x A∗
B

x −x

. Note

x A − x A∗
x aux − x aux∗



B∗

A> A

=


.

By Claim 4.6 and Claim 4.7,
2
x A∗ A> A

A

x −
Then, we lower bound ΠA c A

2

A



1
x B − x B∗
≤ 1+
α

2
B >B

.

. By Claim 4.2,

c −

2
ΠA c A 2


=

1
1+
α



c B − ΠB c B

2
.
2

Since ΠA c A ⊥ c A − ΠA c A , and ΠB c B ⊥ c B − ΠB c B , we have
2
2

ΠA c A

2
2
B 2
c 2

2
2
B 2
ΠB c 2 −

= cA

− c A − ΠA c A

=

− cB −
2
2

= ΠB c B

−

1
c A − ΠA c A
α+1

1
c A − ΠA c A
α+1

2
2

2
.
2

Thus,
ΠB c B

2
2

= ΠA c A

2
2

1
c A − ΠA c A
α+1

+

2
.
2

Since A> c A 6= 0, and because A> c A is integral, by Lemma 4.8, we have
ΠA c A

2
2

1

≥

2 (A)
σmax

.

Thus,
2
ΠB c B 2

≤

2
ΠA c A 2

1
+
cA
α+1

2
2

≤

1+

cA

!
2 2
σ
(A)
2 max
α+1

2
.
2

ΠA c A

Therefore,
A

Ax − ΠA c

A
2



1 1/2
≤ 1+
α

cA

1+

!1/2
2 2
σ (A)
B
2 max


α+1

Finally, by the setting
A

B =
1+


1 1/2
α



1+

2
kc A k22 σmax
(A)
α+1

1/2 ,

we have
Ax A − ΠA c A

2

≤ A ΠA c A

This completes the proof.

23

2

.

ΠA c A

2

.

4.4

Bounding Condition Number of the New Matrix

In this section, we show that the condition number of B is upper bounded by the condition number
of A with a poly(n) multiplicative factor.
 A 
x
B
We first characterize the null space of B. Recall that in Equation (2), we write x =
.
x aux
In the following, we will employ a different indexing of x aux than the one defined at the beginning
of Section 4. We also reorder the columns of B so that Bx B = 0 represents the same equations as
before. For appropriately chosen indices g1 and g2 , we define
1. x aux
1:g1 corresponds to the u-coordinates of the auxiliary variables created in MC2 -gadgets
whose coefficients are nonzero.
2. x aux
g1 +1:g2 corresponds to the v -coordinates of the auxiliary variables created in MC2 -gadgets
whose coefficients are nonzero.
3. x aux
g2 +1:nB −nA corresponds to the coordinates of the auxiliary variables created in MC2 -gadgets
whose coefficients are zero.
Using this ordering, for 0 ≤ i ≤ (g2 − g1 )/6, x aux
g1 +6i+1:g1 +6i+6 corresponds to the v -coordinates of
the auxiliary variables with non-zero coefficients in a single MC2 -gadget.
Given any fixed x A ∈ null(A), we extend x A to a vector in dimension nB :
 A 
x
A def
.
(12)
p(x ) =
x aux
We assign the values of the auxiliary variables of p(x A ) in the order that they are created in
Algorithm 1. In a MC2 -gadget, suppose the values of variables u jk and u lk have already been
assigned. Let u t , u t+i , v t+j , 3 ≤ i ≤ 6, 1 ≤ j ≤ 6 be the auxiliary variables created in this gadget.
We assign values as follows,
u t = (u jk + u lk )/2,
u t+3 = u t .
u t+4 = u jk ,
u t+5 = u t ,
u t+6 = u lk ,
v t+1 = 0,
v t+2 = u t − u jk ,
v t+3 = 0,
v t+4 = u t − u jk ,
v t+5 = u t − u jk ,
v t+6 = 0.
This gives the first g2 entries of x aux , and we set all the rest of the entries to be 0.
Let e i ∈ R(g2 −g1 )/6 be the ith standard basis vector and 1 be the all-one vector in 6 dimensions.
We define p i ∈ RnB to be a vector whose nonzero entries are given by
p inA +g1 +1:nA +g2 = e i ⊗ 1.
24

(13)

Let e j ∈ RnB −nA −g2 is the jth standard basis vector. We define q j ∈ RnB to be a vector whose
nonzero entries are given by
q jnA +g2 +1:nB = e j .
Lemma 4.10. null(B) = Span p(x A ), p i , q j : x A

(14)

1
∈ null(A), 1 ≤ i ≤ g2 −g
6 , 1 ≤ j ≤ nB − g2 .


1
Proof. Let S be the subspace of Span p(x A ), p i , q j : x A ∈ null(A), 1 ≤ i ≤ g2 −g
6 , 1 ≤ j ≤ nB − g2 .
According to definitions of the vectors, we can check that S ⊆ null(B).
It remains to show that S ⊇ null(B). Let x ∈ null(B). By Claim 4.2 with c A = 0, we have
Ax A = 0,
that is, x A ∈ null(A). According to the MC2 -gadget constraints in Algorithm 2, we have
v t+1 = v t+3 = v t+6 = γ, and v t+2 = v t+4 = v t+5 = θ,
where γ − θ = (u jk − u jl )/2. Besides, since all entries in x nA +g2 +1:nB have zero coefficients, they
are free to choose. Thus, x ∈ S , that is, S ⊇ null(B).
This completes the proof.
We upper bound the largest singular value of B in the following claim.



Lemma 4.11. λmax B > B = O α nnz(A)2 kAk2∞ log2 kAk∞ .
Proof. By the Courant-Fischer Theorem,


x > B > Bx
λmax B > B = max
.
x
x >x
We write x as (u 1 , v 2 , . . . , u nB , v nB ), so that each (u i , v i ) corresponds to the two coordinates
of vertex i in the graph. Expanding the right hand side, we get
P
P
P
2
2
2
(i,j)∈E2 w(i,j) (v i − v j ) +
(i,j)∈E1+2 w(i,j) (u i − v i − (u j − v j ))
(i,j)∈E1 w(i,j) (u i − u j ) +
P 2
2
i ui + vi
where w(i,j) ’s are the edge weights. Let
wmax =

max

(i,j)∈E1 ∪E2 ∪E1+2

w(i,j) .

According to our construction of B in Algorithm 1,
wmax ≤ α nnz(B) kAk2∞ .

(15)


Thus, λmax B > B is at most
P
P
P
2
2
2
(i,j)∈E1 (u i − u j ) +
(i,j)∈E2 (v i − v j ) +
(i,j)∈E1+2 (u i − v i − (u j − v j ))
P 2
.
wmax max
2
u
i ui + vi
We upper bound each term of the numerator. Let d1 , d2 , d1+2 be the maximum vertex degree of
the graphs G1 , G2 , G1+2 , respectively. By Cauchy-Schwarz inequality,
!
X
X
X
u 2i + u 2j ≤ 2d1
u 2i + v 2i ,
(u i − u j )2 ≤ 2
(i,j)∈E1

i

(i,j)∈E1

25

!
X

2

(v i − v j ) ≤ 2d2

X

u 2i

+

v 2i

,

i

(i,j)∈E2

and
X

X

(u i − v i − (u j − v j ))2 ≤ 4

u 2i + v 2i + u 2j + v 2j

(i,j)∈E1+2

(i,j)∈E1+2

!
X

≤ 4d1+2

u 2i

+

v 2i

.

i

Plugging the above inequalities to the expansion, we have


λmax B > B ≤ wmax · 8 max{d1 , d2 , d1+2 }
≤ 8wmax nnz(B).
By the upper bound of nnz(B) in Lemma 4.1 and the upper bound of wmax in Equation (15), we
have




λmax B > B ≤ O α nnz(A)2 kAk2∞ log2 kAk∞ .
This completes the proof.
Recall that, for a vector x ∈ RnB , x nA +1:nA +g1 corresponds to the auxiliary u-variables with
non-zero coefficients in MC2 -gadgets, x nA +1+g1 +6i:nA +6+g1 +6i for 0 ≤ i ≤ (g2 − g1 )/6 corresponds
to the auxiliary v -variables with non-zero coefficients in a single MC2 -gadget, and x n+g2 +1:nB
corresponds to the variables with zero coefficient.
Lemma 4.12. Let 0 ≤  ≤ (100nB log kAk∞ )−1 . Let x B ∈ RnB satisfying
1. x B
nA +g2 +1:nB = 0,
2. x B
nA +g1 +6i+1 = 0, ∀0 ≤ i ≤ (g2 − g1 )/6, and
3.

Bx B

2

≤  xB

2

.

Then,
Ax

A
2



1 1/2 1/2
≤4 1+
nB  x A
α

2

.

Proof. We first show that under the conditions of the Lemma,
xA

∞

≥

1
xB
4

def

∞

.

(16)

Let δ =  x B 2 . By the 3rd condition, each entry of Bx B has absolute value at most δ.
We first show that all u-variables in x aux cannot be large. According to Algorithm 1, for each
row of B, all nonzero coefficients have same absolute value, which is at least 1. Based on this fact

26

and the MC2 -gadget constructed in Algorithm 2, we have
u jk + u lk
− 5δ
2
u jk + u lk
− 6δ
2
u jk − δ
u jk + u lk
− 6δ
2
u lk − δ

≤ ut ≤
≤ u t+3
≤ u t+4
≤ u t+5
≤ u t+6

u jk + u lk
+ 5δ
2
u j + u lk
≤ k
+ 6δ
2
≤ u jk + δ
u j + u lk
≤ k
+ 6δ
2
≤ u lk + δ

where u jk , u lk being paired-and-replaced, and all others are entries of x aux . By the triangle inequality,
1
|u t | , |u t+k | ≤ (|u jk | + |u lk |) + 6δ, ∀3 ≤ k ≤ 6.
2
Note that the sum of coefficients on both sides are equal. We can repeat this type of substitution
on the right hand side until u jk and u lk are variables of x A . At the ith iteration of Algorithm 1
line 20, we have
X
|u t | , |u t+k | ≤
αj x A
∀3 ≤ k ≤ 6,
j + 6rδ,
j

where αj ≥ 0 and

P

j

αj = 1. By the Hölder inequality,
|u t | , |u t+k | ≤ x A

∞

∀3 ≤ k ≤ 6.

+ 6rδ,

(17)

We then argue that all v -variables in x aux cannot be large. Note that by the 2nd condition, in
the MC2 -gadget, we have v t+1 = 0. According to the equations in Algorithm 2, we have
u t+3 − u t+4 − 3δ ≤ v t+2 ≤ u t+3 − u t+4 + 3δ
−δ ≤ v t+3 ≤ δ
u t+3 − u t+4 − 2δ ≤ v t+4 ≤ u t+3 − u t+4 + 2δ
u t+3 − u t+4 − 4δ ≤ v t+5 ≤ u t+3 − u t+4 + 4δ
−δ ≤ v t+6 ≤ δ
By the triangle inequality,
|v t+k | ≤ |u t+3 | + |u t+4 | + 4δ,

∀2 ≤ k ≤ 6.

By (17), at the rth iteration of Algorithm 1 line 20,
|v t+k | ≤ 2 x A

∞

+ 12rδ,

∀2 ≤ k ≤ 6.

Since there are at most log kAk∞ iterations, the above inequality together with Equation (17)
implies
kx aux k∞ ≤ 2 x A ∞ + 12δ log kAk∞ .
Adding x A

∞

on both sides and substituting δ =  x B
xB

∞

≤ 3 xA
≤ 3 xA

2

gives

+ 12 x B 2 log kAk∞
√
+ 12 nB x B ∞ log kAk∞ .
∞
∞

27

−1
√
Given  ≤ 100 nB log kAk∞
, we have Equation (16). This further implies
xB

2

√
≤ 4 nB x A

2

.

(18)

By Claim 4.2 with c B = 0, we have
r
Ax

A
2

≤

α+1
Bx B
α

2

.

By the 3rd condition and the bound (18), we have
Ax

A
2



1 1/2 1/2
nB  x A
≤4 1+
α

2

.

This completes the proof.
We use the above lemma to lower bound the smallest nonzero eigenvalue of B > B.



(A> A)}
.
Lemma 4.13. λmin B > B = Ω min{1,λmin
n2
B

Proof. Let
def

δ = min{1, λmin (A> A)},

(19)

and
def

 =

δ
.
1616(1 + α−1 )n2B

(20)

The goal is to prove that λmin (B > B) ≥ . Assume by contradiction, there exists an x ⊥ null(B)
such that
(x B )> B > Bx B < (x B )> x B .

(21)

We show a contradiction by case analysis according to whether x A is orthogonal to the null space
of A.
Case 1: Suppose x A ⊥ null(A).
Recall that x B
nA +1:nA +g1 corresponds to the auxiliary u-variables in the constraints,
B
x nA +g1 +1:nA +g2 corresponds to the auxiliary v -variables in the constraints, and x B
nA +g2 +1:nB correB
j
sponds to the variables with zero coefficient. By Lemma 4.10, we have x ⊥ q , ∀1 ≤ j ≤ nA − g2 ,
where q j is defined in (14). Thus, x nA +g2 +1:nB = 0.
def P
B
i
i
Let r =
0≤i≤(g2 −g1 )/6 x nA +g1 +6i+1 p , where p is defined in (13). By Lemma 4.10, r ∈
def

null(B). Let y B = x B − r . Similarly, we write
 A 
y
B
y =
,
y aux
where y A corresponds to original variables and y aux corresponds to auxiliary variables. Note
y A = x A . Since r ∈ null(B), we have
(y B )> B > By B = (x B )> B > Bx B .
28

Since x B ⊥ r , we have
(y B )> y B = (x B )> x B + r > r ≥ (x B )> x B .
By assumption (21), we have
(y B )> B > By B < (y B )> y B .
By our definition of y B , y B satisfies the conditions of Lemma 4.12. By Lemma 4.12 and the
definition of  in Equation (20), we have
xA

>

A> Ax A < δ x A

>

x A.

Since x A ⊥ null(A) and the definition of δ in Equation 19, we have


n
o
λmin A> A < δ = min 1, λmin (A> A) ,
which is a contradiction.
Case 2: Suppose x A ∈ null(A).
A ) is the extension of vector x A defined in Equation (12).
P Recall that p(x
i
0≤i≤(g2 −g1 )/6 αi p be a vector such that

Let r =

def

y B = x B − p(x A ) − r
satisfying y B
nA +g1 +6i+1 = 0, ∀0 ≤ i < (g2 − g1 )/6. By this definition, we have
y A = 0.

(22)

Since p(x A ) + r is in the null space of B, by assumption (21),
(y B )> B > By B < (y B )> y B .
Note y B satisfies the conditions of Lemma 4.12. By Equations (16) and (22), we have
yB

∞

= 0.

This implies that x B = p(x A ) + r is in the null space of B, which contradicts that x B ⊥ null(B).
Case 3: Suppose x A = z̃ + ẑ , where z̃ 6= 0 is in null(A) and ẑ 6= 0 is orthogonal to null(A).
P
def
Similarly, let r 1 = i αi p i such that y B = x B − r 1 satisfying y B
nA +g1 +6i+1 = 0, ∀i. By this
definition, we have
y A = x A.
Then by assumption (21),
(y B )> B > By B < (y B )> y B .
y B satisfies the conditions of Lemma 4.12. By Lemma 4.12, we have


1
A > >
A
nB (y A )> y A .
(y ) A Ay ≤ 16 1 +
α
Since y A = z̃ + ẑ with z̃ ∈ null(A) and ẑ ⊥ null(A), we have




1
> >
ẑ A Aẑ ≤ 16 1 +
nB  z̃ > z̃ + ẑ > ẑ .
α
29

If ẑ > ẑ > z̃ > z̃ /(100nB ), then we have a contradiction with (19) and we have done. Otherwise, we
have
kz̃ k22
.
kẑ k22 ≤
100nB
P
Let r 2 = i βi p i such that
def

z B = y B − p(z̃ ) − r 2
satisfies z B
nA +g1 +6i+1 = 0, ∀i. By assumption (21),
(z B )> B > Bz B < (z B )> z B .
Vector z B satisfies the conditions of Lemma 4.12. By Equation (16), we have
zB

2

√
≤ 4 nB kẑ k2 .

Since ẑ > ẑ ≤ z̃ > z̃ /(100nB ), we have
zB

2

≤ kz̃ k2 .

On the other hand,
2
2

zB

= xB

2
2

+ kp(z̃ ) + r 1 + r 2 k22

> kp(z̃ ) + r 1 + r 2 k22
≥ xA

2
2

= kz̃ k22 + kẑ k22
> kz̃ k22 .
The first equality is due to x B ⊥ (p(z̃ ) + r 1 + r 2 ), the third inequality is due to the x A part of
r 1 + r 2 is 0, and the fourth equality is due to that z̃ ⊥ ẑ . Thus, we get a contradiction.
This completes the proof.
Lemma 4.11 and Lemma 4.13 immediately imply the following lemma.


nnz(A)2 kAk∞ log2 kAk∞
.
Lemma 4.14. κ(B) = O
min{1,σmin (A)}

4.5

Putting it All Together

Proof of Lemma 3.5. We set α = 1 in Algorithm 1.
By Lemma 4.1 and Claim B.1, we have
nnz(B) = O (s log(sU )) .
According to our reduction in Algorithm 1, the largest entry and the smallest nonzero entry of the
right hand side vector does not change. Besides, all nonzero entries of B have absolute value at
least 1. The largest entry of B is upper bounded by Equation (15)
p

√
wmax = O
nnz(B) kAk∞ .
By Claim B.1, we have

√



wmax = O s3/2 U log1/2 (sU ) .
30

By Lemma 4.14 and Lemma B.1, B.3, we have



nnz(A)2 kAk∞ log2 kAk∞
= O s4 U 2 K log2 (sU ) .
κ(B) = O
min{1, σmin (A)}
By Lemma 4.9, we have
(B )−1 = (A )−1 O 1 + c A



A −1
2
σ
(A)
=
(
)
O
sU
.
max
2

This completes the proof.

5

MC2 Efficiently Reducible to MC2>0

The matrices generated by interior point methods (see Section 9 for details) are more restrictive
than MC2 : for every edge present, the weights of all three types of matrices are non-zero. Formally,
the class of matrices MC2>0 consists of matrices of the form






1 −1
0 0
1 0
2
1
1+2
,
+L
⊗
+L ⊗
L ⊗
−1 1
0 1
0 0
where the non-zero support of all three matrices, L1 , L2 , and L1+2 are the same (but the weights
may vary greatly). On the other hand, the matrices that we generate in Section 4 can be transformed
into such a matrix by adding a small value, δ to the weight of all edges where one of the types
of edges have non-zero support. We will describe this construction in Section 5.1, and bound its
condition number in Section 5.2.

5.1

Construction

In this section, we show the construction from an instance of MC2 to an instance of MC2>0 . The
strategy is to add extra edges with a sufficiently small weight, such that L1 , L2 , L1+2 have identical
nonzero stricture and the solution of the linear system does not change much.
The reduction from MC2 to MC2>0 , with pseudo-code in Algorithm 4, simply adds edges with
weight δ to all the missing edges. The transformation of solutions of the corresponding instance
of MC2>0 back to a solution of the original MC2 instance in Algorithm 5 simply returns the same
vector x .

31

Algorithm 4 ReduceMC2 ToMC2>0
Input: (B, c B , 1 ) where B ∈ MC2 is an m × n matrix, c B ∈ Rm , and 0 < 1 < 1.
>0
>0
0
Output: (B >0 , c B , 2 ) where B >0 ∈ MC2>0 is an m0 × n0 matrix, c B ∈ Rm , and 0 < 2 < 1.

1: δ ← 100κ(B)σ 1 (B)kc B k
max
2
b← ∅
2: B
{New equations for MC2>0 }
3: for each pair of vertices i, j whose blocks are involved in some equation in B do
4:
if B does not contain a type 1 equation between i and j then
5:
Bb ← Bb ∪ {u i − u j = 0}.
6:
end if
7:
if B does not contain a type 2 equation betwen i and j then
8:
Bb ← Bb ∪ {v i − v j = 0}.
9:
end if
10:
if B does not contain a type 1 + 2 equation betwen i and j then
11:
Bb ← Bb ∪ {u i + v i − (u j + v j ) = 0}.
12:
end if
13: end for
e be the coefficient matrix of equations in B.
b
14: Let B
15: return

 B 

B
1
c
B>0
>0
,
2 =
,
c
=
B =
.
e
0
δB
100

Algorithm 5 MapSolnMC2>0 ToMC2
Input: m × n matrix B ∈ MC2 , m0 × n0 matrix B >0 ∈ MC2>0 , vector c B ∈ Rm vector x ∈ Rn .
Output: Vector y ∈ Rn .
1: if B > c B = 0 then
2:
return y ← 0
3: else
4:
return y ← x
5: end if
>0

Let B >0 x = c B be the linear system returned by a call to ReduceMC2 ToMC2>0 (B, c B , 1 ).
As we only add up to two edges per original edge in B, and do not introduce any new variables,
>0
the size of B >0 is immediate from this routine. According to the algorithms, c B is simply
 B   A 
c
c
B>0
c
=
=
.
0
0
Note the two zero vectors in the above equation have different dimensions.
For a given matrix B, we will choose the additive term to ensure non-zeros to be:
def

δ =

1
.
100κ(B)σmax (B) kc B k2

(23)

Note that we can bound κ(B), σmax (B) by condition number of A (i.e, the linear system instance
of G).
Lemma 5.1. nnz(B >0 ) = O (nnz(B)). By our setting of δ in (23), the largest entry of B >0 does
not change, the smallest entry of B >0 is at least δ.
32

Note that the addition of all three types of edges means that the null space of B >0 is now given by
the connected components in its graph theoretic structure, and is likely significantly different from
the null space of B. In order to solve the Linear System Approximation Problem(lsa) problems
min Bx − c B
x

,

2

min B >0 x − c B

>0

x

2

,

we need to solve the two linear systems:
B > Bx = B > c B ,

(24)

(B >0 )> B >0 x = (B >0 )> c

B>0

.

(25)

First note that the RHS the two equations are the same because:

  cB 
>
>0 > B>0
>
e
(B ) c
= B
= B >cB.
δB
0
This means that for a sufficiently small choice of δ, the differences between the solutions of
these two linear systems is small.
>0

Lemma 5.2. Let B >0 x
=
>0
B
ReduceMC2 ToMC2 (B, c , 1 ).
>0 ∗

Let x B

cB
be the linear system returned by a call to
B∗
Let x
∈ arg minx Bx − c B 2 and x B∗ ⊥ null(B).

>0

∈ arg minx B >0 x − c B

2

. Then we have:

>0 ∗

x B∗ − x B

>

B B

≤ O δ(1 + δ)κ(B) c B


2

.

Proof. The desired distance can be written as
>0 ∗

x B∗ − x B

B >B



>0
= B x B∗ − x B ∗

2

.

Also, the optimality of x B∗ means that Bx B∗ − c B is perpendicular to anything in the rank space
of B, in particular,


>0
B x B∗ − x B ∗ ⊥ Bx B∗ − c B ,
which in turn gives:


>0
B x B∗ − x B ∗

2
2

>0 ∗

= Bx B

>0 ∗

We now bound the right hand side. Since x B
B >0 x B

>0 ∗

>0

− cB

2
2

− cB

2
2

− Bx B∗ − c B

2
.
2

(26)

is a minimizer of B >0 x − c B
>0

≤ B >0 x B∗ − c B

>0

2
2

,

e term in B >0 separately to get:
where we can extract out the δ B
>0 ∗

Bx B

− cB

2
2

e B>0 ∗
+ δ 2 Bx

2
2

≤ Bx B∗ − c B

33

2
2

e B∗
+ δ 2 Bx

2
2

.

2

, we have

Together with Equation (26), this then gives
>0 ∗

x B∗ − x B
e B>0 ∗
It remains to upper bound Bx

2

B >B

e B∗
≤ δ Bx

.

2

. By the assumption of x B∗ ⊥ null(B), we get:
2

 > 
e B
e
λmax B
2
.
≤
λmin B > B
kBx B∗ k22
e B∗
Bx


By a proof similar to Lemma 4.11, which upper bounds λmax B > B , we have:
 
e = O(1 + δ)σmax (B),
σmax B
which implies
e B∗
Bx

2

≤ O(1 + δ)κ(B) Bx B∗

Since x B∗ is a minimizer of minx Bx − c B
Bx B∗

2

2

2

.

, we have

= ΠB c B

2

≤ cB

2

.

Therefore,
>0 ∗

x B∗ − x B

≤ O (δ(1 + δ)κ(B)) kck2 ,

B >B

which completes the proof.
We now check that the approximate solutions of the two linear systems (24) and (25) are also
close to each other.
>0

Lemma 5.3. Let B >0 x
=
cB
be the linear system returned by a call to
>0
B
ReduceMC2 ToMC2 (B, c , 1 ), and let 2 be the error parameter returned by this call.
Let x be a vector such that
>0

>0

B >0 x − ΠB >0 c B

2

≤ 2 ΠB >0 c B

2

.

Let y be the vector returned by a call to MapSolnMC2>0 ToMC2 (B, B >0 , c B , x ). Then
By − ΠB c B

2

≤ 1 ΠB c B

2

.

Proof. In Algorithm 5, if the condition B > c B = 0 is true, then y = 0. This implies
By − ΠB c B

2

= 0.

In the following, we assume that B > c B 6= 0, in which case y = x . We first show that the
>0
>0
construction implies Bx − ΠB c B 2 ≤ B >0 x − ΠB >0 c B
. Once again, let x B∗ and x B ∗ be
2

the vectors that minimize Bx − c B

2

and B >0 x − c B

Bx − ΠB c B

2

>0

= x − x B∗
34

2

respectively. This choice gives:
B >B

.

As k·kB > B is a norm, by triangle inequality we have:
x − x B∗

>0 ∗

≤ x − xB

B >B

>

B B

+ xB

>0 ∗

− x B∗

B >B

.

We will bound these two terms separately.
e >B
e = (B >0 )> B >0 , the first term is less than its norm in the
Since B > B 4 B > B + δ 2 B
>0 > >0
(B ) B norm:
>0 ∗

>0 ∗

x − xB

>

B B

≤ x − xB



>

>

e B
e
B B+δ 2 B

= B >0 x − ΠB >0 c B



>0

2

,

while the second term is precisely the distances between the two optimums, which we just bounded
in Lemma 5.2. Combining these bounds then gives:
Bx − ΠB c B

>0

2

≤ B >0 x − ΠB >0 c B

2

+ O δ(1 + δ)κ(B) c B


2

.

(27)

As the equations in B >0 is a superset of the ones in B, we have
>0 ∗

B >0 x B
Substituting Bx B∗ − c B

2
2

>0

>0 ∗

− cB

2
2
2

= cA

≥ Bx B

2
2

− ΠB c B
>0

ΠB >0 c B

2

− cB

2

≥ Bx B∗ − c B

2

.

and its equivalent in B >0 gives

≤ ΠB c B

2

.

Together with the condition of the lemma,
>0

B >0 x − ΠB >0 c B

2

≤ 2 ΠB c B

2

.

Plugging this into Equation (27), we have
Bx − ΠB c B
It remains to upper bound c B
and assumption,

2
Bc B 2

2

2

= O 2 ΠB c B

+ δκ (B) c B
2

by a function of ΠB c B

2


2

.

. By our construction in Algorithm 1

is a positive integer. By Lemma 4.8,
ΠB c B

2

≥

1
.
σmax (B)

Thus,
Bx − ΠB c B

2

= O 2 + δκ(B)σmax (B) c B

By our setting of δ in Equation (23), we have
Bx − ΠB c B

2

≤ 1 ΠB c B

This completes the proof.

35

2

.


2

ΠB c B

2

.

5.2

Bounding Condition Number of the New System in MC2>0

We now establish bounds on the numerical quantities related to B >0 . By a proof similar to the
upper bound on λmax B > B in Lemma 4.11, we have:





λmax (B >0 )> B >0 = O λmax B > B .
As a result, we focus on the lower bound here:
Lemma 5.4. The matrix B >0 from the linear system returned by a
ReduceMC2 ToMC2>0 (B, c B , 1 , δ), where δ is set according to Equation 23, satisfies



1
>0
σmin B
=Ω
.
κ(B)σmax (B)nB

call

to

Proof. Let G be a unit-edge weight graph whose vertex set and edge set are same as the underlying
graph of B >0 . Let LG be the associated Laplacian matrix of G, and let M := LG ⊗ C , where


2 −1
C =
−1 2
is a symmetric PSD matrix. Note






1 −1
0 0
1 0
,
+ LG ⊗
+ LG ⊗
M = LG ⊗
−1 1
0 1
0 0
which means M has the same null space as (B >0 )> B >0 .


λmin (B >0 )> B >0 =

x > (B >0 )> B >0 x
x >M x
2
≥
δ
min
,
x >x
x ⊥null(M ) x > x
x ⊥null(B >0 )
min

the last inequality is due to that δ is the minimum edge weight.
Note that the sum of the 3 types of blocks is positive definite and has eigenvalue at least 1: the
type 1 and 2 blocks already sum to I . Formally:
λmin (M ) = λmin (LG ) · λmin (C ) = λmin (LG ).
The result then follows from the folklore bound that the minimum non-zero eigenvalue of a
unit weighted graph is at least n12 . One way to see this is via Cheeger’s inequality (see e.g. [Spi07]
applied to each block: this decomposition is equivalent to spearating the matrix into its diagonal
blocks based on the connected components, and then invoking the fact that the minimum weight
of
 
1
a cut is at least 1, and there are at most n vertices. Together these tools imply λmin (LG ) = Ω n2 ,
B
and in turn:


 

δ
1
σmin B >0 ≥ Ω
=Ω
.
n
κ(B)σmax (B)nB
This also implies a bound on the condition number of B >0 :

2
Lemma 5.5. κ(B >0 ) = O −1
1 σmax (B)κ(B)nB .

36

5.3

Putting it All Together

Proof of Lemma 3.6. By Lemma 5.1, we have
nnz(B >0 ) = O(s),
The largest entry is U , and the smallest nonzero entry is δ. By Equation (23) and Lemma B.2
and B.3,
  
.
δ=Ω
KU 2
By Lemma 5.5 and Lemma B.2 and B.3, the condition number is
 2 2 
s U K
>0
κ(B ) = O
.

By Lemma 5.3, the accuracy error for B >0 is O().

6

Rounding and Scaling Weights to Integers

In this section, we show a reduction from the linear system with strict 2-commodity matrix to the
linear system with integral strict 2-commodity matrix.
>0
Algorithm 6 presents the pseudo-code for the algorithm ReduceMC2>0 ToMC2,Z
. Given an
>0

>0

>0
(B >0 , c B , 1 ) reinstance (B >0 , c B , 1 ) where B >0 ∈ MC2>0 , the call ReduceMC2>0 ToMC2,Z
>0,Z

>0
. Algorithm 7 provides the (trivial)
turns an instance (B >0,Z , c B
, 2 ) where B >0,Z ∈ MC2,Z
>0
>0
to a
ToMC2>0 which maps a solution of an instance over MC2,Z
pseudo-code for MapSolnMC2,Z
>0
solution of an instance over MC2 .
>0
Algorithm 6 ReduceMC2>0 ToMC2,Z
>0

>0

Input: (B >0 , c B , 1 ) where B >0 ∈ MC2>0 is an m × n matrix, c B ∈ Rm , and 0 < 1 < 1
>0,Z
>0,Z
0
>0
is an m0 × n0 matrix, c B
∈ Rm , and 0 <
Output: (B >0,Z , c B
, 2 ) where B >0,Z ∈ MC2,Z
2 < l1.
m
nnz(B >0 )
1: k ← log2
+2
1
2:
3:
4:
5:
6:
7:
8:
9:
10:

B >0,Z ← B >0
>0,Z
>0
cB
← cB
>0,Z
for each entry (B
 >0,Z)ij do k 
>0,Z
(B
)ij ← (B
)ij · 2
end for
>0,Z
for each entry c B
do
i
>0,Z
>0,Z
B
B
k
ci
← ci
·2
end for
>0,Z
return (B >0,Z , c B
, 1 /3).
>0

For simplicity, we analyze an intermediate linear system B int y = c B , where
def

(B int )ij = 2−k · (B >0,Z )ij .
Recall the definition of B >0,Z in Algorithm 6 line 5, we have
l
m
(B int )ij = 2−k · (B >0 )ij · 2k .
37

(28)

>0
Algorithm 7 MapSolnMC2,Z
ToMC2>0
0

Input: m × n matrix B >0 ∈ MC2 , m0 × n0 matrix B >0,Z ∈ MC2>0 , vector x ∈ Rn .
Output: Vector y ∈ Rn .
1: return x
>0,Z

>0

The linear system B int y = c B is exactly the linear system B >0,Z y = c B
multiplying a factor
2−k on both sides.
Note the condition number and the eigen-space of a matrix, the optimal solutions of the
>0,Z
projection problems are invariant under scaling. Let y ∗ ∈ arg miny B >0,Z y − c B
=
2

arg miny B

int

y−

>0
cB

2

. The following two inequalities:
ky − y ∗ k(B >0,Z )> B >0,Z ≤  ΠB >0,Z c B

and

>0,Z

2

>0

ky − y ∗ k(B int )> B int ≤  ΠB int c B

2
>0

are equivalent. Thus, it suffices to analyze the linear system B int y = c B .
Solving the two projection problems
>0

min B >0 x − c B
x

min B int y − c B

2
>0

y

2

,
,

is equivalent to solving the following two linear systems:
>0

(B >0 )> B >0 x = (B >0 )> c B ,
(B int )> B int y = (B int )> c

B>0

(29)

.

>0

Note c B = (c A ; 0), and all entries of the rows of B >0 corresponding to the original rows of A
are integers. Thus,
>0
>0
(B int )> c B = (B >0 )> c B .
That is, the 2nd linear system is equivalent to
>0

(B int )> B int y = (B >0 )> c B .
Let

(30)

def
c def
M = (B >0 )> B >0 and M
= (B int )> B int .

c by the following lemma.
We bound the eigenvalues of M
c ) ≤ (1 + 2−k+2 )λmax (M ) and λmin (M
c ) ≥ λmin (M ). Furthermore, κ(M
c) ≤
Lemma 6.1. λmax (M
−k+2
(1 + 2
)κ(M ).
Proof. Note M can be written as B > W B, where B is the incidence-structured block matrix with
unit nonzero entries and W is the diagonal matrix with edge weights. Let x be an arbitrary vector,
def
and let y = Bx .
X
x >M x = y >W y =
W ii y 2i ,
i

38

and
cx =
x >M

X

c iy 2.
W
i

i

Let
def

α = min
i

c ii
c ii
W
W
def
and β = max
.
i W ii
W ii

We have
α

X

W i y 2i ≤

X

i

c iy 2 ≤ β
W
i

X

i

W i y 2i .

i

This implies that
c 4 βM .
αM 4 M

(31)

Now we bound the values of α, β. According to Equation (28),
2

−k · W 1/2 · 2k
2
c
ii
W ii
α = min
≥
= 1.
i W ii
W ii

(32)

Similarly,

β = max
i

c ii
W
≤
W ii

2−k



1/2
W ii

·

2k

+1

2
≤ 1 + 2−k+2 .

W ii

(33)

The last inequality is due to W ii ≥ 1. Thus,
c ) ≤ (1 + 2−k+2 )λmax (M ) and λmin (M
c ) ≥ λmin (M ).
λmax (M
c ) ≤ (1 + 2−k+2 )κ(M ) immediately follows the above two
The bound on the condition number κ(M
inequalities.
We show the exact solutions of the two linear systems in Equation (29) and (30) are close, by
the following lemma.
Lemma 6.2. Let

def

x ∗ = M † (B >0 )> c B

>0

>0
def c †
and y ∗ = M
(B >0 )> c B .

x ∗ and y ∗ are exact solutions of linear system (29) and (30) respectively.
>0

ky ∗ − x ∗ kM ≤ 2−k+2 (B >0 )> c B

M†

.

c are both symmetric. Expanding the left hand side,
Proof. Note M , M
 †

 †

>0
c − M † (B >0 )> c B>0
c −M† M M
ky ∗ − x ∗ k22 = (c B )> B >0 M

2
>0
c † M 1/2 − I M †1/2 (B >0 )> c B>0
= (c B )> B >0 M †1/2 M 1/2 M
†

c M 1/2 − I
≤ M 1/2 M

2
2

>0

M †1/2 (B >0 )> c B

The last inequality is by the Courant-Fischer theorem.
39

2
2

.

c have same null space,
By Equation (31) and the fact that M , M


c † M 1/2 − I 4 α−1 − 1 I .
β −1 − 1 I 4 M 1/2 M
By Equation (32) and (33),
†

c M 1/2 − I
M 1/2 M

2

≤ 2−k+2 .

Therefore,
>0

ky ∗ − x ∗ kM ≤ 2−k+2 (B >0 )> c B

M†

.

We then show the approximate solutions of the two linear systems in Equation (29) and (30)
are close.
Lemma 6.3. Let x be a vector such that
>0

B int x − ΠB int c B

>0

2

≤ 2 ΠB int c B

2

.

Then,


>0
≤ 2 + 2−k+2 (1 + 2 ) ΠB >0 c B

>0

B >0 x − ΠB >0 c B

2

2

.

Proof. Let
def

>0

x ∗ = M † (B >0 )> c B

>0
def c †
and y ∗ = M
(B >0 )> c B .

Expanding the left hand side,
B >0 x − ΠB >0 c B

>0

2

= kx − x ∗ kM
≤ kx − y ∗ kM + ky ∗ − x ∗ kM .

The last inequality is due to the triangle inequality. By Equation (31), the first term can be upper
bounded by
kx − y ∗ kM ≤ α−1/2 kx − y ∗ kM
c.
The second term is upper bounded by Lemma 6.2. Thus, we have
>0

B >0 x − ΠB >0 c B

>0

≤ α−1/2 2 ΠB int c B

2

>0

+ 2−k+2 ΠB >0 c B

2

2

.

Note that
ΠB int c B

>0

>0

2

= B >0 c B

c
M

†

and

ΠB >0 c B

>0

>0

2

= B >0 c B

M†

Again, by Equation (31),
B >0 c B

>0

>0

c†
M

≤ α−1/2 B >0 c B

M†

.

By Equation (32) and (33),
>0

B >0 x − ΠB >0 c B

2



>0
≤ 2 + 2−k+2 (1 + 2 ) ΠB >0 c B

This completes the proof.
40

2

.

.

Proof of Lemma 3.7. We set

l
sm
def
k = log2
+ 2.

After rounding and scaling, the number of nonzero entries does not change. The value of the
smallest nonzero entry does not decrease. The value of the largest entry is multiplied by 2k = s−1 ,
which is upper bounded by s−1 U .
By Lemma 6.1,
κ(B >0,Z ) ≤ (1 + s−1/2 1/2 )κ(B >0 ).
By Lemma 6.3, the accuracy is bounded by /3.

G Efficiently Reducible to Gz,2

7

In this section, we prove Lemma 3.4.
Definition 7.1. We let Gz denote the class of all matrices with integer valued entries s.t. there is
at least one non-zero entry in every row and column, and every row has zero row sum.

7.1

G ≤f Gz

In this section, we show the reduction from G to Gz . Given an instance of linear system (A, c A , 1 )
with A ∈ G, the goal is to construct an instance of linear system (AZ , c Z , 2 ) with AZ ∈ Gz such
that, there exists a map between the solutions of the two linear systems.
Algorithm 8 shows the construction of (AZ , c Z , 2 ), and Algorithm 9 shows the transform from
a solution of (AZ , c Z , 2 ) to a solution of (A, c A , 1 ). Recall that A has nA columns. According to
the algorithm, AZ has (nA + 1) columns, and c Z = c.
Algorithm 8 Reduce GtoGz
Input: (A, c A , 1 ) where A1 ∈ G is an m × n matrix, c A ∈ Rm , and 0 < 1 < 1.
0
Z
Output: (AZ , c Z , 2 ) where
∈ Gz 
is an m0 × n0 matrix, c Z ∈ Rm , and 0 < 2 < 1.
 A
A −A1 , c A , O n1 .
1: return

Algorithm 9 MapSoln Gz toG
0

Input: m × n matrix A ∈ G, m0 × n0 matrix AZ ∈ Gz , vector x Z ∈ Rn .
Output: Vector x ∈ Rn .
Z
1: return x ← x Z
1:n − x n+1 1
Lemma 7.2. Let AZ x Z = c Z be the linear system returned by a call to Reduce GtoGz (A, c A , 1 ).
Then,
nnz(AZ ) = O (nnz(A)) ,
and the largest entry of AZ is at most kAk∞ .
We show the relation between the exact solvers of the two linear systems.
Claim 7.3. Let x Z∗ ∈ arg minx AZ x − c Z

2

. Let x ∗ be the output vector of Algorithm 9, that is,

Z∗
x ∗ = x Z∗
1:nA − x nA +1 1.

41

Then,
x ∗ ∈ arg min Ax − c A

.

2

x

Proof. Note 1 ∈ N ull(AZ ), thus,

AZ x Z∗ = AZ x Z∗ − x Z∗
nA +1 1 .
Expanding it gives
AZ x Z∗ =



A −A1

Note c Z = c. Therefore, x ∗ ∈ arg minx Ax − c A

2



x∗
0



= Ax ∗ .

.

We then show the relation between the approximate solvers.
Lemma 7.4. Let x Z be a vector such that AZ x Z − ΠAZ c A
vector of Algorithm 9, that is, x = x Z1:nA − x ZnA +1 1. Then,
Ax − ΠA c

2

√

A
2

nA σmin (A)
σmin (AZ )

≤O



≤ 2 ΠAZ c A

2 ΠA c A

2

2

. Let x be the output

.

Proof. Let x Z∗ ∈ arg minx AZ x − c Z whose last entry is 0. Expanding the left hand side norm,
AZ x Z − ΠAZ c A
By Claim 7.3, x ∗ ∈ arg minx Ax − c A

2

2

= AZ x Z − x Z∗


2

.

. Thus,

kA (x − x ∗ )k2 = AZ (x Z − x Z∗ )

2

.

Note that
Ax = AZ x Z and Ax ∗ = AZ x Z∗ .
Together with the Lemma condition, this gives
Ax − ΠA c A

= AZ x Z − ΠAZ c Z

2

It remains to upper bound ΠAZ c A
ΠA c A

2

= A> c A

It suffices to work on A> c A

†

have
(AZ )> c A

by a function of ΠA c A

†

and

ΠAZ c A

and (AZ )> c A
2
2

≤ 2 ΠAZ c A

2

(A> A)

(A> A)

2

= A> c A

2
2

2

2

†

†

. Since AZ =

+ 1> A> c A

2
2

2

≤ k1k2 A> c A

42

2

=

√

(34)

((AZ )> AZ )

.

By Cauchy-Schwarz inequality,
1> A > c A

.

. Note

= (AZ )> c A

((AZ )> AZ )

2

nA A > c A

2

.

.

(35)


A −A1 , we

Plugging this gives
(AZ )> c A

2

≤ (nA + 1) A> c A

2

2
2

.

By Courant-Fischer theorem,
Z > A

(A ) c

2
†
(AZ )> AZ

(

)

≤ λmax



Z >

(A ) A

Z

† 

2

(AZ )> c A

and since A> c A is in the eigenspace of A> A,

† 
2
> A
>
A c
A A
A> c A
† ≥ λmin
(A> A)

2
2

2



Z > Z
(AZ )> c A
≤ λ−1
(A
)
A
min



>
= λ−1
A> c A
max A A

2
2

2
2

.

Combining the above two inequalities,
Z > A

(A ) c

2
†
(AZ )> AZ

(

)

≤

λmax A> A
λmin



(AZ )> AZ

Given equations (34) and (35), and condition AZ x Z − c A
Ax − ΠA c A

2

≤ 2

2

 (nA + 1) A> c A
2

(A> A)

≤  ΠAZ c A

σmax (A) √
nA + 1 ΠA c A
σmin (AZ )

2

2

†

.

, we have

.

This completes the proof.
We compute the nonzero condition number of AZ . Note


M
−M 1
Z > Z
(A ) A =
,
−1> M 1> M 1

(36)

where M = A> A.

Claim 7.5. λmax (AZ )> AZ ≤ 2nA λmax (A> A).
Proof. Let λ1 be the largest eigenvalue of (AZ )> AZ , and y = (ỹ ; α) be the associated eigenvector
of unit length. Let µ1 be the largest eigenvalue of A> A = M .




M
−M 1
ỹ
>
Z > Z
>
λ1 = y (A ) A y = ỹ
α
α
−1> M 1> M 1
= ỹ > M ỹ − 2α1> M ỹ + α2 1> M 1.
By Courant-Fischer Theorem,
ỹ > M ỹ ≤ µ1 ỹ > ỹ , and 1> M 1 ≤ µ1 nA .
By Cauchy-Schwarz inequality,
1> M ỹ ≤ k1k2 kM ỹ k2 ≤

√

nA µ1 kỹ k2 .

Putting all together,


√
λ1 ≤ µ1 kỹ k22 + 2 |α| nA kỹ k2 + α2 nA
√
= µ1 (kỹ k2 + |α| nA )2 .
Since kỹ k2 , |α| ≤ 1, we have
λ1 ≤ 2µ1 nA .
This completes the proof.
43

Before lower bounding the smallest nonzero singular value of AZ , we characterize the null space
of AZ by the following Claim.
Claim 7.6.

 


z
null(A ) = Span 1,
: z ∈ null(A) .
0
 


z
Proof. Let S be the subspace of Span 1,
: z ∈ null(A) .
0
We first show that S ⊆ null(AZ ). Clearly,

AZ 1 = A −A1 1 = 0.
Z

For each z ∈ null(A), we have


z
Z
A
=
0

A −A1





z
0


= Az = 0.

Thus, S ⊆ null(AZ ).
We then show that S ⊇ null(AZ ). For any z 0 ∈ null(AZ ), z 0 − z 0nA +1 1 ∈ null(AZ ). Let z ∈ Rn
such that the ith entry of z is z 0i − z 0nA +1 . Thus,
0=A

Z

0

z −

z 0nA +1 1



=

A −A1





z
0


= Az .


That is, any vector in S can be written as a linear combination of 1 and

z
0


where z ∈ null(A).

Thus, S ⊇ null(AZ ).
Therefore, null(AZ ) = S .
By the above claim, we know that A and AZ have same rank. We bound the smallest nonzero
singular value of AZ by the following Claim.

Claim 7.7. λmin (AZ )> AZ ≥ λmin (A> A)/(nA + 1).
Proof. Let λk be the smallest nonzero eigenvalue of (AZ )> AZ , and y = (ỹ ; α) be the associated
eigenvector of unit length. Let µk be the smallest nonzero eigenvalue of A> A = M , and x be the
associated eigenvector of unit length.
λk = ỹ > M ỹ − 2α1> M ỹ + α2 1> M 1.
Since M is symmetric PSD, it can be written as M = M 1/2 M 1/2 .
λk = M 1/2 (ỹ − α1)

2
2

.

We decompose ỹ −α1 := z 1 +z 2 , where z 1 ∈ null(M ) and z 2 ⊥ null(M ). Then, λk = M 1/2 z 2
By Courant-Fischer Theorem
λk ≥ µk kz 2 k22 .
To lower bound λk , it suffices to lower bound kz 2 k2 .
44

2
2

.

Since y = (ỹ ; α),
kỹ −

α1k22


= y −α

1
1



2

.
2

Since y ⊥ 1, which is in null(AZ ), we have
kỹ − α1k22 = ky k22 + α2 (1 + nA ) = kỹ k22 + α2 (2 + nA ).
Vectors ỹ , −α1, ỹ − α1 form a triangle. Let θ be the angle between ỹ and α1. See Figure 1.

Figure 1: The relation between vectors ỹ , −α1 and ỹ − α1.
Then,
cos(θ) =

kỹ k22 + kα1k22 − kỹ − α1k22
|α|
=−
√ .
2 kỹ k2 kα1k2
kỹ k2 nA

Thus,
kz 2 k2 = kỹ k2 + kα1k2 cos(π − θ) = kỹ k2 +

α2
≥ kỹ k2 .
kỹ k2

On the other hand,
!2
1 = α2 +

kỹ k22

=

X

ỹ i

i

+

X

ỹ 2i ≤ (nA + 1)

i

X

ỹ 2i = (nA + 1) kỹ k22 .

i

Therefore,
λk ≥ µk kỹ k22 ≥

µk
.
nA + 1

This completes the proof.
The above lemmas give an upper bound on the condition number of AZ .
3/2

Lemma 7.8. κ(AZ ) ≤ O(nA )κ(A).

45

7.2

Gz ≤f Gz,2

In this section, we show the reduction from Gz to Gz,2 . Given an instance of linear system (AZ , c Z , 1 )
with AZ ∈ Gz , the goal is to construct an instance of linear system (AZ,2 , c Z,2 , 2 ) with AZ,2 ∈ Gz,2
such that, there is a map between the solutions of these two linear systems.
Algorithm 10 shows the construction of (AZ,2 , c Z,2 , 2 ), and Algorithm 11 shows the transform
from a solution of (AZ,2 , c Z,2 , 2 ) to a solution of (AZ , c Z , 1 ). Note c Z,2 = (c Z ; 0).
Algorithm 10 Reduce Gz toGz,2
Input: (AZ , c Z , 1 ) where AZ ∈ Gz is an m × n matrix, c Z ∈ Rm , and 0 < 1 < 1.
0
Output: (AZ,2 , c Z,2 , 2 ) where AZ,2 ∈ Gz,2 is an m0 × n0 matrix, c Z,2 ∈ Rm , and 0 < 2 < 1.
√
−1
mσmax (AZ ) AZ ∞ c Z 2
1: w ← 31

2: k ∗ ← min k ∈ Z : 2k ≥ AZ ∞
3: Let a ∈ Rn
4: for i ← 1 to n do
∗
5:
a i ← 2k − AZi 1 /2
6: end for

AZ a −aa
Z,2
7: A
=
0 w −w
 Z 
c
8: c Z,2 ←
0



9: return
AZ,2 , c Z,2 , O σ (AZ1 )kc Z k
.
max

2

Algorithm 11 MapSoln Gz,2 toGz
0

Input: m × n matrix AZ ∈ Gz , m0 × n0 matrix AZ,2 ∈ Gz,2 , vector c Z ∈ Rm , vector x ∈ Rn .
Output: Vector y ∈ Rn .
1: if (AZ )> c Z = 0 then
2:
return y ← 0
3: else
4:
return y ← x 1:n
5: end if
Lemma 7.9. Let AZ,2 x Z,2 =
duce Gz toGz,2 (AZ , c Z , 1 ). Then,

c Z,2 be the linear system returned by a call to Re-


nnz(AZ,2 ) = O nnz(AZ ) .
√
and the largest entry of AZ,2 is at most −1
mσmax (AZ ) AZ ∞ c Z 2 .
1
  


0
x
Claim 7.10. null(AZ,2 ) = Span  1  ,  0  : x ∈ null(AZ ).
1
0
Proof. Let y = (x ; α; β) ∈ RnA +2 such that y ∈ null(AZ,2 ), that is,



 x
 Z

A a −aa
 α  = A x + (α − β)aa
AZ,2 y =
= 0.
0 w −w
(α − β)w
β
46

Since w > 0, we have
α = β.
and thus
AZ x = 0.
Therefore,


 


0
x
null(AZ,2 ) = Span  1  ,  0  : x ∈ null(AZ ) .
1
0
This completes the proof.
We write x Z,2 as (x ; α, β). By Claim 7.10, the vector (0; 1; 1) is in the null space of AZ,2 . Thus,

 


x
0
α + β  
.
1
AZ,2 x Z,2 = AZ,2 x Z,2 −
= AZ,2  α−β
2
2
β−α
1
2
Without loss of generality, we assume α = β, that is, x Z,2 = (x ; α; −α).
We first show that the exact solutions of the two linear systems are close.
Lemma 7.11. Let x Z,2∗ ∈ arg minx AZ,2 x − c Z,2

2

. Write x Z,2∗ = (s ∗ ; α; −α). Then,

s ∗ ∈ arg min AZ x − c Z
x

P

,

aa >
.
a k22
w2 +ka

where P = I −

Proof. Without loss of generality, we write x Z,2 = (x ; α; −α). We expand AZ,2 x Z,2 − c Z,2 ,
min AZ,2 x Z,2 − c Z,2

2
2

Z



x Z,2



A
0

= min
x ,α

a −aa
w −w

Z

Z

= min A x − c + 2αaa




 Z 
x
 α − c
0
−α
2
2

2

2

2 2

+ 4w α

2


Z
>
Z
a (A x − c ) 

= min min 4 w2 + kaa k22 α + 
+ AZ x − c Z
2
x
α
2 w2 + kaa k2
x ,α

2
2

−

2
a > (AZ x − c Z )
w2 + kaa k22

.

Since the minimization over x and α is independent of each other, for any fixed x , the above value
is minimized when
a > (AZ x − c Z )
.
α=− 
2 w2 + kaa k22
Plugging this value of α gives
min A
x Z,2

Z,2

x

Z,2

−

2
c Z,2 2

Z

= min A x −
x

Z

2
cZ 2

= min A x − c

−

x

x

47

2
P

w2 + kaa k22

aa >
I−
w2 + kaa k22


Z >

= min AZ x − c Z

2
a > (AZ x − c Z )

,

!
AZ x − c Z



where P = I −

aa >
.
a k22
w2 +ka

The fact x Z,2∗ ∈ arg minx AZ,2 x − c Z,2
s ∗ ∈ arg min AZ x − c Z
x

P

2

implies that

,

which completes the proof.
Lemma 7.12. Let x Z,2∗ ∈ arg minx AZ,2 x − c Z,2 2 such that x Z,2∗ has the form (s ∗ ; α; −α) and
x Z,2∗ 2 is minimized. Let x ∗ ∈ arg minx AZ x − c Z 2 such that kx ∗ k2 is minimized. Then,
AZ (s ∗ − x ∗ )

def

Proof. Let γ =

a k22
ka
2
a k22
w +ka

2

≤

kaa k2
σmax (AZ ) c Z
w

ΠAZ c Z

2

2

.

. By the definition of P in Lemma 7.11,
(1 − γ) I 4 P 4 I .

Thus,
(1 − γ) AZ s ∗ − c Z

2
2

≤ AZ s ∗ − c Z

2
P

2
P

≤ AZ x ∗ − c Z

≤ AZ x ∗ − c Z

2
.
2

(37)

The second inequality is due to s ∗ ∈ arg minx AZ x − c Z P . Note that AZ x ∗ − c Z is orthogonal
to the column space of AZ , in particular, it is orthogonal to AZ (s ∗ − x ∗ ).
AZ (s ∗ − x ∗ )

2
2

2
2

≤

= AZ s ∗ − c Z

2
2

− AZ x ∗ − c Z

2
.
2

By Equation (37),
AZ (s ∗ − x ∗ )

γ
AZ x ∗ − c Z
1−γ

2
2

≤

γ
cZ
1−γ

2
.
2

If ΠAZ c Z = 0, then s ∗ = x ∗ = 0 and the claim holds. Otherwise, by Lemma 4.8,
r
γ
Z ∗
∗
A (s − x ) 2 ≤
σmax (AZ ) c Z 2 ΠAZ c Z 2 .
1−γ
Plugging the value of γ completes the proof.
We then show the approximate solvers of the two linear systems are close.
Lemma 7.13. Let 2 be the error parameter returned by a call to Reduce Gz toGz,2 (AZ , c Z , 1 )
(Algorithm 10). Let x Z,2 be a vector such that AZ,2 x Z,2 − ΠAZ,2 c Z,2 2 ≤ 2 ΠAZ,2 c Z,2 2 . Let x
be the vector returned by Algorithm 11. Then,
AZ x − ΠAZ c Z

2

48

≤ 1 ΠAZ c Z

2

.

Proof. If ΠAZ c Z = 0, then x = 0 and the statement holds. If ΠAZ c Z 6= 0, then x = x Z,2
1:nA +1 . Let
Z
∗
Z
x ∈ arg minx A x − c 2 . We now bound the difference between our solution x and x ∗ . Let
x Z,2∗ ∈ arg minx AZ,2 x − c Z,2 2 of the form (s ∗ ; α; −α). By the triangle inequality,
AZ (x − x ∗ )

2

≤ AZ (x − s ∗ )

2

+ AZ (s ∗ − x ∗ )

2

.

(38)

The second term is upper bound by Lemma 7.12. It remains to upper bound the first term. Let
def

δ = 2 ΠAZ,2 c Z,2

2

.

(39)

Without loss of generality, we write x Z,2 as (x ; β; −β), where x is the output of Algorithm 11.
AZ,2 x Z,2 − ΠAZ,2 c Z,2

2
2

= AZ,2 x Z,2 − x Z,2∗



2
2


x − s∗
A a −aa
 β−α 
=
0 w −w
α−β

 Z
A (x − s ∗ ) + 2 (β − α) a
=
2w (α − β)
Z







2

2
2
2

= kA (x − s ∗ ) + 2 (β − α) a k22 + 4w2 (α − β)2 .
Since AZ,2 x Z,2 − ΠAZ,2 c Z,2

2

≤ δ, we have
4w2 (α − β)2 ≤ δ 2 ,

that is,
(α − β)2 ≤

δ2
.
4w2

(40)

Similarly, we have
AZ (x − s ∗ ) + 2 (β − α) a

2

≤ δ.

By the triangle inequality,
AZ (x − s ∗ )

2

− 2 k(β − α) a k2 ≤ δ.

Plugging (40) and (39) into the above inequality, and rearranging it,


kaa k2
Z
∗
2 c Z,2 2 .
A (x − s ) 2 ≤ 1 +
w
By Lemma 4.8,
Z

∗

A (x − s )

2



kaa k2
≤ 1+
2 σmax (AZ ) c Z,2
w

2

ΠAZ c Z

Together with Lemma 7.12 and Equation (38), we have


 
kaa k2
kaa k2
Z
∗
A (x − x ) 2 ≤
+ 1+
2 σmax (AZ ) c Z
w
w
49

2

2

.

ΠAZ c Z

2

.

According to our setting of w and a in line 1 and line 5 of Algorithm 10,
w≥

kaa k2
.
2

This implies that
AZ x − ΠAZ c Z

2

≤ 1 ΠAZ c Z

2

,

which completes the proof.
We then bound the singular values of AZ,2 .


Z
Z > Z
Claim 7.14. λmax (AZ,2 )> AZ,2 ≤ O −2
1 mλmax ((A ) A ) A

2
∞

cZ

2
2



.

Proof. Let λ1 bet he largest eigenvalue of (AZ,2 )> AZ,2 , and y = (ỹ ; α; −α) be the associated
eigenvector of unit length. Let µ1 be the largest eigenvalue of (AZ )> AZ .

λ1 =

AZ a −aa
0 w −w

2


y

2

 2
0
0
0
=
y
y+
0 w − 1 −w + 1
2


>
2
>
Z > Z
> Z
2
= ỹ (A ) A ỹ + 4αaa A ỹ + 4α a a + 1 + 8α (w − 1) + 4α2 (w − 1)2 .


Z

A
0

a −aa
1 −1





By the Courant-Fischer Theorem,
ỹ > (AZ )> AZ ỹ ≤ µ1 ỹ > ỹ ≤ µ1 .
By Cauchy-Schwarz inequality,
a > AZ ỹ ≤ kaa k2 AZ ỹ

2

≤ kAk∞

√ √
m µ1 .

Thus,
λ1 ≤ µ1 + 4 |α| AZ
√
≤
µ1 + 2 A Z
≤ 2µ1 + 8m AZ


√ √
2
m
µ
+
4α
AZ
1
∞
√ 2
m + 12α2 w2
∞
2
∞

2
m
∞

+ 12α2 w2 .

Our setting of w in line 1 of Algorithm 10 gives


λmax (AZ,2 )> AZ,2 ≤ 2µ1 + 8m AZ

2
∞

Z
+ 108−2
1 mµ1 A

This completes the proof.
Claim 7.15. λmin ((AZ,2 )> AZ,2 ) ≥


+ 1 + 12α2 w2

Z > Z
2λmin ((A

 ) A )
.
Z 2
2 kA k∞ m+1 +λmin ((AZ )> AZ )

Proof. Let

C =

AZ a −aa
0 1 −1
50


.

2
∞

cZ

2
.
2

Note
A

Z,2


=C +

0
0
0
0 w − 1 −w + 1


.

We can check that


(AZ,2 )> AZ,2


0
0
0
w2 − 1
−(w2 − 1)  .
= C >C +  0
2
0 −(w − 1)
w2 − 1

By our setting, w2 − 1 ≥ 0. The second matrix is a rank-one PSD matrix. Since AZ,2 and C has
same null space,




λmin (AZ,2 )> AZ,2 ≥ λmin C > C .
It suffices to lower bound the smallest nonzero eigenvalue of C > C .

def
Let λk = λmin C > C , and y = (ỹ ; α; −α) be the associated eigenvector of unit length. Let

def
µk = λmin (AZ )> AZ , and x be the associated eigenvector of unit length.


λk = ỹ > (AZ )> AZ ỹ + 4αaa > AZ ỹ + 4α2 a >a + 1
= AZ ỹ + 2αaa

2
2

(41)

+ 4α2

≥ 4α2 .
On the other hand,
2
 p
ỹ > (AZ )>a a > AZ ỹ
a > AZ ỹ
>
+ ỹ > (AZ )> AZ ỹ −
λk = 2α a a + 1 + √
a >a + 1
a >a + 1


aa >
≥ ỹ > (AZ )> I − >
AZ ỹ .
a a +1
Take eigen-decomposition of the matrix in the middle,
λk ≥ ỹ > (AZ )> QDQ > AZ ỹ ,


a >a
where D := diag 1 − a >
,
1,
.
.
.
,
1
. By Claim 7.10, y ⊥ null(AZ,2 ) implies ỹ ⊥ null(AZ ). By
a +1
the Courant-Fischer Theorem,
λk ≥

1
a >a

+1

AZ ỹ

2
2

≥

1
a >a

+1

µk kỹ k22 .

Together with (41),



µk
2
λk ≥ max 4α , >
kỹ k2 .
a a +1
2

Since ky k22 = ỹ 22 + 2α2 = 1,
λk ≥

2µk
2µk

≥ 
.
2
2 (aa >a + 1) + µk
2 kAk∞ m + 1 + µk

This completes the proof.
51

The above two lemmas give the following bound on the condition number of AZ,2 .

√
√
Lemma 7.16. κ(AZ,2 ) = O −1
m AZ ∞ c Z 2 σmax (AZ ) + κ(AZ ) m AZ ∞ .
1
Proof of Lemma 3.4. By Lemma 7.2 and 7.9, we have
nnz(AZ,2 ) = O(s).
The smallest nonzero entry does not change. The largest entry of AZ,2 is at most

√
max{O(kAk∞ ), w} = O −1 n2 mσmax (A) kAk∞ c A 2 .
By Claim B.1 and B.2, the largest entry is upper bounded by


O −1 s9/2 U 3 .
By Lemma 7.8 and 7.16, and Claim B.1 and B.2

κ(AZ,2 ) = O −1 s8 U 3 K .
By Lemma 7.4 and 7.13, we have


5/2 2
−1 .
−1
=
O
s
U
2
This completes the proof.

8

2D Trusses

In this section, we show that the matrix B constructed by Reduce Gz,2 toMC2 (Algorithm 1), is a
2D Truss Incidence Matrix defined in Definition 1.4. It follows that for any function f , G ≤f MC2
implies G ≤f T2 . We assume that the algorithm Reduce Gz,2 toMC2 is called on a matrix A
with no two identical rows: if there are identical rows, these rows can be collapsed into one row
without changing the associated normal equations, by reweighting the resulting row, similar to the
technique used in the proof of Claim 9.4. Details are left to the reader. A key step is to show
that a 2-commodity gadget in the reduction corresponds to a 2D truss subgraph, which we call the
2D-truss gadget.
Without loss of generality, we let u-variables correspond to the horizonal axis and v -variables
to the vertical axis of the 2D plane. According to Definition 1.2 and 1.4:
1. an equation u i − u j = 0 in a 2-commodity linear system corresponds to a horizontal edge in
the 2D plane;
2. an equation v i − v j = 0 in a 2-commodity linear system corresponds to a vertical edge in the
2D plane;
3. an equation u i − v i − (u j − v j ) = 0 in a 2-commodity linear system corresponds to a diagonal
edge in the 2D plane.
Note that our reduction here heavily relies on the ability to choose arbitrary weights. In particular,
the weights on the elements are not related at all with the distances between the corresponding
vertices.

52

Our strategy for picking the coordinates of the vertices of the constructed 2D truss is the
following: we first pick the coordinates of the original n vertices randomly, and then determine the
coordinates of the new vertices constructed in the reduction to satisfy all the truss equations.
For the n original vertices, we pick their u-coordinates arbitrarily and pick their v -coordinates
randomly. We pick an n-dimensional random vector y uniformly distributed on the n-dimensional
sphere centered at the origin and with radius R = kAk1 n10 . We then round each entry of y
to have precision δ = 10−10 , so that the total number of bits used to store an entry is at most
O(log(n kAk1 )). Let ỹ be the vector after rounding. We assign the v -coordinate of the ith vertex
to be the ith entry of ỹ .
We then pick the coordinates of the new vertices in the order they are created. Note that
each time we replace two vertices in the current equations, say s j1 , s j2 , whose coordinates have
already been determined, we create a 2D truss gadget with 7 new vertices, say s t , s t+1 , . . . , s t+6
(See Algorithm 2 MC2 Gadget for the construction.). According to the construction of this gadget,
the new vertices s t+1 , . . . , s t+6 only appear in this single gadget, whose coordinates do not affect
other vertices. Figure 2 is the corresponding subgraph which satisfies all the equations in the 2D
truss gadget. Note the two triangles (s t+3 , s t+5 , s t+6 ) and (s t+3 , s t+4 , s t+5 ) need to be isosceles
right triangles, which implies v t = (v j1 + v j2 )/2. Note also that we can assign u-coordinates
to the new vertices which are not between the u-coordinates of s j1 and s j2 . In fact, using an
appropriate choice of u-coordinates and edge weights, we can always place s t , s t+1 , . . . , s t+6 to get
the desired equations, provided v j1 6= v j2 , which we later will argue holds with high probability
using Lemma 8.1. First, however, we state and prove Lemma 8.1.
s j2

st

s t+6

s t+2

s t+3

s t+5

s t+1

s t+4

s j1

Figure 2: Geometric realization of the mutlicommodity flow gadget generated by as a truss matrix
Lemma 8.1. Let a ∈ Rn be a fixed vector such that −2 ≤ a i ≤ 2, ∀i ∈ [n] and a > 1 = 0, and a 6= 0.
Let ỹ be a vector picked as above. Then,


2δn2
.
Pr a > ỹ = 0 ≤
kaa k2 R
def

Proof. Let ∆ = ỹ − y . Clearly, −δ ≤ ∆i ≤ δ, ∀i ∈ [n].
X
a > ỹ = a > (y + ∆) ≤ a > y +
|aa i | |∆i | ≤ a > y + 2δn.
i∈[n]

53

Similarly, a > ỹ ≥ a > y − 2δn. Thus,




Pr a > ỹ = 0 ≤ Pr a > y ≤ 2δn .
Since the distribution of y is rotation invariant, we assume without loss of generality a =
(kaa k2 , 0, . . . , 0). Let A be the area of the n-dimensional sphere.
A=

2π (n+1)/2
 · Rn .
Γ n+1
2

Let Aδ be the area of {y : ky k2 = R, |y 1 | ≤ 2δn/ kaa k2 }. Then,
Aδ ≤

2π n/2
2δn
n−1

·
.
n ·R
kaa k2
Γ 2

Thus, (assume n is even, the case of odd n can be checked similarly)
Pr



 A
δ
a >y ≤ δ =
A
Γ( n+1
2δn
2 )
√
·
n
kaa k2 πR Γ( 2 )
√
2δn
n! π
=
·
kaa k2 R n!!(n − 2)!!
2δn2
.
≤
kaa k2 R
≤

This completes the proof.
By our construction of the truss, for each vertex, its v -coordinate can be written as a fixed
convex combination of ỹ , say c > ỹ in which c > 1 = 1 and c i ≥ 0, ∀i ∈ [n]. Note that when
algorithm Reduce Gz,2 toMC2 is applied to a matrix A, it processes the jth row in k iterations
def

where by Lemma 4.1 we have k ≤ log kAj k1 . Let p = log kAj k1 . Given how the convex combination
specified by c is formed, it follows that 2p c is an integer vector.
Next we argue that when two variables are chosen for pairing by the Reduce Gz,2 toMC2
algorithm as it processes some row Aj , these two variables will have their v -coordinates represented
as convex combinations c > ỹ and d > ỹ where crucially c 6= d . This ensures that 2p (c − d ) is a
non-zero integer vector, and hence kc − d k2 ≥ 2−p = 1/ kAj k1 . This follows from stronger claim
stated below, which we prove later.
Claim 8.2. Suppose algorithm Reduce Gz,2 toMC2 is processing some row Aj . Let V l be the set
of variables with non-zero coefficients in the main equation Aj at the lth iteration of the while-loop
in Line 20 of Reduce Gz,2 toMC2 . Let S l be the set of associated vertices. Consider two arbitrary
vertices s, t ∈ S l , and let c, d ∈ Rn be the non-negative vectors with c > 1 = d > 1 = 1 s.t. the
v -coordinates of s and t represented as convex combinations are c > ỹ and d > ỹ respectively. For
every i ∈ [n], we view entries c i and d i as fixed point binary numbers: c i = α0i .α1i α2i . . . αpi and
P
P
d i = β0i .β1i β2i . . . βpi , or equivalently c i = pk=0 αki 2−k and d i = pk=0 βki 2−k . Then there is no
index k s.t. 1 = αki = βki , i.e. there is no index where the kth bit is 1 in both strings.

54

def

These two vertices have same v -coordinate if and only if (c − d )> ỹ = 0. Let a = c − d . Then,
−2 ≤ a i ≤ 2, ∀i ∈ [n], a > 1 = 0, and
kaa k2 ≥ 1/ kAj k1 .
By Lemma 8.1,
 2δn2 kA k

j 1
Pr c > ỹ = d > ỹ ≤
.
R
By Lemma 4.1, the total number of the vertices in the truss is at most

O n2 log n .
By a union bound, the probability that there exist two different vertices with same v -coordinate is
at most
 2 
2
2δn2 kAk1
log n
.
· O n2 log n = O
R
n4
Proof of Lemma 3.8. Since the linear system for 2D trusses is the same as the linear system for
2-commodity, all complexity parameters of these two linear systems are the same.
Proof of Claim 8.2. We prove the claim by induction. Our induction hypothesis is simply that the
claim holds in round l.
Note that by Lemma 4.1 every time a new vertex is created (during the lth iteration of whileloop in Line 20), it is always paired in the following iteration (iteration l + 1) and then disappears
(i.e. has zero coefficient) in the main equation in all following iterations.
Observe that the convex combination vector a for each vertex that corresponds to an original
variable i is has a i = 1 and a h = 0 for all h 6= i. This proves the induction hypothesis for the base
case of the variables in the main equation Aj before the first iteration of the while-loop (i.e. l = 0).
Suppose c,d are the convex combination vectors for two variables that exist in some round l.
Assume the induction hypothesis holds for round l − 1. Write the binary strings
ith of
Ppfor the
−k
i
i
i
i
i
i
i
i
i
both vectors as c i = α0 .α1 α2 . . . αp and d i = β0 .β1 β2 . . . βp , or equivalently c i = k=0 αk 2 and
P
d i = pk=0 βki 2−k . Suppose for the sake of contradiction that there exists some k s.t. αki = βki = 1.
Trivially, it cannot be the case that the such a collision occurs if either variable is an original
variable. So both variables must be new variables. The bit string for entry c i is created by
averaging two bit strings of variables from the main equation in round l − 1, say η0i .η1i η2i . . . ηpi and
γ0i .γ1i γ2i . . . γpi . Similarly, bit string for entry d i is created by averaging two bit strings of variables
from the main equation in round l − 1, say θ0i .θ1i θ2i . . . θpi and σ0i .σ1i σ2i . . . σpi . Note that each variable
can only be paired once in each iteration, so the four bit strings must come from distinct variables
in round l.
αki = 1 requires that exactly one of the following conditions is true:
i
1. γk−1
=1
i
2. ηk−1
=1
i
i
3. A “carry” occurred when adding strings γki γk+1
. . . γpi and ηki ηk+1
. . . ηpi .

But, Case 3 immediately leads to a contradiction, as a carry can only occur when there exists some
bit position g s.t. γgi = ηgi . But this is false by the induction hypothesis. Thus we must be in
either Case 1 or Case 2. By similar logic, we can conclude from βki = 1 that exactly one of the
i
i
following must be true: either θk−1
= 1 or σk−1
= 1. All together, we have concluded that exactly
i
i
i
i
two of the four bits ηk−1 , γk−1 , θk−1 , and σk−1 must be set to 1. This contradicts the induction
hypothesis. Having established a contradiction whenever the induction hypothesis fails at step l,
we have shown that it holds at this step.
55

9

Connections with Interior Point Methods

In this section, we discuss how applications in scientific computing and combinatorial optimization
produce the linear systems that we show are hard to solve. We first give a brief overview of interior
point methods, with focus on how they generate linear systems in Section 9.1. Then we formalize
the matrices that interior point methods produce when run on 2-commodity flow matrices 9.2 and
isotropic total variation minimization 9.3.

9.1

Brief Overview of Interior Point Methods

Interior point methods [Wri97, Ye11, Nem, BV04, DS08, LS14, LS15] can be viewed as ways of
solving convex optimization problems via a sequence of linear systems. For simplicity, we choose
the log-barrier based interpretation from Chapter 11 of the book by Boyd and Vandenberghe as our
starting point. The main idea is to represent a convex optimization problem with constraints as a
sequence of linear programs with terms called barrier functions added to the objective. Specifically,
turning the problem:
min

f (y )

y

subject to:

My = b
y ≥0

into the equality-constrained optimization problem:
min
y

t · f (y ) −

X

log (y i )

i

subject to:

My = b

for a parameter t that is gradually increased (by factors of about 1 + n−1/2 throughout the course
of the algorithm). The solution of the above optimization problem converges to an optimal solution
of the original linear programming, as t goes to infinity. Between these increase steps, the algorithm performs Newton steps on this log barrier objective, which when combined with the equality
constraint M y = b requires solving the problem:
min
∆y

subject to:

− g (y )> ∆y
M ∆y = 0
k∆y kH (y ) ≤ 0.1

which is to maximize the projection along the gradient subject to the second order term being at
most 0.1 and staying in the null space. This can in turn be interpreted as a least squares problem,
and solving the linear system:
M H (y )−1 M > x = −tM H (y )−1 g (y ).

9.2

(42)

2-Commodity Flow

We now show that solving 2-commodity flow problems using interior point methods as described
in Subsection 9.1 can lead to any system in the class MC2>0 . There are also many variants of the
multicommodity flow problem [Mad10], and we work with the minimum cost version due to it being
the most general.
56

Definition 9.1 (Min-cost 2-commodity flow problem). Given a directed graph G = (V, E) with n
vertices and m edges, a positive edge-capacity vector z ∈ Rm , a positive edge-cost vector c ∈ R2m ,
and two vertex-demand vectors d 1 , d 2 ∈ Rn . The goal is to compute two flows y 1 , y 2 such that,
1. y 1 satisfies the demand d 1 , and y 2 satisfies the demand d 2 ,
2. for each edge e ∈ E, the sum of the two flows on e is no larger than the edge capacity z e , and
3. the total cost of the two flows is minimized.

To formulate this as a linear program, we let y = y 1 ; y 2 be the two flows, and N be the
edge-vertex incidence matrix of graph G.
min c > y
s.t. N > y 1 = d 1
N >y 2 = d 2
y 1, y 2, z − y 1 − y 2 ≥ 0
Write this linear programming as the following minimization problem with a logarithmic barrier
function,
min c > y −

1 X
log y 1i + log y 2i + log(z i − y 1i − y 2i )
t
i∈[m]

s.t.

> 1

N y = d1
N >y 2 = d 2

where t > 0 is a parameter.
Definition 9.2. We say a linear system B > Bx = B > c is a Minimum Cost 2-commodity Flow
IPM Linear System if it can be obtained from as an instance of the Newton-Step Linear System of
Equation (42) for some Min-cost 2-commodity flow problem.
The class of Minimum Cost 2-commodity Flow IPM Linear Systems is appears more restrictive
than MC2>0 , but as we will see, it is essentially equivalent, and still sufficiently expressive that any
linear system can be reduced to it. The main result that we will sketch in this section is:
Lemma 9.3. For any linear system Ax = c (with error parameter ) with polynomially bounded
sparse parameter complexity, there exists an efficient reduction to a Minimum Cost 2-commodity
Flow IPM Linear System.
Our reduction from a general linear system (A, c, ), where A ∈ G, to a Minimum Cost 2commodity Flow IPM Linear System is the same as our reduction to MC2>0 , except that before
our chain of reductions, we first multiply the linear system by a diagonal matrix S with diagonal
entries that are ±1. We will later specify how to choose these signs. This gives us the linear system
S Ax = S c with the same error parameter as before. It is easy to verify that this system has the
same sparse parameter complexity as the original linear system, and an approximate solution to
this system is an approximate solution to the original system with the same .
We now apply our usual chain of reductions to get a linear system in over a matrix in MC2>0 .
We write this system as
>0
(B >0 )> B >0 x = (B >0 )> c B .
57

The remainder of this Section is dedicated to showing that this linear system is a Minimum Cost
2-commodity Flow IPM Linear System, thus proving Lemma 9.3.
b > W B,
b where B
b is the
We can pull out the edge weights from B >0 by writing (B >0 )> B >0 = B
unweighted 2-commodity edge-vertex incidence matrix with the same edge structure as B >0 , and
W is the diagonal matrix of the edge weights. Then, the linear system in MC2>0 can be written as
b > W Bx
b =B
b > W 1/2 c B>0 .
B

(43)

Before we prove Lemma 9.3, we explore some properties of the above linear system.
e a diagonal matrix W
f , and
Claim 9.4. There exist a 2-commodity edge-vertex incidence matrix B,
a vector c̃ such that
e >W
fB
e =B
b >W B
b and B
b > W 1/2 c B>0 = B
e > W 1/2 c̃,
1. B
e are distinct.
2. all the rows of B
b are distinct, then we set B
e = B,
b W
f = W and c̃ = c B>0 . Otherwise,
Proof. If all rows of B
b are the same. Note
assume the ith row and the jth row of B
X
b >W B
b =
b >B
b k + (W ii + W jj )B
b >B
b i,
B
W kk B
k
i
k6=i,j

and

b > W 1/2 c B>0 =
B

X

>0
1/2
b > + (W 1/2 c B>0 + W 1/2 c B>0 )B
b >.
W kk c B
B
k
i
j
i
k
ii
jj

k6=i,j

e and W
f by removing the jth row of B,
b W and set the corresponding edge
We can construct B
>0
f
weight W ii to be W ii + W jj . We construct c̃ by removing the jth entry of c B and set c̃ i to be
>0
>0
1/2
1/2
(W ii c B
+ W jj c B
)/(W ii + W jj )1/2 . We repeat the above process to merge the identical
i
j
rows until we get the desired matrices and vector.
>0

b − cB
Remark. Let x ∗ be a minimizer of W 1/2 Bx

2

. We know that x ∗ is a solution of the

b > W Bx
b = B
b > W 1/2 c B>0 . By Claim 9.4, x ∗ is a solution of B
e >W
fB
e >x =
normal equations B
e >W
f 1/2 c̃, which in turn gives that x ∗ ∈ arg minx W
f 1/2 Bx
e − c̃ . Let x be an approximate
B
2

solution such that
kx − x ∗ k b >

b
B WB

b > W 1/2 c B>0
≤ B

b > W B)
b †
(B

.

By Claim 9.4, the above equation is equivalent to
e >W
f 1/2 c̃
kx − x ∗ k e > f e ≤  B
B WB

e >W
f B)
e †
(B

f 1/2 Bx
e − c̃
It means that x is an approximate solution of minx W

.

with -accuracy, vice versa.
2

Thus, the two minimization problems before and after merging identical rows are equivalent.
b are distinct. We define the
Without the loss of generality, we assume that all the rows of B
b
underlying graph of B to be the simple graph over n vertices where vertices i and j are connected
b has at least one type of edges between i and j. According to the constructions
if and only if B
58

in Reduce Gz,2 toMC2 (Algorithm 1) and ReduceMC2 ToMC2>0 (Algorithm 4), each edge of the
b (that is, type 1, type 2 and type 1 + 2), and
underlying graph has exactly 3 types of edges in B
>0
the edge weights of B
are either at least 1 or equal to a tiny number δ assigned in line 1 of
Algorithm 4. We call edge weights which are at least 1 as large weights.
b exactly one of its type 1, type 2, and type
Claim 9.5. For each edge in the underlying graph of B,
1 + 2 edges has large weight.
Proof. Note that the large-weight edges all appear in B ∈ MC2 constructed in Algorithm 1. An
edge in B is either a gadget edge in B or a non-gadget edge in A.
1. By Algorithm 2, all edges in a single MC2 -gadget are distinct.
2. We show that each gadget edge is distinct from all other edges. By Algorithm 2, in a 2commodity gadget, each “old” vertex is connected to a “new” vertex, and all new vertices
except x t are independent of the new vertices created in all other gadgets. However, in each
gadget including x t (as a new vertex or an odd vertex), x t is connected to new vertices which
are created in that gadget and do not appear in any other edge outside that gadget.
3. It is possible that two identical non-gadget edges exist when the original linear system instance
has redundant or inconsistent constraints. But note all non-gadget edges are type 1 edges, it
b By Claim 9.4, we can
implies that two identical edges correspond to two identical rows in B.
always merge such identical rows so that all non-gadget edges are distinct.
This completes the proof.
b W and the vector c B>0 are given by the reductions
Proof of Lemma 9.3. Note the matrices B,
from the general linear system instance Ax = c A . The goal of proving Lemma 9.3 is to determine
the edge-vertex incidence matrix N , the cost vector c, the demand vectors d 1 , d 2 , the edge capacity
vector z , and the flows y 1 , y 2 such that the intermediate linear systems that arise in an interior
point method (42) for solving the Min-cost 2-commodity Flow Problem is exactly the above linear
system (43).
We first make the coefficient matrices of the two linear systems in (42) and (43) to be the same,
by choosing the matrix N and the vectors y 1 , y 2 , d 1 , d 2 , z properly. The coefficient matrix of the
equality constraints of the 2-commodity linear programming is


N>
M =
.
(44)
N>
b By denoting
We choose N to be the edge-vertex incidence matrix of the underlying graph of B.
the residue flow amount along an edge as
def

y r = z − y 1 − y 2,
we can rewrite the Hessian as:
1 X 1
1
1
H =
e 2i−1 e >
e 2i e >
(e 2i−1 + e 2i ) (e 2i−1 + e 2i )> ,
2i−1 +
2i +
r
1
2
2
2
2
t
(y
)
(y i )
(y i )
i
i∈[m]

where e i ∈ R2m is the standard basis vector.
We can rearrange the rows and columns of H such that, row 2i−1 and column 2i−1 correspond
to the y 1 -flow of the ith edge, and row 2i and column 2i correspond to the y 2 -flow of the ith edge.
59

After this rearrangement, H becomes a block diagonal matrix, where each edge i corresponds to a
one 2 × 2 block given by
!
1
1
1
+
r
r
1
2
2
2
def
(y i )
(y i )
(y i )
.
(45)
H 2i−1:2i =
1
1
+ (y1r )2
r
2
2
2
(y )
(y )
i

i

i

The block diagonal structure of H then gives:

−1
H −1 = diag H −1
1:2 , . . . , H 2m−1:2m ,
and for each edge i we have:




1
(y 1i )2 (y ri )2
0
1 −1
−1
1 2
2 2
H 2i−1:2i =
,
+ (y i ) (y i )
0
(y 2i )2 (y ri )2
−1 1
αi

(46)

where αi = (y 1i )2 + (y 2i )2 + (y ri )2 . It means that
1. the type 1 edge has weight (y 1i )2 (y ri )2 /αi ,
2. the type 2 edge has weight (y 2i )2 (y ri )2 /αi , and
3. the type 1 + 2 edge has weight (y 1i )2 (y 2i )2 /αi .
Note these three types of weights are symmetric. Let w(i,1) , w(i,2) , w(i,1+2) be the edge weights
of type 1, type 2 and type 1 + 2 edge of the ith edge, respectively. These edge weights are given by
the corresponding diagonals of W . We set the values of y 1i , y 2i , y ri such that
(y 1i )2 (y ri )2 /αi = w(i,1) ,
(y 2i )2 (y ri )2 /αi = w(i,2) ,
(y 1i )2 (y 2i )2 /αi = w(i,1+2) .
Solving the above equations gives:
w(i,1) w(i,1+2)
+ w(i,1) + w(i,1+2) ,
w(i,2)
w(i,2) w(i,1+2)
(y 2i )2 =
+ w(i,2) + w(i,1+2) ,
w(i,1)
w(i,1) w(i,2)
(y ri )2 =
+ w(i,1) + w(i,2) .
w(i,1+2)
(y 1i )2 =

(47)

The y 1 , y 2 and y r determine the edge-capacity vector z , and they together with N determine the
vertex-demand vectors d 1 and d 2 .
Since we rearranged the rows and columns of H , we need do the same rearrangement for M .
After rearranging row and columns of M , M is of the form that, column 2i − 1 and column 2i
correspond to the y 1 -flow and y 2 -flow of the ith edge, respectively, and row 2i − 1 and row 2i
correspond to the u-coordinate and v -coordinate of vertex i, respectively. In M , each edge has a
2 × 2 identity matrix for one endpoint, and a negative 2 × 2 identity matrix for the other endpoint.
By the above setting, we can check that
>

b W B.
b
M H −1 M > = B
60

(48)

b has dimension 3m × 2n.
Note that M has dimension 2n × 2m, and B
Secondly, we make the right hand side vectors of the two linear systems (42) and (43) to be the
same, that is,
b > W 1/2 c B>0 ,
−tM H −1 g = B

(49)

by choosing the cost vector c properly. The gradient g is given by:
g =c−

1
1
1 X 1
e 2i−1 + 2 e 2i − r (e 2i−1 + e 2i ) .
1
t
yi
yi
yi
i∈[m]

For simplicity, let
def

f (y ) = t (c − g ) .
f (y ) is a vector in 2m-dimension such that, the ith edge corresponds to the 2 × 1 sub-vector
!
1
1
−
r
1
y
yi
i
.
1
1
−
2
yr
y
i

i

Then we have
−tM H −1 g = M H −1 (f − tc) .
b > W 1/2 c B>0 . According to our
Recall the right hand side vector of Equation (49) is B
construction in Reduce Gz,2 toMC2 (Algorithm 1) and ReduceMC2 ToMC2>0 (Algorithm 4),
>0
c B = (c A ; 0), where c A is the right hand side vector of the general linear system instance.
>0
Note all the nonzero entries of c B correspond to the type 1 edges in the main constraint set A,
see line 16 of Algorithm 1.
Recall the structure of M in Equation (44), the edge-vertex incidence matrix N contains all
the m edges of the underlying graphs, which are exactly the set of edges in B constructed in
Reduce Gz,2 toMC2 (Algorithm 1). Note for each edge of N , its large-weight edge can be any
of the corresponding type 1, type 2 and type 1 + 2 edges. By Claim 9.5, we can partition the
underlying graph edges into E1 ∪ E2 ∪ E1+2 such that Ek contains all the underlying edges whose
>0
large-weight edges are the corresponding type k edges, k ∈ {1, 2, 1 + 2}. Note c B
(i,1) 6= 0 only if edge
i is in E1 .
We define a vector c 0 ∈ R2m such that the entry of c 0 which corresponds to the y 1 -flow of the
>0
1/2
ith edge equals to w(i,1) c B
(i,1) . This definition gives:
b > W 1/2 c B>0 = M c 0 .
B

(50)

By Equation (49),
M H −1 (f − tc) = M c 0 .
Rearranging it,

M H −1 f − tc − H c 0 = 0.
According to MC2 Gadget (Algorithm 2), each gadget has a cycle containing both the 2 edges
in E1+2 : (t + 3, t + 4, t + 2, t + 5, t + 6, t + 1, t + 3). Let E 0 be the set of all the edges in such
gadget cycles. Let p ∈ {0, 1}m such that each entry corresponding to an edge in E 0 has value 1.
def

Let q = (p; p). According to Reduce Gz,2 toMC2 (Algorithm 1), every gadget edge i has c 0i = 0.
Thus, we are free to choose the sign of a column of N > which corresponds to a gadget edge, which
61

7

does not change Equations (48) and (50).
sign flipping, M q = 0. It gives that

Without loss of generality, we assume that after proper


M H −1 f − tc − H c 0 + H q = 0.
We set tc = f − H c 0 + λH q , where λ > 0 is a sufficiently large constant to be determined later.
For each edge i, let Ii be an indicator whose value is 1 if edge i in the set E 0 and 0 otherwise.
For the ith edge, its corresponding 2 × 1 sub-vector in tc is
!
!


1
1
1
1
−
+
r
r
1
1
2
2
ci
y
(y i )
yi
(y i )
i
t
=
− c 0i
+ Ii λs i ,
1
1
1
−
c m+i
r
r
2
2
y
(y )
y
i

i

i

where

1
(y 1i )2
1
(y 2i )2

i def

s =

+
+

2
(y ri )2
2
(y ri )2

!
.

Note that s i ≥ 0 (entry-wise). Plugging the solution of y from Equation (47) gives:
q

 q
!
w(i,2)
w(i,1+2)


w(i,2) +w(i,1+2)
−
ci
γ
γ
0
γ
 − ci
q
t
+ Ii λs i ,
=  qw
w(i,1+2)
w(i,1+2)
(i,1)
c m+i
−
γ
γ

γ

def

where γ = w(i,1) w(i,1+2) + w(i,1) w(i,2) + w(i,2) w(i,1+2) .
By Claim 9.5, each edge of the underlying graph has exactly one type of edges with large weight.
We deal with each of these cases separately under the condition that two of the edge’s weights are
much smaller than that of the third one. We will also denote this ratio using  > 0:
1. If w(i,2) = w(i,1+2) = w(i,1) , then

t

ci



c m+i

≥ −c 0i

2
µ

µ

!
+

0√

1− 
√
µ

!
,

where µ = (2+2 )w(i,1) . This corresponds to an edge of E1 . If c 0i ≤ 0, then tc i , tc m+i ≥ 0. If
1/2

>0

>0

B
A
c 0i > 0, then c 0i = w(i,1) c B
(i,1) where c (i,1) comes from an entry of c . We multiply -1 on both
sides of the corresponding equation in the general linear system instance Ax = c A which
creates this edge i, and will get c 0i < 0 instead. 8

2. If w(i,1) = w(i,1+2) = w(i,2) , then

t

ci
c m+i


≥

−c 0i

1+
µ

µ

!
+

√
1− 
√
µ

0

!
,

where µ = (2 + 2 )w(i,2) . This corresponds to an edge of E2 , in which c 0i = 0. We can check
that tc i , tc m+i ≥ 0.
7

Note the demand vectors d 1 and d 2 change after we change N .
We first flip signs of equations of Ax = c A so that c A ≤ 0 (entry-wise), which determines N , M and c 0 . We
then flip signs of columns of N , which we guarantee does not violate Equation (48) and (50).
8

62

3. If w(i,1) = w(i,2) = w(i,1+2) , then

t

ci
c m+i



1+
µ
1
µ

= −c 0i

!
+

√
−1
√
√ µ
−1
√
µ

!
+ Ii λs i ,

where µ = (2 + 2 )w(i,1+2) . This corresponds to an edge of E1+2 , in which c 0i = 0 and Ii = 1.
We choose λ such that, for every edge i in this case, tc i , tc m+i ≥ 0.

Note in IPMs, the values of t is a sequence of increasing values. Whenever t ≥ m/0 , the
optimality gap is smaller than 0 . This ensures that the cost vector c’s entries are in a reasonable
range.

9.3

Isotropic Total Variation Minimization

For the isotropic total variation problem, we follow the formulations given in [CMMP13a], namely
given a graph G = (V, E, w ), we partition the sets into S1 . . . Sk , and minimize the objective
s X
X
2
ky − sk2 +
wi
w uv (y u − y v )2 .
1≤i≤k

(u,v)∈Si

where s is the input signal vector. Let N be the edge-vertex incidence matrix of the graph G, then
the dual of the grouped flow problem is:
max
subject to:

s >N >f
X
−1
2
w −1
e f e ≤ wi

∀1 ≤ i ≤ k

e∈Si

Both the primal and dual problems can be formulated into log barriers via second order
cones [GY04]. In the case with
P minimizing vertex labels, we introduce a variable y i for each
cluster i, and instead minimize i w i y i subject to the constraint
X
2
y 2i ≥
w −1
e (y u − y v ) ,
(u,v)∈Si

which in turn leads to the logarithmic barrier function


X
2
log y 2i −
w −1
,
e (y u − y v )
(u,v)∈Si

while in the flow case the barrier functions for each set of edges Si is


X
2
log w −1
w −1
.
e fe
i −
e∈Si

Due to the connections with the multicommodity flow problems, we will only state the connections through the flow version here. Recall that linear systems in MC2 have the form:






1 0
0 0
1 −1
1
2
1+2
L ⊗
+L ⊗
+L
⊗
,
0 0
0 1
−1 1
63

where E1 , E2 , E1+2 denotes the edges in L1 , L2 , L1+2 , respectively. Edges in E1 ∪ E2 correspond to
individual edges, so it remains to show that edges in E1+2 can correspond to the Hessian matrices
of clusters containing pairs of edges. In a cluster with two edges f 1 and f 2 , the gradient of the
function

−1 2
2
log α − w −1
1 f 1 − w2 f 2
is
2
r
where we treat



−w −1
1 f1
−w −1
2 f2


,

def

−1 2
2
r = α − w −1
1 f 1 − w2 f 2

as a new free variable because we are free to choose α in the IPM instance. Differentiating again
(via the product rule) then gives the Hessian matrix:




−2
2
−4
−2 w −1
w −2
w −1
0
1 f1
1 w 2 f 1f 2
1
+
.
−2
2
w −1
w −2
0
w −1
r2
r
1 w 2 f 1f 2
2 f2
2
As we are free to choose w 1 and w 2 , the second matrix can be any diagonal matrix. Then the
freedom in choosing f 1 and f 2 then means the first matrix can be any rank 1 object. Therefore,
this Hessian is equivalent to the block Hessian for two commodity flows generated in Equation (45)
(up to a change of sign, since we are dealing with a maximization problem here), and its inverse as
given in Equation (46) provides the characterization from Definition 1.5. We can also check more
directly that the 1 + 2 edge can be represented as one of these Hessian inverse blocks.
Claim 9.6. For each edge (i, j) in E1+2 , there exist an edge-vertex incidence matrix N , a diagonal
matrix W and a vector r such that W < r r > and




1 −1
1+2
= N > W − rr> N .
Lij ⊗
−1 1
Proof. For simplicity, we remove all zero rows and columns of L1+2
so that the 2-commodity matrix
ij
for edge (i, j) only has dimension 4 × 4. We pick


1
0


 
 −1 0 
2 0
1


N =
,W =
and r =
.

0
1
0 2
1
0 −1
We check the PSD condition,
W − rr

>


=

1 −1
−1 1


< 0.

Then, we check that this decomposition equals to the 2-commodity matrix


1
0





 −1 0 
1 −1
1 −1 0 0

N > W − rr> N = 
 0
1  −1 1
0 0 1 −1
0 −1




1
1 −1 −1 1


 −1 
  −1 1
1 −1 
1 −1
1+2




1 −1 −1 1 = 
=
= Lij ⊗
.
−1 
1 −1 −1 1 
−1 1
1
−1 1
1 −1

64

Proof of Lemma 3.9. Since the linear system related to the isotropic total variation minimization
matrix is the same as the linear system for 2-commodity, all complexity parameters of these two
linear systems are the same.

10

A Natural Decision Problem: Vector in the Image of a Matrix

When we are given a linear system Ax = c, a natural question is whether the system in fact has
a solution. This is true whenever c lies in the image of A, or equivalently, when ΠA c = c.
We consider a slightly relaxed version of the question of whether c lies in the image of A. We
define the Linear System Decision Problem, abbreviated lsd, as follows:
Definition 10.1 (Linear System Decision Problem, lsd). Given a linear system (A, c) and an
approximation parameter  > 0, we refer to the lsd problem for the triple (A, c, ) as the problem
of outputting
1. yes if there exists x s.t. Ax = c.
2. no if for all x , we have kAx − ck >  kck.
3. Any output is acceptable if neither of the above two cases hold.
Algorithm 12 Reduce lsd to lsa.
1: Find x which is a solution to lsa instance (A, c, ).
2: if kAx − ck ≤  kck then
3:
return yes.
4: else
5:
return no.
6: end if
Lemma 10.2. Algorithm 12 solves lsd instance (A, c, ).
Proof. First, suppose we are in Case 1, i.e. there exists x s.t. Ax = c, or equivalently ΠA c = c.
From this and x being an -approximate solution to (A, c), it follows that
kAx − ck = kAx − ΠA ck ≤  kΠA ck =  kck .
So, the algorithm will return yes.
Second, suppose we are in Case 2, then for all x we have kAx − ck >  kck, so the algorithm
must return no.

10.1

Why Not Solve the Exact Decision Problem?

Earlier in this section, we showed how to solve an approximate linear system decision problem, lsd
by a reduction to an lsa problem. However, given matrix A and a vector c with polynomially
bounded integer entries and condition number, it is natural to ask whether solving lsa to reasonably
high accuracy could in fact let us solve exactly the decision problem of whether there exists an x
s.t. Ax = c. Note that this is equivalent to (I − ΠA ) c = 0. Lemma 4.8 implies that when
A and c have integer valued entries, either ΠA c = 0 or kΠA ck ≥ 1/σmax (A). If a similar
separation result was true for k(I − ΠA ) ck, we could use this to solve the exact decision problem.
65

However, as we will show below, there exists a matrix and vector pair A ∈ Z(n+1)×n and a vector
c ∈ Zn+1 whose entries are polynomially bounded (in fact bounded between -1 and 2), and with
polynomially bounded σmin (A) and σmax (A) s.t. k(I − ΠA ) ck is exponentially small in n. This
suggests that solving the exact decision problem of whether (I − ΠA ) c = 0 would require solving
lsa to exponentially high precision, or adopting some radically different approach.
We now outline the construction of A and c. Each column i of A has two nonzero entries:
Ai,i = 2 and Ai,i+1 = −1. The rank of A is n. Let r be the only unit vector of null(A> ). Then,
I − ΠA = r r > .
We compute r . Since A> r = 0, we have
r i+1 = 2r i , ∀1 ≤ i ≤ n.
Since kr k2 = 1, we have r 1 =
are 0. Then,

1
.
2n+2 −1

Let c be a vector whose first entry is 1 and all other entries

k(I − ΠA ) ck2 = r > c

2

=

1
2n+2

−1

,

which is exponentially small.
We show that the nonzero singular values of A are bounded by constants. Note that A> A is
a tridiagonal matrix such that,
1. all entries on the main diagonal are 5,
2. all entries on the first diagonal below the main diagonal are -2, and
3. all entries on the first diagonal below the main diagonal are -2.
A> A has full rank. A> A can be written as L + D, where L is the graph Laplacian matrix for
an undirected path with identical edge weight 2, and D is a diagonal matrix whose first and last
diagonal entries are 3 and all other diagonal entries are 1. Thus,
L + I 4 A> A 4 L + 3I .
Since λmax (L) = 4, we have for all eigenvalues of λ(A> A) that


1 ≤ λ A> A ≤ 11,
that is,
1 ≤ σmin (A) ≤ σmax (A) ≤

√

11.

Acknowledgements
We thank Richard Peng and Daniel Spielman for helpful comments and discussions and we thank
Richard Peng for suggesting the example in Section 10.1 which demonstrates that the projection
of a vector onto the null space of a matrix can be exponentially small despite well-conditionedness
assumptions on the matrix-vector pair.

66

A

Linear Algebra Background

Proof of Lemma 2.4. First observe
A> Ax − A> c

2
(A> A)†

= (A> Ax − A> c)> (A> A)† (A> Ax − A> c)
= (Ax − c)> Π A (Ax − c)
= (Π A Ax − Π A c)> (Π A Ax − Π A c)
= kAx − Π A ck22 .

Taking x = 0 it follows that kΠ A ck = A> c (A> A)† , and combining our observations immediately
gives the second part of the Fact.
A similar argument gives claims for kx − x ∗ kA> A .
kx − x ∗ k2A> A = (x − x ∗ )> A> A(x − x ∗ )
= x > A> Ax − 2x > A> Ax ∗ + (x ∗ )> A> Ax ∗
= x > A> Ax − 2x > A> Π A c + c > Π >
AΠ Ac
= kAx − Π A ck22 .
The third equality uses the fact that Ax ∗ = Π A c. This completes the proof.

B

Using Complexity Parameters

In this section, we bound the norms, eigenvalues of matrix A ∈ Rn×m and the norms of vector
c ∈ Rn , using the complexity parameters defined in Definition 2.5:
1. s : nnz(A),


2. U : max kAkmax , kckmax , min+1 (A) , min1+ (c) ,
3. K : κ(A).
Note n, m ≤ s.
Claim B.1. kAk∞ ≤ sU .
√
Claim B.2. kck2 ≤ sU .
Proof. This follows the relations of norms,
kck2 ≤

Claim B.3.

1
sU 2

√

n kck∞ ≤

√

≤ λmax (A> A) ≤ sU 2 and λmin (A> A) ≥

sU.

1
.
sK 2 U 2

Proof. By the definition of Frobenius norm,
min{m,n}

kAk2F

=

X

A2ij

=

i,j

X
i=1

67

λi (A> A).

Thus,
λmax (A> A) ≤ kAk2F ≤ sU 2 ,
and
λmax (A> A) ≥

kAk2F
1
≥
.
min{m, n}
sU 2

By the definition of κ(A> A),
λmin (A> A) =

1
λmax (A> A)
≥
.
>
sK 2 U 2
κ(A A)

This completes the proof.

References
[ABFM14]

Boris Alexeev, Afonso S Bandeira, Matthew Fickus, and Dustin G Mixon. Phase
retrieval with polarization. SIAM Journal on Imaging Sciences, 7(1):35–66, 2014.

[ANKKS+ 12] Mica Arie-Nachimson, Shahar Z Kovalsky, Ira Kemelmacher-Shlizerman, Amit
Singer, and Ronen Basri. Global motion estimation from point matches. In 3D
Imaging, Modeling, Processing, Visualization and Transmission (3DIMPVT), 2012
Second International Conference on, pages 81–88. IEEE, 2012.
[BCPT05]

Erik G Boman, Doron Chen, Ojas Parekh, and Sivan Toledo. On factor width and
symmetric h-matrices. Linear algebra and its applications, 405:239–248, 2005.

[BV04]

S. Boyd and L. Vandenberghe. Convex Optimization. Camebridge University Press,
2004. Available at https://web.stanford.edu/˜boyd/cvxbook/bv cvxbook.pdf.

[CFM+ 14]

Michael B. Cohen, Brittany Terese Fasy, Gary L. Miller, Amir Nayyeri, Richard
Peng, and Noel Walkington. Solving 1-laplacians in nearly linear time: Collapsing
and expanding a topological ball. In Proceedings of the Twenty-Fifth Annual ACMSIAM Symposium on Discrete Algorithms, SODA 2014, Portland, Oregon, USA,
January 5-7, 2014, pages 204–216, 2014.

[CKM+ 14]

Michael B. Cohen, Rasmus Kyng, Gary L. Miller, Jakub W. Pachocki, Richard Peng,
Anup B. Rao, and Shen Chen Xu. Solving sdd linear systems in nearly mlog1/2n
time. In Proceedings of the Forty-sixth Annual ACM Symposium on Theory of Computing, STOC ’14, pages 343–352, New York, NY, USA, 2014. ACM.

[CLM+ 16]

Michael B Cohen, Yin Tat Lee, Gary Miller, Jakub Pachocki, and Aaron Sidford.
Geometric median in nearly linear time. In Proceedings of the 48th Annual ACM
SIGACT Symposium on Theory of Computing, pages 9–21. ACM, 2016. Available
at: https://arxiv.org/abs/1606.05225.

[CMMP13a]

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, ITCS ’13, pages 269–282, New York, NY, USA,
2013. ACM.

68

[CMMP13b] 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.
[CMSV17]

Michael B. Cohen, Aleksander Madry, Piotr Sankowski, and Adrian Vladu. Negativeweight shortest paths and unit capacity minimum cost flow in Õ(m10/7 log w) time:
(extended abstract). In Proceedings of the Twenty-Eighth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA ’17, pages 752–771, Philadelphia, PA, USA,
2017. Society for Industrial and Applied Mathematics.

[CS05]

Tony Chan and Jianhong Shen. Image Processing and Analysis: Variational, PDE,
Wavelet, and Stochastic Methods. Society for Industrial and Applied Mathematics,
Philadelphia, PA, USA, 2005.

[DS07]

Samuel I. Daitch and Daniel A. Spielman. Support-graph preconditioners for 2dimensional trusses. CoRR, abs/cs/0703119, 2007.

[DS08]

Samuel I. Daitch and Daniel A. Spielman. Faster approximate lossy generalized flow
via interior point algorithms. In Proceedings of the 40th annual ACM symposium on
Theory of computing, STOC ’08, pages 451–460, New York, NY, USA, 2008. ACM.
Available at http://arxiv.org/abs/0803.0988.

[Fle00]

Lisa K. Fleischer. Approximating fractional multicommodity flow independent of the
number of commodities. SIAM Journal on Discrete Mathematics, 13:505–520, 2000.

[GK98]

Naveen Garg and Jochen Könemann. Faster and simpler algorithms for multicommodity flow and other fractional packing problems. In In Proceedings of the 39th
Annual Symposium on Foundations of Computer Science, pages 300–309, 1998.

[GY04]

Donald Goldfarb and Wotao Yin. Second-order cone programming methods for total
variation-based image restoration. SIAM J. Sci. Comput, 27:622–645, 2004.

[GY05]

Donald Goldfarb and Wotao Yin. Second-order cone programming methods for
total variation-based image restoration. SIAM Journal on Scientific Computing,
27(2):622–645, 2005.

[KLOS14]

Jonathan A. Kelner, Yin Tat Lee, Lorenzo Orecchia, and Aaron Sidford. An almostlinear-time algorithm for approximate max flow in undirected graphs, and its multicommodity generalizations. In Proceedings of the Twenty-fifth Annual ACM-SIAM
Symposium on Discrete Algorithms, SODA ’14, pages 217–226, Philadelphia, PA,
USA, 2014. Society for Industrial and Applied Mathematics.

[KLP+ 16]

Rasmus Kyng, Yin Tat Lee, Richard Peng, Sushant Sachdeva, and Daniel A. Spielman. Sparsified cholesky and multigrid solvers for connection laplacians. In Proceedings of the Forty-eighth Annual ACM Symposium on Theory of Computing, STOC
’16, pages 842–850, New York, NY, USA, 2016. ACM.

[KMP12]

Jonathan A. Kelner, Gary L. Miller, and Richard Peng. Faster approximate multicommodity flow using quadratically coupled flows. In Proceedings of the 44th symposium on Theory of Computing, STOC ’12, pages 1–18, New York, NY, USA, 2012.
ACM. Available at http://arxiv.org/abs/1202.3367.

69

[KMT11]

Ioannis Koutis, Gary L Miller, and David Tolliver. Combinatorial preconditioners
and multilevel solvers for problems in computer vision and image processing. Computer Vision and Image Understanding, 115(12):1638–1646, 2011.

[LMP+ 91]

Tom Leighton, Fillia Makedon, Serge Plotkin, Clifford Stein, Eva Tardos, and Spyros
Tragoudas. Fast approximation algorithms for multicommodity flow problems. In
JOURNAL OF COMPUTER AND SYSTEM SCIENCES, pages 487–496, 1991.

[LS14]

Yin Tat Lee and Aaron Sidford.
Path finding methods for linear programming: Solv√
ing linear programs in Õ( rank) iterations and faster algorithms for maximum flow.
In Foundations of Computer Science (FOCS), 2014 IEEE 55th Annual Symposium
on, pages 424–433. IEEE, 2014. Available at http://arxiv.org/abs/1312.6677 and
http://arxiv.org/abs/1312.6713.

[LS15]

Yin Tat Lee and Aaron Sidford. Efficient inverse maintenance and faster algorithms for linear programming. In Foundations of Computer Science (FOCS),
2015 IEEE 56th Annual Symposium on, pages 230–249. IEEE, 2015. Available at:
https://arxiv.org/abs/1503.01752.

[Mad10]

Aleksander Madry. Faster approximation schemes for fractional multicommodity
flow problems via dynamic graph algorithms. In STOC ’10: Proceedings of the 42nd
ACM symposium on Theory of computing, pages 121–130, New York, NY, USA,
2010. ACM.

[Mad13]

Aleksander Madry. Navigating central path with electrical flows: From flows to
matchings, and back. In Proceedings of the 2013 IEEE 54th Annual Symposium on
Foundations of Computer Science, FOCS ’13, pages 253–262, Washington, DC, USA,
2013. IEEE Computer Society.

[Mad16]

Aleksander Madry. Computing maximum flow with augmenting electrical flows. In
Proceedings of the 2016 IEEE 57th Annual Symposium on Foundations of Computer
Science, FOCS ’16, Washington, DC, USA, 2016. IEEE Computer Society.

[MTW14]

Stefano Marchesini, Yu-Chao Tu, and Hau-tieng Wu. Alternating projection, ptychographic imaging and phase synchronization. arXiv preprint arXiv:1402.0550, 2014.

[Nem]

Arkadi Nemirovski. Interior point polynomial time methods in convex programming.

[OSB15]

Onur Ozyesil, Amit Singer, and Ronen Basri. Stable camera motion estimation using
convex programming. SIAM Journal on Imaging Sciences, 8(2):1220–1262, 2015.

[Pen16]

Richard Peng. Approximate undirected maximum flows in o(mpoly log(n)) time.
In Proceedings of the Twenty-seventh Annual ACM-SIAM Symposium on Discrete
Algorithms, SODA ’16, pages 1862–1867, Philadelphia, PA, USA, 2016. Society for
Industrial and Applied Mathematics.

[Saa03]

Y. Saad. Iterative Methods for Sparse Linear Systems. Society for Industrial and
Applied Mathematics, Philadelphia, PA, USA, 2nd edition, 2003. Available at
http://www-users.cs.umn.edu/˜saad/toc.pdf.

[She13]

Jonah Sherman. Nearly maximum flows in nearly linear time. In Proceedings of the
2013 IEEE 54th Annual Symposium on Foundations of Computer Science, FOCS
’13, pages 263–269, Washington, DC, USA, 2013. IEEE Computer Society.
70

[Spi07]

Daniel A Spielman. Spectral graph theory and its applications. In Foundations of
Computer Science, 2007. FOCS’07. 48th Annual IEEE Symposium on, pages 29–38.
IEEE, 2007.

[Spi16]

Daniel A. Spielman. Nsf award 1562041: Generalized algebraic graph theory: Algorithms and analysis. ALGORITHMIC FOUNDATIONS, 2016.

[SS11]

Amit Singer and Yoel Shkolnisky. Three-dimensional structure determination from
common lines in cryo-em by eigenvectors and semidefinite programming. SIAM
journal on imaging sciences, 4(2):543–572, 2011.

[SS12]

Yoel Shkolnisky and Amit Singer. Viewing direction estimation in cryo-em using
synchronization. SIAM journal on imaging sciences, 5(3):1088–1110, 2012.

[ST08]

Gil Shklarski and Sivan Toledo. Rigidity in finite-element matrices: Sufficient conditions for the rigidity of structures and substructures. SIAM J. Matrix Analysis
Applications, 30(1):7–40, 2008.

[ST14]

D. Spielman and S. Teng.
Nearly linear time algorithms for preconditioning and solving symmetric, diagonally dominant linear systems. SIAM Journal on Matrix Analysis and Applications, 35(3):835–885, 2014. Available at
http://arxiv.org/abs/cs/0607105.

[Str69]

Volker Strassen. Gaussian elimination is not optimal. Numer. Math., 13(4):354–356,
August 1969.

[SZS+ 08]

Richard Szeliski, Ramin Zabih, Daniel Scharstein, Olga Veksler, Vladimir Kolmogorov, Aseem Agarwala, Marshall Tappen, and Carsten Rother. A comparative
study of energy minimization methods for markov random fields with smoothnessbased priors. IEEE Transactions on Pattern Analysis and Machine Intelligence,
30:1068–1080, 2008.

[TBI97]

Lloyd N Trefethen and David Bau III. Numerical linear algebra, volume 50. Siam,
1997.

[Ten10]

Shang-Hua Teng. The Laplacian Paradigm: Emerging Algorithms for Massive
Graphs. In Theory and Applications of Models of Computation, pages 2–14, 2010.

[Vai89]

P. M. Vaidya. Speeding-up linear programming using fast matrix multiplication. In
Proceedings of the 30th Annual Symposium on Foundations of Computer Science,
pages 332–337, Washington, DC, USA, 1989. IEEE Computer Society.

[Wil12]

Virginia Vassilevska Williams. Multiplying matrices faster than coppersmithwinograd. In Proceedings of the Forty-fourth Annual ACM Symposium on Theory of
Computing, STOC ’12, pages 887–898, New York, NY, USA, 2012. ACM.

[Woo14]

David P Woodruff. Sketching as a tool for numerical linear algebra. Theoretical Computer Science, 10(1-2):1–157, 2014. Available at http://arxiv.org/abs/1411.4357.

[WR07]

B. Wohlberg and P. Rodriguez. An iteratively reweighted norm algorithm for minimization of total variation functionals. Signal Processing Letters, IEEE, 14(12):948
–951, dec. 2007.
71

[Wri97]

Stephen J Wright. Primal-dual interior-point methods. SIAM, 1997.

[WW10]

Virginia Vassilevska Williams and Ryan Williams. Subcubic equivalences between
path, matrix and triangle problems. In Proceedings of the 2010 IEEE 51st Annual
Symposium on Foundations of Computer Science, FOCS ’10, pages 645–654, Washington, DC, USA, 2010. IEEE Computer Society.

[WYYZ08]

Yilun Wang, Junfeng Yang, Wotao Yin, and Yin Zhang. A new alternating minimization algorithm for total variation image reconstruction. SIAM Journal on Imaging
Sciences, 1(3):248–272, 2008.

[Ye11]

Yinyu Ye. Interior point algorithms: theory and analysis, volume 44. John Wiley &
Sons, 2011. Available at: http://web.stanford.edu/˜yyye/main.ps.

[ZS14]

Zhizhen Zhao and Amit Singer. Rotationally invariant image representation for
viewing direction classification in cryo-em. Journal of structural biology, 186(1):153–
166, 2014.

72

