Iterative Methods via Locally Evolving Set Process

Authors: baojian  zhou
yifan  sun
reza  babanezhad harikandeh
xingzhi  guo
deqing  yang
yanghua  xiao
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-10 10:04:40

Abstract: Given the damping factor $\alpha$ and precision tolerance $\epsilon$, \citet{andersen2006local} introduced Approximate Personalized PageRank (APPR), the \textit{de facto local method} for approximating the PPR vector, with runtime bounded by $\Theta(1/(\alpha\epsilon))$ independent of the graph size. Recently,  Fountoulakis \& Yang asked whether faster local algorithms could be developed using $\tilde{\mathcal{O}}(1/(\sqrt{\alpha}\epsilon))$ operations. By noticing that APPR is a local variant of Gauss-Seidel, this paper explores the question of *whether standard iterative solvers can be effectively localized*. We propose to use the *locally evolving set process*, a novel framework to characterize the algorithm locality, and demonstrate that many standard solvers can be effectively localized. Let $\overline{\operatorname{vol}}{ (\mathcal S_t)}$ and $\overline{\gamma_t}$ be the running average of volume and the residual ratio of active nodes $\textstyle \mathcal{S_t}$ during the process. We show $\overline{\operatorname{vol}}{ (\mathcal S_t)}/\overline{\gamma_t} \leq 1/\epsilon$ and prove APPR admits a new runtime bound $\tilde{\mathcal{O}}(\overline{\operatorname{vol}}(\mathcal S_t)/(\alpha\overline{\gamma_t}))$ mirroring the actual performance. Furthermore, when the geometric mean of residual reduction is $\Theta(\sqrt{\alpha})$, then there exists $c \in (0,2)$ such that the local Chebyshev method has runtime $\tilde{\mathcal{O}}(\overline{\operatorname{vol}}(\mathcal{S_t})/(\sqrt{\alpha}(2-c)))$ without the monotonicity assumption. Numerical results confirm the efficiency of this novel framework and show up to a hundredfold speedup over corresponding standard solvers on real-world graphs.

Main paper: 
 introduction: 
 Personalized PageRank (PPR) vectors are key tools for graph problems such as clustering [2, 3, 30, 36, 54, 57] , diffusion [10, 14, 15, 29] , random walks [25, 32, 44] , neural net training [7, 27, 20, 21] , and many others [17, 48] . The Approximate PPR (APPR) [2] and its many variants [6, 9, 13, 37] efficiently approximate PPR vectors by exploring the neighbors of a specific node at each time, only requiring access to a tiny part of the graph -hence the number of operations needed is independent of graph size. These local solvers are well-suited for large-scale graphs in modern graph data analysis. Specifically, let A and D be the adjacency and degree matrices of a graph G, respectively. Given a source node s and the damping factor α ∈ (0, 1), this paper studies local solvers for the linear systemEQUATIONwhere e s is the standard basis of s and π is the PPR vector [2, 12, 37] . Given the error tolerance ϵ, a local solver needs to find π such that ∥D -1 ( ππ)∥ ∞ ≤ ϵ without accessing the entire graph G. 2 Andersen et al. [2] proposed the local APPR algorithm, which pushes large residuals to neighboring nodes until all residuals are small. Its runtime is upper bounded by Θ(1/(αϵ)) independent of graph size. Based on a variational characterization of Equ. (1), Fountoulakis et al. [13] reformulated the problem as optimizing a quadratic objective plus ℓ 1 -regularization and later asked [12] whether there exists a local solver with runtime Õ(1/( √ αϵ)). This corresponds to an accelerated rate since α is the strongly convex parameter. Recently, Martínez-Rubio et al. [37] provided a method based on a nested subspace pursuit strategy, and the corresponding iteration complexity is bounded byÕ(|S * |/ √ α)where S * is the support of the optimal solution. This bound deteriorates to Õ(n/ √ α) when the solution is dense, with n representing the number of nodes in G, which could be less favorable than that of standard solvers under similar conditions. Moreover, the nested computational structure provides a constant factor overhead, which could be significant in practice.The bound analysis of the above local methods critically depends on the monotonicity properties of the designed algorithms. These requirements may hinder the development of simpler and faster local linear solvers that lack such monotonicity properties. Specifically, the runtime analysis of APPR relies on the non-negativity and decreasing monotonicity of residuals. Conversely, the runtime bounds developed in Fountoulakis et al. [13] and Martínez-Rubio et al. [37] depend on the monotonicity of variable updates, ensuring that the sparsity of intermediate variables increases monotonically.Our contributions. Based on a refined analysis of APPR, our starting point is to demonstrate that APPR is a local variant of Gauss-Seidel Successive Overrelaxation (GS-SOR) that can be treated as an evolving set process. 3 This insight leads us to explore whether standard solvers can be effectively localized to solve Equ. (1) . To develop faster local methods with improved local bounds, we propose a novel locally evolving set process framework inspired by its stochastic counterpart [38] . This framework enables the development of faster local methods and circumvents the monotonicity requirement barrier in runtime complexity analysis in the existing literature. For example, our analysis of the local Chebyshev method does not depend on the monotonicity of residual or the active node sets processed. Specifically,• As a core tool, we propose an algorithm framework based on the locally evolving set process. We show that APPR is a local variant of GS-SOR using this process. This framework is powerful enough to facilitate the development of new local solvers. Specifically, standard gradient descent (GD) can be effectively localized for solving this problem and admits Θ(1/(αϵ)) runtime bound. • This local evolving set process provides a novel way to characterize the algorithm locality; hence, new runtime bounds can be derived. Let vol(S t ) and γ t be the running average of volume and the residual ratio of active nodes S t during the process; we prove the ratio vol(S t )/γ t serving as a lower bound of 1/ϵ. We further show both APPR and local GD have Õ(vol(S t )/(αγ t )) runtime bound mirroring the actual performance of these two methods. • Using our framework, we show there exists c ∈ (0, 2) such that both the localized Chebyshev and Heavy-Ball methods admit runtime bound Õ(vol(S t )/( √ α(2c))) with the assumption that the geometric mean of active ratio factors is Θ( √ α). Importantly, our analysis does not require any monotonicity property. The technical novelty is that we effectively characterize residuals of these two methods by using second-order difference equations with parameterized coefficients. • We demonstrate, over 17 large graphs, that these localized methods can significantly accelerate their standard counterparts by a large margin. Furthermore, our proposed LOCSOR, LOCCH, and LOCHB are significantly faster than APPR and ℓ 1 -based solvers on two huge-scale graphs.Paper structure. We begin by clarifying notations and reviewing APPR in Sec. 2. Sec. 3 introduces the locally evolving set process. Sec. 4 presents localized Chebyshev and Heavy-Ball methods along with our novel techniques. We discuss open questions in Sec. 5 . Experiments and conclusions are covered in Sec. 6 and 7, respectively. Detailed related works and all missing proofs are included in the Appendix. Our code is available at https: // github. com/ baojian/ LocalCH .In the rest of the appendix, we use the following notations:Description GAn undirected connected simple graph with unit weights. 
notations and preliminaries: 
 Notations. We consider an undirected simple graph G(V, E) where V = {1, 2, . . . , n} and E ⊆ V ×V with |E| = m are the node and edge sets, respectively. The set of neighbors of v is denoted as N (v) ⊆ V. The adjacency matrix A of G assigns unit weight a u,v = 1 if (u, v) ∈ E and 0 otherwise. The v-th entry of the degree matrix D is d v = |N (v)|. Given S ⊆ V, we define the volume of S as vol(S) ≜ v∈S d v . The support of x ∈ R n is the set of nonzero indices supp(x) ≜ {v : x v ̸ = 0, v ∈ V}. The eigendecomposition of D -1/2 AD -1/2 = V ΛV ⊤ where each column of V is an eigenvector and Λ = diag(λ 1 , λ 2 , . . . , λ n ) with 1 = λ 1 ≥ λ 2 ≥ • • • ≥ λ n ≥ -1.We use (p, z) for solving Equ. (1) while use (x, r) for solving Equ. (3) or Equ. (4) . With the initial setting p ← 0, z ← e s , APPR obtains a local estimate of π denoted as p by using a sequence of PUSH operations defined as APPR(α, ϵ, s, G) :EQUATIONAlgo. 1 PUSH(u, α, p, z)1: ν = z u 2: p u ← p u + α • ν 3: z u ← (1 -α) • ν/2 4: for v ∈ N (u) do 5: z v ← z v + (1-α)ν 2du 6: Return (p, z)At each repeat step (p, z) ← PUSH(u, α, p, z), it synchronously updates both p and residual z whenever there exists an active node u ∈ V (a node with a large residual, i.e., z u ≥ ϵd u ). Specifically, for each active u, it updates p u , z u , and z v for v ∈ N (u) by using a PUSH operator illustrated on the left. It stops when no active nodes are left. APPR can be implemented locally so that the runtime is independent of G. In particular, vol(supp(p)) is locally bounded, demonstrating the sparsity effect. We restate the existing main results as follows. Lemma 2.1 (Runtime bound of APPR [2] ). Given α ∈ (0, 1) and the precision ϵ ≤ 1/d s for node s ∈ V with p ← 0, z ← e s at the initial, APPR(α, ϵ, s, G) defined in (2) returns an estimate p of π. There exists a real implementation of (2) (e.g., Algo. 2) such that the runtime T APPR satisfies T APPR ≤ Θ (1/(αϵ)) .Furthermore, the estimate π := p satisfies ∥D -1 ( ππ)∥ ∞ ≤ ϵ and vol(supp( π)) ≤ 2/((1α)ϵ).The main argument for proving Lemma 2.1 is critically based on: 1) z ≥ 0 and ∥z∥ 1 decreases during the updates; 2) for each active u, z u ≥ ϵd u , implying that ∥z∥ 1 is decreased by at least αϵd u , consecutively leading to u d u ≤ 1/(αϵ).To approximate the PPR vector π, the original linear system in Equ. (1) can be reformulated as an equivalent symmetric version defined asEQUATIONwhere again e s is the standard basis of s, and Q is a symmetric positive-definite M -matrix with all eigenvalues in [ 2α 1+α , 2 1+α ]. To solve Equ. ( 3) is equivalent to solving a quadratic problemx * = arg minEQUATIONwhere f is strongly convex with condition number 1/α. Indeed, Equ. ( 3) is a symmetrized version of Equ. (1) and has a unique solution x * = Q -1 b. The PPR vector π can be recovered from x * by π = D 1/2 x * . It is convenient to denote estimate of π as π (t) ≜ D 1/2 x (t) . Given x (t) , we define the residual r (t) ≜ b -Qx (t) . If x (t) is returned by a local solver for solving either Equ. (3) or Equ. (4), we then equivalently require ∥D -1/2 (x (t)x * )∥ ∞ ≤ ϵ. Hence, it is enough to have a stop condition ∥D -1/2 r (t) ∥ ∞ ≤ 2αϵ/(1 + α) for local solvers of Equ. (3) and Equ. (4). 4Fountoulakis et al. [13] demonstrated that APPR is equivalent to a coordinate descent solver for minimizing f in Equ. (4) and introduced an ISTA-style solver by minimizing f (x) + ϵα∥D 1/2 x∥ 1 , which provides a method with runtime bound Õ(1/(ϵα)) for achieving the same estimation guarantee of APPR. On one hand, one may note that the runtime bound Θ(1/(αϵ)) provided in Lemma 2.1 becomes less valuable when ϵ ≤ 1/m; on the other hand, all previous local variants [6, 9, 13, 37] of APPR are critically based on some monotonicity property. This limitation could impede the development of faster local methods that might violate the monotonicity assumption. The following two sections present the techniques and tools to address these challenges. 
