On Solving Asymmetric Diagonally Dominant Linear Systems in Sublinear Time

For $n\in {\mathbb{Z}}^{+}$, we define $[n]:=\{1,2,\mathrm{\dots },n\}$. We call a matrix $𝐌\in {ℝ}^{n\times n}$ RDD (row diagonally dominant) if it satisfies $𝐌(j,j)\ge {\sum }_{k\ne j}\left|𝐌(j,k)\right|$ for all $j\in [n]$ and its diagonal entries are positive. We call a matrix CDD (column diagonally dominant) if its transpose is RDD, and call a matrix SDD (symmetric diagonally dominant) if it is symmetric and RDD. It is well-known that any SDD matrix is PSD (positive semidefinite). We call a square matrix Z-matrix if its off-diagonal entries are nonpositive. We use ${𝒆}_k$ to denote the $k$-th canonical unit vector and $\mathrm{𝟏}$ to denote the all-one vector. For any matrix or vector, we use $|\cdot |$ to denote taking the entrywise absolute value.

We call two real numbers $p,q>1$ Hölder conjugates (or $q$ is conjugate to $p$) if they satisfy $1/p+1/q=1$. By convention, we also formally let $1/\mathrm{\infty }:=0$ and view $\mathrm{\infty }$ and $1$ as Hölder conjugates. For any $p\in [1,\mathrm{\infty }]$, we use ${\Vert 𝒙\Vert }_p$ to denote the $p$-norm of a vector $𝒙$, and ${\Vert 𝐌\Vert }_p$ to denote the matrix norm induced by vector $p$-norm of a matrix $𝐌\in {ℝ}^{n\times n}$.

We consider a linear system $𝐌𝒙=𝒃$, where $𝐌\in {ℝ}^{n\times n}$ is an RDD/CDD matrix and $𝒃\in \mathrm{range}(𝐌)$, and a coefficient vector $𝒕\in {ℝ}^n$. We assume that $𝒃,𝒕\ne \mathrm{𝟎}$ and all nonzero entries in $𝐌$, $𝒃$, and $𝒕$ have absolute values in $[1/\mathrm{poly}(n),\mathrm{poly}(n)]$. We assume that the algorithms are given the dimension $n$ and have oracle access to $𝐌$, $𝒃$, and $𝒕$ via the following basic queries:

• $\mathrm{■}$

  Diagonal queries for $𝐌$: return $𝐌(k,k)$ in $O(1)$ time for a given index $k\in [n]$;

• $\mathrm{■}$

  Row/column queries for $𝐌$: return the indices and corresponding values of nonzero entries for a specified row/column of $𝐌$, in time linear in the number of returned indices;

• $\mathrm{■}$

  Entrywise queries for $𝒃$ and $𝒕$: return $𝒃(k)$ or $𝒕(k)$ in $O(1)$ time for a given index $k\in [n]$.

Our results will assume additional access operations, which will be specified in the statements.

Following the concept of local computation algorithms [36] and the previous work [4], we consider a fixed solution ${𝒙}^{\ast }$ that is determined by $𝐌$ and $𝒃$, and require invoking the algorithm with different $𝒕$ and accuracy parameters returns estimates of ${𝒕}^{\top }{𝒙}^{\ast }$ that are all consistent with the “global” solution ${𝒙}^{\ast }$. Our choice of ${𝒙}^{\ast }$ will be given in Theorem 1.

We also consider the problems of computing (Personalized) PageRank [9] and effective resistance [17] on graphs, viewing them as special cases of our formulation of solving RDD/CDD systems. For these graph problems, we assume the standard adjacency-list model [22], where each degree query takes $O(1)$ time and each neighbor query returns a neighbor index along with the edge weight in $O(1)$ time.

For Personalized PageRank (PPR), we consider a directed graph $G$, a decay factor $\alpha \in (0,1)$, and a source distribution $𝒔\in \{𝒚\in {ℝ}_{\ge 0}^n:{\Vert 𝒚\Vert }_1=1\}$. To ensure that PPR is well-defined, we assume that ${\delta }_G^{+}>0$. The PPR vector ${𝝅}_{G,\alpha ,𝒔}$ is defined as the unique solution to the following two equivalent forms of the PPR equation:

