Massively Parallel Single-Source SimRanks in O(log N) Rounds

1 Introduction

Evaluating the structural similarity between two nodes in a graph is fundamental in plenty of data mining and machine learning tasks. Examples include recommendation systems [Chen et al., 2022; Mo and Luo, 2021], avoiding customer churn [Luo et al., 2019], spam detection [Badola and Gupta, 2021], link prediction [Zhou et al., 2019; Xu et al.,

* The authors are ordered alphabetically.

2017; Hours et al., 2016] and graph mining [Zhu et al., 2024; Jin et al., 2011; Liao et al., 2022; Zhu et al., 2022]. Among many similarity measures, SimRank [Jeh and Widom, 2002] is one of the most widely adopted measures over graphs. SimRank is defined based on the intuition that two nodes are similar only when their neighboring nodes are similar. Formally, it uses the following recursive equation to compute the SimRank between two nodes u and v, where c ∈ (0, 1) and I(u) denotes the in-neighbor set of node u.

Since it was proposed by Jeh and Widom [Jeh and Widom, 2002], SimRank has gained increasing popularity in various application domains, e.g., social analysis [Zheng et al., 2013], nearest neighbor search [Lee et al., 2012], as well as recommender systems [Antonellis et al., 2008; Chuanyan and Xiaoguang, 2021; Symeonidis et al., 2021].

Distributed SimRank Computation. We focus on singlesource SimRank computation, whose importance has been uncovered in a plethora of recent studies [Wang et al., 2020b; Li et al., 2015; Wang et al., 2020a; Shi et al., 2020; Wang et al., 2021; Maehara et al., 2014; Tian and Xiao, 2016; Kusumoto et al., 2014]. Given a graph G, a single-source SimRank for source node u evaluates the SimRank values between u and all other nodes in the graph. Single-source SimRank computation is widely used in applications where a ranking of the other objects with respect to an object is required. For example, it can be applied in search engines to locate the most similar web pages to a given one [Fogaras and R´acz, 2005], or in social network services to recommend new friends to a given user [He et al., 2010; Nguyen et al., 2015], or act as a ranking measurement to cluster objects [Cai et al., 2008].

Given a graph of n nodes, computing single-source SimRank is challenging for large graphs because it inherently involves O(n) times of pairwise SimRank evaluations, each of which can already be too costly. Particularly, following the recursive form in Equation 1, computing s(u, v) requires accessing many pairs of nodes in the graph, leading to O(n2 ) complexity. To address the efficiency issue, recent works [Jiang et al., 2017; Zhang et al., 2017; Shao et al., 2015] employ a random-walk-based approach to approximate the SimRank values. The main idea is to translate the SimRank computation into estimating the meeting probability of two decay-based random walks from the two source nodes (See Section 2.2 for more details). The computational complexity is then dependent on the number of random walks sampled, achieving a significant speed-up. However, even with the improved approach, it is still challenging to compute single-source SimRanks when n is large.

Therefore, it is increasingly popular to apply distributed computation [Wang et al., 2020b; Li et al., 2015; Song et al., 2018] to SimRank computation, which involves multiple machines to compute SimRank values in a collaborative manner, ultimately scaling up the computation to large graphs. In modeling distributed computation, the local computation cost within each machine is typically dominated by the synchronization cost among the machines. Hence, by convention, the local computation cost is omitted in MPC analysis and the focus is on reducing communication rounds. Furthermore, there is an intrinsic trade-off between the number of communication rounds and space-per-machine.

Open Problem. Existing distributed SimRank algorithms (e.g., [Li et al., 2015; Song et al., 2018; Wang et al., 2020b]) mostly focus on empirical evaluation, and nontrivial theoretical analysis with provable approximation result guarantees is rarely given. We analyze existing representative distributed SimRank algorithms and summarize their requirements regarding communication rounds and machine space in Table 1*. In a nutshell, they require either Ω(n) per-machine space (UniWalk [Song et al., 2018]) or Ω(log n) communication rounds (DISK [Wang et al., 2020b] and CloudWalker [Li et al., 2015]). We remark that overcoming the Θ(logn) round complexity barrier in MPC for natural problems is usually challenging and it has attracted tremendous interest [Luo et al., 2022; Lacki et al., 2020; Luo, 2019; Behnezhad et al., 2019; Czumaj et al., 2021; Luo, 2020] to improve the distributed computation down to O(poly(log log n)) rounds. Hence, a natural open theoretical problem is raised:

Given an error threshold ϵ ∈ (0, 1), let an approximate single-source SimRank algorithm be an algorithm that outputs SimRank values with ϵ absolute errors. Is there a distributed approximate single-source SimRank algorithm over a graph of n nodes that can be finished in o(logn) rounds, while each machine only needs o(n) space?

OurMainResults. In this paper, we give a positive answer to the aforementioned open problem and present a distributed single-source SimRank algorithm that suits the Massively Parallel Computing (MPC) model [Karloff et al., 2010]. Particularly, for a graph of n nodes and m edges, we focus on the MPC model that involves a set of machines, each having a sub-linear space S = nα for some α ∈ (0, 1). Our main results can be stated by the following theorem. To the best of our knowledge, this is the first distributed single-source SimRank algorithm that achieves sub-logn communication rounds while only requiring strongly sub-linear per-machine space. We leave the proofs of lemmas and theorems in the Appendix of our technique report [Luo and Zhu, 2024].

* For the limitation of space, we leave the related work analysis in Appendix C. As the authors did not give an analysis based on the MPC model, we analyze them on our own.

Theorem 1. Given a source node u in a graph of n nodes and m edges, there is an algorithm that computes ϵ-absolute-error guaranteed SimRank values between u and all the other nodes using M machines in O (log2 log n) communication rounds with high probability † . This algorithm only requires the space per machine is S = nα for some α < 1, and M is some value of ˜O((m+ n)/S).

2 Preliminaries

2.1 Distributed ComputationModel

Distributed computation is one of the most important techniques to address various computation tasks [Luo et al., 2023; Wang et al., 2020b; Li et al., 2015; Song et al., 2018; Luo et al., 2022; Luo et al., 2014; Luo et al., 2012; Klauck et al., 2014]. In recent years, Massively Parallel Computation (MPC) [Karloff et al., 2010] becomes a popular theoretical framework in modeling the complexity of a distributed algorithm, because it closely simulates the situation of general distributed computation. An MPC model has three important parameters: the input data size F , the number of involved machines M, and the space capacity (words) S on each machine. The frequently used notations are listed in Table 2 in the Appendix. It is required that when given F and S, the number of machines M should be of ˜O( FS ), where ˜O(·) hides a polylogarithmic factor. Our main focus on the MPC model derives from several perspectives as follows:

Space. Consider an input graph of n nodes. MPC for graph algorithms can be categorized into three types: strongly super-linear space (S = n1+ω ) for some constant ω > 0, near linear space (S = Θ(n)), and strongly sub-linear space (S = nα) for some constant α ∈ (0, 1). A super-linear model can typically be employed with a local algorithm, which loses the generality to be employed on a large scale of data. Hence, many studies focus on sub-linear settings that better capture the scalability of a distributed system. In this paper, we focus on the strongly sub-linear setting where each machine cannot even store the whole set of graph nodes.

Communication Rounds. The computation in the MPC model is based on communication rounds. Initially, each edge is randomly assigned to a machine. We assume each node has an integer ID from 1 to n, and each machine has an integer ID from 1 to M. The computation proceeds in synchronized rounds. At the beginning of a round, each machine may receive the messages sent from some other machines in the previous rounds. During a round, every machine conducts local computation based on its local data or messages received. Then each machine will send the computed results, packed as messages, to target machines. Each machine creates message packages to be routed onto the network and hence the size of messages sent/received per machine in one round shall not exceed the space capacity S. A new round starts only after the end of the previous round. The term round complexity [Ghaffari et al., 2019] refers to the number of rounds required for program execution and serves as a measure of distributed computation cost, primarily driven by the expensive communication among machines. The main principle of designing algorithms in the MPC model is to achieve low round complexity subject to the space constraint on each machine.

† We say an event happens with high probability, if there exists a constant τ > 0 such that the event happens with probability at least 1- 1 nτ , where n is the number of graph nodes.

Table 1: Communication rounds of distributed single-source SimRank Computation.