local methods via evolving set process: 
 Our investigation begins with the locally evolving set process, as inspired by the stochastic counterpart [38] . The process reveals that APPR is essentially a local variant of GS-SOR. We then show how to use this process to build faster local solvers based on GS-SOR. We further develop a local parallelizable gradient descent with runtime Θ(1/(αϵ)).Given α, ϵ, s, and G, a local solver for Equ. (3) keeps track of an active set S t ⊂ V at each iteration t.That is, only nodes in S t are used to update x or r. The next set S t+1 is determined by current S t and an associated local solver A. We define this process as the following local evolving set system.Definition 3.1 (Locally evolving set process). Given a parameter configuration θ ≜ (α, ϵ, s, G), and a local iterative method A, the locally evolving set process generates a sequence of (S t , x (t) , r (t) ) representing as the following dynamic system S t+1 , x (t+1) , r (t+1) = Φ θ S t , x (t) , r (t) , A , ∀t ≥ 0,where S t+1 ⊆ S t ∪ (∪ u∈St N (u)) and we denote the active set S t = {u 1 , u 2 , . . . , u |St| }. The set S t is maintained via a queue data structure. We say this process converges when the last set S T = ∅ if there exists such T ; the generated sequence of active nodes are(S 0 , x (0) , r (0) ) → (S 1 , x (1) , r (1) ) → (S 2 , x (2) , r (2) ) → . . . → (S T = ∅, x (T ) , r (T ) ).The runtime of the local solver, A for this whole local process, is then defined as EQUATIONwhere ∆ i is a smaller time magnitude. We define ∆ i = (i -1)/|S t | for the analysis of APPR and LOCSOR while ∆ i = 0 for LOCGD in our later analysis. In the rest, we denote I T = supp(r (T ) ). To demonstrate this local process, Fig. 1 shows vol(S t ) of APPR peaks at the early stage, and the active ratio decreases as the active volume diminishes. The quantity vol(S T )/γ T is strictly smaller than 1/ϵ, indicating that it could serve as a better factor in the runtime analysis.×10 6 1/ vol(St) γ t • ln 2 (1 -α)|It|We first demonstrate how this locally evolving set process can represent APPR. For solving Equ. (1), the set S 0 = {s} and the queue-based of APPR (see Algo. 2 in Appendix A) naturally forms a sequence of active sets from S 0 = {s} to S T = ∅, hence converging. Active nodes u in queue satisfy z u ≥ ϵd u . To delineate successive iterations S t and S t+1 , one can insert * at the beginning of S t . After processing S t , it serves as an indicator for the next iteration. The star * is reinserted into the queue iteratively until the queue is empty. We use a slightly different notation for presenting tuple (S t , p (t) , z (t) ) to consistent with Sec. 2.1 and write out such evolving process as followsEQUATIONThe following lemma establishes the equivalence between APPR and the local variant of GS-SOR method (see Appendix B.1) and provides a new evolving-based bound. Lemma 3.2 (New local evolving-based bound for APPR). Let M = α -1 I -1-α 2 I + AD -1 and s = e s . The linear system M π = s is equivalent to Equ. (1) . Given p (0) = 0, z (0) = e s with ω ∈ (0, 2), the local variant of GS-SOR (15) for M π = s can be formulated asp (t+∆i+1) ← p (t+∆i) + ωz (t+∆i) ui M uiui e ui , z (t+∆i+1) ← z (t+∆i) - ωz (t+∆i) ui M uiui M e ui ,where u i is an active node in S t satisfying z ui ≥ ϵd ui and ∆ i = (i -1)/|S t |. Furthermore, when ω = 1+α 2 , this method reduces to APPR given in (7) , and there exists a real implementation (Aglo. 2) of APPR such that the runtime T APPR is bounded, that isTAPPR ≤ vol(ST ) αγT ln CT ϵ , where vol(ST ) γT ≤ 1 ϵ , CT = 2 (1 -α)|IT | , γT ≜ 1 T T -1 t=0 |S t | i=1 |z (t+∆ i ) u i | ∥z (t) ∥1.Lemma 3.2 points to the sub-optimality of APPR, as GS-SOR allows for a larger ω. For solving Equ. (3), since APPR essentially serves as a local variant of GS-SOR, we can develop a faster local variant based SOR. To extend this method to solve Equ. (3), we propose a local GS-SOR based on an evolving set process, namely LOCSOR, as the followingEQUATIONWhen ω ∈ (0, 1], the residual r is still nonnegative and monotonically decreasing, we establish the convergence of LOCSOR stated in the following theorem. Theorem 3.3 (Runtime bound of LOCSOR (ω = 1)). Given the configuration θ = (α, ϵ, s, G) with α ∈ (0, 1) and ϵ ≤ 1/d s and let r (T ) and x (T ) be returned by LOCSOR defined in (8) for solving Equ. (3). There exists a real implementation of (8) such that the runtime T LOCSOR is bounded by1 + α 2 • vol(S T ) αγ T 1 - ∥D 1/2 r (T ) ∥ 1 ∥D 1/2 r (0) ∥ 1 ≤ T LOCSOR ≤ 1 + α 2 • min 1 αϵ , vol(S T ) αγ T ln C ϵwhere vol(S T ) and γ T are defined in (6) and C = 1+α (1-α)|I T | with I T = supp(r (T ) ). Furthermore, vol(S T )/γ T ≤ 1/ϵ and the local estimate Our new evolving bound Õ(vol(S T )/(αγ T )) mirroring the actual performance of APPR and empirically much smaller than Θ(1/(αϵ)) as illustrated in Fig. 2 . Our lower bounds are quite effective when ϵ is relatively large, while our upper bound is better than Anderson's when ϵ is small. When ϵ ≪ Θ(1/m), this new bound is superior to both O(1/(αϵ)) and Õ(1/( √ αϵ)). This superiority is evident when compared to algorithms like ISTA or FISTA [5] to minimize the ℓ 1 -regularization of f for obtaining an approximate solution of Equ. (3) . Additionally, when ω ∈ (1, 2) and recalling that Q is an M -matrix, the standard analysis of SOR shows that the spectral norm of the iteration matrix must be larger than |ω -1|. Hence, 0 < ω < 2 if and only if global SOR converges [55] . When ω * is optimal (the point that the spectral radius of the iteration matrix is minimized), we have the following result.π := D 1/2 x (T ) satisfies ∥D -1 ( π -π)∥ ∞ ≤ ϵ.Corollary 3.4. Let ω = ω * ≜ 2/(1 + 1 -(1 -α) 2 /(1 + α) 2) and S t = V, ∀t ≥ 0 during the updates, the global version of LOCSOR has the following convergence bound∥r (t) ∥ 2 ≤ 2 (1 + α) √ d s 1 - √ α 1 + √ α + ϵ t t ,where ϵ t are small positive numbers with lim t→∞ ϵ t = 0.Asymptotically, when ϵ t = o( √ α), then the runtime of global LOCSOR is Õ(m/ √ α) where Õ hides log 1/ϵ. The main difficulty of analyzing the optimal local LOCSOR is that the nonnegativity and monotonicity of r (t) do not hold. Instead, by using a parameterized second-order difference equation, we develop new techniques based on the Chebyshev method detailed in Sec. 4.One disadvantage of LOCSOR is its limited potential for parallelization. The standard GD x (t+1) = x (t) -∇f (x (t) ) (step size = 1), in contrast, is easy to parallelize across the coordinates of the update. Instead of updating r and x synchronously per-coordinate, we propose the followingΦ θ S t , x (t) , r (t) , A = LOCGD : x (t+1) ← x (t) + r (t) St , r (t+1) ← r (t) -Qr (t) St (9)Every coordinate in S t is updated in parallel at iteration t. Interestingly, LOCGD exhibits nonnegativity and monotonicity properties, and its runtime complexity is similar to that of LOCSOR, as stated in the following theorem (To remind, ∆ i = 0 for γ T of LOCGD in Equ. ( 6) ). Theorem 3.5 (Runtime bound of LOCGD). Given the configuration θ = (α, ϵ, s, G) with α ∈ (0, 1) and ϵ ≤ 1/d s and let r (T ) and x (T ) be returned by LOCGD defined in (9) for solving Equ. (4). There exists a real implementation of (9) such that the runtime T LOCGD is bounded by1 + α 2 • vol(S T ) αγ T 1 - ∥D 1/2 r (T ) ∥ 1 ∥D 1/2 r (0) ∥ 1 ≤ T LOCGD ≤ 1 + α 2 • min 1 αϵ , vol(S T ) αγ T ln C ϵ , where C = (1 + α)/((1 -α)|I T |), I T = supp(r (T ) ). Furthermore, vol(S T )/γ T ≤ 1/ϵ and the estimate π := D 1/2 x (T ) satisfies ∥D -1 ( π -π)∥ ∞ ≤ ϵ.Note that γ T of LOCGD is empirically smaller than that of LOCSOR. Hence, LOCGD is empirically slower than LOCSOR by only a small constant factor (e.g., twice as slow), a finding consistent with observations of their standard counterparts [19] . Nonetheless, LOCGD is much simpler and more amenable to parallelization on platforms such as GPUs compared to APPR. 
accelerated local iterative methods: 
 This section presents our key contributions where we propose faster local methods based on the Chebyshev method for solving Equ. (3) and the Heavy-Ball (HB) method for Equ. (4).Compared with GS and GD, the standard Chebyshev method offers optimal acceleration in solving Equ. (3) . Following existing techniques (e.g., see d'Aspremont et al. [11] ), we show there exists an upper runtime bound Õ(m/ √ α) to meeting the stopping condition where Õ hides log 1/ϵ (we presented it in Theorem C.6). Hence, the Chebyshev method is one of the optimal first-order linear solvers for solving Equ. (3). However, localizing Chebyshev poses greater challenges due to the additional momentum vector involved in updating x (t) . Our key observation is that if a substantial reduction in the magnitudes of r (t) is required within a subset of S t , then the corresponding momentum coordinates are likely to possess significant acceleration energy. Intuitively, a viable strategy involves localizing both the residual and momentum vectors. For t ≥ 1, denote the "momentum" vector as ∆ (t) := x (t)x (t-1) and δ t:t+1 = δ t δ t+1 , we propose the localized Chebyshev as the followingEQUATIONwhere t ≥ 1 with the initials x (0) = 0, x (1) = r (0) , δ 0 = 0, δ 1 = (1α)/(1 + α), and W = D -1/2 AD -1/2 is normalized adjacency matrix. Our key strategy for analyzing (10) is to rewrite the updates of r (t) as a nonhomogeneous second-order difference equation (see details in Lemma C.8)EQUATIONwhere we denote S j,t = S j ∩ • • • ∩ S t-1 ∩ S t given t ≥ j ≥ 0 where S t = V\S t . In the rest, we define α = (1 -√ α)/(1 + √ α) and recall the eigendecomposition ofD -1/2 AD -1/2 = V ΛV ⊤ .Based on the above Equ. (11) , we have the following key lemma.Lemma 4.1. Given t ≥ 1, x (0) = 0, x (1) = r (0) . The residual r (t) of LOCCH defined in (10) can be expressed as the followingV ⊤ r (t) = δ 1:t Z t V ⊤ r (0) + δ 1:t tu 0,t + 2 t-1 k=1 δ k+1:t (t -k)u k,t ,where Z t is a diagonal matrix such that ∥Z t ∥ 2 ≤ 1 andu k,t =      t-1 j=1 δ2:j t H j,t I -1-α 1+α Λ V ⊤ r (0) S0,j if k = 0 t-1 j=k δ k+1:j (t-k) H j,t 1+α 1-α I -Λ V ⊤ r (k) S k,j if k ≥ 1, where H k,t is a diagonal matrix such that ∥H k,t ∥ 2 ≤ t -k.This key lemma essentially captures the process of residual reduction r (t) of LOCCH. Specifically, given current iteration t, we define the running residual reduction rate for r (k) with k = 0, 1, 2, . . . , t -1 of step t as β k,t , that is,EQUATIONNote thatβ k,t ≤ 2∥r (k) S k ∥ 2 (1 -α)∥r (k) ∥ 2 ≈O(ϵ) + t-1 j=k+1 4α j-k (t -j) (1 -α)(t -k) ∥r (k) S k,j ∥ 2 ∥r (k) ∥ 2 ≤ 4 α (1-α)(1-α), where whether the last term can be even smaller depends on ∥r(k) S k,j ∥ 2 for S k,j = S k ∩ • • • ∩ S j-1 ∩ S j .However, we notice that the running geometric mean β t ≜ ( t-1 j=0 (1 + β j )) 1/t is even smaller in practice. Based on these observations and the assumption on β t , we establish the following theorem. Theorem 4.2 (Runtime bound of LOCCH). Given the configuration θ = (α, ϵ, s, G) with α ∈ (0, 1) and ϵ ≤ 1/d s and let r (T ) and x (T ) be returned by LOCCH defined in (10) for solving Equ. (3). For t ≥ 1, the residual magnitude ∥r (t) ∥ 2 has the following convergence bound∥r (t) ∥ 2 ≤ δ 1:t t-1 j=0 (1 + β j )y t ,where y t is a sequence of positive numbers solving y t+1 -2y t + y t-1 /((1 + β t-1 )(1 + β t )) = 0 with y 0 = y 1 = ∥r (0) ∥ 2 . Suppose the geometric mean β t ≜ ( t-1 j=0 (1 + β j )) 1/t of β t be such thatβ t = 1 + c √ α 1-√ α where c ∈ [0, 2). There exists a real implementation of (10) such that the runtime T LOCCH is bounded by Golub & Overton [18] considered the approximate Chebyshev method by assuming that the inexact residual is sufficiently smaller than ϵ∥r (t) ∥ 2 , where ϵ must be small enough to ensure convergence. However, this assumption is overly stringent for our case. The novelty of our analysis lies in a more elegant treatment of a parameterized second-order difference equation, allowing us to circumvent this assumption. The nested APGD(ε), namely ASPR proposed in Martínez-Rubio et al. [37] has runtime complexityT LOCCH ≤ Θ (1 + √ α)vol(S T ) √ α(2 -c) ln 2y T ϵ .Õ(|S * | vol (S * ) / √ α + |S * | vol (S * ))where S * is the optimal support of arg min x f (x) + ϵα∥D 1/2 x∥ 1 and vol(S * ) = nnz (Q S * ,S * ). Although it is difficult to compare our bound to this, one limitation of ASPR is that it assumes to call APGD(ε) O(|S * |) times to finish in the worst case. However, our iteration complexity isÕ(1/( √ α(2 -c))).Asymptotically, c = o( √ α) (ϵ → 0), our complexity is Õ(1/ √ α) could be better than Õ(|S * |/ √ α).Fig. 3 presents a preliminary study on ASPR, indicating that it requires more operations than APPR.We conclude our analysis by presenting a similar result for the local Heavy Ball (HB). Note the HB method is the one when δ t δ t+1 → α2 where α = (1 -√ α)/(1 + √ α). Hence, it has similar convergence analyses as to LOCCH shown in Theorem D.8. The LOCHB has the following updatesEQUATION5 Generalization and Open Problems T LOCA = Θ vol(S T ) √ αγ T ln C ϵ, where C a graph-independent constant.Additionally, this work primarily focuses on using first-order neighbors at each iteration. An area for future exploration is generalizing to higher-order neighbors to determine if this leads to faster or more efficient methodologies, which remains an open question. 
experiments: 
 We conduct experiments over 17 graphs to solve (3) and explore the local clustering task. We address the following questions: 1) Can iterative solvers be effectively localized? 2) How does the performance of accelerated local methods compare to non-accelerated ones? 3) Can our proposed methods reduce the number of operations required for local clustering? 6Baselines. We consider four baselines: 1) Conjugate Gradient Method (CGM) as a benchmark to compare local and non-local methods; 2) ISTA, the local method proposed by Fountoulakis et al. [13] ; 3) FISTA, the momentum-based local algorithm proposed by Hu [22] ; and 4) APPR, the classic local method proposed by Andersen et al. [2] . All methods are implemented in Python 3.10 with the numba library [33] .Efficiency of localized algorithms. To compare local solvers to their standard counterparts, we set α = 0.1, randomly select 50 nodes from each graph to serve as e s in (3), and run standard GD, SOR, HB, and CH solvers along with their local counterparts: LOCGD, LOCSOR, LOCHB, and LOCCH. We measure the efficiency by the speedup, defined as the ratio between the runtime of the standard and local solver. The range of ϵ is ϵ ∈ [ α 2(1+α)ds , 10 -4 /n]. The results, presented in Fig. 4 , clearly indicate that our design demonstrates significant speedup, especially around ϵ = 1/n. Remarkably, they still show better performance even when ϵ ≈ 10 -4 /n (Fig. 5 ). These results suggest that local solvers are preferred over non-local ones when the precision requirement is in this range. Comparison with local baselines and CGM. We next compare our three accelerated methods with four baselines. Fig. 5 presents the ℓ 1 -estimation error in terms of the number of operations (quantified as t • vol(S t )) executed. It is evident that our three solvers use significantly fewer operations compared to CGM and the other three local methods. Again, due to maintaining a nondecreasing set of active nodes, ISTA and FISTA require more operations than the locally evolving set process. Ours are more efficient than APPR, where r (0) = e s is used. Efficiency in terms of α and huge-graph tests. We demonstrate the performance of local solvers in terms of different α ranging from 0.005 to 0.25. Interestingly, in Fig. 13 , LOCGD show faster convergence when α is small; this may be because of the advantages of monotonicity properties, which is not present in the accelerated methods. However, in other regions of α, accelerated methods are faster. We also tested local solvers on two large-scale graphs where pa-pers100M has 111M nodes and 1.6B edges while com-friendster has 65M nodes with 1.8B edges. Results are shown in Fig. 6 ; compared with current default local methods, it is several times faster, especially on ogbn-papers100M.Following the experimental setup in Fountoulakis et al. [13] , we consider the task of local clustering on 15 graphs. As partially demonstrated in Tab. 1, compared with APPR and FISTA, LOCSOR uses the least operations and is the fastest, demonstrating the advantages of our proposed local solvers.  
