Faster Local Solvers for Graph Diffusion Equations

Authors: jiahe  bai
baojian  zhou
deqing  yang
yanghua  xiao
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 08:01:02

Abstract: Efficient computation of graph diffusion equations (GDEs), such as Personalized PageRank, Katz centrality, and the Heat kernel, is crucial for clustering, training neural networks, and many other graph-related problems. Standard iterative methods require accessing the whole graph per iteration, making them time-consuming for large-scale graphs. While existing local solvers approximate diffusion vectors through heuristic local updates, they often operate sequentially and are typically designed for specific diffusion types, limiting their applicability. Given that diffusion vectors are highly localizable, as measured by the participation ratio, this paper introduces a novel framework for approximately solving GDEs using a local diffusion process. This framework reveals the suboptimality of existing local solvers. Furthermore, our approach effectively localizes standard iterative solvers by designing simple and provably sublinear time algorithms. These new local solvers are highly parallelizable, making them well-suited for implementation on GPUs. We demonstrate the effectiveness of our framework in quickly obtaining approximate diffusion vectors, achieving up to a hundred-fold speed improvement, and its applicability to large-scale dynamic graphs. Our framework could also facilitate more efficient local message-passing mechanisms for GNNs.

Main paper: 
 introduction: 
 Graph diffusion equations (GDEs), such as Personalized PageRank (PPR) [41, 44, 62] , Katz centrality (Katz) [46] , Heat kernel (HK) [21] , and Inverse PageRank (IPR) [51] , are fundamental tools for modeling graph data. These score vectors for nodes in a graph capture various aspects of their importance or influence. They have been successfully applied to many graph learning tasks including local clustering [2, 81] , detecting communities [49] , semi-supervised learning [78, 80] , node embeddings [63, 67] , training graph neural networks (GNNs) [10, 16, 19, 31, 75] , and many other applications [32] . Specifically, given a propagation matrix M associated with an undirected graph G(V, E), a general graph diffusion equation is defined asEQUATION)where f is the diffusion vector computed from a source vector s, and the sequence of coefficients c k satisfies c k ≥ 0. Equation (1) represents a system of linear, constant-coefficient ordinary differential equation, ẋ(t) = M x(t), with an initial condition x(0) and t ≥ 0. Standard solvers [34, 57] for computing f require access to matrix-vector product operations M x, typically involving O(m) operations, where m is the total number of edges in G. This could be time-consuming when dealing with large-scale graphs [42] .A key property of f is the high localization of its entry magnitudes, which reside in a small portion of G. Figure 1 demonstrates this localization property of f on PPR, Katz, and HK. Leveraging this locality property allows for more efficient approximation using local iterative solvers, which heuristically ) of example diffusion vectors f over 18 graphs, ordered from small (cora) to large (ogbn-papers100M). The ratio p(f ) is normalized by the number of nodes n.avoid O(m) operations. Local push-based methods [2, 7, 49] or their variants [4, 12] are known to be closely related to Gauss-Seidel [49, 18] . However, current local push-based methods are fundamentally sequential iterative solvers and focused on specific types such as PPR or HK. These disadvantages limit their applicability to modern GPU architectures and their generalizability.By leveraging the locality of f , we propose, for the first time, a general local iterative framework for solving GDEs using a local diffusion process. A novel component of our framework is to model local diffusion as a locally evolving set process inspired by its stochastic counterpart [58] . We use this framework to localize commonly used standard solvers, demonstrating faster local methods for approximating f . For example, the local gradient descent is simple, provably sublinear, and highly parallelizable. It can be accelerated further by using local momentum. Our contributions are -By demonstrating that popular diffusion vectors, such as PPR, Katz, and HK, have strong localization properties using the participation ratio, we propose a novel graph diffusion framework via a local diffusion process for efficiently approximating GDEs. This framework effectively tracks most energy during diffusion while maintaining local computation.-To demonstrate the power of our proposed framework, we prove that APPR [2] , a widely used local push-based algorithm, can be treated as a special case. We provide better diffusion-based bounds Θ(vol(S t )/(α • γ t )) where vol(S t )/γ t is a lower bound of 1/ϵ. This bound is effective for both PPR and Katz and is empirically smaller than Θ(1/(αϵ)), previously well-known for APPR.-We design simple and fast local methods based on standard gradient descent for APPR and Katz, which admit runtime bounds of min(vol(S t )/(α • γ t ), 1/(αϵ)) for both cases. These methods are GPU-friendly, and we demonstrate that this iterative solution is significantly faster on GPU architecture compared to APPR. When the propagation matrix M is symmetric, we show that this local method can be accelerated further using the local Chebyshev method for PPR and Katz.-Experimental results on GDE approximation of PPR, HK, and Katz demonstrate that these local solvers significantly accelerate their standard counterparts. We further show that they can be naturally adopted to approximate dynamic diffusion vectors. Our experiments indicate that these local methods for training are twice as fast as standard PPR-based GNNs. These results may suggest a novel local message-passing approach to train commonly used GNNs.All proofs, detailed experimental settings, and related works are postponed to the appendix. Our code is publicly available at https://github.com/JiaheBai/Faster-Local-Solver-for-GDEs. 