Existence of an MPC Algorithm. With a per-machine space of S = nα, we acknowledge the existence of an MPC algorithm such that when n ≥ n0 (with n0 being a constant), a distributed algorithm can operate within this space constraint using M machines, where M = ˜O( m+n S ). For ease of analyzing the MPC algorithms, we first give a preparation lemma regarding machine space expansion, as follows.

Lemma 1. If there exists an O(f(n))-round distributed graph algorithm that works for per-machine space S = Θ(nα ) for any 0 < α < 1 using M = ˜O( m+n nα ) machines, and f(n) is a function not related to α (i.e., α only contributes a constant factor to the round complexity), then there exists an MPC algorithm that works for per-machine space S = nα with the same round complexity.

This lemma eliminates the obstacle of analyzing round complexity when there is a constant factor expansion of the machine space. For example, consider that we have an algorithm with a certain distributed algorithm such that (1) its round complexity hides α factors; (2) applies to any α ∈ (0, 1) with M = ˜O( m+n nα ), and (3) the space per machine is Cnα for some constant C. Then Lemma 1 guarantees the existence of an MPC algorithm with nα machine space.

2.2 Approximate SimRank Computation

Given a directed graph G = (V,E), and let n = |V | and m = |E|. Following [Wang et al., 2020a; Shi et al., 2020], we aim to compute approximate SimRank values with constant errors. In particular, given any source node u ∈ G and a constant error ϵ, we aim to compute the SimRank values between u and any other node v ∈ G, such that |s(u, v)-˜ s(u, v)| ≤ ϵ, where s(u, v) denotes the true SimRank value and ˜ s(u, v) denotes the estimated value.

Calculating the SimRank value iteratively according to Equation 1 may occupy large memory space and incur a high computation cost. Therefore, state-of-the-art approaches employ the following √c-decay walk-based computation, first proposed in [Tian and Xiao, 2016]:

Lemma 2. For any two nodes u, v ∈ G, the SimRank between u and v is equal to the meeting probability of two √c-decay walks starting at u and v on the reverse graph ¯G, where two walks meet if there exists an integer i ≥ 0, such that the i-th step of the two walks visit the same node.

Here, a √c-decay walk on ¯G from a node u is a traversal on ¯G such that at each step of the walk, it has 1- √c probability to stop at the current node, and otherwise jumps to the next node that is a uniformly chosen out-neighbor of the current node. Here we define a length-l walk as a path that includes l + 1 nodes and l edges. We also say a length-l walk has l steps. Particularly, we have (√c)i · (1 - √c) probability to generate a length-i walk from a given starting node.

Monte Carlo Method. By sampling N pairs of √c-walks from u and v, and if H pairs of the walks meet, then H N is an estimator of s(u, v). The sampling number N controls the estimation accuracy since more random walks approximate the value of s(u, v) more accurately.

3 O(log2 logn)-Round Algorithm

It is not too challenging to design O(logn)-round MPC algorithm (See Appendix E for details), but further improving the round complexity can be difficult. To explain, we outline challenges of direct adaptation of √c-walk based approach. First, the maximum length of a √c-decay walk can be infinite; a straightforward method to compute length-l in a distributed environment easily entails l communication rounds because one step may require one communication round when some neighbors are located in a different machine. Such design easily leads to Ω(logn)-round algorithms. Second, evaluating single-source SimRank values for node u requires running random walks from both u and all the other nodes, resulting in a large number of random walks being conducted. This can lead to a high space cost per machine because some hub nodes are prone to be passed by many random walks. Therefore, a careful design of the MPC algorithm is required to guarantee a small number of rounds and a low space cost for each machine. Third, the random walks, once computed in the MPC model, are stored in different machines. As such, detecting whether two walks meet may overload the machine regarding space cost.

Algorithm 1: Overall algorithm

Input : Graph G = (V,E); source node u; decay factor c; Montel-Carlo failure probability δ; accuracy error ϵ Output: SimRank scores ˜ s(u, v) for each v ∈ V

• 1 Store G across multiple machines;
• 2 Reverse the edges of G within each machine to form ¯G;
• 3 l ← log 1 √c n;