limitations and conclusion: 
 Our proposed algorithms may have the following limitations: 1) When α is small, the acceleration effect partially disappears, as observed in Fig. 13 . This may be due to the limitations of global counterparts, where the residual may not decrease early; 2) Our new accelerated bound for LocCH depends on an empirically reasonable assumption of residual reduction but lacks theoretical justification.We propose using a new locally evolving set process framework to characterize algorithm locality and demonstrate that several standard iterative solvers can be effectively localized, significantly speeding up current local solvers. Our local methods could be efficiently implemented into GPU architecture to accelerate the training of GNNs such as APPNP [27] and PPRGo [7] . We also offer open problems in developing faster local methods. It is worth exploring whether subsampling active nodes stochastically or using different queue strategies (priority rather than FIFO) could help speed up the framework further. It also remains interesting to see how to design local algorithms for conjugate direction-based methods such as CGM.  
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments? Answer: [Yes] Justification: For our experiment, we used a server powered by an Intel(R) Xeon(R) Gold 5218R CPU, which features 40 cores (80 threads). The system is equipped with 256 GB of RAM. Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines? Answer: [Yes] Justification: The research conducted in the paper conforms in every respect with the NeurIPS Code of Ethics. The authors have thoroughly reviewed and adhered to the guidelines, ensuring that all aspects of their work align with ethical standards. Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 10. Broader Impacts Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed? Answer: [Yes] Justification: The advancements in accelerated algorithms like LOCCH and LOCHB can significantly enhance computational efficiency in various applications, contributing to faster and more effective solutions in fields such as data analysis, machine learning, and optimization. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: The adjacency matrix of G.The diagonal degree matrix of G.The normalized Laplacian matrixW ≜ D -1/2 AD -1/2 . V ΛV ⊤ The eigendecomposition of W is given by W = D -1/2 AD -1/2 = V ΛV ⊤ , where Λ = diag(λ 1 , λ 2 , . . . , λ n ), with 1 = λ 1 ≥ λ 2 ≥ • • • ≥ λ n ≥ -1. When λ n = -1, G is a bipartite graph. Q The underlying matrix of Equ. (3) is Q = I -1-α 1+α D -1/2 AD -1/2 . r (t)Given an estimate x (t) , the residual is defined asr (t) ≜ b -Qx (t) . r(t)The D 1/2 -shifted residual is defined asr(t) = D 1/2 r (t) . αThe damping factor α, which lies in the interval (0,1).α A pre-defined constant α = (1 - √ α)/(1 + √ α).T A The total runtime of a local algorithm A.T t (x)For t ≥ 1, T t denotes the Chebyshev polynomial of the first kind, defined asT t+1 (x) = 2xT t (x) -T t-1 (x), with T 0 (x) = 1, and T 1 (x) = x. δ tThe ratio of T t-1 and T t , i.e.,δ t = T t-1 ( 1+α 1-α )/T t ( 1+α 1-α ), and δ t+1 = (2 1+α 1-α -δ t ) -1 , with δ 1 = 1-α 1+α . δ 1:tThe product of all δ t , i.e., δ 1:t ≜ t i=1 δ i . By default, we set δ 0:1 = 0. S 1:tThe intersection of all t-ordered sets, i.e., S 1:t = S 1 ∩ S 2 ∩ • • • ∩ S t . By default, we set S t:t-1 = V. S tThe complement of S t , i.e., S t = V\S t .S j,t S j,t = S j ∩ S j+1 ∩ • • • ∩ S t-1 ∩ S t = S j:t-1 ∩ S t . By default, we set S t,t = S t .A.2 Proof of Lemma 2.1 Lemma 2.1 (Runtime bound of APPR [2] ). Given α ∈ (0, 1) and the precision ϵ ≤ 1/d s for node s ∈ V with p ← 0, z ← e s at the initial, APPR(G, ϵ, α, s) defined in (2) returns an estimate p of π.There exists a real implementation of (2) (e.g., Algo. 2) such that the runtime T APPR satisfiesT APPR ≤ Θ (1/(αϵ)) .Furthermore, the estimate π := p satisfies ∥D -1 ( π -π)∥ ∞ ≤ ϵ and vol(supp( π)) ≤ 2/((1 -α)ϵ).Proof. To find an upper bound of T APPR , we add a time index for all active nodes u 1 , u 2 , . . . , u t processed in APPR of Algo. 2, a real implementation of (2). The parameter t is the number of active nodes processed, and u i is the node dequeued at time i. So, updates of p and z are from (0, e s ) = (p (0) , z (0) ) to (p (t) , z (t) ) as follows:(p (0) , z (0) ) u1 -→ (p (1) , z (1) ) u2 -→ (p (2) , z (2) ) • • • ut -→ (z (t) , p (t) ).For each active u i , the updates of (p, z) by definition can be represented asp (i) = p (i-1) + αz (i-1) ui • e ui , z (i) = z (i-1) - (1 + α)z (i-1) ui 2 e ui + (1 -α)z (i-1) ui 2 AD -1 e ui .Since z (0) = e s ≥ 0, then z (i) ≥ 0 and p (i) ≥ 0 for all i by induction. Note that ∥AD -1 e ui ∥ 1 = 1, we have the following relation from the updates of z∥z (i) ∥ 1 = ∥z (i-1) ∥ 1 -αz (i-1) ui ⇐⇒ z (i-1) ui = ∥z (i-1) ∥ 1 -∥z (i) ∥ 1 α . Note ϵd ui ≤ z (i-1) uifor each active u i . Summing the above equation over i = 1, 2, . . . , t, we havet i=1 ϵd ui ≤ t i=1 z (i-1) ui = t i=1 ∥z (i-1) ∥ 1 -∥z (i) ∥ 1 α = ∥z (0) ∥ 1 -∥z (t) ∥ 1 α ≤ ∥z (0) ∥ 1 α = 1 α ,where note ∥z (0) ∥ 1 = ∥e s ∥ 1 = 1 by the initial condition. Since t i=1 d ui exactly captures the number of operations needed, the runtime of Algo. 2 is then bounded asT APPR = Θ t i=1 d ui ≤ Θ 1 αϵ .To check the quality of estimate p, using the updates of p (i) and summing over all i, we havep (t) = α t i=1 z (i-1) ui e ui = α (1 + α) 2 I - (1 -α) 2 AD -1 -1 Π t i=1 z (i-1) -z (i) = π -Πz (t) ,where Π is the PPR matrix. The above gives us π -p (t) = Π • z (t) . Since G is undirected, the Π matrix satisfies π v [u] = (d u /d v )π u [v] where π v [u] is the u-th element of PPR vector of sourcing node v. Consider each u-th element of Π • z (t) (Π • z (t) ) u = v∈V z (t) v • π v [u] = v∈V z (t) v • d u d v π u [v] ≤ ϵd u v∈V π u [v] = ϵd u ,where the last equality is due tov∈V π u [v] = 1. Hence, (π -p (t) ) u = (Π • z (t) ) u ≤ ϵd u , which indicates ∥D -1 (p (t) -π)∥ ∞ ≤ ϵ.To see the bound of vol(supp(p (t) )), note for any u ∈ supp(p (t) ), it was an active node and there was at least zu (1α)/2 remain in u-th entry of z (t) where we denote zu as the residual before the last push operation of node u; hence u∈supp(p) d u ≤ u∈supp(p) zu ϵ = u∈supp(p) zu (1 -α)/2 ϵ(1 -α)/2 ≤ u∈supp(p) z (t) u ϵ(1 -α)/2 ≤ 2 (1 -α)ϵ . Algo. 2 APPR(α, ϵ, s, G) via FIFO Queue 1: Initialize: p ← 0, z ← e s , Q ← { * , s}, t = -1 2: while true do 3: u ← Q.dequeue() 4: if u == * then 5: if Q = ∅ then 6:p u ← p u + α • z 12: z u ← z • (1 -α)/2 13:for v ∈ N (u) do 14:z v ← z v + (1-α) 2 • z du 15: if z v ≥ ϵd v and v / ∈ Q then 16: Q.enqueue(v) 17:if z u ≥ ϵd u and u / ∈ Q then 18:Q.enqueue(u) 19: return p Indeed, Lemma 2.1 is a special case of Theorem 1 in [2] . The proof outlined above adheres to the key strategy demonstrated in that theorem, which involves exploring the monotonicity and nonnegativity of z. The real implementation of APPR, as shown in Algo. 2, presents a typical queue-based method. It has monotonic properties during the updates of p and z (Lines 10-16 of Algo. 2). It also holds element-wise that p ≥ 0 and z ≥ 0. The operations of Q.enqueue(u), Q.dequeue(), and v / ∈ Q are all in O(1). Line 4 to Line 9 is to design the marker for distinguishing between S t and S t+1 . If all active nodes are processed and no more active nodes are added into Q, then Q will be empty, and finally, the algorithm returns an estimate p of π.Justification of an equivalent condition. We make the justification of an equivalent stop condition for solving (3) . Note we require a local solver to return an estimate π satisfiesEQUATIONSince we define Qx = b asI -1-α 1+α D -1/2 AD -1/2 x = 2α 1+α D -1/2 e s . With r (t) = b -Qx (t)and the stop condition (14) . To see this, since π = D 1/2 x * , we have∥D -1/2 r (t) ∥ ∞ ≤ 2αϵ 1 + α ensures the estimate π = D 1/2 x (t) satisfies∥D -1 ( π -π) ∥ ∞ = ∥D -1/2 (x (t) -x * )∥ ∞ = ∥D -1/2 (Q -1 b -Q -1 r (t) -Q -1 b)∥ ∞ = ∥D -1/2 Q -1 D 1/2 D -1/2 r (t) ∥ ∞ ≤ ∥D -1/2 Q -1 D 1/2 ∥ ∞ • ∥D -1/2 r (t) ∥ ∞ ≤ ∥D -1/2 Q -1 D 1/2 ∥ ∞ • 2αϵ 1 + α where D -1/2 Q -1 D 1/2 = (I -1-α 1+α D -1 A) -1 = ∞ i=0 ( 1-α 1+α D -1 A) i . This leads to ∥D -1 ( π -π) ∥ ∞ ≤ ∥ ∞ i=0 ( 1-α 1+α D -1 A) i ∥ ∞ • 2αϵ 1 + α ≤ 1 + α 2α • 2αϵ 1 + α = ϵ.The Gauss-Seidel Successive Over-Relaxation (GS-SOR) solver (see Section 11.2.7 of Golub & Van Loan [19] ) for the linear system M π = s via the following forward substitution for i in V := {1, 2, . . . , n} do :p (t+1) i = ω   s i - i-1 j=1 M ij p (t+1) j - n j=i+1 M ij p (t) j   /M ii + (1 -ω)p (t) i ,where p is updated from p (t) to p (t+1) . When the relaxation parameter ω = 1, GS-SOR reduces to the standard GS method. Equivalently, let ∆ i = (i -1)/n for i = 1, 2, . . . , n, then GS-SOR updates can be sequentially represented asfor i in V := {1, 2, . . . , n} do : p (t+∆i+1) ← p (t+∆i) + ω M ii   s i - i-1 j=1 M ij p (t+∆i) j - n j=i M ij p (t+∆i) j   • e i .Therefore, it is natural to define the following local variant of GS-SOR. Definition B.1 (Local variant of GS-SOR). Consider the linear system M π = s. For t ≥ 0, we are given an active node set S t = {u 1 , u 2 , . . . , u |St| } and let∆ i = (i -1)/|S t | for i = 1, 2, . . . , |S t |, providing ω ∈ (0, 2), it is natural to define the local variant of GS-SOR as follows:EQUATION)where S t ⊆ V. When ω = 1 and S t = V = {1, 2, . . . , n}, it reduces to the standard GS.We use the above definition to show APPR is a local variant of GS-SOR as the following. Lemma 3.2 (New local evolving-based bound for APPR). Let M = α -1 I -1-α 2 I + AD -1 and s = e s . The linear system M π = s is equivalent to Equ. (1) . Given p (0) = 0, z (0) = e s with ω ∈ (0, 2), the local variant of GS-SOR (15) for M π = s can be formulated asp (t+∆i+1) ← p (t+∆i) + ωz (t+∆i) ui M uiui e ui , z (t+∆i+1) ← z (t+∆i) - ωz (t+∆i) ui M uiui M e ui ,where u i is an active node in S t satisfying z ui ≥ ϵd ui and ∆ i = (i -1)/|S t |. Furthermore, when ω = 1+α 2 , this method reduces to APPR given in (7) , and there exists a real implementation (Aglo. 2) of APPR such that the runtime T APPR is bounded byTAPPR ≤ vol(ST ) αγT ln CT ϵ , where vol(ST ) γT ≤ 1 ϵ , CT = 2 (1 -α)|IT | , γT ≜ 1 T T -1 t=0 |S t | i=1 |z (t+∆ i ) u i | ∥z (t) ∥1 .Proof. Note that M π = s is equivalent to Equ. (1) . We first rewrite p (t+∆i+1) in terms of the residual z. The residual z at time t + ∆ i can be written asz (t+∆i) = s -M p (t+∆i) . Note s ui - i-1 j=1 M uiuj p (t+∆i) uj - n j=i M uiuj p (t+∆i) uj = (s -M p (t+∆i) ) ui = z (t+∆i) ui. Then, the updates of local GS-SOR defined in ( 15) can be rewritten asp (t+∆i+1) = p (t+∆i) + ω Mu i u i z (t+∆i) ui• e ui . Hence, the updates of z (t+∆i) can be written asz (t+∆i+1) = s -M p (t+∆i+1) = s -M p (t+∆i) + ωz (t+∆i) ui M uiui • e ui = z (t+∆i) - ωz (t+∆i) ui M uiui M e ui .Note the diagonal element M uiui = (1 + α)/(2α). Hence, when ω = 1+α 2 , we havep (t+∆i+1) = p (t+∆i) + αz (t+∆i) ui • e ui z (t+∆i+1) = z (t+∆i) -αz (t+∆i) ui M e ui = z (t+∆i) -z (t+∆i) ui ( 1 + α 2 I - 1 -α 2 AD -1 )e ui .The above updates match APPR's evolving set process formulation in (7) . The rest is to show a new runtime bound. Adding ℓ 1 -norm on both sides of the above equation, then note ∥z (t+∆i+1) ∥ 1 = ∥z (t+∆i) ∥ 1αz (t+∆i) ui for i = 1, 2, . . . , |S t |. We have∥z (t+1) ∥ 1 = 1 - α |St| i=1 z (t+∆i) ui ∥z (t) ∥ 1 ∥z (t) ∥ 1 = (1 -αβ t ) ∥z (t) ∥ 1 = t i=0 (1 -αβ t ) ∥z (0) ∥ 1 ,where we defineβ t := |S t | i=1 |z (t+∆ i ) u i | ∥z (t) ∥1 . Let t = T -1, we have ln ∥z (T ) ∥ 1 ∥z (0) ∥ 1 = T -1 t=0 ln (1 -αβ t ) ≤ - T -1 t=0 αβ t ⇒ T ≤ 1 αγ T ln ∥z (0) ∥ 1 ∥z (T ) ∥ 1 ,where the first inequality is due to ln(1 + x) ≤ x for x > -1. For each nonzero node u ∈ I T = {z (T ) u : z (T ) u ̸ = 0, u ∈ V},consider the last time t ′ that it was altered. Then, either the alteration came from u being an active node, with z(t ′ ) u≥ d u ϵ, and after the PUSH operation it became z ϵ; or the alteration came from a neighboring node v u ∈ N (u) pushing its mass onto u, which ensures that z(t ′′ ) u ≥ (1-α)du(t ′′ ) u ≥ (1-α) 2dv u z (t ′ ) vu ≥ (1-α) 2 ϵ. These two cases provide a lower bound of 1-α 2 ϵ. Hence, ∥z (T ) ∥ 1 ≥ ϵ(1-α)|I T | 2, which leads to the corresponding constant C T .To see the lower bound of 1/ϵ, note ϵd ui ≤ z (t+∆i) ui for all i = 1, 2, . . . , |S t |. Then we haveϵ vol(S t ) ≤ |St| i=1 z (t+∆i) ui = β t ∥z (t) ∥ 1 ≤ β t ,where we definedβ t := |S t | i=1 |z (t+∆ i ) u i | ∥z (t) ∥1and the last inequality is due to the monotonic decreasing of∥z (t) ∥ 1 , i.e., 1 ≥ ∥z (0) ∥ 1 ≥ • • • ≥ ∥z (T ) ∥ 1 . Applying the above inequality for all t = 0, 1, 2 . . . , T - 1, it leads to ϵ vol(S t ) ≤ β t ⇒ ϵ T -1 t=0 vol(S t ) ≤ T -1 t=0 β t ⇒ vol(S T ) γT ≤ 1 ϵ ,where the last derivation is from the fact that γT = 1T T -1 t=0 β t := |S t | i=1 |z (t+∆ i ) u i | ∥z (t) ∥1.Remark B.2. The connection between APPR and the Gauss-Seidel is not new [28, 29, 16, 9] . Our work is the first work that has linked APPR and the Gauss-Seidel with a locally evolving set process.In this subsection, recall we defined r(t) = D 1/2 r (t) . Lemma B.3 (Local iteration complexity of LOCSOR (ω ≤ 1)). Denote S t = {u 1 , u 2 , . . . , u |St| } as the active node set at the t-th iteration. When ω ∈ (0, 1], all vectors r(t) ≥ 0 are nonnegative and magnitudes are decreasing ∥r (t+1) ∥ 1 < ∥r (t) ∥ 1 . Let T be the total number of iterations needed. Then, at iteration T , we haveEQUATIONwhere γ t = t -1 t-1 τ =0 γ τ is the mean of active ratio factors defined in Equ. (6).Proof. Recall u i ∈ S t = {u 1 , . . . , u |St| } and∆ i = i-1 |St| , LOCSOR in Algo. 3 updates x (t+∆i+1) = x (t+∆i) + ωr (t+∆i) ui • e ui r (t+∆i+1) = r (t+∆i) -ωr (t+∆i) ui • e ui + (1-α)ω 1+α r (t+∆i) ui • D -1/2 AD -1/2 e uiNote r ≥ 0 during updates when ω ∈ (0, 1] and recall r(t) = D 1/2 r (t) , we have∥r (t+∆i+1) + ωr (t+∆i) ui • e u ∥ 1 = ∥r (t+∆i) + (1-α)ω 1+α r(t+∆i) ui • AD -1 e ui ∥ 1 ∥r (t+∆i+1) ∥ 1 + ωr (t+∆i) ui = ∥r (t+∆i) ∥ 1 + (1-α)ω 1+α r(t+∆i) ui .Summing over the above equations over u i , we haveEQUATIONGiven {x i } T -1 i=0 and x i ∈ (0, 1), the Weierstrass product inequality provides 1 -T -1 i=0 x i ≤ T -1 i=0 (1 -x i ). By using this inequality, we continue to have a lower bound of T as the following1 - T -1 t=0 2αωγ t 1 + α ≤ T -1 t=0 1 - 2αωγ t 1 + α = r(T ) 1 r(0) 1 ⇒ (1 + α) 1 -∥r (T ) ∥ 1 /∥r (0) ∥ 1 2αωγ T ≤ T.To get upper bound of γ t , note each active residual r(t+∆i) ui pushes at most (1-α)ω (1+α) times magnitude to rui+1 , rui+2 , and ru |S t | ; hence, |St| j=i r(t+∆j) uj will increase by at most r(t+∆i) ui• (1-α)ω (1+α) ≤ r(t+∆i)ui in total. Hence, overall u i , we haver(t) St 1 = |St| i=1 r(t) ui ≤ |St| i=1 r(t+∆i) ui ≤ 2 r(t) St 1 .We reach the following lower and upper bounds of γ t , ∥ r(t)S t ∥1 ∥ r(t) ∥1 ≤ γ t := |St| i=1 r(t+∆ i ) u i ∥ r(t) ∥1 ≤ 2∥ r(t) S t ∥1 ∥ r(t) ∥1 .To check the upper bound of T , from Equ. (17) ,∥r (T ) ∥ 1 = T -1 t=0 1 -2αωγt 1+α ∥r (0) ∥ 1 and ln ∥r (T ) ∥ 1 ∥r (0) ∥ 1 = T -1 t=0 ln 1 - 2αωγ t 1 + α ≤ - T -1 t=0 2αωγ t 1 + α ⇒ T ≤ (1 + α) 2αωγ T ln ∥r (0) ∥ 1 ∥r (T ) ∥ 1 ,where the first inequality is due to ln(1 + x) ≤ x for x > -1.Theorem 3.3 (Runtime bound of LOCSOR (ω = 1)). Given the configuration θ = (α, ϵ, s, G) with α ∈ (0, 1) and ϵ ≤ 1/d s and let r (T ) and x (T ) be returned by LOCSOR defined in (8) for solving Equ. (3). There exists a real implementation of (8) such that the runtime T LOCSOR is bounded by1 + α 2 • vol(S T ) αγ T 1 - ∥D 1/2 r (T ) ∥ 1 ∥D 1/2 r (0) ∥ 1 ≤ T LOCSOR ≤ 1 + α 2 • min 1 αϵ , vol(S T ) αγ T ln C ϵwhere vol(S T ) and γ T are defined in (6) andC = 1+α (1-α)|I T | with I T = supp(r (T ) ). Furthermore, vol(S T )/γ T ≤ 1/ϵ and the local estimate π := D 1/2 x (T ) satisfies ∥D -1 ( π -π)∥ ∞ ≤ ϵ.Proof. After the last iteration T , for each nonzero residual r(T ) u ̸ = 0, u ∈ I T , there must be at least one update that happened at node u: Node u has a neighbor v u ∈ N (u), which was active. This neighbor v u pushed some residual(1-α)r (t ′ ) vu (1+α)dv uto u where t ′ < T . Hence, for all u ∈ I T , we have∥r (T ) ∥ 1 = u∈I T r(T ) u ≥ u∈I T (1 -α)r (t ′ ) vu (1 + α)d vu ≥ u∈I T (1 -α)2αϵd vu /(1 + α) (1 + α)d vu = ϵ|I T | 2α(1 -α) (1 + α) 2 ,where the second inequality is because r(t ′ ) vu was active before the push operation. Applying the above lower bound of ∥r (T ) ∥ 1 to Equ. ( 16) of Lemma B.3 and note ∥r (0) ∥ 1 = 2α/(1 + α), we obtain∥r (0) ∥ 1 ∥r (T ) ∥ 1 ≤ ∥r (0) ∥ 1 ϵ|I T | • 2α(1-α) (1+α) 2 = 1 + α ϵ(1 -α)|I T | := C 1 ϵ .The rest is to prove an upper bound 1/(αϵ) of T LOCSOR . Recall that for any active node u, we have residual updates from Algo. 3 as the followingD 1/2 r (t+1) = D 1/2 r (t) -ωr (t) u D 1/2 e u + (1 -α)ωr (t) u 1 + α AD -1 D 1/2 e u .Move -ωr(t)u D 1/2 e u to the left and note ∥AD -1 D 1/2 e u ∥ 1 = √ d u , we then obtain∥D 1/2 r (t+1) ∥ 1 + ω d u r (t) u = ∥D 1/2 r (t) ∥ 1 + (1 -α)ω 1 + α d u r (t) u .Hence, for each active u, we have2αω √ dur (t) u 1+α = ∥D 1/2 r (t) ∥ 1 -∥D 1/2 r (t+1) ∥ 1 .Summing them over all active nodes u and noticing r(t) u ≥ 2αϵ √ d u /(1 + α) by the active condition. Note ω = 1 and ||D 1/2 r (0) || 1 = 2α1+α , we have run time bounded byT LOCSOR = u d u ≤ 1 + α 2α 2 t (∥D 1/2 r (t) ∥ 1 -∥D 1/2 r (t+1) ∥ 1 ) ωϵ ≤ (1 + α) 2αϵ .Combining the above bound and the bound T shown in Lemma B.3, we prove the lower and upper bound of T LocSOR . To check the lower bound of 1/ϵ,i.e., vol(S T )/γ T ≤ 1/ϵ, note2αϵdu i 1+α ≤ r(t+∆i) ui for all i = 1, 2, . . . , |S t |. Then we have 2αϵ 1 + α vol(S t ) ≤ |St| i=1 r(t+∆i) ui = γ t ∥D 1/2 r (t) ∥ 1 ≤ γ t ∥D 1/2 r (0) ∥ 1 = 2αγ t 1 + α ,where the last inequality is due to the monotonic decreasing of ∥D 1/2 r (t) ∥ 1 , i.e., 2α 1+α ≥ ∥D 1/2 r (0) ∥ 1 ≥ • • • ≥ ∥r (T ) ∥ 1 . Applying the above inequality over all t = 0, 1, 2 . . . , T -1, it leads to ϵ vol(S t ) ≤ γ t ⇒ ϵ T -1 t=0 vol(S t ) ≤ T -1 t=0 γ t ⇒ vol(S T ) γ T ≤ 1 ϵ . Algo. 3 LOCSOR(α, ϵ, s, G, ω) via FIFO Queue 1: Initialize: r ← ce s , x ← 0, c = 2α 1+α , t = -1 2: Q ← { * , s} //x u ← x u + ω • r 15: r u ← r u -ω • r 16:for v ∈ N (u) do 17:r v ← r v + (1-α)ω (1+α) • r du 18: if |r v | ≥ c • ϵd v and v / ∈ Q then 19: Q.enqueue(v) 20: if |r u | ≥ c • ϵd u and u / ∈ Q then 21:Q.enqueue(u) 22: return xThe real queue-based implementation of LOC-SOR is presented in Algo. 3. It has monotonic and nonnegative properties during the updates of r ≥ 0 and x ≥ 0 when ω ∈ (0, 1]. Same as APPR, the operations of Q.enqueue(u), Q.dequeue(), and v / ∈ Q are all in O(1). During the updates, one should note that the real vector r presents D 1/2 r (t) while the vector x is D 1/2 x (t) . In this case, the original active node condition is implicitly shifting from|r u | ≥ 2αϵ √ du 1+α to √ d u |r u | ≥ 2αϵdu 1+α .We use this shifted active condition in Lines 11 and 13 and inactive condition in Line 5. When ω ∈ (1, 2), it is possible |r u | < c • ϵd u and LOCSOR will ignore this inactive node u during the updates. This step makes sure S t = {u i : |r(t+∆ i ) u i | ≥ 2αϵ √ du1+α } during the updates.We introduce the following standard result. Lemma B.4 (Young [55] , Section 12.2, Theorem 2.1 ). Given the GS-SOR method for solving Qx = b, if the underlying matrix Q is a Stieltjes matrix and set relaxation parameter ω asEQUATION)where ρ(B) is the largest eigenvalue (in magnitude) ofEQUATIONwhereL ω := (diag(Q) -ωQ L ) -1 (ωQ U -(ω -1)diag(Q)) with Q = diag(Q) -Q U -Q L . Corollary 3.4. Let ω = ω * ≜ 2/(1 + 1 -(1 -α) 2 /(1 + α) 2) and S t = V, ∀t ≥ 0, the global version of LOCSOR has the following convergence boundEQUATIONwhere ϵ t are small positive numbers with lim t→∞ ϵ t = 0.Proof. Recall Q = I -1-α 1+α D -1/2 AD -1/2 and we consider the underlying graph as simple which means A has 0 diagonal. Hence, diag(Q) = I and B is defined asB = 1 -α 1 + α D -1/2 AD -1/2 , ρ(B) = 1 -α 1 + α .Since Q is a Stieltjes matrix, then Lemma B.4 gives a bound on the spectral radius of L ω asρ(L ω * ) ≤ 2 1 + 1 -ρ(B) 2 -1 1/2 = 2(1 + α) 1 + α + 2 √ α -1 1/2 = 1 - √ α 1 + √ α .Recall that Gelfand's formula states [52] :Given spectral radius ρ(L ω * ) := max i∈[n] |λ i (L ω * )|,where λ i (•) is the i-th eigenvalue, there exists a sequence{ϵ t } ≥0 such that ∥L t ω * ∥ 2 = (ρ(L ω * ) + ϵ t ) tand lim t→∞ ϵ t = 0. The standard SOR method is defined as the followingx (t+1) = (diag(Q) -ω * Q L ) -1 ω * b + (ω * Q U -(ω * -1)diag(Q)) x (t) = L ω * x (t) + ω * (diag(Q) -ω * Q L ) -1 b, Note r (t) = Qe (t) = QL t ω * Q -1 Qe (0) = QL t ω * Q -1 r (0) . We have ∥r (t) ∥ 2 = ∥QL t ω * Q -1 r (0) ∥ 2 ≤ ∥Q∥ 2 ∥L t ω * ∥ 2 ∥Q -1 ∥ 2 ∥r (0) ∥ 2 ≤ 2 1 + α • ∥L t ω * ∥ 2 • 1 + α 2α • 2α (1 + α) √ d s = 2∥L t ω * ∥ 2 (1 + α) √ d s .To meet the stop condition, we require |r(t) u | ≤ 2αϵ √ du 1+α . It is enough to make sure ∥r (t) ∥ 2 ≤ 2αϵ (1+α) √ds . This leads to find t such that∥r (t) ∥ 2 ≤ 2∥L t ω * ∥ 2 (1 + α) √ d s ≤ 2αϵ (1 + α) √ d s ⇔ 1 - √ α 1 + √ α + ϵ t t ≤ αϵ.When ϵ t = o( √ α), then the runtime of global LOCSOR is Õ(m/ √ α) where Õ hides log 1 ϵ .B.4 LOCGD and Proof of Theorem 3.5The local gradient descent, namely LOCGD is to usex (t+1) = x (t) + r (t)St and r (t+1) = r (t) -Qr(t)St , where S t = {u i : |r(t+∆i) ui | ≥ 2αϵ √ d u /(1 + α)}where ∆ i = 0. Algo. 4 presents our actual implementation of LOCGD via FIFO Queue.Algo. 4 LOCGD(α, ϵ, s, G) via FIFO Queue 1: Initialize: r ← ce s , x ← 0, c = 2α 1+α 2: Q ← {s} // Assume ϵ ≤ 1 ds 3: t = 0 4: while Q ̸ = ∅ do 5: S t ← [] 6: while Q ̸ = ∅ do 7: u ← Q.dequeue() 8:S t .append((u, r u ))9:x u ← x u + r u 10:r u ← 0 11:for (u, r) ∈ S t do 12:for v ∈ N (u) do 13:r v ← r v + (1-α)r (1+α)du 14: if |r v | ≥ c • ϵd v and v / ∈ Q then 15:Q.enqueue(v)16:t ← t + 1 17: return x, rAlgo. 4 presents LOCGD similar to the real queue-based implementation of LOCSOR. It has monotonic and nonnegative properties during the updates of r ≥ 0 and x ≥ 0. Again, the operations of Q.enqueue(u), Q.dequeue(), and v / ∈ Q are all in O(1). During the updates, one should note that r presents D 1/2 r (t) while x is D 1/2 x (t) . All shifted conditions are the same as of LOCSOR. The key advantage of LOCGD is that it is highly parallelizable, while LOCSOR is truly an online update, so it is hard to parallelize.With the initial x (0) = 0, r (0) = b, S 0 = supp(r (0) ), denote r(t) = D 1/2 r (t) . LOCGD defined in (9) has the following properties: 1) x (t) ≥ 0, r (t) ≥ 0 and ∥r (t) ∥ ≥ ∥r (t+1) ∥ 1 ; 2) The residual and estimation error satisfies∥r (t+1) ∥ 1 = 1 - 2αγ t 1 + α ∥r (t) ∥ 1 , γ t = |St| i=1 r(t+∆i) ui ∥r (t) ∥ 1, where ∆ i = 0.Proof. We first show x (t) ≥ 0, r (t) ≥ 0 are all nonnegative vectors during the updates when b ≥ 0. This can be seen from the induction. At the initial stage, x (0) ≥ 0 and r (0) = b ≥ 0. Now assume that for any t ≥ 0, x (t) ≥ 0 and r (t) ≥ 0. Thenx (t+1) = x (t) + r (t)S k ≥ 0, andr (t+1) = r (t) St + 1-α 1+α D -1/2 AD -1/2 r (t)St ≥ 0. Therefore, x (t) ≥ 0 and r (t) ≥ 0 for all t. Note r(t+1) = r(t)St + 1-α 1+α AD -1 • r(t) St and since ∥AD -1 r(t) St ∥ 1 = ∥r (t)St ∥ 1 , we will have∥r (t+1) ∥ 1 = 1 - 2α 1 + α ∥r (t) St ∥ 1 ∥r (t) ∥ 1 ∥r (t) ∥ 1 , where γ t := ∥r (t) St ∥ 1 /∥r (t) ∥ 1 .Then, we can bound the total residual as the following theorem. Theorem 3.5 (Runtime bound of LOCGD). Given the configuration θ = (α, ϵ, s, G) with α ∈ (0, 1) and ϵ ≤ 1/d s and let r (T ) and x (T ) be returned by LOCGD defined in (9) for solving Equ. (4). There exists a real implementation of (9) such that the runtime T LOCGD is bounded by1 + α 2 • vol(S T ) αγ T 1 - ∥r (T ) ∥ 1 ∥r (0) ∥ 1 ≤ T LOCGD ≤ 1 + α 2 • min 1 αϵ , vol(S T ) αγ T ln C ϵ ,where C = (1 + α)/((1α)|I T |), I T = supp(r (T ) ). Furthermore, vol(S T )/γ T ≤ 1/ϵ and the estimate π := D 1/2 x (T ) satisfies ∥D -1 ( ππ)∥ ∞ ≤ ϵ.Proof. We first show bound 1/(αϵ). We first rearrange r (t+1) = r (t) -Qr(t) St into D 1/2 r (t+1) + D 1/2 r (t) St = D 1/2 r (t) + 1 -α 1 + α AD -1 D 1/2 r (t)St .Note r (t) ≥ 0 and∥AD -1 D 1/2 r (t) St ∥ 1 = ∥D 1/2 r (t)St ∥ 1 . Hence, it leads to∥D 1/2 r (t) St ∥ 1 = 1 + α 2α ∥D 1/2 r (t) ∥ 1 -∥D 1/2 r (t+1) ∥ 1 .At each local iterative t, by the active node condition 2αϵ√ d u /(1 + α) ≤ r (t) u , we have ϵ vol(S t ) = u∈St ϵd u ≤ u∈St (1 + α) √ d u r (t) u 2α = 1 + α 2α D 1/2 r (t) St 1 .Then the total run time of LOCGD presented in Algo. 4 isT -1 t=0 vol(S t ) ≤ 1 ϵ 1 + α 2α 2 ∥D 1/2 r (0) ∥ 1 -∥D 1/2 r (T ) ∥ 1 ≤ 1 + α 2αϵ .Therefore, the total run time is at mostT LOCGD := T -1 t=0 vol(S t ) ≤ 1+α 2αϵ .For estimating the bounds of T , by the Weierstrass product inequality [26] and Lemma B.5, we use the similar argument made in Lemma B.3 and continue to have(1 + α) 2αγ T 1 - ∥r (T ) ∥ 1 ∥r (0) ∥ 1 ≤ T ≤ (1 + α) 2αγ T ln ∥r (0) ∥ 1 ∥r (T ) ∥ 1 ,Note that each nonzero r(T ) u has at least part of the magnitude from the push operation of an active node, say v u at time t ′ < T . This means each nonzero of D 1/2 r (T ) satisfiesr(T ) u ≥ (1 -α)r (t ′ ) vu (1 + α)d vu ≥ (1 -α) • 2αϵd vu /(1 + α) (1 + α)d vu = 2α(1 -α)ϵ (1 + α) 2 , for u ∈ I T .Hence, we have ∥r(T ) ∥ 1 ≥ 2α(1-α)ϵ|I T | (1+α) 2and T is further bounded asT ≤ (1 + α) 2αγ T ln ∥r (0) ∥ 1 2α(1-α)ϵ|I T | (1+α) 2 := (1 + α) 2αγ T ln C T ϵ , where C T = (1 + α) (1 -α)|I T | .The lower bound of 1/ϵ, i.e., vol(S t )/γ T ≤ 1/ϵ, directly follows a similar strategy of previous proof by noticing that ∥D 1/2 r (t) ∥ 1 is monotonically decreasing.Remark B.6. One may consider designing local methods based on Jacobi and Richardson's iterations. Indeed, these two methods have the same updates as standard GD. Recall the standard GD method to solve (4) is x (t+1) = x (t) + r (t) , r (t+1) = r (t) -Qr (t) . The Richardson's iteration is x (t+1) = (I -ωW )x (t) + ωb, i.e., x (t+1) = x (t)ω(W x (t)b). The optimal ω * = 2/(λ min + λ max ) where λ min = 2α/(1 + α) and λ max ≤ 2/(1 + α). Hence one can choose ω = 1 ≤ ω * [19] . It leads to x (t+1) = x (t) + r (t) . One can get the same result for the Jacobi method.We begin by providing the solutions of the second-order nonhomogeneous equation as the following Lemma C.1 (Stević [47] ). The solution of the second-order nonhomogeneous difference equationEQUATION)is characterized by the following two casesx t =      1 λ2-λ1 λt 1 λ2 x 0 -x 1 - t k=1 f k λk 1 + λt 2 x 1 -λ1 x 0 + t k=1 f k λk 2 p 2 ̸ = 4q (-p 2 ) t x 0 - t k=1 kf k (-p/2) k+1 + t(-p 2 ) t-1 x 1 --p 2 x 0 + t k=1 f k (-p/2) k p 2 = 4q, where λ1 , λ2 are two roots of λ 2 + pλ + q = 0, and the summation follows convention0 k=1 • = 0.Based on the above lemma, we have the following corollary Corollary C.2 (Second-order nonhomogeneous equation). Given |a| ≤ 1, the second-order nonhomogeneous equationx t+1 -2ax t + x t-1 = f thas the following solutionEQUATIONProof. The first and last two cases can be directly followed from Lemma C.1. Since λ1 and λ2 are the two complex roots of λ 2 -2aλ + 1 = 0, we write λ1 = re iθ = r(cos θ + i sin θ) and λ2 = re -iθ = r(cos θi sin θ). It indicatesλ 2 -2aλ + 1 = λ -re iθ λ -re -iθ = λ 2 -r(e iθ + e -iθ )λ + r 2 = 0. Since 1 > a 2 , λ1 = a -i √ 1 -a 2, and Re( λ1 ) 2 + Im( λ1 ) 2 = a 2 + (1a 2 ) = 1, then r = 1. Then θ = arccos(Re( λ1 )) = arccos(a), and sin(θ) = Im( λ1 ) = √ 1a 2 . Finally, λt 1 = e itθ = cos(tθ) + i sin(tθ), andλ1 = cos θ + i sin θ, λ2 = cos θ -i sin θ λ1 λ2 = e itθ e -itθ = 1 λ2 -λ1 = -2i sin θ λt 1 = cos(θt) + i sin(θt), λt 2 = cos(θt) -i sin(θt) λt 2 -λt 1 = e -itθ -e itθ = -2i sin(θt). λt-1 1 -λt-1 2 = 2i sin(θ(t -1)) λ1 λ2=1 = λt 1 λ2 -λt 2 λ1Based on these, we haveEQUATIONi -sin(θi(t-1))yEQUATIONFurthermore, when y= λ i y (0) i for i = 1, 2, . . . , n, then solutions can be simplified asEQUATIONProof. The first part is a consequence of Corollary C.2 by letting f t = 0. To see the second identity of this lemma, note that(λ i sin(θ i t) -sin(θ i (t -1))) sin θ i = cos(θ i t).Indeed, by expanding sin(θ i (t -1)), we havesin(θ i (t -1)) = sin(θ i t -θ i ) = sin(θ i t) cos(-θ i ) + cos(θ i t) sin(-θ i ) = sin(θ i t) cos(θ i ) -cos(θ i t) sin(θ i ) = λ i sin(θ i t) -cos(θ i t) sin(θ i )Hence, when yEQUATIONi -sin(θi(t-1))yEQUATIONFurthermore, with initial conditions yEQUATIONi can be simplified asEQUATIONProof. The first part is a consequence of Corollary C.2 by letting f t = h i,t . We prove the second part by considering three cases of λ i as• Case 1. When λ i = 1, we have yEQUATIONi + y (t-1) i = h i,t . For t ≥ 1, the solution the above isy (t) i = y (0) i + t(y (1) i -y (0) i ) + t-1 k=1 (t -k)h i,k = y (0) i + t-1 k=1 (t -k)h i,k ,where the second equation is due to y(1) i = λ 1 y (0) i = y (0) i . • Case 2. When λ n = -1 (G is a bipartite graph), we have y (t+1) i + 2y (t) i + y (t-1) i = h i,t . For t ≥ 1, the solution is y (t) i = (-1) t y (0) i + (-1) t t-1 k=1 (t -k)h i,k (-1) k+1 = (-1) t y (0) i + t-1 k=1 (-1) t-k-1 (t -k)h i,k .• Case 3. When |λ i | < 1, and define θ i = arccos(λ i ). We use a similar argument in Lemma C.3 and continue to havey (t) i = sin(θ i t)y (1) i -sin(θ i (t -1))y (0) i sin θ i + t-1 k=1 sin(θ i (t -k))h i,k sin θ i = (λ i sin(θ i t) -sin(θ i (t -1))) y (0) i sin θ i + t-1 k=1 sin(θ i (t -k))h i,k sin θ i = cos(θ i t)y (0) i + t-1 k=1 sin(θ i (t -k)) sin θ i h i,k .The next lemma presents the properties of Chebyshev polynomials. Lemma C.5 (Chebyshev polynomial bound). For t ≥ 1, the Chebyshev polynomial of the first kind is defined recursively asT t+1 (x) = 2xT t (x) -T t-1 (x) with T 0 (x) = 1, T 1 (x) = x.For t ≥ 1, define δ t = T t-1 ( 1+α 1-α )/T t ( 1+α 1-α ), then 1. T t (x = 1+α 1-α ) and δ t defines the following sequenceδ t+1 = 2 1 + α 1 -α -δ t -1, where δ 1 = 1α 1 + α .2. The closed-form δ 1:t can be upper bounded asδ 1:t = 1 T t ( 1+α 1-α ) = 2 αt + α-t ≤ 2 1 - √ α 1 + √ α t . 3. Note δ 1 = T 0 /T 1 = 1/x = 1-α 1+α, the sequence {δ t } satisfies δ t < 1, ∀t ≥ 1 and 1 = 2δ t+1 xδ t δ t+1 , t = 1, 2, . . . . 4. Denote δ j:t = t i=j δ i for t ≥ j ≥ 0 and set the default value δ j:j-1 = 1 for j ≥ 0, thenδ 1:t /(δ 1:k ) = δ k+1:t ≤ 2 αt-k , for t ≥ k ≥ 0.Proof. For the first item, let x = 1+α 1-α , use the Chebyshev equation, we have1 = 2 1 + α 1 -α T t T t+1 - T t-1 T t+1 = 2 1 + α 1 -α δ t+1 -δ t δ t+1 ⇒ δ -1 t+1 = 2 1 + α 1 -α -δ t .For the second item, for all t ≥ 0, if ξ = x+x -1 2 ̸ = 0, it is well known that the T t can be rewritten asT t ξ = x + x -1 2 = x t + x -t 2 .For our problem, recall we defined α = 1- √ α 1+ √ α and using x = 1+α 1-α ̸ = 0, one can verify thatx = 1 + α 1 -α = α + α-1 2 ⇐⇒ T t x = 1 + α 1 -α = T t α + α-1 2 = αt + α-t 2 and t j=1 δ j = T 0 T 1 • T 1 T 2 • T 2 T 3 • • • T t-1 T t = 1 T t = 2 αt + α-t = 2α t α2t + 1 ≤ 2α t = 2 1 - √ α 1 + √ α t .For the third item, it is sufficient to show that δ t ≤ 1 x for all t ≥ 1, for x = 1+α 1-α . This can be done recursively, since δ 1 = x andδ t+1 = 1 2x -δ t ≤ 1 2x -1 x x≥ 1 x ≤ 1 x .For the last item, note when t ≥ 1 and k ≥ 1, we have the following inequalitiest j=1 δ j • k j=1 δ -1 j = 2 αt + α-t • αk + α-k 2 = αk + α-k αt + α-t = αt-k α2k + 1 α2t + 1 ≤ 2α t-k ,where noteα2k +1 α2t +1 ∈ [1, 2] for t ≥ k.This subsection introduces the standard Chebyshev algorithm. Our following theorem is to prove the runtime complexity of the Chebyshev polynomial iteration for solving Equ. (3). Theorem C.6 (Standard CH). For t ≥ 1, consider the Chebyshev polynomials to solve Equ. (3) asx (t+1) = x (t) + (1 + δ t:t+1 )r (t) + δ t:t+1 x (t)x (t-1) , r (t+1) = 2δ t+1 W r (t)δ t:t+1 r (t-1) , where x (0) = 0, x (1) = x (0) + r (0) and δ t+1 = 2 1+α 1-αδ t -1 with δ 1 = 1-α 1+α . Assume ϵ < 1/d s , then the residual has the following convergence boundr (t) 2 ≤ 2 1 - √ α 1 + √ α t ∥b∥ 2 .Let the estimate be π = D 1/2 x (t) , the the runtime of CH for reaching∥D -1/2 r (t) ∥ ∞ ≤ 2αϵ 1+α with ∥D -1 ( π -π)∥ ∞ ≤ ϵ guarantee is at most T CH ≤ Θ m 1 + √ α 2 √ α ln 2 ϵ = Θ m √ α . Proof. Recall eigendecomposition of W = V ΛV ⊤ where V = [v 1 , v 2 , . . . , v n ]and each v i is the eigenvector. For t ≥ 1, the residual r (t) can be written as n second-order difference equations asV ⊤ r (t+1) -2δ t+1 ΛV ⊤ r (t) + δ t:t+1 V ⊤ r (t-1) = 0,where each i-th element-wise equation of the above can be written as the followingEQUATIONi where v ⊤ i r (0) := y (0) i by default. The above can be rewrittenEQUATIONwhere ⇐⇒ follows from δ 1:t+1 ̸ = 0. Note V ⊤ r (1) = 1-α 1+α V ⊤ V ΛV ⊤ r (0) = δ 1 ΛV ⊤ r (0) , we haveV ⊤ r (0) = y (0) ⇒ V ⊤ r (1) = δ 1 y (1) = δ 1 ΛV ⊤ r (0) = δ 1 Λy (0) ,where it follows from δ 1 ̸ = 0. As y (1) = Λy (0) , follow Equ. ( 24) of Lemma C.3, y (t) i has the solutiony (t) i = y (0) i cos(θ i t) |λ i | < 1 y (0) i λ t i |λ i | = 1 ≤ |y (0) i | |λ i | < 1 |y (0) i | |λ i | = 1, where θ i = arccos(λ i ). We can write down r (t) in terms of y (t) V ⊤ r (t) = δ 1:t y (t) = δ 1:t Z t y (0) , where Z t has two possible formsZ t =diag (1, . . . , cos(θ i t), . . . , (-1) t ) for bipartite graphs diag (1, . . . , cos(θ i t), . . . , cos(θ n t)) for non-bipartite graphs.Hence, ∥Z t ∥ 2 ≤ 1 and ∥Z t y (0) ∥ 2 ≤ ∥y (0) ∥ 2 . We have∥r (t) ∥ 2 = ∥V ⊤ r (t) ∥ 2 ≤ δ 1:t ∥y (0) ∥ 2 ≤ 2 1 - √ α 1 + √ α t ∥b∥ 2 ,where the last inequality follows Lemma C.5 and note z (0) = V ⊤ r (0) . To meet the stop condition of |r(t) u | < 2αϵ √ d u /(1 + α), u ∈ V = ∅, it is sufficient to choose a minimal integer t such that 2 1 - √ α 1 + √ α t ∥b∥ 2 < 2αϵ (1 + α) √ d s .To see this, if the above inequality is satisfied, then note for any node u, we have|r (t) u | ≤ ∥r (t) ∥ ∞ ≤ ∥r (t) ∥ 2 ≤ 2 1 - √ α 1 + √ α t ∥b∥ 2 < 2αϵ (1 + α) √ d s ≤ 2αϵ √ d u (1 + α) ,where all nodes are inactive. So, it gives t ln 1-√ α 1+ √ α ≤ ln ϵ 2 by noticing ∥b∥ 2 = 2α/((1 + α) √ d s ). It indicates t ≥ ln 2 ϵ ln 1+ √ α 1- √ α. As1 ln(1+x) ≤ 1+x x for x > 0, it is sufficient to choose t t = 1 + √ α 2 √ α ln 2 ϵ .We propose the following local Chebyshev iteration procedurex (t+1) = x (t) + (1 + δ t:t+1 )r (t) St + δ t:t+1 x (t) -x (t-1) St .Our next Lemma is to expanding x (t)x (t-1)Lemma C.7. For t ≥ 1 with initials x (0) = 0 and x (1) = x (0) + r (0) S0 , the local Chebyshev iterative is the followingx (t+1) = x (t) + (1 + δ t:t+1 )r (t) St + δ t:t+1 x (t) -x (t-1)St .Denote ∆ (t) = x (t)x (t-1) , then ∆ (t) = t-1 j=0 (1 + δ j:j+1 )t-1 r=j+1 δ r:r+1 r (j) Sj:t-1 , where δ 0:1 = 0, S 0:t = S 0 ∩ S 1 ∩ • • • ∩ S t and δ j:j+1 = δ j δ j+1 . We have the following(1 + δ t:t+1 )r (t) St + δ t:t+1 x (t) -x (t-1) St = t j=0   (1 + δ j:j+1 ) t r=j+1 δ r:r+1 r (j) Sj,t   ,where S j,t ≜ S j:t-1 ∩ S t .Proof. Recall we defined δ 0 = 0 so that δ 0:1 = 0. We prove this lemma by induction.For t = 1, note δ 0:1 = 0 and the support of r (0) is S 0 = supp(r (0) ), thenx (1)x (0) = (1 + δ 0:1 )r (0) .For t = 2, note S 1:1 = S 1 and S 0:1 = S 0 ∩ S 1 by our notation, thenEQUATIONS1:1 + δ 1:2 (x (1) -x (0) ) S1 = (1 + δ 1:2 )r (1) S1:1 + (1 + δ 0:1 )δ 1:2 r (0) S0:1 .For t = 3, one can build x (3)x (2) based on x (2)x (1) and recallS 2:2 = S 2 , S 1:2 = S 1 ∩ S 2 x (3) -x (2) = (1 + δ 2:3 )r (2) S2:2 + δ 2:3 (x (2) -x (1) ) S2 = (1 + δ 2:3 )r (2) S2:2 + δ 2:3 ((1 + δ 1:2 )r (1) S1:1 + (1 + δ 0:1 )δ 1:2 r (0) S0:1 ) S2 = (1 + δ 2:3 )r (2) S2:2 + (1 + δ 1:2 )δ 2:3 r (1) S1:2 + (1 + δ 0:1 )δ 1:2 δ 2:3 r (0) S0:2 .For t = 4, we continue to have x (4) -x (3) = (1 + δ 3:4 )r (3) S3:3 + δ 3:4 (x (3) -x (2) ) S3:3 = (1 + δ 3:4 )r(x (4) -x (3) ) S4 = (1 + δ 3:4 )r (3) S3:3 + δ 3:4 (x (3) -x (2) ) S3:3 S4 . By induction t ≥ 1, x (t) -x (t-1) = t-1 j=0 (1 + δ j:j+1 ) t-1 r=j+1 δ r:r+1 r (j) Sj:t-1 ,where the convention notation 0 i=1 • = 0 and i j=i+1 • = 1. To verify the inductive step, consider for t + 1, we havex (t+1) -x (t) = (1 + δ t:t+1 )r (t) St + δ t:t+1 (x (t) -x (t-1) ) St = (1 + δ t:t+1 )r (t) St + δ t:t+1 (x (t) -x (t-1) ) St = (1 + δ t:t+1 )r (t) St + δ t:t+1   t-1 j=0   (1 + δ j:j+1 ) t-1 r=j+1 δ r:r+1 r (j) Sj:t-1     St = t j=0   (1 + δ j:j+1 ) t r=j+1 δ r:r+1 r (j) Sj:t   .To see the second equation, note(1 + δ t:t+1 )r (t) St + δ t:t+1 x (t) -x (t-1) St = (1 + δ t:t+1 )r (t) St + t-1 j=0   (1 + δ j:j+1 ) t r=j+1 δ r:r+1 r (j) Sj:t-1∩St   = t j=0   (1 + δ j:j+1 ) t r=j+1 δ r:r+1 r (j) Sj,t   ,where recall we denote S j,t ≜ S j:t-1 ∩ S t .Lemma C.8 (Local Chebyshev updates). Given the updates of x (t+1) as defined by LOCCH in (10), we have the following local updatesEQUATIONProof. We only need to show the second equation of (28) . The residual of LOCCH updates isr (t+1) = b -Qx (t+1) r (t+1) = b -Q x (t) + (1 + δ t:t+1 )r (t) St + δ t:t+1 (x (t) -x (t-1) ) St = r (t) -(1 + δ t:t+1 )Qr (t) St -δ t:t+1 Q(x (t) -x (t-1) ) St r (t+1) + (1 + δ t:t+1 )Qr (t) + δ t:t+1 Q(x (t) -x (t-1) ) -r (t) u = (1 + δ t:t+1 )Qr (t) St + δ t:t+1 Q x (t) -x (t-1) St small noisy part . Note Q(x (t) -x (t-1) ) = b -Qx (t-1) -(b -Q(x (t) ) = r (t-1) -r (t) and then u becomes u = r (t+1) + (1 + δ t:t+1 )Qr (t) + δ t:t+1 (r (t-1) -r (t) ) -r (t) = r (t+1) -2δ t+1 W r (t) + δ t:t+1 r (t-1) ,where the last equality is due to (1 + δ t:t+1 )Qr (t) = (1 + δ t:t+1 )r (t) -2δ t:t+1 W r (t) by noticing (1 + δ t δ t+1 ) 1-α 1+α = 2δ t+1 in Lemma C.5. Hence, we have the second equation. To see the noisy part, note by Lemma C.7EQUATIONS0 . The residual r (t) of LOCCH defined in Equ. (28) satisfiesEQUATIONwhereEQUATION1, . . . , cos(θ i t), . . . , (-1) t ) for bipartite graphs diag (1, . . . , cos(θ i t), . . . , cos(θ n t)) for non-bipartite graphs,H k,t =        diag t -k, . . . , sin(θi(t-k)) sin θi , . . . , (-1) t-k-1 (t -k)for bipartite graphs diag tk, . . . , sin(θi(t-k)) sin θi , . . . , sin(θn(t-k))sin θn for non-bipartite graphs.Proof. We first decompose the residual equation in (28) asV ⊤ r (t+1) -2δ t+1 ΛV ⊤ r (t) + δ t:t+1 V ⊤ r (t-1) = t j=0   (1 + δ j:j+1 ) t r=j+1 (δ r:r+1 ) I - 1 -α 1 + α Λ V ⊤ r (j) Sj,t   f (t). Define V ⊤ r (t) = δ 1:t y (t) and V ⊤ r (0) = δ 1:0 y (0) = y (0) by default. Then we haveδ 1:t+1 y (t+1) -2δ 1:t+1 Λy (t) + δ 1:t+1 y (t-1) = f (t) ⇒ y (t+1) -2Λy (t) + y (t-1) = f (t) δ1:t+1 . Note V ⊤ r (1) = δ 1 y (1) = V ⊤ δ 1 W r (0) = δ 1 Λy (0), which indicates y (1) = Λy (0) . Then, by the Lemma C.4, each y (t) i has the solutiony (t) i =          y (0) i + t-1 k=1 (t -k)f (k) i /(δ 1:k+1 ) if λ i = 1 cos(θ i t)y (0) i + t-1 k=1 sin(θi(t-k)) sin θi f (k) i /(δ 1:k+1 ) if |λ i | < 1 (-1) t y (0) i + t-1 k=1 (-1) t-k-1 (t -k)f (k) i /(δ 1:k+1 ) if λ i = -1.Use Z t and H k,t , we write the solution of the second-order difference equation asy (t) = Z t y (0) + t-1 k=1 H k,t f (k) /(δ 1:k+1 ) V ⊤ r (t) = δ 1:t y (t) = δ 1:t Z t V ⊤ r (0) + δ 1:t t-1 k=1 H k,t f (k) /(δ 1:k+1 ) V ⊤ r (t) = δ 1:t Z t V ⊤ r (0) + t-1 k=1 δ k+2:t H k,t k j=0   (1 + δ j:j+1 ) k r=j+1 (δ r:r+1 ) I - 1 -α 1 + α Λ V ⊤ r (j) S j,k   . Note 1 + δ j:j+1 = 2δ j+1 1+α 1-α , j ≥ 1, then (1 + δ j:j+1 ) k r=j+1 (δ r:r+1 ) = 2 1+α1-α δ j+1:k δ j+1:k+1 . Then, we haveV ⊤ r (t) = δ 1:t Z t V ⊤ r (0) + t-1 k=1 δ k+2:t H k,t δ 2:k δ 1:k+1 I - 1 -α 1 + α Λ V ⊤ r (0) S 0,k + 2 t-1 k=1 δ k+2:t H k,t k j=1 δ j+1:k δ j+1:k+1 1 + α 1 -α I -Λ V ⊤ r (j) S j,k = δ 1:t Z t V ⊤ r (0) + δ 1:t t-1 k=1 δ 2:k H k,t I - 1 -α 1 + α Λ V ⊤ r (0) S 0,k + 2 t-1 k=1 δ k+2:t H k,t k j=1 δ j+1:k δ j+1:k+1 1 + α 1 -α I -Λ V ⊤ r (j) S j,k ,where the last term can be expanded as t-1 k=1 δ k+2:t H k,t k j=1 δ j+1:k δ j+1:k+1 1 + α 1 -α I -Λ V ⊤ r (j) S j,k w (j) S j,k = δ 3:t H 1,t δ 2:1 δ 2:2 w (1) S1,1 + δ 4:t H 2,t δ 2:2 δ 2:3 w (1) S1,2 + δ 4:t H 2,t δ 3:2 δ 3:3 u (2) S2,2 + δ 5:t H 3,δ k+1:j H j,t 1 + α 1 -α I -Λ V ⊤ r (k) S k,j .Here, we denote δ k+1:k = 1. The final iterative update is (1 + β j )y t , where y t+1 -2y t + y t-1V ⊤ r (t) = δ 1:t Z t V ⊤ r (0) + δ 1:t t t-1 j=1 δ 2:j H j,t I - 1 -α 1 + α Λ V ⊤ r (0) S0,j t u0,t + 2 t-1 k=1 δ k+1:t (t -k) t-1 j=k δ k+1:j H j,t 1 + α 1 -α I -Λ V ⊤ r (k) S k,j (t -k) u k,t . C.EQUATION)wherey 0 = y 1 = ∥r 0 ∥ 2 .Proof. Since ∥u k,t ∥ 2 ≤ β k ∥r (k) ∥ 2 , the final iterative updates ( 29) can be bounded asEQUATION)where t = 0, 1, . . . , T . These T + 1 (including a trivial one where ∥r (0) ∥ 2 ≤ ∥r (0) ∥ 2 ) inequalities shown in Equ. (31) form a system of linear inequality matrix as            1 0 0 • • • 0 -z 21 1 0 • • • 0 -z 31 -z 32 1 • • • 0 . . . . . . . . . . . . . . . -z T 1 -z T 2 -z T 3 • • • 1             I-Z L             ∥r (1) ∥ 2 ∥r (2) ∥ 2 ∥r (3) ∥ 2 . . . ∥r (T -1) ∥ 2             ≤             δ 1:1 (1 + 1β 0 )∥r (0) ∥ 2 δ 1:2 (1 + 2β 0 )∥r (0) ∥ 2 δ 1:3 (1 + 3β 0 )∥r (0) ∥ 2 . . . δ 1:T (1 + T β 0 t)∥r (0) ∥ 2             := c,where (Z L ) tk = 2δ k+1:t (tk)β k for t = 2, 3, . . . , T and k = 1, 2, . . . , t -1. Assume that N ∈ R T ×T is a strictly lower triangular matrix, then we know the established formula(I + N ) -1 = I + T -1 k=1 (-1) k N k .Hence, we have the following(I -Z L ) -1 = I + T -1 k=1 Z k L .Given that (I -Z L ) -1 ≥ 0, then we obtain an upper-bound of         ∥r (1) ∥ 2 ∥r (2) ∥ 2 . . . ∥r (T ) ∥ 2          ≤ z ≜ (I -Z L ) -1 c.It leads to the following new second-order difference equationz t -2 t-1 k=1 δ k+1:t (t -k)β k z k = δ 1:t (1 + tβ 0 )z 0 , for t = 1, 2, . . . , T,where the initial value of z 0 = ∥r (0) ∥ 2 . Following the argument in Theorem 1 of Golub & Overton [18] , we construct a second-order homogeneous equation for z t asEQUATIONwhere Equ. ( 32) is obtained by the difference between the second equation and the third equation. Similarly,EQUATIONwhere Equ. ( 33) is obtained by the difference of the first two. Hence, Equ. ( 32) - (33) gives usz t+1 -2(1 + β t )δ t+1 z t + δ t:t+1 z t-1 = 0, for t = 1, 2, . . . , T,where two initials are z 0 = ∥r (0) ∥ 0 andz 1 = δ 1 (1 + β 0 )∥r (0) ∥ 2 . Let z t = δ 1:t ẑt , then ẑt+1 -2(1 + β t )ẑ t + ẑt-1 = 0,where two initials are ẑ0 = z 0 = ∥r (0) ∥ 0 and ẑ1 = (1 + β 0 )∥r (0) ∥ 2 . We finish the proof by setting t-1 j=0 (1 + β j )y t = ẑt . Remark C.10. Key points of the above proof strategy largely follow from Golub & Overton [18] . However, different from the original technique, we generalize the strategy to a parameterized version. Lemma C.11. Given β j ≥ 0, the following second-order difference equationx t+1 -2(1 + β t )x t + x t-1 = 0.has the following solutionx t = t-1 j=0 (1 + β j )y t ,where y t+1 -2y t + yt-1 (1+βt-1)(1+βt) = 0 with y 0 = x 0 andy 1 = x 1 /(1 + β 0 ). Proof. Assume x t = -1 2 t t-1 j=0 (-2(1 + β j )) y t .Then, following the equation, we have- 1 2 t+1 t j=0 (-2(1 + β j )) y t+1 -2(1 + β t ) - 1 2 t t-1 j=0 (-2(1 + β j )) y t + - 1 2 t-1 t-2 j=0 (-2(1 + β j )) y t-1 = 0.Since β j ≥ 0, the term t j=0 (-2(1 + β j )) ̸ = 0, we divide it on both sides to have- 1 2 t+1 y t+1 + - 1 2 t y t + - 1 2 t-1 1 4(1 + β t-1 )(1 + β t ) y t-1 = 0.Hence, it is simplified into y t+1 -2y t + yt-1(1+βt-1)(1+βt) = 0.To make a simplification on x t , we prove the lemma. Theorem 4.2 (Runtime bound of LOCCH). Given the configuration θ = (α, ϵ, s, G) with α ∈ (0, 1) and ϵ ≤ 1/d s and let r (T ) and x (T ) be returned by LOCCH defined in (10) for solving Equ. (3). For t ≥ 1, the residual magnitude ∥r (t) ∥ 2 has the following convergence bound∥r (t) ∥ 2 ≤ δ 1:t t-1 j=0 (1 + β j )y t ,where y t is a sequence of positive numbers solving y t+1 -2y t + y t-1 /((1+ β t-1 )(1 + β t )) = 0 with y 0 = y 1 = ∥r (0) ∥ 2 . Suppose the geometric mean β t ≜ ( t-1 j=0 (1 + β j )) 1/t of β t be such that β t = 1 + c √ α 1-√ α where c ∈ [0, 2). There exists a real implementation of (9) such that the runtime T LOCCH is bounded byT LOCCH ≤ Θ (1 + √ α)vol(S T ) √ α(2 -c) ln 2y T ϵ .Proof. The convergence bound of r (t) directly follows from Corollary C.9. Since we assume that there exists c ∈ [0, 2) such thatt-1 j=0 (1 + β j ) ≤ 1 + c √ α 1- √ α t. Then multiplying both sides by αt , we haveαt t-1 j=0 (1 + β j ) ≤ 1 - (2 -c) √ α 1 + √ α t .Then we have∥r (t) ∥ 2 ≤ δ 1:t t-1 j=0 (1 + β j )y t δ1:t≤2 αt ≤ 2 αt βt t y t ≤ ϵ t ln   1 - √ α 1 + √ α t-1 j=0 (1 + β j ) 1/t   ≤ ln ϵ 2y t t ≥ ln 2y t ϵ ln 1 + √ α 1 - √ α β t Since β t = t-1 j=0 (1 + β j ) 1/t ,x ≥ 1 ln(1+x) and lettingx = 1+ √ α 1- √ α β t -1 > 0,then t can be lower bounded further byt ≥ 1 + √ α 1 + √ α -(1 - √ α)β t ln 2y t ϵ ≥ ln 2y t ϵ ln 1 + √ α 1 - √ α β t .Since we assumedβ t = (1 + c √ α 1- √ α ), which means 1 ≤ βt = (1 + c √ α 1- √ α ), so β t ∈ 1, 1+ √ α 1-√ α . Then, we find such an upper bound of t so that LOCCH converges.t = 1 + √ α 1 + √ α -(1 - √ α)β t ln 2y t ϵ = 1 + √ α (2 -c) √ α ln 2y t ϵ .We present the implementation of LOCCH as follows: Recall the sequenceδ t+1 = 2 1+α 1-α - δ t -1 , t = 1, 2, . . . with δ 1 = 1-α 1+α . Denote x(t) ≜ x (t) -x (t-1) , ∆ (t) := (1 + δ t:t+1 )r (t) + δ t:t+1x(t) , we havex (t+1) = x (t) + (1 + δ t:t+1 )r (t) St + δ t:t+1 x(t) St = x (t) + ∆ (t) St r (t+1) = b -Q x (t) + (1 + δ t:t+1 )r (t) St + δ t:t+1 x(t) St = r (t) -Q∆ (t) St x(t+1) = x(t) + ∆ (t) St -∆ (t-1)St-1 .• When t = 0, we haveEQUATION• When t ≥ 1, we can recursively calculate the following vectorsx (t+1) = x (t) + ∆ (t) St r (t+1) = r (t) -∆ (t) St + 1 -α 1 + α W ∆ (t) St x(t+1) = x(t) + ∆ (t) St -∆ (t-1)St-1 .Therefore, at per-iteration, we only need to save sub-vectors ∆ St and ∆ St-1 and update x locally. Lemma D.1 (The standard HB updates). The updates x (t) and r (t) of the HB method for solving Equ. (4) can be written as 1) .x (t+1) = x (t) + (1 + α2 )r (t) + α2 x (t) -x (t-1) r (t+1) = 2αW r (t) -α2 r (t-The residual updates can be rewritten as a second-order homogeneous equationy (t+1) -2Λy (t) + y (t-1) = 0, ∀t = 1, 2, 3, . . . where y (t) is such that r (t) = αt V y (t) , t ≥ 0 with y (0) = V ⊤ r (0) = V ⊤ b.Proof. We follow the standard Polyak's heavy-ball method [40] asx (t+1) = x (t) -η α ∇f (x (t) ) + η β (x (t) -x (t-1) ),where ∇f (x (t) ) = Qx (t) -b andη α = 4/( 2/(1 + α)+ 2α/(1 + α)) 2 = 2(1+α)/(1+ √ α) 2 = 1 + α2 and η β = ( 2/(1 + α) -2α/(1 + α)) 2 /( 2/(1 + α) + 2α/(1 + α)) 2 = α2. Hence, it leads to the following updates 1) .x (t+1) = x (t) + (1 + α2 )r (t) + α2 x (t) -x (t-Insertingr (t) = Q(x * -x (t) ) = b -I - 1 -α 1 + α W x (t) = b -I - 2 α 1 + α2 W x (t)then x (t+1) = 2αW x (t) -α2 x (t-1) + (1 + α2 )b and sinceQW = I - 1 -α 1 + α W W = W I - 1 -α 1 + α W = W Q So r (t+1) = -Q(x (t+1) -x * ) = -2αQW x (t) + α2 Qx (t-1) + (I -(1 + α2 )Q)b = -2αQW x (t) + α2 Qx (t-1) + -α 2 + 2 αW b = 2αW r (t) -α2 r (t-1)Using 1) .r (t) = αt V y (t) , t ≥ 0 αt+1 V y (t+1) = 2αW αt V y (t) -α2 αt-1 V y (t-1) ⇒ V y (t+1) = 2W V y (t) -V y (t-As W = V ΛV ⊤ and V ⊤ = V -1 is orthogonal matrix, we continue to haveV ⊤ V y (t+1) -2V ⊤ V ΛV ⊤ V y (t) + V ⊤ V y (t-1) = 0 ⇒ y (t+1) -2Λy (t) + y (t-1) = 0.Theorem D.2 (Convergence analysis of Heavy-Ball (HB)). To solve the minimization problem in Equ. (4), we propose the following standard HB updates asx (t+1) = x (t) + (1 + α2 )r (t) + α2 x (t) -x (t-1) , r (t+1) = 2αW r (t) -α2 r (t-1) ,where the initial condition is x (0) = 0, r (0) = b, x (1) = x (0) + Γr (0) , r (1) = b -Qx (1) . Then there exists a constant τ such that the total iteration complexity to reach the stop condition {u :|r u | ≤ ϵd u , u ∈ V} = ∅ is t = 1 + √ α 2 √ α ln C t ∥r (0) ∥ 2 ϵ ,where 1) . By Lemma D.1, we have y (t+1) -2Λy (t) + y (t-1) = 0, where we obtained n second-order difference equationsC t = 1 if Γ = Q -1 (I -αW ) (ideal case); C t = max 1+ α-1 √ 1-λ 2 2 , 1 + (1 + α-1 )t if Γ = 0 (practical case); and C t = 2 √ 1-λ 2 2 if Γ = (1-α)(1+α) 2 I (G is not bi-partite graph). Proof. Recall W = V ΛV ⊤ and then V ⊤ r (t+1) = 2αΛV ⊤ r (t) -α2 V ⊤ r (t-y (t+1) i -2λ i y (t) i + y (t-1) i = 0, ∀i = 1, 2, . . . , n.Follow the Lemma C.3, Equ. ( 23) has the solutionEQUATIONi -sin(θi(t-1))yEQUATIONi + (yiλ i yEQUATIONwhere in the case of |λ i | < 1. We consider the three cases of Γ• Ideal Case: We can eliminate t in (34), when y (1) = Λy (0) , we get yEQUATIONi , and then y (t) i can be simplified intoy (t) i = (λi sin(θit)-sin(θi(t-1)))y (0) i sin θi |λ i | < 1 y (0) i λ t i |λ i | = 1 = y (0) i cos(θ i t) |λ i | < 1 y (0) i λ t i |λ i | = 1 ≤ |y (0) i | |λ i | < 1 |y (0) i | |λ i | = 1 .In this case, Γ needs to be Γ = Q -1 (I -αW ). Therefore, we have∥V ⊤ r (t) ∥ 2 = ∥r (t) ∥ 2 = αt ∥y (t) ∥ 2 ≤ αt ∥y (0) ∥ 2 = αt ∥V ⊤ r (0) ∥ 2 = αt ∥r (0) ∥ 2 .• Practical Case: Just letting x (1) = x (0) = 0, we have αy (1) = y (0) , theny (t) i =    α-1 sin(θit)-sin(θi(t-1)) sin(θi) y (0) i |λ i | < 1 (1 + ( α-1 -λ i )t)y (0) i λ t i |λ i | = 1 ≤ max 1 + α-1 1 -λ 2 2 , 1 + (1 + α-1 )t |y (0) i |,where θ i = arccos(λ i ).• Non-bipartite graph Case: When the graph is non-bipartite, we can eliminate t, as the following:We choose Γ = τ I, we havex (1) = x (0) + Γr (0) = Γr (0) , r (1) = b -Qx (1) = r (0) -τ Qr (0) V ⊤ r (1) = (1 -τ )V ⊤ r (0) + (1 -α)τ 1 + α ΛV ⊤ r (0) , v ⊤ i r (1) = (1 -τ + (1 -α)τ 1 + α λ i )v ⊤ i r (0)We have the following relationsv ⊤ i r (0) = y (0) i v ⊤ i r (1) = αy (1) i = (1 -τ + (1 -α)τ 1 + α λ i )v ⊤ i r (0) = (1 -τ + (1 -α)τ λ i 1 + α )y (0) i ,To make t disappear when λ i = 1, we need yEQUATIONi , or1 -τ + (1 -α)τ 1 + α = 1 - √ α 1 + √ α ⇐⇒ τ = 1 + α α + √ α .In this case, we have|y (t) i | ≤ 2 1 -λ 2 2 |y (0) i |To make sure the algorithm stops when the stop condition is met, it is enough for∥r (t) ∥ 2 = αt ∥y (t) ∥ 2 ≤ αt C t ∥y (0) ∥ 2 = αt C t ∥V ⊤ r (0) ∥ 2 = αt C t ∥r (0) ∥ 2 ≤ ϵ.This means C t ∥r (0) ∥ 2 αt ≤ ϵ, which leads to the followingt = 1 + √ α 2 √ α ln C t ∥r (0) ∥ 2 ϵ .Remark D.3. The constant that appears in the bound involves the second largest eigenvalue λ 2 of AD -1 . It is deeply related to the mixing time of random walk [8] where the second largest eigenvalue determines the mixing time of the walk. A smaller absolute value of the second largest eigenvalue indicates that a random walk on the graph will mix (i.e., approach its steady-state distribution) more quickly. Our proof is partially inspired by d'Aspremont et al. [11] where we directly bound r (t) instead of providing bound for e (t) .Lemma D.4. Let the local heavy-ball method be defined as 1) , where x (0) = 0, x (1) = Γr (0) and Γ = diag(Γ 1 , Γ 2 , . . . , Γ n ) is initial step size matrix. We have the following expanding sequencex (t+1) = x (t) + ∆ (t) St , r (t+1) = r (t) -Q∆ (t) St , ∆ (t) = (1 + α2 )r (t) + α2 x (t) -x (t-α2 (x (t) -x (t-1) ) St = (1 + α2 ) t-1 i=1 α2(t-i) r (i) Si:t-1∩St + α2t Γr (0) S0:t-1∩St , ∀t ≥ 1 ∆ (t) St = (1 + α2 ) t i=1 α2(t-i) r (i) Si:t-1∩St + α2t Γr (0) S0:t-1∩St , ∀t ≥ 1.Furthermore, we have the following sequenceEQUATIONwhere we denote S i,t ≜ S i:t-1 ∩ S t .Proof. We assume all nonzeros in b are active nodes at time t = 0 and t = 1, i.e., S 0 = r (0) = supp(b). The local updates can be expressed asx (t+1) = x (t) + (1 + α2 )r (t) St + α2 x (t) -x (t-1) St r (t+1) = b -Qx (t+1) = r (t) -(1 + α2 )Qr (t) -α2 Q x (t) -x (t-1)original updates+ (1 + α2 )Qr (t) St + α2 Q x (t) -x (t-1)St noisy with small magnitudes = 2αW r (t) -α2 r (t-1) + (1 + α2 )Qr(t) St + α2 Q x (t) -x (t-1)St noisy with small magnitudes= 2αW r (t) -α2 r (t-1) + (1 + α2 ) I - 1 -α 1 + α W r (t) St + α2 Q x (t) -x (t-1) St r (t+1) -2αW r (t) + α2 r (t-1) = (1 + α2 )r (t) St -2αW r (t) St + α2 Q x (t) -x (t-1)For t ≥ 1, we can expand x (t+1)x (t) as the followingx (t+1) = x (t) + (1 + α2 )r (t) St + α2 x (t) -x (t-1) St ∆ (t) = x (t+1) -x (t) = (1 + α2 )r (t) St + α2 (1 + α2 )r (t-1) St-1 + α2 x (t-1) -x (t-2) St-1 St = (1 + α2 )r (t) St + α2 (1 + α2 )r (t-1)St-1:t + α4 x (t-1)x (t-2)St-1:t= (1 + α2 )r (t) St + α2 (1 + α2 )r (t-1)St-1:t + α4 (1 + α2 )rEQUATIONEQUATIONS0:t , we continue to haveEQUATIONEQUATIONS0:t-1 , ∀t ≥ 1 α2 (x (t) -x (t-1) ) St = (1 + α2 ) t-1 i=1 α2(t-i) r (i) Si:t-1∩St + Γα 2t r (0) S0:t-1∩St , ∀t ≥ 1The rest follows readily.Lemma D.5 (The nonhomogeneous difference equation). Given y (1) = Λy (0) , equationsy (t+1) -2Λy (t) + y (t-1) := f (t) .have the following solutionsy (t) = Z t y (0) + t-1 k=1 H k,t f (k) ,whereZ t = diag (1, . . . , cos(θ i t), . . . , (-1) t )for bipartite graphs diag (1, . . . , cos(θ i t), . . . , cos(θ n t)) for non-bipartite graphs,H k,t =      diag t -k, . . . , sin(θi(t-k)) sin θi , . . . , (-1) t-k-1 (t -k) for bipartite graphs diag t -k, . . . , sin(θi(t-k))sin θi , . . . , sin(θn(t-k))sin θn for non-bipartite graphs.Proof. This directly follows from Lemma C.4.Theorem D.6 (Representation of r (t) for LOCHB). Given t ≥ 1, x (0) = 0 and x (1) = Γr (0) S0 . The residual of r (t) of LOCHB satisfiesEQUATIONProof. Follow Lemma D.4, we haver (t+1) -Q∆ (t) St = r (t) -Q∆ (t) = r (t) -(1 + α2 )Qr (t) -α2 Q(x (t) -x (t-1) ) = -α 2 r (t) + 2 αW r (t) -α2 Q(x (t) -x (t-1) ) = 2αW r (t) -α2 r (t-1)So we can write the updates of LOCHB asEQUATIONWrite V ⊤ r (t) = αt y (t) , and noteEQUATIONSt = αt V y (t) St ⇒ V ⊤ r (t) St = αt V ⊤ V y (t) St .Hence,V ⊤ r (t+1) -2αV ⊤ W r (t) + α2 V ⊤ r (t-1) = (1 + α2 ) t i=1 α2(t-i) V ⊤ Qr (i) Si,t + α2t V ⊤ QΓr (0) S0,t ⇐⇒ αt+1 y (t+1) -2α αt W y (t) + α2 αt-1 y (t-1) = (1 + α2 ) t i=1 α2(t-i) αi QV ⊤ (V ⊤ y (i) ) Si,t + α2t QΓV ⊤ (V y (0) ) S0,t ⇐⇒ y (t+1) -2W y (t) + y (t-1) = f (t) t i=1 (1 + α2 ) αt-1 α-i QV ⊤ (V y (i) ) Si,t f (t) i + αt-1 QΓV ⊤ (V y (0) ) S0,t f (t) 0 .Then, from Lemma D.5y (t) = Z t y (0) + t-1 k=1 H k,t f (k) ⇐⇒ V ⊤ r (t) = αt Z t V ⊤ r (0) + αt t-1 k=1 H k,t f (k)so expanding the error term EQUATIONβ k,t ≜ ∥u k,t ∥ 2 /∥r (k) ∥ 2 , β k ≜ max t β k,t .Then the upper bound of ∥r (t) ∥ 2 can be characterized asEQUATIONwhere y t+1 -2y t + y t-1 /((1 + β t-1 )(1 + β t )) = 0 where y 0 = y 1 = ∥r (0) ∥ 2 .Proof. Since ∥u k,t ∥ 2 ≤ β k ∥r (k) ∥ 2 , then given the final iterative updates (29)V ⊤ r (t) = αt Z t V ⊤ r (0) + αt tu 0,t + 2 t-1 k=1 αt-k (t -k)u k,tand since ∥Z t ∥ 2 ≤ 1 we can boundEQUATIONwhere t = 0, 1, . . . , T . The rest just follows a similar strategy shown in Corollary C.9. We have the following inequalitiesI -V L :=             1 0 0 • • • 0 -v 21 1 0 • • • 0 -v 31 -v 32 1 • • • 0 . . . . . . . . . . . . . . . -v n1 -v n2 -v n3 • • • 1                         ∥r (1) ∥ 2 ∥r (2) ∥ 2 ∥r (3) ∥ 2 . . . ∥r (T ) ∥ 2             ≤             α1 (1 + β 0 )∥r (0) ∥ 2 α2 (1 + 2β 0 )∥r (0) ∥ 2 α3 (1 + 3β 0 )∥r (0) ∥ 2 . . . αT (1 + T β 0 )∥r (0) ∥ 2             := c,where (V L ) tk = 2α t-k (tk)β k . Denote each upper bound as τ t = ∥r (t) ∥ 2 , we will haveτ t = c t + t-1 k=1 v t,k = αt (1 + tβ 0 )τ 0 + 2 t-1 k=1 (t -k) αt-k β k τ k τ t+1 = c t+1 + t k=1 v t+1,k = αt+1 (1 + (t + 1)β 0 )τ 0 + 2 t k=1 (t + 1 -k) αt-k+1 β k τ k = αt+1 (1 + (t + 1)β 0 )τ 0 + 2 t-1 k=1 (t -k) αt-k+1 β k τ k + 2 t k=1 αt-k+1 β k τ k ατ t = αt+1 (1 + tβ 0 )τ 0 + 2 t-1 k=1 (t -k)α t-k+1 β k τ k τ t+1 -ατ t = αt+1 β 0 τ 0 + 2 t k=1 αt-k+1 β k τ k τ t-1 = αt-1 (1 + (t -1)β 0 )τ 0 + 2 t-2 k=1 (t -k -1)α t-k-1 β k τ k α2 τ t-1 = αt+1 (1 + (t -1)β 0 )τ 0 + 2 t-1 k=1 (t -k) αt-k+1 β k τ k -2 t-1 k=1 αt-k+1 β k τ k ατ t -α2 τ t-1 = αt+1 β 0 τ 0 + 2 t-1 k=1 αt-k+1 β k τ k . τ t+1 -2ατ t + α2 τ t-1 = 2αβ t τ tThe above analysis finally leads to τ t+1 -2(1 + β t )ατ t + α2 τ t-1 = 0. So for αt t-1 j=0(1 + β j )y t = τ t we have y t+1 -2y t + y t (1 + β t )(1 + β t-1 ) = 0.Following the same strategy in Corollary C.9, we obtain the upper bound.Theorem D.8 (Convergence of LOCHB). Let the geometric mean of β t be β t ≜ t-1 j=0 (1 + β j ) 1/t . Then the upper bound of ∥r (t) ∥ 2 can be characterized asEQUATIONwhere y t+1 -2y t + y t-1 /((1 + β t-1 )(1 + β t )) = 0 where y 0 = y 1 = ∥r (0) ∥ 2 . Assume that there exists a constant c ∈ [0, 2) such thatβ t ≤ 1 + c √ α 1- √ α .Then the total number of iterations can be bounded asT ≤ 1 + √ α (2 -c) √ α ln y t ϵ .Then the total runtime T is bounded byT ≤ Θ (1 + √ α)vol(S T ) (2 -c) √ α ln y t ϵ = O vol(S T ) (2 -c) √ α ,where O hides ln(y t /ϵ).Proof. = b Q = Λ -σD -1/2 AD -1/2 [µ, L] Ref. I -σW x = αe s 1. Λ = I, σ = 1 -α b = αe s [α, 2 -α] [27] [53] αI + 1-α 2 L x = αD -1/2 e s 2. L = I -D -1/2 AD -1/2 Λ = I, σ = 1-α 1+α b = 2αD -1/2 s/(1 + α) 2α 1+α , 2 1+α [2] [13] [12] [37] I -(1 -α)AD -1 y = αe s 3. Λ = I, σ = 1 -α b = αD -1/2 e s , α ∈ (0, 1) x = D -1/2 y [α, 2 -α] [7] [56] λ n I n + D -A y = 2λe s 4. Λ = λ n D -1 + I n , σ = 1 b = 2λD -1/2 e s , λ ∈ [1, n] x = D 1/2 y λ ndmax , λ+2n n [41] [57]We present the implementation of LOCHB as follows: Recall the updates of LOCHB isx (t+1) = x (t) + (1 + α2 )r (t) St + α2 x (t) -x (t-1)St , r (t+1) = 2αW r (t) -α2 r (t-1) . The corresponding local updates are∆ (t) = (1 + α2 )r (t) + α2 x(t) x (t+1) = x (t) + ∆ (t) St r (t+1) = r (t) -∆ (t) + 1 -α 1 + α W ∆ (t)x(t+1) = x(t) + ∆ (t) -∆ (t- 1) , where if we choose x (0) = x (1) = x(1) = 0, r (0) = r (1) = b and ∆ (0) = 0.This section presents most commonly used graph-induced linear system as the followingΛ -σD -1/2 AD -1/2 Q x = b,where Q is the generalized version of the perturbed normalized graph Laplacian matrix with perturbation parameter σ > 0, and b is a sparse vector. A typical example ofQ = I -1-α 1+α D -1/2 AD -1/2 with Λ = I and b = 2αD -1/2 e s /(1 + α)The detailed parameters are:• 1. W = D-1/2 Ã D-1/2 and Ã = I + A is the adjacency matrix defined on G(V, E) by adding self-loops for all nodes, and D = I + D is defined as the augmented degree matrix by adding self-loops. α ∈ (0, 1) and usually α < 0.5, The I -(1α)W is the perturbed augmented normalized Laplacian with perturbed parameter α.• 2. Q = D -1/2 D -1-α 2 (D + A) D -1/2 = αI + 1-α 2 L > 0 and Q = 1-α 1+α D -1/2 AD -1/2. This is known as the lazy random-walk version of PPR vectors.• 3. x = α I -(1α)AD -1 -1 e s This is the standard Personalized PageRank vectors widely used for graph embeddings and graph neural network designing [7] . It is also used for decoupling for large-scale GNNs [56] . Graph kernel computation for online learning. Each computed vector serves as semi-supervised learning feature vectors [24] or as online node labeling learning vectors [41, 57] . Note the target linear system when σ = 1D -1/2 λ n I n + D -A D -1/2 D 1/2 y = 2λD -1/2 e s λ n D -1 + I n -D -1/2 AD -1/2 D 1/2 y = 2λD -1/2 e s .Hence, we have Λ = λ n D -1 + I n , σ = 1, b = 2λD -1/2 e s .Following Leskovec et al. [34] , we treat all 17 graphs as undirected with unit weights. We remove selfloops and keep the largest connected component when the graph is disconnected. After preprocessing, the graphs range from 169, 343 nodes in ogbn-proteins to 111, 059, 433 in ogbn-papers100M, as presented in Table 3 .For solving Equation (3), we randomly select 50 source nodes s from each graph. The damping factor is fixed at 0.1, i.e., α = 0.1 for all experiments, and it varies within the range {0.005, 0.01, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3} for others. The ϵ is chosen from the range 2α/((1 + α)d s ), 10 -4 /n .For solving the local clustering problem, we follow the greedy strategy from Andersen et al. [2] , where a local cluster is identified by examining the top magnitudes in PPR vectors. Specifically, we denote the boundary of S as ∂(S) = {(u, v) ∈ E : u ∈ S, v / ∈ S}. The conductance of S is defined asΦ(S) ≜ |∂(S)| min(vol(S), 2m -vol(V\S)). The goal of local clustering is to obtain PPR vectors using these local methods and then apply clustering algorithms to find clusters with low conductance. For the sorting process, given the approximate PPR vector π, we sort D -1/2 π in decreasing order of magnitudes. Let the ordered nodes be v 1 , v 2 , . . . , v t ; the local clustering algorithm iteratively checks the conductance reduction by v 1 , v 2 , • • • , v k where k = 1, 2, . . . , t, and after completing all checks, it returns a subset v 1 , v 2 , . . . , v k ′ that has the minimal conductance among all examined subsets.Parameter settings of baselines. For the local ISTA method [13] , the precision parameter is set to ε = 0.5 for all experiments. According to the algorithm's description of ISTA, the corresponding ρ value is given by ϵ/(1 + ε). For LOCSOR, the parameter ω is calculated as2(1 + α)/(1 + √ α) 2 .For the local FISTA, as demonstrated in [22] , we adopt the same settings as for ISTA and follow its implementation guidelines. We also include preliminary results on ASPR [37] . The algorithm incorporates a parameter, ε, to control the number of iterations in the nested Accelerated Projected Gradient Descent (APGD). We adjust ε from low precision, ε = 0.1/n, to high precision, ε = 10 -4 /n, to ensure the identification of a good approximation.For our experiment, we used a server powered by an Intel(R) Xeon(R) Gold 5218R CPU, which features 40 cores (80 threads). The system is equipped with 256 GB of RAM.In all 15 graphs, we set α = 0.1 and ϵ = 0.1/n. For each of the testing graphs, we randomly select 50 nodes and run LOCSOR and CGM. Fig. 9 presents the missing results on the estimation error reduction for 15 datasets. Compared with the global solver CGM, all local methods show significant speedup in the early stages. To compare our three local solvers, we zoom in on our results and present them in Fig. 10 . Empirically, LOCSOR is the fastest algorithm when the parameter ω is chosen optimally. Many graph applications [2, 7, 29, 15, 27, 36, 30, 44, 25, 45, 51, 48] Local algorithms. Local solvers, in contrast to standard counterparts, leverage the fact that the energy of π lives in a small portion of the graph and hence do not require O(m) or O(n) per iteration. They are advantageous for huge-scale graphs demonstrated in Fig. 15 . Andersen et al. [2] used APPR to obtain an approximate of π for local clustering. Quite similar algorithms were developed in Berkhin [6] (bookmark-coloring algorithm) and Kloster & Gleich [28] (Gauss-Southwell procedure).Under the same stopping condition as APPR, Fountoulakis et al. [13] demonstrated that APPR is equivalent to coordinate descent via variational characterization, with a runtime of Õ(1/(αϵ)) using ISTA where the monotonicity and conservation properties remain. Hence, it is nature to ask whether Õ(1/( √ αϵ)) could be achieved by FISTA [5] in Fountoulakis & Yang [12] . However, the difficulty is that FISTA violates the monotonicity property where the volume accessed of per-iteration cannot be bounded properly. To overcome this, Martínez-Rubio et al. [37] proposed a nested accelerated projected gradient descent (APGD) and gradually expanding solutions so that the monotonicity property still holds. However, nested APGD requires solving subproblems accurately, which in practice may be cumbersome if the precision requirement of the inner problem is too stringent. All current local methods rely on some monotonicity property of variables to guarantee locality, which does not exist in most accelerated frameworks; thus, developing an accelerated method that is guaranteed to preserve intermediate variable sparsity remains challenging.It is worth mentioning that local methods are also closely related to sublinear time and local computational algorithms [42, 1] . From the optimization perspective, the equivalence between Gauss-Seidel and coordinate descent has been considered [49, 35, 39, 50] but does not focus on local analysis.• If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.The local variant of GS-SOR is defined in Appendix B.1See a justification in Appendix B.In practice, TA := T -1 t=0 (vol(St) + |St|) where we ignore |St| for simplicity as vol(St) dominates |St|.Additional experimental results, setups, and algorithm parameters are provided in Appendix F. 


