Iterative Methods via Locally Evolving Set Process
Authors: Baojian Zhou
Yifan Sun
Reza Babanezhad Harikandeh
Xingzhi Guo
Deqing Yang
Yanghua Xiao
Venue: NeurIPS 2024 poster (2024)
Keywords: local computation
Personalized PageRank
graph clustering
Primary area: optimization
TL;DR: This paper introduces a novel approach to efficiently solve sparse linear systems by localizing standard iterative solvers using a locally evolving set process.

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.

Iterative Methods via Locally Evolving Set Process
Baojian Zhou 1.2
Yifan Sun
Reza Babanezhad Harikandeh"
Xingzhi Guo
Deqing Yang!
Yanghua Xiao²
the School of Data Science, Fudan University, Shanghai Key Laboratory of Data Science, School of Computer Science, Fudan University 3 Department of Computer Science, Stony Brook University, 4 Samsung SAIT Al Lab.
Abstract
Given the damping factor and precision tolerance Andersen et al. [2] intro- duced Approximate Personalized PageRank (APPR), the de facto local method for approximating the PPR vector, with runtime bounded by ((1/(ae)) indepen- dent of the graph size. Recently, Fountoulakis & Yang [12] asked whether faster local algorithms could be developed using ((1/(Vac) operations. By noticing that APPR is local variant of Gauss-Seider this paper explores the question whether standard iterative solvers can be effectively localized. We propose the locally evolving set process, novel framework to characterize the al- gorithm locality, and demonstrate that many standard solvers can be effectively localized. Let vol(S,) and 3, the running average of volume and the residual ratio of active nodes S, during the process. We show vol(S)) S 1/c and prove APPR admits a new runtime bound (vvo((())) (ar7,) mirroring the actual perfor- mance. Furthermore, when the geometric mean of residual reduction is @(va), then there exists E (0, 2) such that the local Chebyshev method has runtime (vvol(S)) (//a//a(2222/)) - c))) without the monotonicity assumption. Numerical results confirm the efficiency of this novel framework and show up a hundredfold speedup over corresponding standard solvers on real-world graphs.
Introduction
Personalized PageRank (PPR) vectors are key tools for graph problems such 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 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 graph G. respectively. Given source node s and the damping factor a E (0, 1). this paper studies local solvers for the linear system
(I- (I AD- /2) ae.,
where e, is the standard basis of and is the PPR vector [2, 12, 37]. Given the tolerance 6, local solver needs to find such that D - न T)) T) S without accessing the entire graph G.²
"Corresponding Baojian Zhou, bjzhou @ add.n.eduu.cn Local methods for solving Equ. (1) can be naturally extended to other linear systems defined on G.
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

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/(ae)) independent of graph size. Based on a variational characterization of Equ. (1), Fountoulakis et al. [13] reformulated the problem as optimizing quadratic objective plus l errreuuaarizttion and later asked [12] whether there exists a local solver with runtime ((1/(vac)) This corresponds to an accelerated rate since a is the strongly convex parameter. Recently, Martínez-Rubio et al. [37] provided method based on a nested subspace pursuit strategy, and the corresponding iteration complexity is bounded by [[IS"I/va) where S is the support of the optimal solution. This bound deteriorates to (n/ vα) 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 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- Seiissl Successive Overrelaxation (GS-SOR) that can be treated as an evolving set process." This insight leads us to explore whether standard solvers can be effectively localized solve Equ. (I). To develop faster local methods with improved local bounds, we propose novel locally evolving set process framework inspired by its stochastic counterpart [38]. This framework enables the development of faster local methods and eircumvents 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 core tool, 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/(ac)) runtime bound.
This local evolving process provides a novel way to characterize the algorithm locality; hence, new runtime bounds can be derived. Let vol(S,) and 7. be the running average of volume and the residual ratio of active nodes S, during the process;
prove the ratio vol(S,) T, serving lower bound of 1 1/c. We further show both APPR and local GD have (vvo((S))/(en) (ar))) runtime bound mirroring the actual performance of these two methods.
Using our framework, we show there exists E (0, 2) such that both the localized Chebyshev and Heavy-Ball methods admit runtime bound ((vol(S))/(Va(2 c))) with the assumption that the geometric mean of active ratio factors is O (( 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 -base solvers on two huge-scale graphs.
Paper structure. We begin by clarifying notations and reviewing APPR in Sec. 2. Sec. 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: // g gighttu thub. com/ hom/baojtan baojian, LocalcaCCH
2 Notations and Preliminaries
Notations. We consider an undirected simple graph ((V, g) where V == = n.. and E C VXV with |8| m are the node and edge sets, respectively. The set of neighbors of v is denoted as
The local variant of GS-SOR is defined in Appendix B.1

---

(v) S V. The adjacency matrix A of G assigns unit weight au,e = if (u, u) & and 0 otherwise. The v-th entry of the degree matrix D is d, [W())). Given S S V, we define the volume of S as vol(S) & Enes dy. The support of 0, V). The eigendecompossition of D ADD 1/2 = VAVT where each column of V is an eigenvector and A diag( (1, ,2,...,,.. An) with X 2 do 2 ... An
IR" the set of nonzero indices supp(æ) & { 2e o
2.1 Revisiting Anderson's 's APPR and its local runtime bound
We use (p, 2) for solving Equ. (1) while use (æ, r) for solving Equ. (3) or Equ. (4). With the initial setting p 0, Ex, APPR obtains a local estimate of denoted as p by using a sequence of PUSH operations defined as
APPR (a, 6, 6, 9) Repeat (p. z) < PUSH(u, a, Until Vo, Fv. edo; Return p.
At each repeat step (p.*) PUSH( u, a, p, x), it synchronously up- dates both p and residual æ whenever there exists an active node u E V (a node with large residual, i..., 2u 2 ed,). Specifically, for each active u, it updates Pu. 2u, and 2. for v E (N()) by using a PUSH operator illustrated 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.
Algo. PUSH( u, Ot, P,
1: V ===
2: Pa 6 Pu +V.
3: f- (1-a) (1-a)./2
the left. It stops when active nodes are left.
4: EN(u) do (1-a)u
5: Zy f- Zu + 2d..
6: Return (p.z) 2)
Lemma 2.1 (Runtime bound APPR (2]). Given a E (0, 1) and the precision S 1/ds for node E V with p & 0, (2) returns an estimate p of π. There exists a real implementation of (2) (e.g,, Algo. 2) such that the runtime TAPPR satisfies
at the initial, APPR(a e, (a,,,,,,,,, 8, defined
TAPPR S (1/(aE))
Furthermore, the estimate क := p satisfies D TT)l S and vol(supp (+)) 2/((1-a)e
The main argument for proving Lemma 2.1 is critically based on: 1) 2 0 and |||||| decreases during the updates; 2) for each active u, Zu 2 edn, implying that ||=||] is decreased by least avedu consecutively leading to Σn d, I 1/(a6).
2.2 Problem reformulation
To approximate the PPR vector æ, the original linear system in Equ. (1) can reformulated as an equivalent symmetric version defined as
Qx b, with Q - I
D 1/2 A/2AD-11 11/2 and A
20
1/2
(1+a)
where again eg is the standard basis of 8, and Q is a symmetric positive-definite M-matrix with all
2a eigenvalues in 140 I+al To solve Equ. (3) is equivalent to solving quadratic problem
arg min =ER=
{(#))
2 TQx x b].
where f is strongly convex with condition number 1/a. Indeed, Equ. (3) is a symmetrized version be recovered from x" by 1/222" It is convenient to denote estimate of π (!) (1)22(())) Given (+()) we define the residual (() & (x If ((() is returned by a local solver for solving either Equ. (3) or Equ. (4), we then equivalently require D-11/2
of Equ. (1) and has unique solution xt Q The PPR vector
1l S 6. Hence, it is enough to have stop condition D (t) I 2 /11 + a:) for local solvers of Equ. (3) and Equ. (4).4
Fountoulakis et al. [13] demonstrated that APPR is equivalent to coordinate descent solver for minimizing f in Equ. (4) and introduced an ISTA-style solver by minimizing f(x) + ea D 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/(ae)) provided in Lemma 2.1 becomes less valuable when e < 1/m; on the other hand, all previous local variants [6, 9, 13, 37]
'See justification in Appendix B.
3

---

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.
3 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 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 (11/(ac)).
3.1 Locally evolving set process
Given a; 6, 8, and g, a local solver for Equ. (3) keeps track of active set St C V at each iteration t. That is, only nodes in St are used to update a: r. The next set S1+1 is determined by current S 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 ê (a, 6, s, 9), and a local iterative method A, the locally evolving set process generates a sequence of (S,, x (0) p())) representing the following dynamic system
(t+1),
,,((11) = tb (0), .A). 1I 2 0,
(5)
where S S S, U (Unes,N(u)) and we denote the active set S, {ui, U2, "Is, J. The set S, is maintained via a queue data structure. We say this process converges when the last set ST = 0 if there exists such T; the generated sequence of active nodes are
(So, (0) (0))
(8,,x1) (()) - x2)
(2))
(S, 0, a
((
The runtime of the local solver, A for this whole local process, is then defined as
TA
vol(S.).
0=
The framework of this set process provides new way design local methods. Furthermore, it helps to analyze the convergence and runtime bound of local solvers by characterizing the sequences (v))(.:)), and [[]r()) [} generated by o To analyze new runtime bound, for T 2 1, we define the average of the volume of active node sets {vol(S.)} and active ratio sequence {7t as
T-1 vol,S,)
7-1
Σ |S
IV du, Tw;
((+A))
(Sr
& T 1=0
(6)
T
t=0
where A; is smaller time magnitude. We define A (i
1)/1Sel for the analysis of APPR and LocSOR while A, 0 for LocGD in our later analysis. In the rest, denote Iy (T))(r(T)).
These two metrics vol(Sr) and Tr characterize the 06
1.0
x10"
1,00
as
CLS
locality of local methods. To demonstrate this local process, Fig. 1 shows vol(S;) of
0.50
1.0
0.4
0.75
0.5
20
40
20
10
70
APPR peaks at the early stage, and the active ratio decreases the active volume diminishes. The quantity vol(Sr)/7) is strictly smaller than 1/e, indicating that it could serve as better factor in the runtime analysis.
Figure I: Runtime of APPR in the locally evolving set process the com-dblp graph with = O, a = 0.1, and e = 1/m. The red region the left figure is TAPPR- The right two figures show active ratios and vol(Sr) /r I 1/e.
'In practice, Ta ET= (vol(S.) + (Se]) where ignore [S.] for simplicity as vol(Se) dominates [St]-

---

3.2 APPR via locally evolving set process
We first demonstrate how this locally evolving set process can represent APPR. For solving Equ. (1), the set So = {s} and the queue-based of APPR (see Algo. 2 in Appendix A) naturally forms sequence of active sets from Su {s} to ST (, hence converging. Active nodes u in queue satisfy Su 2 <dg. To delineate successive iterations S and S++1, one can insert at the beginning of S. After processing St, 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 to consistent with Sec. 2.1 and write out such evolving process as follows
(St, (5, z (()
To (S., 200), A APPR): for th in St == {wI, W2,.... "|s: do p((+Am)) t p (++Ai) + a24
(H+A)
Cu,
D = (i - 1)/(5a]
2((+Ai11) t- ((+A.)) (1+a) (1+A.)
Ca; (1-a) (l+b.) AD-leus
2th
The following lemma establishes the equivalence between APPR and the local variant of GS-SOR method (see Appendix B. and provides a new evolving-based bound.
Lemma 3.2 (New local evolving-based bound for APPR). Let M a a1 +I+ AD-1) and 8 = Cs. The linear system MT π is equivalen to Equ. (1). Given (0)) 0. (0) E, with
E (0, the local variant of GS-SOR (15) for MT can be formulated as
(t+A:)
(++A,)
pl(+Am+) t p(t+A.))
Whn
(1+A4441) + (1+
Win:
M Me
Maga
Bu,
M
where ui is an active node in St satisfying Zu, 2 edu,, and A, = (i - 1)/1S;). Furthermore, when
w = L+a this method reduces to APPR given in (7), and there exists a real implementation (Aglo. 2)
of APPR such that the runtime TAPPR bounded, that is
vol(Sr) aryr
OT
vol(Sr) 3rr
2 [1 (I-a)
TAPPR <
where
S 7: Cr Or=
Pr
T
[+)
3.3 Faster local variant of GS-SOR
Lemma 3.2 points to the sub-optimality of APPR, as GS-SOR allows for a larger w. For solving Equ. (3), since APPR essentially serves local variant of GS-SOR, we can develop a faster local variant based SOR. To extend this method solve Equ. (3), we propose a local GS-SOR based on an evolving set process, namely LocSOR, as the following
Da (8,, p(t), A LOCSOR) for uy in Se {ui, .u.,..., l|s.i and do ((t+A + (t+Ai) + wTu, ((+A.) ((++:)1) t- ((t+ wr (++
Cus A, = (i - 1)
(8)
(1-a)w (t+a:) D D-1/2AD-1/22 1/22e 140 rus
Cu. +
When w E (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 (w = 1)). Given the configuration 0 = (ax, 6, 8, G) with 10 and e S 1/d, and let (T) and (T) be returned by LoCSOR defined in (8) for solving
Equ. (3). There exists real implementation of (8) such that the runtime TLocSOR is bounded by
1 a vol(ST) (1-
D (1-1))(----I))
++a
min { (VE'
vol(ST)
< TLocSOR <
In
CTT
2
atT
where vol(ST) and Fr are defined in (6) and C [I-a)]<1 with IT supp(r(())) Furthermore, vol(Sr)/77 $ 1/c and the local estimate ît (T) satisfies || D D T))lbo Sc..
140
5

---

Our new evolving bound (xvol(Sr)/(777)) mirrr-
ing the actual performance of APPR and empirically much smaller than ((1/(xx) as illustrated in Fig. 2. Our lower bounds are quite effective when e rel- atively large, while our upper bound is better than Anderson's when e is small. When < ((1/m) this new bound is superior to both (()((ae) and (L)((aa)) This superiority is evident when com- pared to algorithms like ISTA or FISTA [5] to min- imize the f-rregularization of for obtaining ap- proximate solution of Equ. (3). Additionally, when
Comparison of runtime between APPR and LOCSOR (left) and runtime bounds (right) as a function of e. We used the same setting as in Fig. 1.
Figure 2:
(1, and recalling that Q is M-matrix, the standard analysis of SOR shows that the spectral norm of the iteration matrix must be larger than mW 1/. Hence, 0 < w < 2 if and only if global SOR converges [55]. When w optimal (the point that the spectral radius of the iteration matrix is minimized), we have the following result.
Corollary 3.4. Let = w² A 2/(1 + -1 a)2/(1 2) and S = V. 2 0 during the updates, the global version of LOCSOR has the following convergence bound
2
1 1- va
|r S (1 (1+a)Vd,
1+ Va + Ec
where e are small positive numbers with lim = 0.
Asymptotically, when = a)) then the runtime of global LOCSOR is ((m/ Va) where Ô hides log 1/c. The main difficulty of analyzing the optimal local LocSOR is that the nonnegativity and monotonicity of (() 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.
3.4 Parallelizable local gradient descent
One disadvantage of LocSOR is its limited potential for parallelization. The standard GD (+++1) x fxx (step size 1), in contrast, is easy to parallelize across the coordinates of the update.
Instead of updating r and æ synchronously per-coordinate, we propose the following
E (0,, A LocGD) ((++1)
(*) Trs. ((+ t- (() Qr (9)
(t)
(9)
Every coordinate in S, is updated in parallel at iteration t. Interestingly, LocGD exhibits nonnegativ- ity and monotonicity properties, and its runtime complexity is similar to that of LocSOR, as stated in the following theorem (To remind, A, = 0 for Tr of LocGD in Equ. (6)
Theorem 3.5 (Runtime bound of LOcGD). Given the configuration 0 = (a, (,,,, with E (0, and < l/d, and let (T)) and (T) be returned by LocGD defined in (9) solving Equ. (4). There exists a real implementation of (9) such that the runtime TLoCGD is bounded by
1 vol(ST)
(1- (1
12222
1+a 2
1 ax
vol(ST)
}.
2 oTr
1
< TLOGGD S
min
o²T
where C = + ++aa/((1 (1-aa)Izi I), Ir supp(r(()), Furthermore, vol(Sr) & 1/e and the estimate मे = (1)22(T) satisfies | (f T 1l S.
Note that 7 of LocGD empirically smaller than that of LOCSOR. Hence, LocGD is empirically slower than LocSOR by only a small constant factor (c.g.. twice as slow), 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.
4 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).

---

4.1 Local Chebyshev method
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/ va) to meeting the stopping condition where O hides log 1/c (we presented 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 xx()) Our key observation is that if a substantial reduction in the magnitudes of r(" is required within a subset St, then the corresponding momentum coordinates are likely possess significant acceleration energy. Intuitively, viable strategy involves localizing both the residual and momentum vectors. For 2 1, denote the "momentum" " vector as (t) = (t (t--1) and 8c1+1 = 6,68+1, we propose the localized Chebyshev as the following
Do (St, (11), (()) A = LOCCH):
a(1) (1+ 1) (1) A (t)
5,1 6c41 = t2 6t) 6,)=
(10)
1-a
(1+1)) -- (1) =(1,, (((+1) p(() x(t)
1-a W i)
Ita
where L 1 with the initials (0)) (1)) (0), 8n 0,8 (1 a) a), and W
D D-1/2 AD 1/2 is normalized adjacency matrix. Our key strategy for analyzing (10) is to rewrite the updates of (t) as nonhomogeneous second-order difference equation (see details in Lemma C.8)
1/2
((+2) - 28 Wr + 1
(1 (1+8)++11 r=j+1
II 8p Qr
S S..
(II)
f=0
where we denote S3.1 = Sj n..... n n nS given 2 0 where S, V.S.. In the rest, we
define à = (1 - /a)/(1 (1 + (a) and recall the cigendeconposition of 1011/2AD 1/2 VAVT. Based on the above Equ. (11), we have the following key lemma.
Lemma 4.1. Given 2 1, x0)
x( -()), The residual r(()) of LoCCH defined in (10) can
be expressed as the following
V Z.V (0))) + d1.etuot +
>
Set1: k))KR
k=1
where Z, is a diagonal matrix such that Zell2 S and
Σ 625
j.t
---
1- A
Ita
V
Sa,
(0)
ifk=
tgg
Serra Σ j
It I A
1
(T-))
I-a
Sky
where Hy is a diagonal matrix such that Ha,lle S tt k.
This key lemma essentially captures the process of residual reduction (()) of LocCH. Specif- ically, given current iteration t, define the running residual reduction rate for (k) with 0, 0,1,2,.... t-1 of step t as 32,e. that is,
3x,1 A |
Bg a max Br.e-
(12)
|r(()] 112
Note that
2/|7 (A)
3x 1l2
45 1āi-k(1 (t 5) || | ll2
Se, < (1 a) (A) 11
(1 - a) - k)) [|((k)) 12
=0(e)
(TTaall-ar
where whether the last term can be even smaller depends on ||| (k)
Sz llz for Sk.j = =S --
However, we notice that the running geometric mean B, A (T= + 8/)) is even smaller in practice. Based these observations and the assumption 3,, we establish the following theorem.
7

---

Theorem 4.2 (Runtime bound of LocCH), Given the configuration 0 = (a, 9) with E (0, 1) and e S 1/d, and let r!T) and (T)) be returned by LocCH defined in (10) for solving Equ. (3). For t 2 1, the residual magnitude || (2l la has the following convergence bound
|r(0 2 S 81:1 II (1+33)))
where Yt is sequence of positive numbers solving M++1 2ye + 1t-1 + 3 B.))
with yo = 111 || Suppose the geometric mean 3, - (T + B,)) of 3e be such that
= cva vvح where c (0. There exists a real implementation of (10) such that the runtime
I-va
Locch is bounded by
(1+ ( va)vol(ST) )
(a(2-c) In 2yT
Toccu Se
Golub & Overton [18] considered the approximate Chebyshev method by assuming that the inexact resid- ual is sufficiently smaller than <||>(*)||, where c must be small enough to ensure convergence. How- ever, this assumption overly stringent for our case. The novelty of our analysis lies in more elegant treatment parameterized second-order difference equation, allowing to circumvent this assumption. The nested APGD(E), namely ASPR proposed in Mubboo et al. [37] has runtime complexity vol"))vo //a + |S*] vol (8*)) where S
ASP% (1w00)
-35
20000 3003) 000XS Tx (
Figure 3: Comparison of runtime between APPR and ASPR. The setting the same as in Fig. 1. Left e = 10-4 while J for right.
the optimal support of arg min {(()) + eall ca||_\ and vol(S*) = nnz (Qs-.s-) Although it difficult compare our bound to this, one limitation of ASPR is that it assumes call APGD(E)
((S*)) times to finish in the worst case. However, our iteration complexity is ((1/(V/(( c))). Asymptotically, c aa (E + 0), our complexity ((1/ could be better than I((*\///) Fig. 3 presents preliminary study on ASPR, indicating that it requires more operations than APPR.
We conclude analysis by presenting, similar result for the local Heavy Ball (HB). Note the HB method is the one when 5,0,41 - d² where ā (1 (a)/(1 + va). Hence, it has similar convergence analyses as LOCCH shown in Theorem D.8. The LoCHB has the following updates
(S a, LoCHBB 1(1) (1 + 1 xr) (1)
à²
A
(13)
S.
(1++1) (1) 410, (1+1)
(())
200) I-a ( T+a
5 Generalization and Open Problems
Our framework can be applied to various local methods for large-scale linear systems. Extensions of this framework other linear systems detailed in Tab. 2 of Appendix E. More broadly, we consider the fcasibility of local methods for solving Qx = b, where b is a sparse vector () supp(bb)] < n) and Q is a positive definite, graph-induced matrix with bounded eigenvalues. This leads to question whether all standard iterative methods can be effectively localized, raising two key questions
1. Given graph-induced matrix Q and its spectral radius p(Q) 1, standard solver A, and the corresponding local evolving process (S, (1),, ((t) Loc.), does localized version of A (over S2) converge and have local runtime bounds?
2. Based on current analysis, Theorem 4.2 relies on the geometric mean of residual reduction on [|=(( being small. How feasible is acceleration within locality constraints? Specifically, stronger bound could be established for solving Equ. (3) via LOCHB and LocCH, with graph-independent bound of
( (oo(ST) VOTT
In ) where C' graph-indeeendent constant.
TooA =
8

---

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?"
Baselines. 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 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 a = 0.1, randomly select 50 nodes from each graph to serve as e, 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 the standard and local solver. The range of c is [2(1+ajd, 100411// +/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 x 10 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.
ogt4
ogen-anly
ogbn-mag
cobn producrs
oginn proteim
30
JUD
30
100
Figure 4: The speedup of local solvers as a function of c. The vertical line is c = 1/n.
Comparison with local baselines and CGM. We next compare our three accelerated methods with four baselines. Fig. 5 presents the l f-eestimation ror in terms of the number of operations (quanti- fied as t vol(S,)) executed. It is evident that our three solvers use significantly fewer operations com- pared to CGM and the other three local methods. Again, due maintaining nondecreasing set of ac- tive nodes, ISTA and FISTA require more operations
oglaxptodrets
Inlgtt
0.
65
L8
L.
20
Operitiuss
Openatians
Figure 5: Estimation error function of operations required. (e = 10 1001/n/ than the locally evolving set process. Ours are more efficient than APPR, where (0) e, used.
Efficiency in terms of a and huge-graph tests. We demonstrate the performance of local solvers terms of different a ranging from 0.005 to 0.25. Interest- ingly, in Fig. 13, LocGD show faster convergence when a is small; this may be because of the advan- tages of monotonicity properties, which is not present in the accelerated methods. However, in other regions of at, accelerated methods are faster. We also tested local solvers on two large-scale graphs where pa- pers100M 100M 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
Ps
SeciOo
tn Bait:
Openstions
Operstoon
Figure 6: Performance on large-scale graphs.
Additiooal experimental results, setups, and algorithm parameters provided in Appendix F.

---

Table 1: Operations/runtime comparison on local clustering.
Case study on local clustering. Following the experimental setup
Operations
Run Time (Seconds)
in Fountoulakis et al. [13]. we consider the task of local cluster- ing on 15 graphs. As partially demonstrated in Tab. I, compared with APPR and FISTA, LOCSOR uses the least operations and the fastest, demonstrating the advan- tages of our proposed local solvers.
APPR
LocSOR
FISTA
APPR
LocSOR
FISTA
91 6.9e+05
6.5e+04
5.7e+05
0.127
0.043
0.093
Ø2 6.7e+05
8.9c+04
4.4e+05
0.362
0.125
0.308
G3 4.3e+05
3.5e+04
2.9e+05
0.069
0.014
0.042
94 5.7e+05
7.6e+04
4.4e+05
0.357
0.175
0.229
95 5.4c+05
9.0e+04
5.0c+05
0.072
0.055
0.084
7 Limitations and Conclusion
Our proposed algorithms may have the following limitations: When a is small, the acceleration effect partially disappears, as observed in Fig. 13. This may be due to the limitations of global coun- terparts, 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.
locally evolving set process framework 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 developing faster local methods. It 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.
We propose using
10

---

Acknowledgments and Disclosure of Funding
The authors would like thank the anonymous reviewers for their helpful comments. The work of Baojian Zhou is sponsored by Shanghai Pujiang Program (No. 22PJ1401 and the National Natural Science Foundation of China (No. KRH2305047). The work of Deqing Yang is supported by Chinese NSF Major Research Plan No.92270121. The computations in this research were performed using the CFFF platform Fudan University.
References
[1] Alon, N., Rubinfeld, R., Vardi, S., and Xie, N. Space-efficient local computation algorithms. In Proceedings of the twenty-rhird amnual ACM-SIAM symposium on Discrete Algorithms (SODA), pp. 1132- 1139. SIAM, 2012.
[2] Andersen, R., Chung, F., and Lang, K. Local graph partitioning using PageRank vectors. In 47th Annual IEEE Symposium Foundations of Computer Science (FOCS), 2006.
[3] Andersen, R., Gharan, S. O., Peres, Y., and Trevisan, L. Almost optimal local graph clustering using evolving sets. Journal of the ACM (JACM), 63(2):1-31, 2016.
[4] Anikin, A., Gasnikov, A., Gornov, A., Kamzolov, D., Maximov, Y., and Nesterov, Y. Efficient numerical methods solve sparse linear equations with application to PageRank. Optimization Methods and Software, pp. 1-29, 2020.
[5] Beck, A. and Teboulle, M. A fast iterative brrikkag--hhrehoodii algorithm for linear inverse problems. SIAM Imaging Sci., 2:183-202, 2009. 2009.
[6] Berkhin, P. Bookmark-coloring algorithm for personalized pagerank computing. Internet Mathematics, 3(1):41-62, 2006.
[7] Bojchevski, A., Klicpera, J., Perozzi, B., Kapoor, A., Blais, M., Rózemberezki, B., Lukasik, M., and Günnemann, S. Scaling graph neural networks with approximate pagerank. In Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining (KDD), 2020.
[8] Boyd, S., Diaconis, P., and Xiao, L. Fastest mixing markov chain graph. SIAM review, 46 (4):667-689, 2004.
[9] Chen, Z., Guo, X., Zhou, B., Yang, D., and Skiena, S. Accelerating personalized PageRank vector computation. In Proceedings of the 29th ACM SIGKDD Conference on Knowledge Discovery and Data Mining (KDD), pp. 262-273, 2023.
[10] Chung, F. The heat kernel as the pagerank of a graph. Proceedings of the National Academy of Sciences (PNAS), (04(50):17335-19744) 2007.
[11] Aspremont, A., Scicur, D., Taylor, A., et al. Acceleration methods. Foundations and Trends in Optimization, 5(1-2):1-225, 2021.
[12] Fountoulakis, K. and Yang. S. Open problem: Running time complexity of accelerated 41-
regularized PageRank. In Conference on Learning Theory (COLT), 2022.
[13] Fountoulakis, K., Roosta-Khorasani F., Shun, J., Cheng, X., and Mahoney, M. W. Variational perspective on local graph clustering. Mathematical Programming, 174(1):553-573, 2019.
[14] Fountoulakis, K., Wang, D., and Yang, S. p-norm flow diffusion for local graph clustering. In ICML, 2020.
[15] Gasteiger, J.. Weißenbergerr S., and Günnemann, Diffusson improves graph learning. In Advances in neural information processing systems (NeurIPS), 2019.
[16] Gleich, D. and Mahoney, M. Anti-differentiating approximation algorithms: A case study with min-cuts, spectral, and flow. International Conference on Machine Learning, pp. 1018-1025. PMLR, 2014.

---

[17] Gleich, D. F. PageRank beyond the web. siam REVIEW, 57(3):321-300 -3333,,133,,, 2015.
[18] Golub, G. H. and Overton, M. L. The convergence of inexact chebyshev and richardson iterative methods for solving linear systems. Numerische Mathematik, 53(5):571-33 1988.
[19] Golub, G. H. and Van Loan, C. F. Matrix computations (4th Edition). JHU press, 2013.
[20] Guo, X., Zhou, B., and Skiena, S. Subset node representation learning over large dynamic graphs. In Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining, pp. 516 516-526, 526, 2021.
[21] Guo, X., Zhou, B., and Skiena, S. Subset node anomaly tracking over large dynamic graphs. In Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, pp. 475,485, 2022.
[22] Hu, C. Local graph clustering using 11-regularized PageRank algorithms. Master's thesis, University of Waterloo, 2020.
[23] Hu, W., Fey, M.. Zitnik, M., Dong, Y., Ren. H., Liu, B., Catasta, M., and Leskovec, J. Open graph benchmark: Datasets for machine learning on graphs. Advances neural information processing systems, 33:221 18-2213, 2020.
[24] Johnson, R. and Zhang, T. Graph-based semi-supert supervise learning and spectral kernel design. IEEE Transactions on Information Theory, 54(1):275-285 2008
[25] Kapralov, M., Lattanzi, S., Nouri, N., and Tardos, J. Efficient and local parallel random walks. Advances in Neural Information Processing Systems (NeurIPS), 2021.
[26] Klamkin, M. and Newman, D. J. Extensions of the weierstrass product inequalities. Mathematics Magazine, 137-)113 1970.
[27] Klicpera, J., Bojchevski, A., and Günnemann, S. Predict then propagate: Graph neural networks meet personalized pagerank. In International Conference on Learning Representations (ICLR), 2019.
[28] Kloster, K. and Gleich, D. F. A nearly-sublinear method for approximating column of the matrix exponential for matrices from large, sparse networks. In Algorithms and Models for the Web Graph: 10th International Workshop. WAW 2013, Cambridge, MA, USA, December 14-15, 2013, Proceedings 10, pp. 68-79. Springer, 2013.
[29] Kloster, K. and Gleich, D. F. Heat kernel based community detection. In Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining (KDD), 2014.
[30] Kloumann, I. M. and Kleinberg, J. M. Community membership identification from small seed sets. In Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining (KDD), 2014.
[31] Koutis, I., Miller, G. L., and Peng, R. A nearly-m log time solver for sdd linear systems. In 2011 IEEE 52nd Annual Symposium on Foundations of Computer Science (FOCS), pp. 590-598. IEEF, 2011.
[32] Lącki, J., Mitrovićć S., Onak, K., and Sankowski, P. Walking randomly, massively, and efficiently. In Proceedings of the 52nd Annual ACM SIGACT Symposium on Theory of Computing (STOC), pp. 364-377, 2020.
[33] Lam, S. K., Pitrou, A., and Seibert, S. Numba: A llym-based python jit compiler: In Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC, pp. 1-6, 2015.
[34] Leskovec, J., Lang, K. J., Dasgupta, A., and Mahoney, M. W. Community structure in large net- works: Natural cluster sizes and the absence of large well-defined clusters. Internet Mathematics, 6(1):29-123, 2009.
[35] Leventhal, D. and Lewis, A. S. Randomized methods for linear constraints: convergence rates and conditioning. Mathematics of Operations Research, 35(3):641-654, 2010.
12

---

[36] Mahoney, M. W., Orecchia, L., and Vishnoi, N. K. A local spectral method for graphs: With applications improving graph partitions and exploring data graphs locally. The Journal of Machine Learning Research (JMLR), 13(1):2339-233) 2012.
[37] Moatinez-Rubic D., Wirth, E., and Pokutta, S. Accelerated and sparse algorithms for ap- proximate personalized PageRank and beyond. In Proceedings of Thirty Sixth Conference on Learning Theory (COLT), volume 195 of Proceedings of Machine Learning Research, pp. 2852-2876. PMLR, 2023.
[38] Morris, B. and Peres, Y. Evolving sets and mixing. In Proceedings of the Thirty- Annual ACM Symposium on Theory of Computing (STOC), pp. 279-286, New York, NY, USA, 2003. Association for Computing Machinery.
[39] Nutini, J., Schmidt, M., Laradji, I., Friedlander, M., and Koepke, H. Coordinate descent converges faster with the gauss-southwell rule than random selection. In ICML, pp. 1632-1641. PMLR, 2015.
[40] Polyak, T. Introduction to optimization. New York, Optimization Software, 1987.
[41] Rakhlin, A. and Sridharan, K. Efficient online multiclass prediction on graphs via surrogate losses. In Artificial Intelligence and Statistics, pp. 1403-1411, PMLR, 2017.
[42] Rubinfeld, R. and Shapira, A. Sublinear time algorithms. SIAM Journal on Discrete Mathemat- ics, 25(4):1522-1588, 2011.
[43] Saad, Iterative methods for sparse linear systems. SIAM, 2003.
[44] Schaub, M. T., Benson, R., Horn, P., Lippner, G., and Jadbabaie, A. Random walks on simplicial complexes and the normalized hodge 1-laplacian. SIAM Review, 62(2):353-391, 2020.
[45] Spielman, D. Algorithms, graph theory, and linear equations in laplacian matrices. In Proceedings of the International Congress of Mathematicians 2010 (ICM 2010) (In 4 Volumes) Vol. I: Plenary Lectures and Ceremonies Vols. II-IV: Invited Lectures, pp. 2698-2722. World Scientific, 2010.
[46] Spielman, D. A. and Teng. S.-H. Nearly linear time algorithms for preconditioning and solving symmetric, diagonally dominant linear systems. SIAM Journal on Matrix Analysis and Applications, 35(3):835-8855 2014.
[47] Stević, S. Bounded solutions to nonhomogeneous linear second-order difference equations. Symmetry, 9(1))227, 2017.
[48] Teng, S.-H. et al. Scalable algorithms for data and network analysis. Foundations and Trends in Theoretical Computer Science, 12(1-2):1-27,, 2016.
[49] Tseng, P. and Yun, S. A coordinate gradient descent method for nonsmooth separable minimiza- tion. Mathematical Programming, 117:387-423, 2009.
[50] Tu, S., Venkataraman, S., Wilson, A. C.,, Gittens, A., Jordan, M. I., and Recht, B. Breaking locality accelerates block gauss-seidel. In ICMI, pp. 1482-3491. PMLR, 2017.
[51] Vishnoi, N. K. et al. Lx= b. Foundations and Trends in Theoretical Computer Science, 8(1-2): 1-141, 2013.
[52] Wang, J..K., Lin, C.-H., and Abernethy, D. A modular analysis of provable acceleration via Polyak's 's momentum: Training wide ReLU network and a deep linear network. In ICML, pp. 10816-10827. PMLR, 2021.
[53] Wu, F., Souza, A., Zhang, T., Fifty, C., Yu, T., and Weinberger, K. Simplifying graph convolu- tional networks. In ICML, pp. 6861-6871, PMLR, 2019.
[54] Yin, H., Benson, A. R., Leskovec, J., and Gleich, D. F. Local higher-order graph clustering. In Proceedings of the 23rd ACM SIGKDD international conference on knowledge discovery and data mining (KDD), 2017.
13

---

[55] Young, D. M. Iterative solution of large linear systems. Elsevier, 2014.
[56] Zeng, H., Zhang, M., Xia, Y., Srivastava, A., Malevich, A., Kannan, R., Prasanna, V., Jin, L.,. and Chen, R. Decoupling the depth and scope graph neural networks. Advances in Neural Information Processing Systems, 34: 1965--9667,, 2021.
[57] Zhou, B., Sun, Y., and Harikandeh, R. B. Fast online node labeling for very large graphs. In ICML, pp. 42658 42697. PMLR, 2023.
14

---

Appendix supplemental material
Notations and Proof of Lemma 2.1
16
A.1 List of Notations
16
A.2 Proof of Lemma 2.1
16
Local Iterative Methods via Evolving Set Process
18
B.1 Local Variant of GS-SOR and Proof of Lemma 3.2
18
B.2 LocSOR and Proof of Theorem 3.3
20
B.3 Optimal GS-SOR and Proof of Corollary 3.4
22
B.4 LocGD and Proof of Theorem 3.5
23
C Local Chebyshev Method LoCCH
25
C.1 Nonhomogeneous of Second-order Difference Equation
25
C.2 Properties on Ratio of Chebyshev Polynomials
28
C,3 Standard Chebyshev (CH) Method and Proof of Theorem C.6
29
C,4 Residual Updates LocCH and Proof of Lemma 4.1
30
C.5 Convergence of LOcCH and Proof of Theorem 4.2
34
C.6
Implementation of LoCCH
37
Local Heavy-Ball Method LoCHB
37
D.I Standard HB and Proof Theorem D.2
37
D.2 Residual Updates of LOCHB and Proof of Theorem D.6
................................
40
D.3 Convergence of LOCHB of Proof of Theorem D.8
43
D.4 Implementation of LoCHB
45
E Instances of Sparse Linear Systems
45
E. Table of Popular Graph-induced Linear Systems
45
F Experimental Details and Missing Results
46
F.I
Datasets and Preprocessing
46
F22 Problems Settings and Baseline Methods
46
F3 Full results of Fig. 15456
47
F.4 Results local clustering
52
Related work
53
15

---

A Notations and Proof of Lemma 2.1
A.1 List of Notations
In the rest of the appendix, we use the following notations:
Description
An undirected connected simple graph with unit weights.
A
The adjacency matrix of g.
D
The diagonal degree matrix of g.
W
The normalized Laplacian matrix W D- AD-1/2,
VAV
The eigendecompossitinn of W is given by W =D AD 1/2 VAVT, where
= A2. An), with = A 2 X2 2 ... 2 An 2 -1 When Xn -1,9 g
bipartite graph.
Q
The underlying matrix of Equ. (3) I 1-a D 1/2 AD 1/2 Ita
(0
Given estimate (!)), the residual is defined as (() A b Qx!!),
F(0)
The Di/2-shifted residual is defined as f(t) (1/2
The damping factor a, which lies in the interval (0,1).
&
A pre-defined constant á = (1 -1- / Va).
Ta
The total runtime of a local algorithm A.
T((x)
For t 2 1, T. denotes the Chebyshev polynomial of the first kind, defined as 1xt 2xTT (2) T (x), with (x 1, and Ti (x:) w.
di
The ratio of T,-1 and Tú, i.e., dc Ti-1 T (1 (T(((tt) and 81+1 (2]+a - -5,)=1 1-o
with 61
T4o
61.t
The product of all 8, i... 81 - I 8, By default, we set d0:1 = O.
S
The intersection of all t-ordered sets, i.e., Sr. = S1 O S2 ... Sr. By default, we set S = V.
St
The complement of S,, i.e., S, V.,s,
SJ, = S nS
n3, Spa-1 n Sr By default, we set St, = S,.
A.2 Proof of Lemma 2.1
Lemma 2.1 (Runtime bound of APPR [2]). Given (0,1) and the precision S. L/ds for node E V with p + 0, + e. at the initial APPR(G, a, defined in (2) returns an estimate of æ. There exists a real implementation of (2) (e.g., Algo. 2) such that the runtime TAPPR satisfies
TAPPR θ (1/(EE))
Furthermore, the estimate ît p satisfies D- (f--r))lle c cnd vol(supp(ft)) 2/((11 aa))..
Proof. To find upper bound of TAPPR, add time index for all active nodes u1, U2,.... 04 processed in APPR of Algo. 2, a real implementation of (2). The parameter is the number of active nodes processed, and uy the node dequeued time So, updates of p and z are from ez) (p) (0)) to (p(t), (()) follows:
(p
(1) (())
For each active the updates of (p, x) by definition can be represented as
(1 +
a)zu,
(i-1)
(5-1)
(1 - a (i-1)
p = =p (i-) + (i-i)) eu. (1)) = z(( (i-1)
ev, +
AD*
eu:
2
2
16

---

Since (0)) = e, 2 O, then z 2 0 and p 2 0 for all by induction. Note that ||A -eeee
Cu, |l
=
we have the following relation from the updates of
12
|| 11
||
(i-1)
(2-1)
Note edu, S Zz:
(1-1) for each active u4- Summing the above equation over i 1. ..,,,
,t. we have
|| 1) - 11)
|| ||
||x00//2
| 11 1
edu
(1-1) uu
= Σ
a'
1=1
where note ||
||e/l1 1 by the initial condition. Since Ei= du, exactly captures the number of operations needed, the runtime Algo. is then bounded
|1
Σ dus
TAPPR = +
aE
To check the quality of estimate p. using the updates of (()) and summing over all i, we have
eu. =a ( (1++))
(1-o) a
(())
a
(1-1)
AD
2 (i-1) 2(1) = π II,(1),
2
11
where is the PPR matrix. The above gives us -p(t) = II (() Since G is undirected, the II matrix satisfies Te (d,/ /d, ) [0] where Tu[u| is the u-th element PPR vector of sourcing node v. Consider each u-th element of II 2(1)
(1) du Tulu]
III 2() = Σ
(c)
To [u]
1 edu Σ Tu [u] = edus
de
VEV
NEV
vEV
where the last equality is due to Evey Tu 1. Hence, ( ( (II ())) S ed,, which indicates [p-)(p(()-TTT)
I e. To see the bound of vol(supp(p note for any u (t))
it was an active node and there was at least (1 (1-1a)) a)/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
(0) Σ Euesupp(p) 2u
Σ du kEsupp(p)
Σ
Σ
2. a2/2
2
<11 a)/2
c(1 - a2/2
(1-a)e
[Esupp(p)
@Esupp(p)
Algo. 2 APPR(a, e, 9) via FIFO Queue
1: Initialize: p < 0, f- C., Q + {*, s,,,,, s)
Indeed, Lemma 2.1 is special case of The- orem 1 in [2]. The proof outlined above ad- heres to the key strategy demonstrated in that theorem, which involves exploring the mono- tonicity and nonnegativity of z. The real implementation of APPR, shown in Algo. 2, presents typical queue-based method. It has monotonie properties dur- ing the updates of and z (Lines 10-16 of Algo. 2). It also holds clement-wise that p Z and z > O. The operations of Q.nnqueee(u), u), Q.dequeue), and Q are all in (1). Line to Line 9 is to design the marker for distinguishing between S, and S.+1- all active nodes are processed and more active nodes added into Q then Q will be empty, and finally, the algorithm returns an estimate p of π.
2: while true do
3: u < Q.dequeue()
4: if == then
5:
if i Q then break
6:
7:
1 1 ++1 Q.enqueue(*) continue
// Starting time of S, // Marker for next S
8:
9:
3 2u
11:
Pu + Pa + a
zu x (1-1a)/2
13:
for E ENu do
2v + 2v + (1-a)
d.
15:
2, 2 ed, and $ Q then
16:
Qenqueue(v)
17:
if z. 2 ed,, and u $ Q then @enqueue(u)
18:
19:
return
17

---

B Local Iterative Methods via Evolving Set Process
Justification of an equivalent condition. We make the justification of equivalent stop condition for solving (3). Note we require a local solver to return an estimate म satisfies
|
( π) b0 I
(14)
Since we define Qx as (I- and the stop condition
1-a D Ita
1/DAD 1/2) 1/2) x
D
e-1/ With (()) b- (x((())
I+a
2cre
[ 1/-1/2 1 1 1+ a
ensures the estimate # = (1) satisfies (14). To see this, since π = D1/22", we have
D 1 (# ④ 1oo
| (x x² ll)
| 1/2 ( Q (t) Q-1b))lcc
|| 1-1'/QQ' D 1/2 12,111))1m
| 1DD11/2-11
D-1122((0) oo
20re
D
1 a
where D D-1/2Q-1D D1/2
1-a D D-1A) 1--1 Dizo
D- A)' This leads
Ita
Ita
2ne
| (+ - T) 1l < ||
L- D
A)
T+o
I+a
2=
1 1+
2
2a
1 a
=6
B.1 Local Variant GS-SOR and Proof of Lemma 3.2
The Gauss-Seide Successive Ever-Relaxation (GS-SOR) solver (see Section 11.2.7 of Golub & Van Loan [19]) for the linear system MT via the following forward substitution
for i V = f1,2...., n) do
" Σ s=/+1
(1+1) P.
w
MyP;
(1+1)
Mypp
/Ma + (1 (11-w)) w)p
where p is updated from (t) to p (t+1), When the relaxation parameter w 1. GS-SOR reduces to the standard GS method. Equivalently, let A, = (i - 1)/n for i = 1,2,...,7. then GS-SOR updates can be sequentially represented as
for in = f1. do:
pittAm+) (HtA)) +
W
8; Σ
i-1
Mapp
((+.))
((+i,)
Σ
a
Map (ist (+2+)) e,
Therefore, is natural to define the following local variant of GS-SOR.
Definition B.1 (Local variant of GS-SOR). Consider the linear system MT For t 0, we given an active node set S; = {u1, U2 ..., U|s,I } and let A; = (i 1)/(S, for i = 1,2..... SS), providing w e (0, 2), is natural to define the local variant of GS-SOR as follows: for nn in S = {uI, u2, ) do:
((+tbi+)) + ((+ti)) +
w Mu,a,
Su,
Munu
plt (++
Σ
Mun Mu,u, plet
eu. (15)
j=l
where S, C V. When w = and S === 2,, ...., n). it reduces to the standard GS.
18

---

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 = - (I 1-a and $- = l.s. The linear system MT- = 8 is equivalent to Equ. (1). Given pl)) = 0, 20) = e, with
2 (I+ AD A))
(0, the local variant of GS-SOR (or M can be formulated as
(t+A:)
(t+A,) Meus
p p(tAmt)) < -pltti))
Wwy.
Buc (1++ c- (1+D.)) Whys
Maros
Muca.
where ug is an active node in St satisfying Zu, 2 cdu, and A; = (i - 1)/[Sa]. Furthermore, when w Lta, this method reduces to APPR given in (7), and there exists a real implementation (Aglo. 2)
of APPR such that the runtime TAPPR is bounded by
vol(Sr) abyr
In Cr where
vol(Sr) frr
- 1 E' (I II
2
TAPER <
CT =
fr
T
Proof. Note that MT 8 is equivalent to Equ. (I). We first rewrite p 11)+)) in terms of the residual z. The residual at time + A; can be written as 2((+Ai)= = MM (p(++A)) Note Su, Pi; Pu, (++A.) -j= Mu.u,P Ph; Mu. (++A.))
= (s Mp((+A:))a )a: = 2n; (1+A.) Then, the updates of local GS-SOR defined in (15) can be rewritten as p p(++Am+)) = (((+ii) + Musus Sue
Ej= Mu,,,
((+A,)
Bu.. Hence,
the updates of x ((++A)) can be written
(++A,)
(++A,) Mew;
2 (t+A++1)) = -8- (p 8 M (+tta)))
Whise
(((+A)))
Whys
Mapu
-uus
Musu
Note the diagonal element Mars. = (1 + (1+a)/(2a). Hence, when w ((++A++)) ((++Ai) + ahuu
2 2
1aa
we have
(I+Ai)
Eu:
(1+A4))) =2( (++))) 0:
((+A)) Me., =((+A
(1-- 1 1+a I 1-a 2
AD AD- )eas
2
The above updates match APPR's evolving set process formulation in (7). The rest is to show a new runtime bound. Adding fn-ormm on both sides of the above equation, then note [| =
||((t+)))/11 az2u; (1+A.) for i = 1,2,..., [.e... We have
.
Σ
|Sc ((+A,)
||2 (1+1)111
Ei=1
Zu:
11= 112
II (1 aB)) ||= ==0
1210
where we define Be =
Lett= 11,,, 1, we have
||
In
1
T-1
Σ
In aBB)) < Σ 03
T-1
T
1 oryr
In
| 11 [z(
[z
|1
where the first incquality is due to In(1 + x) S x for x > -1. For each nonzero node u E In [1T)) = * 0, V}. consider the last time l' that it was altered. Then, either the alteration came from u being an active node, with (!') 2 duc, and after the PUSH operation it became (()) (1 (1-a)d= e; or the alteration came from a neighboring node v. E (N) pushing its mass onto
(#u
zu
u. which ensures that 2((") (1-a) (t') (1-a) These two cases provide lower bound of 1-a
2d 20
2
Hence, || 2 <(1-a) Irl which leads to the corresponding constant CT.
2
To see the lower bound of I/e, note edu, I Zus
(++Ai) for all i 1,2,.... ... |S,].
Then we have
IS.
vol((S)) <
(++A))
i=1
= 8.11= (111
< Br,
19

---

where we defined B, ;- El +
lzu,
and the last inequality is due to the monotonic decreasing of
[1=(1)]
12, 1 > ||=(0||| 2 ... 2 || ||.. Applying the above inequality for all t 0, 1, 2....... T
1, leads to
(vol(Si) S B
T-1
vot -
vol(ST)
-
where the last derivation is from the fact that fr
1 T
ΣT- 1== B
== Σ |2
Remark B.2. The connection between APPR and the Gauss-Seide is not new [28, 29, 16, 9]. Our work is the first work that has linked APPR and the Gtaus--Seide with locally evolving set process.
B.2 LOCSOR and Proof of Theorem 3.3
In this subsection, recall we defined f(t) = (1)2
Lemma B.3 (Local iteration complexity of LOcSOR (W S 1)). Denote S {un, U2..... Us 3
as the active node set at the t-th iteration. When w E (0, 1), all vectors f(() 2 0 are nonnegative and magnitudes are decreasing [[=(t+)1|2 Then, at iteration T, we have
< [1. Let be the total number of iterations needed.
(1 + a) TE YT
(T))
(
T
18d (++A:)
In
|
7
F(
(16)
where T, = 1-- Y/ is the mean of active ratio factors defined in Equ. (6).
Proof. Recall E St {u1,,-.. .......:0... U[s, and A,
LOCSOR in Algo. 3 updates [8,1
i-1
x (+A441) = (++.i) + wr (++1)) Cu:
(++2i)_ eu, + (1-a)w (L+2.) (Lta)) D-1/2AD AD 1/2 (+41+1) = (1+1)) wr a;
Ita 04
Note
0 during updates when w E (0, 1] and recall zi(1) (12 we have
|| + (+(++A.) ealli = ||((++/))) + (1-a)w p(t+Aa)) AD- eu,ll 11
Ita "i
||ا+)م)) 1 + wf(++A:) = || + (1-a)e ((t+A.)
1+a m
Summing over the above equations over u4, we have
[S ((++.) = 1
|s (t++a,)
p(++1) = || 1
2ow
2ow
Tut
|
(17)
1+ 1=1
Given {ai)/D and x, E (0, the Weierstrass product inequality provides 1- -T x S IT=-0 (1 - 2t). By using this inequality, we continue to have a lower bound of T as the following
T-1 2ow7t < II (1 (1- 2awne I+a) 0=
T-1
(((T)
(1+ - (17l1 11000111/)
1
< T.
I+a
2awhr
0
To get upper bound of 7h note each active residual F
((LAA.) pushes at most
(1-a)w
((+a) times magnitude to
Pultr Fusp and Pujs. hence, Ej=r Tiy
([t+A)))
will increase by most fuu
(.(++A))
(1-a)w ((++,,)
(1+o)
in total. Hence, overall u, we have
|5
|S
((+A))
11
S,
i=1
20

---

We reach the following lower and upper bounds of %/.
2 11 |
S 7t
+((+a,)
Fuc
21
To
1|1
|
| 1|1
check the upper bound of T, from Equ. (17), |[
11
II:=0
(1
20
|| 1 and
I+a
|
In 1 = Σ
T-1
Yo
T-1 Yown
T Ta
(1+a)
|
||001
1 +
1=0 l+a
|
t=0
where the first inequality is due to In(1 + S a for x -.
Theorem 3.3 (Runtime bound of LoCSOR (w 1)). Given the configuration 0 (ax, c, 8, 9) with 10 and S 1/d, and let r(T) and (T) be returned by LocSOR defined in (8) for solving
Equ. (3). There exists real implementation of (8) such that the runtime TLoSOR is bounded by
1 ++a vol(Sr) (1 - [D1/2 S TToCOOR I
1+a
1 vol(ST) In C )
min
OTr where vol(Sr) and "r are defined (6) and C (1-a) Trl with IT supp(r(T)), Furthermore,
2
O6
1+a
vol(Sr)/7rr S 1/e and the local estimate == 11/2 satisfies || 1 1)) S..
Proof. After the last iteration T, for each nonzero residual ((T) * 0. E Ir, there must be at least update that happened node Node has neighbor uu E (u),, which was active. This
(I=a)=(t)
neighbor pushed some residual (I+o)d to u where t' < Hence, for all E Ir, we have
(1-a) a)
(1 - aa)aaddd. (1 (1+ =[[T]]
2a(1 a)
|1
Σ u
Σ
(1 + a)do
(1 add
(1 + a)²
MEZT
MEIT
where the second inequality is because p((")
was active before the push operation. Applying the above
lower bound of [|
[a Equ. (16) of Lemma B.3 and note || ||
20x(1 (1 + a:), we obtain
| 111
||((0]
1+a =
[
<(Irl
20(1-a) (I+a)²
c(1-a) Irl
The rest is to prove upper bound 1/(ae) of TOOSORR Recall that for any active node we have residual updates from Algo. 3 as the following
(1- (1-a)wr!" a:wrr DI
D1/2 (1)2,(()) - wr(()) 12/2) eu
+ a
(1) 1I2 eu to the left and note AI DD12 1/122e.l1 e/2lll1 Move -w
= Vdu, we then obtain
|D (1) 111 (1 -lwaww B (d
|D | (e+ + + du /d
1+a a
Hence, for each active u, we have 2awn || 12, Ita
|D1/22((++1) Summing them
over all active nodes and noticing r' > 2n /da/ (1 + a) by the active condition. Note
and (02 20 we have run time bounded by
1+a)
122
D D1/ (e+1)))) <
(1+a)
LLCSORR Σ d S
2a
WE
20E
Combining the above bound and the bound T shown in Lemma B.3, we prove the lower and upper bound of 'LocSOR. To check the lower bound of 1/c,i... vol(Sr)/Tr S 1/e, note Then we have
20edu, Ita
((+A.)
< Fa;
for all i = 1,2...., |Sy].
20E
vol(S.) S
Σ
((+A.)
+ a
i=1
= = (0)
5 Sr 1/2,
20%
1 a
21

---

where the last inequality is due to the monotonic decreasing of |D [1, i...
2a
2
1+a
11)22 2 ... ||((T)|1- Applying the above inequality over all 0, .., ,
it leads to
7 e evol(S)) S
T-1 Σ
vol
-
1=0
vol(Sr)
5
PT
Algo. LocSORR(a, s, via FIFO Queue
1: Initialize: r +- ce It f- 0, c I-a t
2a
2: Q <- { {*,8)
// As we assume < 1/ds 3: while true do
The real queue-based implementation of Loc- SOR is presented in Algo. 3. It has mono- tonic and nonnegative properties during the up- dates of > 0 and æ > when w (0, 1J. Same as APPR, the operations of G.enqueeewu). Q.dequeue(), and & Q are all in O(1). During the updates, one should note that the real vector r presents (122,(() while the vector x is (1) this case, the original active node condition is implicitly shifting from |ru] 2 20cyd_ to / 2 20edz,
4: tr --= Q.dqqueee))
5: if
then
6: if O = 0 then break t < ++1
7:
8:
// Starting time of S // Marker for next S141
9:
Q.enqueue(*) continue
10:
F Tu
12:
if |ral <e.ed, then continue
I+a We use this shifted active condition Lines 11 and and inactive condition in Line 5. When E (1, 2). it is possible |r, < e- eded, and LocSOR will ignore this inactive node during the updates. This
1+a
14:
Iu Iu +
15:
Tu 4 Tu w
16:
for E N(u) do
To f To + (1-a)w (1+a) du
18:
if |ro| 2 ed and & Q then Q.nnqueee(v)
step makes sure S. {u : |reg | 2 1+a during the updates.
20cyd
20:
if Irl 2 c. ed,, and a o Q then
Q enquueee(r)
22:
return
B.3
Optimal GS-SOR and Proof of Corollary 3.4
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 Stielijes matrix and set relaxation parameter as
((B)
11
(1+
(18)
1+ V1 (BB
11 -
(B)
where p(B) is the largest eigenvalue (in magnitude) of = diag(Q) ()) (QQQ then
w -1115 1111spp p(Lw < -w--1
(19)
where L_ = (diag Q) w - (w
(w 1)diag Q)) with Q = diag Q) - Qu Q
Corollary 3.4. Let w & 2/(11 VI - (1 (1-1((11//11111 + a)²) and V, Vt 0, the global
version of LOCSOR has the following convergence bound
2
(1
1 11 Va
( 11
(1 Va +a)
(20)
(1+a)
where e, small positive numbers with lim = 0.
22

---

Proof. Recall -
D ADD 1/2 and we consider the underlying graph as simple which I+a
1-
means A has 0 diagonal. Hence, diag( Q) I and B is defined as
B
1- 1+a
D 1/2 1 ATAAD-11 1/2, (
1-a a
Since Q is Stieltjes matrix, then Lemma B.4 gives a bound on the spectral radius of Lu
VI (B) -1 1 - -
2
1/2
2(1 +a)
1/2
Vo
p(
2/a
1t va
Recall that Gelfand's formula states [52]: Given spectral radius p( maxie(n) [>s(Lw))] where x,(-) is the i-th eigenvalue, there exists a sequence { such that || (p(Lwr) +G) and lim Et 0. The standard SOR method defined as the following
x (++1 (diag( Q) - w -Q)) (w'b + (w Qu (w- - 1)diag((a))))) x!!))
w (t) w" (diag Q) w Qt) b.
Note (()) Qe QL Qee QL Q--iQ--p0 (0), We have
||r(() 112 = IQQ Q (
1|2 -
[Q||2|| 112lllQ11000 1allQ-11 (12(l1 1|2
2
| ll2 1
20
2 [ 11
S 1
(1 (1+a)Vd,
(1 + dT
To meet the stop condition, we require |r(0 This leads to find such that (1talvd,"
S
20e/d. 1
It is enough make sure |||(1) S
200
(n(0) 12 < 2| L 1l2
20E
1- va
t,, +.) ) s ac.
(1 + ad
(1+a)vd,
(1
1 Vo
When Er = v() then the runtime of global LOCSOR is Ö(m/ (a) where Ô hides log -
B.4 LoCGD and Proof of Theorem 3.5
The local gradient descent, namely LocGD is to use (f+1) = (t) tr's, and (++1) = ((t) Qrs, where S. = (u, Irth ((++)))] 2 2or Aowe/dy/t + a}} where A; = 0. Algo presents our actual
implementation of LocGD via FIFO Queue.
Algo. 4 Lo, 9) via FIFO Queue
1: Initialize: r' f- ces, æ 4- 0, C
2a
1 T+a 2: Q {s) Assume < d
3: t=
Algo. 4 presents LOCGD similar to the real queue-based implementation of LocSOR. It has monotonic and nonnegative properties during the updates of r 2 and x 2 0. Again, the operations of Q.enqueue(u), Q.dequeue(), and & Q are all in O(1). During the updates, one should note that r presents 11/2,(()) while æ is 1 All shifted conditions are the same as of LOCSOR. The key advantage of LocGD is that it highly parallelizable, while LOCSOR is truly an online update, it is hard to parallelize.
4: while 2 do
5: S. -
6: while Q + 0 do
7: u + d.dequee
8:
S. rw)) En
9:
+ xn + ra
10:
Fu 0
11:
for (u, f) E S, do
12:
for v E NN(u do
To + To + (1-a)F
13:
(1+a)d.
14:
if |r, > c- ed, and
Q then
15:
(.enqueue())
16:
1+1+1
17:
return Z,
23

---

Lemma B.5 (Iterations of LocGD). With the initial x0)) = 0, (() = b, b. supp(r(0)) (0)), denote
f(() = 11/2,((), LOcGD defined in (9) has the following properties: 1) x+) 2 0, (t) 2 0 and
||(0|| 2 | ||2: 2) The residual and estimation error satisfies
S ((+A.)
| 1 (
-11-1
111+a
207
|1
7t-
where A O.
Proof. We first show x01) 2 0 (()))) 2 O are all nonnegative vectors during the updates when
2 O. This can be seen from the induction. At the initial stage, x(()) 2 0 and r(0) = 2 0.
Now assume that for any 1 2 0, ((()) > 0 and r()) 2 0. Then (+(+1) a0 + +rs. S 2 0, and
(t)
((+1)) = =r (t)
+ 1-a D
1/2 AD 1/2
(1) S. 2 0. Therefore, x(() 2 0 and (()) 2 0 for all t. Note
S.
p((+1) ".. +
1-a
AD-1 11- rs. and since || | - ll = | ll. we will have
I+a
| (1
(1) 11 1 + a [F(0]]
2a
|| where 70 == || (1/
Then, we can bound the total residual as the following theorem.
Theorem 3.5 (Runtime bound of LocGD). Given the configuration 0 (a, 6, 9) with (0, 1) and S 1/ds and let (T) and (T) be returned by LocGD defined in for solving Equ. (4). There exists a real implementation of (9) such that the runtime TLOGGD is bounded by
1 vol(ST)
(1
S TLCCD I 1
1 are
(ol(ST)
Cl
|
min
In
2
O5T
0JT
where C= (1 + a)/((1 - [I l). Ir = supp(r(T)) Furthermore, vol(Sr) Fr S 1/c and the - 1l Se
estimate क == 1122 satisfies
Proof. We first show bound 1/(ac). We first rearrange ((++) p()) Qr.. into
11 (t+1) + DD222 T
(1) = (1),,((() + 1+a a
1 a AD- DD 12
(!) || Hence, it leads to Note (+)) 0 and |AD-1D' 1 s. |11 |2
(1)
IDd 1
s.
11
1+
(D1/2,2(0) 111 - (
2a
At each local iterative t. by the active node condition 20c// (1 + a) S r() we have
cvll
Σ cd, S
(1+
1
D1/2
D
(t) S.
2a
2a
MES,
Then the total run time of LocGD presented in Algo. is
(S, < 1 (11(1111 2a )
(0) - 11122,10011)) (1 <
1 - a
2ac
Therefore, the total run time is at most TLOGGD == ΣT t= (S, of by the Weierstrass product inequality [26] and Lemma B.5, we use the similar argument made in Lemma B.3 and continue to have
< 2nc
Lia For estimating the bounds
(1 + a) (1- |
(1 + a)
ST
In
2077
20%
|
24

---

Note that each nonzero ((T) has at least part of the magnitude from the push operation of an active node, say tu at time t' < T. This means each nonzero of D D1/2,(TT) satisfies
a 2 (1 a (t)
20cde, /11 + (1 a) a)
2x(1 -ac for u E IT.
add.
a)d..
(1 a1+
Hence, we have |
2n(1 aadd- and T is further bounded as (Ita)²
(1+a)
||
(1+a) Cr where Cr
(1 + a)
T
207T
20(1-a)e[Zrl (1+a)²
207T
(1 a) Lrl
The lower bound of 1/c, Le., (o) / I 1/c, directly follows a similar strategy of previous proof by noticing that || (1), (t))) 1|1 is monotonically decreasing.
Remark B.6. One may consider designing local methods based on Jacobi and Richardson's 's iterations.
Indeed, these two methods have the same updates as standard GD. Recall the standard GD method
to solve (4) is (t-1) (t) + p(()) (t+1) p(()) Qt()), The Richardson's is iteration x
(I - wW)x(() + wb, i.e., x (f+11
(t) w( W x
b). The optimal w' 2/ //min + Amax)
where Amin 2a/( (1 a)) and Amax S 2/(1 + a). Hence one can choose w 1 S w" [19]. It leads
x (++1) (1) + ((),, One can get the same result for the Jacobi method.
C Local Chebyshev Method LOcCH
C.1 Nonhomogeneous of Second-order Difference Equation
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 nonhomogeeeenn difference equation
Jt+1 + pxt + qIt-1 = ft. t= =
(21)
characterized by the following two cases
-2-d1 (i1 (72200 - I1 Σi=1 X
+ \/. I1 ~1æ0 + Σi= #t
p² # 49
xt-
(-2)² ao Ex-n
kSr (-D/2]0F]
+ t(
- (I)
(- () xo + Σx=
p² = 44
kml
2/
where 2i d2 are two roots of X2 + pd + q 0, and the summation follows convention Eb=1 0.
Based on the above lemma, we have the following corollary
Corollary C.2 (Second-orde) nonhomogeneous equation). Given |a| I 1, the second-ore -odder nonho-
mogeneous equation
T++1 22x + It-1 fe
has the following solution
20 + +((( - xo) + Σi= f)fx
if = 1
sin(O1) -in(0--)))) + Σi= sin(0(t-k))ff
ifla < 1 where 0 = arccos(a)
8
8in
- (=10
x1)) + (-1)² (1)= k))sk)
= -1
(22)
Proof. The first and last two cases can be directly followed from Lemma C. Since Ån and Å
are the two complex roots of X 2ad + 1 O, we write X₃ r r(cos i sin 0) and
A = re 10 = r(cos 0 - sin 0). indicates
x - 2ad + (1 (0-re-)) () - re r0- x2 r( - + x = 0.
25

---

Since > a², "1 = a a -iiii-iia a², and Re(A))² + Im())) = (1- = L, then r = 1.
Then 0 = arccos( arccos(Re(M)) Re({))) arccos(a), and sin(0) Im()1) -1- a² Finally, \f E
cos(to) + siin(t)), and
\ι cos issiin 8, 22 = e it0 e e-it =1
d2 = cos isn
A2 dι -iisin 0
1 = cos(Ot) (t) + i sin(ft), X33 = cos(Ot) - sin((
N/2 \ e e-t ite eit sin (0t).
1 1/ 1t-1 = si - 1)) 2122=1 iř A9L1
Based on these, we have
Sk
Dt 22 - it 2220 211-
Σ
+
X/
\1z0
+
Σ
24)
1m
R=
ii 112:0000 x1) + X2(000 - -1100)
Σ
+
A22 - A11
A2 A
k=1
k=
(2isin(0(t 1))xr 24 sin(()211))
1
X }
-2i sin
-2i sin 0
X1
sin sin(Ot)æ1 sin sin(O(t 1))x0
1
(1/ - \2- x11+)) fN.
sin 0
22 sin 0
sin(0t)æ1 - sin sin(O(t - 1))20
1
24 sin(((t - k))
sin 0
-2 sin 0
k=1
sin(Ot)æ1 - sin 1))
D kkm sin(t k)),f
sin 0
sin 9
Lemma C.3. For 1adl S 1,i 1,2, .n. equations y following solutions
(++1) 22,8 (1) + Vi (1-1)
0 have the
+ +(yi Vi (0)
if
A =1
Vi
(8)
sln(8,0w,)) sin(t)((--
if
|& < 1 where 0; arccos(A;)
(23)
sin(H.)
(-1) t(y/ yi (1))
(0)
if
x; -11
Furthermore, when (2) Xavi (0) for = 1, ...... then solutions can be simplified
(U) Vi
if \;- =1
yi
(0)
cos (0+t) if |>j where 0; = arccos(A.)
(24)
21,(-11)² if A -1
Proof. The first part a consequence of Corollary C.2 by letting fr = 0. To see the second identity of this lemma, note that
(A, (0,() sin(O - 1))) cos(0.t).
sin 0
26

---

Indeed, by expanding sin(0,(t 1)), we have
sin(0, 1)) sin(0,t-00) sin(0,t) (o0c)))
- 0;) + cos( (in)s sin(O:))coos(.) cos (0,1) sin(O:)
= X sin(0,t) - (0,( sin(0;)
Hence, when vft \ (0) we have
sin(0,1),,, sin(Os(t - 1))a (0)
sin(())) sin Di
cos(0,t)v (0)
sin(O;)
sin 0g
Lemma C.4. Given t 2 1, |A S 1, the second-order difference equations
V1 (1+1) 2 (0) v (t-1) ha,0.
i= 1,2, ........
have the following solutions
Vi (0)
+tu: )) + Ek=L k)hi,k
if d =1
yi
(1)
sin(O.t))) (in().(te1))1)) (0)
- sin k))
Aal < 1 where 0, = arcoos(A:) if
sin(e.)
(-1) - llv (0) yi () + Σ (-1)----111 b))e
if =
(25)
Furthermore, with initial conditions Vi (1) >eti Vi
(0)
(0)
be simplified
V/ + - (t- k:)
if
A;
V/
cOs(O;t)); (0) + - (0n (t-b)) T
if
Dal < 1 where 0; = arccos(A;)
(26)
sin 0:
(-1)'vt
(0) C-
if
X -1
Proof. The first part is consequence of Corollary C.2 by letting fo = ha, We prove the second part by considering three cases of X as
Case 1. When X; = 1, we have ((+1) above
2y/ + +gi
(t)
(1-1) = hig For 2 1, the solution the
1-1
Vi v; (0)
(1)
+
(1)
vi (0)) + Σ (t k)hi,k
Vi (0)
k))a,k
where the second equation is due to y (1)
\1 (0)
(0)
Case 2. When A = -1 (G is a bipartite graph), we have (1+1)
+ 2y/;
+ 11
(t-1) higt-
For i 1, the solution is
k)bEA (-1)++] = (-1)'y²0)
4-1 1-) (t -ki)hi..
(t) y: = (-1)'y (0) (-1)*
+
k=1
Case 3. When [Agl < and define 0, arccos(A:). We similar argument in Lemma C.3 and continue to have
V1 sin(O:)) sin(O:(1 1))9) (0)
D sin s0n(O. k))
sin 0:
sin 0:
(A,sin(Ost) sin(O,(t - 1))) Vi
km sin - )
sin 0:
sin 0,
sin(O,(t k))
cos():t) (0)
li.k-
sin 0;
27

---

C.2 Properties on Ratio of Chebyshev Polynomials
The next lemma presents the properties of Chebyshev polynomials.
Lemma C,5 (Chebyshev polynomial bound). For t 2 1, the Chebysher polynomial of the first kind defined recursively as
Ti x. T (x) 2x - (1)) with To(x) 1,
Fort 1. 2 1, define óz -T- / (T(()=)) then
1. Tix
I-a and ò, defines the following sequence
82+1 = 2 1+a a
1 1-a a 1
1 - de
where ō1
2. The flosed-form dn can be upper bounded as
2 ā' a-
1 Va
8L T((
S2
I va
3. Note 61 = To/TT = 1/2 =
I-a Ita' the sequence {ōe} satisfies 8, < 1, 14 2 1 and
20 - bedr+1, t 1,2,
4. Denote 6/23 = IT for t 2 2 and set the default value d5:511 1 for 2 then
ó 5x+12 S 26'-k fort 2 k 2 0.
Proof. For the first item, let x T-a
use the Chebyshev equation, we have
1=2 + 1-a T--1
T. T--
2
1 + a
8.41 - 850141
o+1 =
+ (2(11))) 8c.
T++1
1 1- a
1 1-a
For the second item, for all t 2 0, i
xtz-
00, is well known that the Tj can be rewritten as
2
+
+
x²
+
T {
2
2
For our problem, recall we defined ā 1-va and using
1+a I-a
f 0, one can verify that
Itva
1 +
à+ā-
1+a
à
-
ã
-ä-t
36
To x
= T₁
2
1-
2
2
and
To
T2
Ti
T2
T
1
28
1 Va
2n
=2
T
T2
T3
T,
T,
ā' + δ-t
2
+
For the third item, is sufficient to show that 8, < for all 1 2 for x recursively, since 5, = x and
1ta I-a
This can be done
1
81+1
2x 8,
2x
For the last item, note when / > 1 and k > 1, we have the following inequalities
k + -
ãh +
= 5L-A 2k +
II
5
< 2 1-
ā² + ā
2 ā a'a
ā +1
1=1
where note क
[1,2] for t > k.
28

---

C.3 Standard Chebyshev (CH) Method and Proof of Theorem C.6
This subsection introduces the standard Chebyshev algorithm. Our following theorem prove the runtime complexity of the Chebyshev polynomial iteration for solving Equ. (3).
Theorem C.6 (Standard CH). For 2 1, consider the Chebyshev polynomials to solve Equ. (3) as
x (1+1) xx + (1 8((+1)r 61141(1 - xl(-1)),, p((+1) -28++1 1 11
where x( (0) 0, x0)) (0) + ( and de+1 (2 - 51) with 61 1-a Assume e 1/dg.
then the residual has the following convergence bound
|
1 1-vaa
Va (+
2
||b|l2-
Let the estimate be 21/22 the the runtime of CH for reaching ||D 1/ S
20e with Ita
| тT) To)o S e guarantee is at most
TCh 0 s mm
[ 1 Va 2/a
2 2) = m Vo
In
Proof. Recall eigendecomposition of W VAVT where V (01, U2..... Un and each v, the eigenvector. For 1, the residual (t) can be written as n second-order difference equations as
T (++1) - 26, AV (1) + s V (t-1) 0,
where each i-th element-wise equation of the above can be written as the following
i p (t+1) 28,+1A2 (() + Set1d1e, (t-1)) 0, i= 1, ....
Define VT, ( (1) Each component u. (0) is w.. p()) S1 (c) where uf xTT10) x=0)== v/ by
(0)
default. The above can be rewritten
6144196
((+1)
28 + 66418101= (t-1) =
yi ((+1) 2>13/ (1)
(1-1) 0.
+ bi
(27)
where have
follows from 61:441 # O. Note V (1))= = 1-a V VAV
= &AA (0), we (0))) (VAVT.100
14a
(T, (0)
V vT k1y = TT = (y)
where follows from on + O. As (1) = Ay (0), follow Equ. (24) of Lemma C.3, y has the solution
y/ cos(O,1)
(0)
|& < <
10)
[& <
V. 10) X
[A = (
[Ag] A
where 8; arccos(())) A,). We can write down (t) in terms of (t)
VT, = (1 = (.
where Z, has two possible forms
diag (L...., (ostt..,...............
for bipartite graphs
Z,
diag (L....,cos((Et)) ...
cost))).......(((..))) (0x)))
for non-bipartite graphs
Hence, 12,12 < and Z,y ( (1)
||2- We have
1 Va 1 va
(1((() 2| = T < 01 (0)
2
|bbll2,
29

---

where the last inequality follows Lemma C.5 and note 2 = VTr()) To meet the stop condition of
{1rl) < 2ae/ (1 at), u E v} 0. it is sufficient to choose minimal integer é such that
20e
2
||||
(1+a)vd,
To see this, if the above inequality is satisfied, then note for any node we have
Ir S 1lx - | 122 s 2
(a))
||b||2 (1+a)vd, (1 + a)
2ae S
20ev
where all nodes are inactive. So, it gives t
l-vā Itva
< In f by noticing |b1l 2a/ ((1 + ad.d))
It indicates t 2 In 2 In
1+va [-
As
in(1=x) I
1tz for > O, it is sufficient to choose Ł
+ Va 2
==
In
2/a E
C.4 Residual Updates of LOCCH and Proof of Lemma 4.1
We propose the following local Chebyshev iteration procedure
(++1) = (1) + (1 + 01s41) + se:41
- (1044((00=>x11111100
Our next Lemma is to expanding (xx(1) (1(11)) 8.
Lemma C.7. Fort 2 1 with initials x the following
(0) 0 and x( (0) + r.. the local Chebyshev iterative
x (++1) = ((()) + (1 + dea+1)r (t) + Seeei1 (x(() x (0((-1))ss
( Denote A (1) = (t)) - (t--1), then A (1) = Σ)=o (1 + S5g+1) Ir S 11 Spa--1) where
50.1 = 0, So.t = So Si 0 ..... S. and 0jz+1 = 8,d5t1. We have the following
(1 (1+0 (t) + 8 (t) x: (1-11) 8,
Σ
(1+053+)) II
( be
1=0
r=f+1
where Sz e S
Proof. Recall we defined 8y 0 so that dg:1 = 0. We prove this lemma by induction. For L = 1, note dn:1 0 and the support of (0) is Sn (u)p(r(()), then
ax -r) (1 do. (0)
For t = 2, note S = S and So.1 = =G n S₁ by our notation, then
(2) (1) (1 + 5122) (1)
Su.. + 61.2 (x (0) ss
(1) + (1+ + 10:1) 2 (0) (1+6)1))11
So.1
For = 3, one can build (3) (2) based on (2) (1) and recall S2.2 = S2, S1:2 = S1 n
x. (3)
x (2)
= (1 + 82:3 )2:3)r (2)
Sun + 823 n (2) - x (1))s)
= (1 + 82.33 Suz
(2)
+ 62 + )1.2)r (1)
S. + (1 + 8p )1:27 Su.. )s
= (1 + 62.3) (2)
S2a
+ (1 + 61.2) 62.3
S.12
(1)
+ (1 +
Sm
(0)
30

---

Fort = 4, we continue to have
(1) xx()
(1 + 83.4 I
S43 + 63.4 ( (x) (3) x 12(2))ss.
(3)
(1 d3,4 |r (3)
So. + 63.1 + (r S2. (1 + 61.2 02.5 (1)
S1z2
+ 0o.1 061222233
So Sss
(0)
(1 63.4 Sz3
(8)
+ + 523)6
(2) S23
+ (1 + 05:200 Sia 8
(1)
+ (1 + 60.1 011222:3033311 (0)
Son
(xx(4 )
((1 + )3.4 S2.3
(3)
+ (z) - xz )sa)) Sa¹
By induction 1 1,
1-1
1-1
(x() 11 (1 Σ
(1+ 853+1)
II
1=0
r=1+1
where the convention notation Σi= 0 and IT |= 1. To verify the inductive step. consider
for 1. + 1, we have
(1+1) (t) (1 + d (0) + 6 - x (1(11))
(1 1) + 6 x(1-1)) )
1-1 Σ
1-1
1) Spa-1 St
(1 + d (9) + drat1
(1+ )z
8 8 1
r=
Σ
(1+ 81j+11))
II
(1)
Syit
==0
r=
To see the second equation, note
(1 + S1(41) (1) 3.
+ Seet1 -x -x(1-1)) 3
=(1+ + (1+00111)) 01
(1)
+
Σ
(1+8
II
8 8p ()
Si
S5c-188)
r=t+1
Σ j=0
(1 +1 (zj+1)
II
0 (1)
Spee
r=tt1
where recall we denote S₁ P Sp S..
Lemma C.8 (Local Chebyshev updates). Given the updates of (1+1) as defined by LocCH in (10).
we have the following local updates
(a(t+1) =x(() + +(1 l1 + Tr (1) S. + 81.41(2((1) x(t-1)) s.
(28)
p((+1) - -20+1 (r + (1111 Ej (1 dzز+1)) IT- & Qr!)).
31

---

Proof. We only need to show the second equation of (28). The residual of LocCH updates is
p(++) = b - Qz((+1)
((+1 = - (1)((( + (1 dest Ir + 6 (((-1))s))
=r(() - -(1 1o+tt1 1 (t) b 1(=11) -x _1(1-1))ss
(1+2) + (1 + 1 1Qqr d Q(x!!t) -x(1 7()
-1+60 1Qr+ (1) + seat Q(x"1) - x111111)
small noisy part
Note ((x!+) - (1(-1))- = Qx(t-1)) (b Q(x(!1) = (1 (()) and then u becomes
= ((++1) + (1 + ((r+1)Q1l1) + 81:41(1(tt)) - p(t)) p(t
= ((+1) - 25.+1 Wr (1) (x(-1)
where the last equality is due to (1 + 8 1
QQr(1) = (1 + des )r(()
20 Wr(() by noticing
(1 + 5)0+()) 1-a
I-a 28,41 in Lemma C.5. Hence, we have the second equation. To see the noisy part,
note by Lemma C.7
(8) (1+ 1 (1) + de.t1 (x(()) x x(1-1))3 (1 + 6ea4
1)
Σ
(1 (1+ 03g+41 r=
II
+ 1n (1) 15a
S.
1=0
(1 + Q (t) + + - x((-1)) 3. = Σ
(1 835411
II
8p
(1) Qr
j=0
r=j+1
Lemma 4.1 (Residual updates LOCCH). Given 1,) 0, (1) p(0) of LocCH defined in Equ. (28) satisfies
x
(0)
The residual
1-1
vT, 81. V.V (0) + Sintuat + Σ Sptrelt k))kR.
(29)
where
1-1 ΣØ=1 j=1 82.5 Hjs I
1-a A V Lta
Tou /t
ifk 0
ug.t
Ejsk (kk+1.) Hjgt
1+a I A V 1-a
(k) Sn.
/(t k))
2 1,
(-1))) diag (1, cos((.)) ... co...,,.))..........o ((s()))).....c(((..))
diag (1, (1....,cos((::))) (iiag(()...coo())))........
for bipartite graphs
Z,
for non-bipartite graphs
diag - k k..., (t-
(-1)-=h-1 R))
for bipartite graphs
sin 0;
Hy
diag
(1 k ....,
sin(P.(t-R) sin 0,
sini0. (t-A))
for non-bipartite graphs.
On
Proof.
We first decompase the residual equation in (28) as
V ((+1) - 16 T (1) + S (t-1)
(I- 1- 1-1--1) A v 1+a
(1 + (18414))
II (Brart)
vTrg) r
0=0
r=it1
32

---

Define VT = "1 and V (0) = )1.08(0) =y by default. Then we have
11 (++1) 28 Ay 11411y(111) = f(()) y y 2Ay (t) (t
(
Note VTr(U) = =dyy(1 = V 61 Wr = 5n AyAy( which indicates y = Ay Then, by the
Lemma C.4, each yi (1) has the solution
y + 1t), (b) (01.t+1)
if A =
cos (0:1))) (0)
- sin (8)
sin R; Ji / (01:2+1)
if A1 <1
if Use Ze and Hu,t, we write the solution of the second-order difference equation as
(-1)(yt)
(0) (-) k),, / (61
A = -1.
yd Z,y (0) + Σ
HR H(h) 1(011(1))+1)
k=1
V (+) 01,1311) S1 Z. -d1
(1)
1-1
k
V (+) 8.. (,())
2 H
E
(1 (1+8g3+1)
II (8rot1) I 1 1 -0-1 A
V vT) ()
1 1 a
j=0
r=
Note 1 + 85.5+1 = 28,+1 1ta j 2 1, then (1 + 09.5+1)
11 (6pr41) = 2 22
Then, we have
1-1
V (0) 51:1 Z.V (.
+
Σ
6k+2 Has
I 1-1-am A V VTr
S
(0)
1 1+a
k=1
8n 2: 2 Hust E Ö5+1K05+H1AH
1+a 1 A V
()) 15
())
+ 2
1 1-a
=1: D.V. (0) + 51:1
82n H
I (I-1-
1 a
1+a A
V
(0)
k=1
+ 2
2k k Σ 0541K00+11AHH
1-a a
1 1 a I 1- A V
T8
(1)
where the last term can be expanded as
k
2 HRS
05+11805+11441
1+a I
A
T ()
b,k
63. H H2, 12:22:12:220
(1)
S
w
S12
+ 6 H2 13220
S2.2 +
85 H Hs,d220d
w
Sia
(1)
+ 6 H 1,100
+ 85 H 04 w
Saa" +
(3)
8 H H4.224022
(1)
+ 66 H4
+ 66, H 04
S
(3)
+ db. H S5. Ss.s (4)
u (4)
87 Sn.H.,022502000
(1) + 87,
52, + 87 Hr.HH, , S0.6
+ 67, 34.5
SL.
+ 87 10
(5)
5o,5
1-1
8k+1:1
8kt1y
(1 a I A V
(b))
(k)
11 a
j=k
33

---

Here, we denote Ok+1: = 1. The final iterative update is
= =, (0)) + dirl
o2.5 Hjs
1
1--1-a
1 +
1-
A
V
T (0)
r
So.3
16..2
+
Σ
1-1
Σ
1-1
1 + i-a
I -11
V VTr
/ /(-
jmn
C.5 Convergence of LOCCH and Proof of Theorem 4.2
Corollary C.9. Let By be lower bound of residual reduction satisfies |ux,cl22 S BB (x then the upper bound of [|r(0)]]2 can be characterized as
[-1
|||00/2 S 51 II (1 + B,)Ue, where ye+1 2ye + 1=0
yt-1
0,
(30)
(1 + 1)-1) +())
where yo = 1/1 ||r0
Proof. Since ||ux,el22 S Bx [|)) the final iterative updates (29) can be bounded as
(Ty = -.. (0) + 51 + 2
Σ
S na(l k)uue
k=
||r 122 < 81se|lr(0)11 1|2 + 61.1180l111011 (0) +
d - k:)32l1 (
t-1
(|r(()
112
6x- k) S 61 1z0
(31)
where 0, 1...., T. These T + 1 (including a trivial one where
||:00)|2 S ||r(0))|2) inequalities
shown in Equ. (31) form system of linear inequality matrix
1
0
||>(2) ll2
81.1 + +10))l1
221
0
0
||r(3)
61:(11 + +220)11((01
231
232
1
0
||r63)
S
61.(1 136)))1((011
C,
2T1
225
1
||<(T))
Ghr(1+T3h) + (r)0))
I 7
where (Zz) = 26x+1a (t k) for = ....... T and k = 1.2...., 1. Assume that N E R TxT is strictly lower triangular matrix, then we know the established formula (I ++N)-== I - Hence, we have the following
(I- 1-Zz))--14 -I+
Σ z
k=
34

---

||n00
||r(2)112
Given that (I- -τ) 2 0. then we obtain an upper-bound of
SI Zz)
||
leads to the following new second-order difference equation
2 2
da+1a - k. 81:1(1 + to:
for t = I......
1=
where the initial value of 20 || 2|.. Following the argument in Theorem of Golub & Overton [18], we construct second-order homogeneous equation for 24 as
= (1 +
1a k)
2+1 61:41(1 (1((1+1)))))))))
8x+1 + 3k
81112 = -111+1 + Σ
8u+ 1e+1144(1 B4
k=1
2++1 6e+12 = 61:4+11832 + 2
1
(32)
k=1
where Equ. (32) is obtained by the difference between the second equation and the third equation. Similarly,
žt- = (1 +2
8at 1:tl 11 11- - 1)32zz
2-2
8 (1- (1) +
Sx+1 - - -1)3822
1-1
d1+12 - 6144121-1 = 81.4+10020 +
>
8x -x-1:1+1 1:0+1 ByZk.
(33)
k=1
where Equ. (33) is obtained by the difference of the first two. Hence, Equ. (32) (33) gives
_-- 2(1 + 3:)6 žt +
=0 for = 1,2,...,
where two initials are =n = [0) l1o and 21 61(1 + (B0)ll((0))) ||2. Let 2, 81 then 111 2(1 + ++) + 1 = 0,
where two initials are ån = 20 ||r and 21 = + (0: We finish the proof by setting (1 + By)y= 2.
Remark 1. 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 B3 O, the following second-order difference equation
Ft+1 2(1 + Br)az, + It-1 =0.
has the following solution
(-1 Ft
II (1 + 8,)y,, 1=
where U1++1 2y + 1148-1014
Mt-1
0 with yo and 91 I1 11/(1 + (0).
35

---

Proof. Assume (--2) [[}= (-2(1 + B,)) yr. Then, following the equation, we have
[+
(-2)
II (-2(1 + B,)) ye+1
2
j=U
1-1
1-1
- -(11+ + +())
(-2(1 + 3,)) w +
(-2
II
2
j=0
j=0
Since 8j 2 0, the term II, == (-2(1 + 83)) # 0, we divide it on both sides to have
(-2
1 6-1
(-22)
1
(-1
36+1 +
2 Vt +
4(1 11 1111+00 3,) 31 0.
Hence, is simplified into yt+1 2tlt + prove the lemma.
0. To make a simplification on 2t, we Vt-1
(1+P--1)(1-8
Theorem 4.2 (Runtime bound of LOCCH). Given the configuration 0 = s, 9) with E (0, and S 1/d, and let (T) and (T) be returned by LoCCH defined in (10) for solving Equ. (3). For t 2 1, the residual magnitude ||r(() [l2 has the following convergence bound
|r(0) 12 S E1
II [1 (1+8,)88r. 1=0
where yt is a sequence of positive numbers solving V++1 2yr + yt-1 + 9(-1) + 3:)) =
with you = y1 = -1|r(0) [|2. Suppose the geometric mean 3, & IT + 3;)) of 3e be such that
evaa 3, 1+ I-va where e E (0,2). There exists a real implementation of (9) such that the runtime
TLoCCH is bounded by
(u+ (L+va]vo((yp) In Va((2-
2yr
TLACC ~0
Proof. The convergence bound of (+) directly follows from Corollary C.9. Since we assume that
evaa there exists E 2) such that 11 - + B₂)) S (1 + ( Then multiplying both sides by at, we have
(1 11 (1+))) (1
(2 - c)
ā'
1 va
Then we have
[-1 II + Bi)u
Bue52al
|| l S d1s
2ā 2ā 3.y. S S
==
1 1- Va
aI B,) 11)+ )
E
Aln
< In
1 Va
23
1=0
2y/e
1 va
2 In
In
(1- Va
Since 3,- = (II)=b(1 (1+ + ++,)) 1/t
and by using 1+x 2
1 and letting
= Itya (1~v0)
1 > 0.
In(1+x)
then can be lower bounded further by
+ Va
2t
2y
1 a
t 1+ a (1 v
In
>
In
e
(1 va)3,
36

---

Since we assumed B (1 eva Then, we find such an upper bound of so that LocCH converges.
) which means 1 S 3
(1
1-val so B E
eva
1, 1
I-va
1 vā 1+ 1+va (1 va)3,
20
1 V
2y ]
In
In
C.6 Implementation of LocCH
We present the implementation of LocCH as follows: Recall the sequence ói-1 = (2
6.) 1 = 1, with 61
1-a 14a"
Denote a(t) & xl1) (((-1) (c) == (1 + des+1))(())
Sus+18 a(1), have
x (++)) = (1) + (1 + 81s+1) 10 + 1 = (9) + + (4)
S.
(1+1) = =b Q (() (1 + :1 )rs? 11e4 (1) p(() -QA (8)
((+2) = (1()
(t) (t-1)
+
5,-1
When t = 0, we have x (0)
= 0,
r (0) b, a(0) 0
(0) = (0)
When 1, we have (!))
T
(0)
(())
W
S.
(0)
a(1)
T So
(0)
A (1
61)2)((1) (1))
When t 2 1, we can recursively calculate the following vectors
(++1) (1
(t)
S.
r (1+1) p())
(4)
1
WA (1) S.
+
(141) :0
(t-1)
S
Therefore, at per-iteration, we only need to save sub-vectors As, and As and update æ locally,
D Local Heavy-Ball Method - LOCHB
D.1 Standard HB and Proof Theorem D.2
Lemma D.1 (The standard HB updates). The updates (! and (()) of the HB method for solving Equ. (4) can be written as
(1+1) xx) + (1+ (1) + ā² (x! - x((11))
(1+1) 2 (t-1t
The residual updates can be rewritten as second-order homogeneous equation
y (++1) 2Ay 10) (t-1) = Vt 1,2,3.....
where (0) is such that (1) à' V (t) t with (0)) VT, T
Proof. We follow the standard Polyak's heavy-ball method [40] as
x (1) The V (xx + 7/ (x. (0(t))))))))))
where f((((t)) = Qx b and na = 4/ 20/ 2
1 à and 773 = 2/ a) 2o (1 + (1)) + + / 20 a))² = Hence,
leads to the following updates
((+1) = (t) + +1+ (0) + (x x (t-1)
37

---

Inserting
r()) = Q(x" b (1 b- 1- w 2 = )
2ā
W ax(1)
I+a*
1 a²
then a (1+1) 22WW (1)) -2x(11111 (1 a²2) ) and since
QW 1 - 1-Itaa 1-1--w) w = w (1 1- 1 (1-1- W = WQ
+a
So
((+1)
Q(x((+1) -xx)
2aQ (1) a² Qz((11) + (I (1 Q)6
-20QW (') à² 2²Qu (1-1) (-a² 26W)
26Wr(( (t-11
Using (() = (1,
ā Vg (1+1) 2åW a' - -2d
Vy 2WV3 (y((-1),
As W = VAVT and V V is orthogonal matrix, we continue to have
-2 V -VV AV Vy + + (t- =0 ((++1) 2Ay (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 as
(++1) = (1) + (1+ ((()) + â² (x - (1(-11),
(t+1) 22W (.t)) - (t-11t-1)
where the initial condition is x (0) = 0, (0) b, (1)) (0)) + ()) p(L) = Qx
Then there exists a constant J such that the total iteration complexity to reach the stop condition {u: (rul & edg, u E V} = is
[ Va Cr t
In
2
e
where C. = -W) aW) (ideal case); Ct max
1ta=
1+ 1f
VT-x
(practical case); and C, VI-xg iT (1-8)(++)) I (G is not bi-purtite graph).
Proof. Recall W VAVT and then V (1+11) 22AV 1(0 T (1-11) By Lemma D.I,
we have
y (+11 (1) + y (t-1) = 0,
where we obtained second-order difference equations
(1+1)
2A:1) (t) + y (t-1) = 0, vi= 1;,,,...,,
yi
Follow the Lemma C.3, Equ. (23) has the solution
sin(A, () ( (t (1),
A& < where = arccos(\;)
w/ (1)
sin(e,)
(34)
(1( (0) + (y) Xa X DE)A
A
where in the case of |As] < 1. We consider the three cases of r
38

---

Ideal Case: We can eliminate / in (34), when (0 Ay we get v/ = Xw,, (0) and then y (0) can be simplified into
(ii (t-1))) (0)
[1d 1 y cos():t) [1d = U 1
11/< <1
1
Pa. <1
yi
sind.
I
X
| =
In this case, I needs to be r = Q--Q--111 - a.). Therefore, we have
10 1|2 = || = a' 1|2 S a" (0) 2|2 = a'll V T 11 d'(|r100//2
Practical Case: Just letting x!!1)
x
0. we have ay = (0)), then
Vi (0)
1 (1- 13)
(0m(
(0)
| <
1+a 1 +1+ 11+t111
)) 11:0%
I max
(1 + (ā );) (0) X
A|=1
VV- X3"
where 02 arccos Az).
Non-bipartite graph Case: When the graph is non-bipartite, we can eliminate t, as the following:
We choose r = -I, we have
(1) (0) (r) = (r) ( = b -Qz!1) = (0))= =Qr (0)
VT, (1) = (1- (0) (1 (1-oj) O)T AV (0))
vT.(1 = (1 + (1 a)T x,)er (0)
+
1
We have the following relations
v.. (0) 1/ of r1) ànر²) (1 T1 (1- a)T (1-a)T
(1 --))dd (0)
11),, (0) = (1 7
13;
1 Ita a
1 1+a a
To make disappear when A; = L, we need w/0) y (1)
1 (1 -a)
1 Va
T
1+a 1 a
a + Va
In this case, we have
19 <
2
(0)
yi
V1 X
To make sure the algorithm stops when the stop condition is met, it is enough for
|| 1|2 ā' || 19 1l2 S a' y ( 1|2 = a' V'C, , (01|2 = a' Call ||2 S e.
This means (1))0' e, which leads to the following
1+ Va
2/ In C/ll
t =
Remark D.3. The constant that appears in the bound involves the second largest eigenvalue X2 of AD-1 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 directly bound (() instead of providing bound for (t))
39

---

D.2 Residual Updates of LOCHB and Proof of Theorem D.6
Lemma D.4. Let the local heavy-ball method be defined as
(1-1)) x1(0)
(9) Se ((+1))
-QA
(1) S.
A 10) (1 (1+ + a²(x!!)) xll-11),
where (0) 0, x(1) (r) and r 1, F') is initial step size matrix. We have
the following expanding sequence
a² (x x! (((-)))ss =(1
Σ
á
2(
p (1)
Tr (0)
1/ 2 1
1n5.
i=1
A (9) (1+
2(
(1)
ā Tr (0)
Vt 2 1.
Sun- 15
Furthermore, we have the following sequence
@²(TT a AV (xx(1) - =0(1-1))3.
I+a
(L+ (V
1-1
1- AV
(0)
+Ta2
1- AV "'."""
1+a
S.,
1+a
where we denote Si,, - S n St
Proof. We assume all nonzeros in b are active nodes time t and 1, i.e., So =
supp(b). The local updates can be expressed as
(++1) (t) + (1 + (
+ .2 (2)(((() (1(1-1)) St.
((+1) Ox((11)
(t) (1 + ā 2 Qr (t) - a² (z(((=(1) -c
original updates
+ + 8² Q
8 Q ( (1 11)
noisy with small magnitudes
2åW (1) - (1-1) + (1+ a² Qr 19) + ((( - (1)
3.
noisy with small magnitudes
(1 (1+ (111 -
1
1+a -/ + Q(x)(!) (2(1-1))
W (8)
-
(((+)) -2a + (1-1) (1 + _ (t) 26
(
aa² (1(-1)))
s,
40

---

Fort 2 1, we can expand x (++1) x(() as the following
x (0) + (1 +1+ + (2((0)) -x(1)) 8,
(t) (f+1) x (1) = (1 [1+ (1+ (1++))) +0 -2(-)) 5 (
(1+ (
- 211
St-
= =(1+ + S. + (t-2)
S1-2
+ +
S
= =1+ 10 + 0² (1 + ) S.m + + + - (t-2)
((-1
+ t x
2) x (-3))
S1221
= (1+
Σ
à 20 rr Site + ā (1 - (1(-33) S1-2.
1=l-2
(1 + a²)
2( (1-1) (1)
Sue + =2((2((1) =x(")) Si
-1++a²
21
(i) (0)
Tr
Sun SI.C
Note S. = (u) then (0)
(0) Sozs
S..
S continue to have
1-1
(1) (t- = (1+
1-1 ( (1)
S
+ Taa 1), (0)
r Soat11
Vt 2
-2(x(0) - 1(1-1))33 = =(++
(t- (1)
+ r (0)
r
14 2
Sa 1
Sea 155
The rest follows readily.
Lemma D.5 (The nonhomogeneous difference equation). Given y Ay equations
(1+1) - -2Ay (t-1)) = ((1)
have the following solutions
( = Z,y
H (())
where
diag (0s..,,, (-1)9)
for bipartite graphs
Z.
diag (Os()),, coss(An!)
for non-bipariit graphs,
diag (1 - k ...,
sin(8,
(-1)a-t -- k))
for bipartite graphs
sin 0;
HE
diag
- ky..,,
(t
sin (t-A))
for non-bipartite graphs.
sin 0,
siu 6
Proof. This directly follows from Lemma C.4.
41

---

Theorem D.6 (Representation of r() for LoCHB), Given 1 2 1, xx = 0 and x+ residual of (()) of LOCHB satisfies
Trrs. [0) The
VT,(0) = a'Z T (0)) + à't
á H V
Qr
So
no,6
2
a'- k)
2 H
1 a
1+a a
A V
/ (l-k) k)-
Su.s
(k)
Proof. Follow Lemma D.4, we have
((+1) Q4 (t)
1( QA
(0 -r - - -x01-1)
-a'p(1) +2WW, (1) - ā² 2()(2(1)) - -x((-1)
20W, (1-1)
So we can write the updates of LOCHB as
(l++) - -2äW (() + (t-1( QA S.
(1 + 02)
à
2(- 1) Qr
+ -r (0)
S.,
So..
1=1
Write (T) = (") and note
(00 = a' V'V (1)
( S.
= ā² (Vyg2))
V (t)
= =à'VT (Vyli))
=
T.
3.
Hence,
VT ((+1) - -2āv (7(() (1) + ==aV (T
=(1 V Qr (1)
s.. +ā V @rr 0
So.
(1t1))) -2āā W (0) + aadat-y(1--
= = 0²) 2 y()),,, + +a"'QrT((VyV) ) so.
= y + ty (-1)
Σ (1+a2ja' (V S. + QTVTQTVVVV y (0) So.
$6
Then, from Lemma D.5
(0)) - ()) ()) +
Σ (k) (k)
- (T = T' 10) + ā
H H (1)
k=1
k=
42

---

so expanding the error term
(H Σ
H ( Hz.Ho
1=1
= =11
Σ H QV (Vy
- QTV (Vy S5n
= (1+
Qv
10 S. + @ H H
(0)
So.
==1
2
t+ k- 2i H
1ta
QV
S. +a
(i)
+ H2 (0)
Su.k
1=1
Σ a' Ha.f (k)
t-1
H
1+
QU
(7)
+
HIV (0) Tu.
1-16-
(-1
2
サトナ 21
H
+ a QV
Sk,)
(R)
+ Σ
gltk-1 2
H QrV T
Sp.k
a
A= j=k
which when simplified, gives the relation proposed.
D.3
Convergence of LOCHB of Proof of Theorem D.8
Corollary D.7. Define
Bk, - ||ux,ell//|1 1|,,
Bg a max Br,e.
Then the upper bound of ||)(() 1l2 can be characterized as
f-1
||~ ( 22 S á' II
(1 + B,)yt,
(35)
5=
where Uo+11 -2ye + 3t-1/((1+ B + 3.)) where yo y1 ||>((0) 1|2.
Proof. Since 1|ux,el2 s Bx (x) 2. then given the final iterative updates (29)
(T) G'z.V (0) + ā āttuo, +
Σ ā
(t k)uE
k=1
and since |Ze < 1 we can bound
[-1
(t 2 S a² (0) a't +
- k- k - k))x1l1((0011 (10
+ [|r(()
Σ
S ā' + t || (0) 1|2,
(36)
km
where 0, .,.... T. following inequalities
The rest just follows similar strategy shown in Corollary C.9. We have the
0
[|x
a' + 36) (0)
121
1
||x(2) ll2
6 236) (0)) (002
14l + 3/6) 33)))(1((0))
I Vz
U31
U32
0
||
C,
:
||
ā7 (1 + T3)) (3 (001
Unl
Un2
Un3
43

---

where (VL)x* = 2
k) km Denote each upper bound as = ||r(0//22 we will have
Cx +
Ur - a1 + 2
k)ā "a)"""BBn BETTE
k=1
k=1
U1+1, at (t 1),000 + Σ (t +1+1 - k))atktt 3K
7++1 G+1
1-1
- (1 + (t+ 11)0000 +2 Σ
(1 - k)a² k)aa-
k+ 36TH + 2 Σ at- k+1 BETH
k=1
1-1
αTg = ā ++1 (1 + TO + Σ (1 - -k)ā -++++=+111 BK
k=1
T(+1 ÔT; = á
g++ p0 2
á
xt-++1 BKTh
k=1
2-2
T-- = =--1=
a'-1 10 +2 Σ
(t k -1)
1)aa=== 1kBTT-
k=1
t-1
[-1
T' T-- T
at (1+ (t 1) + Σ
(t - k) a 1 36TH - 2
Σ
k=1
1-1
ãT, à T--1 = ā g++ To + 2
Σ á zt-k+1
BKT..
kml
T++1 2hT, + à T(- B
The above analysis finally leads T++1
2(1 + +.):aT. ā²T-11 = 0. So for (-1
ā'
II (1 + By) Tt
we have
yt
yt+1 - 2yr + (1+ + B,) + 3,-1)
Following the same strategy in Corollary C.9, we obtain the upper bound.
Theorem D.8 (Convergence of LOCHB). Let the geometric mean of 3r be B, A (1 (1 + B,
Then the upper bound of [|r(() l| can be characterized as
1-1
||00(0) 1|2 S a' II + (1++)))
(37)
3=0
where Yt+1 - 2yt + yt-1/((1 + 11 + 3()) = 0 where yo = y1 = (1)(0)) 1|2. Assume that there
exists constant e. (0,2) such that Br < 1 v bounded as
Then the total number of iterations can be
l-va
1 Vo
y/
T
In
(2 c) Va
Then the total runtime T is bounded by
( T
(1+ (1+voool(Srr
In Vt
vol(Sr)
(2-c)va
22 clvo
(2
c)vo
where Õ hides In(yt/e).
Proof. The first part follows from Corollary of D.7. The rest follows the same strategy of LocCH as in Theorem 4.2. (See also Theorem 4.2.)
44

---

Table 2: Examples of sparse linear systems
Our target Qu b
Original Linear system
[u, L]
Ref.
Q A oD
AD
A 1,0 -aa
[27]
(1-oW) x ae,
[a,2- of
b a63
[53]
[2]
-0 (al + -- C) aD 1/2e,
A 1, 14a
1-a
[13]
2a 2
C C= C-I-D-11 AD
b 20 D (
Ita I+a
[12]
[37]
A 1,0 1
(I- ae2 3.
b = aD aD-1//2e., e,, a E (0,
[a, 2 a a|
[7]
[56]
a D 0-1/2y
A "
D D In, 0=
[41]
I + D A) y 2
MAD-1/22.0
X E [I
A 2+
Fallaax
[57]
x D1/2y
D.4 Implementation LOCHB
We present the implementation of LoCHB as follows: Recall the updates of LoCHB is (++1) (*) (1 + ā² rr. (t)
+ à² (1)
x x((-1))ss) ((t1) = 2äWr 2
The corresponding local updates are
(4)
(1+ a² ))) ())( +a²
= (1) + (1)
x
(1-1) = (((
(9) 1- a (A())
1+a
(t-1) (t) + (7) A (t-1)
where if we choose (0) = x (1) ã
0, r (0) (1) b and A (0) 0.
E
Instances of Sparse Linear Systems
E.1 Table of Popular Graph-induced Linear Systems
This section presents most commonly used graph-inden induued linear system as the following
D 1/2 122 AD-1/D11122 1/2
b,
where Q is the generalized version of the perturbed normalized graph Laplacian matrix with perturba- tion parameter > 0, and b is sparse vector. A typical example of Q 1-12 AD-1/2 1/2 with A- and = D-/D2 e./ ar)
I-a
T+a
The detailed parameters are:
W = --1-2AD-111 and Ā + 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 aa is the perturbed augmented normalized Laplacian with perturbed parameter a.
2. Q = 112 (D
I- (D + A)) 1/2 al + 1-a 2
and Q = 1-a 14a
1/2 AD 1/2
2
This known as the lazy random-walk version of PPR vectors.
3. (I - (1)-a)AD AD-
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].
45

---

Dataset ID
Dataset Name
1
m
91
as-skitter
1,694,616
11,094,209
92
cit-patent
3,764,117
16,511,740
G
com-dblp
317,080
1,049,866
G4
com-lj
3,997,962
34,681,189
95
com-orkut
3,072,441
117,185,083
com-youtube
1,134,890
2,987,624
97
ogbn-arxiv
169.343
1,157,799
Ge
ogbn-mag
1,939,743
21,091,072
Gg
ogbn- boon-products
2,385,902
61,806,303
910
ogbn-proteins
132,534
39,561,252
G11
soc-Ij1
4,843,953
42,845,684
912
soc-pokec
1,632,803
22,301,964
913
wiki-talk
2,388,953
4,656,682
G14
ogbl-ppa 576,039
21,231,776
91s
wiki-en21 6,216,199
160,823,797
91
com-friendster
65,608,366
1,806,067,135
917
ogbn-papers100m
111,059,433
1,614,061,934
Table 3: Dataset Statistics
4. Graph kernel computation for online learning. Each computed vector serves semi-supervised learning feature vectors [24] or as online node labeling learning vectors [41, 57]. Note the target linear system when =
D D-1/2 1/2
In D A) D 1/2 112 2
\ D D-11/2 1/2 11/AAD111 1/2 Di 2XD 2XD 1/2e...
11
Hence, we have A = D + In, = 1, = 2/2-1/2e..
F Experimental Details and Missing Results
F1 Datasets and Preprocessing
Following Leskovec et al. [34], we treat graphs as undirected with unit weights. remove self- loops and keep the largest connected component when the graph is disconnected. After preprocessing. the graphs range from 169, 343 nodes in ogbn-proteins 111, 059, 433 ogbn-papers (0M,, as presented in Table 3.
F.2 Problems Settings and Baseline Methods
For solving Equation (3), we randomly select 50 source nodes s from each graph. The damping factor is fixed at 0.1, i.e., range (0.005, 0.01, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3) for others. The chosen from the range [2a/((1 + a.d.), 10_4/n].
0.1 for all experiments, and it varies within the
For solving the local clustering problem, we follow the greedy strategy from Andersen et al. [2]. where local cluster is identified by examining the top magnituees in PPR vectors. Specifically, denote the boundary of S as O(S) {(u, E u E S, v & S}. The conductance of S is defined
as
[A(S)
((S) min(vol(SS, 2m vol(V))))
46

---

dipelet
Oat Seedil
oml)
-11
-1P
XBF
x30
XDF
con yunuu
ogh roghro
agh prevolo
-00
-1.
-1A
xay
x
x1e
xbr
xIpe
ophrpeddats
scold
whben21
wibbtak
-13
-05
25
18
opeziiiss
Expecctios ktt
opecrioss kiv
Figure 7: The LOcSOR method compared with CGM over 18 graphs.
patent
oom ctip
comii
orkul
500
000
comn-ynuuube
ogiloops
ogin--rch
oglnnrrag
ogbh-procucts
60
500
1.90
40
50
25
wik-on21
Spaedup
100
10
000
30
Figure 8: The speedup of local solvers compared with their standard counterparts.
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 - in decreasing order of magnitudes. Let the ordered nodes be v1, .., Uz; the local clustering algorithm iteratively checks the conductance reduc- tion by 01; U2 ... Ut where & 1, ... I, and after completing all checks, returns subset U1 02.-... Uh that has the minimal conductance among all examined subsets.
Parameter settings baselines. For the local ISTA method [13], the precision parameter is set to é é= 0.5 for all experiments. According to the algorithm's 's description of ISTA, the corresponding value is given by c/(1 + é). For LoCSOR, the parameter is calculated as 2(1 + a1 +1+ va)² 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, to ensure the identification of good approximation.
0.1/n, to high precision, é 10-1/n,
For our experiment, we used 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.
F.3 Full results of Fig. 15 5
In 15 graphs, we set a 0.1 and 0.1/n. For each of the testing graphs, we randomly select 50 nodes and run LocSOR and CGM.
47

---

(e 12
paleez =1.13-00)
dil-lbl 12e00)
cotrli 11.12000
colr
LocSOR LocHi LocCIl
-2
-2
2
0.0
0.3
1.0
2
x10
X10
x10
X10
X10
comcyoulaale 12600)
ug agbippsn (t =1 13e-00))
udd beraaiii (e =1 =1.12-00)
uglu-unn (e =1, 12-00)
geberpeodd (e =1 -2.120
-2
In Mllll
0.0
2.5
5.0
0.0
0.5
1.0
X107
107
10P
x
X
oplmn-pmttiiee (e =1 126-t00)
soc-Ij1 (e =1 126-20)
socepoloe (e =1 12-00)))))
wiki-en21 (e =1 12h401)
w MlC-talk aik de m1 126-20
-3.5
-2
In
-4.0
-.5
Operations 10
Operatioos 10
Operations
Operations x
Operations x107
Figure 10: Comparison of three local solvers 15 graphs.
Fig. 8 presents speedup tests on datasets. evident that these standard linear solvers can be localized effectively.
-7.11s-001
coza-dlip
Lockon
Looks
X127
xtdf
ousepostale is -2.36-00)
ca24415e (s (00)000)
1alis-endl (6.600))
oplasans (s -i (860000
rgloggroiats -1012-00
calospediita -2022048)
mooli m1
=1.4e-01)
wilo-e222 (e m4.32-30)
valo-tak m1.
-55
-00
-33
4.5
0.2
65
1.0
1.5
0.0
255
1.0
Operatins
Operstsos
Operatione
Openatine
Operation
Figure 9: The estimation error reduction tests on 7 solvers including our LOCSOR, LoCHB, and LOCCH. The experiments were conducted on 15 datasets.
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 Fig. 10. Empirically, LoCSOR is the fastest algorithm when the parameter is chosen optimally.
48

---

(Nocl
(Nodie
1/4))
-70
LocitB
2.25
-2.4
200
LocCH
LocSOR
-200
In M0lla.1
-3.0
-270
-2.
-275
-30
-0.5
107
x10
x107
x100
x10
(Node 0. 1/m)
[Node 1, 1/m)
(Node 2, 1/m)
(Node 3. 1/m)
[Node 4, 1//e)
-20
-2.5
-30
3.0
3.0
-5.0
In
-3.5
-35
-5.
4.0
1
2
Opraalios xII
Opentiooa
X10
Operatios
x10
Opemntion
X3I
Openitions
XII
Figure 11: Estimation error as function of the number of operations on nom-friendster. We randomly select different nodes and use = 1./ and = 1/m.
(Node 0, 1/m)
(Node 1. 1/m)
(Node 1/m)
(Node 3, 1/m0
[Node 4, 1/me)
-15
APPR
LoeCH LocSOR
-20
-2.5
2.8
-7.
-2.
-500
-3.0
-8.0
-3.2
0.0
2.5
5.0
x10
x107
X307
X207
X10
(Node it, 1/bm)
(Noce 1, - 1/m]
(Node 2, c- 1/m)
(Node 3, c 1/mm)
(Node 4, 1/m))
2.5
-25
-30
-3.
8.5
-35
-.0
0.0
0.5
1.0
0
2
0.0
0.5
1.0
Operations
x10
Operations
x100
Opentions x10
Operations x.I
Operations
x10
Figure 12: Estimation error as a function of the number of operations on agbn-pappes10000 100m. We randomly select different nodes and use 1./n and
1/m.
Run time (seconds)
Number of operations
Dotaset
LoCSOR
LocCH
CGM
LocSOR
LocCH
as-skitter
0,350 0.054
0.567 0.095
2.626 + 0.551
7,827e+05
1.026e406
1.524c+08
cit-patent
0,966 0.120
1.609 0.189
14.660 +: 2.249
1.804=006
2.298e406
2.873e+08
com-dblp
0.068 0.058
0.104 0.059
0.469 + 0.112
8.222c+04
1.166es05
1,562e+07
com-friendster
15.29 1.89
26.54 3.52
508.50 99.12
7.027c+07
8.063c+07
1.442c+10
com-lj
2.361
4.122e+08
3.221e+06
9.220e+08
com-youtube
0.290 0.070
0.501 0.099
1.314 + 0.257
7.617c+05
9.323e+05
3.561c+07
mgh-mag240m
85.14 1: 16.05
108.86 = 7.994
549.51 367.2
4.541e+07
5.820e#07
9.426e+09
oghl-ppa
0.116:00000
0.202 0.040
5.624 + 1.108
6.117c405
6.445e405
1.682c408
ugbn-arxiv
0.039 +: 0.060
0.058 0.063
0.239 + 0.104
1.158e=05
1.354e+05
1.473e+07
ugbn
0.520 1: 0.079
0.921 1 0.128
6.136 1: 0.974
1.804e+06
1 1,947e+06
2,050e+08
ogbn-papers 100M 100M
27.20 3.94
41.99 6.17
750.96+ 10045
8.893c-07
1. 1,055c488
1.604c+10
ogbn-products
0.695 0.236
1059 + 0.249
31,907 + 4.961
1.777c+06
2.385c406
7,071c+08
oghn-proteins
0,021 0.057
0.025 0.057
0.910 +: 0.306
7.941e+04
6.610ex04
1.488c+08
soc-Ij1
0.487
7.102 +
3.263c406
4.045e+06
4.833c008
soc pokee
1.568 0.207
2.020ee006
2.225e+016
2.160c+08
wiki-en21
1,436 2.708
0.215
19.658 4,576
5. 5.996c=06
6.659c406
1.329c+09
wiki-talk
0.251 0.049
+ 0.091
0.642 + 0.071
1.090c+06
1.290c406
4.284c+07
Table 4: Summary of runtime and operations for 15 datasets. (e 10 10-6)
49

---

ogôl (a 0.005)
agblippa ppa (a 0.1)
agbi ppa (o 0.25)
-21
LocGD
LoxSOR
-25
LocHB
-23
-3.0
-2.4
2
2
0.0
0.5
1.0
1.5
X107
x10
x10
ogbn-products (ir = 0.005)
ogbn-products (a 0.11
ogbn-proddcts (er = 0.26)
-5
-5
-6
0.0
0.5
1.0
Operations
x10
Operations
110"
Operations
x10
Figure 13: Estimation error as a function operations needed. For a 0.005, a 0.1, and 0.25.
puoee
Clellcceas
-03
Lee'l
-2
1000
100
0/
-05
-0.
-1.
359
600
154
300
2003
4000
23
0.0
0.0
-0.5
-2
-10
200
1000
1331
50
100
200
2540
5009
7501
1000
2030
3001
Support
Suspect
Support
Suppoct Sto
Figure 14: The graph conductance found by local graph clustering method using different local approximate methods. Experiments ran 15 graphs. 10-6)
50

---

Dataset
APPR
LoCCH
LoCHB
LocSOR
ISTA
FISTA
as-skitter
5.90e-04
5.90e-04
5.90e-04
5.90e-04
5.90e-04
5.90e-04
cit-patent
4.23e-02
4.23e-02
4.23e-02
4.23e-02
4.23e-02
4.23e-02
com-dblp
4.12e-03
4.12c-03
4.12c-03
4.12e-03
4.12c-03
4.12e-03
com-lj
2.94e-04
2.94e-04
2.94e-04
2.94e-04
2.94e-04
2.94e-04
com-orkut
1.75e-01
1.76e-01
1.76e-01
1.75e-01
1.76c-01
1.75e-01
com-youtube
5.46e-03
5.46c-03
5.46c-03
5.46e-03
5.46e-03
5,46c-03
ogbn-arxiv
2.11e-02
2.11e-02
2.11e-02
2.12e-02
2.14c-02
2.11e-02
ogbn-mag
3.18e-03
1.59e-03
3.18e-03
3.18e-03
4.74e-03
3.18e-03
ogbn-products
7.26e-02
1.02e-01
9.68e-02
9.65e-02
9.51e-02
7.24e-02
ogbn-proteins
5.92e-03
5.92e-03
5.92e-03
5.92e-03
5.92e-03
5.92e-03
soc-Ij1
5.07e-02
1.18e-01
1.18e-01
1.18e-01
5.25e-02
5.07e-02
soc-pokec
7.80e-05
7.80e-05
7.80e-05
7.80e-05
7.80e-05
7.80e-05
wiki-talk
1.64e-02
1.64e-02
1.64e-02
1.64e-02
1.64e-02
1.64e-02
ogbl-ppa
5.13e-02
5.13e-02
5.13e-02
5.13e-02
5.13e-02
5.13e-02
wiki-en21
1.31e-01
1.31e-01
1.31e-01
1.31e-01
1.31e-01
1.31e-01
Table 5: The local conductance for six local solvers tested on 15 graphs datasets. (c= 10
Dataset
APPR
LoCCH
LoCHB
LocSOR
ISTA
FISTA
as-skitter
0.127
0.147
0.144
0.043
0.323
0.093
cit-patent
0.362
0.516
0.457
0.125
0.939
0.308
com-dblp
0.069
0.033
0.028
0.014
0.297
0.042
com-lj
0.357
0.440
0.664
0.175
0.493
0.229
com-orkut
0.072
0.141
0.139
0.055
0.108
0.084
com-youtube
0.128
0.176
0.131
0.040
0.682
0.102
ogbl-ppa
0.102
0.047
0.091
0.027
0.146
0.102
ogbn-arxiv
0.042
0.013
0.014
0.006
0.237
0.032
oghn-mag
0.068
0.137
0.091
0.035
0.253
0.090
produtroducts
0.072
0.121
0.117
0.045
0.135
0.090
ogbn-proteins
0.005
0.005
0.005
0.002
0.005
0.004
soc-Ij1
0.239
0.376
0.350
0.103
0.512
0.194
soc-pokec
0.067
0.096
0.059
0.033
0.222
0.051
wiki-talk
0.197
0.314
0.274
0.109
0.508
0.176
wiki-en21
0.121
0,290
0.279
0.106
0.197
0.147
Table 6: Runtime (seconds) for six local solvers tested
15 graphs datasets.
10 10-6)
51

---

E4
Results on local clustering
Dataset
APPR
LoCCH
LoCHB
LocSOR
ISTA
FISTA
as-skitter
6.9e+05
7.6e+04
8.1e+04
6.5e+04
2.9e+06
5.7e+05
cit-patent
6.7e+05
1.0e+05
1.le+05
8.9e+04
2.3e+06
4.4e+05
com-dblp
4.3e+05
4.8e+04
5.0e+04
3.5e+04
1.9e+06
2.9e+05
com-lj
5,7c+05
8.9c+04
1.0c+05
7.6c+04
2.8c+06
4.4c+05
com-orkut
5.4e+05
8.9e+04
1.le+05
9.0e+04
1.3e+06
5.0e+05
com-youtube
5.3c+05
6.5c+04
6.6e+04
5.7e+04
3,8c+06
4.4c+05
ogbl-ppa
6.7e+05
9.6e+04
1.le+05
9.3e+04
1,5e+06
6.0e+05
ogbn-arxiv
6.8c+05
8.7c+04
9.7c+04
7.6e+04
2.7c+06
5,2c+05
ogbn-mag
5.2e+05
8.2e+04
9.4c+04
8.3e+04
2.4e+06
4.4e+05
ogbn-producte
9.0c+05
1.1c+05
1.3e+05
9.8c+04
2.1c+06
7.6c+05
ogbn-proteins
5.3e+05
4.0e+04
5.5e+04
5.0e+04
6.8e+05
6.2e+05
soc-Ij1
5.6e+05
8.9e+04
1.0e+05
7.7e+04
2.8e+06
4.4e+05
soc-pokec
5.9e+05
1.le+05
1.3e+05
1.1e+05
2.5e+06
5.0e+05
wiki-talk
4.4e+05
5.5e+04
5.2e+04
4.6e+04
1.8e+06
3.4e+05
wiki-en21
5.2c+05
8.0c+04
9.7e+04
8.1e+04
1,8e+06
4.9c+05
Table 7: Operations Needed for six local solvers tested
15 graphs datasets.
= 10-6)
52

---

G
Related work
Many graph applications [2, 7, 29, 15, 27, 36, 30, 44, 25, 45, 51, 48] only require solving Equ. (1) approxi- mately. The reasons could be either the most energies of π are among small set of nodes forming small subgraphs, or one wants to study large graphs by checking them locally. Given graph G with nodes
0.0
LocCH CGM
-00
-1.
-1.5
-2.0
and m edges, there are two main types of iterative solvers for Equ. (I) as follows:
Standard iterative methods. Methods for solving linear systems have been well-established over the past decades (see textbooks of Saad [43], Golub Van Loan [19], Young [55]). The fastest linear solver for solving the symmetrized version of Equ. (1) the Conjugate Gradient Method (CGM) with runtime
0.00
0.25
Number of operations needed
0.50
0.75
1.00
1.25
1.00
1.75
Figure 15: Comparison of the error reduction between the proposed LocCH and the stan- dard CGM on the papers 100M dataset [23], in terms of the number of operations required. complexity Omm/v ((m where is the number of edges in the graph. It costs (m) to access the entire graph each iteration; hence, it much slower than local solvers, as demonstrated in Fig. 15. The symmetric diagonally dominant (SDD) solvers advance CGM further to have complexity Ô (in log" n (1)) [31, 46). Anikin et al. [4] considered the PageRank problem and proposed an algorithm with runtime depending on (n). This paper focuses on local algorithms where the goal is to avoid the dominant factor m by avoiding the full Qæ: operation.
Local algorithms. Local solvers, contrast standard counterparts, leverage the fact that the energy of π lives in small portion of the graph and hence do not require (mm) or (nz) per iteration. They are advantageous for huge-scale graphs demonstrated in Fig. 15. Andersen al. [2] used APPR to obtain an approximate of π for local clustering. Quite similar algorithms were developed in Berkhin [6] (bookmark-colol coloring algorithm) and Kloster & Gleich [28] Trauss-Southwedt 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 ((x)) using ISTA where the monotonicity and conservation properties remain. Hence, is nature to ask whether ((1((Vac) 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 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 guarantee locality, which does not exist in most accelerated frameworks; thus, developing an accelerated method that is guaranteed preserve intermediate variable sparsity remains challenging.
is worth mentioning that local methods are also closely related to sublinear time and local computa- tional al 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.
53

---

NeurIPS Paper Checklist
1. Claims
Question:
Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The abstract and introduction clearly state the main claims made by the paper, including its key contributions and scope.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contribu- tions made in the paper and important assumptions and limitations. A No NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: The accelerated algorithms, such as LocCH and LOCHB, may exhibit instability when a is small. This limitation arises due to the inherent constraints of global methods.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to vi- olations of these assumptions (e.g.. independence assumptions, noiseless settings, model well-speciffication asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low images are taken in low lighting. Or a speech-to-texx system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might be used by review- grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [Yes]
54

---

Justification:
All assumptions are clearly stated.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross-referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proof' can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition.
Inversely, any informal proof provided in the of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main experi- mental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: We have provided our code for the review process and will make publicly available upon publication.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, may be necessary to either make it possible for others replicate the model with the same dataset, provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., g., in the case a large language model), releasing of a model checkpoint, other means that appropriate the research performed.
While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution primarily a new algorithm, the paper should make clear how to reproduce that algorithm.
(b) If the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution be a way to access this model for reproducing the results or a way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
new model (e.g., a large language model), then there should either
(d) We recognize that reproducibility may be tricky in some cases, in which case authors welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification:
All datasets used in this study are publicly available.
55

---

Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines gui guides/ quidee/oodeeee Seta//ooeSuubisseeeeee ssiooilccy for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (c.g., for a new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https:// ssi cessionooligy cy) for more details.
nips. httpii/c/giidee/eeeeeeee
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted.
Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Answer: [Yes]
Justification:
All parameter settings our methods and baselines are included.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them.
TTee full details can be provided either with the code, in appendix, or supplemental material.
7. Experiment Statistical Significance
Question: Does the paper report eπror bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: For most of our results, we report the standard error over 50 random sampling nodes. Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
It is OK report I-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show tables or figures symmetric error bars that would yield results that out of range (c.g, negative error rates).
56

---

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.
8. Experiments Compute Resourees
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 experi- ments?
Answer: [Yes]
Justification: For our experiment, we used a server powered by an Intel(R) Xcon(R) Gold 5218R CPU, which features 40 cores (80 threads). The system 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 experi- mental 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 failed experiments that didn't 't make into the paper).
Code Of Ethies
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code Ethics https://neurips cc/public/Etiiicstuudll.....
Answer: [Yes]
Justification: The research conducted 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 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
preserve anonymity (e.g., if there 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 signifi- cantly 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 unintended uses (e.g.. disinformation, generating fake profiles, surveillance). faimeess considerations (c.g.. deploy- ment 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 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.
57

---

The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology heing used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of MI.).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, scraped datasets)?
Question:
Answer: [NA]
Justification: None.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions access the model implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do require this, but we encourage authors to take this into account and make a best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [NA]
Justification: None.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include a URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of derived asset (if it has changed) should be provided.
If this information is not available online, the authors are encouraged reach out to the asset's creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: None.
Guidelines:
The answer NA means that the paper does not release new assets.
58

---

Researchers should communicate the details of the dataset/ccoeemmodl 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.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and Research with Human Subjects
For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, applicable, as well as details about compensation (if any)?
Question:
Answer: [NA]
Justification: No human subjects involved.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Ethics, workers involved in data collection, curation, other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Sub- jects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/revieww based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
No human subjects involved.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
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.
59


Reviews:
[{"review_id": "HQkYPlroli", "reviewer": "Reviewer_cyZp", "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.", "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": "NA", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper addresses an important problem, provides deeper insights into an existing algorithm, provides a new algorithm and also reanalyzes the algorithm in a more fine-grained way. All of this is done via connection to GSSOR which seems to be new.\n\nI find the result quite interesting. However, I am not very familiar with recent work on personalized page rank.  For this reason, I recommend accepting but with a low confidence.", "weaknesses": "NA", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for taking the time and effort to review our paper carefully. Your positive perspective on our work is inspiring. The main concern on the assumption we made is addressed as follows:\n\n**Q1.** Comments on the assumptions we made in analyzing the local Chebyshev method.\n\n**A:** (Ignore this answer if you read the general response section). Let us recall that our Theorem C.6 (Standard CH) provides a runtime bound for the standard Chebyshev method. The total number of operations is bounded by\n\n$$\n\\mathcal{T}_{\\mathrm{CH}} \\leq m \\cdot \\left\\lceil\\frac{1+\\sqrt{\\alpha}}{2 \\sqrt{\\alpha}} \\ln \\frac{2}{\\epsilon}\\right\\rceil,\n$$\n\nwhere $m$ is the number of edges. Our key lemma, Lemma 4.1 (Line 246), captures the evolving process of the local Chebyshev method. By considering the geometric mean reduction on $\\beta_t$ and using Lemma 4.1, LocCH has the following runtime bound:\n\n$$\n\\mathcal{T} _{\\text{LocCH}} \\leq \\overline{\\operatorname{vol}}({\\mathcal S } _T) \\cdot \\left\\lceil\\frac{1+\\sqrt{\\alpha}}{(2-c) \\sqrt{\\alpha}} \\ln \\frac{2 y_T}{\\epsilon} \\right\\rceil\n$$\n\n\nwhere $\\overline{\\operatorname{vol}} ({\\mathcal S } _T)$ is the averaged volume and $y_T$ is defined by $y _ {t+1}-2 y _t+y _{t-1} / ( (1+\\beta _{t-1}) (1+\\beta _t ) )=0$. Ideally, when $\\epsilon \\rightarrow 0$, it roughly indicates the local solver is getting closer to the global one. This leads to $y_T \\rightarrow 1$, $\\overline{\\operatorname{vol}} ({\\mathcal S } _T) \\rightarrow m$, and $c \\rightarrow 0$ where $c$ is the parameter defined via the geometric mean of $\\beta_t$. However, verifying how strong or weak this assumption is can be difficult. The main reason is the convergence of the second-order difference equation itself is complicated and makes the analysis harder. We are investigating and exploring alternative directions now. A more reasonable approach might be to consider a typical example on a chain graph, where one can obtain a much simpler formulation for $y_t$. We will develop a refined analysis for specific graph types in future work.\n\nWe are happy to have further discussions if needed.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "dbM68Snnqf", "reviewer": "Reviewer_oXbQ", "summary": "This paper uses the evolving set procedure to give a local PageRank algorithm whose dependence on \\alpha (the reset probability) is \\sqrt{\\alpha}.\n\nIt 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.\n\nThe 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.", "questions": "Is the dependence on \\epsilon optimal? Aka. have methods with \\sqrt{\\eps} (or even \\log(1 / \\eps)) dependences been ruled out?", "limitations": "yes, limitations have been addressed, and are entirely theoretical w.r.t. some graph parameter regimes.", "rating": 9, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "Local algorithms are widely used in graph analytics. The question studied is natural, and has been proposed before.\n\nThe method is theoretically well-founded, and has significant technical depth.\n\nThe experiments are thorough and well documents, and clearly demonstrate the advantages of this method in multiple parameter regimes.", "weaknesses": "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.\n\nIdeally 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for taking the time and effort to review our paper carefully. Your positive perspective on our work is so inspiring. We also believe our work is novel and some new interesting problems are worth to explore. Your main concerns and our responses are as follows:\n\n---\n\n**Q1.** The gains only kick in at a relatively large number of steps: it's unclear that these are the parameter regimes in which local algorithms get used.\n\n**A:** Thank you for this excellent point. We assume you mean in what range of $\\epsilon$ is useful in practice. Let's take typical examples to see when $\\epsilon$ is enough/useful in practice. For the local clustering algorithm, in Fountoulakis's paper (see Fountoulakis' work in [1]), they used $\\epsilon = 10^{-5}$ where the number of nodes in graphs is in the range $[2\\times 10^6, 3\\times 10^6]$. This is around $\\epsilon \\approx 1/n$ (roughly corresponding to the largest speedup as shown in Figure 4, dashed vertical lines are $\\epsilon = 1/n$). For learning GNN models, PPRs are used to train the PPRGo model (see Bojchevski's work in [2]); the parameter $\\epsilon$ used is $10^{-4}$ where the number of nodes in graphs is in the range $[1.87\\times 10^4, 1.05 \\times 10^7]$. This means $\\epsilon \\leq 1/n$ in all graphs. We found similar parameter settings in graph embedding and online graph node classification. The effective range of $\\epsilon$ is $\\epsilon \\leq 10^{-4} / n$ where the speedup is significant, as presented in Fig. 4. This clearly indicates that our design demonstrates significant speedup, especially around $\\epsilon=1 / n$.\n\nIf you mean the number of local iterations needed, then it is true that the number of local iterations will be slightly larger than that of standard solvers. The key point is the combination of both the number of iterations and local volumes, and then the runtime of local solvers is much smaller than that of standard solvers. Whether there is an optimal tradeoff between a local number of iterations and local volumes is an interesting problem.\n\n---\n\n**Q2.** 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.\n\n**A:** Yes, the downstream tasks, such as training good GNN models that have decent F1 scores, are typical examples of the methods that our proposed methods could apply. Any promising results will make our paper even stronger.  However, due to space limits, we will consider these tasks in our future work. And we are looking forward to downstream applications.\n\n---\n\n**Q3.** Is the dependence on $\\epsilon$ or $\\alpha$ optimal?\n\n**A:** Let us recall that the runtime bound is optimal for the standard Chebyshev method. Specifically, the *first-order optimal methods (such as Chebyshev) for linear system* have the runtime bound $\\mathcal O (m/\\sqrt{\\alpha} \\cdot \\log 1/\\epsilon)$. This runtime bound is optimal for finding a solution of a sparse linear system defined on a chain graph (see page Theorem 3.15 of Bubeck's work in [1] for more details). We will discuss both parameters $\\alpha$ and $\\epsilon$.\n\n---\n**References**\n\n- [1] Fountoulakis, K., Roosta-Khorasani, F., Shun, J., Cheng, X., \\& Mahoney, M. W. (2019). Variational perspective on local graph clustering. Mathematical Programming, 174, 553-573.\n\n- [2] Bojchevski, Aleksandar, Johannes Gasteiger, Bryan Perozzi, Amol Kapoor, Martin Blais, Benedek Rózemberczki, Michal Lukasik, and Stephan Günnemann. \"Scaling graph neural networks with approximate PageRank.\" In Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery \\& Data Mining, pp. 2464-2473. 2020.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Q5X1JcsEJ3", "reviewer": "Reviewer_N8Gt", "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.", "questions": "One question is that what is the core reason for not obtaining convergence result for accelerated methods for all $\\epsilon>0$.", "limitations": "Authors have addressed limitations.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The main strengths of this paper are to develop a new algorithmic framework that can not only encompass existing algorithms but lead to the development of better ones that overcome previously known limitations for designing local graph clustering algorithms. They also back their theoretical analysis with the practical implementation of their method which is also shown to be superior to previous algorithms.", "weaknesses": "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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for taking the time and effort to review our paper carefully. We appreciate the positive perspective on our work. Your concern on the range of $\\epsilon$ can be effectively addressed as follows:\n\n**Q1.** What is the core reason for not obtaining convergence results for accelerated methods for all $\\epsilon > 0$?\n\n**A:** To clarify, our results cover all cases for   $0 \\leq \\epsilon \\leq 1/d_s$, including Theorems 3.3, 3.5, and 4.2. When $\\epsilon > 1/d_s$, the algorithm terminates before taking a single step, e.g. $T = 0$. In other words, we assume that the precision parameter satisfies $\\epsilon \\leq 1/d_s$ in our theorems because $\\epsilon \\leq 1/d_s$ ensures all local solvers run at least one local iteration, so that $T \\geq 1$. All theorems hold when $1/d_s < \\epsilon$. We thank the reviewer for pointing out if this was unclear in the manuscript and will revise it. \n\nWe are happy to have further discussions if needed!", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
**General Responses**

We thank all reviewers for their time and effort in carefully reading our paper. We are very happy that you like our work. Some general concerns are worth discussing as follows:

---

**Q1.** Comments on the assumption of the local Chebyshev (LocCH) method.

**A:** Let us recall that our Theorem C.6 (Standard CH) provides a runtime bound for the standard Chebyshev method. The total number of operations is bounded by

$$
\mathcal{T}_{\mathrm{CH}} \leq m \cdot \left\lceil\frac{1+\sqrt{\alpha}}{2 \sqrt{\alpha}} \ln \frac{2}{\epsilon}\right\rceil,
$$

where $m$ is the number of edges. Our key lemma, Lemma 4.1 (Line 246), captures the evolving process of the local Chebyshev method. By considering the geometric mean reduction on $\beta_t$ and using Lemma 4.1, LocCH has the following runtime bound:

$$
\mathcal{T} _{\text{LocCH}} \leq \overline{\operatorname{vol}}({\mathcal S } _T) \cdot \left\lceil\frac{1+\sqrt{\alpha}}{(2-c) \sqrt{\alpha}} \ln \frac{2 y_T}{\epsilon} \right\rceil
$$


where $\overline{\operatorname{vol}} ({\mathcal S } _T)$ is the averaged volume and $y_T$ is defined by $y _ {t+1}-2 y _t+y _{t-1} / ( (1+\beta _{t-1}) (1+\beta _t ) )=0$. Ideally, when $\epsilon \rightarrow 0$, it roughly indicates the local solver is getting closer to the global one. This leads to $y_T \rightarrow 1$, $\overline{\operatorname{vol}} ({\mathcal S } _T) \rightarrow m$, and $c \rightarrow 0$ where $c$ is the parameter defined via the geometric mean of $\beta_t$. However, verifying how strong or weak this assumption is can be difficult. The main reason is the convergence of the second-order difference equation itself is complicated and makes the analysis harder. We are investigating and exploring alternative directions now. A more reasonable approach might be to consider a typical example on a chain graph, where one can obtain a much simpler formulation for $y_t$. We will develop a refined analysis for specific graph types in future work.

---

**Q2.** Comments on the parameter regimes in which local algorithms get used.


**A:** The range of $\epsilon$ that is useful in real-world applications depends on the specific task. However, as far as we know, the speedup over standard solvers can be applied to a wide range of downstream tasks. Let's look at typical examples to see when $\epsilon$ is sufficient or useful in practice. For the local clustering algorithm, in Fountoulakis's paper (see Fountoulakis's work in [1]), they used $\epsilon = 10^{-5}$ where the number of nodes in graphs is in the range $[2\times 10^6, 3\times 10^6]$. This is around $\epsilon \approx 1/n$ (roughly corresponding to the largest speedup shown in Figure 4, where the dashed vertical lines are $\epsilon = 1/n$). For learning GNN models, PPRs are used to train the PPRGo model (see Bojchevski's work in [2]); the parameter $\epsilon$ used is $10^{-4}$ where the number of nodes in graphs is in the range $[1.87\times 10^4, 1.05 \times 10^7]$. This means $\epsilon \leq 1/n$ in all graphs. We found similar parameter settings in graph embedding and online graph node classification. 

Overall, the effective range of $\epsilon$ is $\epsilon \leq 10^{-4} / n$, where the speedup is significant, as presented in Fig. 4. This should cover most interesting downstream tasks where PPR is a crucial tool.

---

**Q3.** Is the dependence on $\epsilon$ optimal?

**A:** Let us recall that the runtime bound is optimal for the standard Chebyshev method. Specifically, the *first-order optimal methods (such as Chebyshev) for linear systems* have the runtime bound $\mathcal{O}(m/\sqrt{\alpha} \cdot \log 1/\epsilon)$. This runtime bound is optimal for finding a solution of a sparse linear system defined on a chain graph (see Theorem 3.15 of Bubeck's work in [3] for more details).

We believe this is true in terms of the **optimal local first-order methods** for local solvers. However, it remains interesting to see whether the lower bound can be identified, where one may find that it matches the bound conjectured in Line 295 of our manuscript.


---
**References**

- [1] Fountoulakis, K., Roosta-Khorasani, F., Shun, J., Cheng, X., & Mahoney, M. W. (2019). Variational perspective on local graph clustering. Mathematical Programming, 174, 553-573.

- [2] Bojchevski, Aleksandar, Johannes Gasteiger, Bryan Perozzi, Amol Kapoor, Martin Blais, Benedek Rózemberczki, Michal Lukasik, and Stephan Günnemann. "Scaling graph neural networks with approximate pagerank." In Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery \& Data Mining, pp. 2464-2473. 2020.

- [3] Sébastien Bubeck, Convex Optimization: Algorithms and
Complexity, https://arxiv.org/pdf/1405.4980.

Meta-review: 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.
Decision: Accept (poster)