• 8 Shuffle and decompose random walks (See Sections 4.2 and 4.3);
• 9 Run Algorithm 2 to calculate ˜ s(u, v) for each v ∈ V ;

3.1 New Interpretation of Walk-based Approach

To address the aforementioned challenges, we reinterpret the SimRank computation between node u and node v in a batch manner using the following three operations, which are more MPC-friendly:

(a) Random walk generation. Instead of generating √cdecay random walks whose lengths are non-deterministic, we generate random walks with deterministic length distribution. Particularly, to generate N √c-walks from u, we will generate N · √ci (1 - √c) length-i walks for i ≥ 0, following the corresponding geometric distribution of walk lengths. Similar operations are conducted for random walks starting from v. Fixing the random walk lengths is more MPC-friendly. As we will show in Section 4.1, we can generate sufficient such random walks both from the source node u and each node v ∈ V in only O (log2 log n) MPC rounds.

• (b) Random walk shuffling and matching. Since the random walks from u (or v) are generated in the order of increasing lengths (due to Step (a)), pairing up the random walks from u and v directly for meeting detection is not valid for estimating SimRanks as the i-th walks from u and v are correlated (they have the same length). In order to maintain the randomness, we shuffle the random walks above.
• (c) SimRank computation. Last, we compute the SimRanks by detecting how many pairs of walks meet and estimating the SimRank value s(u, v) by Lemma 2.

3.2 Overview of Main Algorithmic Steps

Based on the new interpretation, we give our main algorithmic steps in Algorithm 1, which consists of the following five main stages. For ease of presentation, we first outline the main idea in this section and defer the detailed MPC-related operations in Section 4.

Initial State. The SimRank algorithms are based on the topology information of the input graph. Initially, the whole graph should be partitioned across different machines, as assumed by the MPC model. We assume that each machine holds a random partition of edges of the graph (Line 1 in Algorithm 1), and hence each machine holds roughly m M edges.

Algorithm 2: Detect meeting and calculate SimRank

Input : Tuple set DG from Section 4.3; source node u Output: Estimated SimRank score ˜ s(u, v) for each v ∈ V

We also reverse the edges so that the graph represents ¯G (Line 2 of Algorithm 1).

Parallel Random Walks Generation. This stage corresponds to Operation (a) mentioned earlier. As we need to compute single-source SimRank values from any source node u to all other nodes, we generate N √c-walks from each node. Further, we truncate them at log1/√c n length (Line 3 of Algorithm 1). By setting N = log 2n 2(ϵ-3/n)2 , we show that the truncated walks only cause negligible influence on the final accuracy (see Section 4), and we can guarantee the ϵ error bound of SimRank values. Another problem is when generating Ni = N ·(√c)i (1-√c) length-i walks, Ni may not be an integer. In Section 5 we give a rounding technique and prove that the rounding still guarantees the error bound. We also prove that these random walks can be obtained in the MPC model using O(log2 log n) communication rounds (Lines 4-7 of Algorithm 1 and details in Section 4.1).

Shuffling Random Walks. It's important to note that in Lemma 2, the lengths of the random walks originating from u and v may differ, as they are independently and randomly generated. To guarantee the randomness of each pair of random walks in our Monte Carlo simulation, the generated random walks have to be shuffled (Line 8 of Algorithm 1). The shuffling operation in the MPC model is not as trivial as in a single machine because the generated walks from the same node can be located in different machines. Shuffling these walks incurs communication between machines, which may violate the space bound in each machine.

Decomposing Random Walks. After shuffling, we pair the j-th walk from the source node u with the j-th walks from other nodes respectively, and detect whether two walks in each pair meet. Detecting whether two walks meet is relatively simple in the single-machine setting because all information can be loaded locally. However, in a distributed environment, we need to detect whether two paired-up walks meet while the walks are stored at different machines. This poses drastically different challenges. To address the challenges, we decompose the generated random walks into tuples, each containing one visited node in the walk (except the starting node). Each tuple contains five elements to convey all the information needed to detect whether the walk intersects another walk in the specified node. As shown in Figure 1, there is a walk {u, v1, v2, v3} decomposed into three tuples (j, 1, v1, 3, u), (j, 2, v2, 3, u), (j, 3, v3, 3, u), where j implies that the tuple is decomposed from the j-th walk of those starting from u. The remaining elements in the tuple are the walk step, visited node at the particular walk step, the total steps in the walk, and source node, respectively.

