Strictly Balancing Matrices in Polynomial Time Using Osborne's Iteration

TL;DR

This paper proves that Osborne's iteration can strictly balance matrices in polynomial time across any finite Lp norm, providing the first such guarantee and improving upon previous weak balancing results.

Contribution

The paper introduces a polynomial-time implementation of Osborne's iteration for strict matrix balancing in any Lp norm, a novel theoretical guarantee.

Findings

Converges to a strictly epsilon-balanced matrix in polynomial time.

First proof of polynomial-time strict balancing in Lp norms.

Improves upon previous weak balancing results.

Abstract

Osborne's iteration is a method for balancing $n\times n$ matrices which is widely used in linear algebra packages, as balancing preserves eigenvalues and stabilizes their numeral computation. The iteration can be implemented in any norm over $\mathbb{R}^n$, but it is normally used in the $L_2$ norm. The choice of norm not only affects the desired balance condition, but also defines the iterated balancing step itself. In this paper we focus on Osborne's iteration in any $L_p$ norm, where $p < \infty$. We design a specific implementation of Osborne's iteration in any $L_p$ norm that converges to a strictly $\epsilon$-balanced matrix in $\tilde{O}(\epsilon^{-2}n^{9} K)$ iterations, where $K$ measures, roughly, the {\em number of bits} required to represent the entries of the input matrix. This is the first result that proves that Osborne's iteration in the $L_2$ norm (or any $L_p$ norm, $p < \infty$) strictly balances matrices in polynomial time. This is a substantial improvement over our recent result (in SODA 2017) that showed weak balancing in $L_p$ norms. Previously, Schulman and Sinclair (STOC 2015) showed strong balancing of Osborne's iteration in the $L_\infty$ norm. Their result does not imply any bounds on strict balancing in other norms.

Equations

$$\frac{\max \left\{\Vert a_{.,i}{\Vert }_1,\Vert a_{i,.}{\Vert }_1\right\}}{\min \left\{\Vert a_{.,i}{\Vert }_1,\Vert a_{i,.}{\Vert }_1\right\}}\le 1+ϵ.$$

$$\frac{\partial f(\mathbf{x})}{\partial x_i}=\sum _{j=1}^n{a}_{ij}\cdot e^{x_i-x_j}-\sum _{j=1}^n{a}_{ji}\cdot e^{x_j-x_i},$$

$$f_A(\mathbf{x})-f_A({\mathbf{x}}^{\prime })$$

$$\Vert b_{.,i}{\Vert }_1+\Vert b_{i,.}{\Vert }_1<\frac{1}{8n^3}\cdot f_A(\vec{0}),$$

$$O\left({ϵ}^{-2}{n}^9\log (wn/ϵ)\log w/\log n\right)$$

$$f^{({\mathcal{B}}_s)}({\mathbf{x}}^{(\mathbf{t})})\le f^{({\mathcal{B}}_s)}({\mathbf{x}}^{({\mathbf{t}}_{\mathbf{s},\mathbf{r}})})\le \left(n-\left|{\mathcal{B}}_{s,t_{s,r}}\right|\right)\cdot {\tau }_s=\left(n-\left|{\mathcal{B}}_{s,t}\right|\right)\cdot {\tau }_s.$$

$$f^{({\mathcal{B}}_s)}({\mathbf{x}}^{({\mathbf{t}}_{\mathbf{s},\mathbf{r}+1})})\le \left(n-\left|{\mathcal{B}}_{s,t_{s,r}}\right|\right)\cdot {\tau }_s+k\cdot {\tau }_s=\left(n-\left|{\mathcal{B}}_{s,t_{s,r+1}}\right|\right)\cdot {\tau }_s.$$

$$f^{({\mathcal{B}}_s)}({\mathbf{x}}^{({\mathbf{t}}_{\mathbf{s}+1})})\le n\cdot {\tau }_s=\frac{1}{4n^2}\cdot f^{({\mathcal{B}}_{s-1})}({\mathbf{x}}^{({\mathbf{t}}_{\mathbf{s}})}),$$

$${\tau }_s=\frac{1}{4n^3}\cdot f^{({\mathcal{B}}_{s-1})}({\mathbf{x}}^{({\mathbf{t}}_{\mathbf{s}})})\le \frac{1}{4n^3}\cdot n{\tau }_{s-1}=\frac{1}{4n^2}\cdot {\tau }_{s-1},$$

$${\tau }_s-\sum _{j=s+1}^{s^{\prime }}n{\tau }_j\ge {\tau }_s\cdot \left(1-n\cdot \sum _{k=1}^{s^{\prime }-s}(2n{)}^{-2k}\right)\ge \frac{{\tau }_s}{2}.$$

$$t^{\prime }-t=O\left({ϵ}^{-2}{n}^7\log (wn/ϵ)\right).$$

$$f^{({\mathcal{B}}_s)}({\mathbf{x}}^{({\mathbf{t}}^{\prime \prime })})-f^{({\mathcal{B}}_s)}({\mathbf{x}}^{({\mathbf{t}}^{\prime \prime }+1)})$$

$$\frac{\Vert \nabla f^{({\mathcal{B}}_s)}({\mathbf{x}}^{({\mathbf{t}}^{\prime \prime })}){\Vert }_1}{f^{({\mathcal{B}}_s)}({\mathbf{x}}^{({\mathbf{t}}^{\prime \prime })})}>{ϵ}^{\prime }=\frac{{ϵ}^2}{64n^4}.$$

