Efficient Estimation of Heat Kernel PageRank for Local Clustering

TL;DR

This paper introduces TEA and TEA+ algorithms for local graph clustering that efficiently estimate heat kernel PageRank with strong theoretical guarantees, significantly outperforming existing methods on large-scale graphs.

Contribution

The paper proposes novel algorithms TEA and TEA+ that provide accurate HKPR estimation with theoretical guarantees and improved efficiency for large-scale graph clustering.

Findings

TEA+ is over four times faster than previous algorithms on benchmark datasets.

TEA+ achieves an order of magnitude speedup on large graphs like Twitter and Friendster.

The algorithms maintain high clustering quality while reducing computational time.

Abstract

Given an undirected graph G and a seed node s, the local clustering problem aims to identify a high-quality cluster containing s in time roughly proportional to the size of the cluster, regardless of the size of G. This problem finds numerous applications on large-scale graphs. Recently, heat kernel PageRank (HKPR), which is a measure of the proximity of nodes in graphs, is applied to this problem and found to be more efficient compared with prior methods. However, existing solutions for computing HKPR either are prohibitively expensive or provide unsatisfactory error approximation on HKPR values, rendering them impractical especially on billion-edge graphs. In this paper, we present TEA and TEA+, two novel local graph clustering algorithms based on HKPR, to address the aforementioned limitations. Specifically, these algorithms provide non-trivial theoretical guarantees in relative error of HKPR values and the time complexity. The basic idea is to utilize deterministic graph traversal to produce a rough estimation of exact HKPR vector, and then exploit Monte-Carlo random walks to refine the results in an optimized and non-trivial way. In particular, TEA+ offers practical efficiency and effectiveness due to non-trivial optimizations. Extensive experiments on real-world datasets demonstrate that TEA+ outperforms the state-of-the-art algorithm by more than four times on most benchmark datasets in terms of computational time when achieving the same clustering quality, and in particular, is an order of magnitude faster on large graphs including the widely studied Twitter and Friendster datasets.

Tables

Table 1. Table 1. Theoretical guarantee of our solution against that of the state-of-the-art solutions.

Algorithm	Accuracy Guarantee	Time Complexity
$\mathrm{𝖢𝗅𝗎𝗌𝗍𝖾𝗋𝖧𝖪𝖯𝖱}$ (chung2014hklocal, )	with probability at least $1-ϵ$,	$O\left(\frac{t\log \left(n\right)}{{ϵ}^3}\right)$
$\mathrm{𝖧𝖪}\text{-}\mathrm{𝖱𝖾𝗅𝖺𝗑}$ (kloster2014heat, )		$O\left(\frac{t{e}^t\log \left(1/{ϵ}_a\right)}{{ϵ}_a}\right)$
Our solutions	with probability at least $1-p_f$,	$O\left(\frac{t\log \left(n/p_f\right)}{{ϵ}_r^2\cdot \delta }\right)$

Table 2. Table 2. Frequently used notations.

Notation	Description
$G$=$(V,E)$	An undirected graph with node set $V$ and edge set $E$
$n,m$	The numbers of nodes and edges in $G$, respectively
$N\left(v\right)$	The set of neighbors of node $v$
$d\left(v\right)$	The degree of node $v$
$\overline{d}$	The average degree of the graph, i.e., $\frac{2m}{n}$
$𝐀,𝐃,𝐏$	The adjacency, diagonal degree, and transition matrices of $G$
$t$	The heat constant of HKPR
$\eta \left(k\right),\psi \left(k\right)$	See Equation (1) and Equation (3), respectively
${𝝆}_s\left[v\right]$	HKPR of $v$ w.r.t. $s$, defined by Equation (2)
${ϵ}_r,\delta ,p_f$	Parameters of an approximate HKPR, as in Section 3
${𝐫}_s^{\left(k\right)}\left[v\right]$	The $k$-hop residue of $v$ during performing push operations from $s$
${𝐪}_s\left[v\right]$	The reserve of $v$ during performing push operations from $s$
$K$	The maximum number of hops during performing push operations from the seed node

Table 3. Table 3. An example for 𝖳𝖤𝖠 + limit-from 𝖳𝖤𝖠 \mathsf{TEA+}

	$k\le 2$	$k=3$	$k=4$
${\sum }_{v\in V}{𝐫}_s^{\left(k\right)}\left[v\right]$	$0$	$0.1$	$0.3$
${\max }_{v\in V}\frac{{𝐫}_s^{\left(k\right)}\left[v\right]}{d\left(v\right)}$	$0$	$\frac{{𝐫}_s^{\left(k\right)}\left[v_1\right]}{d\left(v_1\right)}=0.0025$	$\frac{{𝐫}_s^{\left(k\right)}\left[v_2\right]}{d\left(v_2\right)}=0.0076$
${\max }_{v\in V}{𝐫}_s^{\left(k\right)}\left[v\right]$	$0$	${𝐫}_s^{\left(k\right)}\left[v_1\right]=0.0025$	${𝐫}_s^{\left(k\right)}\left[v_2\right]=0.076$

Table 4. Table 4. The reserves and residues after the first round of push operations.

	$𝒔$	${𝒗}_{\mathrm{𝟏}}$	${𝒗}_{\mathrm{𝟐}}$	${𝒗}_{\mathrm{𝟑}}$	${𝒗}_{\mathrm{𝟒}}$	${𝒗}_{\mathrm{𝟓}}$	${𝒗}_{\mathrm{𝟔}}$	${𝒗}_{\mathrm{𝟕}}$
${𝐪}_s$	$\frac{1}{e^3}$	0	0	0	0	0	0	0
${𝐫}_s^{\left(1\right)}$	0	$\frac{e^3-1}{2e^3}$	$\frac{e^3-1}{2e^3}$	0	0	0	0	0

Table 5. Table 5. The reserves and residues after the second round of push operations.

	$𝒔$	${𝒗}_{\mathrm{𝟏}}$	${𝒗}_{\mathrm{𝟐}}$	${𝒗}_{\mathrm{𝟑}}$	${𝒗}_{\mathrm{𝟒}}$	${𝒗}_{\mathrm{𝟓}}$	${𝒗}_{\mathrm{𝟔}}$	${𝒗}_{\mathrm{𝟕}}$
${𝐪}_s$	$\frac{1}{e^3}$	$\frac{3}{2e^3}$	0	0	0	0	0	0
${𝐫}_s^{\left(1\right)}$	0	0	$\frac{e^3-1}{2e^3}$	0	0	0	0	0
${𝐫}_s^{\left(2\right)}$	$\frac{\tau }{6}$	0	$\frac{\tau }{6}$	$\frac{\tau }{6}$	0	0	0	0

Table 6. Table 6. The reserves and residues after the third round of push operations.

	$𝒔$	${𝒗}_{\mathrm{𝟏}}$	${𝒗}_{\mathrm{𝟐}}$	${𝒗}_{\mathrm{𝟑}}$	${𝒗}_{\mathrm{𝟒}}$	${𝒗}_{\mathrm{𝟓}}$	${𝒗}_{\mathrm{𝟔}}$	${𝒗}_{\mathrm{𝟕}}$
${𝐪}_s$	$\frac{1}{e^3}$	$\frac{3}{2e^3}$	$\frac{3}{2e^3}$	0	0	0	0	0
${𝐫}_s^{\left(1\right)}$	0	0	0	0	0	0	0	0
${𝐫}_s^{\left(2\right)}$	$\frac{\tau }{4}$	$\frac{\tau }{12}$	$\frac{\tau }{6}$	$\frac{\tau }{6}$	$\frac{\tau }{12}$	$\frac{\tau }{12}$	$\frac{\tau }{12}$	$\frac{\tau }{12}$

Table 7. Table 7. Statistics of graph datasets.

Dataset	$𝒏$	$𝒎$	$\overline{𝒅}$
DBLP	317,080	1,049,866	6.62
Youtube	1,134,890	2,987,624	5.27
PLC	2,000,000	9,999,961	9.99
Orkut	3,072,441	117,185,083	76.28
LiveJournal	3,997,962	34,681,189	17.35
3D-grid	9,938,375	29,676,450	5.97
Twitter	41,652,231	1,202,513,046	57.74
Friendster	65,608,366	1,806,067,135	55.06

Table 8. Table 8. The result of evaluating all algorithms on finding real-world communities.

Data	$F_1$-measure					Running Time (ms)
	$\mathrm{𝖢𝗅𝗎𝗌𝗍𝖾𝗋𝖧𝖪𝖯𝖱}$	$\mathrm{𝖬𝗈𝗇𝗍𝖾}\text{-}\mathrm{𝖢𝖺𝗋𝗅𝗈}$	$\mathrm{𝖧𝖪}\text{-}\mathrm{𝖱𝖾𝗅𝖺𝗑}$	$\mathrm{𝖳𝖤𝖠}$	$\mathrm{𝖳𝖤𝖠}+$	$\mathrm{𝖢𝗅𝗎𝗌𝗍𝖾𝗋𝖧𝖪𝖯𝖱}$	$\mathrm{𝖬𝗈𝗇𝗍𝖾}\text{-}\mathrm{𝖢𝖺𝗋𝗅𝗈}$	$\mathrm{𝖧𝖪}\text{-}\mathrm{𝖱𝖾𝗅𝖺𝗑}$	$\mathrm{𝖳𝖤𝖠}$	$\mathrm{𝖳𝖤𝖠}+$
DBLP	0.13655	0.13631	0.13592	0.13679	0.136699	3053.95	2891.64	297.78	176	109.66
Youtube	0.10113	0.10097	0.09858	0.10133	0.10334	7.76	7.2	8.11	2.49	2
LiveJournal	0.64644	0.65105	0.64516	0.64959	0.67	1.34665	1.2	3.57	0.55	0.29
Orkut	0.18497	0.18464	0.19375	0.19333	0.19636	29.95	29.35	62.17	24.78	14.86

Equations

$$\Phi (S)=\frac{|\mathrm{cut}(S)|}{\min \{(\mathrm{vol}(S),\mathrm{vol}(V\setminus S))\}},$$

$$\eta (k)=\frac{e^{-t}{t}^k}{k!}.$$

$${\mathbf{\rho }}_s[v]=\sum _{k=0}^{\infty }\eta (k)\cdot {\mathbf{P}}^k[s,v].$$

$$\left|\frac{{\hat{\mathbf{\rho }}}_s[v]}{d(v)}-\frac{{\mathbf{\rho }}_s[v]}{d(v)}\right|\le {ϵ}_r\cdot \frac{{\mathbf{\rho }}_s[v]}{d(v)};$$

$$\left|\frac{{\hat{\mathbf{\rho }}}_s[v]}{d(v)}-\frac{{\mathbf{\rho }}_s[v]}{d(v)}\right|\le {ϵ}_r\cdot \delta .$$

$$1-n\cdot \exp \left(-\frac{n_r\cdot {ϵ}_r^2\cdot \delta }{2(1+{ϵ}_r/3)}\right).$$

$$\psi (k)=\sum _{\ell =k}^{\infty }\eta (\ell ),$$

$${\mathbf{\rho }}_s[v]={\mathbf{q}}_s[v]+{\sum }_{u\in V}{\sum }_{k=0}^K{\mathbf{r}}_s^{(k)}[u]\cdot {\mathbf{h}}_u^{(k)}[v],$$

$${\mathbf{h}}_u^{(k)}[v]={\sum }_{\ell =0}^{\infty }\frac{\eta (k+\ell )}{\psi (k)}\cdot {\mathbf{P}}^{\ell }[u,v],$$

$$\omega =\frac{2(1+{ϵ}_r/3)\log (1/p_f^{\prime })}{{ϵ}_r^2\cdot \delta },$$

