Matrix Balancing in Lp Norms:
A New Analysis of Osborne’s Iteration

arXiv:1606.08083v1 [] 26 Jun 2016

Rafail Ostrovsky
UCLA
rafail@cs.ucla.edu∗

Yuval Rabani
The Hebrew University of Jerusalem
yrabani@cs.huji.ac.il†

Arman Yousefi
UCLA
armany@cs.ucla.edu∗

June 28, 2016

Abstract
We study an iterative matrix conditioning algorithm due to Osborne (1960). The goal of the algorithm
is to convert a square matrix into a balanced matrix where every row and corresponding column have
the same norm. The original algorithm was proposed for balancing rows and columns in the L2 norm,
and it works by iterating over balancing a row-column pair in fixed round-robin order. Variants of the
algorithm for other norms have been heavily studied and are implemented as standard preconditioners
in many numerical linear algebra packages. Recently, Schulman and Sinclair (2015), in a first result of
its kind for any norm, analyzed the rate of convergence of a variant of Osborne’s algorithm that uses the
L∞ norm and a different order of choosing row-column pairs. In this paper we study matrix balancing in
the L1 norm and other Lp norms. We show the following results for any matrix A = (aij )ni,j=1 , resolving
in particular a main open problem mentioned by Schulman and Sinclair.
1. We analyze the iteration for the L1 norm under a greedy order of balancing. We show that it
converges to an ǫ-balanced matrix in K = O(min{ǫ−2 log w, ǫ−1 n3/2 log(w/ǫ)}) iterations that cost
a total of O(m + Kn log n) arithmetic operations
over O(n log w)-bit numbers. Here m is the
P
number of non-zero entries of A, and w = i,j |aij |/amin with amin = min{|aij | : aij 6= 0}.

2. We show that the original round-robin implementation converges to an ǫ-balanced matrix in
O(ǫ−2 n2 log w) iterations totalling O(ǫ−2 mn log w) arithmetic operations over O(n log w)-bit numbers.

3. We show that a random implementation of the iteration converges to an ǫ-balanced matrix in
O(ǫ−2 log w) iterations using O(m+ǫ−2 n log w) arithmetric operations over O(log(wn/ǫ))-bit numbers.
√
4. We demonstrate a lower bound of Ω(1/ ǫ) on the convergence rate of any implementation of the
iteration.
5. We observe, through a known trivial reduction, that our results for L1 balancing apply to any Lp
norm for all finite p, at the cost of increasing the number of iterations by only a factor of p.
We note that our techniques are very different from those used by Schulman and Sinclair.

∗

Research supported in part by NSF grants 1065276, 1118126 and 1136174, US-Israel BSF grants, OKAWA Foundation
Research Award, IBM Faculty Research Award, Xerox Faculty Research Award, B. John Garrick Foundation Award, Teradata
Research Award, and Lockheed-Martin Corporation Research Award. This material is also based upon work supported in
part by DARPA Safeware program. The views expressed are those of the authors and do not reflect the official policy or
position of the Department of Defense or the U.S. Government.
†
Research supported in part by ISF grant 956-15, by BSF grant 2012333, and by I-CORE Algo.

1

Introduction

Let A = (aij )n×n be a square matrix with real entries, and let k·k be a given norm. For an index i ∈ [n], let
kai,. k and ka.,i k, respectively, denote the norms of the ith row and the ith column of A, respectively. The
matrix A is balanced in k · k iff ka.,i k = kai,. k for all i. An invertible diagonal matrix D = diag(d1 , . . . , dn )
is said to balance a matrix A iff DAD −1 is balanced. A matrix A is balanceable in k · k iff there exists a
diagonal matrix D that balances it.
Osborne [8] studied the above problem in the L2 norm and considered its application in preconditioning
a given matrix in order to increase the accuracy of the computation of its eigenvalues. The motivation is
that standard linear algebra algorithms that are used to compute eigenvalues are numerically unstable for
unbalanced matrices; diagonal balancing addresses this issue by obtaining a balanced matrix that has the
same eigenvalues as the original matrix, as DAD −1 and A have the same eigenvalues. Osborne suggested
an iterative algorithm for finding a diagonal matrix D that balances a matrix A, and also proved that his
algorithm converges in the limit. He also observed that if a diagonal matrix D = diag(d1 , . . . , dn ) balances
a matrix A, then the diagonal vector d = (d1 , . . . , dn ) minimizes the Frobenius norm of the matrix DAD −1 .
Osborne’s classic algorithm is an iteration that at each step balances a row and its corresponding column
by scaling them appropriately. More specifically the algorithm balances row-columnppairs in a fixed cyclic
order. In orderpto balance row and column i, the algorithm scales the ith row by ka.,i k/kai,. k and the
ith column by kai,. k/ka.,i k. Osborne’s algorithm converges to a unique balanced matrix, but there have
been no upper bounds on the converges rate of Osborne’s algorithm for the L2 norm prior to our work.
Parlett and Reinsch [9] generalized Osborne’s algorithm to other norms without proving convergence.
The L1 version of the algorithm has been studied extensively. The convergence in the limit of the L1
version was proved by Grad [4], uniqueness of the balanced matrix by Hartfiel [5], and a characterization
of balanceable matrices was given by Eaves et al. [3]. Again, there have been no upper bounds on the
running time of the L1 version of the iteration. The first polynomial time algorithm for balancing a
matrix in the L1 norm was given by Kalantari, Khachiyan, and Shokoufandeh [6]. Their approach is
different from the iterative algorithm of Osborne-Parlett-Reinsch. They reduce the balancing problem to
a convex optimization problem and then solve that problem approximately using
P the ellipsoid algorithm.
4
Their algorithm runs in O(n log(n log w/ǫ)) arithemtic operations where w = i,j |ai,j |/amin for amin =
min{|aij | : aij 6= 0} and ǫ is the relative imbalance of the output matrix (see Definition 1).
For matrix balancing in the L∞ norm, Schneider and Schneider [11] gave an O(n4 )-time non-iterative
algorithm. This running time was improved to O(mn + n2 log n) by Young, Tarjan, and Orlin [14]. Despite
the existence of polynomial time algorithms for balancing in the L1 and L∞ norms, and the lack of any
theoretical bounds on the running time of the Osborne-Parlett-Reinsch (OPR) iterative algorithm, the
latter is favored in practice, and the Parlett and Reinsch variant [9] is implemented as a standard in almost
all linear algebra packages (see Chen [2, Section 3.1], also the book [10, Chapter 11] and the code in [1]).
One reason is that iterative methods usually perform well in practice and run for far fewer iterations than
are needed in the worst case. Another advantage of iterative algorithms is that they are simple, they
provide steady partial progress, and they can always generate a matrix that is sufficiently balanced for the
subsequent linear algebra computation.
Motivated by the impact of the OPR algorithm and the lack of any theoretical bounds on its running
time, Schulman and Sinclair [12] recently showed the first bound on the convergence rate of a modified version of this algorithm in the L∞ norm. They prove that their modified algorithm converges in
O(n3 log(ρn/ǫ)) balancing steps where ρ measures the initial imbalance of A and ǫ is the target imbalance
of the output matrix. Their algorithm differs from the original algorithm only in the order of choosing
row-column pairs to balance (we will use the term variant to indicate a deviation from the original roundrobin order). Schulman and Sinclair do not prove any bounds on the running time of the algorithm for
other Lp norms; this was explicitly mentioned as an open problem. Notice that when changing the norm,
1

not only the target balancing condition changes but also the iteration itself, so we cannot deduce an upper
bound on the rate of convergence in the Lp norm from the rate of convergence in the L∞ norm.
In this paper we resolve the open question of [12], and upper bound the convergence rate of the OPR
iteration in any Lp norm.1 Specifically, we show the following bounds for the L1 norm. They imply the
same bounds with an extra factor of p for the Lp norm, by using them on the matrix with entries raised
to the power of p. (Below, the Õ(·) notation hides factors that are logarithmic in various parameters of
the problem. Exact bounds await the statements of the theorems in the following sections.) We show that
the original algorithm (with no modification) converges to an ǫ-balanced matrix in Õ(n2 /ǫ2 ) balancing
steps, using Õ(mn/ǫ2 ) arithmetic operations. We also show that a greedy variant converges in Õ(1/ǫ2 )
balancing steps, using O(m) + Õ(n/ǫ2 ) arithmetic operations; or alternatively in Õ(n3/2 /ǫ) iterations,
using Õ(n5/2 /ǫ) arithmetic operations. Thus, the number of arithmetic operations needed by our greedy
variant is nearly linear in m or nearly linear in 1/ǫ. The near linear dependence on m is significantly
better than the Kalantari-Khachiyan-Shokoufandeh algorithm that uses O(n4 log(n log w/ǫ)) arithmetic
operations (and also the Schulman and Sinclair version with a stricter, yet L∞ , guarantee). For an accurate
comparison we should note that we may need to maintain Õ(n) bits of precision, so the running time is
actually O(m + n2 log n log w/ǫ2 ) (the Kalantari et al. algorithm maintains O(log(wn/ǫ))-bit numbers).
We improve this with yet another, randomized, variant that has similar convergence rate (nearly linear in
m), but needs only O(log(wn/ǫ)) bits of precision. Finally, we show that the dependence on ǫ given by our
√
analyses is within the right ballpark—we demonstrate a lower bound of Ω(1/ ǫ) on the convergence rate
of any variant of the algorithm to an ǫ-balanced matrix. Notice the contrast with the Schulman-Sinclair
upper bound for balancing in the L∞ norm that has O(log(1/ǫ)) dependence on ǫ (this lower bound is for
the Kalantari et al. notion of balancing so it naturally applies also to strict balancing).
Osborne observed that a diagonal matrix D = diag(d1 , . . . , dn ) that balances a matrix A in the L2 norm
also minimizes the Frobenius norm of the matrix DAD −1 . Thus, the balancing problem can be reduced to
minimizing a convex function. Kalantari et al. [6] gave a convex program for balancing in the L1 norm. Our
analysis is based on their convex program. We relate the OPR balancing step to the coordinate descent
method in convex programming. We show that each step reduces the value of the objective function. Our
various bounds are derived through analyzing the progress made in each step. In particular, one of the main
tools in our analysis is an upper bound on the distance to optimality (measured by the convex objective
function) in terms of the the L1 norm of the gradient, which we prove using network flow arguments.
For lack of space, many proofs are missing inline. They appear in Section 7.