gdes, localization, and existing local solvers: 
 Table 1 : Example GDEs with their corresponding propagation matrix M , coefficients c k , and source s.Equ. M c k s PPR [62] AD -1 α(1α) k e s Katz [46] A α k e s HK [21] AD -1 e -τ τ k /k! e sAD -1θ k (θ k +θ 10 ) 2 e s APPNP [30] D -1 2 AD -1 2 α(1 -α) k xNotations. We consider an undirected graph G(V, E) where V = {1, 2, . . . , n} is the set of nodes, and E is the edge set with |E| = m. The degree matrix is D = diag(d 1 , d 2 , . . . , d n ), and A is the adjacency matrix of G. The volume of subset nodes S is vol(S) = v∈S d v . The set of neighbors of node v is denoted by N (v). The standard basis for node s is e s . The support of x is defined as supp(x) = {u : x u ̸ = 0}.Many graph learning tools can be represented as diffusion vectors. Table 1 presents examples widely used as graph learning tools. The coefficient c k usually exhibits exponential decay, so most of the energy of f is related to only the first few c k . We revisit these GDEs and existing local solvers.Personalized PageRank. Given the weight decaying strategy c k = α(1α) k and a predefined damping factor α ∈ (0, 1) with M = AD -1 , the analytic solution for PPR isEQUATIONThese vectors can be used to find a local graph cut or train GNNs [8, 10, 25, 79] . The well-known approximate PPR (APPR) method [2] locally updates the estimate-residual pair (x, r) as follows APPR :x ← x + αr u • e u , r ← r -r u • e u + (1 -α)r u • AD -1 • e u .Here, each active node u has a large residual r u ≥ ϵαd u with initial values x ← 0 and r ← αe s .The runtime bound for reaching r u < ϵαd u for all nodes is graph-independent and is Θ(1/(αϵ)), leveraging the monotonicity property. Previous analyses [49, 18] suggest that APPR is a localized version of the Gauss-Seidel (GS) iteration and thus has fundamental sequential limitations.Katz centrality. Given the weight c k = α k+1 with α ∈ (0, 1/∥A∥ 2 ) and a different propagation matrix M = A. The solution for Katz centrality can then be rewritten as:EQUATION)The goal is to approximate the Katz centrality ((I -αA) -1 -I)e s . Similar to APPR, a local solver [11] for Katz centrality, denoted as AKatz, can be designed as follows AKatz :x ← x + r u e u , r ← rr u e u + αr u Ae u .AKatz is a coordinate descent method [11] . Under the nonnegativity assumption of r (i.e., α ≤ 1/d max ), a convergence bound for the residual ∥r t+1 ∥ 1 is provided.Heat Kernel. The heat kernel (HK) [21] is useful for finding smaller but more precise local graph cuts. It uses a different weight decaying c k = τ k e -τ /k! and the vector is then defined asf HK = exp {-τ (I -M )} e s ,where τ is the temperature parameter of the HK equation and M = AD -1 . Unlike the previous two, this equation does not admit an analytic solution. Following the technique developed in [48, 49] ,given an initial vector s and temperature τ , the HK vector can be approximated using the first N terms of the Taylor polynomial approximation: definex N = N k=0 v k with v 0 = e s and v k+1 = M v k /k for k = 0, . . . , N . It is known that ∥f HK -x N ∥ 1 ≤ 1/(N !N ). This is equivalent to solving the linear system (I N +1 ⊗ I n -S N +1 ⊗ M ) v = e 1 ⊗ e s , where ⊗ denotes the Kronecker product. Then, the push-based method, namely AHK, has the following updatesAHK : v ← v + r u (e k ⊗ e u ) , r ← r -r u (I N +1 ⊗ I n -S N +1 ⊗ M ) (e k ⊗ e u ) .There are many other types of GDEs, such as Inverse PageRank [51] , which generalize PPR. Many GNN propagation layers, such as SGC [69] and APPNP [30] , can be formulated as GDEs. For all these f , we use the participation ratio [54] to measure its localization ability, defined asp(f ) = n i=1 |f i | 2 2 / n • n i=1 |f i | 4 .When the entries in f are uniformly distributed, such thatf i ∼ O(1/n), then p(f ) = O(1). However, in the extremely sparse case where f is a standard basis vector, p(f ) = 1/n indicates the sparsity effect. Figure 1 illustrates the participation ratios, showing that almost all vectors have ratios below 0.01. Additionally, larger graphs tend to have smaller participation ratios.APPR, AKatz, and AHK all have fundamental limitations due to their reliance on sequential Gauss-Seidel-style updates, which limit their parallelization ability. In the following sections, we will develop faster local methods using the local diffusion framework for solving PPR and Katz.(0) ∥ 1 ≥ • • • ∥r (t+ti) ∥ 1 ≥ ∥r (t+ti+1) ∥ 1 • • • .If the local diffusion process converges (i.e., S T = ∅), then T is bounded byT ≤ 1 ωγ T (1 -βP max ) ln ∥r (0) ∥ 1 ∥r (T ) ∥ 1 , where γ T ≜ 1 T T -1 t=0 γ t ≜ |St| i=1 r (t+ti) ui ∥r (t) ∥ 1 .Based on Theorem 3.2, we establish the following sublinear time bounds of LocalSOR for PPR.Theorem 3.3 (Sublinear runtime bound of LocalSOR for PPR). Let I T = supp(r (T ) ). Given an undirected graph G and a target source node s with α ∈ (0, 1), ω = 1, and provided 0 < ϵ ≤ 1/d s , the run time of LocalSOR in Equ. (6) for solving (I -(1α)AD -1 )f PPR = αe s with the stop condition ∥D -1 r (T ) ∥ ∞ ≤ αϵ and initials x (0) = 0 and r (0) = αe s is bounded as the followingEQUATION)where C PPR = 1/((1 -α)|I T |). The estimate x (T ) satisfies ∥D -1 (x (T ) -f PPR )∥ ∞ ≤ ϵ.I T = supp(r (T ) ) and C Katz = 1/((1 -α)|I T |).Given an undirected graph G and a target source node s with α ∈ (0, 1/d max ), ω = 1, and provided 0 < ϵ ≤ 1/d s , the run time of LocalSOR in Equ. (6) for solving (I -αA)x = e s with the stop condition ∥D -1 r (T ) ∥ ∞ ≤ ϵ and initials x (0) = 0 and r (0) = e s is bounded as the followingEQUATION)The estimate fKatz = x (T )e s satisfies ∥ fKatz - f Katz ∥ 2 ≤ ∥(I -αA) -1 D∥ 1 • ϵ. Accerlation when Q is Stieltjes matrix (ω * = 2/(1 + 1 -(α -1) 2 )).(0) ∥ 1 ≥ • • • ∥r (t+ti) ∥ 1 ≥ ∥r (t+ti+1) ∥ 1 • • • .If the local diffusion process converges (i.e., S T = ∅), then T is bounded byT ≤ 1 ωγ T (1 -βP max ) ln ∥r (0) ∥ 1 ∥r (T ) ∥ 1 , where γ T ≜ 1 T T -1 t=0 γ t ≜ |St| i=1 r (t+ti) ui ∥r (t) ∥ 1 .Proof. At each step t and t i = (i -1)/|S t | where i = 1, 2, . . . , |S t | , recall LocalSOR for solving Qx = s has the following updatesx (t+ti+1) = x (t+ti) + ω • r (t+ti) ui q uiui e ui , r (t+ti+1) = r (t+ti) - ω • r (t+ti) ui q uiui Q • e ui .Since we define Q as Q = I -βP , where P uu = 0 for all u ∈ V, and using q uiui = 1, we continue to have the updates of r as followsr (t+ti+1) = r (t+ti) - ω • r (t+ti) ui q uiui • Qe ui = r (t+ti) -ω • r (t+ti) ui e ui + ωβr (t+ti)(0) ∥ 1 ≥ • • • ∥r (t) ∥ 1 ≥ ∥r (t+1) ∥ 1 • • • .If the local diffusion process converges (i.e., S T = ∅), then T is bounded byT ≤ 1 γ T (1 -βP max ) ln ∥r (0) ∥ 1 ∥r (T ) ∥ 1 , where γ T ≜ 1 T T -1 t=0 γ t ≜ ∥r (t) St ∥ 1 ∥r (t) ∥ 1 .Proof. At each step t, recall LocalGD for solving Qx = (I -βP )x = s has the following updatesx (t+1) = x (t) + r (t) St , r (t+1) = r (t) -Qr (t)St . Since Q = I -βP , we continue to have the updates of r as the followingr (t+1) = r (t) -(I -βP )r (t) St = r (t) -r (t) St + βP r (t) St .By using induction, we show the nonnegativity of r (t) . First of all, r (0) ≥ 0 by our assumption. Let us assume r (t) ≥ 0. Then, the above equation gives r (t+1) ≥ 0 since β ≥ 0 and P ≥ 0 n×n by our assumption. The above equation can be written asr (t+1) + r (t) St = r (t) + βP r (t) StTaking ∥ • ∥ 1 on both sides, we have∥r (t+1) ∥ 1 + ∥r (t) St ∥ 1 = ∥r (t) ∥ 1 + β∥P r (t) St ∥ 1 .To make an effective reduction, β should be such that β∥P r(t) St ∥ 1 < ∥r (t)St ∥ 1 for all S t ⊆ V. For this, we then have the followingEQUATIONwhere note we definedγ t = ∥r (t) S t∥1 ∥r (t) ∥1 and assumed ∥P r(t) St ∥ 1 ≤ P max ∥r (t)St ∥ 1 . Apply ( 21) from t = 0 to T , it leads to the following bound∥r (T ) ∥ 1 ≤ T -1 t=0 (1 -γ t (1 -βP max )) ∥r (0) ∥ 1Note each of the term 1γ t (1 -βP max ) ≥ 0 during the updates, then γ t (1 -βP max ) ∈ (0, 1). To check the upper bound of T , we haveln ∥r (T ) ∥ 1 ∥r (0) ∥ 1 ≤ T -1 t=0 ln (1 -γ t (1 -βP max )) ≤ - T -1 t=0 γ t (1 -βP max ) ,which leads toT ≤ 1 γ T (1 -βP max ) ln ∥r (0) ∥ 1 ∥r (T ) ∥ 1 .Corollary 3.6 (Convergence of LocalGD for PPR and Katz). Let I T = supp(r (T ) ) and C = 1(1-α)|I T | . Use LocalGD to approximate PPR or Katz by using iterative procedure (10) . Denote T PPR and T Katz as the total number of operations needed by using LocalGD, they can then be bounded byEQUATION)for a stop condition ∥D -1 r (t) ∥ ∞ ≤ α PPR • ϵ. For solving KATZ, then the toal runtime is bounded byEQUATION)for a stop condition ∥D -1 r (t) ∥ ∞ ≤ ϵd u . The estimate equality is the same as of LocalSOR.Proof. We first show graph-independent bound O(1/(αϵ)) for PPR computation. Since r (t+1) = r (t) -Qr(t)St , then rearrange it to6 r (t+1) + r(t) St = r (t) + (1 -α)AD -1 r (t)St .Note that entries in r (t) are nonnegative. It leads to∥r (t+1) ∥ 1 + ∥r (t) St ∥ 1 = ∥r (t) ∥ 1 + ∥(1 -α)AD -1 r (t) St ∥ 1 , where note ∥(1 -α)AD -1 r (t) St ∥ 1 = (1 -α)∥r (t)St ∥ 1 , then we have ∥r(t) St ∥ 1 = (∥r (t) ∥ 1 -∥r (t+1) ∥ 1 )/α.At step t, LocalGD accesses the indices in S t . By the active node condition ϵαd ui ≤ r (t) ui , we havevol(S t ) = k i=1 d ui ≤ k i=1 r (t) ui ϵα = ∥r (t) St ∥ 1 ϵα = ∥r (t) ∥ 1 -∥r (t+1) ∥ 1 ϵα 2 .Then the total run time of LocalGD isT PPR := T -1 t=0 vol(S t ) ≤ 1 ϵα 2 T -1 t=0 ∥r (t) ∥ 1 -∥r (t+1) ∥ 1 = ∥r (0) ∥ 1 -∥r (T ) ∥ 1 ϵα 2 ≤ 1 ϵα ,where the last inequality is due to ∥r (t) ∥ 1 = α. Therefore, the total run time is O(1/(αϵ)). Follow the similar technique, we have sublinear runtime bound for Katz, i.e., T Kata ≤ 1 (1-αKatz•dmax)ϵ . Next, we show the local diffusion bound for PPR. Recall LocalGD has the initial x (0) = 0, r (0) = αe s and the following updatesx (t+1) = x (t) + r (t) St , r (t+1) = r (t) -Qr (t) St , S t = u : r (t)u ≥ ϵαd u , u ∈ V .Since β = 1α and P = AD -1 , by applying Theorem 3.5, we have the upper bound of T∥r (T ) ∥ 1 = T -1 t=0 (1 -αγ t ) ∥r (0) ∥ 1 ⇒ T ≤ 1 αγ T ln ∥r (0) ∥ 1 ∥r (T ) ∥ 1 .Note that each nonzero r (T ) u has at least part of the magnitude from the push operation of an active node. This means each nonzero of r (T ) satisfiesr (T ) u ≥ (1 -α)r ( t) v d v ≥ (1 -α)αϵd v d v = (1 -α)αϵ, for u ∈ I T ⇒ ∥r (T ) ∥ 1 ≥ (1 -α)αϵ|I T |,where t ≤ T . Hence, we haveT ≤ 1 αγ T ln ∥r (0) ∥ 1 ∥r (T ) ∥ 1 ≤ 1 αγ T ln ∥r (0) ∥ 1 αϵ(1 -α)|I T | := 1 αγ T ln C ϵ , where C = 1 (1-α)|I T | .To see the additive error, note that r (T ) := b -Qx (T ) = αe s -I -(1α)AD -1 x (T ) . It has the following equalityf s -x (T ) = I -(1 -α)AD -1 -1 r (T )To estimate the bound |f s [u]-x (T ) u |, we need to know ( I -(1α)AD -1 -1 r (T ) ) u . Specifically, the u-th entry of the above isf s [u] -x (T ) u = 1 α v∈V f v [u]r (T ) v By Proposition B.1, we know d u f u [v] = d v f v [u],then we continue to havef s [u] -x (T ) u = 1 α v∈V d u f u [v] d v r (T ) v ≤ 1 α v∈V d u f u [v]ϵα = ϵd u v∈V f u [v] = ϵd u ,where the first inequality is due to r (T ) v < ϵαd v , and the last equality follows fromv∈V f u [v] = 1.Combining all inequalities for u ∈ V, we have the estimate|D -1 ( f -f PPR )| 1 ≤ ϵ.We omit the proof of the sublinear runtime bound of LocalGD for Katz, as it largely follows the proof technique in Corollary 3.4. For the two lower bounds of 1/ϵ, i.e., vol(S T )/γ T ≤ 1 ϵ , it directly follows from the monotonicity and nonnegativity properties stated in Theorem 3.2 and Theorem 3.5.Remark B.3. The part of the theorem shares the similar proof strategy we provided in [77] . Here, we use a slightly different formulation of the linear system. 