Both equations can be viewed as linear systems of the form $𝐌𝒙=𝒃$, where the coefficient matrices $𝐈-(1-\alpha ){𝐀}_G^{\top }{𝐃}_G^{-1}$ and ${𝐃}_G-(1-\alpha ){𝐀}_G^{\top }$ are both CDDZ and invertible. Note that for the second form, the solution to the corresponding system is ${𝐃}_G^{-1}{𝝅}_{G,\alpha ,𝒔}$, an outdegree-scaled version of the PPR vector. We define the PPR value from $s$ to $t$ as ${𝝅}_{G,\alpha }(s,t):={𝝅}_{G,\alpha ,{𝒆}_s}(t)$, and the PageRank vector ${𝝅}_{G,\alpha }$ as ${𝝅}_{G,\alpha }:={𝝅}_{G,\alpha ,1/n\cdot \mathrm{𝟏}}$. It holds that ${𝝅}_{G,\alpha }(t)=\frac{1}{n}{\sum }_{s\in V}{𝝅}_{G,\alpha }(s,t)$ for all $t\in V$. We also consider the following two equivalent forms of the PageRank contribution equation:

where $t\in V$ is a specified target node and ${𝝅}_{G,\alpha ,t}^{-1}$ is called the PageRank contribution vector to $t$. It holds that ${𝝅}_{G,\alpha ,t}^{-1}(s)={𝝅}_{G,\alpha }(s,t)$ for all $s\in V$. Similarly, both equations can be viewed as linear systems of the form $𝐌𝒙=𝒃$, where the coefficient matrices $𝐈-(1-\alpha ){𝐃}_G^{-1}{𝐀}_G$ and ${𝐃}_G-(1-\alpha ){𝐀}_G$ are both RDDZ and invertible. Note that for the second form, the corresponding vector $𝒃$ is $\alpha {𝐃}_G{𝒆}_t$.

For effective resistance, we consider a connected undirected graph $G$ and two distinct nodes $s,t\in V$. The effective resistance (a.k.a. resistance distance) between $s$ and $t$, denoted by $R_G(s,t)$, is defined as the equivalent resistance between $s,t$ if the graph is thought of as an electrical network with each edge $(u,v)\in E$ having resistance $1/{𝐀}_G(u,v)$. Algebraically, $R_G(s,t)={({𝒆}_s-{𝒆}_t)}^{\top }{𝐋}_G^{+}({𝒆}_s-{𝒆}_t)$. As we will establish in Lemma 27, setting $𝐌={𝐋}_G$ and $𝒃=𝒕={𝒆}_s-{𝒆}_t$ in our formulation yields ${𝒕}^{\top }{𝒙}^{\ast }=R_G(s,t)$. Here, $𝐌={𝐋}_G$ is SDDZ and $𝒃={𝒆}_s-{𝒆}_t\in \mathrm{range}(𝐌)$.

[4] studies the case when the linear system is $𝐒𝒙=𝒃$ for some SDD matrix $𝐒$. Define ${𝐃}_{𝐒}$ as the diagonal matrix that satisfies $𝐃(k,k)=𝐒(k,k)$ for each $k\in [n]$ and $\widetilde{𝐒}:={𝐃}_{𝐒}^{-1/2}{\mathrm{𝐒𝐃}}_{𝐒}^{-1/2}$. They formulate the fixed solution as ${𝒙}^{\ast }:={𝐃}_{𝐒}^{-1/2}{\widetilde{𝐒}}^{+}{𝐃}_{𝐒}^{-1/2}𝒃$ and give a Neumann series expansion of ${𝒙}^{\ast }$. For the algorithmic results, they assume that the algorithm is given $\kappa$, an upper bound on $\kappa (\widetilde{𝐒})$ (alternatively, $\gamma$ as a lower bound on $\gamma (𝐒)$) and set a truncation parameter for the Neumann series as $L:=\mathrm{\Theta }\left(\kappa \log \left(\kappa \cdot \kappa ({𝐃}_{𝐒}){\Vert 𝒃\Vert }_0/\epsilon \right)\right)=\tilde{\mathrm{\Theta }}(\kappa )$.