2

Preliminaries

In this section we introduce notation and definitions, we discuss some previously known facts and results,
and we prove a couple of useful lemmas.
The problem. Let A = (aij )n×n be a square real matrix, and let k · k be a norm on Rn . For an
index i ∈ [n], let kai,. k and ka.,i k, respectively, denote the norms of the ith row and the ith column of
A, respectively. A matrix A is balanced in k · k iff ka.,i k = kai,. k for all i. An invertible diagonal matrix
D = diag(d1 , . . . , dn ) is said to balance a matrix A iff DAD −1 is balanced. A matrix A is balanceable in
k · k iff there exists a diagonal matrix D that balances it.
For balancing a matrix A in the Lp norm only the absolute values of the entries of A matter, so we
may assume without loss of generality that A is non-negative. Furthermore, balancing a matrix does not
1
It should be noted that the definition of target imbalance ǫ in [12] is stricter than the definition used by [6]. We use the
definition in [6]. This is justified by the fact that the numerical stability of eigenvalue calculations depends on the Frobenius
norm of the balanced matrix, see [9].

2

change its diagonal entries, so if a diagonal matrix D balances A with its diagonal entries replaced by
zeroes, then D balances A too. Thus, for the rest of the paper, we assume without loss of generality that
the given n × n matrix A = (aij ) is non-negative and its diagonal entries are all 0.
A diagonal matrix D = diag(d1 , . . . , dn ) balances A = (aij ) in the Lp norm if and only if D p =
diag(d1 p , . . . , dn p ) balances the matrix A′ = (aij p ) in the L1 norm. Thus, the problem of balancing
matrices in the Lp norm (for any finitie p) reduces to the problem of balancing matrices in the L1 norm;
for the rest of the paper we focus on balancing matrices in the L1 norm.
For an n × n matrix A, we use GA = (V, E, w) to denote the weighted directed graph whose adjacency
matrix is A. More formally, GA is defined as follows. Put V = {1, . . . , n}, put E = {(i, j) : aij 6= 0}, and
put w(i, j) = aij for every (i, j) ∈ E. We use an index i ∈ [n] to refer to both the ith row or column of A,
and to the node i of the digraph GA . Thus, the non-zero entries of the ith column (the ith row, respectively)
correspond to the arcs into (out of, respectively) node i. In the L1 norm it is useful to think of the weight of
an arc as a flow being carried by that arc. Thus, ka.,i k1 is the total flow into vertex i and kai,. k1 is the total
flow out of it. Note that if a matrix A is not balanced then for some nodes i, ka.,i k1 6= kai,. k1 , and thus the
flow on the arcs does not constitute a valid circulation because flow conservation is not maintained. Thus,
the goal of balancing in the L1 norm can be stated as applying diagonal scaling to find a flow function on
the arcs of the graph GA that forms a valid circulation. We use both views of the graph (with arc weights
or flow), and also the matrix terminology, throughout this paper, as convenient.
Without loss of generality we may assume that the undirected graph underlying GA is connected.
Otherwise, after permuting V = {1, . . . , n}, the given matrix A can be replaced by diag(A1 , . . . , Ar ) where
each of A1 , . . . , Ar is a square matrix whose corresponding directed graph is connected. Thus, balancing
A is equivalent to balancing each of A1 , . . . , Ar .
The goal of the iterative algorithm is to balance approximately a matrix A, up to an error term ǫ. We
define the error here.
Definition 1 (approximate balancing). Let ǫ > 0.
√Pn
k −kai,. k1 )2
i=1 (ka
P .,i 1
≤ ǫ.
1. A matrix A is ǫ-balanced iff
ai,j
i,j

2. A diagonal matrix D with positive diagonal entries is said to ǫ-balance A iff DAD −1 is ǫ-balanced.
The algorithms. Kalantari et al. [6] introduced the above definition of ǫ-balancing, and showed that
their algorithm for ǫ-balancing a matrix in the L1 norm uses O(n4 ln((n/ǫ) ln w)) arithmetic operations. In
their recent work, Schulman and Sinclair [12] use, in the context of balancing in the L∞ norm, a stronger
notion of strict balancing (that requires even very low weight row-column pairs to be nearly balanced).
Their iterative algorithm strictly ǫ-balances a matrix in the L∞ norm in O(n3 log(nρ/ǫ)) iterations where
ρ measures the inital imbalance of the matrix. In this paper, we prove upper bounds on the convergence
rate of the Osborne-Parlett-Reinsch (OPR) balancing.
The OPR iterative algorithm balances indices in a fixed round-robin order. Schulman and Sinclair
considered a variant that uses a different rule to choose the next index to balance. We consider in this paper
several alternative implementations of OPR balancing (including the original round-robin implementation)
that differ only in the rule by which an index to balance is chosen at each step. For all rules that we
consider, the iteration generates a sequence A = A(1) , A(2) , . . . , A(t) , . . . of n × n of matrices that converges
to a unique balanced matrix A∗ (see Grad [4] and Hartfiel [5]). The matrix A(t+1) is obtained by balancing
(t)
(t+1) = D (t) A(t) D (t) −1 where D (t) is a
an index of A(t) . If the ith index
q of A is chosen, we get that A
(t)

diagonal matrix with dii =

(t)

(t)

(t)

(t)

(t)

ka.,i k1 /kai,. k1 and djj = 1 for j 6= i. Note that ai,. (a.,i , respectively)

denotes the ith row (ith column, respectively) of A(t) . Also, putting D̄ (1) = In×n and D̄ (t) = D (t−1) · · · D (1)
for t > 1, we get that A(t) = D̄ (t) A(D̄ (t) )−1 .
3

The following lemma shows that each balancing step reduces the sum of entries of the matrix.
Lemma 1. Balancing the ith index
matrix B = (bij )n×n (with bii = 0) decreases the
p of a non-negative
p
total sum of the entries of B by ( kb.,i k1 − kbi,. k1 )2 .

Proof. Before balancing, the total sum of entries
k1 +kb.,i k1 . Balp in the ith row and in the ith column is kbi,.p
ancing scales the entries of the ith column by kbi,. k1 /kb.,i k1 and entries of the ith row by kb.,i k1 /kbi,. k1 .
Thus, after p
balancing the sum of entries in the ith column, which equals the sum of entries in the ith row,
is equal to kbi,. k1 · kb.,i k1 . The entries that are not in the
P balanced row and columnpare not changed.
Therefore, keeping in mind that bii = 0, balancing decreases i,j bij by kb.,i k1 +kbi,. k1 −2 kbi,. k1 · kb.,i k1 =
p
p
( kb.,i k1 − kbi,. k1 )2 .

A reduction to convex optimization. Kalantari et al. [6], as part of their algorithm, reduce matrix
balancing to a convex optimization problem. We overview their reduction here. Our starting point is
Osborne’s observation that if a diagonal matrix D = diag(d1 , . . . , dn ) balances a matrix A in the L2 norm,
then the diagonal vector d = (d1 , . . . , dn ) minimizes the Frobenius norm of the matrix DAD −1 . The
analogous claim for the L1 norm is that if a diagonal matrix D = diag(d1 , . . . , dn ) balances
A
P a matrix
di
in the L1 norm, then the diagonal vector d = (d1 , . . . , dn ) minimizes the function F (d) = i,j aij dj . On
the other hand, Eaves et al. [3] observed that a matrix A can be balanced if and only if the digraph GA is
strongly connected. The following theorem [6, Theorem 1] summarizes the above discussion.
Theorem 1 (Kalantari et al.). Let A = (aij )n×n be a real non-negative matrix, aii = 0, for all i = 1, . . . n,
such that the undirected graph underlying GA is connected. Then, the following statements are equivalent.
(i) A is balanceable (i.e., there exists a diagonal matrix D such that DAD −1 is balanced).
(ii) GA is strongly connected.
P
(iii) Let F (d) = (i,j)∈E aij ddji . There is a point d∗ ∈ Ω = {d ∈ Rn : di > 0, i = 1, . . . , n} such that
F (d∗ ) = inf{F (d) : d ∈ Ω}.
We refer the reader to [6, Theorem 1] for a proof. We have the following corollary.
Corollary 1. d∗ minimizes F over Ω if and only if D ∗ = diag(d∗1 , . . . , d∗n ) balances A.
∗

∂F (d )
∗
∗
Proof. P
As F attains its infimum
Pn at d ∈∗ Ω,∗its gradient ∇F satisfies ∇F (d ) = 0. ∗Also, ∂di = 0 if and
n
∗
∗
only if j=1 aij · (di /dj ) = j=1 aji · (dj /di ) for all i ∈ [n]. In other words, ∇F (d ) = 0 if and only if the
matrix D ∗ AD ∗−1 is balanced where D ∗ = diag(d∗1 , . . . , d∗n ). Thus, d∗ minimizes F over Ω if and only if
D ∗ = diag(d∗1 , . . . , d∗n ) balances A.

It can also be shown that under the assumption of Theorem 1, the balancing matrix D ∗ is unique up to
a scalar factor (see Osborne [8] and Eaves et al. [3]). Therefore, the problem of balancing matrix A can be
reduced to optimizing the function F . Since we are optimizing over the set Ω of strictly positive vectors,
we can apply a change of variables d = (ex1 , . . . , exn ) ∈ Rn to obtain a convex objective function:
f (x) = fA (x) =

n
X

aij exi −xj .

(1)

i,j=1

Kalantari et al. [6] use the convex function f because it can be minimized using the ellipsoid algorithm.
We do not need the convexity of f , and use f instead of F only because it is more convenient to work
4

with, and it adds some intuition. Notice that the partial derivative of f with respect to xi is
n

n

j=1

j=1

X
∂f (x) X
aji · exj −xi ,
aij · exi −xj −
=
∂xi

(2)

which is precisely the difference between the L1 norms of the ith row and the ith column of the matrix
DAD −1 , where D = diag(ex1 , . . . , exn ). Also, by definition, the diagonal matrix diag(ex1 , . . . , exn ) ǫbalances A iff
r
2
Pn
Pn Pn
xi −xj −
xj −xi
a
e
a
e
ij
ji
j=1
j=1
i=1
k∇f (x)k2
Pn
=
≤ ǫ.
(3)
x
−x
i
j
f (x)
i,j=1 aij e

We now state and prove a key lemma that our analysis uses. The lemma uses combinatorial flow and
circulation arguments to measure progress by bounding f (x) − f (x∗ ) in terms of k∇f (x)k1 which is a
global measure of imbalances of all vertices.
Lemma 2. Let f be the function defined in Equation (1), and let x∗ be a global minimum of f . Then, for
all x ∈ Rn , f (x) − f (x∗ ) ≤ n2 · k∇f (x)k1 .

Proof. Recall that f (x) = fA (x) is the sum of entries of a matrix B = (bij ) defined by bij = aij · exi −xj .
Notice that f (x) = fB (~0), and f (x∗ ) = fB (x∗∗ ), where x∗∗ = x∗ − x. Alternatively, f (x) is the sum of
flows (or weights) of the arcs of GB , and f (x∗ ) is the sum of flows of the arcs of a graph G∗ (an arc ij
∗
∗
∗
of G∗ carries a flow of aij · exi −xj ). Notice
Pnthat GB and G have the same set of arcs, but with different
weights. By Equation (2), k∇fA (x)k1 = i=1 kb.,i k1 − kbi,. k1 , i.e., it is the sum over all the nodes of GB
of the difference between the flow into the node and flow out of it. Also notice that GB is unbalanced (else
the statement of the lemma is trivial), however G∗ is balanced. Therefore, the arc flows in G∗ , but not
those in GB , form a valid circulation.
Our proof now proceeds in two main steps. In the first step we show a way of reducing the flow on
some arcs of GB , such that the revised flows make every node balanced (and thus form a valid circulation).
We also make sure that the total flow reduction is at most n2 · k∇fA (x)k1 . In the second step we show that
sum of revised flows of all the arcs is a lower bound on f (x∗ ). These two steps together prove the lemma.
We start with the first step. The nodes of GB are not balanced. Let S and T be a partition of the
unbalanced nodes of GB , with S = {i ∈ [n] : kb.,i k1 > kbi,. k1 } and T = {i ∈ [n] : kb.,i k1 < kbi,. k1 }. That
is, the flow into a node in S exceeds the flow out of it, and the flow into a node in T is less than the flow
out of it. We have that
X
X
X
(kb.,i k1 − kbi,. k1 ) −
(kbi,. k1 − kb.,i k1 ) =
(kb.,i k1 − kbi,. k1 ) = 0.
i∈S

i∈T

i∈[n]

Thus, we can view each node i ∈ S as a source with supply kb.,i k1 − kbi,. k1 , and each node i ∈ T as a sink
with demand kbi,. k1 − kb.,i k1 , and the total supply equals the total demand. We now add some weighted
arcs connecting the nodes in S to the nodes in T . These arcs carry the supply at the nodes in S to the
demand at the nodes in T . Note that we may add arcs that are parallel to some existing arcs in GB . Such
arcs can be replaced by adding flow to the parallel existing arcs of GB . In more detail, to compute the
flows of the added arcs (or the added flow to existing arcs), we add arcs inductively as follows. We start
with any pair of nodes i ∈ S and j ∈ T , and add an arc from i to j carrying flow equal to the minimum
between the supply at i and the demand at j. Adding this arc will balance one of its endpoints, but in the
new graph the sum of supplies at the nodes of S is still equal to the sum of demands at the nodes of T , so
we can repeat the process. (Notice that either S or T or both lose one node.) Each additional arc balances

5

at least one unbalanced node, so GB gets balanced by adding
P at most n additional 1arcs from nodes in S
to nodes in T . The total flow on the added arcs is exactly i∈S (kb.,i k1 − kbi,. k1 ) = 2 · k∇f (x)k1 .
Let E ′ be the set of newly added arcs, and let GB ′ be the new graph with arc weights given by B ′ = (b′ij ).
Since GB ′ is balanced, the arc flows form a valid circulation. We next decompose the total flow of arcs into
cycles. Consider a cycle C in GB ′ that contains at least one arc from E ′ (i.e., C ∩ E ′ 6= ∅). Reduce the
flow on all arcs in C by α = minij∈C b′ij . This can be viewed as peeling off from GB ′ a circulation carrying
flow α. This reduces the flow on at least one arc to zero, and the remaining flow on arcs is still a valid
circulation, so we can repeat the process. It can be repeated as long as there is positive flow on some arc
in E ′ . Eliminating the flow on all arcs in E ′ using cycles reduces the total flow on the arcs by at most n
times the total initial flow on the arcs in E ′ (i.e., n2 · k∇f (x(1) )k1 ), because each cycle contains at most n
arcs and its flow α that is peeled off reduces the flow on at least one arc in E ′ by α. After peeling off all
the flow on all arcs in E ′ , all the arcs with positive flow are original arcs of GB . Let GB ′′ be the graph
with the remaining arcs and their flows which are given by B ′′ = (b′′ij ). The total flow on the arcs of GB ′′
is at least f (x) + 12 · k∇f (x)k1 − n2 · k∇f (x)k1 ≥ f (x) − n2 · k∇f (x)k1 .
Next we show that the total flow on the arcs of GB ′′ is a lower bound on f (x∗ ). Our key tool for
this is the fact that balancing operations preserve the product of arc flows on any cycle in the original graph GB , because balancing a node i multiplies the flow on the arcs into i by some factor r and
the flow on the arcs out of i by 1r . Thus, the geometric mean of the flows of the arcs on any cycle is
not changed by a balancing operation. The arc flows in GB ′′ form a valid circulation, and thus can be
decomposed into flow cycles C1 , . . . , Cq by a similar peeling-off process that was described earlier. Let
n1 , . . . , nq P
be the lengths of cycles, and let α1 , . . . , αq be their flows. The total flow on arcs in GB ′′ is,
therefore, qk=1 nk ·αk . Notice that, by construction, b′′ij ≤ bij , and the decomposition into cycles gives that
P
P
P
P
P
∗∗
∗∗
∗∗
∗∗
∗∗
∗∗
b′′ij = k:ij∈Ck αk . Thus, f (x∗ ) = ni,j=1 bij exi −xj ≥ ni,j=1 b′′ij exi −xj = ni,j=1 k:ij∈Ck αk exi −xj =

Q
P
P
Pq
Pq P
∗∗ 1/nk
∗∗
x∗∗
x∗∗
i −xj
i −xj
= qk=1 nk αk = ni,j=1 b′′ij , where the last
α
e
α
e
≥
n
ij∈Ck k
ij∈Ck k
k=1 k
k=1
inequality uses the arithmetic-geometric mean inequality. Notice that the right-hand side is the total flow
on the arcs of GB ′′ , which is at least f (x) − n2 · k∇f (x(1) )k1 . Thus, f (x∗ ) ≥ f (x) − n2 · k∇f (x)k1 , and this
completes the proof of the lemma.

3

Greedy Balancing

Here we present and analyze a greedy variant of the OPR iteration. Instead of balancing indices in a fixed
round-robin order, the greedy modification chooses at iteration t an index it of A(t) such that balancing
the chosen index results in the largest decrease in the sum of entries of A(t) . In other words, we pick it
such that the following equation holds.
2
q
q
(t)
(t)
ka.,i k1 − kai,. k1
(4)
it = arg max
i∈[n]

We give two analyses of this variant, one that shows that the number of operations is nearly linear
in the size of GA , and another that shows that the number of operations is nearly linear in 1/ǫ. More
specifically, we prove the following theorem.
Theorem 2. Given an n × n matrix A, let m = |E(GA )|, the greedy implementation of the OPR iterative
algorithm outputs an ǫ-balanced matrix in K iterations which
 cost a total of O(m + Knlog n) arithmetic
operations over O(n log w)-bit numbers, where K = O min ǫ−2 log w, ǫ−1 n3/2 log(w/ǫ) .

The proof uses the convex optimization framework introduced in Section 2. Recall that A(t) =
(t)
(t)
D̄ (t) A(D̄ (t) )−1 . If we let D̄ (t) = diag(ex1 , . . . , exn ), the iterative sequence can be viewed as generating a
6

(t)

(t)

sequence of points x(1) , x(2) , . . . , x(t) , . . . in Rn , where x(t) = (x1 , . . . , xn ) and A(t) = D̄ (t) A(D̄ (t) )−1 =
(t)

(t)

(t)

(aij exi −xj )n×n . Initially, x(1) = (0, . . . , 0), and x(t+1) = x(t) + αt ei , where αt = ln(dii ) and ei is the ith
vector of the standard basis for Rn . By Equation (1), the value f (x(t) ) is sum of the entries of the matrix
A(t) . The following key lemma allows us to lower bound the decrease in the value of f (x(t) ) in terms of a
value that can be later related to the stopping condition.
Lemma 3. If index it defined in Equation (4) is picked to balance A(t) , then f (x(t) ) − f (x(t+1)) ) ≥
k∇f (x(t) )k22
.
4f (x(t) )

