On Solving Linear Systems in Sublinear Time

We first present our simpler algorithm for linear systems in Laplacians with a bounded condition number.

A few extensions of the theorem follow easily from our proof. First, if the algorithm is given also an upper bound $B_{up}$ on $||b||_{0}$, then the expression for $s$ can be refined by replacing $n$ with $B_{up}\leq n$. Second, we can improve the runtime to $O(\epsilon^{-2}s^{3}\log s)$ whenever the representation of $G$ allows to sample a uniformly random neighbor of a vertex in constant time. Third, the algorithm has an (essentially) cubic dependence on the condition number $\kappa(L_{G})$, which can be improved to quadratic if we allow a preprocessing of $G$ (or, equivalently if we only count the number of probes into $b$). Later we show that this quadratic dependence is near-optimal.

We further design an algorithm for SDD matrices with bounded condition number. The formal statement appears in Theorem 5.1 and is a natural generalization of Theorem 1.1 with two differences. One difference is that a natural solution to the system $Sx=b$ is $x=S^{+}b$, but our method requires $S$ to have normalized diagonal entries, and thus we aim at another solution $x^{*}$, construed as follows. Define

then our linear system can be written as $\tilde{S}(D^{1/2}x)=D^{-1/2}b$, which has a solution

that is expressed using the pseudo-inverse of $\tilde{S}$, rather than of $S$.

A second difference is that Theorem 5.1 makes no assumptions about the multiplicity of the eigenvalue $0$ of $\tilde{S}$, e.g., if $S$ is a graph Laplacian, then the graph need not be connected. The assumptions needed to achieve a polylogarithmic time, beyond $\tilde{S}$ having a bounded condition number, are only that a random “neighbor” in the graph corresponding to $S$ can be sampled quickly, and that $\frac{\max_{i\in[n]}D_{ii}}{\min_{i\in[n]}D_{ii}}\leq\mbox{poly}(n)$, which holds if $S$ has polynomially-bounded entries.

Our first lower bound shows that the above guarantees cannot be obtained for a general PSD matrix, even if we are allowed to preprocess the matrix $S$, and only count probes into $b$. The proof employs a PSD matrix $S$ that is invertible (i.e., positive definite), in which case the linear system $Sx=b$ has a unique solution $x=S^{-1}b$.

The second lower bound shows that our SDD algorithm has a near-optimal dependence on the condition number of $S$, even if we are allowed to preprocess the matrix S, and only count probes into $b$. The lower bound holds even for Laplacian matrices.

An example application of our algorithmic results is computing the effective resistance between a pair of vertices $u,v$ in a graph $G$ (given $u$,$v$ and $G$ as input). It is well known that the effective resistance, denoted $\Reff(u,v)$, can be expressed as $x_{u}-x_{v}$, where $x$ solves $L_{G}x=e_{u}-e_{v}$. The spectral-sparsification algorithm of Spielman and Srivastava [SS11] relies on a near-linear time algorithm (that they devise) for approximating the effective resistances of all edges in $G$. For unweighted graphs, there is also a faster algorithm [Lee14] that runs in time $\tilde{O}(n)$, which is sublinear in the number of edges, and approximates effective resistances within a larger factor $\polylog(n)$. In a $d$-regular expander $G$, the effective resistance between every two vertices is $\Theta(1/d)$, and our algorithm in Theorem 1.1 can quickly compute an arbitrarily good approximation (factor $1+\epsilon$). Indeed, observe that we can use $B_{up}=2$, hence the runtime is $O(\frac{1}{\epsilon^{2}}\polylog\frac{1}{\epsilon})$, independently of $n$. The additive accuracy is $\epsilon\cdot\max_{ij\in E(G)}|x_{i}-x_{j}|$; in fact, each $x_{i}-x_{j}$ is the potential difference between $i$ and $j$ when a potential difference of $\Reff(u,v)$ is imposed between $u$ and $v$, thus $\max_{ij\in E(G)}|{x_{i}-x_{j}}|\leq{x_{u}-x_{v}}=\Reff(u,v)$, and hence with high probability the output actually achieves a multiplicative guarantee $\hatReff(u,v)\in(1\pm\epsilon)\Reff(u,v)$.

Our basic technique relies on a classic idea of von Neumann and Ulam [FL50, Was52] for estimating a matrix inverse by a power series; see Section 1.3 for a discussion of related work. Our starting point is the identity