Computing SimRank Values. Evaluating the SimRank value between u and v is based on the probability that two walks from u and v meet. This can be detected by sorting the decomposed tuples (Line 9 of Algorithm 1→ Line 1 of Algorithm 2). Figure 1 shows the j-th walks from u and v respectively. The two walks meet because they share the same node v2 in their third walk steps. Corresponding to v2, there are two decomposed tuples from the two walks, denoted by (j, 2, v2, 3, u) and (j, 2, v2, 4, v). We note that they share the first three elements. In general, it is easy to verify that if two walks meet, there must be two decomposed tuples, each from one walk, that share the first three elements in the tuple. More formally, suppose there are two j-th walks {v0, v1, . . . , vl1 } and {u0, u1, . . . , ul2 }(0 ≤ l1, l2 ≤ log1/√c n) starting from v0 and u0. Walk {v0, v1, . . . , vl1 } is decomposed into l1 tuples, denoted by (j, i, vi, l1, v0) (1 ≤ i ≤ l1); Walk {u0, u1, . . . , ul2 } is decomposed into l2 tuples, denoted by (j, k, uk, l2, u0) (1 ≤ k ≤ l2). Then, if the two walks meet, or equivalently, they share the same node vi in the same step i, there must be vi = ui and the corresponding two tuples (for vi and ui) share the first three elements. Again, the challenges are in a proper rearrangement of the tuples across the machines, and the collection and aggregation of the walkmeeting cases within each machine.

With the above five stages, we highlight the major differences between our approach and the classic √c-walk-based algorithms: (a) the lengths of the walks are regularized to geometric distributions, and longer walks are truncated and we show that this would not violate the error tolerance in Section 5; (b) most steps are redesigned non-trivially so that all the steps can be efficiently implemented in the MPC model.

4 DetailedMPCOperations

4.1 Parallel RandomWalks Generation

Performing random walks in the MPC model will incur communication rounds because the next sampled neighbor for each node can be located in different machines. Let us consider the situation of a node v0 sampling its neighboring node v1 in the walk. In the MPC model, as v0 and v1 can be in different machines, forming a walk of v0 → v1 may need one communication round. In general, forming an l-step random walk in the MPC model typically incurs l communication rounds. For the purpose of reducing communication rounds, we first introduce a result provided by WRME [Lacki et al., 2020], and then give our extension in Theorem 3:

Theorem 2. Let G be a directed graph. Let N and l be positive integers such that l = o(S)/ log3 n, where S is the available space per machine. For the task that samples N independent random walks of length l starting from each node v in G, there exists an MPC algorithm that runs in O (log2 log n + log2 l) rounds and uses O (m+n1+o(1) l3.5 +Nnl2+o(1)) total space and strongly sub-linear space per machine S = nα (0 < α < 1). The algorithm is an imperfect sampler that does not fail with probability 1 - O (n-1 ).

Essentially, Theorem 2 states that the generation of a length-l random walk from every node can be roundefficient in the MPC model. Particularly, it takes only O (log2 logn+ log2 l) communication rounds to generate N length-l random walks, each with different source node. However, even with Theorem 2, we cannot directly give a reasonable round-complexity for SimRank evaluations if we use the original √c-decay walk-based method. The reason is that the length l of a √c-decay walk can be infinite. To address this issue, as we show in Algorithm 1 Line 3 and Line 5, we carefully design the truncated walk length (Line 3) and the number of random walks to be sampled for each walk length (Line 5). The truncated length guarantees a reasonable bound of O (log2 logn) communication rounds when l = log1/√c n, and the number of random-walk samples ensures the SimRank estimation accuracy (we will formally prove it in Section 5).

Parallel RandomWalks in O(log2 log n) Rounds. To effectively apply Theorem 2, we let Task-i be the generation of Ni length-i walks from every node. We let all Task-i (1 ≤ i ≤ log1/√c n) be conducted in MPC in parallel. We apply log1/√c n WRME algorithm [Lacki et al., 2020] instances concurrently in the MPC model, where each instance corresponds to one task.