$$f^{({\mathcal{B}}_s)}({\mathbf{x}}^{({\mathbf{t}}^{\prime \prime })})-f^{({\mathcal{B}}_s)}({\mathbf{x}}^{({\mathbf{t}}^{\prime \prime }+1))})$$

$$f^{({\mathcal{B}}_s)}({\mathbf{x}}^{({\mathbf{t}}^{\prime \prime }+1)})-f^{({\mathcal{B}}_s)}({\mathbf{x}}^{\ast })\le \left(1-\frac{{ϵ}^{\prime }}{8n^2}\right)\cdot (f^{({\mathcal{B}}_s)}({\mathbf{x}}^{({\mathbf{t}}^{\prime \prime })})-f^{({\mathcal{B}}_s)}({\mathbf{x}}^{\ast })).$$

$$f^{({\mathcal{B}}_s)}({\mathbf{x}}^{(\mathbf{t}+\mathbf{T})})-f^{({\mathcal{B}}_s)}({\mathbf{x}}^{\ast })\le {\left(1-\frac{{ϵ}^{\prime }}{8n^2}\right)}^T\cdot (f^{({\mathcal{B}}_s)}({\mathbf{x}}^{(\mathbf{t})})-f^{({\mathcal{B}}_s)}({\mathbf{x}}^{\ast })).$$

$$f^{({\mathcal{B}}_s)}({\mathbf{x}}^{({\mathbf{t}}^{\prime }-1)})-f^{({\mathcal{B}}_s)}({\mathbf{x}}^{\ast })\le {\left(\frac{{ϵ}^{\prime }}{4\sqrt{n}}\right)}^2\cdot \frac{1}{w^n}\cdot \sum _{i,j=1}^n{a}_{ij}\le {\left(\frac{{ϵ}^{\prime }}{4\sqrt{n}}\right)}^2\cdot f^{({\mathcal{B}}_s)}({\mathbf{x}}^{({\mathbf{t}}^{\prime }-1)}).$$