Based on the truncated Neumann series, they present a randomized algorithm that, given a coordinate $t\in [n]$, computes an estimate ${\widehat{x}}_t$ satisfying $\left|{\widehat{x}}_t-{𝒙}^{\ast }(t)\right|\le \epsilon {\Vert {𝐃}_{𝐒}^{-1}𝒃\Vert }_{\mathrm{\infty }}$ with probability at least $3/4$. The algorithm runs in time $O\left(f(𝐒)L^3\log L/{\epsilon }^2\right)=\tilde{O}\left(f(𝐒){\kappa }^3{\epsilon }^{-2}\right)=\tilde{O}\left(f(𝐒){\gamma }^{-3}{\epsilon }^{-2}\right)$, where $f(𝐒)$ is the maximum time cost to simulate one step in the random walk defined by $𝐒$. This result is implicit in the proof of [3, Theorem 5.1].

On the negative side, [4] proves an $\mathrm{\Omega }\left(\kappa {(𝐒)}^2/{\log }^{3}n\right)=\tilde{\mathrm{\Omega }}\left(\kappa {(𝐒)}^2\right)$ query lower bound (in terms of probing $𝒃$) for achieving a weaker absolute error bound of $\epsilon {\Vert {𝒙}^{\ast }\Vert }_{\mathrm{\infty }}$, for $\kappa (𝐒)=O(\sqrt{n}/\log n)$ and $\epsilon =\mathrm{\Theta }(1/\log n)$. The matrix $𝐒$ in their hard instance is a Laplacian matrix of a fixed unweighted undirected graph with maximum degree $4$ and thus satisfies $\kappa (\widetilde{𝐒})=\mathrm{\Theta }\left(\kappa (𝐒)\right)$. Therefore, this lower bound can also be written as $\tilde{\mathrm{\Omega }}\left(1/\gamma {(𝐒)}^2\right)$. To our knowledge, no other work has explicitly studied sublinear-time SDD/RDD/CDD solvers.

Our first contribution is a Neumann-series-based characterization of a solution ${𝒙}^{\ast }$ to $𝐌𝒙=𝒃$, which is consistent with the solution considered by [4] for SDD systems.

We decompose $𝐌$ uniquely as $𝐌={𝐃}_{𝐌}-{𝐀}_{𝐌}^{\top }$, where ${𝐃}_{𝐌}$ is a diagonal matrix and all diagonal entries of ${𝐀}_{𝐌}^{\top }$ are $0$. Define $\widetilde{𝐌}:={𝐃}_{𝐌}^{-1/2}{\mathrm{𝐌𝐃}}_{𝐌}^{-1/2}$ and ${\widetilde{𝐀}}_{𝐌}^{\top }:={𝐃}_{𝐌}^{-1/2}{𝐀}_{𝐌}^{\top }{𝐃}_{𝐌}^{-1/2}$.

The next theorem gives the definition of ${𝒙}^{\ast }$ and its properties.

Next, we study a truncated version of ${𝒙}^{\ast }$ and upper bound the truncation error. As previous analysis based on eigendecomposition is not directly applicable to asymmetric matrices, we introduce a novel concept called the $p$-norm gap of $𝐌$: for any $p\in [1,\mathrm{\infty }]$, we define the $p$-norm gap of $𝐌$ as

where $q$ is conjugate to $p$. We further define the maximum $p$-norm gap of $𝐌$ as ${\gamma }_{\max }(𝐌):={\max }_{p\in [1,\mathrm{\infty }]}{\gamma }_p(𝐌)$.³³3One could alternatively define ${\gamma }_{\max }(𝐌):={\max }_{p\in \{1,2,\mathrm{\infty }\}}{\gamma }_p(𝐌)$, which yields a possibly smaller quantity but our main results continue to hold with this definition. To our knowledge, no prior work has explicitly studied these quantities.

The following theorem guarantees that for any RDD/CDD matrix $𝐌$, the maximum $p$-norm gap ${\gamma }_{\max }(𝐌)$ lies in $(0,1]$, and when $𝐌$ is SDD, ${\gamma }_{\max }(𝐌)$ coincides with the spectral gap $\gamma (𝐌)$. Thus, it is a natural generalization of the spectral gap to asymmetric matrices.

We will devise algorithms to estimate the truncated version of ${𝒙}^{\ast }$, defined as

for an integer truncation parameter $L$. The next theorem upper bounds the truncation error in terms of a given lower bound on the maximum $p$-norm gap ${\gamma }_{\max }(𝐌)$.