Unfortunately, concurrently conducting multiple WRME algorithms will incur a space cost higher than nα in each machine because by default each algorithm instance can incur a local cost up to nα. To address this issue, we let β satisfy that nβ = nα/ log1/√c n, and we apply the WRME algorithm concurrently for each Task-i with S = nβ . Particularly, each machine space nα is divided into subspaces of size nβ to compute Task-i. By careful analysis of the round complexity and total space cost, we give the following theorem.

Theorem 3. The random walk generation in Algorithm 1 can be done in O (log2 log n) rounds with high probability, with machine space S = nα such that nα > log5 n/(log 1 √c ) and a total space O (m log n + n1+o(1) log4.5 n + n log3+o(1) n ).

We note that the WRME algorithm is an imperfect sampler that fails with at most a probability O(n-1 ). We only apply O(log n) times the sampler, which can still easily guarantee that our algorithm is successful with high probability (because O(log n) = o(nτ ) for any constant τ > 0).

4.2 Shuffling RandomWalks

We need to shuffle the generated N random walks sourced at each node to fully mimic the behavior of generating a √cwalk. The reason for shuffling is that when we select two random walks that start from two nodes and calculate their meeting probability, the lengths of these two walks are not bound to be the same. This shuffling operation helps us guarantee the randomness in our Monte Carlo simulation of the generated √c-decay walks.

Shuffling in O(1) Rounds. We note that sorting in MPC can be done in O(1) rounds [Goodrich et al., 2011] as long as the number of items is O(nη ) for some constant η > 0. Observe that N = O( log n ϵ2 ), and hence there exists η such that nN = O(nη ). Hence, sorting the nN walks all together can be finished in constant rounds in MPC, where the comparison of two walks is based on the comparison of the corresponding nodes in the walks. As such, all the walks sourced at the same node will be clustered.

For sufficiently large n we have nα > 2N log1/√c n, and the walks sourced at the same node, occupying a space of N log1/√c n, can be held in one machine. Since all the walks are sorted and nα N log1/√c n may not be an integer, each set of N walks sourced at the same node may go across two machines. Let Wu be the set of walks sourced at node u. We discuss two cases to shuffle the walks of the same source node.

Case 1: If the walks inWu are fully located in a machine, shuffling is done locally. We shuffle all sets of walks, and for each set we assign numbers from 1 to N to each walk in Wu after shuffling.

Case 2: If the walks sourced at the same node are located at two consecutive machines. Without loss of generality, we denote the source node as u, and there are N0 walks in the first machine and N -N0 in the second. We then swap those N -N0 walks in the second machine with the last N -N0 walks sourced at node u - 1, as illustrated in Figure 2 in the Appendix. Since nα > 2N log1/√c n for sufficiently large n, it is guaranteed that those N -N0 walks sourced at u-1 are located in the first machine and have been shuffled locally in Case 1. Once swapped, the shuffling within Wu can be done locally. We note that by guaranteeing nα > 2N log1/√c n we can make sure the swaps are not conflicting and thus can be done in one round.

Any set Wu can be shuffled in Case 1 or Case 2. During shuffling, each walk is labeled by 1 to N sequentially, and the j-th walk refers to the walk assigned with a number j.

4.3 Decomposing RandomWalks

Decomposition in O(1) Rounds. Recall that each length-t walk {v0, v1, . . . , vt} is to be decomposed into t tuples each of which has 5 elements: (j, i, vi, t, v0), for 1 ≤ i ≤ t. Tuple (j, i, vi, t, v0) implies that in the j-th random walk of the walks starting from v0, the i-th step of the walk visits node vi and the walk has overall t steps. In the MPC model, the decomposition of a random walk can be done within the machine holding the walk. While each walk will generate multiple tuples, the total space is only amplified by a constant factor (i.e., 5 times) because each tuple corresponds to one visited node along the walk. Hence, this operation will not incur any communication among the machines, and only has a constant expansion of the space cost, which has a negligible effect due to Lemma 1. We denote the set of all the decomposed tuples by DG.

4.4 DetectingMeeting-Walks

As the decomposed tuples are stored at different machines, challenges exist if we use these tuples to detect whether some pairs of random walks meet. Algorithm 2 shows the pseudocode of detecting walk-meetings and computing SimRank.