Corollary 2. If matrix A(t) is not ǫ-balanced, by balancing index it at iteration t, we have f (x(t) ) −
2
f (x(t+1)) ) ≥ ǫ4 · f (x(t) ).
Proof of Theorem 2. By Corollary 2, while A(t) is not ǫ-balanced,
there

 exists an index it to balance such
2
ǫ
2
that f (x(t) ) − f (x(t+1)) ) ≥ ǫ4 · f (x(t) ). Thus, f (x(t+1) ) ≤ 1 −
· f (x(t) ). Iterating for t steps yields
4


P
2 t
f (x(t+1) ) ≤ 1 − ǫ4 ·f (x(1) ). So, on the one hand, f (x(1) ) = ni,j=1 aij since f (x(1) ) is the sum of entries

in A(1) . On the other hand, we argue that the value of f (x(t+1) ) is at least min(i,j)∈E aij . To see this,
consider a directed cycle in the graph GA . It’s easy to see that balancing operations preserve the product of
weights of the arcs on any cycle. Thus, the weight
one arcin the cycle
least
t Pis at least its weight in the
 of at
t
ǫ2
ǫ2
(t+1)
(1)
input matrix A. Therefore, amin ≤ f (x
) ≤ 1 − 4 ·f (x ) = 1 − 4 · ni,j=1 aij . Thus, t ≤ ǫ42 ·ln w
and this is an upper bound on the number of balancing operations before an ǫ-balanced matrix is obtained.
The algorithm initially computes ka.,i k1 and kai,. k1 for all i ∈ [n] in O(m) time. Also the algorithm initially
p
p
2
kai,. k1 − ka.,i k1 for all i in O(m) time and inserts the values in a priority
computes the value of
2
q
q
(t)
(t)
(t)
(t)
are updated
kai,. k1 − ka.,i k1
queue in O(n log n) time. The values of kai,. k1 , ka.,i k1 for all i and
after each balancing operation. In each iteration the weights of at most n arcs change. Updating the
2
q
q
(t)
(t)
(t)
(t)
kai,. k1 − ka.,i k1
involves
values of kai,. k1 and ka.,i k1 takes O(n) time and updating the values of