For our algorithmic results, we assume that the algorithm is given a quantity $\gamma >0$ as a lower bound on ${\gamma }_{\max }(𝐌)$ and an accuracy parameter $\epsilon >0$. We use $\gamma$, $\epsilon$, and the specific terms in the accuracy guarantee to set the truncation parameter $L$ according to Equation (6), where we assume that suitable upper bounds on the quantities in the logarithmic factor are known.

The statements of our results will use the following definition. For RDD $𝐌$, we use $f_{\mathrm{row}}(𝐌)$ to denote the maximum time cost to simulate a single step in the random walk defined by the row substochastic matrix $\frac{1}{2}\left(𝐈+{𝐃}_{𝐌}^{-1}\left|{𝐀}_{𝐌}^{\top }\right|\right)$. This quantity depends on the structure and representation of $𝐌$. For instance, if each row of $𝐌$ has at most $d$ nonzero entries, then $f_{\mathrm{row}}(𝐌)=O(d)$; if the nonzero entries in ${𝐀}_{𝐌}$ have equal absolute values and we are allowed to sample a uniformly random index of the nonzero entries in each column of ${𝐀}_{𝐌}$ in $O(1)$ time, then $f_{\mathrm{row}}(𝐌)=O(1)$. For CDD $𝐌$, we define $f_{\mathrm{col}}(𝐌):=f_{\mathrm{row}}\left({𝐌}^{\top }\right)$.

Each of our algorithmic results excels under different system types, access models, and parameter regimes. We present a selection of our results here, with additional results given in the full version of this paper [28]. Furthermore, a remark at the end of this subsection establishes that most of our results have “symmetric” counterparts (e.g., for CDD systems) that can be obtained by exchanging the roles of certain quantities.

We first present our results of using random-walk sampling for RDD systems.

This theorem subsumes the main algorithmic result of [4] for SDD systems and extends it to the RDD case while achieving a mild $(\log L)$-factor improvement. Their original result is recovered as a special case when $𝐌$ is SDD and $𝒕$ is a canonical unit vector. Our algorithm is similar to theirs, but we achieve the improved complexity by adopting a different random-walk sampling scheme.

We also show that for RDDZ $𝐌$ along with nonnegative vectors $𝒃$ and $𝒕$, if we allow the relaxed error bound $\epsilon {\Vert {𝒙}^{\ast }\Vert }_{\mathrm{\infty }}$, the complexity can be improved to depend quadratically on $L$ via a variance analysis of the random-walk sampling process. We adopt this error measure since it provides a natural accuracy guarantee and has been previously studied in [4]. It is shown in [4] that ${\Vert {𝐃}_{𝐌}^{-1}𝒃\Vert }_{\mathrm{\infty }}\le 2{\Vert {𝒙}^{\ast }\Vert }_{\mathrm{\infty }}$ for RDD systems.

The following theorem further considers the relative error guarantee of $\epsilon \cdot {𝒕}^{\top }{𝒙}^{\ast }$. The complexity depends quadratically on $L$ and linearly on ${\Vert 𝒕\Vert }_1{\Vert {𝐃}_{𝐌}^{-1}𝒃\Vert }_{\mathrm{\infty }}/{𝒕}^{\top }{𝒙}^{\ast }$ and is also achieved using random-walk sampling.

Our next result leverages the local push method to derive a deterministic algorithm for special RCDD systems, whose complexity depends linearly on $1/\epsilon$.

Lastly, applying the bidirectional method to special RCDD systems yields complexity bounds with improved dependence on $L$ and $\epsilon$, achieving either $L^{7/3}{\epsilon }^{-2/3}$ or $L^{5/2}{\epsilon }^{-1}$ using different parameter settings.

Our results can be directly applied to the PPR or PageRank contribution equations and the single-pair effective resistance problem to yield complexity bounds for different graph types under various access models and accuracy guarantees. Notably, as we shall see, for PageRank computation, half the decay factor serves as a lower bound on the maximum $p$-norm gap of the corresponding systems; for effective resistance computation, the spectral gap of the graph Laplacian equals the maximum $p$-norm gap of the system. Rather than exhaustively presenting all applications of our results, here we highlight selected results that generalize and improve upon the previously best complexity bounds.

The following theorem is derived by applying Theorem 6 to the PageRank equation (2) combined with tighter lower bounds on ${𝝅}_{G,\alpha }(t)$ that we establish for Eulerian graphs.