(Recall that $||X||$ denotes the spectral norm of a matrix $X$.) Now given a Laplacian $L=L_{G}$ of a $d$-regular graph $G$, observe that $\frac{1}{d}L=I-\frac{1}{d}A$, where $A$ is the adjacency matrix of $G$. Assume for a moment that $||\frac{1}{d}A||<1$; then by the above identity, $(\frac{1}{d}L)^{-1}=(I-\frac{1}{d}A)^{-1}=\sum_{t=0}^{\infty}(\frac{1}{d}A)^{t}$, and the solution of the linear system $Lx=b$ would be $x^{*}=L^{-1}b=\frac{1}{d}\sum_{t=0}^{\infty}(\frac{1}{d}A)^{t}b$. The point is that the summands decay exponentially because $||(\frac{1}{d}A)^{t}b||_{2}\leq||(\frac{1}{d}A)^{t}||\cdot||b||_{2}\leq||(\frac{1}{d}A)||^{t}\cdot||b||_{2}$. Therefore, we can estimate $x^{*}_{u}$ using the first $t_{0}$ terms, i.e., $\hat{x}_{u}=e_{u}^{\mathsf{T}}\frac{1}{d}\sum_{t=0}^{t_{0}}(\frac{1}{d}A)^{t}b$, where $t_{0}$ is logarithmic (with base $\|{\frac{1}{d}A}\|^{-1}>1$). In order to compute each term $e_{u}^{\mathsf{T}}\frac{1}{d}(\frac{1}{d}A)^{t}b$, observe that $e_{u}^{\mathsf{T}}(\frac{1}{d}A)^{t}e_{w}$ is exactly the probability that a random walk of length $t$ starting at $u$ will end at vertex $w$. Thus, if we perform a random walk of length $t$ starting at $u$, and let $z$ be its (random) end vertex, then

If we perform several random walks (specifically, $\poly(t_{0},\frac{1}{\epsilon})$ walk suffice), average the resulting $b_{z}$’s, and then multiply by $\frac{1}{d}$, then with high probability, we will obtain a good approximation to $e_{u}^{\mathsf{T}}\frac{1}{d}(\frac{1}{d}A)^{t}b$.

As a matter of fact, we have a non-strict inequality $||\frac{1}{d}A||\leq 1$, because of the all-ones vector $\vec{1}\in\mathbb{R}^{n}$. Nevertheless, we can still get a meaningful result if all eigenvalues of $A$ except for the largest one are smaller than $d$ (equivalently, the graph $G$ is connected). First, we get rid of any negative eigenvalues by the standard trick of considering $(dI+A)/2$ instead of $A$, which is equivalent to adding $d$ self-loops at every vertex. We may assume $b$ is orthogonal to $\vec{1}$ (otherwise the linear system has no solution), hence the linear system $Lx=b$ has infinitely many solutions, and since $L$ is PSD, we can aim to estimate the specific solution $x^{*}\stackrel{{\scriptstyle{\rm def}}}{{=}}L^{+}b$ by $\frac{1}{d}\sum_{t=0}^{t_{0}}(\frac{1}{d}A)^{t}b$. Indeed, the idealized analysis above still applies by restricting all our calculations to the subspace orthogonal to $\vec{1}$. This is carried out in Theorem 1.1.

To generalize the above approach to SDD matrices, we face three issues. First, due to the irregularity of general SDD matrices, it is harder to properly define the equivalent random walk matrix. We resolve this by normalizing the SDD matrix $S$ into $\tilde{S}$ defined in (2), and solving the equivalent (normalized) system $\tilde{S}(D^{1/2}x)=D^{-1/2}b$. Second, general SDD matrices can have positive off-diagonal elements, in constrast to the Laplacians. To address this, we interpret such entries as negative-weight edges, and employ random walks that “remember” the signs of the traversed edges. Third, diagonal elements may strictly dominate their row, which we address by terminating the random walk early with some positive probability.

We first discuss our lower bound for PSD matrices, which is one of the main contributions of our work. The goal is to exhibit a matrix $S$ for which estimating a coordinate $x^{*}_{u}$ of the solution $x^{*}=S^{-1}b$ requires $n^{\Omega(1)}$ probes into the input $b$.