faster solvers via local diffusion process: 
 We introduce a local diffusion process framework, which allows standard iterative solvers to be effectively localized. By incorporating this framework, we show that the computation of PPR and Katz defined in ( 2) and ( 3) can be locally approximated sequentially or in parallel on GPUs.To compute f PPR and f KATZ , it is equivalent to solving the linear systems, which can be written asEQUATIONwhere Q is a (symmetric) positive definite matrix with eigenvalues bounded by µ and L, i.e., µ ≤ λ(Q) ≤ L. For PPR, we solve(I -(1 -α)AD -1 )x = αe s , which is equivalent to solving (I -(1 -α)D -1/2 AD -1/2 )D -1/2 x = αD -1/2 e s .For Katz, we compute (I -αA)x = e s . The vector s is sparse, with supp(s) ≪ n. We define local diffusion process, a locally evolving set procedure inspired by its stochastic counterpart [58] as follows Definition 3.1 (Local diffusion process). Given an input graph G, a source distribution s, precision tolerance ϵ, and a local iterative method A with parameter θ for solving (4), the local diffusion process is defined as a process of updates {(x (t) , r (t) , S t )} 0≤t≤T . Specifically, it follows the dynamic systemEQUATION)where ϕ is a mapping via some iterative solver A θ . For all three diffusion processes (PPR, Katz, and HK), we set S 0 = {s}. We say this process converges when S T = ∅ if there exists such T ; the generated sequence of active nodes are S t . The total number of operations of the local solver A θ isT A θ = T -1 t=0 vol(S t ) = T • vol(S T ),where we denote the average of active volume as vol(ST ) = T -1 T -1 t=0 vol(S t ).Intuitively, we can treat this local diffusion process as a random walk on a Markov Chain defined on the space of all subsets of V, where the transition probability represents the operation cost. More efficient local solvers try to find a shorter random walk from S 0 to S T . The following two subsections demonstrate the power of this process in designing local methods for PPR and Katz.Given the estimate x (t) at t-th iteration, we define residual r (t) = s -Qx (t) and Q be the matrix induced by A and D. The typical local Gauss-Seidel with Successive Overrelaxation has the following online estimate-residual updates (See Section 11.2 of [34] ): at each time t = 0, 1, . . . , T -1, for each active node u i ∈ S t = {u 1 , u 2 , . . . , u |St| } with i = 1, 2, . . . , |S t |, we update each u i -th entry of x and corresponding residual r as the following LocalSOR :EQUATIONwhere t i ≜ (i -1)/|S t | is the current time, and update unit vector isẽ(t+ti) ui ≜ r (t+ti) uie ui /q uiui . Note that each update of ( 6) is Θ(d ui ). If S t = V, it reduces to the standard GS-SOR [34] . Therefore, APPR is a special case of (6), a local variant of GS-SOR with ω = 1. Figure 2 illustrates this procedure. APPR updates x at some entries and keeps track of large magnitudes of r per iteration, while LocalSOR allows for a better choice of ω; hence, the total number of operations is reduced. We establish the following fundamental property of LocalSOR. Theorem 3.2 (Properties of local diffusion process via LocalSOR). Let Q ≜ I -βP where P ≥ 0 n×n and P uv ̸ = 0 if (u, v) ∈ E; 0 otherwise. Define maximal value P max = max u∈V ∥P e u ∥ 1 . Assume that r (0) ≥ 0 is nonnegative and P max , β are such that βP max < 1, given the updates of (6), then the local diffusion process of ϕ x (t) , r (t) , S t ; s, ϵ, G, A θ = (LocalSOR, ω) with ω ∈ (0, 1) has the following properties 1. Nonnegativity. r (t+ti) ≥ 0 for all t ≥ 0 and t i = (i -1)/|S t | with i = 1, 2, . . . , |S t |. LocalSOR uses fewer local iterations, costs less total active volume, and obtains better approximate solutions. We choose the source node s = 0 with ϵ = 0.02 and α = 0.25.The fundamental limitation of LocalSOR is its reliance on essentially sequential online updates, which may be challenging to utilize with GPU-type computing resources. Interestingly, we developed a local iterative method that is embarrassingly simpler, highly parallelizable, and provably sublinear for approximating PPR and Katz. First, one can reformulate the equation of PPR and Katz as2 x * t = arg minEQUATIONA natural idea for solving the above is to use standard GD. Hence, following a similar idea of local updates, we simultaneously update solutions for S t . We propose the following local GD.LocalGD :EQUATIONTheorem 3.5 (Properties of local diffusion process via LocalGD). Let Q ≜ I -βP where P ≥ 0 n×n and P uv ̸ = 0 if (u, v) ∈ E; 0 otherwise. Define maximal valueP max = max St⊆V ∥P r St ∥ 1 /∥r St ∥ 1 .Assume that r (0) = s ≥ 0 and P max , β are such that βP max < 1, given the updates of (10), then the local diffusion process of ϕ S t , x (t) , r (t) ; G, A θ = (LocalGD, µ, L) has the following properties 1. Nonnegativity. r (t) ≥ 0 for all t ≥ 0.2. Monotonicity property.∥r (0) ∥ 1 ≥ • • • ∥r (t) ∥ 1 ≥ ∥r (t+1) ∥ 1 • • • .If the local diffusion process converges (i.e., S T = ∅), then T is bounded byT ≤ 1 γ T (1 -βP max ) ln ∥r (0) ∥ 1 ∥r (T ) ∥ 1 , where γ T ≜ 1 T T -1 t=0 γ t ≜ ∥r (t) St ∥ 1 ∥r (t) ∥ 1 .Indeed, the above local updates have properties that are quite similar to those of SOR. Based on Theorem 3.5, we establish the sublinear runtime bounds of LocalGD for solving PPR and Katz.Corollary 3.6 (Convergence of LocalGD for PPR and Katz). Let I T = supp(r (T ) ) and C = 1 (1-α)|I T | . Use LocalGD to approximate PPR or Katz by using iterative procedure (10) . Denote T PPR and T Katz as the total number of operations needed by using LocalGD, they can then be bounded byEQUATIONfor a stop condition ∥D -1 r (t) ∥ ∞ ≤ α PPR • ϵ. For solving KATZ, then the toal runtime is bounded byEQUATION)for a stop condition ∥D -1 r (t) ∥ ∞ ≤ ϵd u . The estimate of equality is the same as that of LocalSOR.Remark 3.7. Note that LocalGD is quite different from iterative hard-thresholding methods [43] where the time complexity of the thresholding operator is O(n) at best, hence not a sublinear algorithm.Accerlerated local Chevbyshev. One can extend the Chebyshev method [39] , an optimal iterative solver for (9) . Following [24] , Chebyshev polynomials T n are defined via following recursionsEQUATIONGiven δ 1 = L-µ L+µ , x 1 = x 0 -2 L+µ ∇f (x 0 ), the standard Chevyshev method is defined asx k = x k-1 - 4δ k L -µ ∇f (x k-1 ) + 1 -2δ k L + µ L -µ (x k-2 -x k-1 ) , δ k = 1 2 L+µ L-µ -δ k-1, where µ ≤ λ(Q) ≤ L. For example, for approximating PPR, we propose the following local updatesLocalCH : π (t+1) = π (t) + 2δt+1 1-α D 1/2 r (t) St + δ t:t+1 (π (t) -π (t-1) ) St , δ t+1 = 2 1-α -δ t -1 D 1/2 r (t+1) = D 1/2 r (t) -(π (t+1) -π (t) ) + (1 -α)AD -1 (π (t+1) -π (t) ).The sublinear runtime analysis for LocalCH is complicated since it does not follow the monotonicity property during the updates. Whether it admits, an accelerated rate remains an open problem. 