This improves over the previously best upper bounds of $O\left(\frac{1}{{\epsilon }^2{\delta }_G}\cdot \min (\frac{d_G(t)}{{\alpha }^2},\frac{\sqrt{m}}{{\alpha }^2})\right)$, which was given by [45] and stated for unweighted undirected graphs. Our algorithm is essentially the same as that in [45], which generates random walks from the target node $t$, and the improvement comes from our tighter lower bounds on ${𝝅}_{G,\alpha }(t)$ for Eulerian graphs.

For effective resistance computation, recall that we consider connected undirected graphs $G$. Lemma 27 establishes that the value ${𝒕}^{\top }{𝒙}^{\ast }$ that our algorithms approximate equals $R_G(s,t)$ when setting $𝐌={𝐋}_G$ and $𝒃=𝒕={𝒆}_s-{𝒆}_t$. Thus, Theorems 4 and 8 directly imply the following complexity bounds for estimating effective resistance.

This subsumes the previous bounds of $O\left(\min (\frac{L^3\log L}{{\epsilon }^2\cdot \min {(d_G(s),d_G(t))}^2},\frac{L^{7/3}\log L}{{\epsilon }^{2/3}})\right)$ given by [14] with essentially the same setting of $L$. Moreover, since $R_G(s,t)$ can be lower bounded by $1/2/\min (d_G(s),d_G(t))$ [34, Corollary 3.3], by setting the absolute error parameter $\epsilon$ to be ${\epsilon }_r\cdot 1/2/\min (d_G(s),d_G(t))$, the last bound in our result implies that an estimate of $R_G(s,t)$ within relative error ${\epsilon }_r$ can be computed in time $O\left(\min {(d_G(s),d_G(t))}^{1/2}{L}^{5/2}{\epsilon }_r^{-1}\right)$. This improves over the previous bound of $O(\min {(d_G(s),d_G(t))}^{1/2}{L}^3\log L\cdot {\epsilon }_r^{-1})$ given by [52]. Our algorithms are essentially the same as those in [14, 52], i.e., using random-walk sampling and the bidirectional method, and the improvements stem from our simpler sampling scheme that avoids using extra data structures and a refined analysis.

On the other hand, known hardness results for local PageRank and effective resistance computation can potentially yield lower bounds for sublinear-time solvers. We highlight the following result, derived by establishing a reduction from estimating single-node PageRank on undirected graphs to solving SDD systems and applying the lower bound for PageRank computation from [45].

This result gives an $\mathrm{\Omega }(1/\epsilon )$ lower bound for local SDD solvers with accuracy guarantee $\epsilon {\Vert {𝒙}^{\ast }\Vert }_{\mathrm{\infty }}$, demonstrating the necessity of a linear dependence on $1/\epsilon$ in our Theorem 4. This lower bound on the accuracy parameter complements the lower bound on the spectral gap given by [4].

Inspired by the local push algorithms ForwardPush and BackwardPush, we formulate our Push algorithm as a unified primitive for estimating a summation of matrix powers applied to a vector, which is closely related to our approach to solving RDD/CDD systems. This abstraction reveals that ForwardPush and BackwardPush, despite appearing as distinct algorithms for two problems with different propagation strategies, are equivalent to applying Push to different linear systems, modulo variable scaling. The apparent differences arise from the scaling and the distinct behaviors that Push exhibits on different types of linear systems.

Specifically, for PageRank computation, ForwardPush from node $s$ corresponds to applying Push to Equation (2) with $𝒔={𝒆}_s$ and outdegree scaling, while BackwardPush from node $t$ applies Push to the contribution equations (3) or (4). Our analysis demonstrates that Push provides closed-form complexity bounds on certain CDD systems and accuracy guarantees on RDD systems, which explains the known computational properties of ForwardPush and BackwardPush on directed graphs.

For RCDD systems, Push inherits both complexity and accuracy advantages, clarifying why ForwardPush and BackwardPush perform well on undirected graphs. This perspective further reveals that on Eulerian graphs, ForwardPush from any node is equivalent to BackwardPush from that node on the transpose graph, establishing a basic connection previously unrecognized.

Our characterizations of the solution ${𝒙}^{\ast }$ and the $p$-norm gaps are based on fundamental properties of Neumann series and restricted linear maps. Since asymmetric matrices may be non-diagonalizable, the eigendecomposition techniques used for symmetric matrices become inapplicable. We address this challenge by analyzing the operator norms of restricted linear maps to establish series convergence and derive truncation error bounds.