Without the sparsity constraint on $S$, one can deduce such a lower bound via a reduction to the communication complexity of the Vector in Subspace Problem (VSP), in which Alice has an $n/2$-dimensional subspace $H\subset\mathbb{R}^{n}$, Bob has a vector $b\in\mathbb{R}^{n}$, and their goal is to determine whether $b\in H$ or $b\in H^{\perp}$. The randomized communication complexity of this promise problem is between $\Omega(n^{1/3})$ [KR11] and $O(\sqrt{n})$ [Raz99] (while for quantum communication it is $O(\log n)$). To reduce this problem to linear-system solvers, let $P_{H}\in R^{n\times n}$ be the projection operator onto the subspace $H$, and set $S=I+P_{H}$. Consider the system $Sx=b$, and note that Alice knows $S$ and Bob has $b$. It is easy to see that the unique solution $x^{*}$ is either $b$ or $\tfrac{1}{2}b$, depending on whether $b\in H^{\perp}$ or $b\in H$. Alice and Bob could use a solver that makes few probes to $b$, as follows. Bob would pick an index $u\in[n]$ that maximizes $|{b_{u}}|$ (and thus also $|{x_{u}}|$), and send it to Alice. She would then apply the solver, receiving from Bob only a few entries of $b$, to estimate $x_{u}$ within additive error $\tfrac{1}{2}\|{x}\|_{\infty}$, which suffices to distinguish the two cases. This matrix $S$ is PSD with condition number $\kappa(S)\leq 2$. However it is dense.

We thus revert to a different approach of proving it from basic principles. Our high-level idea is to take a $2d$-regular expander and assign its edges with signs ($\pm 1$) that are random but balanced everywhere (namely, at every vertex, the incident edges are split evenly between positive and negative). The signed adjacency matrix $A\in\set{-1,0,+1}^{n\times n}$ should have spectral norm $\mu\stackrel{{\scriptstyle{\rm def}}}{{=}}\|{A}\|=O(\sqrt{d})$, and then instead of the (signed) Laplacian $L=(2d)I-A$, we consider $S=2\mu I-A$, which is PSD with condition number $\kappa(S)\leq 3$, as well as invertible and sparse. Now following similar arguments as in our algorithm, we can write $S^{-1}$ as a power series of the matrix $A$, and express coordinate $x^{*}_{u}$ of the solution $x^{*}=S^{-1}b$ via $\EX_{z}[b_{z}]$ where $z$ is the (random) end vertex of a random walk that starts at $u$ and its length is bounded by some $t_{0}$ (performed in the “signed” graph corresponding to $A$). Now if the graph around $u$ looks like a tree (e.g., it has high girth), then not-too-long walks are highly symmetric and easy to count. We now let $b_{v}$ be non-zero only at vertices $v$ at distance exactly $t_{0}$ from $u$, and for these vertices $b_{v}$ is set to $+1$ or $-1$ at random but with a small bias $\delta$ towards one of the values. Some calculations show that $\sgn(\EX_{z}[b_{z}])$, and consequently $\sgn(x^{*}_{u})$, will be according to our bias (with high probability), however discovering this $\sgn(x^{*}_{u})$ via probes to $b$ is essentially the problem of learning a biased coin, which requires $\Omega(\delta^{-2})$ coin observations. An additional technical obstacle is to prove that the solution $x^{*}$ has a small $\ell_{\infty}$-norm, so that we can argue that an $\tfrac{1}{5}\|{x^{*}}\|_{\infty}$-additive error to $x^{*}_{u}$ will not change its sign. Overall, we show we can set $t_{0}=\Omega(\log n)$ and $\delta\approx((2d-1)^{t_{0}})^{-1/2}$, thus concluding that the algorithm must observe $\Omega(\delta^{-2})=n^{\Omega(1)}$ entries of $b$.

It is instructive to ask where in the above argument is it crucial to have $\mu=O(\sqrt{d})$, because if it were also valid for $\mu=d$, in which case the matrix $S=2\mu I-A$ is SDD, then it would contradict our own algorithm for SDD matrices. The answer is that $\mu\ll 2d$ is required to bound $\|{x^{*}}\|_{\infty}$, specifically in the analysis that follows immediately after Eqn. (10).

We now outline the ideas to prove the $\tilde{\Omega}(\kappa^{2})$ lower bound even for Laplacian systems with condition number $\kappa$. First we note that it is relatively straight-forward to prove that a linear dependence on the condition number is necessary. Namely, consider a dumbbell graph (two 3-regular expanders connected by a bridge edge $(u,v)$), for which we need to estimate $x^{*}_{u}-x^{*}_{v}$. For $b$ defined as $b=e_{i}-e_{j}$, the value of $x^{*}_{u}-x^{*}_{v}$ will be non-zero iff vertices $i,j$ are on the opposite sides of the bridge. To determine the latter, one requires $\Omega(n)$ queries into $b$. Since this graph has a condition number of $O(n)$, we obtain an $\Omega(\kappa)$ lower bound.