at most n updates of values in the priority queue, each taking time O(log n). Thus, the first iteration takes
O(m) operations and each iteration after that takes O(n log n) operations, so the total running time of the
algorithm in terms of arithmetic operations is O(m + (n log n log w)/ǫ2 ).
√
An alternative analysis completes the proof. Notice that k∇f (x(t) )k2 ≤ k∇f (x(t) )k1 ≤ n·k∇f (x(t) )k2 .
k∇f (x(t) )k2

(t)

(t)

(x )k2
(x )k2
√
Therefore, f (x(t) )−f (x(t+1)) ) ≥ 4f (x(t) ) 2 ≥ k∇f
·k∇f (x(t) )k1 ≥ 2n13/2 · k∇f
·(f (x(t) )−f (x∗ )),
4 n·f (x(t) )
f (x(t) )
where the first inequality follows from Lemma 3, and the last inequality follows from Lemma 2. Therefore,
(x(t) )k2
> ǫ), we have that f (x(t) ) − f (x(t+1)) ≥ 2nǫ3/2 · (f (x(t) ) − f (x∗ )).
while At is not ǫ-balanced (so k∇f
f (x(t) )


Rearranging the terms, we get f (x(t+1) ) − f (x∗ ) ≤ 1 − 2nǫ3/2 · (f (x(t) ) − f (x∗ )). Therefore, f (x(t+1) ) −

t
f (x∗ ) ≤ 1 − 2nǫ3/2 · (f (x(1) ) − f (x∗ )). Notice that by Lemma 3, f (x(t+1) ) − f (x∗ ) ≥ f (x(t+1) ) −
2
2


(t+1) )k
(x(t+1) )k2
2
(t+1) ) ≥ k∇f (x
·
f
(x
· amin . On the other hand, f (x(1) ) − f (x∗ ) ≤
f (x(t+2) ) ≥ k∇f
(t+1)
(t+1)
2f (x
)
2f (x
)
P
(x(t+1) )k2
f (x(1) ) ≤ ni,j=1 aij . Thus, for t = 2ǫ−1 · n3/2 ln(4w/ǫ2 ), we have that k∇f
≤ ǫ, so the matrix is
f (x(t+1) )
ǫ-balanced.

7

4

Round-Robin Balancing (the original algorithm)

Recall that original Osborne-Parlett-Reinsch algorithm balances indices in a fixed round-robin order. Although the greedy variant of the OPR iteration is a simple modification of the implementation, the convergence rate of the original algorithm (with no change) is interesting. This is important because the original
algorithm has a slightly simpler implementation, and also because this is the implementation used in almost all numerical linear algebra software including MATLAB, LAPACK and EISPACK (refer to [13, 7]
for further background). We give some answer to this question in the following theorem.
Theorem 3. Given an n × n matrix A, the original implementation of the OPR iteration outputs an ǫbalanced matrix in O(ǫ−2 n2 log w) iterations totalling O(ǫ−2 mn log w) arithmetic operations over O(n log w)bit numbers (m is the number of non-zero entries of A).

5

Randomized Balancing

In Theorem 2 the arithmetic operations were applied to O(n ln w)-bit numbers. This will cause an additional
factor of O(n ln w) in the running time of the algorithm. In this section we fix this issue by presenting a
randomized variant of the algorithm that applies arithmetic operations to numbers of O(ln(wn/ǫ)) bits.
Thus, we obtain a algorithm for balancing that
in nearly linear time. While the greedy algorithm
p runs p
works by picking the node i that maximizes ( kai,. k − ka.,i k)2 , the key idea of the randomized algorithm
is sampling a node for balancing using sampling probabilities that do not depend on the difference in arc
weights (the algorithm uses low-precision rounded weights, so this can affect significantly the difference).
Instead, our sampling probabilities depend on the sum of weights of the arcs incident on a node.
We first introduce some notation. We use O(ln(wn/ǫ)) bits of precision to approximate xi -s with x̂i -s.
(t) (t)
(t)
b(t) = (x̂1 , x̂2 , . . . , x̂n ) at every time t.
Thus, xi − 2−O(ln(nw/ǫ)) ≤ x̂i ≤ xi . In addition to maintaining x
(t)

(t)

(t)

(t)

The algorithm also maintains for every i and j the value of b
aij which is aij = aij ex̂i −x̂j truncated to
O(ln(wn/ǫ)) bits of precision. We set the hidden constant to give a truncation error of r = (ǫ/wn)10 amin , so
P
P
(t)
(t)
(t)
(t)
(t)
(t)
(t)
aij − r ≤ b
aij ≤ aij . The algorithm also maintains for every i, kb
ai,. k = nj=1 b
aij and kb
a.,i k = nj=1 b
aji .
P
P
(t)
(t)
(t)
For every i, we use the notation kai,. k = nj=1 aij and ka.,i k = nj=1 aji . Note that the algorithm does
(t)

(t)

(t)

not maintain the values aij , ka.,i k or kai,. k.
The algorithm works as follows (see the pseudo-code of Algorithm 1 that appears in Section 7). In each
(t)

(t)

kb
ai,. k+kb
a.,i k
P
(t) .
2 i,j b
aij

iteration it samples an index i with probability pi =

If i is sampled, a balancing operation

is applied to index i only if the arcs incident on i have significant weight, and i’s imbalance is sufficiently
(t)
(t)
(t)
(t)
large. Put M̂i = max{kb
ai,. k, kb
a.,i k} and put m̂i = min{kb
ai,. k, kb
a.,i k}. The imbalance is considered large if
m̂i = 0 (this can happen because of the low precision), or if m̂i 6= 0 and
(t)

(t)

(t)

M̂i
m̂i

≥ 1 + nǫ . A balancing operation

bi , where α = 12 ln(kb
is done by adding α to x
a.,i k/kb
ai,. k), unless m̂i = 0, in which case we replace the 0
value by nr. This updates the weights of the arcs incident on i. Also, the L1 norms of changed rows and
columns are updated. (For convenience we use in this section k · k instead of k · k1 to denote the L1 norm.)
Note that in the pseudo-code, ← indicates an assignment where the value on the right-hand side is
computed to O(ln(wn/ǫ)) bits of precision. Thus, we have
(t+1)

and

α − (ǫ/wn)10 ≤ x
bi

(t+1)

aij ex̂i

(t+1)

−x̂j

(t+1)

−r ≤b
aij
8

(t)

−x
bi ≤ α
(t+1)

≤ aij ex̂i

(5)
(t+1)

−x̂j

,

(t+1)

aji ex̂i

Theorem 4. With probability at least
an ǫ-balanced matrix.

(t+1)

−x̂j
9
10 ,

(t+1)

−r ≤b
aji

(t+1)

≤ aji ex̂i

(t+1)

−x̂j

.

Algorithm 1 returns in time O(m ln

P

ij

aij + ǫ−2 n ln(wn/ǫ) ln w)

The idea of proof is to show that in every iteration of the algorithm we reduce f (.) by at least a factor
of 1 − Ω(ǫ2 ). Before we prove the theorem, we state and prove a couple of useful lemmas.
(t)
(t)
Fix an iteration t, and define three sets of indices as follows: A = {i : kb
ai,. k + kb
a.,i k ≥ ǫamin /10wn},
B = {i : m̂i 6= 0 ∧ M̂i /m̂i ≥ 1 + ǫ/n}, and C = {i : m̂i = 0}. If the random index i satisfies i ∈
/ A or
i ∈ A \ (B ∪ C), the algorithm does not perform any balancing operation on i. The following lemma states
that the expected decrease due to balancing such indices is small, and thus skipping them does not affect
the speed of convergence substantially.
q
2
q
P
2
(t)
(t)
< 2ǫn · f (x(t) ), where p is the
pi ·
kai,. k − ka.,i k
Lemma 4. For every iteration t, i∈A∩(B∪C)
/
probability distribution over indices at time t.

We now show a lower bound on the decrease in f (·), if a node i ∈ A ∩ (B ∪ C) is balanced.
2
q
q
(t)
(t)
1
(t+1)
(t)
Lemma 5. If i ∈ A∩(B∪C) is balanced in iteration t, then f (b
x
)−f (b
x ) ≥ 10 ·
kai,. k − ka.,i k .
Proof of Theorem 4. By Lemma 5, the expected decrease in f (.) in iteration t is lower bounded as follows.
E[f (b
x(t) ) − f (b
x(t+1) )] ≥

X

i∈A∩(B∪C)

=

pi ·

1
10

q
2
q
(t)
(t)
kai,. k − ka.,i k


q
2
q
n
1 X
(t)
(t)
pi ·
·
kai,. k − ka.,i k −
10
i=1

The second term can be bounded, using Lemma 4, by
For the first term, we can write
n
X
i=1

pi ·

q

(t)
kai,. k −

q

(t)
ka.,i k

2

≥
=

≥
=

where the penultimate inequality holds because

M̂i
Mi

P

n
X
i=1

X

i∈A∩(B∪C)
/

i∈A∩(B∪C)
/

q

(t)

pi ·

pi ·

(t)
kai,. k

−

q
q

(t)

kai,. k −
2

(t)
ka.,i k

q

≤

(t)

ka.,i k

2ǫ2
n

2



.

· f (b
x(t) ).

(t)

(kai,. k − ka.,i k)2
(t)

(t)

2(kai,. k + ka.,i k)

(t)
(t)
(t)
(t)
n
X
kb
ai,. k + kb
a.,i k (kai,. k − ka.,i k)2
·
P (t)
(t)
(t)
2 ij b
aij
2(kai,. k + ka.,i k)
i=1
(t)
(t) 2
n
1 X (kai,. k − ka.,i k)
·
P (t)
16
a
i=1

ij

k∇f (b
x(t) )k22
16f (b
x(t) )

≥ 12 , so

9

≥

ǫ2

16

(t)

(t)

(t)

(t)

kb
ai,. k+kb
a.,i k
kai,. k+ka.,i k

ij

· f (b
x(t) ),
≥

M̂i
2Mi

≥ 14 , and the last inequality

holds as long as the matrix is not ǫ-balanced, so

E[f (b
x(t) ) − f (b
x(t+1) )] ≥
≥

k∇f (b
x(t) )k2
f (b
x(t) )

≥ ǫ. Combining everything together, we get



q
q
2
2
q
q
n
X
1 X
(t)
(t)
(t)
(t)
·
pi ·
pi ·
kai,. k − ka.,i k −
kai,. k − ka.,i k 
10
i=1
i∈A∩(B∪C)
/
 2

ǫ
2ǫ2
ǫ2
1
·
· f (b
x(t) ) −
· f (b
x((t)) ) ≥
· f (b
x(t) ),
10
16
n
320

where the last inequality assumes n ≥ 64. This implies that the expected number of iterations to obtain
9
an ǫ-balanced matrix is O(ǫ−2 ln w). Markov’s inequality implies that with probability 10
an ǫ-balanced
−2
matrix is obtained in O(ǫ ln w) iterations. It P
is easy to see that each iteration of the algorithm takes
O(n ln(wn/ǫ)) time. Initializations take O(m ln ij aij ) time. So the total running time of the algorithm
P
is O(m ln ij aij + ǫ−2 n ln(wn/ǫ) ln w).

6

A Lower Bound on the Rate of Convergence

In this section we prove the following lower bound.
Theorem 5. There are matrices for which all variants of the Osborne-Parlett-Reinsch iteration (i.e.,
√
regardless of the order of indices chosen to balance) require Ω(1/ ǫ) iterations to balance the matrix to the
relative error of ǫ.
Before proving this theorem, we present the claimed construction. Let A be the following 4 × 4 matrix,
and let A∗ denote the corresponding fully-balanced matrix.




0
1
0
0
0 1
0
0
p
1 0 β + ǫ 0
1
ǫ(β + ǫ) 0
0
∗




p
A=
,
A
=
0
0 ǫ
0
1
ǫ(β + ǫ)
0
1
0 0
1
0
0
0
1
0
Here ǫ > 0 is arbitrarily small, and β = 100ǫ. It’s easy to see that A∗ = D ∗ AD ∗ −1 where
!
r
r
β+ǫ
β+ǫ
D = diag 1, 1,
.
,
ǫ
ǫ

√
To prove Theorem 5, we show that balancing A to the relative error of ǫ requires Ω(1/ ǫ) iterations,
regardless of the order of balancing operations. Notice that in order to fully balance A, we simply need to
replace a23 and a32 by their geometric mean. We measure the rate of convergence using the ratio a32 /a23 .
1
ǫ
= 101
. When the matrix is fully balanced, the ratio becomes 1. We show that
This ratio is initially β+ǫ
this ratio increases by a small factor in each iteration, and that it has to increase sufficiently for the matrix
to be ǫ-balanced. This is summarized in the following two lemmas.
(t+1)
 1 + 7√β  a(t)
a32
· 32
≤
.
Lemma 6 (change in ratio). (t+1)
(t)
1
+
ǫ
a23
a23
(t)

Lemma 7 (stopping condition). If A(t) is ǫ-balanced, then

a32

(t)
a23

>

1
.
100

Before proving the two lemmas we show how they lead to the proof of Theorem 5.
10

 √ t
 √ t
(t+1)
a
1+7 β
1+7 β
a32
ǫ
Proof of Theorem 5. By Lemma 6, 32
≤
·
· β+ǫ
=
. By Lemma 7, if A(t+1) is
(t+1)
1+ǫ
a23
1+ǫ
a23
 √ t
(t+1)
√ t ǫ
a
1+7 β
ǫ
1
·
< 32
≤
1
+
7
ǫ-balanced, then 100
β · β+ǫ . Using β = 100ǫ, we get the condition
(t+1) ≤
1+ǫ
β+ǫ
√ t 101 a23
√
that (1 + 7 β) > 100 , which implies that t = Ω(1/ ǫ).
P
Proof of Lemma 6. Using the notation we defined earlier, we have that f (x(1) ) = 4i,j=1 aij = 4 + 2ǫ + β
p
P
and f (x∗ ) = 4i,j=1 a∗ij = 4 + 2 ǫ(β + ǫ), so f (x(1) ) − f (x∗ ) < β. We observe that at each iteration t,
(t) (t)

(t) (t)

(t) (t)

a12 a21 = a34 a43 = 1 and a23 a32 = ǫ(β + ǫ) because the product of weights of arcs on any cycle in GA is
preserved (for instance, arcs (1, 2) and (2, 1) form a cycle and initially a12 a21 = 1).
(t) (t)
The ratio a32 /a23 is only affected in iterations that balance index 2 or 3. Let’s assume a balancing
operation at index 2, a similar analysis applies to balancing at index 3. By balancing at index 2 at time t
we have
v
u (t)
(t+1)
(t)
u a + a(t)
a32
a23
23
t 21
=
=
.
(6)
(t)
(t+1)
(t)
(t)
a32
a23
a12 + a32

Thus, to prove Lemma 6, it suffices to show that
(t)

(t+1)

a32

(t)

a32

·

a23

(t+1)

a23

(t) (t)

=

(t)

(t)

(t)

(t)

a21 + a23
a12 + a32

(t)

√
1+7 β
≤
.
1+ǫ
(t)

(7)
(t) (t)

By our previous observation, a12 a21 = 1, so if a21 = y, then a12 = 1/y. Similarly a23 a32 = ǫ(β + ǫ) implies
(t)
(t)
that there exists z such that a23 = (β + ǫ)z and a32 = ǫ/z. Therefore:
(t)

(t)

(t)
a12

(t)
a32

a21 + a23
+

=

y + (β + ǫ)z
(1/y) + (ǫ/z)

(8)

We bound
√ the right hand side of Equation (8) by proving upper bounds on y and z. We first show that
y < 1 + 2 β. To see this notice that on the one hand,
f (x(t) ) =

4
X

i,j=1

(t)

(t)

(t)

(t)

(t)

(t)

(t)

aij = a12 + a21 + a23 + a32 + a34 + a43 ≥ y +

p
1
+ 2 ǫ(β + ǫ) + 2,
y

(9)

p
(t)
(t)
(t) (t)
where we used a34 + a43 ≥ 2 and a34 a43 ≥ 2 ǫ(β + ǫ), both implied by the arithmetic-geometric mean
inequality. On the other hand,
p
f (x(t) ) ≤ f (x(1) ) ≤ f (x∗ ) + β = 4 + 2 ǫ(β + ǫ) + β.
(10)

Combining Equations (9) and (10) together, we have y + (1/y) − 2 ≤ β. For sufficiently small ǫ, the last
inequality√implies, in particular, that y < 2. Thus, we have (y − 1)2 ≤ yβ < 2β, and this implies that
y < 1 + 2 β.
Next we show that z ≤ 1. Assume for contradiction that z > 1. By the arithmetic-geometric mean
(t)
(t)
(t)
(t)
inequality a12 + a21 ≥ 2 and a34 + a43 ≥ 2. Thus,
f (x(t) ) =

4
X

i,j=1

(t)

aij ≥ 2 + (β + ǫ)z +



ǫ
1
+ 2 = 4 + βz + ǫ z +
> 4 + β + 2ǫ = f (x(1) ),
z
z

11

where the last inequality follows because z > 1, and z + 1/z > 2. But this is a contradiction, because each
balancing iteration reduces the value of F , so f (x(t) ) ≤ f (x(1) ).
(t)
(t)
(t)
(t)
We can now bound (a21 + a23 )/(a12 + a32 ). By Equation (8), and using our bounds for y and z,
√
√
√
(t)
(t)
1+4 β
1+7 β
a21 + a23
y + (β + ǫ)z
(1 + 2 β) + (β + ǫ)
≤
≤
=
≤
.
(t)
(t)
ǫ
1
1
(1/y) + (ǫ/z)
1+ǫ
a12 + a32
√ +ǫ
√ +
√
1+2 β
1+2 β 1+2 β
√
The last line uses the fact that β ≫ β = 100ǫ ≥ ǫ, which holds if ǫ is sufficiently small.
Proof of Lemma 7. Let t − 1 be the last iteration before an ǫ-balanced matrix is obtained. We argued
(t)
(t)
that there is z ≤ 1 such that a23 = (β + ǫ)z and a32 = ǫ/z. Assume for the sake of contradiction that
(t) (t)
a32 /a23 < 1/100. This implies that (ǫ/z)/((β + ǫ)z) < 1/100, and thus z 2 > 100/101. So, we get
r
q
q 2
q

1 2
(t)
(t)
(t)
(t)
(t)
(t)
(t)
≥ a23 1 −
a23 − a32
f (x(t) ) − f (x∗ ) ≥ a23 + a32 − 2 a23 a32 =
100
r
100
= 0.81 · (β + ǫ)z ≥ 0.81 · (β + ǫ) ·
≥ 81 · ǫ.
(11)
101
By Lemma 2, the left hand side the of above can be bounded as follows.
f (x(t) ) − f (x∗ ) ≤ nk∇f (x(t) )k1 ≤ n2 k∇f (x(t) )k2

(12)

Note that for sufficiently small ǫ, f (x(t) ) ≤ f (x(1) ) ≤ 5. Combining Equations (11) and (12), and using
n = 4 and f (x(t) ) ≥ 5, we get that
81
k∇f (x(t) )k2
>
· ǫ > ǫ.
(13)
(t)
80
f (x )
By Equation (3), this contradicts our assumption that t − 1 is the last iteration.

7

Proofs

Proof of Lemma 3. The value f (x(t) ) is the sum of the entries of A(t) . By Lemma 1, balancing the i-th
q
2
q
(t)
(t)
(t)
(t)
index of A reduces the value of f (x ) by
ka.,i k1 − kai,. k1 . To simplify notation, we drop the
superscript t in the following equations. We have
q
2
q
(ka.,i k1 − kai,. k1 )2
(ka.,i k1 − kai,. k1 )2
ka.,i k1 − kai,. k1
.
= p
p
2 ≥
2 (ka.,i k1 + kai,. k1 )
ka.,i k1 + kai,. k1

It is easy to see that

Pn
(ka.,i k1 − kai,. k1 )2
(ka.,i k1 − kai,. k1 )2
≥ Pi=1
.
max
n
i∈[n] (ka.,i k1 + kai,. k1 )
i=1 (ka.,i k1 + kai,. k1 )

(14)

(15)

But the right hand side of the above
 resuming the use of the superscript t) equals
 inequality (after
k∇f (x(t) )k22
(t)
(t)
. This is because for all i, kai,. k1 − ka.,i k1 is by Equation (2) the i-th coordinate of ∇f (x(t) ),
2f (x(t) )

Pn  (t)
(t)
= 2f (x(t) ). Together with Equations (14) and (15),
ka
k
+
ka
k
and in the denominator
i,. 1
.,i 1
i=1
(q
2 )
q
(t)
(t)
ka.,i k1 − kai,. k1
decreases f (x(t) ) by the claimed
this implies that balancing it = arg maxi∈[n]
value.

12

Proof of Corollary 2. From Equation (3), we know that the diagonal matrix diag(ex1 , . . . , exn ) balances
k∇f (x)k2
k∇f (x(t) )k2
A with relative error ǫ if and only if
> ǫ. By
≤ ǫ. Thus, if A(t) is not ǫ-balanced,
f (x)
f (x(t) )


2
2
k∇f (x(t) )k2
(x(t) )k2
· f (x(t) ) ≥ ǫ4 · f (x(t) ).
Lemma 3, f (x(t) ) − f (x(t+1)) ) ≥ 4f (x(t) ) 2 = 41 · k∇f
f (x(t) )

Proof of Theorem 3. In the original Osborne-Parlett-Reinsch algorithm, the indices are balanced in a fixed
round-robin order. A round of balancing is a sequence of n balancing operations where each index is
balanced exactly once. Thus, in the OPR algorithm all n indices are balanced in the same order every
round. We prove a more general statement that any algorithm that balances indices in rounds (even if
the indices are not balanced in the same order every round) obtains an ǫ-balanced matrix in at most
O((n log w)/ǫ2 ) rounds. To this end, we show that applying a round of balancing to a matrix that is not
ǫ-balanced reduces the value of function f at least by a factor of 1 − ǫ2 /16n.
To simplify notation, we consider applying a round of balancing to the initial matrix A(1) = A. The
argument clearly holds for any time-t matrix A(t) . If A is not ǫ-balanced, by Lemma 3 and Corollary 2,
there exists an index i such that by balancing i the value of f is reduced by:
(1)

f (x

(2)

) − f (x

)=

q

2
q
ǫ2
ka.,i k1 − kai,. k1
≥ f (x(1) ).
4

(16)

If i is the first index to balance in the next round of balancing, then in that round the value of f is
reduced at least by a factor of 1−ǫ2 /4 ≥ 1−ǫ2 /16n, and we are done. Consider the graph GA corresponding
to the matrix A. If node i is not the first node in GA to be balanced, then some of its neighbors in the
graph GA might be balanced before i. The main problem is that balancing neighbors of i before i may
reduce the imbalance of i significantly, so we cannot argue that when we reach i and balance it the value of
f reduces significantly. Nevertheless, we show that balancing i and its neighbors in this round will reduce
the value of f by at least the desired amount. Let t denote the time that i is balanced in the round. For
(t)
(t)
every arc (j, i) into i, let δj = |aji − aji |, and for every arc (i, j) out of i let σj = |aij − aij |. These values
measure the weight change of these arcs due to balancing a neighbor of i at any time since the beginning
of the round. The next lemma shows if the weight of an arc incident on i has changed since the beginning
of the round, it must have reduced the value of f .
Claim 1. If balancing node j changes aji to aji + δ, then the balancing reduces the value of f by at least
δ2 /aji . Similarly if balancing node j changes aij to aij + δ, then the balancing reduces the value of f by at
least δ2 /aij .
Proof. To simplicity notation we assume that j is balanced in the first iteration of the round. If balancing
j changes aji to aji + δ, then by the definition of balancing,
s
ka.,j k1
aji + δ
=
.
(17)
aji
kaj,. k1
Thus, by Lemma 1 the value of f reduces by
q

ka.,j k1 −

q

kaj,. k1

2

=

s

!2
 2
2

δ
δ2
aji + δ
ka.,j k1
kaj,. k1 ≥
− 1 kaj,. k1 =
− 1 kaj,. k1 =
kaj,. k1
aji
aji
aji

The proof for the second part of the claim is similar.

13

Going back to the proof of Theorem 3, let t denote the iteration in the round that i is balanced. By
Claim 1, balancing neighbors of i has already reduced the value of f by
δj2

X

j:(j,i)∈E

Balancing i reduces value of f by an additional
round is reduced by at least:
R=

X

j:(j,i)∈E

δj2
+
aji

X

aji

j:(i,j)∈E

q

j:(i,j)∈E

σj2

X

+

(t)
ka.,i k1

−

aij
q

.

(18)

(t)
kai,. k1

2

, so the value of f in the current

2
q
q
σj2
(t)
(t)
ka.,i k1 − kai,. k1
+
aij

Assume without loss of generality that kai,. k1 > ka.,i k1 . To lower bound R, we consider two cases:
X
X
1
case (i)
δj +
σj ≥ (kai,. k1 − ka.,i k1 ). In this case,
2
j:(j,i)∈E

j:(i,j)∈E

R≥

X

j:(j,i)∈E

δj2
+
aji

X

j:(i,j)∈E

σj2
1
≥
aij
ka.,i k1

X

δj2 +

j:(j,i)∈E

X

1
kai,. k1

σj2

j:(i,j)∈E

X
X
1
1
≥
(
δj )2 +
(
nka.,i k1
nkai,. k1
j:(j,i)∈E

σj )2 ,

j:(i,j)∈E

where the last inequality follows by Cauchy-Schwarz inequality. By assumption of case (i),
X
X
1
max(
δj ,
σj ) ≥ (kai,. k1 − ka.,i k1 )
4
j:(j,i)∈E

(19)

(20)

j:(i,j)∈E

Equations (19) and (20) together imply that
P
P
( j:(j,i)∈E δj )2 + ( j:(i,j)∈E σj )2
1 (kai,. k1 − ka.,i k1 )2
R≥
≥
n max(ka.,i k1 , kai,. k1 )
16n max(ka.,i k1 , kai,. k1 )
p
p
p
2 p
2
2
q
q
ka.,i k1 − kai,. k1
ka.,i k1 + kai,. k1
1
≥
=
ka.,i k1 − kai,. k1 .
16n max(ka.,i k1 , kai,. k1 )
16n
case (ii)

X

δj +

j:(j,i)∈E

X

j:(i,j)∈E

1
σj < (kai,. k1 − ka.,i k1 ). By definition of δj ’s and σj ’s:
2
ka.,i k1 −
kai,. k1 −

X

j:(j,i)∈E

X

j:(i,j)∈E

X

(t)

δj ≤ ka.,i k1 ≤ ka.,i k1 +
(t)

σj ≤ kai,. k1 ≤ kai,. k1 +

δj

(21)

σj .

(22)

j:(j,i)∈E

X

j:(i,j)∈E

Combining Equations (21) and (22), and the assumption of case (ii) gives:
X
X
(t)
(t)
kai,. k1 + ka.,i k1 ≤ kai,. k1 + ka.,i k1 +
σj +
δj ≤ 2 (kai,. k1 + ka.,i k1 )
(t)

(t)

kai,. k1 − ka.,i k1 ≥ kai,. k1 − ka.,i k1 −

j:(i,j)∈E

j:(j,i)∈E

X

X

j:(i,j)∈E

14

σj −

j:(j,i)∈E

δj ≥

1
(kai,. k1 − ka.,i k1 ) .
2

(23)
(24)

Using Equations (23) and (24), we can write:

R≥

q

(t)

kai,. k1 −

q

(t)

ka.,i k1

2

(t)
(t)
ka.,i k1 − kai,. k1
(kai,. k1 − ka.,i k1 )2

= q
2 ≥  (t)
q
(t)
(t)
(t)
k
k
+
ka
8
ka
.,i 1
i,. 1
ka.,i k1 + kai,. k1

2

(kai,. k1 − ka.,i k1 )2
1
≥
≥
16 (kai,. k1 + ka.,i k1 )
16

q

kai,. k1 −

q

ka.,i k1

2

.

Thus, we have shown in both cases that in one round the balancing operations on node i and its neighbors
reduces the value of f by at least
q
2
q
1
ka.,i k1 − kai,. k1 ,
(25)
16n
2

which in turn is at least Ω( ǫn f (x(1) )) by Equation (16). Thus, we have shown that if A is not ǫ-balanced,
one round of balancing (where
index is balanced exactly once) reduces the objective function f by
 2 each 
ǫ
(1)
a factor of at least 1 − Ω n f (x ) . By an argument similar to the one in the proof of Theorem 2,

we get that the algorithm obtains an ǫ-balanced matrix in at most O(ǫ−2 n log w) rounds. The number of
balancing iterations in each round is n, and the number of arithmetic operations in each round is O(m), so
the original OPR algorithm obtains an ǫ-balanced matrix using O(ǫ−2 mn log w) arithmetic operations.
Proof of Lemma 4. Notice that for every i,
(t)

(t)

(kb
ai,. k + kb
a.,i k) ·
(t)

q

(t)

kai,. k −

(t)

q

(t)

(t)

ka.,i k

2

(t)

(t)

≤ (kb
ai,. k + kb
a.,i k) ·


2
(t)
(t)
kai,. k − ka.,i k
(t)

(t)

kai,. k + ka.,i k


2
(t)
(t)
≤ kai,. k − ka.,i k ,

(t)

because kb
ai,. k + kb
a.,i k ≤ kai,. k + ka.,i k. We first bound the sum over i ∈
/ A.
X
i∈A
/

pi ·

q

(t)
kai,. k

−

q

2

(t)
ka.,i k

2
q
(t)
(t)
q
X kb
ai,. k + kb
a.,i k
(t)
(t)
kai,. k − ka.,i k
·
P
(t)
2 i,j b
aij
i∈A
/
2
X  (t)
1
(t)
·
ka
k
−
ka
k
≤
P
i,.
.,i
(t)
2 i,j b
aij i∈A
/
2
X  (t)
1
(t)
≤
·
kb
a
k
+
kb
a
k
+
2nr
P
i,.
.,i
(t)
2 i,j b
aij i∈A
/
X
1
·
(2ǫamin /10wn)2
≤
P
(t)
2 i,j b
aij i∈A
/

=

≤

2

P

1

(t)

i,j

b
aij

· n · (ǫamin /5wn)2 ≤
(t)

ǫ2
ǫ2
· amin ≤
· f (x(t) )
25n
25n

(t)

(t)

(t)

where the second inequality follows because, for every j, aij ≤ b
aij + r and aji ≤ b
aji + r, and the third
(t)

(t)

inequality follows because kb
ai,. k + kb
a.,i k < ǫamin /10wn and nr < ǫamin /10wn.
(t)

(t)

(t)

(t)

Next, we bound the sum over i ∈ A\(B∪C). Recall M̂i = max{kb
ai,. k, kb
a.,i k} and m̂i = min{kb
ai,. k, kb
a.,i k}.
(t)

(t)

(t)

(t)

Put Mi = max{kai,. k, ka.,i k} and mi = min{kai,. k, ka.,i k}. Let k = arg maxi∈A\(B∪C) (Mi − mi )2 . We
15

Algorithm 1 RandomBalance(A, ǫ)
Input: Matrix A ∈ Rn×n , ǫ
Output: An ǫ-balanced matrix
r = amin · (ǫ/wn)10
(1)
(1)
2: Let âij = aij for all i and j
1:

(t)

(t)

(t)

(t)

Let kb
ai,. k = kai,. k and kb
a.,i k = ka.,i k for all i
−2
4: for t = 1 to O(ǫ ln w) do
3:

(t)

5:

Pick i randomly with probability pi =

6:

if kb
ai,. k + kb
a.,i k ≥ ǫamin /10wn then

(t)

(t)

(t)

8:
9:

(t)

(t)

12:
13:
14:
15:

(t+1)

b
aij

16:

bits)