$$p_f^{\prime }=\{\begin{array}{ll}{p}_f,& \text{if }{\sum }_{v\in V}{p}_f^{d(v)-1}\le 1\\ \frac{p_f}{{\sum }_{v\in V}{p}_f^{d(v)-1}},& \text{otherwise}.\end{array}$$

$${\mathbf{\rho }}_s[v]={\mathbf{q}}_s[v]+{\mathbf{a}}_s[v],$$

$${\mathbf{a}}_s[v]=\alpha \cdot {\sum }_{k=0}^K{\sum }_{u\in V}\frac{{\mathbf{r}}_s^{(k)}[u]}{\alpha }\cdot {\mathbf{h}}_u^{(k)}[v].$$

$$\mathbb{E}[Y_i]={\sum }_{u\in V}{\sum }_{k=0}^K\frac{{\mathbf{r}}_s^{(k)}[u]}{\alpha }\cdot {\mathbf{h}}_u^{(k)}[v].$$

$$\mathbb{E}\left[{\sum }_{i=1}^{n_r}{Y}_i\cdot \frac{\alpha }{n_r}\right]={\mathbf{a}}_s[v],$$

$$\mathbb{P}[|X-\mathbb{E}[X]|\ge \lambda ]\le \exp \left(-\frac{{\lambda }^2}{2\mathbb{E}[X]+2\lambda /3}\right).$$

$$O\left(\frac{1}{r_{max}}+\alpha \cdot \omega t\right)=O\left(\frac{t\log (1/p_f^{\prime })}{{ϵ}_r^2\cdot \delta }\right)=O\left(\frac{t\log (n/p_f)}{{ϵ}_r^2\cdot \delta }\right),$$

$$\sum _{\ell =0}^K\underset{v\in V}{\max }\left\{\frac{{\mathbf{r}}_s^{(\ell )}[v]}{d(v)}\right\}\le {ϵ}_a,$$

$$\left|\frac{{\mathbf{q}}_s[v]}{d(v)}-\frac{{\mathbf{\rho }}_s[v]}{d(v)}\right|\le {ϵ}_a.$$

$$K=c\cdot \frac{\log (\frac{1}{{ϵ}_r\cdot \delta })}{\log \left(\bar{d}\right)},$$

$${\sum }_{\ell =0}^4{\max }_{v\in V}\left\{\frac{{\mathbf{r}}_s^{(\ell )}[v]}{d(v)}\right\}=0.0025+0.0076=0.0101,$$

$$\alpha ={\sum }_{k=0}^K{\sum }_{u\in V}{\mathbf{r}}_s^{(k)}[u]=\mathrm{0.4.}$$

$${\beta }_3\cdot {ϵ}_r\delta \cdot d(v)=0.0025\cdot d(v),$$

$${\beta }_4\cdot {ϵ}_r\delta \cdot d(v)=0.0075\cdot d(v).$$

$$\omega =\frac{8(1+{ϵ}_r/6)\log (1/p_f^{\prime })}{{ϵ}_r^2\cdot \delta },$$

$${\beta }_k=\frac{{\sum }_{u\in V}{\mathbf{r}}_s^{(k)}[u]}{{\sum }_{\ell =0}^K{\sum }_{u\in V}{\mathbf{r}}_s^{(\ell )}[u]}.$$

$$\frac{{\mathbf{r}}_s^{(2)}[s]}{d(s)}=\frac{\tau }{8},\frac{{\mathbf{r}}_s^{(2)}[v_1]}{d(v_1)}=\frac{\tau }{36},\frac{{\mathbf{r}}_s^{(2)}[v_2]}{d(v_2)}=\frac{\tau }{36},\frac{{\mathbf{r}}_s^{(2)}[v_3]}{d(v_3)}=\frac{\tau }{6},$$

$$\frac{{\mathbf{r}}_s^{(2)}[v_4]}{d(v_4)}=\frac{{\mathbf{r}}_s^{(2)}[v_5]}{d(v_5)}=\frac{{\mathbf{r}}_s^{(2)}[v_6]}{d(v_6)}=\frac{{\mathbf{r}}_s^{(2)}[v_6]}{d(v_6)}=\frac{\tau }{12}.$$

$${\mathbf{r}}_s^{(2)}[v_3]={\mathbf{r}}_s^{(2)}[v_3]-{ϵ}_r\cdot \delta =\frac{\tau }{18},$$

$${\mathbf{r}}_s^{(2)}[s]={\mathbf{r}}_s^{(2)}[s]-2\cdot {ϵ}_r\delta =\frac{\tau }{36},$$

$$n_r=\alpha \cdot \omega =\sum _{k=0}^2\sum _{u\in V}{\mathbf{r}}_s^{(k)}[u]\cdot \omega =\frac{\tau }{12}\cdot \omega \approx 81.$$

Full text

In this paper, we present $\mathsf{TEA}$ and $\mathsf{TEA+}$, two novel local graph clustering algorithms based on HKPR, to address the aforementioned limitations. Specifically, these algorithms provide non-trivial theoretical guarantees in relative error of HKPR values and the time complexity. The basic idea is to utilize deterministic graph traversal to produce a rough estimation of exact HKPR vector, and then exploit Monte-Carlo random walks to refine the results in an optimized and non-trivial way. In particular, $\mathsf{TEA+}$ offers practical efficiency and effectiveness due to non-trivial optimizations. Extensive experiments on real-world datasets demonstrate that $\mathsf{TEA+}$ outperforms the state-of-the-art algorithm by more than four times on most benchmark datasets in terms of computational time when achieving the same clustering quality, and in particular, is an order of magnitude faster on large graphs including the widely studied Twitter and Friendster datasets.

heat kernel PageRank; local clustering

††journalyear: 2019††copyright: acmcopyright††conference: 2019 International Conference on Management of Data; June 30-July 5, 2019; Amsterdam, Netherlands††booktitle: 2019 International Conference on Management of Data (SIGMOD ’19), June 30-July 5, 2019, Amsterdam, Netherlands††price: 15.00††doi: 10.1145/3299869.3319886††isbn: 978-1-4503-5643-5/19/06††ccs: Mathematics of computing Graph algorithms

1. Introduction

Graph clustering is a fundamental problem that finds numerous applications, e.g., community detection  (fortunato2010community, ; leskovec2010empirical, ; wang2015community, ), image segmentation (felzenszwalb2004efficient, ; tolliver2006graph, ), and protein grouping in biological networks (liao2009isorankn, ; voevodski2009finding, ). The problem has been studied extensively in the literature, and yet, clustering massive graphs remains a challenge in terms of computation efficiency. This motivates a series of algorithms (spielman2004nearly, ; andersen2006local, ; andersen2009finding, ; gharan2012approximating, ; spielman2013local, ; zhu2013local, ; chung2014hklocal, ; orecchia2014flow, ; kloster2014heat, ; avron2015community, ; veldt2016simple, ; wang2017capacity, ) for local clustering, which takes as input an undirected graph $G$ and a seed node $s$, and identifies a cluster (i.e., a set of nodes) containing $s$ in time depending on the size of the cluster, regardless of the size of $G$.

Local clustering algorithms have the potential to underpin interactive exploration of massive graphs. Specifically, they can facilitate exploration of a relatively small region of a large graph that is of interest to a user. For example, consider Bob, a budding entrepreneur, who is interested in exploring the local clusters of visionary entrepreneurs in Twitter. Particularly, he wishes to begin his exploration with the cluster associated with Elon Musk (i.e., seed). Since Bob thinks that Elon Musk is an inspirational entrepreneur, he would like to explore if there are any other notable entrepreneurs (e.g., Kevin Rose) in Elon’s local cluster (e.g., followers, followees) and wishes to further explore the local neighborhoods of these entrepreneurs. In order to ensure a palatable and non-disruptive interactive experience, Bob needs an efficient local clustering framework that can return high quality clusters within few seconds. Which existing local clustering framework can Bob utilize for his exploration?

Spielman and Teng (spielman2004nearly, ; spielman2013local, ) are the first to study the local clustering problem, and they propose a random-walk-based algorithm, $\mathsf{Nibble}$, that optimizes the conductance (bollobas1998modern, ) of the cluster returned. Specifically, the conductance of cluster $S$ is defined as $\Phi(S)=\frac{|\mathrm{cut}(S)|}{\min\{\mathrm{vol}(S),2m-\mathrm{vol}(S)\}}$, where $\mathrm{vol}(S)$ is the sum of the degrees of all nodes in $S$, $m$ is the number of edges in the graph $G$, and $|\mathrm{cut}(C)|$ is number of edges with one endpoint in $S$ and another not in $S$. Intuitively, if a cluster $C$ has a small conductance, then the nodes in $S$ are better connected to each other than to the nodes apart from $S$, in which case $S$ should be considered a good cluster. This algorithm is subsequently improved in a series of work (andersen2006local, ; andersen2009finding, ; gharan2012approximating, ; spielman2013local, ; zhu2013local, ; chung2014hklocal, ; orecchia2014flow, ; wang2017capacity, ; veldt2016simple, ) that aims to either improve the efficiency of local clustering or reduce the conductance of the cluster returned.

The state-of-the-art solutions (chung2014hklocal, ; kloster2014heat, ) for local clustering are based on heat kernel PageRank (HKPR) (chung2007heat, ), which is a measure of the proximity of nodes in graphs. Given a seed node $s$, these solutions first compute a vector ${\widehat{\boldsymbol{\rho}}_{s}}$ where each element ${\widehat{\boldsymbol{\rho}}_{s}}[v]$ approximates the HKPR value of a node $v$ with respect to $s$ (i.e., ${\widehat{\boldsymbol{\rho}}_{s}}[v]$ approximately measures the proximity of $s$ to $v$). Then, they utilize ${\widehat{\boldsymbol{\rho}}_{s}}$ to derive a local cluster $C$ containing $s$. It is shown that the quality of $S$ depends on the accuracy of ${\widehat{\boldsymbol{\rho}}_{s}}$ (kloster2014heat, ; chung2015computing, ), in the sense that the conductance of $S$ tends to decrease with the approximation error in ${\widehat{\boldsymbol{\rho}}_{s}}$. Therefore, existing HKPR-based solutions (chung2014hklocal, ; kloster2014heat, ) all focus on striking a good trade-off between time efficiency and the accuracy of ${\widehat{\boldsymbol{\rho}}_{s}}$. In particular, the current best solution $\mathsf{HK{\hbox{-}}Relax}$ (kloster2014heat, ) ensures that (i) $\frac{1}{d(v)}\left|{\widehat{\boldsymbol{\rho}}_{s}}[v]-{\boldsymbol{\rho}}_{s}[v]\right|<\epsilon_{a}$ for any node $v$, where $\epsilon_{a}$ is a given threshold, $d(v)$ is the degree of node $v$ and ${\boldsymbol{\rho}}_{s}[v]$ is the exact HKPR value of node $v$ with respect to $s$, and (ii) ${\widehat{\boldsymbol{\rho}}_{s}}$ is computed in $O\left(\frac{te^{t}\log{(1/\epsilon_{a})}}{\epsilon_{a}}\right)$ time, where $t$ is constant (referred to as the heat constant) used in the definition of HKPR.

Motivation. The time complexity of $\mathsf{HK{\hbox{-}}Relax}$ has a large factor $e^{t}$, where $t$ (i.e., the heat constant) could be as large as a few dozens (chung2014hklocal, ; kloster2014heat, ; shun2016parallel, ). Consequently, it can be inefficient for several applications. For instance, reconsider Bob’s endeavor to explore the local clusters of Elon Musk and Kevin Rose. $\mathsf{HK{\hbox{-}}Relax}$ consumes around 15s and 48s, respectively, to compute their local clusters. Such performance is disruptive for any interactive graph exploration. In addition, $\mathsf{HK{\hbox{-}}Relax}$ provides an accuracy guarantee on each $\frac{{\widehat{\boldsymbol{\rho}}_{s}}[v]}{d(v)}$ in terms of its absolute error, but as we discuss in Section 3, this guarantee is less than ideal for accurate local clustering. Specifically, HKPR-based local clustering requires ranking each node $v$ in descending order of $\frac{{\widehat{\boldsymbol{\rho}}_{s}}[v]}{d(v)}$, which we refer to as $v$’s normalized HKPR. To optimize this accuracy of this ranking, we observe that it is more effective to minimize the relative errors of normalized HKPR values than their absolute errors. To explain, we note that the normalized HKPR varies significantly from nodes to nodes. For the aforementioned ranking, nodes with large normalized HKPR could tolerate more absolute errors than nodes with small normalized HKPR, and hence, imposing the same absolute error guarantees on all nodes tend to produce sub-optimal results.

Our contributions. Motivated by the deficiency of existing solutions, we present an in-depth study on HKPR-based local clustering, and make the following contributions. First, we formalize the problem of approximate HKPR computation with probabilistic relative error guarantees, and pinpoint why none of the existing techniques could provide an efficient solution to this problem.

Second, based on our problem formulation, we propose two new algorithms, $\mathsf{TEA}$ and $\mathsf{TEA+}$, both of which (i) take as input a seed node $s$, two thresholds $\epsilon_{r},\delta$, and a failure probability $p_{f}$, and (ii) return an approximate HKPR vector ${\widehat{\boldsymbol{\rho}}_{s}}$ where each element ${\widehat{\boldsymbol{\rho}}_{s}}[v]$ with $\frac{{\boldsymbol{\rho}}_{s}[v]}{d(v)}>\delta$ has at most $\epsilon_{r}$ relative error with at least $1-p_{f}$ probability (i.e., all significant HKPR values are accurately approximated with high probability). The core of $\mathsf{TEA}$ is an adaptive method that combines deterministic graph traversal with random walks to estimate normalized HKPR in a cost-effective manner, while $\mathsf{TEA+}$ significantly improves over $\mathsf{TEA}$ in terms of practical efficiency by incorporating a number of non-trivial optimization techniques. Both algorithms have a time complexity of $O\left(\frac{t\log{(n/p_{f})}}{\epsilon^{2}_{r}\cdot\delta}\right)$, which eliminates the exponential term $e^{t}$ in $\mathsf{HK{\hbox{-}}Relax}$’s efficiency bound (see Table 1).

Third, we experimentally evaluate them against the state of the art, using large datasets with up to 65 million nodes and 1.8 billion edges. Our results show that $\mathsf{TEA+}$ is up to an order of magnitude faster than competing methods when achieving the same clustering quality. In particular, it can compute the local clusters of Elon Musk and Kevin Rose within 1.3s and 6.1s, respectively, thereby facilitating interactive exploration.

Paper Organization. The rest of the paper is organized as follows. In Section 2, we introduce background on heat kernel-based local clustering. An overview of our solution framework is presented in Section 3. We present $\mathsf{TEA}$ and $\mathsf{TEA+}$ in Sections 4 and 5, respectively. Related work is reviewed in Section 6. We evaluate the practical efficiency of our algorithms against the competitors in Section 7. Finally, Section 8 concludes the paper. Proofs of theorems and lemmas are given in Appendix B. Table 2 lists the notations that are frequently used in our paper.

2. Preliminaries

2.1. Basic Terminology

Let $G=(V,E)$ be an undirected and unweighted graph, where $V$ and $E$ denote the node and edge sets, respectively. We use $d(v)$ to denote the degree of node $v$, and $\mathbf{A}$ to denote the adjacency matrix of $G$; i.e., $\mathbf{A}[i,j]=\mathbf{A}[j,i]=1$ if and only if $(v_{i},v_{j})\in E$. Let $\boldsymbol{\mathrm{D}}$ be the diagonal degree matrix of $G$, where $\mathbf{D}[i,i]=d(v_{i})$. Then, the probability transition matrix (a.k.a. random walk transition matrix) for $G$ is defined as $\boldsymbol{\mathrm{P}}=\boldsymbol{\mathrm{D}}^{-1}\boldsymbol{\mathrm{A}}$. Accordingly, $\mathbf{P}^{k}[s,v]$ denotes the probability that a $k$-hop ($k\geq 1$) random walk from node $s$ would end at $v$.

A cluster in $G$ is a node set $S\subseteq V$. Intuitively, a good cluster should be both internally cohesive and well separated from the remainder of $G$. We say that $S$ is a high-quality cluster if it has a small conductance (bollobas1998modern, ) $\Phi(S)$, defined as:

[TABLE]

where $\mathrm{vol}(S)$ is the volume of $S$, namely, the sum of the degrees of all nodes in $S$, and $\mathrm{cut}(S)$ is the cut of $S$, i.e., the set of edges with one endpoint in $S$ and another not in $S$.

2.2. Heat Kernel-based Local Clustering

Given a heat constant $t$ and two nodes $u$ and $v$, the HKPR value from $u$ to $v$ is defined as the probability that a random walk of length $k$ starting from $u$ would end at $v$, where $k$ is sampled from the following Poisson distribution:

[TABLE]

Let $s$ be the seed node for local clustering. We define the HKPR vector ${\boldsymbol{\rho}}_{s}$ of $s$ as an $n$-size vector, such that the $i$-th element of ${\boldsymbol{\rho}}_{s}$ equals the HKPR value from $s$ to the $i$-th node in $G$. In addition, we use ${\boldsymbol{\rho}}_{s}[v]$ to denote the HKPR value from $s$ to $v$, which is defined by

[TABLE]

Existing heat-kernel-based algorithms (chung2014hklocal, ; kloster2014heat, ; chung2015distributed, ; shun2016parallel, ) all adopt a two-phase approach. In particular, they first compute an approximate HKPR vector ${\widehat{\boldsymbol{\rho}}_{s}}$ for $s$, and then perform a sweep as follows:

1. (1)

Take the set $S^{*}$ of nodes with non-zero values in $\widehat{\boldsymbol{\rho}}_{s}$. 2. (2)

Sort the nodes $v\in S^{*}$ in descending order of $\frac{{\widehat{\boldsymbol{\rho}}_{s}}[v]}{d(v)}$. Let $S^{*}_{i}$ be a set containing the first $i$ nodes in the sorted sequence. 3. (3)

Inspect $S^{*}_{i}$ in ascending order of $i$. Return the set $S^{*}_{i}$ with the smallest conductance among the ones that have been inspected.

It is shown in (zhu2013local, ; shun2016parallel, ) that the above sweep can be conducted in ${O\left(|S^{*}|\cdot\log{|S^{*}|}\right)}$ time, assuming that $\widehat{\boldsymbol{\rho}}_{s}$ is given in a sparse representation with $O(|S^{*}|)$ entries. In contrast, the computation of $\widehat{\boldsymbol{\rho}}_{s}$ is much more costly, and hence, has been the main subject of research in existing work (chung2014hklocal, ; kloster2014heat, ; chung2015distributed, ; shun2016parallel, ).

3. Solution Overview

Our solution for local clustering is based on heat kernel PageRank, and it follows the same two-phase framework in the existing work (chung2007heat, ; chung2009local, ; chung2014hklocal, ; kloster2014heat, ; chung2015distributed, ; shun2016parallel, ). That is, we also compute an approximate HKPR vector ${\widehat{\boldsymbol{\rho}}_{s}}$ for $s$, and then conduct a sweep on ${\widehat{\boldsymbol{\rho}}_{s}}$. However, we require that ${\widehat{\boldsymbol{\rho}}_{s}}$ should be a $(d,\epsilon_{r},\delta)$-approximate HKPR vector, which is a criterion not considered in any existing work (chung2007heat, ; chung2009local, ; chung2014hklocal, ; kloster2014heat, ; chung2015distributed, ; shun2016parallel, ).

Definition 0.

($(d,\epsilon_{r},\delta)$-approximate HKPR) Let ${\boldsymbol{\rho}}_{s}$ be the HKPR vector for a seed node $s$, and ${\widehat{\boldsymbol{\rho}}_{s}}$ be an approximated version of ${\boldsymbol{\rho}}_{s}$. ${\widehat{\boldsymbol{\rho}}_{s}}$ is $(d,\epsilon_{r},\delta)$-approximate if it satisfies the following conditions:

• •

For every $v\in V$ with $\frac{{\boldsymbol{\rho}}_{s}[v]}{d(v)}>\delta$,

[TABLE]

• •

For every $v\in V$ with $\frac{{\boldsymbol{\rho}}_{s}[v]}{d(v)}\leq\delta$,

[TABLE]

In other words, we require $\frac{{\widehat{\boldsymbol{\rho}}_{s}}[v]}{d(v)}$ to provide a relative error guarantee when $\frac{{\boldsymbol{\rho}}_{s}[v]}{d(v)}>\delta$, and an absolute error guarantee when $\frac{{\boldsymbol{\rho}}_{s}[v]}{d(v)}\leq\delta$. This is to ensure that when we sort the nodes in descending order or $\frac{{\widehat{\boldsymbol{\rho}}_{s}}[v]}{d(v)}$ (which is a crucial step in the sweep for local clustering), the sorted sequence would be close to the one generated based on $\frac{{\boldsymbol{\rho}}_{s}[v]}{d(v)}$. We do not consider relative error guarantees when $\frac{{\boldsymbol{\rho}}_{s}[v]}{d(v)}\leq\delta$, because (i) ensuring a small relative error for such a node $v$ requires an extremely accurate estimation of its normalized HKPR, which would incur significant computation overheads, and (ii) providing such high accuracy for $v$ is unnecessary, since $v$’s tiny normalized HKPR value indicates that it is not relevant to the result of local clustering.

By the definition of HKPR (in Equation (2)), the HKPR value of $v$ w.r.t. $s$ is a weighted sum of $k$-hop random walk transition probabilities from $s$ to $v$, where $k$ is a Poisson distributed length. Thus, a straightforward method to compute $(d,\epsilon_{r},\delta)$-approximate HKPR for seed node $s$ is to conduct Monte-Carlo simulations using a large number of random walks. Specifically, each random walk should start from $s$, and should have a length that is sampled from the Poisson distribution in Equation (1). Let $n_{r}$ be the total number of random walks, and ${\widehat{\boldsymbol{\rho}}_{s}}[v]$ be the fraction of walks that end at a node $v$. Then, we can use ${\widehat{\boldsymbol{\rho}}_{s}}[v]$ as an estimation of ${\boldsymbol{\rho}}_{s}[v]$. By the Chernoff bound (see Lemma 5) and union bound, it can be shown that ${\widehat{\boldsymbol{\rho}}_{s}}$ is $(d,\epsilon_{r},\delta)$-approximate with probability at least

[TABLE]

Therefore, if we are to ensure that the above probability is at least $1-p_{f}$, then we can set $n_{r}=\frac{2(1+\epsilon_{r}/3)\log{(n/p_{f})}}{\epsilon^{2}_{r}\cdot\delta}$. In that case, the time required to generate the random walks is $O\left(\frac{t\log{(n/p_{f})}}{\epsilon^{2}_{r}\cdot\delta}\right)$. The main issue of this Monte-Carlo approach, however, is that it incurs considerable overheads in practice (see our experimental results in Section 7.4). To explain, consider a node $v$ with a small ${\boldsymbol{\rho}}_{s}[v]$ but a relatively large $\frac{{\boldsymbol{\rho}}_{s}[v]}{d(v)}>\delta$. By the requirements of $(d,\epsilon_{r},\delta)$-approximation, we need to ensure that $|{\widehat{\boldsymbol{\rho}}_{s}}[v]-{\boldsymbol{\rho}}_{s}[v]|\leq\epsilon_{r}\cdot{\boldsymbol{\rho}}_{s}[v]$. In turn, this requires that the number $n_{r}$ of random walks should be large; otherwise, the number of walks that end at $v$ would be rather small, in which case the estimation of ${\boldsymbol{\rho}}_{s}[v]$ would be inaccurate.

We also note that none of the existing methods (chung2014hklocal, ; kloster2014heat, ) can be adopted to compute $(d,\epsilon_{r},\delta)$-approximate HKPR efficiently. In particular, $\mathsf{HK{\hbox{-}}Relax}$ (kloster2014heat, ) only ensures that for any node $v\in V$, $\frac{1}{d(v)}\left|{\widehat{\boldsymbol{\rho}}_{s}}[v]-{\boldsymbol{\rho}}_{s}[v]\right|<\epsilon_{a}.$ If we are to use $\mathsf{HK{\hbox{-}}Relax}$ for $(d,\epsilon_{r},\delta)$-approximation, then we need to set $\epsilon_{a}=\epsilon_{r}\cdot\delta$, in which case its complexity would be $O\left(\frac{te^{t}\log{\left(\frac{1}{\epsilon_{r}\cdot\delta}\right)}}{\epsilon_{r}\cdot\delta}\right)$, which is considerably worse than the time complexity of the Monte-Carlo approach, due to the exponential term $e^{t}$. The $\mathsf{ClusterHKPR}$ (chung2014hklocal, ) algorithm suffers from a similar issue, as we point out in Section 6.

To mitigate the deficiencies of the aforementioned methods, we present (in Section 4 and Section 5) two more efficient HKPR algorithms that satisfy the following requirements:

1. (1)

Return a $(d,\epsilon_{r},\delta)$-approximate HKPR vector ${\widehat{\boldsymbol{\rho}}_{s}}$ with at least $1-p_{f}$ probability, where $p_{f}$ is a user-specified parameter; 2. (2)

Run in $O\left(\frac{t\log(n/p_{f})}{\epsilon^{2}_{r}\cdot\delta}\right)$ expected time.

4. The $\boldsymbol{\mathsf{TEA}}$ Algorithm

This section presents $\mathsf{TEA}$ 111Two-Phase Heat Kernel Approximation, our first-cut solution for $(d,\epsilon_{r},\delta)$-approximate HKPR. $\mathsf{TEA}$ is motivated by the inefficiency of the Monte-Carlo approach explained in Section 3, i.e., it requires a large number of random walks to accurately estimate HKPR values. To address this issue, we propose to combine $\mathsf{Monte\textrm{-}Carlo}$ with a secondary algorithm that could help reduce the number of random walks needed. In particular, we first utilize the secondary algorithm to efficiently compute a rough estimation ${\mathbf{q}_{s}}[v]$ of ${\boldsymbol{\rho}}_{s}[v]$, and then perform random walks to refine each ${\mathbf{q}_{s}}[v]$, so as to transform ${\mathbf{q}_{s}}$ into a $(d,\epsilon_{r},\delta)$-approximate HKPR vector ${\widehat{\boldsymbol{\rho}}_{s}}$. Towards this end, there are several issues that we need to address:

1. (1)

How to design a secondary algorithm that could generate a rough approximation of the HKPR vector at a small computation cost? 2. (2)

How to enable $\mathsf{Monte\textrm{-}Carlo}$ to leverage the output of the secondary algorithm for improved efficiency? 3. (3)

How to ensure that the combination of $\mathsf{Monte\textrm{-}Carlo}$ and the secondary algorithm still provides strong theoretical assurance in terms of time complexity and accuracy?

To answer the above questions, we first present our choice of the secondary algorithm, referred to as $\mathsf{HK{\hbox{-}}Push}$, in Section 4.1; after that, we elaborate the integration of $\mathsf{HK{\hbox{-}}Push}$ and $\mathsf{Monte\textrm{-}Carlo}$ in Section 4.2, and then provide a theoretical analysis of the combined algorithm in Section 4.3.

4.1. $\boldsymbol{\mathsf{HK{\hbox{-}}Push}}$

Algorithm 1 shows the pseudo-code of our secondary algorithm, $\mathsf{HK{\hbox{-}}Push}$, for deriving a rough approximation ${\widehat{\boldsymbol{\rho}}_{s}}$ of the HKPR vector. Its basic idea is to begin with a vector ${\widehat{\boldsymbol{\rho}}_{s}}$ where ${\widehat{\boldsymbol{\rho}}_{s}}[s]=1$ and ${\widehat{\boldsymbol{\rho}}_{s}}[v]=0$ for all nodes $v$ except $s$, and then perform a traversal of $G$ starting from $s$, and keep refining ${\widehat{\boldsymbol{\rho}}_{s}}$ during the course of the traversal. In addition, to facilitate its combination with random walks, it not only returns an approximate HKPR vector ${\mathbf{q}_{s}}$, but also outputs $K+1$ auxiliary vectors $\mathbf{r}^{(0)}_{s},\ldots,\mathbf{r}^{(K)}_{s}\in\mathbb{R}^{n}$ that could be used to guide the random walks conducted by $\mathsf{Monte\textrm{-}Carlo}$. We refer to ${\mathbf{q}_{s}}$ as the reserve vector and $\mathbf{r}^{(k)}_{s}$ as the $k$-hop residue vector. Accordingly, for any node $v$, ${\mathbf{q}_{s}}[v]$ and $\mathbf{r}^{(k)}_{s}[v]$ are referred to as the reserve and $k$-hop residue of $v$, respectively.

More specifically, $\mathsf{HK{\hbox{-}}Push}$ takes as input $G$, $s$, and a residue threshold $r_{max}$. It begins by setting all entries in ${\mathbf{q}_{s}}$ and ${\mathbf{r}^{(k)}_{s}}$ to zero, except that $\mathbf{r}^{(0)}_{s}[s]=1$ (Lines 1-2). After that, it starts an iterative process to traverse $G$ from $s$ (Lines 3-7). In particular, in each iteration, it inspects the $k$-hop residue vectors to identify a node $v$ whose $k$-hop residue ${\mathbf{r}^{(k)}_{s}}[v]$ is above $r_{max}\cdot d(v)$. If such a node $v$ exists, then the algorithm updates the reserve and $k$-hop residue of $v$, as well as the $(k+1)$-hop residues of $v$’s neighbors. In particular, it first adds $\frac{\eta(k)}{\psi(k)}$ fraction of $v$’s $k$-hop residue ${\mathbf{r}^{(k)}_{s}}[v]$ to its reserve ${\mathbf{q}_{s}}[v]$, where $\eta(k)$ is as defined in Equation (1) and

[TABLE]

and then evenly distribute the other $1-\frac{\eta(k)}{\psi(k)}$ fraction to the $(k+1)$-hop residues of $v$’s neighbors (Lines 4-6). For convenience, we refer to the operation of distributing a fraction of $v$’s $k$-hop residue to one of its neighbors as a push operation. The rationale of the aforementioned push operations is that, if a random walk from $s$ arrives at $v$ at the $k$-th hop, then it has $\frac{\eta(k)}{\psi(k)}$ probability to terminate at $v$, and has $1-\frac{\eta(k)}{\psi(k)}$ probability to traverse to a randomly selected neighbor of $v$ at the next hop. After that, the algorithm sets ${\mathbf{r}^{(k)}_{s}}[v]=0$ (Line 7), and proceeds to the next iteration. After the iterative process terminates, $\mathsf{HK{\hbox{-}}Push}$ identifies the largest $K$ such that $\mathbf{r}^{(K)}_{s}$ has at least one non-zero entry, and returns ${\mathbf{q}_{s}}$ and $\mathbf{r}^{(0)}_{s},\ldots,\mathbf{r}^{(K)}_{s}$. The following lemma shows a crucial property of these reserve and residue vectors:

Lemma 0.

Consider any iteration in Algorithm 1. Let ${\mathbf{q}_{s}}$ and $\mathbf{r}^{(0)}_{s},\ldots\mathbf{r}^{(K)}_{s}$ be the reserve and residue vectors constructed by the end of the iteration. We have

[TABLE]

where

[TABLE]

i.e., $\mathbf{h}^{(k)}_{u}[v]$ is the probability that a random walk stops at $v$, conditioned on the $k$-hop of the walk is at $u$. $\square$

Intuitively, Lemma 1 indicates that for any node $v$, ${\mathbf{q}_{s}}[v]$ is a lower bound of ${\boldsymbol{\rho}}_{s}[v]$ in any iteration in Algorithm 1. Since each iteration of Algorithm 1 only increases the reserve of a selected node and never decreases any others, it guarantees that the difference between ${\mathbf{q}_{s}}$ and ${\boldsymbol{\rho}}_{s}$ monotonically decreases, i.e., ${\mathbf{q}_{s}}$ becomes a better approximation of ${\boldsymbol{\rho}}_{s}$ as the algorithm progresses. Although $\mathsf{HK{\hbox{-}}Push}$ may produce results that are far from satisfying the requirements of $(d,\epsilon_{r},\delta)$-approximate HKPR, it is sufficient for the integration with $\mathsf{Monte\textrm{-}Carlo}$, as we show in Section 4.2.

4.2. Algorithm

Basic Idea. To incorporate $\mathsf{HK{\hbox{-}}Push}$ into $\mathsf{Monte\textrm{-}Carlo}$, we utilize Equation (4), which shows that the exact HKPR vector ${\boldsymbol{\rho}}_{s}$ can be expressed as a function of ${\mathbf{q}_{s}}$, ${\mathbf{r}^{(k)}_{s}}$, and $\mathbf{h}^{(k)}_{u}[v]$ for any $u,v\in V$, and $k\in[0,K]$. Recall that ${\mathbf{q}_{s}}$ and ${\mathbf{r}^{(k)}_{s}}$ are outputs of $\mathsf{HK{\hbox{-}}Push}$, while $\mathbf{h}^{(k)}_{u}[v]$ is the conditional probability that a random walk terminates at node $v$ given that its $k$-th hop is at node $u$. If we can accurately estimate $\mathbf{h}^{(k)}_{u}[v]$ for any given $u$, $v$, and $k$, then we can combine the estimated values with ${\mathbf{q}_{s}}$ and ${\mathbf{r}^{(k)}_{s}}$ to obtain an approximate version of ${\boldsymbol{\rho}}_{s}$. Towards this end, we conduct random walks starting from $u$, and estimate $\mathbf{h}^{(k)}_{u}[v]$ based on the frequency that $v$ appears at the $k$-th hop of the random walks. Algorithm 2 shows the pseudo-code of our random walk generation method, referred to as $k\textrm{-}\mathsf{RandomWalk}$.

The following lemma proves that $k\textrm{-}\mathsf{RandomWalk}$ samples each node $v$ with probability $\mathbf{h}^{(k)}_{u}[v]$.

Lemma 0.

Given $G$, $u$, and $k$, for any node $v$, Algorithm 2 returns $v$ with probability $\mathbf{h}^{(k)}_{u}[v]$. $\square$

Details. Algorithm 3 illustrates the pseudo-code of $\mathsf{TEA}$, our first-cut solution that (i) incorporates $\mathsf{HK{\hbox{-}}Push}$ and $k\textrm{-}\mathsf{RandomWalk}$ and (ii) computes a $(d,\epsilon_{r},\delta)$-approximate HKPR vector ${\widehat{\boldsymbol{\rho}}_{s}}$ with at least $1-p_{f}$ probability for any given seed node $s$. Given $G$, $\epsilon_{r}$, $\delta$, $r_{max}$, and failure probability $p_{f}$ as inputs, the algorithm starts by invoking $\mathsf{HK{\hbox{-}}Push}$ with three parameters: $G$, $s$, and $r_{max}$ (Line 6). It returns an approximate HKPR vector ${\widehat{\boldsymbol{\rho}}_{s}}$ and $K+1$ residue vectors $\mathbf{r}^{(0)}_{s},\ldots,\mathbf{r}^{(K)}_{s}$ from $\mathsf{HK{\hbox{-}}Push}$. Then, $\mathsf{TEA}$ proceeds to refine ${\widehat{\boldsymbol{\rho}}_{s}}$ using $k\textrm{-}\mathsf{RandomWalk}$ (Lines 7-13). In particular, $\mathsf{TEA}$ first computes the sum $\alpha$ of the residues in $\mathbf{r}^{(0)}_{s},\ldots,\mathbf{r}^{(K)}_{s}$ (Line 7), and computes

[TABLE]

where

[TABLE]

Note that $p^{\prime}_{f}$ can be pre-computed when the graph $G$ is loaded. Given $\omega$, $\mathsf{TEA}$ performs $n_{r}=\alpha\cdot\omega$ random walks using $k\textrm{-}\mathsf{RandomWalk}$ (Lines 9-12), such that the starting point $u$ of each walk is sampled with probability $\frac{{\mathbf{r}^{(k)}_{s}}[u]}{\alpha}$ (Line 10). Note that this sampling procedure can be conducted efficiently by conducting an alias structure (walker1974new, ) on the non-zero elements in $\mathbf{r}^{(0)}_{s},\ldots,\mathbf{r}^{(K)}_{s}$. For each random walk that ends at a node $v$, $\mathsf{TEA}$ increases ${\widehat{\boldsymbol{\rho}}_{s}}[v]$ by $\frac{\alpha}{n_{r}}$ (Line 12).

Observe that the parameter $r_{max}$ in $\mathsf{TEA}$ controls the balance between its two main components: $\mathsf{HK{\hbox{-}}Push}$ and $k\textrm{-}\mathsf{RandomWalk}$. In particular, by Algorithm 1, the entries in $\mathbf{r}^{(0)}_{s},\ldots,\mathbf{r}^{(K)}_{s}$ are bounded by $r_{max}$. Therefore, when $r_{max}$ is small, $\alpha=\sum_{k=0}^{K}{\sum_{u\in V}{{\mathbf{r}^{(k)}_{s}}[u]}}$ would decrease, in which case the total number $\alpha\cdot\omega$ of random walks conducted by $\mathsf{TEA}$ would also be small. As a trade-off, the processing cost of $\mathsf{HK{\hbox{-}}Push}$ would increase, as shown in the following lemma:

Lemma 0.

Given residue threshold $r_{max}$, Algorithm 1 runs in $O\left(\frac{1}{r_{max}}\right)$ time and requires $O\left(\frac{1}{r_{max}}\right)$ space (excluding the space required by the input graph). In addition, in the residue vectors $\mathbf{r}^{(0)}_{s},\ldots,\mathbf{r}^{(K)}_{s}$ returned by Algorithm 1, there are $O\left(\frac{1}{r_{max}}\right)$ non-zero elements in total.

To strike a balance between the costs incurred by $\mathsf{HK{\hbox{-}}Push}$ and $k\textrm{-}\mathsf{RandomWalk}$, we set $r_{max}=O(\frac{1}{\omega\cdot t})$. In that case, the processing cost of $\mathsf{HK{\hbox{-}}Push}$ is $O(\omega\cdot t)$, while $k\textrm{-}\mathsf{RandomWalk}$ incurs $O(\alpha\omega t)$ expected cost, due to the following lemma:

Lemma 0.

The expected cost of each invocation of $k\textrm{-}\mathsf{RandomWalk}$ is $O(t)$.

Hence, setting $r_{max}=O(\frac{1}{\omega\cdot t})$ ensures that the overheads of $\mathsf{HK{\hbox{-}}Push}$ and $k\textrm{-}\mathsf{RandomWalk}$ are roughly comparable.

4.3. Analysis

Correctness. Let ${\mathbf{q}_{s}}$ denote the approximate HKPR vector obtained from $\mathsf{HK{\hbox{-}}Push}$ in Line 6 of $\mathsf{TEA}$, and ${\widehat{\boldsymbol{\rho}}_{s}}$ be the approximate HKPR vector eventually output by $\mathsf{TEA}$. In the following, we show that ${\widehat{\boldsymbol{\rho}}_{s}}$ is a $(d,\epsilon_{r},\delta)$-approximate HKPR vector.

First, by Lemma 1, we have the following equation for any node $v$:

[TABLE]

where

[TABLE]

Consider the $i$-th invocation of $k\textrm{-}\mathsf{RandomWalk}$ in $\mathsf{TEA}$. Let $(u,k)$ be the entry sampled by $\mathsf{TEA}$ (in Line 10) before the invocation, and $v$ be the node returned by $k\textrm{-}\mathsf{RandomWalk}$. Let $Y_{i}$ be a Bernoulli variable that equals $1$ if $v$ is returned, and [math] otherwise. By Lemma 2,

[TABLE]

Combining Equations (8) and (9), we have

[TABLE]

which indicates that $\sum_{i=1}^{n_{r}}{Y_{i}\cdot\frac{\alpha}{n_{r}}}$ is an unbiased estimator of $\mathbf{a}_{s}[v]$. By the Chernoff bound (in Lemma 5), we can prove that this estimator is highly accurate, based on which we obtain Theorem 6.

Lemma 0 (Chernoff Bound (chung2006concentration, )).

Let $X_{1},X_{2},\cdots,X_{n_{r}}\in[0,1]$ be i.i.d. random variables, and $X=\sum_{i=1}^{n_{r}}{X_{i}}$. Then,

[TABLE]

Theorem 6.

$\mathsf{TEA}$* outputs a $(d,\epsilon_{r},\delta)$-approximate HKPR vector ${\widehat{\boldsymbol{\rho}}_{s}}$ with probability at least $1-p_{f}$.*

Time and Space Complexities. Given $r_{max}=O\left(\frac{1}{\omega\cdot t}\right)$, $\mathsf{HK{\hbox{-}}Push}$ runs in $O(\omega\cdot t)$ time and $O(\omega\cdot t)$ space. In addition, the computation of $\alpha$ as well as the construction of alias structure on ${\mathbf{r}^{(k)}_{s}}$ can be done in time and space linear to the total number of non-zero entries in $\mathbf{r}^{(0)}_{s},\ldots,\mathbf{r}^{(K)}_{s}$, which is $O(\omega\cdot t)$ (see Lemma 3). Furthermore, in terms of both space and time, the total expected cost incurred by the random walks in $\mathsf{TEA}$ is $O(\alpha\omega t)$, where $\alpha<1$. Therefore, the time complexity of $\mathsf{TEA}$ is

[TABLE]

and its space complexity is $O\left(n+m+\frac{t\log{(n/p_{f})}}{\epsilon^{2}_{r}\cdot\delta}\right)$, where the $n+m$ term is due to storing of the input graph.

5. The TEA+ Algorithm

Although $\mathsf{TEA}$ provides a strong accuracy guarantee, we observe in our experiments that it often performs a large number of random walks, which degrades its computation efficiency. One may attempt to reduce the cost of random walks by decreasing the residue threshold $r_{max}$ in $\mathsf{TEA}$ (see the discussion in the end of Section 4.2), but this cost reduction would be offset by the fact that $\mathsf{HK{\hbox{-}}Push}$ incurs a larger overhead when $r_{max}$ is small.

In this section, we present $\mathsf{TEA+}$, an algorithm that significantly improves over $\mathsf{TEA}$ in terms of practical efficiency without degrading its theoretical guarantees. $\mathsf{TEA+}$ is similar in spirit to $\mathsf{TEA}$ in that it combines random walks with a variant of $\mathsf{HK{\hbox{-}}Push}$, but there is a crucial difference: after $\mathsf{TEA+}$ obtains the residue vectors $\mathbf{r}^{(0)}_{s},\ldots,\mathbf{r}^{(K)}_{s}$, it may reduce some entries in the residue vectors before performing random walks. That is, for each node $u$ that has a non-zero entry $\mathbf{r}^{(k)}_{s}[u]$ for some $k$, $\mathsf{TEA+}$ may choose to perform a smaller number of random walks from $u$ than $\mathsf{TEA}$ does, which effectively reduces the total cost of random walks. Establishing the correctness of this pruning approach, however, is non-trivial. In what follows, we first discuss in Section 5.1 the extreme case where we can derive $(d,\epsilon_{r},\delta)$-approximate HKPR while ignoring all elements in the residue vectors $\mathbf{r}^{(0)}_{s},\ldots,\mathbf{r}^{(K)}_{s}$; after that, in Section 5.2, we generalize our discussions to the case where we reduce the non-zero entries in the residue vectors instead of completely omitting them.

5.1. The Case without Random Walks

Suppose that we are to let $\mathsf{TEA}$ achieve $(d,\epsilon_{r},\delta)$-approximation without performing any random walks. In that case, we would need to ensure that Line 6 of $\mathsf{TEA}$ obtains a $(d,\epsilon_{r},\delta)$-approximate HKPR vector from $\mathsf{HK{\hbox{-}}Push}$. Towards this end, we present the following theorem:

Theorem 1.

Let ${\mathbf{q}_{s}}$ and $\mathbf{r}^{(0)}_{s},\ldots\mathbf{r}^{(K)}_{s}$ be the reserve and residue vectors returned by $\mathsf{HK{\hbox{-}}Push}$. If

[TABLE]

then, for any node $v$ in $G$, we have

[TABLE]

Theorem 1 provides a sufficient condition (i.e., Inequality (11)) for $\mathsf{HK{\hbox{-}}Push}$ to ensure $\epsilon_{a}$ absolute error in each $\frac{{\mathbf{q}_{s}}[v]}{d(v)}$. By Definition 1, such ${\mathbf{q}_{s}}$ is a $(d,\epsilon_{r},\delta)$-approximate HKPR vector as long as $\epsilon_{a}\leq\epsilon_{r}\cdot\delta$. That said, it is rather inefficient to let $\mathsf{HK{\hbox{-}}Push}$ run until Inequality (11) is satisfied. Instead, we propose to let $\mathsf{HK{\hbox{-}}Push}$ run with a fixed budget of processing cost. If it is able to satisfy Inequality (11) with $\epsilon_{a}=\epsilon_{r}\cdot\delta$ before the budget is depleted, then we return ${\mathbf{q}_{s}}$ as the final result; otherwise, we proceed to refine ${\mathbf{q}_{s}}$ using random walks (see Section 5.2).

Based on the above discussion, we present $\mathsf{HK{\hbox{-}}Push+}$ (in Algorithm 4), which is a revised version of $\mathsf{HK{\hbox{-}}Push}$ with three major changes. First, $\mathsf{HK{\hbox{-}}Push+}$’s input parameters include three thresholds $\epsilon_{r}$, $\delta$, and $n_{p}$, and it has two new termination conditions (in Line 6): (i) Inequality (11) holds with $\epsilon_{a}=\epsilon_{r}\cdot\delta$; (ii) The number of push operations that it performs reaches $n_{p}$. Recall that a push operation refers to the operation of converting part of a node’s $k$-hop residue to one of its neighbor’s $(k+1)$-hop residue. In other words, $\mathsf{HK{\hbox{-}}Push+}$ runs in $O(n_{p})$ time and requires $O(n_{p})$ space, and it returns a $(d,\epsilon_{r},\delta)$-approximate HKPR vector whenever Inequality (11) is satisfied.

Second, $\mathsf{HK{\hbox{-}}Push+}$ judiciously performs push operations only on nodes $v$ with residue ${\mathbf{r}^{(k)}_{s}}[v]>\frac{\epsilon_{r}\cdot\delta}{K}\cdot d(v)$ (Line 4), whereas $\mathsf{HK{\hbox{-}}Push}$ conducts push operations only when ${\mathbf{r}^{(k)}_{s}}[v]$ is larger than an input given threshold $r_{max}\cdot d(v)$. The rationale is that $\mathsf{HK{\hbox{-}}Push+}$ strives to reduce the $k$-hop residue of each node below $\frac{\epsilon_{r}\cdot\delta}{K}\cdot d(v)$, so as to satisfy Inequality (11); in contrast, $\mathsf{HK{\hbox{-}}Push}$ is not guided by Inequality (11), and hence, uses an ad hoc threshold $r_{max}$ instead.

Third, $\mathsf{HK{\hbox{-}}Push+}$ makes the maximum number $K$ of hops be specified as an input parameter, whereas $\mathsf{HK{\hbox{-}}Push}$ does not fix $K$ in advance. We use a fixed $K$ in $\mathsf{HK{\hbox{-}}Push+}$ because (i) as mentioned, Line 4 of $\mathsf{HK{\hbox{-}}Push+}$ requires testing whether there exists a node $v$ with ${\mathbf{r}^{(k)}_{s}}[v]>\frac{\epsilon_{r}\cdot\delta}{K}\cdot d(v)$, and (ii) such a test can be efficiently implemented when $K$ is fixed. Otherwise, whenever $K$ changes, we would need to recheck all nodes’ residues to see if ${\mathbf{r}^{(k)}_{s}}[v]>\frac{\epsilon_{r}\cdot\delta}{K}\cdot d(v)$ holds, which would incur considerable overheads. Meanwhile, $\mathsf{HK{\hbox{-}}Push}$ can afford to let $K$ dynamically change, since it uses a fixed residue threshold $r_{max}$ given as input. In our implementation of $\mathsf{HK{\hbox{-}}Push+}$, we set

[TABLE]

where $\bar{d}$ is the average degree of the nodes in $G$, and $c$ is a constant that we decide based on our experiments in Section 7.2. We refer interested readers to Appendix A for the rationale of this setting of $K$.

5.2. The Case with Random Walks

Suppose that $\mathsf{HK{\hbox{-}}Push+}$ depletes its computation budget $n_{p}$ before it satisfies Inequality (11) with $\epsilon_{a}=\epsilon_{r}\cdot\delta$. In that case, the HKPR vector ${\mathbf{q}_{s}}$ returned by $\mathsf{HK{\hbox{-}}Push+}$ does not ensure $(d,\epsilon_{r},\delta)$-approximation, and we need to refine ${\mathbf{q}_{s}}$ by conducting random walks according to the residue vectors $\mathbf{r}^{(0)}_{s},\ldots,\mathbf{r}^{(K)}_{s}$ returned by $\mathsf{HK{\hbox{-}}Push+}$. To reduce the number of random walks required, we propose to reduce the residues in $\mathbf{r}^{(0)}_{s},\ldots,\mathbf{r}^{(K)}_{s}$ when conducting random walks, based on the following intuition.

First, the reduction of residues would incur some errors in the approximate HKPR vector, as demonstrated in Theorem 1 (where we ignore all residues in $\mathbf{r}^{(k)}_{s}$ and using ${\mathbf{q}_{s}}$ directly as the final approximate HKPR vector). Second, if we only reduce the residues in $\mathbf{r}^{(k)}_{s}$ by a small value, then the absolute errors incurred by the reduction could be so small that they would not jeopardize $(d,\epsilon_{r},\delta)$-approximation. In particular, as we show in Section 5.5, if we reduce every residue $\mathbf{r}^{(k)}_{s}[v]$ by $\beta_{k}\cdot\epsilon_{r}\delta\cdot d(v)$ ($k=0,1,\ldots,K$ and $v\in V$), then the absolute error in $\frac{{\widehat{\boldsymbol{\rho}}_{s}}[v]}{d(v)}$ incurred by the residue reduction is at most $\epsilon_{r}\delta\cdot\sum_{k=0}^{K}\beta_{k}$. In other words, if we choose $\beta_{k}$ such that $\sum_{k=0}^{K}\beta_{k}=1$, then the resulting absolute error in $\frac{{\widehat{\boldsymbol{\rho}}_{s}}[v]}{d(v)}$ is at most $\epsilon_{r}\delta$, which is permissible under $(d,\epsilon_{r},\delta)$-approximation. The following example demonstrates the benefit of this residue reduction method.

Example 0.

Suppose that given a graph $G$, seed node $s$, $K=4$, and $\epsilon_{r}\cdot\delta=0.01$, $\mathsf{HK{\hbox{-}}Push+}$ returns residue vectors $\mathbf{r}^{(0)}_{s},\ldots,\mathbf{r}^{(4)}_{s}$ that have the characteristics in Table 3. Note that $v_{1}$ has the largest $3$-hop residue and degree-normalized $3$-hop residue, while $v_{2}$ has the largest $4$-hop residue and degree-normalized $4$-hop residue. In addition, all of the $3$-hop residues except $v_{1}$’s are less than $0.0025$, and all of the $4$-hop residues except $v_{2}$’s are less than $0.075$. We can observe that

[TABLE]

which is slightly larger than $\epsilon_{r}\cdot\delta$. In this case, according to Lines 7-8 in $\mathsf{TEA}$, we need to perform $\alpha\cdot\omega$ random walks, where

[TABLE]

Now suppose that we apply the residue reduction method, setting $\beta_{3}=1/4$ and $\beta_{4}=3/4$. In that case, we reduce every residue $\mathbf{r}^{(3)}_{s}[v]$ by

[TABLE]

and every residue $\mathbf{r}^{(4)}_{s}[v]$ by

[TABLE]

Then, all residues in $\mathbf{r}^{(3)}_{s}$ and $\mathbf{r}^{(4)}_{s}$ are reduced to [math], except that $\mathbf{r}^{(4)}_{s}[v_{2}]$ is decreased to $0.076-\beta_{4}\cdot\epsilon_{r}\delta\cdot d(v_{2})=0.001$. Accordingly, $\alpha=\sum_{\ell=0}^{4}{\max_{v\in V}{\mathbf{r}^{(\ell)}_{s}[v]}}$ is reduced from 0.4 to 0.001, which implies that the number of random walks required is reduced by 400 times. $\square$

5.3. Details of $\boldsymbol{\mathsf{TEA+}}$

Based on the ideas described in Sections 5.1 and  5.2, we present $\mathsf{TEA+}$, which utilizes $\mathsf{HK{\hbox{-}}Push+}$ and random walks to compute a $(d,\epsilon_{r},\delta)$-approximate HKPR vector ${\widehat{\boldsymbol{\rho}}_{s}}$ with at least $1-p_{f}$ probability for any given seed node $s$. Algorithm 5 illustrates the pseudo-code of $\mathsf{TEA+}$. The input parameters of $\mathsf{TEA+}$ are identical to those of $\mathsf{TEA}$, except that $\mathsf{TEA+}$ takes an additional parameter $c$, which, as mentioned in Section 5.1, is used to decide the maximum number $K$ of hops used in $\mathsf{HK{\hbox{-}}Push+}$.

$\mathsf{TEA+}$ starts by invoking $\mathsf{HK{\hbox{-}}Push+}$ with the following parameters (Lines 5-6): $G$, $\epsilon_{r}$, $\delta$, $K=c\cdot\frac{\log{(\frac{1}{\epsilon_{r}\cdot\delta})}}{\log{(\bar{d})}}$, and $n_{p}=\frac{\omega\cdot t}{2}$, where

[TABLE]

and $p^{\prime}_{f}$ is as defined in Equation (6) and is pre-computed when $G$ is loaded. Then, $\mathsf{HK{\hbox{-}}Push+}$ returns an approximate HKPR vector ${\widehat{\boldsymbol{\rho}}_{s}}$ and $K+1$ residue vectors $\mathbf{r}^{(0)}_{s},\ldots,\mathbf{r}^{(K)}_{s}$.

If $\sum_{k=0}^{K}{\max_{u\in V}{\left\{\frac{{\mathbf{r}^{(k)}_{s}}[u]}{d(u)}\right\}}}\leq\epsilon_{r}\cdot\delta$, then by Theorem 1, ${\widehat{\boldsymbol{\rho}}_{s}}$ is a $(d,\epsilon_{r},\delta)$-approximate HKPR vector. In that case, $\mathsf{TEA+}$ immediately terminates and returns ${\widehat{\boldsymbol{\rho}}_{s}}$ (Line 7). Otherwise, $\mathsf{TEA+}$ proceeds to refine ${\widehat{\boldsymbol{\rho}}_{s}}$ using $k\textrm{-}\mathsf{RandomWalk}$ (Lines 8-20). But before that, $\mathsf{TEA+}$ first applies the residue reduction method discussed in Section 5.2. Specifically, it decreases each residue value ${\mathbf{r}^{(k)}_{s}}[u]$ by $\beta_{k}\cdot\epsilon_{r}\delta\cdot d(u)$ (Lines 8-11), where

[TABLE]

The rationale of this choice of $\beta_{k}$ is as follows. First, $\sum_{k=0}^{K}\beta_{k}=1$, which is crucial for $(d,\epsilon_{r},\delta)$-approximation, as we mention in Section 5.2. Second, we set $\beta_{k}$ to be proportional to $\sum_{u\in V}{{\mathbf{r}^{(k)}_{s}}[u]}$ because, intuitively, when $\sum_{u\in V}{{\mathbf{r}^{(k)}_{s}}[u]}$ is large, the residue values in ${\mathbf{r}^{(k)}_{s}}$ also tend to be large, in which case we need a larger reduction of the residues in ${\mathbf{r}^{(k)}_{s}}$ to effectively reduce the number of random walks needed.

After reducing the residues in $\mathbf{r}^{(0)}_{s},\ldots,\mathbf{r}^{(K)}_{s}$, $\mathsf{TEA+}$ performs random walks according to the reduced residues, in the same way as $\mathsf{TEA}$ does (Lines 12-17). This results in a refined approximate HKPR vector ${\widehat{\boldsymbol{\rho}}_{s}}$. Then, $\mathsf{TEA+}$ gives ${\widehat{\boldsymbol{\rho}}_{s}}[v]$ a final touch by adding $\frac{\epsilon_{r}\cdot\delta}{2}\cdot d(v)$ to each ${\widehat{\boldsymbol{\rho}}_{s}}[v]$ (Lines 18-19). The intuition of adding this offset to each ${\widehat{\boldsymbol{\rho}}_{s}}[v]$ is as follows. The residue reduction method leads to an underestimation of each HKPR value, and amount of underestimation is in $[0,\epsilon_{r}\cdot\delta\cdot d(v)]$. By adding an offset $\frac{\epsilon_{r}\cdot\delta}{2}\cdot d(v)$ to ${\widehat{\boldsymbol{\rho}}_{s}}[v]$, the range of the error in ${\widehat{\boldsymbol{\rho}}_{s}}[v]$ becomes $[-\frac{\epsilon_{r}\cdot\delta}{2}\cdot d(v),\frac{\epsilon_{r}\cdot\delta}{2}\cdot d(v)]$, in which case the maximum absolute error in ${\widehat{\boldsymbol{\rho}}_{s}}[v]$ is reduced by half, which help tightening the accuracy bound of $\mathsf{TEA+}$.

Note that Lines 18-19 in $\mathsf{TEA+}$ can be performed in $O(1)$ time, as we can keep each ${\widehat{\boldsymbol{\rho}}_{s}}[v]$ unchanged but record the value of $\frac{\epsilon_{r}\cdot\delta}{2}$ along with ${\widehat{\boldsymbol{\rho}}_{s}}$. Then, whenever ${\widehat{\boldsymbol{\rho}}_{s}}[v]$ is accessed, we can add $\frac{\epsilon_{r}\cdot\delta}{2}\cdot d(v)$ to ${\widehat{\boldsymbol{\rho}}_{s}}[v]$ on the fly. In addition, for the purpose of local clustering, we can simply ignore this offset of $\frac{\epsilon_{r}\cdot\delta}{2}\cdot d(v)$ since it does not affect the ranking of nodes based on $\frac{{\widehat{\boldsymbol{\rho}}_{s}}[v]}{d(v)}$.

5.4. Example for $\boldsymbol{\mathsf{TEA+}}$

We use the graph $G^{\prime}$ in Figure 1 to illustrate $\mathsf{TEA+}$. For ease of presentation, we set $t=3$, $p_{f}=10^{-2}$, $\epsilon_{r}=0.5$, $\delta=\frac{2\tau}{9}$, and $\tau=1-\frac{4}{e^{3}}\approx 0.8$, and we set $c$ to a constant such that $K=2$. Consider that we invoke Algorithm 5 on $G^{\prime}$ with seed node $s$ and above parameters as inputs. $\mathsf{TEA+}$ first computes $\omega\approx\frac{970}{\tau}$ and $n_{p}\approx\frac{1455}{\tau}$, and then invokes Algorithm 4 with seed node $s$ and parameters $\epsilon_{r}$, $\delta$, $K$, and $n_{p}$ as inputs.

Initially, Algorithm 4 sets $\mathbf{r}^{(0)}_{s}[s]=1$ and $k=0$. In that case, $k<K$, and $\frac{\mathbf{r}^{(0)}_{s}[s]}{d(s)}=0.5>\frac{\epsilon_{r}\delta}{K}\approx 0.0445$. Accordingly, Algorithm 4 performs push operations from $s$. It first converts $\frac{\eta(0)}{\psi(0)}=\frac{1}{e^{3}}$ portion of $s$’s reserve $\mathbf{r}^{(0)}_{s}[s]$ into its reserve $\mathbf{q}_{s}[s]$, and then distributes $\frac{1}{2}(1-\frac{\eta(0)}{\psi(0)})=\frac{e^{3}-1}{2e^{3}}$ portion to each of its two neighbors, i.e., $v_{1}$ and $v_{2}$. Table 4 shows the reserves and residues after the first round of push operations.

As shown in Table 4 and Figure 1, $\frac{\mathbf{r}^{(1)}_{s}[v_{1}]}{d(v_{1})}\approx 0.1584$ and $\frac{\mathbf{r}^{(1)}_{s}[v_{2}]}{d(v_{2})}\approx 0.0792$, both of which are greater than $\frac{\epsilon_{r}\delta}{K}\approx 0.0445$. In addition, $k=1<K$, and the number of push operations performed $i$ is 2, which is smaller than $n_{p}$. Hence, Algorithm 4 starts the second round of push operations on $\mathbf{r}^{(1)}_{s}[v_{1}]$. In that case, it converts $\frac{\eta(1)}{\psi(1)}$ portion of $\mathbf{r}^{(1)}_{s}[v_{1}]$ into $v_{1}$’s reserve, and distributes $\frac{1}{3}(1-\frac{\eta(1)}{\psi(1)})\cdot\mathbf{r}^{(1)}_{s}[v_{1}]=\frac{\tau}{6}$ residue to each of $v_{1}$’s three neighbors $s,v_{2}$ and $v_{3}$. Table 5 lists the reserves and residues after this round of push operations.

In the third round, Algorithm 4 converts $\frac{\eta(1)}{\psi(1)}$ portion of $\mathbf{r}^{(1)}_{s}[v_{2}]$ into the reserve of $v_{2}$. Meanwhile, each of $v_{2}$’s neighbors receives $\frac{1}{6}(1-\frac{\eta(1)}{\psi(1)})\cdot\mathbf{r}^{(1)}_{s}[v_{2}]$ residue. After this round, $k=K=2$, and hence, Algorithm 4 terminates with the reserves in Table 6 as the approximate HKPR values in ${\widehat{\boldsymbol{\rho}}_{s}}$.

Based on Table 5 and Figure 1, we have the following results:

[TABLE]

[TABLE]

Since $\sum_{k=0}^{2}{\max_{u\in V}{\left\{\frac{{\mathbf{r}^{(k)}_{s}}[u]}{d(u)}\right\}}}=\frac{\tau}{6}>\epsilon_{r}\cdot\delta=\frac{\tau}{9}$, Algorithm 5 does not terminate, but would perform random walks to refine ${\widehat{\boldsymbol{\rho}}_{s}}$. It first computes $\beta_{0}=\beta_{1}=0$ and $\beta_{2}=1$, and then reduces the residue of each node $v$ by $\mathbf{r}^{(2)}_{s}[v]=\max\left\{0,\ {\mathbf{r}^{(2)}_{s}[v]}-\beta_{2}\cdot\epsilon_{r}\delta\cdot d(v)\right\}$. This leads to

[TABLE]

and other residues become 0. The number of random walks required by Algorithm 5 is therefore

[TABLE]

Before performing random walks from $s$ and $v_{3}$, Algorithm 5 first builds an alias structure on all nodes, such that entries $(s,2)$ and $(v_{3},2)$ would be sampled with probabilities $\frac{\mathbf{r}^{(2)}_{s}[s]}{\alpha}=\frac{1}{3}$ and $\frac{\mathbf{r}^{(2)}_{s}[v_{3}]}{\alpha}=\frac{2}{3}$, respectively. After that, it invokes Algorithm 2 for $n_{r}=81$ times. In each invocation, it samples the starting entry $(v,2)$ from this alias structure, and at the $\ell$-th step, it terminates at the current node $u$ with probability $\frac{\eta(2+\ell)}{\psi(2+\ell)}$, and traverses to a randomly selected neighbor of $u$ with the other $1-\frac{\eta(2+\ell)}{\psi(2+\ell)}$ probability. Once it terminates at a node $w$, Algorithm 5 increases $w$’s approximate HKPR value ${\widehat{\boldsymbol{\rho}}_{s}}[w]$ by $\frac{\alpha}{n_{r}}=\frac{\tau}{970}$. $\mathsf{TEA+}$ returns ${\widehat{\boldsymbol{\rho}}_{s}}$ as an approximation of the exact HKPR vector ${\boldsymbol{\rho}}_{s}$ after 81 random walks.

Finally, we perform a sweep over ${\widehat{\boldsymbol{\rho}}_{s}}$ to find a cut with the best conductance. That is, we sort the nodes $v$ by $\frac{{\widehat{\boldsymbol{\rho}}_{s}}[v]}{d(v)}$ in descending order. If the sorted order is $s,v_{1},v_{3},v_{2},v_{4}-v_{7}$, then the sweep would produce $\{s,v_{1},v_{3}\}$ as the best cut, as its conductance $\frac{1}{3}$ is smaller than other cuts encountered during the sweep.

5.5. Analysis

Correctness.

Let ${\mathbf{q}_{s}}$ denote the approximate HKPR vector returned by $\mathsf{HK{\hbox{-}}Push+}$ in Line 6 of $\mathsf{TEA+}$, $\mathbf{r}^{(0)}_{s},\ldots,\mathbf{r}^{(K)}_{s}$ be the residue vectors output by $\mathsf{HK{\hbox{-}}Push+}$ at the same step, and ${\widehat{\boldsymbol{\rho}}_{s}}$ be the final version of the HKPR vector output by $\mathsf{TEA+}$. We define a residue vector $\mathbf{rb}^{(k)}_{s}$ as:

[TABLE]

Observe that $\mathbf{rb}^{(k)}_{s}[u]$ equals the amount of residue reduction on ${\mathbf{r}^{(k)}_{s}}[u]$ occurred in Lines 8-11 of $\mathsf{TEA+}$.

Similar to the correctness analysis in Section 4.3, by Lemma 1, we have the following equation for any node $v$:

[TABLE]

where

[TABLE]

Then, the approximation error in each ${\widehat{\boldsymbol{\rho}}_{s}}[v]$ can be regarded as the sum of two approximation errors for $\mathbf{a}_{s}[v]$ and $\mathbf{b}_{s}[v]$, respectively. The error in $\mathbf{a}_{s}[v]$ is caused by sampling errors in $k\textrm{-}\mathsf{RandomWalk}$, and hence, it can be bounded using the Chernoff bound, in a way similar to the analysis in Section 4.3. Meanwhile, the error in $\mathbf{b}_{s}[v]$ is due to the residue reduction procedure in Lines 8-11 of $\mathsf{TEA+}$. In what follows, we present an analysis of the error in $\mathbf{b}_{s}[v]$.

By Equation (13), for any node $u\in V$ and $k\in[0,K]$, the amount of residue reduction on ${\mathbf{r}^{(k)}_{s}}[u]$ satisfies $0\leq\mathbf{rb}^{(k)}_{s}[u]\leq\beta_{k}\cdot\epsilon_{r}\delta\cdot d(u)$. Combining this with Equation (16), we have

[TABLE]

Lemma 0 ((pons2005computing, )).

Let $u$ and $v$ be any two nodes in $G$, and $\mathbf{P}^{k}[u,v]$ (resp. $\mathbf{P}^{k}[v,u]$) be the probability that a length-$k$ random walk from $u$ ends at $v$ (resp. from $v$ ends at $u$). Then,

[TABLE]

By Lemma 3 and the definition of $\mathbf{h}^{(k)}_{u}[v]$ in Equation 5, for any node $v\in V$ and $k\in[0,K]$, we have

[TABLE]

Combining Equations (17) and (5.5), we have

[TABLE]

Therefore, estimating $\mathbf{b}_{s}[v]$ as $\frac{\epsilon_{r}\cdot\delta}{2}\cdot d(v)$ incurs an absolute error of at most $\frac{\epsilon_{r}\cdot\delta}{2}\cdot d(v)$.

Based on the above analysis, we establish the accuracy guarantee of $\mathsf{TEA+}$ as follows:

Theorem 4.

$\mathsf{TEA+}$* outputs a $(d,\epsilon_{r},\delta)$-approximate HKPR vector ${\widehat{\boldsymbol{\rho}}_{s}}$ with probability at least $1-p_{f}$.*

Time and Space Complexities. The time and space complexities of $\mathsf{TEA+}$ depend on its two main components: $\mathsf{HK{\hbox{-}}Push+}$ and $k\textrm{-}\mathsf{RandomWalk}$. As discussed in Section 5.1, both the computation and space overheads of $\mathsf{HK{\hbox{-}}Push+}$ are $O(n_{p})$. Since $\mathsf{TEA+}$ sets $n_{p}=\frac{\omega\cdot t}{2}$, its invocation of $\mathsf{HK{\hbox{-}}Push+}$ incurs $O\left(\frac{t\cdot\log{(n/p_{f})}}{\epsilon^{2}_{r}\cdot\delta}\right)$ time and space costs. Meanwhile, the total number of random walks conducted by $\mathsf{TEA+}$ is no more than that by $\mathsf{TEA}$, and hence, the computational and space costs of generating random walks in $\mathsf{TEA+}$ are both $O\left(\frac{t\cdot\log{(n/p_{f})}}{\epsilon^{2}_{r}\cdot\delta}\right)$ in expectation. Thus, the expected time and space complexities of $\mathsf{TEA+}$ are $O\left(\frac{t\cdot\log{(n/p_{f})}}{\epsilon^{2}_{r}\cdot\delta}\right)$ and $O\left(m+n+\frac{t\cdot\log{(n/p_{f})}}{\epsilon^{2}_{r}\cdot\delta}\right)$, respectively, where the $m+n$ term is due to the space required by the input graph.

6. Related Work

In this section, we first review two HKPR algorithms, $\mathsf{ClusterHKPR}$ and $\mathsf{HK{\hbox{-}}Relax}$, that are most related to our solutions; after that, we review other work related to local clustering and HKPR computation.

$\boldsymbol{\mathsf{ClusterHKPR}}$. $\mathsf{ClusterHKPR}$ (chung2014hklocal, ) is a random-walk-based method for computing approximate HKPR. Given a seed node $s$, it performs $\frac{16\log n}{\epsilon^{3}}$ random walks from $s$, with the constraint that each walk has a length at most $K$; after that, for each node $v$, it uses the fraction of walks that end at $v$ as an estimation ${\widehat{\boldsymbol{\rho}}_{s}}[v]$ of $v$’s HKPR. It is shown that with probability at least $1-\epsilon$, $\mathsf{ClusterHKPR}$ guarantees that

[TABLE]

Note that for the above guarantee to be meaningful, $\epsilon\ll 1$ should hold; otherwise, the successful probability $1-\epsilon$ of $\mathsf{ClusterHKPR}$ would be too small, and there could be too many nodes having a large absolute error up to $\epsilon$. In particular, if we are to ensure that $\mathsf{ClusterHKPR}$ achieves $(d,\epsilon_{r},\delta)$-approximation with probability at least $1-p_{f}$, then we have to set $\epsilon\leq\min\left\{\epsilon_{r}\cdot\delta,p_{f}\right\}$. However, when $\epsilon\ll 1$, $\mathsf{ClusterHKPR}$ incurs excessive computation cost, since its expected time complexity is inversely proportional to $\epsilon^{3}$.

$\boldsymbol{\mathsf{HK{\hbox{-}}Relax}}$. $\mathsf{HK{\hbox{-}}Relax}$ (kloster2014heat, ) is a deterministic algorithm that runs in $O\left(\frac{te^{t}\log{(1/\epsilon_{a})}}{\epsilon_{a}}\right)$ time and returns an approximate HKPR vector ${\widehat{\boldsymbol{\rho}}_{s}}$ such that $\left|\frac{{\widehat{\boldsymbol{\rho}}_{s}}[v]}{d(v)}-\frac{{\boldsymbol{\rho}}_{s}[v]}{d(v)}\right|\leq\epsilon_{a}$ for any $v\in V$. $\mathsf{HK{\hbox{-}}Relax}$ is similar to our $\mathsf{HK{\hbox{-}}Push}$ algorithm in that they both (i) maintain an approximation HKPR vector ${\widehat{\boldsymbol{\rho}}_{s}}$ and a number of residue vectors, and (ii) incrementally refine ${\widehat{\boldsymbol{\rho}}_{s}}$ by applying push operations according to the residue of each node. However, there exist three major differences between $\mathsf{HK{\hbox{-}}Relax}$ and $\mathsf{HK{\hbox{-}}Push}$. First, $\mathsf{HK{\hbox{-}}Relax}$ and $\mathsf{HK{\hbox{-}}Push}$ define the residue of each node in different manners, due to which $\mathsf{HK{\hbox{-}}Relax}$ requires more sophisticated approaches than $\mathsf{HK{\hbox{-}}Push}$ to update ${\widehat{\boldsymbol{\rho}}_{s}}$ and the residue vectors after each push operation. Second, $\mathsf{HK{\hbox{-}}Relax}$ ignores all $k$-hop residues with $k>2t\log{\frac{1}{\epsilon_{a}}}$, whereas $\mathsf{HK{\hbox{-}}Push}$ retains all residues generated for the refinement of ${\widehat{\boldsymbol{\rho}}_{s}}$ via random walks. Third, $\mathsf{HK{\hbox{-}}Relax}$ and $\mathsf{HK{\hbox{-}}Push}$ have different termination conditions. Due to these differences, we are unable combine $\mathsf{HK{\hbox{-}}Relax}$ with random walks to achieve the same performance guarantee provided by our $\mathsf{TEA}$ and $\mathsf{TEA+}$ algorithms.

Other methods for local clustering. The first local graph clustering algorithm, $\mathsf{Nibble}$, is introduced in the seminal work (spielman2004nearly, ; spielman2013local, ) by Spielman and Teng. After that, Anderson et al. (andersen2006local, ) propose $\mathsf{PR{\hbox{-}}Nibble}$, a local clustering algorithm based on personalized PageRank (haveliwala2002topic, ; jeh2003scaling, ), and show that it improves over $\mathsf{Nibble}$ in terms of the theoretical guarantees of both clustering quality and time complexity. In turn, Anderson et al.’s method is improved in subsequent work (andersen2009finding, ; gharan2012approximating, ) based on the volume-biased evolving set process (diaconis1990strong, ). Subsequent work (zhu2013local, ; orecchia2014flow, ; wang2017capacity, ; veldt2016simple, ) achieves further improved guarantees on the quality of local clustering, but the methods proposed are mostly of theoretical interests only, as they are difficult to implement and offer rather poor practical efficiency. As a consequence, $\mathsf{HK{\hbox{-}}Relax}$ remains the state-of-the-art for local clustering in terms of practical performance (kloster2014heat, ; avron2015community, ; chung2015computing, ).

In recent work (shun2016parallel, ), Shun et al. study parallel implementations for $\mathsf{Nibble}$, $\mathsf{PR{\hbox{-}}Nibble}$, $\mathsf{ClusterHKPR}$, and $\mathsf{HK{\hbox{-}}Relax}$, respectively, and are able to achieve significant speedup on a machine with $40$ cores. We believe that our algorithms may also exploit parallelism for higher efficiency, but a thorough treatment of this problem is beyond the scope of this paper.

Methods for personalized PageRank. We note that $\mathsf{TEA}$ and $\mathsf{TEA+}$ are similar in spirit to several recent methods (wang2017fora, ; lofgren2016personalized, ; wei2018topppr, ; wang2016hubppr, ) for computing personalized PageRank (PPR), since they all combine a push-operation-based algorithm with random walks. Hence, at first glance, it may seem that we can simply adopt and extend these techniques to address HKPR computation. Unfortunately, this is not the case as HKPR is inherently more sophisticated than PPR. In particular, even though both HKPR and PPR measure the proximity of a node $v$ with respect to another node $u$ by the probability that a random walk starting from $u$ would end at $v$, they differ significantly in the ways that they define random walks:

• •

PPR’s random walks are Markovian: in each step of a walk, it terminates with a fixed probability $\alpha\in(0,1)$, regardless of the previous steps.

• •

HKPR’s random walks are non-Markovian: the termination probability of a walk at the $i$-th step is a function of $i$, i.e., the walk has to remember the number of steps that it has traversed, so as to decide whether it should terminate.

The Markovianness of PPR random walks is a key property exploited in the methods in (wang2017fora, ; lofgren2016personalized, ; wei2018topppr, ; wang2016hubppr, ). Specifically, when computing the PPR $p(u,v)$ from node $u$ to node $v$, the methods in (lofgren2016personalized, ; wei2018topppr, ; wang2016hubppr, ) require performing a backward search which starts from $v$ and traverses the incoming edges of each node in a backward manner. For each node $w$ encountered and each of its incoming neighbor $x$, the backward search needs to calculate the probability that a random walk hitting $x$ at a certain step would arrive at $w$ at the next step. For PPR random walks, this probability is a constant decided only by $\alpha$ and the number of $x$’s outgoing neighbors. Unfortunately, for HKPR random walks the probability is not a constant, but a variable depending on the number of steps that the walk has taken before reaching $x$. In other words, this variable is not unique even when $w$ and $x$ are fixed, due to which the backward search no longer can be utilized. This issue makes it unpalatable to extend the methods in (lofgren2016personalized, ; wei2018topppr, ; wang2016hubppr, ) to compute HKPR.

Meanwhile, the FORA method in (wang2017fora, ) does not require a backward search; instead, it combines random walks with a forward search from $u$ that is similar to the $\mathsf{HK{\hbox{-}}Push}$ algorithm used in $\mathsf{TEA}$. However, $\mathsf{TEA}$ is more sophisticated than FORA as it needs to account for the non-Markovianness of HKPR, and there are three major differences between the two methods. First, $\mathsf{TEA}$ requires maintaining multiple residue vectors in its invocation of $\mathsf{HK{\hbox{-}}Push}$, since it needs to differentiate the residues generated at different steps of the forward search; otherwise, it would not be able to combine the results of $\mathsf{HK{\hbox{-}}Push}$ with random walks because of HKPR’s non-Markovianness. In contrast, FORA only needs to maintain one residue vector, as the Markovianness of PPR allows it to merge the resides produced at different steps of the forward search. Second, the theoretical analysis of $\mathsf{TEA}$ is more challenging than that of FORA, since it is more complicated to (i) maintain and update multiple residue vectors and (ii) combine random walks with the forward traversal in a way that takes into account the non-Markovianness of HKPR. Third, $\mathsf{TEA}$ provides an accuracy guarantee in terms of each node’s degree-normalized HKPR, whereas FORA’s accuracy guarantee is on each node’s PPR without normalization.

Last but not the least, we note that our $\mathsf{TEA+}$ algorithm, which significantly improves over $\mathsf{TEA}$ in terms of practical efficiency, is based on a new optimization that is specifically designed for HKPR and our notion of $(d,\epsilon_{r},\delta)$-approximation. This optimization is inapplicable for PPR computation, which further differentiates $\mathsf{TEA+}$ from FORA.

7. Experiments

7.1. Experimental Setup

We conduct all experiments on a Linux server with a Intel Xeon(R) E5-2650 v2@2.60GHz CPU and 64GB RAM. For fair comparison, all algorithms are implemented in C++ and compiled by g++ 4.8 with -O3 optimization.

We use 6 undirected real-world graphs and 2 synthetic graphs which are used in recent work (kloster2014heat, ; shun2016parallel, ; chung2015computing, ) as benchmark datasets (Table 7). We obtain DBLP, Youtube, Orkut, LiveJournal, and Friendster from (linksnap, ). PLC is a synthetic graph, and it is generated by Holme and Kim algorithm for generating graphs with powerlaw degree distribution and approximate average clustering. 3D-grid is a synthetic grid graph in 3-dimensional space where every node has six edges, each connecting it to its 2 neighbors in each dimension. Twitter is a symmetrized version of a snapshot of the Twitter network (kwak2010twitter, ). For each dataset, we pick 50 seed nodes uniformly at random as our query sets.

Unless specified otherwise, following previous work (kloster2014heat, ; banerjee2015fast, ), we set heat constant $t=5$. In addition, for all randomized algorithms, we set failure probability $p_{f}=10^{-6}$. We report the average query time (measured in wall-clock time) of each algorithm on each dataset with various parameter settings. Note that the $y$-axis is in log-scale and the measurement unit for running time is millisecond (ms).

7.2. Tuning Parameter $c$ for $\boldsymbol{\mathsf{TEA+}}$

First, we experimentally study how to set the parameter $c$ so as to obtain the best performance for $\mathsf{TEA+}$ in practice. We run $\mathsf{TEA+}$ with parameters $\epsilon_{r}=0.5$, $\delta=\frac{1}{n}$, and varying $c$ from $0.5$ to $5$ on all 8 graphs.

Figure 2 plots the running time of $\mathsf{TEA+}$ on each dataset for different $c$. We omit the results for Twitter and Friendster when $c=0.5$, because it takes several hours to finish one query. We can make the following observations. For each dataset, the running time decreases first as $c$ grows. The reason is that $\mathsf{TEA+}$ degrades to $\mathsf{Monte\textrm{-}Carlo}$ when $c$ is very small, and if we keep increasing $c$, $\mathsf{HK{\hbox{-}}Push+}$ will perform more push operations so as to reduce the number of random walks until $c$ balances the costs incurred for $\mathsf{HK{\hbox{-}}Push+}$ and $k\textrm{-}\mathsf{RandomWalk}$. On the other hand, when $c$ increases further, the overhead of $\mathsf{HK{\hbox{-}}Push+}$ goes up gradually, disrupting the balance between $\mathsf{HK{\hbox{-}}Push+}$ and $k\textrm{-}\mathsf{RandomWalk}$. This leads to higher running time. More specifically, we can see that for graphs with small average degree including DBLP, Youtube, PLC and 3D-grid, the costs are minimized when $c$ is around $2$. On the other hand, for graphs with high average degree (e.g., Orkut, LiveJournal, Twitter, and Friendster), we note that $c=2.5$ achieves the best performance. Based on the above observations, a good value choice for $c$ is $2.5$, when the overheads on most of the graphs are minimized. In the sequel, we set $c=2.5$.

7.3. Comparison between $\boldsymbol{\mathsf{TEA}}$ and $\boldsymbol{\mathsf{TEA+}}$

In this set of experiments, we compare $\mathsf{TEA+}$ with $\mathsf{TEA}$ based on identical theoretical accuracy guarantees. For both $\mathsf{TEA}$ and $\mathsf{TEA+}$, we set the relative error threshold $\epsilon_{r}=0.5$. Since the best values for $r_{max}$ vary largely for different parameter settings and various datasets, we are unable to find a universal optimal value for $r_{max}$. Instead, we tune $r_{max}$ for $\mathsf{TEA}$ with different error thresholds on each dataset separately. That is, we scale $\frac{1}{\omega\cdot t}$ up or down such that the costs for $\mathsf{HK{\hbox{-}}Push}$ and $k\textrm{-}\mathsf{RandomWalk}$ in $\mathsf{TEA}$ are roughly balanced and the total cost is minimized.

Figure 3 reports the computational time of $\mathsf{TEA}$ and $\mathsf{TEA+}$ on all datasets when varying $\epsilon_{r}$ from $0.1$ to $0.9$ and fixing $\delta=10^{-6}$. Observe that $\mathsf{TEA+}$ always outperforms $\mathsf{TEA}$ markedly for all datasets. In particular, when the relative error threshold $\epsilon_{r}$ is large (e.g., $0.5-0.9$), $\mathsf{TEA+}$ is one to two orders of magnitude faster than $\mathsf{TEA}$ on most of datasets, but only around $5$ times faster on 3D-grid. This is caused by the fact that each node in 3D-grid has six neighbors, thus residues will drop below the threshold quickly and both $\mathsf{TEA}$ and $\mathsf{TEA+}$ require very few random walks. As we keep decreasing $\epsilon_{r}$, the gap between $\mathsf{TEA}$ and $\mathsf{TEA+}$ is narrowed. Especially, when $\epsilon_{r}=0.1$, $\mathsf{TEA+}$ achieves around $5\times$ to $10\times$ speedup. The reasons are as follows. When the relative error thresholds are large, $\mathsf{TEA+}$ only needs a small number of push operations and random walks due to the new termination conditions of $\mathsf{HK{\hbox{-}}Push+}$ and residue reduction method compared to $\mathsf{TEA}$. However, when the relative error thresholds are very small, the termination conditions of $\mathsf{HK{\hbox{-}}Push+}$ are harder to satisfy and as a result incurs much higher costs to terminate. Furthermore, the residue reduction method is not able to reduce the number of random walks significantly since the residue sum is already small. Thus, both $\mathsf{TEA}$ and $\mathsf{TEA+}$ perform many push operations and random walks, and as a result, the speedup is modest. The results demonstrate the power of new termination conditions of $\mathsf{HK{\hbox{-}}Push+}$ and residue reducation method in $\mathsf{TEA+}$, especially when the error thresholds are not very small. In summary, $\mathsf{TEA+}$ outperforms $\mathsf{TEA}$ without sacrificing theoretical accuracies of HKPR values.

7.4. Comparisons with Competitors

We compare $\mathsf{TEA}$ and $\mathsf{TEA+}$ against $\mathsf{ClusterHKPR}$, $\mathsf{SimpleLocal}$ (veldt2016simple, ), $\mathsf{CRD}$ (wang2017capacity, ), $\mathsf{Monte\textrm{-}Carlo}$ and $\mathsf{HK{\hbox{-}}Relax}$ in terms of clustering quality and efficiency (or memory overheads). Recall that $\mathsf{Monte\textrm{-}Carlo}$ is a random-walk-based approach as described in Section 3. $\mathsf{Monte\textrm{-}Carlo}$ accepts two thresholds $\epsilon_{r}$ and $\delta$, and a failure probability $p_{f}$ as inputs. It performs $\frac{2(1+\epsilon_{r}/3)\log{(n/p_{f})}}{\epsilon^{2}_{r}\cdot\delta}$ random walks from the seed node $s$ and returns a $(d,\epsilon_{r},\delta)$-approximate HKPR vector for $s$ with probability at least $1-p_{f}$. $\mathsf{HK{\hbox{-}}Relax}$ ensures $\epsilon_{a}$ absolute error in each $\frac{{\widehat{\boldsymbol{\rho}}_{s}}[v]}{d(v)}$, which is incomparable with the accuracy guarantees of other three methods, that is $(d,\epsilon_{r},\delta)$-approximation guarantees. Hence, we do not compare all algorithms under the same theoretical accuracy guarantees. Instead, we evaluate each method by its empirical clustering quality (i.e., conductance) and empirical running time (or memory overheads) with various parameter settings and find the method that achieves the best trade-off in terms of clustering quality and running time (or memory overheads). $\mathsf{HK{\hbox{-}}Relax}$ only has one internal parameter $\epsilon_{a}$, which is the absolute error threshold. We vary it in $\{10^{-8},10^{-7},10^{-6},10^{-5},10^{-4}\}$ in our experiments. Since $\mathsf{Monte\textrm{-}Carlo}$, $\mathsf{TEA}$, and $\mathsf{TEA+}$ have almost the same parameters, we set relative error threshold $\epsilon_{r}=0.5$, and $\delta$ is varied in $\{2\times 10^{-8},2\times 10^{-7},2\times 10^{-6},2\times 10^{-5},2\times 10^{-4}\}$ for all of them. $\mathsf{ClusterHKPR}$ has one internal parameter $\epsilon$, which is varied in $\{0.005,0.01,0.02,0.05,0.1\}$. $\mathsf{SimpleLocal}$ has one locality parameter $\delta$. We vary it in $\{0.005,0.01,0.02,0.05,0.1\}$. In addition, we vary the number of iterations of $\mathsf{CRD}$ in $\{7,10,15,20,30\}$ and keep other parameters default. For fair comparison, we let the $x$-axis be the average conductance of the output clusters and the $y$-axis be the average running time (or memory overheads), depicting the empirical clustering quality and empirical efficiency (or memory overheads), respectively.

Figure 4 shows the average conductance of the output clusters and the the average running time when varying the aforementioned parameters for each algorithm. We can make the following observations. First, for each algorithm, when the error thresholds (i.e., $\epsilon_{a},\epsilon$ and $\delta$) become smaller or the number of iterations increase, the conductance of the output clusters reduces (i.e., the quality of the output clusters improves), as well as the computational time goes up markedly, which accords with their theoretical time complexities. In particular, $\mathsf{SimpleLocal}$ incurs very high running time as well as poor cluster quality due to its high time complexity and the fact that it is mainly devised for recovering the cluster for a subset of nodes of the cluster rather than detecting a cluster for a single seed node. We also note that $\mathsf{CRD}$ shows a better performance than $\mathsf{SimpleLocal}$. However, it is still much slower than HKPR-based methods. Hence, we omit the results of $\mathsf{SimpleLocal}$ on the remaining datasets and that of $\mathsf{CRD}$ on Orkut, LiveJournal, 3D-grid, Twitter and Friendster.

Second, we can see that $\mathsf{Monte\textrm{-}Carlo}$ and $\mathsf{ClusterHKPR}$ take several minutes to finish one local clustering query on all graph datasets in order to find clusters with small conductance. Hence, it is outperformed by $\mathsf{HK{\hbox{-}}Relax}$, $\mathsf{TEA}$ and $\mathsf{TEA+}$ by $1$ to $3$ orders of magnitude when they output clusters with almost the same conductance. This is due to the fact that $\mathsf{Monte\textrm{-}Carlo}$ and $\mathsf{ClusterHKPR}$ require performing a large number of random walks. In fact, this is consistent with the experimental results in prior work (chung2015computing, ; shun2016parallel, ). Moreover, it can be observed that $\mathsf{HK{\hbox{-}}Relax}$ always runs faster than $\mathsf{Monte\textrm{-}Carlo}$ and $\mathsf{ClusterHKPR}$ and achieves more than $10\times$ speedup on all datasets except Orkut, Twitter, and Friendster. To explain this phenomenon, recall that $\mathsf{HK{\hbox{-}}Relax}$ requires iteratively pushing residuals to $1-K$-hop nodes from the seed node and $K$ is very large (see Section 6). However, on graphs with large average degrees (Orkut, Twitter, and Friendster) the number of push operations increases dramatically after several hops from the seed node.

Third, $\mathsf{TEA+}$ outperforms $\mathsf{HK{\hbox{-}}Relax}$ by more than $10\times$ speedup on PLC, Orkut, Twitter, and Friendster, and more than $4\times$ speedup on the rest of graphs. The speedup is achieved by new termination conditions for $\mathsf{HK{\hbox{-}}Push+}$ and the residue reduction method for reducing the number random walks. However, we note that the speedup on DBLP, Youtube, LiveJournal and 3D-grid is not as significant as that on PLC, Orkut, Twitter, and Friendster. The reason is that these graphs either have large clustering coefficients (yang2012def, ) or small average degrees. The first one can also be observed in our experiments. With the same parameters as inputs to all three algorithms, the conductance values of output clusters from PLC, Orkut, Twitter, and Friendster are clearly greater than those from the remaining four graphs (i.e., DBLP, Youtube, LiveJournal and 3D-grid). This implies that nodes in these four graphs are more likely to cluster together, and then residues on these graphs tend to be propagated within a small cluster of nodes when $\mathsf{HK{\hbox{-}}Push+}$ is performed. Now consider the second reason. Recall that $\mathsf{HK{\hbox{-}}Push+}$ iteratively distributes residues to neighbors evenly before any termination condition is satisfied. Since the average degrees are small, the residues that each node receives will be large. Additionally, it needs more iterations to distribute the residues to more nodes, which may not be done before termination. As a result of these two factors, a few nodes will hold large residues rather than many nodes holding small residues. Consequently, the residue reduction method in $\mathsf{TEA+}$ fails to significantly reduce the number of random walks. Note that $\mathsf{HK{\hbox{-}}Relax}$, $\mathsf{TEA}$ and $\mathsf{TEA+}$ all terminate very quickly on 3D-grid (less than 10 milliseconds), which is consistent with the observation in (shun2016parallel, ). This is due to the fact that each node in 3D-grid has six neighbors and the residues will drop below the residue threshold quickly after performing several rounds of push operations.

In addition, we also note that $\mathsf{TEA}$ fails to achieve considerable speedup compared with $\mathsf{HK{\hbox{-}}Relax}$. Especially on DBLP, Youtube, LiveJournal, and 3D-grid, $\mathsf{TEA}$’s performance degrades to the same level as that of $\mathsf{HK{\hbox{-}}Relax}$. This is also caused by aforementioned high clustering coefficients and small average degrees of these graphs. $\mathsf{TEA+}$ is around $4\times$ faster than $\mathsf{TEA}$ on Orkut, Twitter, and Friendster. In summary, our experiments demonstrate the power of new termination condition of $\mathsf{HK{\hbox{-}}Push+}$ and residue reduction method, which reduce many push operations and random walks without sacrificing the cluster qualities.

Figure 5 shows the memory overheads (including the space required to store the input graph) for all datasets by varying the error thresholds. First, we observe that the memory overheads increase with the reduction in error thresholds as more space is required to store residues and HKPR values. In particular, memory overheads on 3D-grid and Friendster remain almost stable for all algorithms. As each node in 3D-grid is connected to six neighbors, $\mathsf{ClusterHKPR}$, $\mathsf{Monte\textrm{-}Carlo}$, $\mathsf{HK{\hbox{-}}Relax}$, $\mathsf{TEA}$ and $\mathsf{TEA+}$ easily detect the large set of nodes around the seed nodes and the size of memory allocated for storing HKPR values for this set of nodes remains stable. In Friendster, the input graph size is too large and memory usage is not impacted remarkably. Second, we observe that the memory overheads of all algorithms are roughly comparable because space usage is dominated by the storage of the input graph, and hence, in terms of memory cost, there is no significant advantage to choose one algorithm over another.

7.5. Ranking Accuracy of Normalized HKPR

In this set of experiments, we evaluate the accuracy and efficiency of each method for computing normalized HKPR values (e.g., $\frac{{\boldsymbol{\rho}}_{s}[v]}{d(v)}$). First, we randomly select $50$ seed nodes and apply the power method (page1999pagerank, ) with $40$ iterations to compute the ground-truth normalized HKPR values (we omit the large datasets due to time and memory limitations). Following the experimental settings in Section 7.1, we run $\mathsf{HK{\hbox{-}}Relax}$, $\mathsf{ClusterHKPR}$, $\mathsf{Monte\textrm{-}Carlo}$, $\mathsf{TEA}$ and $\mathsf{TEA+}$ to generate normalized HKPR values for the selected seed nodes with varied error thresholds. Specifically, we vary $\epsilon$ in $\{10^{-8},10^{-7},10^{-6},10^{-5},10^{-4},10^{-3}\}$ for $\mathsf{HK{\hbox{-}}Relax}$, $\epsilon$ in $\{0.01,0.02,0.05,0.1,0.2,0.3\}$ for $\mathsf{ClusterHKPR}$, and we set $\epsilon_{r}=0.5$ and $\delta$ is varied in $\{2\times 10^{-8},2\times 10^{-7},2\times 10^{-6},2\times 10^{-5},2\times 10^{-4},2\times 10^{-3}\}$ for $\mathsf{Monte\textrm{-}Carlo}$, $\mathsf{TEA}$ and $\mathsf{TEA+}$, respectively. Then, we evaluate the accuracy of each method by using Normalized Discounted Cumulative Gain (NDCG) (jarvelin2000ir, ), which is a classic metric for evaluating ranking results.

Figure 6 reports the performance of each method on four datasets. We can make the following observations. First, as we reduce the error thresholds, both the running time and NDCG of each method increase markedly, which is consistent with their theoretical guarantees. Second, $\mathsf{TEA+}$ consistently incurs least running time while achieving the same NDCG compared to the competing methods. In addition, $\mathsf{TEA}$ is $2\times-8\times$ slower than $\mathsf{TEA+}$ while $\mathsf{HK{\hbox{-}}Relax}$ runs even more slowly. Especially on PLC and Orkut datasets, $\mathsf{HK{\hbox{-}}Relax}$’s performance degrades to the same level of $\mathsf{ClusterHKPR}$ and $\mathsf{Monte\textrm{-}Carlo}$. Although $\mathsf{ClusterHKPR}$ and $\mathsf{Monte\textrm{-}Carlo}$ also provide relative-error guarantees, they still incur the highest overheads because they require a large number of random walks. Third, we note that the efficiency and ranking accuracy results accord with the efficiency and clustering quality results reported in Section 7.4. This demonstrates the relationship between ranking accuracy of normalized HKPR and the quality of HKPR-based clustering algorithms, emphasizing why our methods produce clusters with smaller conductance than the competing ones.

7.6. Clusters Produced vs. Ground-truth

We collect the top 5,000 ground-truth communities in DBLP, Youtube, LiveJournal and Orkut datasets from (linksnap, ). We select 100 seed nodes from 100 known communities of size greater than 100 randomly as the query set. For all algorithms, we vary $t$ from $3$ to $10$ ($t>10$ would give us clusters with substantially lower quality) and their error thresholds respectively to produce clusters with highest average $F_{1}$-measure (i.e., harmonic mean of precision and recall). More specifically, we vary $\epsilon$ from $0.005$ to $0.35$ for $\mathsf{ClusterHKPR}$, $\epsilon$ from $10^{-8}$ to $10^{-1}$ for $\mathsf{HK{\hbox{-}}Relax}$. In addition, we fix $\epsilon_{r}=0.5$ and vary $\delta$ from $10^{-8}$ to $10^{-1}$ for $\mathsf{Monte\textrm{-}Carlo}$, $\mathsf{TEA}$ and $\mathsf{TEA+}$ respectively.

Table 8 reports the highest $F_{1}$ measure of each algorithm and their corresponding running times. $\mathsf{TEA+}$ consistently produces clusters with the best average $F_{1}$-measures and least running times for all datasets except DBLP. On DBLP, $\mathsf{TEA}$ produces clusters with the best average $F_{1}$-measure while $\mathsf{TEA+}$ produces clusters with slightly smaller $F_{1}$ measure but significantly faster. We also observe that $\mathsf{ClusterHKPR}$ and $\mathsf{Monte\textrm{-}Carlo}$ generate very similar results for all datasets. They run significantly slower than $\mathsf{TEA}$ and $\mathsf{TEA+}$ and also produce clusters with slightly smaller average $F_{1}$-measures than our methods. In addition, $\mathsf{HK{\hbox{-}}Relax}$ has the worst performance on most datasets. The only exception occurs on Orkut, where it produces clusters with the second best $F_{1}$-measure but $4\times$ slower than $\mathsf{TEA+}$.

7.7. Sensitivity Analysis to the Subgraph Characteristics

Next, we study the impact of query sets generated from subgraphs of different characteristics on clustering quality and efficiency. First, from each dataset of Youtube, PLC and Orkut, we select 250 subgraphs with different densities (lawler2001combinatorial, ) randomly. Then we sort the subgraphs by their densities in descending order (denoted as $\{SG_{1},SG_{2},\cdots,SG_{250}\}$). We pick 50 nodes from $SG_{1},\cdots,SG_{50}$ respectively to form a query set referred to as high-density seed nodes, 50 nodes from $SG_{100},\cdots,SG_{150}$ respectively as medium-density seed nodes, and 50 nodes from $SG_{200},\cdots,SG_{250}$ respectively as low-density seed nodes. We run $\mathsf{ClusterHKPR}$, $\mathsf{Monte\textrm{-}Carlo}$, $\mathsf{HK{\hbox{-}}Relax}$, $\mathsf{TEA}$ and $\mathsf{TEA+}$ with the same parameter settings as in Section 7.4 on these three query sets.

Figure 7 plots the average conductance of the output clusters and the average running times of all algorithms under different error thresholds for the three query sets. We report the results on DBLP, Youtube, PLC and Orkut here. We can make the following observations. First, $\mathsf{TEA}$ and $\mathsf{TEA+}$ are consistently faster than the existing approaches for all query sets. Second, the conductance values of each graph in Figures 7(i)-(l) are higher than the rest. This is because subgraphs with high densities have low conductance. Also, both $\mathsf{ClusterHKPR}$ and $\mathsf{Monte\textrm{-}Carlo}$ show similar results on all query sets for all datasets whereas $\mathsf{HK{\hbox{-}}Relax}$, $\mathsf{TEA}$ and $\mathsf{TEA+}$ are sensitive to the subgraph densities. Since seed nodes picked from subgraphs with high densities would have many neighbors, the residues in $\mathsf{HK{\hbox{-}}Relax}$, $\mathsf{TEA}$ and $\mathsf{TEA+}$ will drop quickly as push operations are performed, making them terminate quickly.

7.8. Effects of Heat Constant $t$

Lastly, we investigate the impact of the heat constant $t$. Using the same parameter settings and query set in Section 7.4, we run all algorithms on DBLP and PLC datasets by varying $t$ in $\{5,10,20,40\}$. Figures 8 and 9 plot the average running time and average conductance of the output clusters of each algorithm on DBLP and PLC respectively. Observe that the running time of each algorithm increases as we increase $t$, which is consistent with their time complexities. $\mathsf{ClusterHKPR}$ and $\mathsf{Monte\textrm{-}Carlo}$ are the slowest as $t$ changes. We further observe that the advantage of $\mathsf{TEA+}$ over competing methods is more prominent as $t$ becomes larger. More specifically, $\mathsf{TEA+}$ is around $4$ times faster than $\mathsf{HK{\hbox{-}}Relax}$ when $t=5$ and the speedup goes up to more than one order of magnitude when $t=40$ on DBLP. Similar observations can be made on PLC. In addition, we find that the conductance values of clusters produced by each algorithm with larger $t$ are smaller than those with smaller $t$. This shows that we can obtain clusters with small conductance by choose a large $t$. However, our “Clusters Produced vs. Ground Truth” experiment (see Section 7.6) reveals that clusters produced by all algorithms with a large $t$ are very different from the ground-truth. This is because algorithms with a large $t$ tend to give us a cluster of nodes that are far from the seed node. As a result, choosing a good $t$ is paramount for finding high quality clusters.

8. Conclusions

In this paper, we propose two novel heat-kernel-based local clustering algorithms, $\mathsf{TEA}$ and $\mathsf{TEA+}$, for computing approximate HKPR values and local graph clustering efficiently. Our algorithms bridge deterministic graph traversal with Monte-Carlo random walks in a non-trivial way, thereby overcoming the drawbacks of both and achieving significant gain in performance in comparison to the state-of-the-art local clustering techniques. Our experiments demonstrate that $\mathsf{TEA+}$ significantly outperforms the state-of-the-art heat-kernel-based algorithm by at least $4$ times on small graphs and up to one order of magnitude on large graphs in terms of computational time when producing clusters with the same qualities.

9. ACKNOWLEDGEMENTS

This research is supported in part by MOE, Singapore under grant MOE2015-T2-2-069 and grant RG128/18 (M4012086.020), by NUS, Singapore under an SUG, by NTU, Singapore under grant M4082311.020, and by National Natural Science Foundation of China (No. 61832017, No. 61772346 and U1809206).

Appendix A Choosing $K$ for TEA+

Recall that $K$ limits the push operations within $K$ hops from the seed node, which is an important parameter controlling the termination of $\mathsf{HK{\hbox{-}}Push+}$. If $K$ is too small, it will lead to a small overhead for $\mathsf{HK{\hbox{-}}Push+}$ but a significant cost for $k\textrm{-}\mathsf{RandomWalk}$, and vice versa. Consequently, it is crucial to find an appropriate value for $K$ so as to strike a good balance between $\mathsf{HK{\hbox{-}}Push+}$ and $\mathsf{Monte\textrm{-}Carlo}$. In particular, there are two factors we need to take into account when choosing $K$:

1. (1)

$K$ ought to be changed along with $\epsilon_{r}\cdot\delta$; 2. (2)

We need large $K$ for graphs with small degrees, and conversely, small $K$ for large-degree graphs.

This reason is as follows. If $K$ is a fixed constant and small, the number of push operations performed by $\mathsf{HK{\hbox{-}}Push+}$ will be bounded by a maximum number. Recall that in each iteration of $\mathsf{HK{\hbox{-}}Push+}$, we will perform $d(v)$ push operations, where $d(v)$ is the number of neighbors of current node $v$. This implies that the total number of push operations performed within $K$ hops from $s$ on the graph with small average degree will be rather limited compared with that of the graph with high average degree. As a result, if the number of push operations is not large enough, $\mathsf{TEA+}$ will degrade to $\mathsf{Monte\textrm{-}Carlo}$ and the overhead for random walks will be significant especially when the value of $\epsilon_{r}\cdot\delta$ is small. Considering these two factors, we propose to choose $K$ by the following equation:

[TABLE]

where $\bar{d}$ is the average degree of the input graph. The intuition behind Equation (20) is that, after $K$-hops of push operations, the residues tend to approach $\epsilon_{r}\cdot\delta$, and hence, the value for $\alpha$ would not be large.

Appendix B Proofs

B.1. Proof of Lemma 1

Proof.

This proof is based on induction. Given reserve vector ${\mathbf{q}_{s}}$ and residue vectors $\mathbf{r}^{(0)}_{s},\ldots\mathbf{r}^{(K)}_{s}$ constructed by the end of any interation in Algorithm 1, we define $\mathbf{f}_{s}[v]$ as

[TABLE]

First, let us consider the initial condition, in which all entries of $\mathbf{r}^{(0)}_{s}$ are zero, except $\mathbf{r}^{(0)}_{s}[s]=1$, and other vectors are zero vectors. This implies that

[TABLE]

Namely, Equation (4) holds in the initial case.

Furthermore, we assume that after $j$ iterations, reserve vector ${\mathbf{q}_{s}}$ and residue vectors $\mathbf{r}^{(0)}_{s},\ldots\mathbf{r}^{(K)}_{s}$ satisfy Equation (4), namely $\mathbf{f}_{s}[v]={\boldsymbol{\rho}}_{s}[v]$. In the beginning of the $(j+1)$-th iteration, we will conduct push operations with an entry $(w,i)$, where $w\in V$ and $i<K$. Then, we define the change of $\mathbf{f}_{s}[v]$ after this iteration as

[TABLE]

where $\mathbf{\widehat{q}}_{s}$ is the updated reserve vector, $\mathbf{\widehat{r}}^{(i)}_{s}$ and $\mathbf{\widehat{r}}^{(i+1)}_{s}$ are the updated residue vectors after performing push operations in this interation. Then, by Lines 4-7 of Algorithm 1, we have

[TABLE]

[TABLE]

and $\forall{u\in N(w)}$, $\textstyle\mathbf{\widehat{r}}^{(i+1)}_{s}[u]-\mathbf{r}^{(i+1)}_{s}[u]=\left(1-\frac{\eta(i)}{\psi(i)}\right)\cdot\mathbf{r}^{(i)}_{s}[w]\cdot\mathbf{P}[w,u].$ Hence, the change of $\mathbf{f}_{s}[v]$ after this iteration is

[TABLE]

Notice that

[TABLE]

Therefore, $\textstyle\Delta(w,i)=\textstyle\mathbf{\widehat{q}}_{s}[v]-{\mathbf{q}_{s}}[v]-\mathbf{r}^{(i)}_{s}[w]\cdot\mathbf{h}^{(i)}_{w}[v]+\mathbf{r}^{(i)}_{s}[w]\cdot\left(\mathbf{h}^{(i)}_{w}[v]-\frac{\eta(i)}{\psi(i)}\right)=0,$ which implies $\mathbf{f}_{s}[v]$ is invariant and still equals ${\boldsymbol{\rho}}_{s}[v]$ after the $(j+1)$-th iteration. This completes the proof of this lemma. ∎

B.2. Proof of Lemma 2

Proof.

Recall that Algorithm 2 consists of several iterations (see Lines 3-8). In the $\ell$-th ($\ell=0,1,\cdots$) iteration, the algorithm terminates with $\frac{\eta(k+\ell)}{\psi(k+\ell)}$ probability (Line 4); with the other $1-\frac{\eta(k+\ell)}{\psi(k+\ell)}$ probability, it samples a neighbor node of current node and set the sampled node as current node (Lines 6-8). Let $v_{i}$ be the node at the beginning of the $i$-th iteration, and ${p}(v,\ell)$ be the probability that the algorithm terminates at the $\ell$-th iteration with $v$ as the returned node. We will prove that, given $G,s=u,k$ as the input to Algorithm 2, we have

[TABLE]

Note that if Equation (21) holds, then the overall probability that $v$ is sampled as returned node is

[TABLE]

which establishes the lemma.

We prove by induction. For the base case where $\ell=0$, we have $v=s$, and the probability that the [math]-th iteration of Algorithm 2 terminates is $\frac{\eta(k)}{\psi(k)}$. In that case, $p(v,\ell)=\frac{\eta(k)}{\psi(k)}$ if $v=u$; otherwise $p(v,\ell)=0$. Meanwhile, it can be verified that when $\ell=0$, the r.h.s of Equation (21) equals $\frac{\eta(k)}{\psi(k)}$ if $v=u$, and [math] otherwise. Therefore, Equation (21) holds when $\ell=0$.

Assume that Equation (21) holds when $\ell=i$. Then, for any node $w\in V$, the probability that $w$ is at the beginning of $i$-th iteration is

[TABLE]

Now consider the $(i+1)$-th iteration. Since Algorithm 2 does not terminate at $i$-th iteration, for any node $v$, the probability that $v$ is picked as current node is

[TABLE]

In the $(i+1)$-th iteration, Algorithm 2 terminates with probability $\frac{\eta(k+i+1)}{\psi(k+i+1)}$, which implies that $v$ is returned with probability

[TABLE]

Therefore, the lemma is proved. ∎

B.3. Proof of Lemma 3

Proof.

Since the initialization of vectors in Algorithm 1 can be implemented in $O(1)$ cost by using hashmap, the total cost of Algorithm 1 is mainly determined by the number of push operations performed before termination. Recall that in each iteration of Algorithm 1 (Lines 4-7), if there exists an entry $(u,k)$ such that $u$’s current $k$-hop residue is greater than $r_{max}\cdot d(u)$, $\mathsf{HK{\hbox{-}}Push}$ will convert $\frac{\eta(k)}{\psi(k)}$ fraction of the current $k$-hop residue into its reserve. This implies that in each interation occurred on entry $(u,k)$, its reserve will be increased by at least $\frac{\eta(k)}{\psi(k)}\cdot r_{max}\cdot d(u)$. Let ${\mathbf{r}^{(k)}_{s}}[u]$ denote the total $k$-hop residue that $u$ receives and $x_{u,k}$ be the number of interations occurred on entry $(u,k)$. Then, we have

[TABLE]

and thus $d(u)\cdot x_{u,k}\leq\cdot\frac{{\mathbf{r}^{(k)}_{s}}[u]}{r_{max}}$. Since in each iteration occurred on entry $(u,k)$, $\mathsf{HK{\hbox{-}}Push}$ will perform $d(u)$ push operations, the total number of push operations occurred on all $(u,k)$ entries is then bounded by

[TABLE]

Since each push operation on entry $(u,k)$ will increase the value of $\mathbf{r}^{(k)}_{s}[u]$, the number of non-zero elements in residue vectors $\mathbf{r}^{(0)}_{s},\ldots,\mathbf{r}^{(K)}_{s}$ returned by Algorithm 1 is bounded by the total number of push operations, i.e., $O(\frac{1}{r_{max}})$.

Excluding the space required by the input graph, the space cost of Algorithm 1 is incurred by the reserve vector and the non-zero elements in residue vectors, which is $O(\frac{1}{r_{max}})$. Therefore, the lemma is proved. ∎

B.4. Proof of Lemma 4

Proof.

With regard to each invocation of $k\textrm{-}\mathsf{RandomWalk}$, the time complexity relies on the total number of nodes it visits in the random walk. According to Lemma 2, the probability that any node $v\in V$ is returned as the end node at $\ell$-th iteration is $p(v,\ell)$. If an invocation of $k\textrm{-}\mathsf{RandomWalk}$ terminates at $v$ in $\ell$-th iteration, then $\ell$ nodes are visited before it returns $v$. Considering all nodes in the graph, we have the expected time of each invocation of $k\textrm{-}\mathsf{RandomWalk}$ as follows

[TABLE]

Note that $\frac{t}{k}\cdot\eta(k-1)=\eta(k)\leq\sum_{\ell=k}^{\infty}{\eta(\ell)}=\psi(k)$. Then $\textstyle t\cdot\frac{\eta(k-1)}{\psi(k)}\leq k.$ This implies that Equation (22) is no greater than $t$, which finishes our proof.

∎

B.5. Proof of Theorem 6

Proof.

Let $Y_{i}$ be as defined in the context of Equation (9), and let $Y=\sum_{i=1}^{n_{r}}{Y_{i}}$. By Line 12 of $\mathsf{TEA}$, $\textstyle{\widehat{\boldsymbol{\rho}}_{s}}[v]={\mathbf{q}_{s}}[v]+\frac{Y\cdot\alpha}{n_{r}}.$ By Equation (10), the expectation of $Y$ is

[TABLE]

Let $\lambda=\frac{n_{r}\epsilon_{r}}{\alpha}\cdot{\boldsymbol{\rho}}_{s}[v]$ and $n_{r}$ be as defined in Line 8 of $\mathsf{TEA}$. By the Chernoff bound (see Lemma 5), for any node $v$ in $V$ with ${\boldsymbol{\rho}}_{s}[v]>d(v)\cdot\delta$, we have

[TABLE]

On the other hand, let $\lambda=\frac{n_{r}\epsilon_{r}\delta d(v)}{\alpha}$. Then, for any node $v$ in $V$ with ${\boldsymbol{\rho}}_{s}[v]\leq d(v)\cdot\delta$,

[TABLE]

By the union bound, for $V_{1}=\{v|v\in V\ s.t.\ {\boldsymbol{\rho}}_{s}[v]>d(v)\cdot\delta\}$, we have

[TABLE]

and for $V_{2}=\{v|v\in V\ s.t.\ {\boldsymbol{\rho}}_{s}[v]\leq d(v)\cdot\delta\}$, we have

[TABLE]

Therefore, we have the following results, respectively. With probability at least $1-\sum_{v\in V}{(p^{\prime}_{f})^{d(v)}}$, for every node $v$ in $V$ with $\frac{{\boldsymbol{\rho}}_{s}[v]}{d(v)}>\delta$, $\textstyle\left|\frac{{\widehat{\boldsymbol{\rho}}_{s}}[v]}{d(v)}-\frac{{\boldsymbol{\rho}}_{s}[v]}{d(v)}\right|\leq\epsilon_{r}\cdot\frac{{\boldsymbol{\rho}}_{s}[v]}{d(v)},$ and for every node $v$ in $V$ with $\frac{{\boldsymbol{\rho}}_{s}[v]}{d(v)}\leq\delta$, $\textstyle\left|\frac{{\widehat{\boldsymbol{\rho}}_{s}}[v]}{d(v)}-\frac{{\boldsymbol{\rho}}_{s}[v]}{d(v)}\right|\leq\epsilon_{r}\delta.$

By the definition of $p^{\prime}_{f}$ in Equation (6), if $\sum_{v\in V}{p_{f}^{d(v)-1}}\leq 1$, then $p^{\prime}_{f}=p_{f}$, which leads to $\textstyle\sum_{v\in V}{(p^{\prime}_{f})^{d(v)}}=\sum_{v\in V}{p_{f}^{d(v)}}\leq p_{f};$ otherwise, $p^{\prime}_{f}=\frac{p_{f}}{\sum_{v\in V}{p_{f}^{d(v)-1}}}$, hence $\textstyle\sum_{v\in V}{(p^{\prime}_{f})^{d(v)}}<\sum_{v\in V}{\frac{p_{f}^{d(v)}}{\sum_{v\in V}{p_{f}^{d(v)-1}}}}=p_{f}.$ Namely, ${\widehat{\boldsymbol{\rho}}_{s}}$ is a $(d,\epsilon_{r},\delta)$-approximate HKPR vector with probability at least $1-p_{f}$. ∎

B.6. Proof of Theorem 1

Proof.

By Lemma 3 and the definition of $\mathbf{h}^{(k)}_{u}[v]$ in Equation 5, we have $\textstyle\frac{\mathbf{h}^{(k)}_{u}[v]}{d(v)}=\frac{\mathbf{h}^{(k)}_{v}[u]}{d(u)}.$ Then we can rewrite Equation (4) as follows

[TABLE]

Now we prove that $\sum_{u\in V}{\mathbf{h}^{(k)}_{v}[u]}=1$ for any node $v\in V$ and $k\in[0,K]$.

[TABLE]

Hence, $\textstyle{\boldsymbol{\rho}}_{s}[v]-{\mathbf{q}_{s}}[v]\leq d(v)\cdot\sum_{k=0}^{K}{\max_{u\in V}{\left\{\frac{{\mathbf{r}^{(k)}_{s}}[u]}{d(u)}\right\}}}.$ Once Inequality (11) held, we have $\textstyle\frac{{\boldsymbol{\rho}}_{s}[v]}{d(v)}-\frac{{\mathbf{q}_{s}}[v]}{d(v)}\leq\epsilon_{a},$ which completes the proof. ∎

B.7. Proof of Theorem 4

Proof.

First, consider the case where $\mathsf{TEA+}$ terminates at Line 7, i.e., $\sum_{k=0}^{K}{\max_{u\in V}{\left\{\frac{{\mathbf{r}^{(k)}_{s}}[u]}{d(u)}\right\}}}\leq\epsilon_{a}$. In that case, by Theorem 1, for any node $v\in V$, $\textstyle\left|\frac{{\widehat{\boldsymbol{\rho}}_{s}}[v]}{d(v)}-\frac{{\boldsymbol{\rho}}_{s}[v]}{d(v)}\right|\leq\epsilon_{a}=\epsilon_{r}\cdot\delta.$ This indicates that, for any node $v$ with ${\boldsymbol{\rho}}_{s}[v]>d(v)\cdot\delta$, $\textstyle\left|\frac{{\widehat{\boldsymbol{\rho}}_{s}}[v]}{d(v)}-\frac{{\boldsymbol{\rho}}_{s}[v]}{d(v)}\right|\leq\epsilon_{r}\cdot\frac{{\boldsymbol{\rho}}_{s}[v]}{d(v)}.$ In addition, for any node $v$ with ${\boldsymbol{\rho}}_{s}[v]\leq d(v)\cdot\delta$, $\textstyle\left|\frac{{\widehat{\boldsymbol{\rho}}_{s}}[v]}{d(v)}-\frac{{\boldsymbol{\rho}}_{s}[v]}{d(v)}\right|\leq\epsilon_{r}\cdot\delta.$ Thus, ${\widehat{\boldsymbol{\rho}}_{s}}$ is a $(d,\epsilon_{r},\delta)$-approximate HKPR vector.

Now consider the case where $\mathsf{TEA+}$ does not terminate at Line 7. Let $Y_{i}$ be as defined in the context of Equation (9), and let $Y=\sum_{i=1}^{n_{r}}{Y_{i}}$. By Lines 17 and 19 of $\mathsf{TEA+}$,

[TABLE]

By Equation (10), the expectation of $Y$ is

[TABLE]

Let $\lambda=\frac{n_{r}\epsilon_{r}}{2\cdot\alpha}\cdot{\boldsymbol{\rho}}_{s}[v]$ and $n_{r}$ be as defined in Line 13 of $\mathsf{TEA+}$. By the Chernoff bound (see Lemma 5), for any node $v$ in $V$ with ${\boldsymbol{\rho}}_{s}[v]>d(v)\cdot\delta$, we have

[TABLE]

On the other hand, let $\lambda=\frac{n_{r}\epsilon_{r}\delta d(v)}{2\alpha}$. Then, for any node $v$ in $V$ with ${\boldsymbol{\rho}}_{s}[v]\leq d(v)\cdot\delta$, we have

[TABLE]

By union bound, for $V_{1}=\{v|v\in V\ s.t.\ {\boldsymbol{\rho}}_{s}[v]>d(v)\cdot\delta\}$, we have

[TABLE]

and for $V_{2}=\{v|v\in V\ s.t.\ {\boldsymbol{\rho}}_{s}[v]\leq d(v)\cdot\delta\}$, we have

[TABLE]

By Inequality (19), $\left|\frac{\mathbf{b}_{s}[v]}{d(v)}-\frac{\epsilon_{r}\delta}{2}\right|\leq\frac{\epsilon_{r}\delta}{2}$. Then, we have the following results, respectively. With probability at least $1-\sum_{v\in V}{(p^{\prime}_{f})^{d(v)}}$, for every node $v$ in $V$ with ${\boldsymbol{\rho}}_{s}[v]>d(v)\cdot\delta$,

[TABLE]

and for every node $v$ in $V$ with ${\boldsymbol{\rho}}_{s}[v]\leq d(v)\cdot\delta$,

[TABLE]

By the definition of $p^{\prime}_{f}$ in Equation (6), the total failure probability will be at most $\sum_{v\in V}{(p^{\prime}_{f})^{d(v)}}\leq p_{f}$. Therefore, ${\widehat{\boldsymbol{\rho}}_{s}}$ is a $(d,\epsilon_{r},\delta)$-approximate HKPR vector with probability at least $1-p_{f}$. ∎