Decision: This paper introduces a novel approach to efficiently compute approximate personalized PageRank (and more generally for sparse linear systems) based on a locally evolving set process. It further develops localized versions of chebyshev and heavy ball methods that achieve a $1/\sqrt{\alpha}$ dependence in the number of iterations necessary for approximating pagerank with dampening factor $\alpha$, in contrast with the $1/\alpha$ dependence for the classic algorithm by Andersen-Chung-Lang. This dependence is optimal, and the authors also present simulations that suggest an empirical advantage for small enough target accuracy.

The reviewers were positive about the paper, highlighting the novel technical contributions. There were some concerns about the empirical advantage being realized only for very small approximations, but I think this is an overall solid contribution.

Reviews: review: paper_summary: paper_summary: This paper considers the approximate personalized page rank. Classical results for this problem have a runtime that is linear in $1/\alpha\epsilon$ where $\alpha$ is the damping factor and $\epsilon$ is the error parameter. The authors show that APPR is simply a local variant of Gauss-Seidel Successive Overrelaxation. Using this connection, the authors derive new run time bounds for APPR and also propose a new algorithm based on Gradient Descent. The execution time for both these are, in the worst-case, identical to the previous bounds. However, they are more sensitive to the state of execution of the algorithms (depend on the active nodes) and seem to mirror the actual performance of these algorithms. Also, under certain assumptions, they improve the worst-case execution time.
strength_weakness: strength_weakness: NA
questions: questions: Can you explain how weak/strong are the assumptions that you make to achieve the improvement for the local Chebyshev method? I couldn't quite gauge the usefulness of this result.
limitations: limitations: NA
paper_summary: paper_summary: This paper uses the evolving set procedure to give a local PageRank algorithm whose dependence on \alpha (the reset probability) is \sqrt{\alpha}.