applications to dynamic gdes and gnn propagation: 
 Our accelerated local solvers are ready for many applications. We demonstrate here that they can be incorporated to approximate dynamic GDEs and GNN propagation. We consider the discretetime dynamic graph [47] which contains a sequence of snapshots of the underlying graphs, i.e., G 0 , G 1 , . . . , G T . The transition from G t-1 to G t consists of a list of events O t involving edge deletions or insertions. The graph G t is then represented as:G 0 O1 --→ G 1 O2 --→ G 2 O3 --→ • • • G T -1 O T --→ G T .The goal is to calculate an approximate f t for the PPR linear system Q t f t = αe s at time t. That is,EQUATIONThe key advantage of updating the above equation (presented in Algorithm 3) is that the APPR algorithm is used as a primary component for updating, making it much cheaper than computing from scratch. Consequently, we can apply our local solver LocalSOR to the above equation, and we found that it significantly sped up dynamic GDE calculations and dynamic PPR-based GNN training. 
experiments: 
 Datasets. We conduct experiments on 18 graphs ranging from small-scale (cora) to large-scale (papers100M), mainly collected from Stanford SNAP [45] and OGB [42] (see details in Table 3 ). We focus on the following tasks: 1) approximating diffusion vectors f with fixed stop conditions using both CPU and GPU implementations; 2) approximating dynamic PPR using local methods and training dynamic GNN models based on InstantGNN models [75] . 3Baselines and Experimental Setups. We consider four methods with their localized counterparts: Gauss-Seidel (GS)/LocalGS, Successive Overrelaxation (SOR)/LocalSOR, Gradient Descent (GD)/LocalGD, and Chebyshev (CH)/LocalCH. 4 Specifically, we use the stop condition ∥D -1 r (t) ∥ ∞ ≤ ϵα for PPR and ∥D -1 r (t) ∥ ∞ ≤ ϵ for Katz, while we follow the parameter settings used in [49] for HK. For LocalSOR, we use the optimal parameter ω * suggested in Section 3.2. We randomly sample 50 nodes uniformly from lower to higher-degree nodes for all experiments. All results are averaged over these 50 nodes. We conduct experiments using Python 3.10 with CuPy and Numba on a server with 80 cores, 256GB of memory, and two 28GB NVIDIA-4090 GPUs.Local solvers are faster and use fewer operations. We first investigate the efficiency of the proposed local solvers for computing PPR and Katz centrality. We set (α PPR = 0.1, ϵ = 1/n) for PPR and (α Katz = 1/(∥A∥ 2 + 1), ϵ = 1/m) for Katz. We use a temperature of τ = 10 and ϵ = 1/ √ n for HK. All algorithms for HK estimate the number of iterations by considering the truncated error of Taylor approximation. Table 2 presents the speedup ratio of four local methods over their standard counterparts. For five graph datasets, the speedup is more than 100 times in terms of the number of operations in most cases. This strongly demonstrates the efficiency of these local solvers. Figure 4 presents all such results. Key observations are: 1) All local methods significantly speed up their global counterparts on all datasets. This strongly indicates that when ϵ is within a certain range, local solvers for GDEs are much cheaper than their global counterparts. 2) Among all local methods, LocalGS and LocalGD have the best overall performance. This may seem counterintuitive since LocalSOR and LocalCH are more efficient in convergence rate. However, the number of iterations needed for ϵ = 1/n or ϵ = 1/m is much smaller. Hence, their improvements are less significant. Graph Index # of Operations ( EQUATIONThe next critical question is when standard solvers can be effectively localized, meaning under which ϵ conditions we see speedup over standard solvers.We conduct experiments to determine when local GDE solvers show speedup over their global counterparts for approximating Katz, HK, and PPR on the wiki-talk graph dataset for different ϵ values, ranging from lower precision (2 -17 ) to high precision (2 -32 ). As illustrated in Figure 6 , when ϵ is in the lower range, local methods are much more efficient than the standard counterparts.As expected, the speedup decreases and becomes worse than the standard counterparts when high precision is needed. This indicates the broad applicability of our framework; the required precisions in many real-world graph applications are within a range where our framework is effective.We conducted experiments on the efficiency of the GPU implementation of LocalGD and GD, specifically using LocalGD's GPU implementation to compare with other methods. Figure 7 presents the running time of global and local solvers as a function of ϵ on the wiki-talk dataset. LocalGD is the fastest among a wide range of ϵ. This indicates that, when a GPU is available and ϵ is within the effective range, LocalGD (GPU) can be much faster than standard GD (GPU) and other methods based on CPUs. We observed similar patterns for computing Katz in Figure 11 . To test the applicability of the proposed local solvers, we apply LocalSOR to GNN propagation and test its efficiency on ogbn-arxiv and ogbn-products. We use the InstantGNN [75] model (essentially APPNP [30] ), following the experimental settings in [75] . Specifically, we set α = 0.1, ϵ = 10 -2 /n, and ω = ω * to the optimal value for our cases. For ogbnarxiv, we randomly partition the graph into 16 snapshots (each snapshot with 59,375 edges), where the initial graph G 0 contains 17.9% of the edges. With a similar performance on testing accuracy, as illustrated in Figure 8 , the Instant-GNN model with LocalSOR (Algo. 7) has a significantly shorter training time than its local propagation counterpart (Algo. 6). Figure 9 presents the accumulated operations over these 16 snapshots. The faster local solver is LocalSOR (Dynamic), which dynamically updates (x, r) for ( 14) according to Algo. 5 and Algo. 3, whereas LocalSOR (Static) updates all approximate PPR vectors from scratch at the start of each snapshot for (14) . We observed a similar pattern for LocalCH, where the number of operations is significantly reduced compared to static ones.2 -23 2 -31 2 -19 2 -27 
limitations and conclusion: 
 When ϵ is sufficiently small, the speedup is insignificant, and it is unknown whether more efficient local solvers can be designed under this setting. Although we observed acceleration in practice, the accelerated bounds for LocalSOR and LocalCH have not been proven. Another limitation of local solvers is that they inherit the limitations of their standard counterparts. This paper mainly develops local solvers for PPR and Katz, and it remains interesting to consider other types of GDEs.We propose the local diffusion process, a local iterative algorithm framework based on the locally evolving set process. Our framework is powerful in capturing existing local iterative solvers such as APPR for GDEs. We then demonstrate that standard iterative solvers can be effectively localized, achieving sublinear runtime complexity when monotonicity properties hold in these local solvers. Extensive experiments show that local solvers consistently speed up standard solvers by a hundredfold. We also show that these local solvers could help build faster GNN models [31, 13, 14, 20] . We expect many GNN models to benefit from these local solvers using a local message passing strategy, which we are actively investigating. Several open problems are worth exploring, such as whether these empirically accelerated local solvers admit accelerated sublinear runtime bounds without the monotonicity assumption. 