Sorting in O(1) Rounds. Our first step is to sort all the tuples across the machines (Line 1). The sorting is based on the multi-dimensional sorting because the tuple contains five elements, i.e., sorting is first done based on the first element, and then for the elements that share the same first element, sorting is based on the second element, and so on. Elements are sorted based on their node IDs. We also intentionally let the node ID of source node u be smaller than other nodes to guarantee that walks sharing the same first three elements start with the walk sourced at u.

Recall in Section 3 we show that if two walks meet, then there must be two decomposed tuples, each from one walk, sharing the first three elements. After sorting the tuples, the tuples sharing the first three elements will be clustered together and the tuple whose fifth element is u, if exists, will be placed in the first position in the cluster. Particularly, suppose we have tuples (j, i, v, l1, u), (j, i, v, l2, v1), (j, i, v, l3, v2) that share the first three elements, implying that the j-th walk started from u meets the j-th walks started from v1 and v2. These three tuples will be clustered in the order of (j, i, v, l1, u), (j, i, v, l2, v1), (j, i, v, l3, v2) after sorting. The formal result is as follows.

Lemma 3. If the j-thwalkW started from the query source node u meets the j-th walk W ′ started from v, then there must be a decomposed tuple T of W and a decomposed tuple T ′ of W ′ that share the first three elements. Furthermore, after sorting the tuples, T is the closest tuple that is decomposed fromW and comes before T ′ .

Walk Paired-Up in O(1) Rounds. To compute the SimRank between u and v1, we need to pair up (j, i, v, l1, u) and (j, i, v, l2, v1); similarly, we also need to pair up (j, i, v, l1, u) and (j, i, v, l3, v2) to compute the SimRank between u and v2. Pairing-up can be challenging because the tuples are stored at different machines, and hence communication between machines is unavoidable. In the MPC model, a machine that communicates with other machines would require the same space as the size of communication messages. To reasonably bound the communication cost between the machines, we employ the following PREDECESSOR procedure [Behnezhad et al., 2019] to couple the tuples sharing the same first three elements, which can be finished in O(1) communication rounds.

PREDECESSOR: Considered an ordered list of tuples such that each tuple is labeled by 0 or by 1. Then, for each tuple T ′ labeled by 0, PREDECESSOR associates the closest tuple T labeled by 1 such that T comes before T ′ in the ordering. PREDECESSOR can be implemented in O(1) MPC rounds with nα space per machine, for any constant α > 0.

To apply PREDECESSOR for pairing up the tuples, each machine can run a local algorithm to distinguish the tuples that are generated from the random walks starting at node u and others. As shown in Figure 1, all machines assign label 1 to the tuples whose fifth elements are u in parallel. Otherwise, the tuples are labeled by 0. By Lemma 3, we can apply the PREDECESSOR procedure to associate the closest Label-1 tuple (i.e., the tuple generated from the walk sourced at u) to each of its following Label-0 tuples. Here, the association means that the machine holding the Label-0 tuple is aware of its closest Label-1 tuple before it. For each such associated pair T and T ′ , if their first three elements are the same, then the corresponding walks of T ′ and T meet.

4.5 Computing SimRanks

The paired-up tuples in the previous step will contribute to the corresponding SimRank value. For example, suppose two walks from u and v meet, then this pair of walks contributes to the SimRank value s(u, v), because s(u, v) is estimated by the meeting probability of the walks from node u and node v. The challenge is how to concurrently distribute and aggregate these SimRank contributions across machines. Note that it is crucial to let those SimRank contributions corresponding to the same pair of nodes (e.g., (u, v)) be processed in the same machine, to avoid repeat counting of walk-meets. For example, walks {u, u1, u2, w} and {v, u1, u2, k} meet at both node u1 and node u2. To ensure accuracy, the event of each walk-meet should be counted only once. For this purpose, we let each machine specifically handle SimRank evaluations of nM nodes with respect to u. Without loss of generality, we assume the node IDs are from 1 to n. Then, the i-th machine handles the node set with node IDs in [ n(i-1) M +1, ni M ].