$$\frac{1}{16n}\cdot \frac{\Vert \nabla f^{({\mathcal{B}}_s)}({\mathbf{x}}^{({\mathbf{t}}^{\prime }-1)}{\Vert }_1^2}{(f^{({\mathcal{B}}_s)}({\mathbf{x}}^{({\mathbf{t}}^{\prime }-1)}){)}^2}\le \frac{f^{({\mathcal{B}}_s)}({\mathbf{x}}^{({\mathbf{t}}^{\prime }-1)})-f^{({\mathcal{B}}_s)}({\mathbf{x}}^{({\mathbf{t}}^{\prime })})}{f^{(\mathcal{B})}({\mathbf{x}}^{({\mathbf{t}}^{\prime }-1)})}\le \frac{f^{({\mathcal{B}}_s)}({\mathbf{x}}^{({\mathbf{t}}^{\prime }-1)})-f^{({\mathcal{B}}_s)}({\mathbf{x}}^{\ast })}{f^{(\mathcal{B})}({\mathbf{x}}^{({\mathbf{t}}^{\prime }-1)})}\le {\left(\frac{{ϵ}^{\prime }}{4\sqrt{n}}\right)}^2,$$

Full text

In this paper we focus on Osborne’s iteration in any $L_{p}$ norm, where $p<\infty$. We design a specific implementation of Osborne’s iteration in any $L_{p}$ norm that converges to a strictly $\epsilon$-balanced matrix in $\tilde{O}(\epsilon^{-2}n^{9}K)$ iterations, where $K$ measures, roughly, the number of bits required to represent the entries of the input matrix.

This is the first result that proves that Osborne’s iteration in the $L_{2}$ norm (or any $L_{p}$ norm, $p<\infty$) strictly balances matrices in polynomial time. This is a substantial improvement over our recent result (in SODA 2017) that showed weak balancing in $L_{p}$ norms. Previously, Schulman and Sinclair (STOC 2015) showed strong balancing of Osborne’s iteration in the $L_{\infty}$ norm. Their result does not imply any bounds on strict balancing in other norms.

1 Introduction

Problem statement and motivation.

This paper analyzes the convergence properties of Osborne’s celebrated iteration [8] for balancing matrices. Given a norm $\|\cdot\|$ in $\mathbb{R}^{n}$, an $n\times n$ matrix $A$ is balanced if and only if for all $i$, the $i$-th row of $A$ and the $i$-th column of $A$ have the same norm. The problem of balancing a matrix $A$ is to compute a diagonal matrix $D$ such that $DAD^{-1}$ is balanced. The main motivation behind this problem is that balancing a matrix does not affect its eigenvalues, and balancing matrices in the $L_{2}$ norm increases the numerical stability of eigenvalue computations [8, 7]. Balancing also has a positive impact on the computational time needed for computing eigenvalues ([7, section 1.4.3]). In practice, it is sufficient to get a good approximation to the balancing problem. For $\alpha\geq 1$, a matrix $B=DAD^{-1}$ is an $\alpha$-approximation to the problem of balancing $A$ if and only if for all $i$, the ratio between the maximum and minimum of the norms of the $i$-th row and column is bounded by $\alpha$. It is desirable to achieve $\alpha=1+\epsilon$ for some small $\epsilon>0$. A matrix $B$ that satisfies this relaxed balancing condition is also said to be strictly $\epsilon$-balanced.

Osborne’s iteration attempts to compute the diagonal matrix $D$ by repeatedly choosing an index $i$ and balancing the $i$-th row and column (this multiplies the $i$-th diagonal entry of $D$ appropriately). Osborne proposed this iteration in the $L_{2}$ norm, and suggested round-robin choice of index to balance. However, other papers consider the iteration in other norms and propose alternative choices of index to balance [10, 13, 9]. Notice that a change of norm not only changes the target balance condition, but also changes the iteration itself, as in each step a row-column pair is balanced in the given norm. An implementation of Osborne’s iteration is used in most numerical algebra packages, including MATLAB, LAPACK, and EISPACK, and is empirically efficient (see [7, 14] for further background). The main theoretical question about Osborne’s iteration is its rate of convergence. How many rounds of the iteration are provably sufficient to get a strictly $\epsilon$-balanced matrix?

Our results.

We consider Osborne’s iteration in $L_{p}$ norms for finite $p$. We design a new simple choice of the iteration (i.e., a rule to choose the next index to balance), and we prove that this variant provides a polynomial time approximation scheme to the balancing problem. More specifically, we show that in the $L_{1}$ norm, our implementation converges to a strictly $\epsilon$-balanced matrix in $O\left(\epsilon^{-2}n^{9}\log(wn/\epsilon)\log w/\log n\right)$ iterations, where $\log w$ is a lower bound on the number of bits required to represent the entries of $A$ (exact definitions await Section 2). The time complexity of these iterations is $O\left(\epsilon^{-2}n^{10}\log(wn/\epsilon)\log w\right)$ arithmetic operations over $O(n\log w)$-bit numbers. This result implies similar bounds for any $L_{p}$ norm where $p$ is fixed, and in particular the important case of $p=2$. This is because applying Osborne’s iteration in the $L_{p}$ norm to $A=(a_{ij})_{n\times n}$ is equivalent to applying the iteration in the $L_{1}$ norm to $(a_{ij}^{p})_{n\times n}$. Of course, the bit representation complexity of the matrix, and thus the bound on the number of iterations, grows by a factor of $p$.

Our results give the first theoretical analysis that indicates that Osborne’s iteration in the $L_{2}$ norm, or any $L_{p}$ norm for finite $p$, is indeed efficient in the worst case. This resolves the question that has been open since 1960. Previously, such a result was obtained only for the $L_{\infty}$ norm [13]. Concerning the convergence rate for the $L_{p}$ norms discussed here, we recently published a result [9] that considers a much weaker notion of approximation. The previous result only shows the rate of convergence to a matrix that is approximately balanced in an average sense. The matrix might still have row-column pairs that are highly unbalanced. The implementations in the common numerical linear algebra packages use as a stopping condition the strict notion of balancing, and not this weaker notion. We discuss previous work in greater detail below.

Previous work.

Osborne [8] studied the $L_{2}$ norm version of matrix balancing, proved the uniqueness of the $L_{2}$ solution, designed the iterative algorithm discussed above, and proved that it converges in the limit to a balanced matrix (without bounding the convergence rate). Parlett and Reinsch [10] generalized Osborne’s iteration to other norms. Their implementation is the one widely used in practice (see Chen [2, Section 3.1], also the book [11, Chapter 11] and the code in [1]). Grad [4] proved convergence in the limit for the $L_{1}$ version (again without bounding the running time), and Hartfiel [5] showed that the $L_{1}$ solution is unique. Eaves et al. [3] gave a characterization of balanceable matrices. Kalantari et al. [6] gave an algorithm for $\epsilon$-balancing a matrix in the $L_{1}$ norm. The algorithm reduces the problem to unconstrained convex optimization and uses the ellipsoid algorithm to approximate the optimal solution. This generates a weakly $\epsilon$-balanced matrix, which satisfies the following definition. Given $\epsilon>0$, a matrix $A=(a_{ij})_{n\times n}$ is weakly $\epsilon$-balanced if and only if $\sqrt{\sum_{i=1}^{n}(\|a_{.,i}\|-\|a_{i,.}\|)^{2}}\leq\epsilon\cdot\sum_{i,j}|a_{i,j}|$. Compare this with the stronger condition of being strictly $\epsilon$-balanced, which we use in this paper, and numerical linear algebra packages use as a stopping condition. This condition requires that for every $i\in\{1,2,\dots,n\}$, $\max\{\|a_{.,i}\|,\|a_{i,.}\|\}\leq(1+\epsilon)\cdot\min\{\|a_{.,i}\|,\|a_{i,.}\|\}$. In $L_{\infty}$, Schneider and Schneider [12] gave a polynomial time algorithm that exactly balances a matrix. The algorithm does not use Osborne’s iteration. Its running time was improved by Young et al. [15]. Both algorithms rely on iterating over computing a minimum mean cycle in a weighted strongly connected digraph, then contracting the cycle. Schulman and Sinclair [13] were the first to provide a quantitative bound on the running time of Osborne’s iteration. They proposed a carefully designed implementation of Osborne’s iteration in the $L_{\infty}$ norm that strictly $\epsilon$-balances an $n\times n$ matrix $A$ in $O(n^{3}\log(\varrho n/\epsilon))$ iterations, where $\varrho$ measures the initial $L_{\infty}$ imbalance of $A$. Their proof is an intricate case analysis. Finally, in [9] we recently proved that a logarithmic dependence on $1/\epsilon$ is impossible in the $L_{1}$ norm (the lower bound is $\Omega(1/\sqrt{\epsilon})$). In the same paper we also showed that several implementations of Osborne’s iteration in $L_{p}$ norms, including the original implementation, converge to a weakly $\epsilon$-balanced matrix in polynomial time (which, in fact, can be either nearly linear in $n$ or nearly linear in $1/\epsilon$). The result of [9] is derived by observing that Osborne’s iteration can be interpreted as an implementation of coordinate descent to optimize the convex function from [6]. This is the starting point of this paper, but to make the approach guarantee strict balancing, we need to revise substantially previous implementations using novel algorithmic ideas. The main difficulty is the need to handle the different scales of row/column norm values; an index may shift between scales over time as a side-effect of balancing other indices. Moreover, the analysis of the convergence rate is more complicated, and requires additional ideas.

Our contribution.

The main difficulty with respect to previous work is the following. The convergence rate of coordinate descent can be bounded effectively as long as there is a choice of coordinate (i.e., index) for which the drop in the objective function in a single step is non-negligible compared with the current objective value. But if this is not the case, then one can argue only about the balance of each index relative to the sum of norms of all rows and columns. Indices that have relatively heavy weight (row norm $+$ column norm) will indeed be balanced at this point. However, light-weight indices may be highly unbalanced. The naive remedy to this problem is to work down by scales. After balancing the matrix globally, heavy-weight indices are balanced, approximately, so they can be left alone, deactivated. Now there are light-weight indices that have become heavy-weight with respect to the remaining active nodes, so we can continue balancing the active indices until the relatively heavy-weight among them become approximately balanced, and so forth. The problem with the naive solution is that balancing the active indices shifts the weights of both active and inactive indices, and they move out of their initial scale. If the scale sets of indices keep changing, it is hard to argue that the process converges. Shifting between scales is precisely what our algorithm and proof deal with. Light-weight indices that have become heavy-weight are easy to handle. They can keep being active. Heavy-weight indices that have become light-weight cannot continue to be inactive, because they are no longer guaranteed to be approximately balanced. Thus, in order to analyze convergence effectively, we need to bound the number (and global effect on weight) of these reactivation events.

2 Preliminaries

The input is a real square matrix $A=(a_{ij})_{n\times n}$. We denote the $i$-th row of such a matrix by $a_{i,.}$ and the $i$-th column by $a_{.,i}$. We also use the notation $[n]=\{1,2,\dots,n\}$. The matrix $A$ is balanced in the $L_{p}$ norm iff $\|a_{.,i}\|_{p}=\|a_{i,.}\|_{p}$ for every index $i\in[n]$. Since the condition for being balanced depends neither on the signs of the entries of $A$ nor on the diagonal values, we will assume without loss of generality that $A$ is non-negative with zeroes on the diagonal.

An invertible diagonal matrix $D=\mathtt{diag}(d_{1},\cdots,d_{n})$ balances $A$ in the $L_{p}$ norm iff $DAD^{-1}$ is balanced in the $L_{p}$ norm. A matrix $A$ is balanceable iff there exists an invertible diagonal matrix $D$ that balances $A$. Balancing a matrix $A=(a_{ij})_{n\times n}$ in the $L_{p}$ norm is equivalent to balancing the matrix $(a_{ij}^{p})_{n\times n}$ in the $L_{1}$ norm. Therefore, for the rest of the paper we focus on balancing matrices in the $L_{1}$ norm.

We use $a_{\min}$ to denote the minimum non-zero entry of $A$. We also define $w=\frac{1}{a_{\min}}\cdot\sum_{ij}a_{ij}$.

Definition 1**.**

Given $\epsilon>0$ and an $n\times n$ matrix $A$, we say that the index $i$ of $A$ (where $i\in[n]$) is $\epsilon$-balanced iff

[TABLE]

We say that $A$ is strictly $\epsilon$-balanced iff every index $i$ of $A$ is $\epsilon$-balanced.

Any implementation of Osborne’s iteration can be thought of as computing vectors $\mathbf{x^{(t)}}\in\mathbb{R}^{n}$ for $t=1,2,\dots$, where iteration $t$ is applied to the matrix $(a_{ij}^{(t)})=DAD^{-1}$ for $D=\mathtt{diag}(e^{x_{1}^{(t)}},e^{x_{2}^{(t)}},\dots,e^{x_{n}^{(t)}})$. Thus, for all $i,j$, $a_{ij}^{(t)}=a_{ij}\cdot e^{x_{i}^{(t)}-x_{j}^{(t)}}$. Initially, $\mathbf{x}^{(1)}=(0,0,\dots,0)$. A balancing step of the iteration chooses an index $i$, then sets $x_{i}^{(t+1)}=x_{i}^{(t)}+\frac{1}{2}\cdot\left(\ln\|a_{.,i}^{(t)}\|_{1}-\ln\|a_{i,.}^{(t)}\|_{1}\right)$, and for all $j\neq i$, keeps $x_{j}^{(t+1)}=x_{j}^{(t)}$. For $\mathbf{x}\in\mathbb{R}^{n}$, we denote the sum of entries of the matrix $DAD^{-1}$ for $D=\mathtt{diag}(e^{x_{1}},e^{x_{2}},\dots,e^{x_{n}})$ by $f(\mathbf{x})=f_{A}(\mathbf{x})=\sum_{ij}a_{ij}\cdot e^{x_{i}-x_{j}}$. For any $n\times n$ non-negative matrix $B=(b_{ij})$, we denote by $G_{B}$ the weighted directed graph with node set $\{1,2,\dots,n\}$, arc set $\{(i,j):b_{ij}>0\}$, where an arc $(i,j)$ has weight $b_{ij}$. We will assume henceforth that the undirected version of $G_{A}$ is connected, otherwise we can handle each connected component separately. We quote a few useful lemmas. The references contain the proofs.

Lemma 1** (Theorem 1 in Kalantari et al. [6]).**

The input matrix $A$ is balanceable if and only if $G_{A}$ is strongly connected. Moreover, $DAD^{-1}$ is balanced in the $L_{1}$ norm if and only if $D=\mathtt{diag}(e^{x_{1}^{*}},e^{x_{2}^{*}},\dots,e^{x_{n}^{*}})$, where $\mathbf{x}^{*}=(x_{1}^{*},x_{2}^{*},\dots,x_{n}^{*})$ minimizes $f(\mathbf{x})$ over $\mathbf{x}\in\mathbb{R}^{n}$.

Notice that $f$ is a convex function and the gradient $\nabla f(\mathbf{x})$ of $f$ at $\mathbf{x}$ is given by

[TABLE]

the difference between the total weight of arcs leaving node $i$ and the total weights of arcs going into node $i$ in the graph of $DAD^{-1}$ for $D=\mathtt{diag}(e^{x_{1}},e^{x_{2}},\dots,e^{x_{n}})$. If $DAD^{-1}$ is balanced then the arc weights $a_{ij}\cdot e^{x_{i}-x_{j}}$ form a valid circulation in the graph $G_{A}$, since the gradient has to be [math]. Some properties of $f$ are given in the following lemma.

Lemma 2** (Lemmas 2.1 and 2.2 in Ostrovsky et al. [9]).**

If $\mathbf{x}^{\prime}$ is derived from $\mathbf{x}$ by balancing index $i$ of a matrix $B=(b_{ij})_{n\times n}$, then $f(\mathbf{x})-f(\mathbf{x}^{\prime})=(\sqrt{\|b_{.,i}\|_{1}}-\sqrt{\|b_{i,.}\|_{1}})^{2}$. Also, for all $\mathbf{x}\in\mathbb{R}^{n}$, $f(\mathbf{x})-f(\mathbf{x}^{*})\leq\frac{n}{2}\cdot\|\nabla f(\mathbf{x})\|_{1}$.

We also need the following absolute bounds on the arc weights.

Lemma 3** (Lemma 3.2 in Ostrovsky et al. [9]).**

Suppose that a matrix $B$ is derived from a matrix $A$ through a sequence of balancing operations. Then, for every arc $(i,j)$ of $G_{B}$, $\left(\frac{a_{\min}}{\sum_{ij}a_{ij}}\right)^{n}\cdot\sum_{ij}a_{ij}\leq b_{ij}\leq\sum_{ij}a_{ij}$. (Notice that the arcs of $G_{B}$ are identical to the arcs of $G_{A}$.)

Finally, we prove the following global condition on indices being $\epsilon$-balanced.

Lemma 4**.**

Consider a matrix $B=DAD^{-1}=(b_{ij})_{n\times n}$, where $D=\mathtt{diag}(e^{x_{1}},e^{x_{2}},\dots,e^{x_{n}})$, that was derived from $A$ by a sequence of zero or more balancing operations. Let $\epsilon\in(0,1/2]$, and put $\epsilon^{\prime}=\frac{\epsilon^{2}}{64n^{4}}$. Suppose that $\|\nabla f_{A}(\vec{0})\|_{1}\leq\epsilon^{\prime}\cdot f_{A}(\vec{0})$. Then, for every $i\in[n]$ we have the following implication. If $\|b_{.,i}\|_{1}+\|b_{i,.}\|_{1}\geq\frac{1}{8n^{3}}\cdot f_{A}(\mathbf{x})$, then index $i$ is $\epsilon$-balanced in $B$.

Proof.

We will show the contrapositive claim that if a node is not $\epsilon$-balanced then it must have low weight (both with respect to $B$). Let $i$ be an index that is not $\epsilon$-balanced in $B$. Without loss of generality we may assume that the in-weight is larger than the out-weight, so ${\|b_{.,i}\|_{1}}/{\|b_{i,.}\|_{1}}>1+\epsilon$. Consider what would happen if we balance index $i$ in $B$, yielding a vector $\mathbf{x}^{\prime}$ that differs from $\mathbf{x}$ only in the $i$-th coordinate.

[TABLE]

where the equation follows from Lemma 2 and the last inequality uses the fact that $\epsilon\leq\frac{1}{2}$. On the other hand, we have

[TABLE]

where the first inequality follows from the the fact that every balancing step decreases $f_{A}$, the second inequality follows from Lemma 2, the third inequality follows from the assumption on $f_{A}(\vec{0})$, and the last equation follows from the choice of $\epsilon^{\prime}$. Combining the bounds on $f_{A}(\mathbf{x})-f_{A}(\mathbf{x}^{\prime})$ in Equations (1) and (2) gives

[TABLE]

and this completes the proof. ∎

3 Strict Balancing

In this section we present a variant of Osborne’s iteration and prove that it converges in polynomial time to a strictly $\epsilon$-balanced matrix. The algorithm, a procedure named StrictBalance, is defined in pseudocode labeled Algorithm 1 on page 1. Lemma 4 above motivates the main idea of contracting heavy nodes in step 14 of StrictBalance.

Our main theorem is

Theorem 1**.**

StrictBalance($A$, $\epsilon$) returns a strictly $\epsilon$-balanced matrix $B=DAD^{-1}$ after at most

[TABLE]

balancing steps, using $O\left(\epsilon^{-2}n^{10}\log(wn/\epsilon)\log w\right)$ arithmetic operations over $O(n\log w)$-bit numbers.

The proof of Theorem 1 uses a few arguments, given in the following lemmas. A phase of StrictBalance is one iteration of the outer while loop. Notice that in the beginning of this loop the variable $s$ indexes the phase number (i.e., $s-1$ phases were completed thus far). Also in the beginning of the inner while loop the variable $t$ indexes the total iteration number from all phases (i.e., $t-1$ balancing operations from all phases were completed thus far).

We identify outer loop iteration $s$ with an interval $[t_{s},t_{s+1})=\{t_{s},t_{s}+1,\dots,t_{s+1}-1\}$ of the inner loop iterations executed during phase $s$. We denote by ${\cal B}_{s,t}$ the value of ${\cal B}_{s}$ in the beginning of the inner while loop iteration number $t$ (dubbed time $t$). If $t\in[t_{j},t_{j+1})$, then ${\cal B}_{s,t}$ is defined only for $s\leq j$. We also use $G^{(\mathcal{B}_{s,t})}$ to denote the graph that is obtained by contracting the nodes of set ${\cal B}_{s,t}$ in graph $G_{A}$. Also $f^{(\mathcal{B}_{s,t})}$ is the function corresponding to graph $G^{(\mathcal{B}_{s,t})}$ and $f^{(\mathcal{B}_{s,t})}(\mathbf{x^{(t)}})$ denotes the sum of weights of arcs of graph $G^{(\mathcal{B}_{s,t})}$ at time $t$. If set $B_{s}$ is unchanged during an interval and there is no confusion, we may use $G^{(\mathcal{B}_{s})}$ instead of $G^{(\mathcal{B}_{s,t})}$. Particularly we use $f^{(\mathcal{B}_{s})}(\mathbf{x^{(t)}})$ instead of $f^{(\mathcal{B}_{s,t})}(\mathbf{x^{(t)}})$. We refer to the quantity $\|a_{.,i}^{(t)}\|_{1}+\|a_{i,.}^{(t)}\|_{1}$ as the weight of node $i$ at time $t$.

Lemma 5**.**

For every phase $s\geq 1$, for every $t\geq t_{s+1}$, $\mathcal{B}_{s,t}=\mathcal{B}_{s,t_{s+1}}$.

Proof.

The claim follows easily from the fact that any iteration $t\geq t_{s+1}$ belongs to a phase $s^{\prime}>s$, so $\mathcal{B}_{s,t_{s+1}}\cap(\mathcal{B}_{s^{\prime},t}\setminus\mathcal{B}_{s^{\prime}-1,t})=\emptyset$, and by line 8 and 9 of StrictBalance none of the nodes in $\mathcal{B}_{s,t_{s+1}}$ will be removed. ∎

Lemma 6**.**

For all $s>1$, for all $t\in[t_{s},t_{s+1})$, $f^{(\mathcal{B}_{s,t})}(\mathbf{x^{(t)}})\leq\left(n-\left|{\mathcal{B}}_{s,t}\right|\right)\cdot\tau_{s}$.

Proof.

Let $t_{s}=t_{s,1}<t_{s,2}<t_{s,3}<\cdots<t_{s,\ell_{s}}$ denote the time steps before which $\mathcal{B}_{s}$ changes during phase $s$. For simplicity, we abuse notation and use $\mathcal{B}_{s,j}$ instead of $\mathcal{B}_{s,t_{s,j}}$. Clearly $\mathcal{B}_{s,1}\supseteq\mathcal{B}_{s,2}\ldots\supseteq\mathcal{B}_{s,\ell_{s}}$, because we only remove nodes from $\mathcal{B}_{s}$ once it is set. Fix $s>1$. We prove this lemma by induction on $r\in\{1,2,\dots,\ell_{s}\}$. For the basis, let $r=1$. Clearly, by the way the algorithm sets $\mathcal{B}_{s}$ before time $t_{s,1}$, all nodes with weight $\geq\tau_{s}$ are in $\mathcal{B}_{s}$, and therefore every node $i\not\in\mathcal{B}_{s}$ has weight at most $\tau_{s}$, so the lemma follows. Now, assume that the lemma is true for every $t\leq t_{s,r}$, we show that the lemma holds for every $t\leq t_{s,r+1}$. If $t\in[t_{s,r},t_{s,r+1})$, then $\mathcal{B}_{s,t}=\mathcal{B}_{s,t_{s,r}}$, and we have:

[TABLE]

The first inequality holds because balancing operations from time $t_{s,r}$ to time $t$ only reduce the value of $f^{(\mathcal{B}_{s})}$, and the second inequality holds by the induction hypothesis.

Just before iteration $t=t_{s,r+1}$, the set $\mathcal{B}_{s}$ changes, and one or more nodes are removed from it. However, every removed node has weight at most $\tau_{s}$, and its removal does not change the weights of the other nodes in $[n]\setminus\mathcal{B}_{s}$. Therefore, if $k$ nodes are removed from $\mathcal{B}_{s}$,

[TABLE]

This completes the proof. ∎

Corollary 1**.**

For all $s>1$, $f^{(\mathcal{B}_{s})}(\mathbf{x^{(t_{s+1})}})\leq\frac{1}{4n^{2}}\cdot f^{(\mathcal{B}_{s-1})}(\mathbf{x^{(t_{s})}})$. If $s>2$, then $\tau_{s}\leq\frac{\tau_{s-1}}{4n^{2}}$.

Proof.

Notice that

[TABLE]

where the inequality follows from Lemma 6, and the equation follows from line 13 of StrictBalance. This proves the first assertion. As for the second assertion, notice that if $s>2$ then $s-1>1$, so using line 13 of StrictBalance and Lemma 6 again,

[TABLE]

as stipulated. ∎

Lemma 7**.**

For every phase $s>1$, for every $t\geq t_{s}$, all the nodes in $\mathcal{B}_{s,t}$ have weight $\geq\tau_{s}/2$ and are $\epsilon$-balanced at time $t$.

Proof.

Fix $s>1$ and let $i\in\mathcal{B}_{s,t}$. Without loss of generality $i\not\in\mathcal{B}_{s-1,t}$, otherwise we can replace $s$ with $s-1$. (Recall that $\mathcal{B}_{1}=\emptyset$ at all times.) Also note that it must be the case that $i\in\mathcal{B}_{s,t_{s}}$, because $\mathcal{B}_{s}$ does not accumulate additional nodes after being created. If $t\in[t_{s},t_{s+1}]$, then lines 13-14 and 8-9 of StrictBalance guarantee that if $i\in\mathcal{B}_{s,t}\setminus\mathcal{B}_{s-1,t}$, then its weight at time $t$ is at least $\tau_{s}$.

Otherwise, consider $t>t_{s+1}$ and let $s^{\prime}>s$ be the phase containing $t$. Consider a phase $j>s$. By Lemma 6 the total weight of $f^{(\mathcal{B}_{j})}$ during phase $j$ is at most $n\tau_{j}$, and $f^{(\mathcal{B}_{j})}$ never drops below 0. So, the total weight that a node $i\in\mathcal{B}_{j}$ can lose (which is at most the total weight that $f^{(\mathcal{B}_{j})}$ can lose) is at most $n\tau_{j}$. By Corollary 1, for every $j>s$, $\tau_{j+1}\leq\frac{\tau_{j}}{4n^{2}}$. Now, suppose that $t$ is an iteration in phase $s^{\prime}>s$. Then, the weight of $i$ at time $t$ is at least

[TABLE]

Thus we have established that at any time $t\geq t_{s}$, if $i\in\mathcal{B}_{s,t}$ then its weight is at least $\frac{\tau_{s}}{2}=\frac{1}{8n^{3}}f^{(\mathcal{B}_{s-1})}(\mathbf{x^{(t_{s})}})$. By line 4 of StrictBalance, $\|\nabla f^{(\mathcal{B}_{s-1,t_{s}})}(\mathbf{x^{(t_{s})}})\|_{1}\leq\epsilon^{\prime}\cdot f^{(\mathcal{B}_{s-1,t_{s}})}(\mathbf{x^{(t_{s})}})$. By Lemma 5, $\mathcal{B}_{s-1}$ does not change in the interval $[t_{s},t]$. Therefore, we conclude from Lemma 4 that $i$ is $\epsilon$-balanced at time $t$. ∎

Lemma 8**.**

Suppose that $t<t^{\prime}$ satisfies $[t,t^{\prime})\subseteq[t_{s},t_{s+1})$, and furthermore, during the iterations in the interval $[t,t^{\prime})$ the set ${\cal B}_{s}$ does not change (it could change after balancing step $t^{\prime}-1$). Then, the length of the interval

[TABLE]

Proof.

Rename the nodes so that $\mathcal{B}_{s,t}=\mathcal{B}_{s,t^{\prime}-1}=\{p,p+1,\ldots,n\}$. The assumption that $\mathcal{B}_{s}$ does not change during the interval $[t,t^{\prime})$ means that the weights of all the nodes $p,p+1,\dots,n$ remain at least $\tau_{s}$ for the duration of this interval. During the interval $[t,t^{\prime})$, the graph $G^{(\mathcal{B}_{s})}$ (which remains fixed) is obtained by contracting the nodes $p,p+1,\ldots,n$ in $G_{A}$. So $G^{(\mathcal{B}_{s})}$ has $p$ nodes $1,2,\dots,p-1,p$, where the last node $p$ is the contracted node. In each iteration in the interval $[t,t^{\prime})$, one of the nodes $1,2,\dots,p-1$ is balanced. Consider some time step $t^{\prime\prime}\in[t,t^{\prime})$, and let $I_{i}$ and $O_{i}$, respectively, denote the current sums of weights of the arcs of $G^{(\mathcal{B}_{s})}$ into and out of node $i$, respectively. Let $j\in[p-1]$ be the node that maximizes $\frac{(I_{j}-O_{j})^{2}}{I_{j}+O_{j}}$. We have

[TABLE]

The first equation follows from the choice of $i$ in line 5 StrictBalance, and Lemma 2. The third inequality follows from an averaging argument and the choice of $j$. The fourth inequality uses Cauchy-Schwarz. The last inequality holds because $\sum_{i=1}^{p}(I_{i}-O_{i})=0$, so $|I_{p}-O_{p}|=\left|\sum_{i=1}^{p-1}(I_{i}-O_{i})\right|\leq\sum_{i=1}^{p-1}|I_{i}-O_{i}|$, and therefore $\sum_{i=1}^{p}|I_{i}-O_{i}|\leq 2\sum_{i=1}^{p-1}|I_{i}-O_{i}|$.

Since the interval $[t,t^{\prime})$ is contained in phase $s$, the stopping condition for the phase does not hold, so

[TABLE]

Therefore,

[TABLE]

where the last inequality follows from Lemma 2. Rearranging the terms gives

[TABLE]

Iterating for $T$ step yields

[TABLE]

Now, by Lemma 3, we have that $f^{(\mathcal{B}_{s})}(\mathbf{x^{(t)}})-f^{(\mathcal{B}_{s})}(\mathbf{x}^{*})\leq f^{(\mathcal{B}_{s})}(\mathbf{x^{(t)}})\leq\sum_{i,j=1}^{n}a_{ij}$, and for all $t^{\prime\prime}$, $f^{(\mathcal{B}_{s})}(\mathbf{x^{(t^{\prime\prime})}})\geq\frac{1}{w^{n}}\sum_{i,j=1}^{n}a_{ij}$. Therefore, if $t^{\prime}-t\geq\frac{8n^{2}}{\epsilon^{\prime}}\cdot\ln\left(16nw^{n}/(\epsilon^{\prime})^{2}\right)+1$, then

[TABLE]

Therefore,

[TABLE]

where the first inequality follows from (3), and the second inequality holds because $f^{(\mathcal{B}_{s})}(\mathbf{x}^{*})\leq f^{(\mathcal{B}_{s})}(\mathbf{x^{(t^{\prime}-1)}})$. We get that $\frac{\|\nabla f^{(\mathcal{B})}(\mathbf{x^{(t^{\prime}-1)}})\|_{1}}{f^{(\mathcal{B})}(\mathbf{x^{(t^{\prime}-1)}})}\leq\epsilon^{\prime}$, in contradiction to our assumption that the phase does not end before the start of iteration $t^{\prime}$. ∎

Corollary 2**.**

In any phase, the number of balancing steps is at most $O\left(\epsilon^{-2}n^{8}\log(wn/\epsilon)\right)$.

Proof.

In the beginning of phase $s$ the set $\mathcal{B}_{s}$ contains at most $n-1$ nodes. Partition the phase into intervals $[t,t^{\prime})$ where $\mathcal{B}_{s}$ does not change during an interval, but does change between intervals. By Lemma 8, each interval consists of at most $O\left(\epsilon^{-2}n^{7}\log(wn/\epsilon)\right)$ balancing steps. Since nodes that are removed from $\mathcal{B}_{s}$ between intervals are never returned to $\mathcal{B}_{s}$, the number of such intervals is at most $n-1$. Hence, the total number of balancing steps in the phase is at most $O\left(\epsilon^{-2}n^{8}\log(wn/\epsilon)\right)$. ∎

Lemma 9**.**

The total number of phases of the algorithm is $O(n\log w/\log n)$.

Proof.

Let $s>2$ be a phase of the algorithm and $t\in[t_{s},t_{s+1})$. By Lemma 3, $f^{(\mathcal{B}_{s,t})}(\mathbf{x^{(t)}})\geq\frac{1}{w^{n}}\cdot\sum_{ij}a_{ij}$. On the other hand, by Lemma 6 and Corollary 1, $\tau_{s}\leq\frac{1}{(4n^{2})^{s-2}}\cdot\tau_{2}\leq\frac{1}{(4n^{2})^{s-2}}\cdot\sum_{ij}a_{ij}$, and $f^{(\mathcal{B}_{s,t})}(\mathbf{x^{(t)}})\leq n\tau_{s}$. Combining these gives $\frac{1}{w^{n}}\cdot\sum_{ij}a_{ij}\leq n\tau_{s}\leq\frac{n}{(4n^{2})^{s-2}}\cdot\sum_{ij}a_{ij}$ which implies that $s\leq\frac{\log(nw^{n})}{\log(4n^{2})}+2$. ∎

Proof of Theorem 1.

By Lemma 9, for some $s=O(n\log w/\log n)$, StrictBalance terminates, so $\mathcal{B}_{s,t_{s}}=[n]$. By Corollary 2, the number of balancing steps in a phase is at most $O\left(\epsilon^{-2}n^{8}\log(wn/\epsilon)\right)$. Therefore, the total number of balancing steps is at most $O\left(\epsilon^{-2}n^{9}\log(wn/\epsilon)\log w/\log n\right)$. These balancing steps require at most $O\left(\epsilon^{-2}n^{10}\log(wn/\epsilon)\log w\right)$ arithmetic operations over $O(n\log w)$-bit numbers. When the algorithm terminates at time $t_{s}$, all the nodes are in $\mathcal{B}_{s,t_{s}}$, and by Lemma 7 they are all $\epsilon$-balanced, so the matrix is strictly $\epsilon$-balanced. ∎