safeguards: 
 Question: 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, or scraped datasets)?Answer: [NA]Justification: The paper does not involve data or models with a high risk for misuse.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a 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 to usage guidelines or restrictions to access the model or 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 is challenging, and many papers do not 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 of 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: The paper does not involve licenses for existing assets.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 are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: Disclaimer: The local diffusion process introduced in this paper shares a similar conceptual foundation with our concurrent work [77] , where we propose the locally evolving set process. However, these two approaches address distinct problems. In this work, our focus is on investigating a different class of equations and exploring whether the local diffusion process can accelerate the training of Graph Neural Networks (GNNs).Graph diffusion equations (GDEs) and localization. Graph diffusion vectors can either be exactly represented as a linear system or be approximations of an underlying linear system. The general form defined in Equation ( 1) has appeared in various literature [49, 51] (see more examples in [68] ) or can be found in the survey works of [32] . The locality and localization properties of the diffusion vector f , including the PPR vector, have been explored in [64, 2, 59, 33, 60] . Inspired by methods for measuring the dynamic localization of quantum chaos [35, 22] , which are also used for analyzing graph eigenfunctions [11, 26, 54] , we use the inverse participation ratio to measure the localization ability of f . However, there is a lack of a local framework for solving GDEs efficiently.Standard solvers for GDEs. The computation of a graph diffusion equation approximates the exponential of a propagation matrix, either through an inherently Neumann series or as an analytic system of linear equations. Well-established iterative methods, such as those approximating via the Neumann series or Padé approximations, exist (see more details in the surveys [57, 6, 61] and [34] ). Each iteration requires computing the matrix-vector dot product, requiring O(m) operations for all the above methods. This is the computation we aim to avoid.Local methods for GDEs. The most well-known local solver for solving PPR was initially proposed in [2] for local clustering. PPR, as a fundamental tool for graph learning, has been applied to many core problems, including graph neural networks [8, 17, 69, 15] , graph diffusion [31, 74] , and its generalizations [51, 19] for GNNs or clustering [50, 71] . Many variants have been proposed for solving such equations in different forms using local methods that couple local methods with Monte Carlo sampling strategies, achieving sublinear time [5, 53] . The work of [55] attempts to accelerate PPR computation via the Gauss-Southwell method, but it does not work well in practice [48] . A notable finding from [70] is that the local computation of PPR is equivalent to power iteration when the precision is high. There are also works on temporal PPR vector computation [66, 36, 56, 52] , and our work can be immediately applied to these. Some local variants of APPR have been proposed for Katz and HK [11, 49] . However, these local solvers are sequential and lack the ability to be easily implemented on GPUs. Our local framework addresses this challenge.PPR on dynamic graphs. Ranking on dynamic graphs has many applications [65] , including motif clustering [71, 29] , embedding on large-scale graphs [63, 23, 37, 38, 76, 18] , and designing graph neural networks [75, 72] . The local computation via forward push (a.k.a. APPR algorithm [2] ) or reverse push (a.k.a. reverse APPR [1] ) is widely used in GNN propagation and bidirectional propagation algorithms [16, 28] . Our framework helps to design more efficient dynamic GNNs.We list all notations used in our proofs here.• e u : The indicator vector where the u-th entry is 1, and 0 otherwise. • N (u): Set of neighbors of node u.• d u : The degree of node u.• G(V, E): The underlying graph with the set of nodes V and set of edges E.• A: Adjacency matrix of G.• D: Degree matrix of G when G is undirected.• Π: Personalized PPR matrix, defined asΠ = α I -(1 -α)AD -1 -1 .• α PPR : Damping factor α PPR ∈ (0, 1) for the PPR equation. • α Katz : Parameter α Katz ∈ (0, 1) for the Katz equation.• ϵ: The error tolerance.• ∥D -1 r∥ ∞ ≤ α PPR ϵ: The stop condition for local solvers of PPR.• ∥D -1 r∥ ∞ ≤ ϵ: The stop condition for local solvers of Katz.In this section, we restate all theorems and present all missing proofs. We first restate the two target linear systems we aim to solve, as introduced in Section 3. Recall PPR is defined asI -(1 -α PPR )AD -1 f PPR = α PPR e s .It can be equivalently written as PPR:EQUATIONwhere α PPR ∈ (0, 1) and x * = αQ -1 D -1/2 e s . Hence, f PPR = D 1/2 x * . The Katz centrality can be written as Katz:EQUATIONwhere α Katz ∈ (0, 1/∥A∥ 2 ). Then f Katz = x *e s . 5Proposition B.1 ( [3, 53] ). Let G(V, E) be an undirected graph and u and v be two vertices in V. Denote f s := α(I -(1α)AD -1 )e s as the PPR vector of a source node s. Thend u • f u [v] = d v • f v [u].Proof. We directly follow the proof strategy in [53] for completeness. For path P = {s, v 1 , v 2 , . . . , v k , t} in G, we denote its length as ℓ(P ) (here ℓ(P ) = k + 1 ), and define its reverse path to be P = {t, v k , . . . , v 2 , v 1 , s}note that ℓ(P ) = ℓ( P ). Moreover, we know that a random walk starting from s traverses path P with probabilityP[P ] = 1 ds • 1 dv 1 • • • • • 1 dv k, and thus, it is easy to see that we haveP[P ] • d s = P[ P ] • d t .Now let P st denote the paths in G starting at s and terminating at t. Then, we can re-write f s [t] as:f s [t] = P ∈Pst α(1 -α) ℓ(P ) P[P ] = P ∈Pst α(1 -α) ℓ(P ) d t d s P[ P ] = d t d s P ∈Pts α(1 -α) ℓ( P ) P[ P ] = d t d s f t [s],Algo 1 APPR(G, ϵ, α, s) [2] adopted from [8] 1: x ← 0, r ← αe s 2: while ∃u, r u ≥ ϵαd u do 3:r, x ← PUSH(u, r, x) 4: Return x 1: PUSH(u, r, x):2:ru ← r u 3:x u ← x u + ru 4: r u ← 0 5: for v ∈ N (u) do 6: r v ← r v + (1-α)ru du 7: Return (r, x)We adopt the typical implementation of APPR as presented in Algo. 1.Theorem 3.2 (Properties of local diffusion process via LocalSOR). Let Q ≜ I -βP where P ≥ 0 n×n and P uv ̸ = 0 if (u, v) ∈ E; 0 otherwise. Define maximal value P max = max u∈V ∥P e u ∥ 1 . Assume that r (0) ≥ 0 is nonnegative and P max , β are such that βP max < 1, given the updates of (6), then the local diffusion process of ϕ x (t) , r (t) , S t ; s, ϵ, G, A θ = (LocalSOR, ω) with ω ∈ (0, 1) has the following properties 1. Nonnegativity. r (t+ti) ≥ 0 for all t ≥ 0 and t i = (i -1)/|S t | with i = 1, 2, . . . , |S t |.• P e ui .By using induction, we show the nonnegativity of r (t) . First of all, r (0) ≥ 0 by our assumption. Let us assume r (t+ti) ≥ 0, then the above equation gives r (t+ti+1) ≥ 0 since ω ≤ 1 and P ≥ 0 n×n by our assumption. Since we assume ω ≤ 1 and P ≥ 0 n×n , the above equation can be written as r (t+ti+1) + ω • r (t+ti) ui e ui = r (t+ti) + ωβr (t+ti) ui• P e ui . By taking ∥ • ∥ 1 on both sides, we have∥r (t+ti+1) ∥ 1 + ωr (t+ti) ui = ∥r (t+ti) ∥ 1 + ωβr (t+ti) ui ∥P e ui ∥ 1To make an effective reduction, β should be such that β∥P e ui ∥ 1 < 1 for all u i ∈ V. Summing over i = 1, 2, . . . , |S t | of the above equation, we then have the followingEQUATIONNote that we definedγ t = |S t | i=1 r (t+t i ) u i ∥r (t) ∥1and (1 -βP max) ≤ (1 -βP ui ) for all u i ∈ S t . The inequality of (17) leads to the following bound∥r (T ) ∥ 1 ≤ T -1 t=0 (1 -ωγ t (1 -βP max )) ∥r (0) ∥ 1 .To further simply the above upper bound, since each term 1ωγ t (1 -βP max ) ≥ 0 during the updates, it follows that ωγ t (1 -βP max ) ∈ (0, 1). If there exists T such that S T = ∅, then we can obtain an upper bound of T asln ∥r (T ) ∥ 1 ∥r (0) ∥ 1 ≤ T -1 t=0 ln (1 -ωγ t (1 -βP max )) ≤ - T -1 t=0 ωγ t (1 -βP max ) ,which leads toT ≤ 1 ωγ T (1 -βP max ) ln ∥r (0) ∥ 1 ∥r (T ) ∥ 1 .Theorem 3.3 (Sublinear runtime bound of LocalSOR for PPR). Let I T = supp(r (T ) ). Given an undirected graph G and a target source node s with α ∈ (0, 1), ω = 1, and provided 0 < ϵ ≤ 1/d s , the run time of LocalSOR in Equ. (6) for solving (I -(1α)AD -1 )f PPR = αe s with the stop condition ∥D -1 r (T ) ∥ ∞ ≤ αϵ and initials x (0) = 0 and r (0) = αe s is bounded as the followingEQUATION)whereC PPR = 1/((1 -α)|I T |). The estimate x (T ) satisfies ∥D -1 (x (T ) -f PPR )∥ ∞ ≤ ϵ.Proof. Recall S t = {u 1 , u 2 , . . . , u |St| } be the set of active nodes processed in t iteration. For convenience, we denote |S t | = k. By LocalSOR defined in ( 6) for solvingI -(1 -α)AD -1 Q f PPR = αe s s .We have the following online iteration updates for each active node u i ∈ S t (recall q uiui = 1).x(t+ i k ) = x (t+ i-1 k ) + ω • r (t+ i-1 k ) ui • e ui , for i = 1, . . . , k and u i ∈ S t , r (t+ i k ) = r (t+ i-1 k ) -ω • r (t+ i-1 k ) ui • e ui + ω • (1 -α) • r (t+ i-1 k ) ui • AD -1 e ui .Note for each active node u i , we have r(t+ i-1 k ) ui≥ ϵαd ui . The total operations for LOCSOR isT LocalSOR := T -1 t=0 vol(S t ) = T -1 t=0 k i=1 d ui ≤ T -1 t=0 k i=1 r (t+(i-1)/k)) ui ϵα = T -1 t=0 k i=1 ∥r (t+ i-1 k ) ∥ 1 -∥r (t+ i k ) ∥ 1 ωϵα 2 = ∥r (0) ∥ 1 -∥r (T ) ∥ 1 ωϵα 2 ≤ ∥r (0) ∥ 1 ωϵα 2 = 1 ϵα ,where the last equality follows from ω = 1 and ∥r (0) ∥ 1 = α. Next, we focus on the proof of our newly derived bound Õ(vol(S T )/(αγT )). To check the upper bound of T , from Theorem 3.2 by noticing β = (1α) and P max = 1, this leads to the followingT ≤ 1 ωγ T (1 -βP max ) ln ∥r (0) ∥ 1 ∥r (T ) ∥ 1 = 1 αγ T ln ∥r (0) ∥ 1 ∥r (T ) ∥ 1 .Next, we try to give a lower bound of r (T ) . Note that after the last iteration T , for each nonzero residual r (T ) u ̸ = 0, u ∈ I T , there is at least one possible update that happened at node u: 1) Node u has a neighbor v u ∈ N (u), which was active. This neighbor v u pushed some residual (1α)r(t ′ )vu /d vu to u where the time t ′ < T . Hence, for all u ∈ I T , we have∥r (T ) ∥ 1 = u∈I T r (T ) u ≥ u∈I T (1 -α)r (t ′ ) vu d vu ≥ u∈I T (1 -α)ϵαd vu d vu = u∈I T (1 -α)ϵα ≥ αϵ(1 -α)|I T |,where the first equality is due to the nonnegativity of r guaranteeing by Theorem 3.2 and the second inequality is due to the fact that r(t ′ ) uwas active residuals before the push operation. Applying the above lower bound of ∥r (T ) ∥ 1 , we obtain∥r (0) ∥ 1 ∥r (T ) ∥ 1 ≤ ∥r (0) ∥ 1 αϵ(1 -α)|I T | = 1 ϵ(1 -α)|I T | := C PPR ϵ ,where C PPR = 1/((1α)|I T |) and I T = supp(r (T ) ). Combine the above inequality and the upper bound T , we obtain our new local diffusion-based bounds. The rest is to show a lower bound of 1/ϵ. By using the active node condition, for u i ∈ S t , we haveEQUATIONWe continue to havevol(S T ) γ T = T -1 t=0 vol (S t ) T -1 t=0 γ t ≤ T -1 t=0 vol (S t ) T -1 t=0 |S t | i=1 r (t+t i ) u i ∥r (0) ∥1 = α T -1 t=0 vol (S t ) T -1 t=0 |St| i=1 r (t+ti) ui ≤ 1 ϵ ,where the first inequality follows from the fact that|S t | i=1 r (t+t i ) u i ∥r (0) ∥1 ≤ i=1 |S t | r (t+t i ) u i ∥r (t) ∥1for t = 0, 1, . . . , T -1, due to the monotonicity property of ∥r (t) ∥ 1 . The last inequality is from (19) . Combining these, we finish the proof of the sublinear runtime bound. The rest is to prove the estimation quality. Note ∥D -1 r (T ) ∥ ∞ ≤ αϵ directly implies ∥D -1 (x (T )f PPR )∥ ∞ ≤ ϵ by using Proposition B.1 and the corresponding stop condition.Remark B.2. The above theorem shares the proof strategy we provided in [77] . Here, we use a slightly different formulation of the linear system. Corollary 3.4 (Runtime bound of LocalSOR for Katz). Let I T = supp(r (T ) ) and C Katz = 1/((1α)|I T |). Given an undirected graph G and a target source node s with α ∈ (0, 1/d max ), ω = 1, and provided 0 < ϵ ≤ 1/d s , the run time of LocalSOR in Equ. (6) for solving (I -αA)x = e s with the stop condition ∥D -1 r (T ) ∥ ∞ ≤ ϵ and initials x (0) = 0 and r (0) = e s is bounded as the followingEQUATION)The estimate fKatz = x (T )e s satisfies ∥ fKatzf Katz ∥ 2 ≤ ∥(I -αA) -1 D∥ 1 ϵ.Proof. The linear system (I -αA)x = s using LocalSOR updates can be written asx (t+ti+∆t) = x (t+ti) + ωr (t+ti) uiq uiui e ui , r (t+ti+∆t) = r (t+ti) -ωr (t+ti) uiq uiui e ui + ωr (t+ti) ui αA e ui q uiui , where q uiui = 1. The updates of r (t+ti) can be simplified as the following r (t+ti+∆t) = r (t+ti)ωr (t+ti) ui e ui + ωr (t+ti) ui αAe ui .The proof of nonnegativity of r (t) follows the similar induction as in previous theorem by noticing the spectral radius of A is |λ(A)| ≤ d max . Hence, we haver (t+1) = r (t) -ω |St| i=1 r (t+ti) ui e ui + ω |St| i=1 r (t+ti) ui αAe ui ,Move the negative term to left and take ℓ 1 -norm on both sides, we haveω |St| i=1 r (t+ti) ui = ∥r (t) ∥ 1 -∥r (t+1) ∥ 1 + ω |St| i=1 r (t+ti) ui αd ui ∥r (t+1) ∥ 1 = 1 - ω |St| i=1 r (t+ti) ui -ω |St| i=1 r (t+ti) ui αd ui ∥r (t) ∥ 1 ∥r (t) ∥ 1 ≤ 1 - ω(1 -αd max ) |St| i=1 r (t+ti) ui ∥r (t) ∥ 1 ∥r (t) ∥ 1 .It leads to the followingω |St| i=1 r (t+ti) ui -ω |St| i=1 r (t+ti) ui αd ui = ∥r (t) ∥ 1 -∥r (t+1) ∥ 1 Since each r (t+ti) ui ≥ ϵd ui , vol(S t ) ≤ |St| i=1 r (t+ti) ui/ϵ, we continue haveω |St| i=1 r (t+ti) ui -ω |St| i=1 r (t+ti) ui αd ui = ω |St| i=1 r (t+ti) ui (1 -αd ui ) ≥ ω |St| i=1 r (t+ti) ui (1 -αd max ) .Finally, the above inequality gives the following upper runtime bound asT -1 t=0 vol(S t ) ≤ T -1 t=0 ∥r (t) ∥ 1 -∥r (t+1) ∥ 1 ωϵ(1 -αd max ) = ∥r (0) ∥ 1 -∥r (T ) ∥ 1 ωϵ(1 -αd max ) ≤ 1 ωϵ(1 -αd max ).By letting β = α and P = A with P max = d max , we apply Theorem 3.2, to obtain the local diffusion-based bound as∥r (t+1) ∥ 1 ≤ (1 -ω(1 -αd max )β t ) ∥r (t) ∥ 1 ∥r (T ) ∥ 1 ≤ T -1 t=0 (1 -ω(1 -αd max )β t ) ∥r (0) ∥ 1 = T -1 t=0 (1 -ω(1 -αd max )β t ) .Using the similar technique provided in the previous case, and letting ω = 1 we can haveT Katz ≤ vol(S T ) (1 -αd max )γ T ln C Katz ϵ .To check the estimate equality, we haver (T ) = e s -(I -αA)x (T ) =⇒ (I -αA) -1 r (T ) = (I -αA) -1 -I e s + e s -x (T ) = f Katz + e s -x (T ) .Let fKatz := x (T )e s This leads to Assume that r (0) = s ≥ 0 and P max , β are such that βP max < 1, given the updates of (10), then the local diffusion process of ϕ S t , x (t) , r (t) ; G, A θ = (LocalGD, µ, L) has the following properties 1. Nonnegativity. r (t) ≥ 0 for all t ≥ 0.∥ fKatz -f Katz ∥ 2 = ∥(I -αA) -1 DD -1 r (T ) ∥ 2 ≤ ∥(I -αA) -1 D∥ 1 • ∥D -1 r (T ) ∥ ∞ ≤ ∥(I -αA) -1 D∥ 1 ϵ.• Chebyshev for PPR.Let Q = I -(1 -α)D -1/2 AD -1/2 and b = αD -1/2 e s . Then f (t) = D 1/2 x (t) . The eigenvalue of Q is in range [α, 2 -α]. So, we let L = 2 -α and µ = α. 1. δ 1 = 1 -α, x (0) = 0, D 1/2 r (0) = αe s . 2. When t = 1, we have D 1/2 x (1) = D 1/2 x (0) -D 1/2 ∇f (x (0) ) = αe s D 1/2 r (1) = αe s -(I -(1 -α)AD -1 )D 1/2 x (1) = α(1 -α)AD -1 e s .3. When t ≥ 1, we haveδ t+1 = 2 1-α -δ t -1 f (t) = 2δt+1 1-α D 1/2 r (t) + δ t δ t+1 D 1/2 ∆ (t) f (t+1) = f (t) + f (t) D 1/2 r (t+1) = D 1/2 r (t) -I -(1 -α)AD -1 f (t) .For LocalCH, we change the update f (t) to the update f (t) St . The final estimate f := D 1/2 x (T ) . • Chebyshev for Katz. We want to solve (I -αA) x = e s . We assume µ = 1αλ max (A), L = 1αλ min (A). Then, the updates of LocalCH for Katz centrality are 1. When t = 0, x (0) = 0, r (0) = e s . To obtain an initial value δ 1 , we have δ 1 = αd max 7 2. When t = 1, we haveEQUATION= (1 -2 L+µ )e s + 2α L + µ Ae s .3. When t ≥ 1, it updates the estimate-residual pair asδ t+1 = 2 αdmax -δ t -1 x(t) = 4δt+1 L-µ r (t) + δ t δ t+1 ∆ (t) x (t+1) = x (t) + x(t) r (t+1) = r (t) -(I -αA) x(t) .For LocalCH, we change the update x(t) to the update x(t) St . The final estimate is then f := x(T )e s .We omit the details of LocalSOR and LocalGD for PPR and Katz as they can directly follow from ( 6) and (10) , respectively. We implement all our proposed local solvers via the FIFO Queue data structure.We examine the different data structures of local methods and explore two types: the First-In-First-Out (FIFO) Queue, which requires constant time O(1) for updates, and the Priority Queue, which is used to implement the Gauss-Southwell algorithm as described in [48, 7] . We test the number of operations needed using these two data structures on five small graphs, including CORA (n = 19793, m = 126842), CORA-ML (n = 2995, m = 16316), CITESEER (n = 4230, m = 10674), DBLP (n = 17716, m = 105734), and PUBMED (n = 19717, m = 88648) as used in [9] and downloaded from https://github.com/abojchevski/graph2gauss.Figure 10 presents the ratio of the total number of operations needed by APPR-FIFO and APPR-Priority-Queue. We tried four different settings; the performance of one does not dominate the other, and the Priority Queue is suitable for some nodes but not all nodes. The number of operations ratio between the APPR (FIFO-Queue) and Gauss-Southwell (Priority-Queue). The number of operations of APPR is defined asT FIFO-Queue = T -1 t=0 S t , T Priority-Queue = K k=0 d u k + log 2 | supp(r (k) )|where K is the total number of push operations used in Gauss-Southwell iteration and supp(r(k) ) = {u : r (k) u ̸ = 0, u ∈ V} and log 2 | supp(r (k) )| is the number of operations needed by the priority queue for maintaining all residuals in r (k) .The dynamic PPR algorithm based on LocalGS is presented in Algo. 2 and 3, initially proposed in [73] . Specifically, Algo. 2 and 3 are the components of LocalGS(Dynamic) for updating the edge event (u, v) while Algo. 4 and 5 are the components of LocagSOR(Dynamic) for update the edge event (u, v). In practice, we follow the batch updates strategy proposed in InstantGNN [75] , which we call Algo. 3 and 5 for list of edge events then call Algo. 2 and 4 once after this batch update of edge events.Algo 2 LocalGS(G, p s , r s , ϵ, α, s)1: while max u r s [u] ≥ ϵd out [u] do 2: Push(u) 3: while min u r s [u] ≤ -ϵd out [u] do 4:Push(u) 5: return (p s , r s ) 6: procedure Push(u): 7: p s [u] ← p s [u] + r s [u] 8: for v in N out (u) do 9: r s [u] ← 0 10: r s [v] ← r s [v] + (1-α)•rs[u] dout[u] Algo 3 LOCALGS(DYNAMIC)(G, p s , r s , ϵ, α, s, (u, v)) 1: Apply Insert/Delete of (u, v) to (p s , r s ). 2: return LocalGS(G, p s , r s , ϵ, α, s) 3: procedure Insert(u, v) 4: p s [u] ← p s [u] * d out [u]/(d out [u] -1) 5: r s [u] ← r s [u] -p s [u]/d out [u] 6: r s [v] ← r s [v] + (1 -α) • p s [u]/d out [u] 7: procedure Delete(u, v) 8: p s [u] ← p s [u] * d out [u]/(d out [u] + 1) 9: r s [u] ← r s [u] + p s [u]/d out [u] 10: r s [v] ← r s [v] -(1 -α) • p s [u]/d out [u] Algo 4 LocalSOR(G, p s , r s , ϵ, α, s, ω) 1: while max u r s [u] ≥ ϵd out [u] do 2: Push(u) 3: while min u r s [u] ≤ -ϵd out [u]p s [u] ← p s [u] + ωr s [u] 8: r s [u] ← r s [u] -ωr s [u] 9:for v in N out (u) do 10:r s [v] ← r s [v]+ω • (1-α)•rs[u] dout[u] Algo 5 LOCALSOR(DYNAMIC)(G, p s , r s , ϵ, α, ω, s, (u, v)) 1: Apply Insert/Delete of (u, v) to (p s , r s ). 2: return LocalSOR(G, p s , r s , ϵ, α, s, ω) 3: procedure Insert(u, v) 4: p s [u] ← p s [u] * d out [u]/(d out [u] -1) 5: r s [u] ← r s [u] -p s [u]/d out [u] 6: r s [v] ← r s [v] + (1 -α) • p s [u]/d out [u] 7: procedure Delete(u, v) 8: p s [u] ← p s [u] * d out [u]/(d out [u] + 1) 9: r s [u] ← r s [u] + p s [u]/d out [u] 10: r s [v] ← r s [v] -(1 -α) • p s [u]/d out [u]We present InstantGNN(LocalGS) and InstantGNN(LcalSOR) in Algo. 6 and Algo. 7, respectively. Additional parameter β is used; in practice, we choose β = 1/2 for all our experiments. Algo 6 InstantGNN(LocalGS)(G, p, r, ϵ, α, s, β)1: while max u |r[u]| ≥ ϵd 1-β out [u] do 2:Push(u) 3: return (p, r) 4: procedure Push(u):5: p[u] ← p[u] + α • r[u] 6: r[u] ← 0 7: for v in N out (u) do 8: r s [v] ← r[v] + (1-α)•r[u] d 1-β out [u]d β out [v]Algo 7 InstantGNN(LocalSOR)(G, p, r, ϵ, α, s, β, ω)1: while max u |r[u]| ≥ ϵd 1-β out [u] do 2:Push(u) 3: return (p, r) 4: procedure Push(u): 5: p[u] ← p[u] + α • ω • r[u] 6: r[u] ← r[u] -ω • r[u] 7: for v in N out (u) do 8: r[v] ← r[v] + ω • (1-α)•r[u] d 1-β out [u]d β out [v]We collected 18 graph datasets listed in Table 3 . To test LocalSOR on GNN propagation, we use three graphs with feature data from [75] .For all experiments, the maximum runtime limit is set to 14,400 seconds. All methods will terminate when this time limit is reached. Specifically, we use the following parameter settings for drawing Figure 1 and Figure 3 .• Experimental settings for Figure 1 . For computing PPR, we use a high precision parameter of ϵ = 10 -10 /m and set α PPR = 0.1. For calculating Katz centrality, we use α Katz = 1/(|A| 2 + 1). The temperature for the Heat Kernel equation is set to τ = 10.• Experimental settings Figure 3 . We use the wiki-talk graph and α = 0.1 for PPR calculation. = ∞ k=0 α(1 -α) k (D β AD 1-β ) k x.It is the same as APPNP [30] when we set β = 0.5.Figure 11 presents the results of GPU-implemented methods for Katz. Table 5 -7 presents more details of participation ratios for PPR, Katz, and HK. The ratios in tables are without normalization. Justification: There are no broader societal impacts, as it focuses on technical contributions. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • 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 is being 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 to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).For the PPR equation, (I -(1 -α)AD -1 )x = αes, we can symmetrize this linear system by rewriting it as (I -(1 -α)D -1/2 AD -1/2 )D -1/2 x = αD -1/2 es. The solution is then recovered by D 1/2 x (t) .More detailed experimental setups and additional results are in the appendix.Note that the approximate local solvers of APPR, AHK, and AKatz will be referred to as GS.We will arbitrarily write αPPR and αKatz as α when the context is clear.To simplify the proof, here r represents D 1/2 r.This is because |λ(A)| ≤ dmax. 