20:
21:
22:
23:
24:

(t)

1
2
1
2

(t)

ln(nr/kb
ai,. k)
(t)

else if m̂i = kb
ai,. k = 0 then α = ln(kb
a.,i k/nr)
end if
b(t+1) ← x
b(t) + αei (truncated to O(ln(wn/ǫ)) bits of precision)
Let x
for j = 1 to n do
if j is a neighbor of i then

11:

19:

(t)

else if m̂i = kb
a.,i k = 0 then α =

10:

18:

(t)

M̂i = max{kb
ai,. k, kb
a.,i k}, m̂i = min{kb
ai,. k, kb
a.,i k}
if m̂i = 0 or M̂i /m̂i ≥ 1 + ǫ/n then
(t)
(t)
a.,i k/kb
ai,. k)
if m̂i 6= 0 then α = 12 ln(kb

7:

17:

(t)

kb
ai,. k+kb
a.,i k
P
(t)
2 i,j b
aij

(t+1)

← aij ex̂i

(t+1)

(t+1)

−x̂j

(t)

(t)

(t+1)

(t+1)

← aji ex̂j

and b
aji
(t+1)

(t+1)

(t+1)

−x̂i
(t)

, (truncated to O(ln(wn/ǫ))
(t)

(t+1)

kb
aj,. k = kb
aj,. k − b
aji + b
aji
and kb
a.,j k = kb
a.,j k − b
aij + b
aij
end if
end for P
P
(t+1)
(t+1)
(t+1)
(t+1)
kb
ai,. k = nj=1 b
aij
and kb
a.,i k = nj=1 b
aji
end if
end if
end for
return the resulting matrix

have
X

i∈A\(B∪C)

pi ·

q

(t)

kai,. k −

q

2
(t)
≤
ka.,i k
≤
≤

2
2

2

P
P
P

To bound the last quantity, we prove an upper bound on
(t)

(t)

we have M̂k + m̂k = kb
ak,. k + kb
a.,k k ≥

ǫamin
10wn .

Thus, M̂k ≥

16

1
(t)

i,j

1

b
aij

(t)
aij
i,j b

1

(t)

i,j

b
aij

·
·

X

i∈A\(B∪C)

X

i∈A\(B∪C)

·n

· m2k




2
(t)
(t)
kai,. k − ka.,i k
(Mi − mi )2

2
Mk
−1 .
mk

M̂k
Mk
ǫ
mk using the fact that m̂k < 1 + n . As k ∈ A,
M̂k
ǫamin
ǫ
20wn . Combining this with m̂k < 1 + n implies

ǫamin
40wn .

that m̂k > 12 M̂k ≥

Hence,

 ǫ 9 M̂
M̂k + nr
M̂k
M̂k
Mk
nr
40ǫ9
2ǫ
Mk
k
≤
≤
≤
≤
+
=
+ 40n ·
+ 8 ≤1+ .
mk
m̂k
m̂k
m̂k
ǫamin /40wn
m̂k
wn
m̂k
n
n
(Notice that w ≥ 1.) Using the upper bound on
X

i∈A\(B∪C)

pi ·

q

(t)
kai,. k −

q

Mk
mk ,

(t)
ka.,i k

we obtain

2

≤

2

P

1
(t)

i,j

1

b
aij

2
Mk
−1
·n
mk
 2
2ǫ
2
· n · mk
n
· m2k



≤

2

≤

ǫ2
2ǫ2
· mk ≤
· f (x(t) ),
n
n

P

(t)

i,j

b
aij

where the penultimate inequality uses the fact that mk ≤ m̂k + nr ≤ m̂k +
P
(t)
aij .
i,j b
Proof of Lemma 5. We will assume that ǫ <

1
10 .

ǫamin
40wn

< 2m̂k ≤ M̂k + m̂k ≤

We first consider the case that i ∈ A ∩ B (notice that
(t)

(t+1)

B ∩ C = ∅). The update using O(ln(wn/ǫ)) bits of precision gives x
bi + α − (ǫ/nw)10 ≤ x
bi
so
v
v
u (t)
u (t)
u kb
u kb
k
a
(t+1)
(t)
10
.,i
x
b
x
b
−(ǫ/wn)
t a.,i k · exb(t)
t
i
i
i .
≤
≤
e
·
e
(t)
(t)
kb
ai,. k
kb
ai,. k

(t)

≤x
bi + α,

Therefore,

(t+1)

kai,.

k=

n
X

(t+1)

aij exbi

(t+1)

−b
xj

j=1

and
(t+1)

ka.,i

k=

n
X
j=1

(t+1)
(t+1)
x
bj
−b
xi

aji e

v
v
v
u (t)
u (t)
u (t)
n
X
u
u a.,i k
u kb
(t)
(t)
a.,i k
kb
a.,i k
(t) t kb
x
bi −b
xj
t
t
a
e
=
·
·
ka
k
,
≤
ij
i,.
(t)
(t)
(t)
kb
ai,. k j=1
kb
ai,. k
kb
ai,. k

10

≤ e(ǫ/wn)

v
v
u (t)
u (t)
n
X
u kb
u kb
(t)
(t)
a
k
ai,. k
i,.
(t)
x
b
−b
x
10
i
· t (t) ·
aji e j
≤ (1 + 2(ǫ/wn) )) · t (t) · ka.,i k.
kb
a.,i k j=1
kb
a.,i k

We used the fact that ex ≤ 1 + 2x for x ≤ 21 . We will now use the notation M̂i , m̂i , Mi , and mi
(the reader can recall the definitions from the proof of Lemma 4). We also put δ = 2(ǫ/wn)10 , and
(t+1)

M̂i /m̂i
Mi /mi .

Thus, decrease of function f (·) due to balancing i is f (b
x(t) ) − f (b
x(t+1) ) = Mi + mi − ka.,i k −

 q
q
p
√  √
(t+1)
1/σ + σ · Mi mi =
kai,. k ≥ Mi + mi − (1 + δ) Mi m̂i /M̂i + mi M̂i /m̂i = Mi + mi − (1 + δ)
√
√
√
√
√ 2
Mi − mi − ((1 + δ)/ σ + (1 + δ) σ − 2) · Mi mi . We now consider three cases, and in each case
show that
 p
√
√
√ 2
9 p
Mi − mi .
·
(1 + δ)/ σ + (1 + δ) σ − 2 · Mi mi ≤
10
σ=

17

4

case (i): 1 ≤ σ < 1 + nǫ 2 . We first note that Mi ≥ M̂i ≥
mi
m̂i +nr
1
ǫ
1
≤ 1+ǫ/n
. Since ǫ < 10
, we have
+ nr ≤ 1 − 2n
Mi ≤
M̂i

M̂i +m̂i
2

>

ǫamin
20wn .

Also, mi ≤ m̂i + nr, so

M̂i

9
10



1−

r

mi
Mi

2

r

2
ǫ
9
≥
1− 1−
10
2n
4
4ǫ
≥
2
n




ǫ4
≥
(1 + δ) + (1 + δ) · 1 + 2 − 2
n

 r
√
1+δ
mi
√ + (1 + δ) σ − 2 ·
,
≥
Mi
σ
where the third inequality holds by definition of δ, and the last inequality holds because mi /Mi ≤ 1 and
σ ∈ [1, 1 + ǫ4 /n2 ]. By multiplying both sides of the inequality by Mi we obtain the desired bound.
i
σ < 1. We first prove a lower bound on the value of σ, as follows: M̂
m̂i ≥




20ǫ9
9
Mi
Mi
Mi
nr
nr
(1
−
)
≥
·
1
−
·
1
−
≥
, and thus σ ≥ 1 − 20ǫ
8 . So we have
mi
Mi
mi
Mi
mi
n
8
n

 r
√
1+δ
mi
1+δ
√ + (1 + δ) σ − 2 ·
≤ q
+ (1 + δ) − 2
9
Mi
σ
1 − 20ǫ
n8


20ǫ9
+ (1 + δ) − 2
≤ (1 + δ) 1 + 8
n

2
r
mi
4ǫ4
9
24ǫ9
,
· 1−
< 2 ≤
≤
n8
n
10
Mi

case (ii):

proving the desired inequality in this case. The first inequality holds because

know that

≥

Mi −nr
mi

≥

9

≤ 1 and 1 − 20ǫ
n8 ≤ σ ≤ 1.

ǫ4
.
n2

The idea is to show that Mi /mi is large so the desired inequality follows. We
M̂i
Mi
nr
= m̂i ≤ m̂i and therefore m̂i ≤ mσi . On the other hand, m̂i ≥ mi − nr, so mi ≤ 1−1/σ
.
ǫamin /20wn
Mi
nr
n6
ǫ4
1 − 2n2 , so mi < ǫ4 /2n2 . Also, Mi ≥ ǫamin /20wn. Therefore, mi ≥ 2n3 r/ǫ4 ≥ 40ǫ5 . Next,

case (iii): σ > 1 +
σMi
mi

mi
Mi

M̂i
mi

Clearly, 1/σ <
notice that since m̂i > 0 it must be that m̂i ≥ r. Therefore, mi ≤ m̂i + nr ≤ 2nm̂i . This implies that
M̂i
Mi
Mi
m̂i ≤ m̂i ≤ 2n · mi , so σ ≤ 2n. Finally,
r

√
√
Mi
1
1+δ
√ + (1 + δ) σ − 2 ≤ (1 + δ) · 2n ≤≤
·
,
σ
10
mi
q
Mi
i
with room to spare (using the lower bound on M
).
Multiplying
both
sides
by
mi
mi gives




 r
√
Mi
1+δ
1 Mi
9
√ + (1 + δ) σ − 2 ·
≤
≤
mi
10 mi
10
σ

with more room to spare. This completes proof of the case i ∈ A ∩ B.
18

r

!2
Mi
−1 ,
mi

1
2

We now move on to the case i ∈ A ∩ C, so M̂i + m̂i ≥

(t)
ln(nr/kb
ai,. k)

or α =

(t)
1
a.,i k/nr).
2 ln(kb

ǫamin
10wn

and m̂i = 0. In the algorithm, α =

The idea is that we therefore replace m̂i (which is 0) by nr in some
 q
q 
nr
(t)
(t+1)
i
+ mi M̂
of the equations. In particular, f (b
x ) − f (b
x
) ≥ Mi + mi − (1 + δ) Mi
nr . Note that
M̂i

since m̂i = 0 then mi ≤ nr. Therefore,

so

M̂i
nr

≥

ǫamin /20wn
n(ǫ/wn)10 amin

≥

n8
.
20ǫ9

M̂i
nr

Thus we get

≤

Mi
nr

≤

Mi
mi .

On the other hand, since i ∈ A, M̂i ≥ ǫamin /20wn,


≥
≥
where the third inequality holds because

nr

s



M̂i 
nr
M̂i
r
r !
20ǫ9
Mi
+ mi
Mi + mi − (1 + δ) Mi
8
n
mi
r
20ǫ9
Mi + mi − 2(1 + δ)Mi
n8


4
√
1 p
1
20ǫ
≥ Mi ≥ ( Mi − mi )2 ,
Mi 1 − 4
n
10
10
q
q
q
mi
nr
i
=
M
≤
M
mi M
.
i
i
mi
Mi

f (b
x(t) ) − f (b
x(t+1) ) ≥ Mi + mi − (1 + δ) Mi
≥

r

+ mi

M̂i

References
[1] EISPACK implementation. http://www.netlib.org/eispack/balanc.f.
[2] T.-Y. Chen. Balancing sparse matrices for computing eigenvalues. Master’s thesis, UC Berkeley, May
1998.
[3] B. C. Eaves, A. J. Hoffman, U. G. Rothblum, and H. Schneider. Line-sum-symmetric scalings of
square nonnegative matrices. In Mathematical Programming Essays in Honor of George B. Dantzig
Part II, pages 124–141. Springer, 1985.
[4] J. Grad. Matrix balancing. The Computer Journal, 14(3):280–284, 1971.
[5] D. J. Hartfiel. Concerning diagonal similarity of irreducible matrices. In Proceedings of the American
Mathematical Society, pages 419–425, 1971.
[6] B. Kalantari, L. Khachiyan, and A. Shokoufandeh. On the complexity of matrix balancing. SIAM
Journal on Matrix Analysis and Applications, 118(2):450–463, 1997.
[7] D. Kressner. Numerical methods for general and structured eigenvalue problems. Princeton University
Press, 2005.
[8] E. E. Osborne. On pre-conditioning of matrices. Journal of the ACM (JACM), 7(4):338–345, 1960.
[9] B. N. Parlett and C. Reinsch. Balancing a matrix for calculation of eigenvalues and eigenvectors.
Numerische Mathematik, 13(4):293–304, 1969.
[10] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical Recipes: The Art of
Scientific Computing, 3rd Edition. Cambridge University Press, 2007.
19

[11] H. Schneider and M. H. Schneider. Max-balancing weighted directed graphs and matrix scaling.
Mathematics of Operations Research, 16(1):208–222, February 1991.
[12] L. J. Schulman and A. Sinclair. Analysis of a classical matrix preconditioning algorithm. In Proceedings
of the Forty-Seventh Annual ACM on Symposium on Theory of Computing, pages 831–840, 2015.
[13] L. N. Trefethen and M. Embree. Spectra and pseudospectra: The behavior of nonnormal matrices and
operators. Springer, 2005.
[14] N. E. Young, R. E. Tarjan, and J. B. Orlin. Faster parametric shortest path and minimum-balance
algorithms. Networks, 21(2):205–221, 1991.

20