The quadratic lower bound requires both a different graph and a different distribution over $b$. We use the following graph $G$ with condition number $O(k)$: take two 3-regular expanders and connect them with $n/k$ edges (“bridge edges”). The vector $b\in\{-1,+1\}^{n}$ will be dense and in particular it is either: 1) balanced, i.e., $\sum b_{i}$ on each expander is zero, or 2) unbalanced, i.e., each $b_{i}$ is chosen $\pm 1$ with a bias $p\approx 1/k$ towards $+1$ on the first expander, and towards $-1$ on the second one. Now, as above, it is simple to prove that on average over the bridge edges $(u,v)$: 1) in the balanced case, the average of $x^{*}_{u}-x^{*}_{v}$ must be zero, and 2) in the unbalanced case, the average must be $\Omega(1)$. However, the main challenge is that the actual values may differ from the average — e.g., even in the balanced case, each bridge edge $(u,v)$ will likely have non-zero value of $x^{*}_{u}-x^{*}_{v}$. Nonetheless, we manage to prove an upper bound on the maximum value of $|x^{*}_{u}-x^{*}_{v}|$ over all edges $(u,v)$ (as in the previous lower bound, we need to bound $\|x^{*}\|_{\infty}$ as well). For the latter, we need to again analyze $\E_{z}[b_{z}]$ where $z$ is the endpoint of a random walk of some fixed length $i\geq 1$ starting from $u$ in the graph $G$. Since the vector $b$ is not symmetric over the graph $G$, a direct analysis seems hard — instead we estimate $\E_{z}[b_{z}]$ via a coupling of such walks in $G$ with random walks in an expander, which is amenable to a direct analysis.

An example application of our results is computing quickly the PageRank (defined in [BP98]) of a single node in an undirected $d$-regular graph. Recall that the PageRank vector of an $n$-vertex graph with associated transition matrix $P$ is the solution to the linear system $x=\frac{1-\alpha}{n}\vec{1}+\alpha Px$, where $0<\alpha<1$ is a given parameter. In personalized PageRank, one replaces $\frac{1}{n}\vec{1}$ (the uniform distribution) with some $b\in\mathbb{R}^{n}$, e.g., a standard basis vector. Equivalently, $x$ solves the system $Sx=\frac{1-\alpha}{n}\vec{1}$ where $S=I-\alpha P$ is an SDD matrix with $1$’s on the diagonal. As all eigenvalues of $P$ are of magnitude at most 1 (recall $P$ is a transition matrix), all eigenvalues of $I-\tilde{S}=I-S=\alpha P$ are of magnitude at most $\alpha$, and the running time guaranteed by Theorem 5.1 is logarithmic (with base $\frac{2}{\alpha+1}$).

Algorithms for the PageRank model were studied extensively. In particular, the sublinear algorithms of [BPP18] approximate the PageRank of a vertex using $\tilde{O}(n^{2/3})$ queries to an arbitrary graph, or using $\tilde{O}((n\Delta)^{1/2})$ queries when the maximum degree is $\Delta$. Another example is the heavy-hitters algorithm of [BBCT14], which reports all vertices whose approximate PageRank exceeds a threshold $T$ in sublinear time $\tilde{O}(1/\Delta)$, when PageRanks are viewed as probabilities and sum to $1$. Other work explores connections to other graph problems, including for instances using PageRank algorithms to approximate effective resistances [CZ10], the PageRank vector itself, and computing sparse cuts [ACL07].

Our algorithms in Theorems 1.1 and 5.1 are local in the sense that they query a small portion of their input, usually around the input vertex, when viewed as graph algorithms. Local algorithms for graph problems were studied in several contexts, like graph partitioning [ST04, AP09], Web analysis [CGS04, ABC+08], and distributed computing [Suo13]. Rubinfeld, Tamir, Vardi, and Xie [RTVX11] introduced a formal concept of Local Computation Algorithms that requires consistency between the local outputs of multiple executions (namely, these local outputs must all agree with a single global solution). As explained earlier, our problem formulation (1) follows this consistency requirement.