Decision: This paper introduces a set of efficient methods for approximating graph diffusion vectors. Their approach achieves significant acceleration when implemented on GPUs. Empirical results demonstrate that these new methods can achieve speedups of up to one hundred times compared to their global counterparts.

Although the score is borderline, I find this contribution meaningful since it can potentially have an impact in speeding up computation which appears in the industry sector.

Reviews: review: paper_summary: paper_summary: This paper proposes a suite of fast methods to approximately compute graph diffusion vectors such as Personalized PageRank, Katz centrality and the heat kernel. A notable feature of the proposed methods is that they are easily parallelizable and hence can achieve further acceleration on GPU. The authors also provide a running time bound for each method that they introduce. Empirical results show that the new local methods can achieve up to a hundred-fold speedup when compared to their global counterpart.
strength_weakness: strength_weakness: - When compared with APPR, the speedup is not really captured by Theorem 3.3 and Corollary 3.6. The authors only showed that $\overline{\mbox{vol}}(\mathcal{S}_T)/\bar\gamma_T \le 1/\epsilon$, but it requires strict inequality to achieve nontrivial speedup in terms of worse-case running time. It is not clear how tight this bound is in general. If the bound is tight, then in the worst case there is no speedup. The authors should comment on if there are classes of graphs over which there will be a notable gap between the 2 quantities $\overline{\mbox{vol}}(\mathcal{S}_T)/\bar\gamma_T$ and $1/\epsilon$, so that the result provides a meaningful improvement.
- Because the theorems concerning the worst-case running time do not seem to capture a clear improvement over simple baseline local methods, it is unclear where exactly the speedup reported in Table 2 comes from.
- In LocalGD and LocalCH, the authors did not provide an update rule (or even a definition) for $\mathcal{S}_t$. Since $\mathcal{S}_t$ is part of the local diffusion process, the authors should specify how $\mathcal{S}_t$ is updated or defined at each step. I guess that $\mathcal{S}_t$ depends on the termination condition, but I could not find where the authors mention about it in the main paper.
questions: questions: - Line 188: The open problem is recently solved by Martínez-Rubio et al, Accelerated and Sparse Algorithms for Approximate Personalized PageRank and Beyond, 2023.
- When eps goes to 0, do the approximate solutions converge to the exact global solutions?
limitations: limitations: Yes
paper_summary: paper_summary: This paper proposes a new local iterative framework for solving graph diffusion equations (GDEs). Specifically, the framework approximates GDEs through a local diffusion process, leveraging the strong localization properties of diffusion vectors such as personalized PageRank, Katz centrality, and Heat Kernel. The proposed local solvers can achieve sublinear runtime complexity under certain monotonicity assumptions. Empirical results demonstrate that these solvers significantly accelerate their standard counterparts on several large-scale benchmark datasets.
strength_weakness: strength_weakness: 1. It is not clear how graph structures, such as sparsity and spectral properties, impact the runtime complexity of the local solvers. A discussion on the potential influence of graph structures on runtime complexity would be helpful.