Computing SimRanks in O(1) Rounds. The detailed MPC operations of computing SimRanks are performed as follows (also illustrated in Algorithm 2 Lines 8 - 11). Using the PREDECESSOR procedure, for each tuple (i, j, vj, l2, v) with label 0, it can be linked to the closest tuple (i, k, vk, l1, u); if j = k and vj = vk, then that indicates the meet of two walks represented by the two tuples. Then, a walk-meet message (j, u, l1, v, l2) will be sent (by the machine holding tuple (j, i, vi, l1, u)) to the target machine that is responsible to computing SimRank s(u, v). Here, the walk-meet message (j, u, l1, v, l2) indicates that a length-l1 walk sourced at node u and a length-l2 walk sourced at node v meet. Also, assuming u and v are node IDs from 1 to n, the target machine that the message is sent to is the ⌊ uM n ⌋- th machine because this machine will compute s(u, v). Each machine may then receive multiple copies of the walk-meet message (j, u, l1, v, l2), and only one of them will be kept because each walk-meet event should be counted only once. For each unique message (j, u, l1, v, l2) received, we add a value |Zu,l1|·|Zv,l2 | | ˆ Zu,l1 |·| ˆ Zv,l2 | · 1 z to the SimRank value ˜ s(u, v), where z = log 2n 2(ϵ- 3n )2 is the expected number of walk samples, ˆ Zu,l denotes the set of actually generated length-l random walks sourced at u, and Zu,l is the expected set of length-l random walks to be generated sourced at u.

We distinguish actual and expected samples because Zu,l may not be an integer. Particularly, we note that | ˆ Zu,l| = ⌈z · (√c)l · (1 - √c)⌉, and |Zu,l| = z·(√c)l ·(1-√c). When z · (√c)l · (1 - √c) is not an integer, the expected walk set Zu,l is defined based on size-⌈|Zu,l|⌉ set Z1 and size-⌊|Zu,l|⌋ set Z2, such that there is probability |Zu,l|-⌊|Zu,l|⌋ that one walk of Zu,l is selected from set Z1 and otherwise selected from set Z2.

Space Analysis for Walk-Meet Messages. The walk-meet messages received by each machine are at most log1/√c n times the number of the random walks started at the nodes handled by the machine. To see this, in the worst case, for each node of a walk W generated from v handled by the machine, there is a walk W ′ from u that meets W at the node, creating a walk-meet message. Since each machine is responsible to compute SimRank values for n/M nodes, the space per machine for receiving the message is O ( nM log 1 √c n · log 2n 2(ϵ- 3n )2 ) = O ( n log2 n M(ϵ- 3n )2 ). The total space cost for walk-meet messages is therefore O ( n log2 n M(ϵ- 3n )2 · M ) = O ( n log2 n (ϵ- 3n )2 ).

5 Summary of the Results

Round Complexity. The round complexity is dominated by the generation of parallel random walks, which is O(log2 log n) by Theorem 3.

Space Complexity. The total space cost for walkmeet messages is dominated by the random-walk generation

(see Theorem 3). We note that the total space cost is ˜O(m+ n), indicating the existence of M = ˜Θ( m+n

S ) for S = nα. All these results give us Theorem 4 as well as Theorem 1.

Theorem 4. Given a graph of n nodes and m edges, and a constant error ϵ > 0, for sufficiently large n, Algorithm 1 can be run in MPC using O (log2 logn) communication rounds, with a total space per-round as O (m log n + n1+o(1) log4.5 n + n log3+o(1) n 2(ϵ- 3n )2 ). The space per-machine needed is S = nα such that nα > log5 n/(log 1 √c ). Algorithm 1 is an imperfect sampler that fails with probability at most O (n-1 ).

Accuracy Analysis. The accuracy guarantee of our random-walk sampling techniques resembles that of a simple Monte Carlo method, but we also incorporate the techniques of length-truncation and the rounding of the number of walks. We only include the main result as follows. The proof and the methodology for further reducing the threshold 3n are detailed in the Appendix.

Theorem 5. Algorithm 1 outputs SimRank values ˜ s(u, v) with error at most ϵ (ϵ ≥ 3n ) and with probability at least 1 - 1n .

Acknowledgments

This research is supported by the Ministry of Education, Singapore, under its AcRF Tier-2 Grant (T2EP20122-0003).