Our algorithms estimate ${𝒕}^{\top }{𝒙}_L^{\ast }=\frac{1}{2}{𝒕}^{\top }{\sum }_{\mathrm{\ell }=0}^{L-1}{\left(\frac{1}{2}\left(𝐈+{𝐃}_{𝐌}^{-1}{𝐀}_{𝐌}^{\top }\right)\right)}^{\mathrm{\ell }}{𝐃}_{𝐌}^{-1}𝒃$ by adapting three techniques for random-walk probability estimation on graphs: random-walk sampling [39, 19], local push [2, 1], and the bidirectional method [32].

When $𝐌$ is RDD, the matrix $\frac{1}{2}\left(𝐈+{𝐃}_{𝐌}^{-1}\left|{𝐀}_{𝐌}^{\top }\right|\right)$ is row substochastic, enabling us to interpret ${𝒕}^{\top }{𝒙}_L^{\ast }$ as an expectation over random walks that start from the distribution $|𝒕|/{\Vert 𝒕\Vert }_1$ and transition according to $\frac{1}{2}\left(𝐈+{𝐃}_{𝐌}^{-1}\left|{𝐀}_{𝐌}^{\top }\right|\right)$. Such random walks may terminate early, and the algorithm needs to record the signs of the entries in ${𝐀}_{𝐌}$ along the walk. This method extends the approach in [4], and we adopt a different sampling scheme and conduct variance analysis in some special cases to reduce the dependence on $L$ in the complexity.

Based on local push methods, we formulate a Push primitive that estimates the vector ${𝒙}_L^{\ast }={\sum }_{\mathrm{\ell }=0}^{L-1}{\left(\frac{1}{2}\left(𝐈+{𝐃}_{𝐌}^{-1}{𝐀}_{𝐌}^{\top }\right)\right)}^{\mathrm{\ell }}{𝐃}_{𝐌}^{-1}𝒃$ through deterministic local computation. Push maintains coordinate variables initialized to ${𝐃}_{𝐌}^{-1}𝒃$ and iteratively applies the linear operator $\frac{1}{2}\left(𝐈+{𝐃}_{𝐌}^{-1}{𝐀}_{𝐌}^{\top }\right)$ to selected coordinates. Our analysis relies on invariant properties preserved by these operations, including an inequality variant that helps to handle negative entries. This approach yields closed-form accuracy guarantees for RDD systems and complexity bounds for special CDD systems. Combining these two aspects yield our result for special RCDD systems.

The bidirectional method combines random-walk sampling and local push from two directions. We adapt the BiPPR framework [32], performing Push from ${𝐃}_{𝐌}^{-1}𝒃$ and exploiting the invariant property to construct an unbiased estimator for ${𝒕}^{\top }{𝒙}_L^{\ast }$, which can be sampled via random walks from $|𝒕|/{\Vert 𝒕\Vert }_1$ when $𝐌$ is RDD. By balancing the costs of both components, we achieve improved dependence on $L$ and $\epsilon$, particularly for RCDD systems.

Crucially, we can transpose the expression for ${𝒕}^{\top }{𝒙}_L^{\ast }$ to obtain an equivalent summation with ${𝐀}_{𝐌}^{\top }$ replaced by ${𝐀}_{𝐌}$ and the roles of $𝒃$ and $𝒕$ exchanged. This allows us to alternatively apply Push from ${𝐃}_{𝐌}^{-1}𝒕$ and random-walk sampling from $|𝒃|/{\Vert 𝒃\Vert }_1$ when $𝐌$ is CDD, leading to symmetric algorithmic procedures and complexity results.

The remainder of this paper is organized as follows. Section 2 introduces more related work. Section 3 elaborates on our formulation of ${𝒙}^{\ast }$ and $p$-norm gaps and proves Theorems 1 to 3. In Sections 4, 5, and 6, we present our algorithms based on random-walk sampling, local push, and the bidirectional method, respectively, and prove Theorems 4 to 8. After that, Sections 7 and 8 relate our problem and results to local computation of PageRank and effective resistance, respectively. In the full version of this paper [28], we provide more preliminaries, discussions on future directions, detailed explanations of the relationship between our Push primitive and the ForwardPush and BackwardPush algorithms, additional results, and all omitted proofs.