2. The runtime complexity analysis assumes that the updates of local solvers satisfy monotonicity properties. However, LocalCH does not satisfy these properties during updates. Establishing runtime bounds for LocalCH may require different techniques.

3. In Figure 7, when $\epsilon \geq 2^{-31}$, the running time of LocalGD is the same as GD, indicating that LocalGD may not speed up the standard counterparts when $\epsilon$ is sufficiently small.
questions: questions: Please refer to Weaknesses.
limitations: limitations: NA
paper_summary: paper_summary: This paper proposes a novel framework for approximately solving graph diffusion equations using a local diffusion process. In addition,  the proposed method can effectively localizes standard iterative solvers by designing simple and provably sublinear time algorithms.
strength_weakness: strength_weakness: - In this paper, the authors explore 18 different graphs, can the authors show/provide which local GDE solver achieves better performance on which type(s) of graphs? 
- Computational complexity/cost is missing.
questions: questions: See comments in Weaknesses.
limitations: limitations: Not applicable.
paper_summary: paper_summary: The paper, "Faster Local Solvers for Graph Diffusion Equations," addresses the efficiency of computing Graph Diffusion Equations (GDEs) such as Personalized PageRank (PPR), Katz centrality, and the Heat kernel, which are essential for various graph-related problems like clustering and training neural networks. Traditional methods for solving GDEs are computationally intensive for large-scale graphs. This paper introduces a novel framework for approximating GDEs using local diffusion processes, which significantly reduces computational time and improves scalability by leveraging the localization property of diffusion vectors. The proposed local solvers are highly parallelizable and suitable for GPU implementation, offering up to a hundred-fold speed improvement and applicability to large-scale dynamic graphs. The paper also discusses the potential for these methods to enhance local message-passing mechanisms in Graph Neural Networks (GNNs).
strength_weakness: strength_weakness: Precision Limitations: While the paper shows significant speedups for lower precision settings, the performance gain diminishes as higher precision is required. This limitation could affect the applicability of the methods in scenarios where high precision is crucial.