It proposes accelerated local iterative methods with coefficients given by Chebyshev iteration. The convergence of this algorithm in both graph theoretic and general sparse linear systems settings are analyzed in detail. Discussions of the relations between this method and other local iterative algorithms are also given in detail.

The method was implemented and tested on a range of graphs, mostly coming from social networks. This includes two large ones with edges in the billions. On moderate ranges of \alpha (reset probability), the experiments show significant speedups (factor of about 3) and convergences (factor of 10) on most graphs.
strength_weakness: strength_weakness: The gains only kick in at a relatively large number of steps: it's not clear to me that these are the parameter regimes in which local algorithms actually get used.

Ideally for the empirical works I'd also like to see comparisons of downstream tasks and effects on overall accuracies (e.g. F-1 score), but the paper itself has already covered a lot of ground.
questions: questions: Is the dependence on \epsilon optimal? Aka. have methods with \sqrt{\eps} (or even \log(1 / \eps)) dependences been ruled out?
limitations: limitations: yes, limitations have been addressed, and are entirely theoretical w.r.t. some graph parameter regimes.
paper_summary: paper_summary: This paper considers the study of local algorithms for graph clustering which is an important problem in the field of graph data analysis. In particular this paper is considers the task of computing Personalized Page Rank (PPR) vectors for a given graph. In this problem the algorithm is given a graph in the form of its adjacency and degree matrices, the goal is to approximate the Personalized Page Rank vector for a given starting vertex and dampening factor $\alpha$ up to precision $\epsilon$ without accessing the entire graph. The classical algorithm of Andersen, Chung and Lang runs in time $O(1/\alpha \epsilon)$, which independent of the graph size. The central question posed by subsequent works is whether the dependence on $\alpha$ can be improved to $1/\sqrt{\alpha}$. The main contribution of the paper is to propose a new algorithmic framework based on the locally evolving set process. Under this framework they are able to implement existing algorithms such as Andersen et al.'s APPR algorithm as well as localized implementation of standard gradient descent. They are also able to develop localized versions of chebyshev and heavy ball methods that do achieve the $1/\sqrt{\alpha}$ dependence for some fixed constant value of $\epsilon$. Finally they show that on several large scale graphs, their new localized chebyshev and heavy ball methods do outperform APPR and related methods empirically.
strength_weakness: strength_weakness: One weakness is that the paper is only able to obtain a quadratic improvement in the dependence on the parameter $\alpha$, obtained by the local implementation of the Chebyshev and Heavy-ball method, only for a value of $\epsilon$ and not for all.
questions: questions: One question is that what is the core reason for not obtaining convergence result for accelerated methods for all $\epsilon>0$.
limitations: limitations: Authors have addressed limitations.

Comments: title: Thank all reviewers for their time and effort in the discussion and for helping us.
Response to rebuttal
thank you
comment: We sincerely thank all reviewers for their time and effort in the discussion and for helping us; we will keep improving our submission.
Thanks to the authors for answering the question. My score remains the same.
Thank you for the detailed comments / clarifications.

It's good to know the \epsilon regime being considered here. However, my worry is that in such parameter regimes, non-local things with some dependencies on n might start winning. This is, of course, beyond the scope of this work, or even this discussion. So I'll keep my review/score unchanged.