Sequential Nature: (Please correct if I am wrong.) Despite improvements, the reliance on sequential updates in some methods (like LocalSOR) still poses challenges for achieving maximal parallel efficiency.

Complexity of Analysis: The runtime analysis for some of the proposed methods, such as the Local Chebyshev method, is noted as complex and remains an open problem. Maybe simplifying this analysis or providing more intuitive explanations could enhance the paper’s accessibility.
questions: questions: What are the potential trade-offs between precision and computational efficiency, and how can practitioners balance these in practical applications?

The paper includes quite some datasets to test, with biggest of ogbn-papers100M. How does the proposed method handle very large-scale dynamic graphs in real-time applications, and what are the potential bottlenecks?
limitations: limitations: See above

Comments: title: We thank all reviewers for their time and efforting, we will keep improving our submission.
Request for Clarification on Remaining Concerns (We are confused)
Thank you.
comment: We sincerely thank all reviewers for their time and effort in the discussion and for helping us; we will keep improving our submission.
I'd like to thank the authors for their detailed responses. My questions have been properly addressed. I will increase my score.
We sincerely appreciate your response and your careful reading of our rebuttal. However, we are still uncertain whether all your concerns have been sufficiently addressed. Based on your score, it seems there are reasons to reject the paper, such as limited evaluation, that outweigh the reasons to accept it, such as a good evaluation. But, currently, you only have one question and one misunderstanding point. So, we are confused.

To facilitate a valid discussion, please specify any other concerns you might have so we can address them appropriately. Thank you.
I appreciate the author for the detailed response. After carefully reading the rebuttal, I am retaining my score.
Thank you for your response. I continue to recommend acceptance and retain my score.