Adaptive Shrinkage Estimation for Streaming Graphs

TL;DR

This paper introduces an adaptive, single-pass sampling framework with a James-Stein shrinkage estimator for efficient higher-order dependency estimation in massive, streaming networks, improving accuracy and computational efficiency.

Contribution

It proposes a novel adaptive sampling method and a shrinkage estimator tailored for streaming networks, enabling accurate higher-order analysis from partial data.

Findings

Outperforms baseline methods in large network experiments

Efficiently estimates higher-order dependencies with small samples

Incrementally updates in streaming environments

Abstract

Networks are a natural representation of complex systems across the sciences, and higher-order dependencies are central to the understanding and modeling of these systems. However, in many practical applications such as online social networks, networks are massive, dynamic, and naturally streaming, where pairwise interactions among vertices become available one at a time in some arbitrary order. The massive size and streaming nature of these networks allow only partial observation, since it is infeasible to analyze the entire network. Under such scenarios, it is challenging to study the higher-order structural and connectivity patterns of streaming networks. In this work, we consider the fundamental problem of estimating the higher-order dependencies using adaptive sampling. We propose a novel adaptive, single-pass sampling framework and unbiased estimators for higher-order network analysis of large streaming networks. Our algorithms exploit adaptive techniques to identify edges that are highly informative for efficiently estimating the higher-order structure of streaming networks from small sample data. We also introduce a novel James-Stein shrinkage estimator to reduce the estimation error. Our approach is fully analytic, computationally efficient, and can be incrementally updated in a streaming setting. Numerical experiments on large networks show that our approach is superior to baseline methods.

Tables

Table 1. Table 1: Summary of Graph Statistics

graph	$\left|V\right|$	$\left|K\right|$	$T$	$T_{\text{max}}$
soc-flickr	514K	3.2M	58.8M	2236
soc-livejournal	4.03M	27.9M	83.6M	586
soc-youtube	1.13M	2.98M	3.05M	4034
wiki-Talk	2.4M	4.7M	9.2M	1631
web-BerkStan-dir	685K	6.7M	64.7M	45057
cit-Patents	3.8M	16.5M	7.5M	591
soc-orkut-dir	3.07M	117.2M	627.6M	9145

Table 2. Table 2: MSE and Relative Spectral Norm at sampling fraction f = 0.2 𝑓 0.2 f=0.2 . APS : Adaptive Sampling, APS JS : APS with shrinkage Estimation, Unif : Uniform Sampling, Triest : Triest Sampling.

	Mean Squared Error (MSE)				${\parallel A-\hat{A}\parallel }_2/{\parallel A\parallel }_2$
graph	APS	APS JS	Unif	Triest	APS	APS JS	Unif	Triest
soc-flickr	22.30K	295.13	6.3K	7.46K	0.5793	0.0478	0.4321	0.5149
soc-livejournal	214.80	16.11	257.60	293.67	0.0269	0.0089	0.429	0.5092
soc-youtube-snap	11.35	6.68	119.79	145.87	0.0455	0.079	0.4159	0.4982
wiki-Talk	7.70	5.32	589.92	680.67	0.0105	0.0359	0.4315	0.5109
web-BerkStan-dir	7.32K	561.20	10.70K	14.03K	0.1169	0.0557	0.4381	0.6163
cit-Patents	6.02	3.03	10.59	10.91	0.0187	0.0428	0.4325	0.4914
soc-orkut-dir	2.08K	70.79	467.90	613.89	0.1086	0.0726	0.4385	0.4241

Table 3. Table 3: Summary of Notation

Notation	Description
$k_t$ (or just $t$)	Edge arriving at time $t$
${\hat{K}}_t$	Sample set after edge $t$ processed
${\hat{K}}_t^{\prime }$	Edges in reservoir prior to selection at time $t$
$J$	Generic edge subset
$J_t$	Edges from $J$ that have arrived by $t$
$S_{J,t}$	Indicator variable that indicates if all edges in $J$ have arrived by $t$
${\hat{S}}_{J,t}$ (${\hat{S}}_{J,t}^{\prime }$)	Inverse probability estimator of $S_{J,t}$  (estimator without last arriving edge)
$I_{J,t}$ ($I_{J,t}^{\prime }$)	Un-normalized estimator of $S_{J,t}$  (estimator not using last arriving edge)
$w_{i,t}$	Weight of edge $i$ at time $t\ge i$
$u_i$	IID uniform $(0,1]$ variable for edge $i$
$r_{i,t}$	Priority rank variable of edge $i$ at time $t\ge i$
$z_{J,t}$	Minimum priority rank of non-$J$ edges prior to $t$
$z_t$	$z_{\varnothing ,t}$ i.e., unrestricted minimum priority rank
$z_t^{\ast }$	Cumulative maximum of $z_{t^{\prime }}$ for $t^{\prime }\le t$
$H$ ($H_t$) ($H_{k,t}$)	Set of motifs (those with all their edges arrived by $t$)  (also containing edge $k$)
$n_{k,t}$	Total number of members of $H_t$ than contain $k$
${\hat{n}}_{k,t}$	Estimator of $n_{k,t}$
$\eta$	Generic James-Stein estimator for an edge count $n$
$\lambda$	Mixture parameter (i.e., shrinkage coefficient) in $\eta$
$p_{i,t}$	Probability of inclusion of edge $i\in {\hat{K}}_t$ at $t\ge i$
$P_{J,t}$	Probability of inclusion of edges from $J_t$ in ${\hat{K}}_t$ at time $t\ge i$
$t_J$	Minimum time over all edges in $J$, i.e. ${\min }_{i\in J}{t}_i$
${\tau }_J$	Time of the last arriving edge in $J$, ${\max }_{i\in J}{t}_i$

Table 4. Table 4: MSE for Non-Adaptive Sampling ( f = 0.2 𝑓 0.2 f=0.2 )

graph	Non-Adapt	Non-Adapt (JS)
soc-flickr	4907.21	2174.9
soc-livejournal	94.46	69.97
soc-youtube-snap	24.78	31.704
wiki-Talk	78.69	98.765
web-BerkStan-dir	1723.63	1236.3
cit-Patents	6.45	5.67
soc-orkut-dir	405.86	227.65

Equations

$$p_{i,t,J}=\{\begin{array}{ll}1,& \text{if }t<i\\ \min \{p_{i,t-1,J},w_{i,t}/z_{J,t}\},& \text{if }t\ge i\end{array}$$

$$\mathcal{L}(\lambda )=\mathrm{Var}(\hat{\eta })+(\mathbb{E}[\hat{\eta }]-n{)}^2={\lambda }^2\mathrm{Var}(\hat{n})+{\overline{\lambda }}^2\mathrm{Var}(w)+2\lambda \overline{\lambda }\mathrm{Cov}(\hat{n},w)+{\overline{\lambda }}^2\mathbb{E}[\hat{n}-w{]}^2$$

$${\mathcal{L}}^{\prime }(\lambda )/2=\lambda \mathrm{Var}(\hat{n})-\overline{\lambda }\mathrm{Var}(w)+(1-2\lambda )\mathrm{Cov}(\hat{n},w)-\overline{\lambda }(\mathbb{E}[\hat{n}-w]{)}^2$$

$$\lambda =1-\frac{\mathrm{Cov}(\hat{n}-w,\hat{n})}{\mathbb{E}[(\hat{n}-w{)}^2]}=1-\frac{\mathrm{Var}(\hat{n})-\mathrm{Cov}(\hat{n},w)}{\mathbb{E}[(\hat{n}-w{)}^2]}$$

$$\mathrm{Var}({\hat{n}}_{j,t})=\sum _{J,L\in H_{j,t}}\mathrm{Cov}({\hat{S}}_{J,t}^{\prime },{\hat{S}}_{L,t}^{\prime })$$

$$\mathop{\mathrm{Var}}(\widehat{n}_{j,t})=\mathop{\mathrm{Var}}(\widehat{n}_{j,t-1})+\sum_{J\in\Delta_{j,t}}\big{[}\mathop{\mathrm{Var}}(\widehat{S}^{\prime}_{J,t})+2\mathop{\mathrm{Cov}}(\widehat{n}_{j,t-1},\widehat{S}^{\prime}_{J,t})+\sum_{\begin{subarray}{c}L\in\Delta_{j,t}\\ L\neq J\end{subarray}}\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{J,t},\widehat{S}^{\prime}_{L,t})\big{]}$$

$$C_{J,t_1;L,t_2}={\hat{S}}_{J,t_1}^{\prime }{\hat{S}}_{L,t_2}^{\prime }-{\hat{S}}_{J\setminus L,t_1}^{\prime }{\hat{S}}_{L\setminus J,t_2}^{\prime }{\hat{S}}_{J\cap L,t_1\vee t_2}^{\prime }$$

$$\mathbb{E}[I(J_t\subset {\hat{K}}_t)|{\mathcal{Z}}_{J,t},J_{t-1}\subset {\hat{K}}_{t-1}]$$

$$\mathbb{E}[{\tilde{S}}_{J,t}|{\mathcal{Z}}_{J,t},J_{t-1}\subset {\hat{K}}_{t-1}]$$

$$\mathbb{E}[{\tilde{S}}_{J,t}|{\mathcal{Z}}_{J,t-1},J_{t-1}\subset {\hat{K}}_{t-1}]={\tilde{S}}_{J,t-1}$$

$$z_t=\frac{w_{d_t,t}}{u_{d_t}}\ge \frac{w_{d_t,s}}{u_{d_t}}>z_s$$

$$\frac{w_{i,t+1}}{z_{t+1}}\ge \frac{w_{i,t+1}}{z_{t+1}^{\ast }}\ge \frac{w_{i,t}}{z_{t+1}^{\ast }}=\frac{w_{i,t}}{z_t^{\ast }}$$

$$\mathrm{Cov}({\hat{S}}_{J,t_1}^{\prime },{\hat{S}}_{L,t_2}^{\prime })$$

$$=\mathbb{E}[{\hat{S}}_{J,t_1}^{\prime }{\hat{S}}_{L,t_2}^{\prime }]-1$$

$$\mathbb{E}[{\hat{S}}_{J\setminus L,t_1}^{\prime }{\hat{S}}_{L\setminus J,t_2}^{\prime }{\hat{S}}_{J\cap L,t_1\vee t_2}^{\prime }]=1$$

$$\mathbb{E}[{\hat{S}}_{J_1,t_1}^{\prime }{I}_{J_2,t_2}^{\prime }|{\mathcal{Z}}_{J,t_2},J_{t_2-1}\subset {\hat{K}}_{t_2-1}]$$

$$\mathbb{E}[{\hat{S}}_{J_1,t_1}^{\prime }{P}_{J_1\cap J_2,t_2}{I}_{J_2\setminus J_1,t_2}^{\prime }|{\mathcal{Z}}_{J,t_2},J_{t_2-1}\subset {\hat{K}}_{t_2-1}]$$

$${\hat{n}}_i$$

$${\hat{n}}_j$$

$${\hat{n}}_k$$

$$\mathrm{Var}({\hat{n}}_i)$$

$$\mathrm{Var}({\hat{n}}_j)$$

$$\mathrm{Var}({\hat{n}}_k)$$

$$\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{J,t},\widehat{S}^{\prime}_{L,s})=(p_{i,t}p_{j,t})^{-1}\big{(}p_{i,s}^{-1}-1\big{)}p_{j^{\prime},s}^{-1}$$

$$\sum _{s<t}\mathrm{Cov}({\hat{S}}_{J,t}^{\prime },{\hat{S}}_{L,s}^{\prime })$$

$$=(p_{i,t}{p}_{j,t}{)}^{-1}\ast U_{i,t}$$

$$\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{J,t},\widehat{S}^{\prime}_{L,s})=(p_{i,t}p_{j,t})^{-1}\big{(}p_{j,s}^{-1}-1\big{)}p_{i^{\prime},s}^{-1}$$

$$\sum _{s<t}\mathrm{Cov}({\hat{S}}_{J,t}^{\prime },{\hat{S}}_{L,s}^{\prime })$$

$$=(p_{i,t}{p}_{j,t}{)}^{-1}\ast U_{j,t}$$

$$\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{J,t},\widehat{S}^{\prime}_{L,s})=(p_{i,t}p_{j,t})^{-1}\big{(}(p_{i^{\prime},s}p_{j^{\prime},s})^{-1}-(p_{i^{\prime},s}p_{j^{\prime},s})^{-1}\big{)}=0$$

$$\mathrm{Var}({\hat{n}}_i)$$

$$\mathrm{Var}({\hat{n}}_j)$$

$$U_{i,t}$$

$$U_{j,t}$$

$$\mathrm{Cov}({\hat{S}}_{J,t}^{\prime },{\hat{I}}_{L,s}^{\prime })=(p_{i,t}{p}_{j,t}{)}^{-1}(1-p_{i,s})$$

Full text

1 Introduction

Network analysis has been central to the understanding and modeling of large complex systems in various domains, e.g., social, biological, neural, and technological systems [7, 37]. These complex systems are usually represented as a network (graph) where vertices represent the components of the system, and edges represent their direct (observed) interactions over time. The success of network analysis throughout the sciences rests on the ability to describe the complex structure and dynamics of arbitrary systems using only observed pairwise interaction data among the components of the system. Many networked systems exhibit rich structural and connectivity patterns that can be captured at the level of pairwise links (edges) or individual vertices. However, higher-order dependencies that capture complex forms of interactions remain largely unknown, since they are beyond the reach of methods that focus primarily on pairwise links. Recently, there has been a surge of studies on higher-order network analysis [4, 9, 52, 43, 20]. These methods focus on generalizing the analysis and modeling of network data from pairwise relationships (e.g., edges) to more complex forms of relationships such as multi-node (many-body) relationships (e.g., motif patterns, hypergraphs) and higher-order network paths that depend on more history [46]. Higher-order connectivity patterns were shown to change node rankings [46, 57], reshape the community structure [52, 9, 56], reveal the hub structure [4], learn more accurate embeddings [42, 41], and generative network models [16].

Many networks are massive, dynamic, and naturally streaming over time [33, 44, 3], with pairwise interactions (i.e., edges that represent communication in the form of user-to-user, user-to-product interactions) are becoming available one at a time in some arbitrary order (e.g., online social networks, Emails, Twitter data, recommendation engines). The massive size and streaming nature of these networks allow only partial observation, since it is infeasible to analyze the entire network. Under such scenarios, the question of how to study and reveal the higher-order connectivity structure and patterns of streaming networks has remained a challenge. This work is motivated by large-scale streaming network data that are generated by measurement processes (i.e., from online social media, sensors, and communication devices), and we study how to estimate the higher-order connectivity structure of streaming networks under the constraints of partial observation and limited memory. We particularly focus on the estimation of higher-order network patterns captured by small subgraphs, also called network motifs (e.g., triangles or small cliques) [34, 6].

Randomization and sampling techniques are fundamental in the context of graph and matrix approximations in both static and streaming settings; see [33, 29, 26, 5]. The general problem is setup as follows: given a graph $G=(V,K)$ and a budget $m$, find a sampled graph $\widehat{G}$ such that the (expected) number of edges (non-zero entries) is at most $m$ and $\widehat{G}$ is a good proxy for $G$. In the data streaming model, the input graph $G$ is a stream of edges $K=\{k_{1}=(u,v),k_{2}=(v,w)\dots\}$ and is partially observed as the edges stream and become available to the algorithm one at a time in some arbitrary order. The streaming model is fundamental to applications of online social networks, social media, and recommendation systems where network data become available one at a time (e.g., friendship links, emails, Twitter feeds, user-item preferences, purchase transactions, etc). Moreover, the streaming model is also crucial where network data is streaming from disk storage and random accesses of edges are too expensive. However, the theory and algorithms of current graph sampling techniques are mostly well developed for sampling individual edges to estimate global network properties (e.g., total number of edges in a graph) [25, 50]. Here, we consider instead sampling techniques that can capture how edges connect locally to form small network substructures (i.e., network motifs). Designing new sampling algorithms to estimate the local higher-order connectivity patterns of streaming networks has the potential to improve accuracy and efficiency of sampling and knowledge discovery in streaming networks.

Contributions. We propose a novel topologically adaptive, single-pass priority sampling framework for unbiased estimation of higher-order network connectivity structure of large streaming networks, where edges become available one at a time in some arbitrary order. Specifically, we propose unbiased estimators for local counts of subgraphs or motifs containing each edge (Theorem 1) and show how to compute them efficiently for streaming networks (Theorem 2). These estimators are embodied in our proposed adaptive sampling framework (see Algorithm 1).

Our proposed adaptive sampling preferentially selects edges to include in the sample based on their importance weight relative to the variable of interest (i.e., higher-order graph properties), then adapts their weights to allow edges to gain importance during stream processing leading to reduction in estimation variance as compared with static and/or uniform weights.

We also propose a novel shrinkage estimator which we formulate as a convex combination estimator to reduce the mean squared error (MSE) (as shown in Figure 1), and we discuss its computation during stream processing (Section 3). Our approach is fully analytic, computationally efficient, and can be incrementally updated as the edges become available one at a time during stream processing. The proposed methods are also generally applicable to a wide variety of networks, including directed, undirected, weighted, and heterogeneous networks.

2 Adaptive Sampling Framework

2.1 Notation and Problem Definition

Consider an arriving stream $K$ of unique graph edges labelled by the edge identifiers $k\in[|K|]$. Let $G=(V,K)$ denote the undirected graph formed by the edges, where $V$ is the vertex set and $K$ is the edge set. Assume $M$ is a motif (subgraph) pattern of interest, let $H$ denote the class of subgraphs in $G$ that are isomorphic to $M$ (e.g., all triangles or cliques of a given size that appear in $G$). We define the $H$-weighted graph of $G$ as the weighted graph $G_{H}=(V,K,N)$ with edge weights $N=\{n_{k}:k\in K\}$, such that for each edge $k\in K$, $n_{k}$ is the number of subgraphs in $H$ that are isomorphic to motif $M$ and incident to $k$, i.e., $n_{k}=|\{h\in H:h\owns k,h\cong M\}|$. We refer to this graph as the motif-weighted graph, and we denote $\mathrm{A}$ as its motif adjacency matrix [9]. For brevity we will identify a subgraph $h\in H$ with its edge set. Table 3 in the supplementary materials provides a summary of notation. Suppose the edges of $G$ are labelled in some arbitrary order based on their arrival in the stream. Let $G_{t}=(V_{t},K_{t})$ denote the subgraph of $G$ formed by the first $t$ edges in this order, $H_{t}=\{h\in H:h\subset K_{t}\}$ be the set of subgraphs in $H$ all of whose edges have arrived by $t$, and $(V_{t},K_{t},N_{t})$ be the corresponding $H$-weighted graph of $G_{t}$ (with weights $N_{t}=\{n_{k,t}:k\in K_{t}\}$). This paper studies two questions: (1) how to maintain a reservoir sample $\widehat{K}$ of $m$ edges from the unweighted edge stream $K$, and (2) how to obtain an unbiased estimate of the $H$-weighted graph $G_{H}=(V_{t},K_{t},N_{t})$ at any time $t\in[|K|]$. We propose a variable-weight adaptive sampling framework for streaming network/graph data, called adaptive priority sampling. Our proposed framework preferentially selects edges to include in the sample based on their importance weight, where the weights are relative to the role of these edges in the formation of motifs and general subgraphs of interest (e.g., triangles or small cliques) and can adapt to the changing topology during streaming. Next, we describe the proposed framework (Alg. 1), and discuss its theoretical foundation.

2.2 Algorithm Description and Key Intuition

We consider a generic reservoir sample $\widehat{K}$ selected progressively from the edge stream labelled $K=[|K|]=\{1,2,\ldots,|K|\}$. We assume edges are unique, and therefore they can be identified by their arrival positions (i.e., edge ids); nevertheless we will sometimes emphasize their graph or time aspects, denoting by $k_{t}$ the edge arriving at time slot $t$, and by $t_{k}$ the arrival time slot of edge $k$. In Alg. 1, the first $m$ edges are admitted to the sample: $\widehat{K}_{t}=[t]$ for $t\leq m$. Then, each subsequent edge $t$ is provisionally included in the current sample to form $\widehat{K}^{\prime}_{t}=\widehat{K}_{t-1}\cup\{t\}$ (see line 8), from which an edge is discarded to produce the sample $\widehat{K}_{t}$, and maintain the sample size $m=|\widehat{K}_{t}|$ at any time $t$.

In Algorithm 1, each edge $i\in\widehat{K}^{\prime}_{t}$ is assigned a priority rank variable defined as $r_{i,t}=w_{i,t}/u_{i}$, where $w_{i,t}$ is the edge weight at time $t$, and $u_{i}$ is a uniformly distributed random variable on $(0,1]$ assigned to the edge on its first arrival. Then, the edge with minimum rank $z_{t}=\min_{j\in\widehat{K}^{\prime}_{t}}r_{j,t}$ is discarded from $\widehat{K}^{\prime}_{t}$ to obtain the sample $\widehat{K}_{t}$ (see lines 19–22). For each edge $i\in\widehat{K}^{\prime}_{t}$, we compute the weight $w_{i,t}>0$ as a function of its previous weight $w_{i,t-1}$ and the sample set $\widehat{K}^{\prime}_{t}$.

Upon its arrival, a new edge $k$ is assigned an IID edge random variable $u_{k}$ uniformly distributed on $(0,1]$, and an initial (constant) weight $\phi$ (lines 5–7), plus the number of target subgraphs/motifs in $\widehat{K}^{\prime}_{t}$ that contains $k$ (see lines 11–17). An edge $i\in\widehat{K}^{\prime}_{t}$ survives the sampling at time $t$, if and only if there is another edge in $\widehat{K}^{\prime}_{t}$ that has the minimum rank, i.e., $r_{i,t}>z_{t}$.

Conditional on $z_{t}$, the effective sampling probability of an edge $i\in\widehat{K}_{t}$ is: $\mathbb{P}\{r_{i,t}>z_{t}\}=\mathbb{P}\{u_{i}<w_{i,t}/z_{t}\}=\min\{1,w_{i,t}/z_{t}\}$. We note that in the experiments of Section 4, we choose the initial edge weight $\phi=1$ to be comparable with the edge weight increments due to subgraphs incident to each edge (see line 6). This procedure allows edges to have a chance to be included in the sample with a non-zero probability, regardless of the number of subgraphs incident to them, but not so large as to damp out their topological weight. Next, we discuss how the approach in Algorithm 1 leads to unbiased estimators of general subgraphs/motifs.

2.3 Unbiased Estimators of General Subgraphs

Let $S_{i,t}$ denote the arrival of an edge $i$, i.e., $S_{i,t}=I(i\leq t)$. For any subgraph $J\subset K$, where $J$ is a subset of edges (or edge ids), let $S_{J,t}=\prod_{i\in J}S_{i,t}$ indicates whether all edges $i\in J$ have arrived by time $t$, i.e., $S_{J,t}=1$ if $J\subset K_{t}$ and [math] otherwise. We observe the local edge count $n_{i,t}=\sum_{J\in H_{i,t}}S_{J,t}$, and $H_{i,t}=\{h\in H_{t}:h\owns i\}$ is the set of subgraphs (motifs) incident to edge $i$ whose edges have arrived by time $t$.

Theorem 1 establishes unbiased inverse probability estimators [23] for $S_{J,t}$ in the form $\widehat{S}_{J,t}=I(J\subset\widehat{K}_{t})/P_{J,t}$ when $t\geq\tau_{J}:=\max_{i\in J}t_{i}$ (i.e., all edges in $J$ have arrived by time $t$), and $P_{J,t}$ is the sampling probability for the subgraph $J$. For any subgraph $J\subset K$ with $|J|\leq m\leq t$, let $J_{t}=J\cap[t]$, and define the conditional minimum edge rank over the sample $\widehat{K}^{\prime}_{t}$ as $z_{J,t}=\min_{j\in\widehat{K}^{\prime}_{t}\setminus J_{t}}r_{j,t}$. Hence, $z_{t}=z_{\emptyset,t}$ is the unrestricted minimum rank over $\widehat{K}^{\prime}_{t}$. For $i\in J$, we define the edge probabilities $p_{i,t,J}$ to be $1$ when $t<i$ and $\min\{1,\min_{i\leq s\leq t}{w_{i,s}}/{z_{J,t}}\}$ otherwise. This can be expressed in an iterative form as follows,

[TABLE]

We distinguish between $\widetilde{P}_{J,t}$ and $P_{J,t}$. We use $\widetilde{P}_{J,t}=\prod_{i\in J_{t}}p_{i,t,J}$ to denote the sampling probability of subgraph $J$ at time $t$, conditional on the ranks of edges not in $J$ (i.e., using the conditional min rank $z_{J,t}$). We also use $P_{J,t}=\prod_{i\in J_{t}}p_{i,t}$, where $p_{i,t}:=p_{i,t,\emptyset}$, to denote the sampling probability of subgraph $J$ that employs the threshold $z_{t}=z_{\emptyset,t}$, i.e., $z_{t}$ is the unrestricted minimum rank over $\widehat{K}^{\prime}_{t}$.

Set $t_{J}=\min_{i\in J}t_{i}$, then define $\widetilde{S}_{J_{t}}=I(J_{t}\in\widehat{K}_{t})/\widetilde{P}_{J,t}$ and the set of variables ${\mathcal{Z}}_{J,t}=\{z_{J,s}:\ t_{J}\leq s\leq t\}$. In Theorem 1, we establish first that $\widetilde{S}_{J,t}$ is an unbiased estimator of $S_{J,t}$, but that estimates can be computed using $\widehat{S}_{J,t}$. This is preferable since $P_{J,t}$ is computed using the unrestricted threshold $z_{t}$, independent of the subgraph $J$ to be estimated.

Theorem 1** (Unbiased Subgraph Estimation111Proofs of all the theorems are discussed in the supplementary materials.).**

1. (I)

The distributions of the edge random variables $\{u_{i}:i\in J\}$, conditional on $J_{t}\subset\widehat{K}_{t}$ and ${\mathcal{Z}}_{J,t}$, are independent, with each $u_{i}$ being uniformly distributed on $(0,p_{i,J,t}]$. 2. (II)

${\mathbb{E}}[I(J_{t}\subset K_{t})|{\mathcal{Z}}_{J,t},J_{t-1}\subset\widehat{K}_{t-1}]=\widetilde{P}_{J,t}/\widetilde{P}_{J,t-1}$** 3. (III)

${\mathbb{E}}[\widetilde{S}_{J,t}|{\mathcal{Z}}_{J,t-1},J_{t-1}\subset\widehat{K}_{t-1}]=\widetilde{S}_{J,t-1}$, and hence ${\mathbb{E}}[\widetilde{S}_{J,t}]=1$, for $t>t_{J}$. 4. (IV)

$\widetilde{P}_{J,t}=P_{J,t}$* when $J_{t}\in\widehat{K}_{t}$ and hence ${\mathbb{E}}[\widehat{S}_{J,t}]=S_{J,t}$, for all $t$.*

Using Theorem 1, it is straightforward to show that for any edge $i\in\widehat{K}_{t}$, $\widehat{n}_{i,t}=\sum_{J\in H_{i,t}}\widehat{S}_{J,t}$ is an unbiased estimator of $n_{i,t}$, i.e. $\mathbb{E}[\widehat{n}_{i,t}]=n_{i,t}$.

Unbiased Estimation from the Last Arriving Edge.

Recall that $\tau_{J}=\max_{i\in J}t_{i}$ denotes the time of the last arriving edge $k_{\tau_{J}}$ of the subgraph $J\subset K$. Set $J^{(0)}=J\setminus\{k_{\tau_{J}}\},$ and define $\widehat{S}^{\prime}_{J,t}=\widehat{S}_{J^{(0)},\tau_{J}-1}$, where $S^{\prime}_{J,t}$ indicates subgraph $J$ right before the arrival of the last edge $k_{\tau_{J}}$.

In Alg. 1, when a new edge arrives at time $t=\tau_{J}$, Algorithm 1 finds all subgraphs $\Delta\subset H_{t}$ that are completed by the arriving edge and whose edges are in the sample $\widehat{K}^{\prime}_{t}$ (see line 10). For each subgraph $J\in\Delta$ and each edge $i\in J$, we increment the estimate $\widehat{n}_{i,t}$ by the inverse probability $1/P_{J^{(0)},t-1}$, where $P_{J^{(0)},t-1}=\prod_{i\in J^{(0)}}p_{i,t-1}$ is the sampling probability for $S^{\prime}_{J,t}$ (lines 15–16).

Corollary 1 results from Theorem 1 and establishes that ${\mathbb{E}}[\widehat{S}^{\prime}_{J,t}]=1$, hence, $\widehat{n}_{i,t}=\sum_{J\in H_{i,t}}\widehat{S}^{\prime}_{J,t}$ is an unbiased estimator for $n_{i,t}$, for all $i\in K_{t}$. This allows us to update the estimates without risking loss of some edge in $J$ during subsequent sampling (i.e., when the edge with minimum rank is discarded from the sample).

Corollary 1**.**

${\mathbb{E}}[\widehat{S}^{\prime}_{J,t}]=1$* and hence $\widehat{n}_{i,t}=\sum_{J\in H_{i,t}}\widehat{S}^{\prime}_{J,t}$ is an unbiased estimator of the local subgraph count $n_{i,t}$ for all $i\in K_{t}$.*

2.4 Special Case of Non-decreasing Sampling Weights

Computing the probabilities $p_{i,t}$ according to Equation 1 requires an update for each each edge $i\in\widehat{K}_{t}$ at each time step $t$, i.e., $O(m)$ for each arriving edge. We now show that this computational cost can be reduced when $w_{i,t}$ is non-decreasing in $t$. Let $d_{t}\leq t$ denote the edge discarded at time $t>m$, i.e., $\{d_{t}\}=\widehat{K}^{\prime}_{t}\setminus\widehat{K}_{t}$ (line 22 in Alg. 1). We define the sample threshold $z^{*}_{t}$ iteratively by $z^{*}_{m}=0$ and $z^{*}_{t}=\max\{z^{*}_{t-1},z_{t}\}$ , for $t>m$ (see line 21 in Algorithm 1). Define $p^{*}_{i,i}=\min\{1,w_{i,i}/z^{*}_{i}\}$ and $p^{*}_{i,t+1}=\min\{p^{*}_{i,t},w_{i,t+1}/z^{*}_{t+1}\}$, for $t\geq i$, i.e., similar to Equation 1 but with $z_{t}$ replaced by $z^{*}_{t}$ ( as shown in line 13 in Alg. 1).

Theorem 2**.**

When $w_{i,t}$ is non-decreasing in $t$ then (I) $d_{t}\neq t$ implies $z^{*}_{t}=z_{t}$; and (II) $p^{*}_{i,t}=p_{i,t}$ for all $t\geq i$.

We take advantage of Theorem 2 to reduce the number of updates to the probability $p^{*}_{i,t}$, Since $w_{i,t}$ is non-decreasing and $z^{*}_{t}$ is also non-decreasing, $w_{i,t}/z^{*}_{t}$ can only increase when $w_{i,t}$ increases.

During the intervals of constant $w_{i,t}$, $w_{i,t}/z^{*}_{t}$ is non-increasing. Therefore, provided that we update $p^{*}_{i,t}$ at times when $w_{i,t}$ increases, all other updates of $p^{*}_{i,t}$ can be deferred until needed for estimation (see line 13 of Alg. 1).

Complexity Analysis. In Algorithm 1, the sampling reservoir is implemented as a min-heap. Any insertion, deletion, update operation has $O(\log m)$ complexity in the worst case. Retrieving the edge with minimum rank is done in constant time $O(1)$. The complexity of the weight update depends on the target subgraph class, being proportional to the number of edges in new subgraphs created by the arriving edge. In the experiments reported in this paper, the target subgraphs are triangles. For an arriving edge $k=(v_{1},v_{2})$, the third vertex of any new triangle incident to $k$ lies in the set intersection of the sampled neighbors of $v_{1}$ and $v_{2}$ which can be computed in $O(\min\{\textrm{deg}(v_{1}),\textrm{deg}(v_{2})\})$, where $\textrm{deg}(v_{1})$ and $\textrm{deg}(v_{2})$ are the sampled vertex degrees of $v_{1}$ and $v_{2}$ respectively. This complexity can be achieved if a hash table (or Bloom filter) is used for storing and looping over the sampled neighborhood of the vertex with minimum degree and querying the hash table of the other vertex.

3 James-Stein Shrinkage Estimator

It is common in graph sampling to seek unbiased estimators with minimum variance that perform well, e.g., the estimator in Section 2. In this section, we also investigate another desirable estimator, called shrinkage estimator [24, 21], that directly reduces the mean squared error (MSE), which is a direct measure of estimation error. In Figure 1, we demonstrate the bias-variance trade-off in graph sampling, which leads to both biased and unbiased estimators. Unbiased estimators of local subgraph counts are subject to high relative variance when the motif counts are small, because in this case the individual count estimates, scaled by the inverse probabilities, are smoothed less by aggregation.

More generally, James and Stein originated the observation that unbiased estimators do not necessarily minimize the mean squared error [24]. In their study, unbiased estimates of high dimensional Gaussian random variables are adjusted through scaling-based regularization and linear combination with dimensional averages. Shrinkage estimation has been used in other settings such as covariance or affinity matrix estimation [45, 55, 11, 28]. Here, we examine shrinkage for the estimated count $\widehat{n}_{k}$ by convex combination with the observed and un-normalized count provided by the edge sampling weight $w_{k}$. By introducing bias through $w_{k}$, we can obtain further reductions in mean squared error (MSE), additional to the adaptive sampling technique discussed in Section 2.

3.1 Optimizing Shrinkage Coefficients

We define a family of shrinkage estimators $\eta=\lambda\widehat{n}+{\overline{\lambda}}w$, where the shrinkage coefficient $\lambda\in[0,1]$ specifies $\eta$ as a convex combination of the unbiased estimator $\widehat{n}=\widehat{n}_{k}$ and the un-normalized edge weight $w=w_{k}$, for any edge $k$. Let ${\overline{\lambda}}$ denote $1-\lambda$. The loss $\mathcal{L}(\lambda)$ associated with the shrinkage coefficient $\lambda$ is the mean squared error:

[TABLE]

since ${\mathbb{E}}[\widehat{\eta}-n]={\mathbb{E}}[\widehat{\eta}-\widehat{n}]={\mathbb{E}}[\lambda\widehat{n}+{\overline{\lambda}}w-\widehat{n}]={\overline{\lambda}}{\mathbb{E}}[w-\widehat{n}]$.

$\mathcal{L}$ is convex with derivative $\mathcal{L}^{\prime}$ specified by,

[TABLE]

We seek the minimum of $\mathcal{L}$ when $\mathcal{L}^{\prime}(\lambda)=0$, i.e., when

[TABLE]

We truncate ${\overline{\lambda}}$ at $1$ so that the constraint ${\overline{\lambda}}\leq 1$ always holds. Since the optimal $\lambda$ is a function of the unknown true covariances, we follow the practice of [12] by employing a plug-in estimator $\widehat{\lambda}$ for $\lambda$ by substituting $(\widehat{n}-w)^{2}$ in the denominator, and an unbiased estimate for $\mathop{\mathrm{Cov}}(\widehat{n}-w,\widehat{n})={\mathop{\mathrm{Var}}(\widehat{n}_{k})}-{\mathop{\mathrm{Cov}}(\widehat{n}_{k},w_{k})}$, whose computation we describe next.

3.2 Unbiased Estimation of the Variance ${\mathop{\mathrm{Var}}(\widehat{n})}$

Let $\Delta_{j,t}=H_{j,t}\setminus H_{j,t-1}$ denote the set of subgraphs in $K_{t}$ that contain an edge $j$ and are completed by the new edge arrival at time $t$. Similarly, let $\widehat{\Delta}_{j,t}$ denote the (possibly empty) set of subgraphs in $\widehat{K}^{\prime}_{t}$ that contain an edge $j\in K^{\prime}_{t}$ and are completed by the new edge arrival at time $t$. Thus, the estimated count $\widehat{n}_{j,t}$ can be decomposed as: $\widehat{n}_{j,t}=\widehat{n}_{j,t-1}+\sum_{J\in\Delta_{j,t}}\widehat{S}^{\prime}_{J,t}$.

For any pair of subgraphs $J,L\in H_{j,t}$, the variance of $\widehat{n}_{j,t}$ is specified by:

[TABLE]

where $\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{J,t},\widehat{S}^{\prime}_{L,t})$ is the covariance between two subgraph estimators. Furthermore, the variance $\mathop{\mathrm{Var}}(\widehat{n}_{j,t})$ can also be computed incrementally at each time $t$ as follows,

[TABLE]

where the term $\mathop{\mathrm{Cov}}(\widehat{n}_{j,t-1},\widehat{S}^{\prime}_{J,t})=\sum_{s<t}\sum_{L\in\Delta_{j,s}}\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{J,t},\widehat{S}^{\prime}_{L,s})$, for $s<t$.

Theorem 3 is used to establish an unbiased estimator for $\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{J,t},\widehat{S}^{\prime}_{L,s})$ in the form,

[TABLE]

where $t_{1}\geq t_{2}$, and $t_{1}\vee t_{2}=\max\{t_{1},t_{2}\}$.

Theorem 3**.**

$C_{J,t_{1};L,t_{2}}$* is an unbiased estimator of $\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{J,t_{1}},\widehat{S}^{\prime}_{L,t_{2}})$, for some time $t_{1}\geq t_{2}$.*

A special case of Theorem 3 happens when $J=L$ and $t_{1}=t_{2}=t$, which leads to $V(\widehat{S}^{\prime}_{J,t})=\widehat{S}^{\prime}_{J,t}(\widehat{S}^{\prime}_{J,t}-1)$, where $V(\widehat{S}^{\prime}_{J,t})$ is an unbiased estimator of $\mathop{\mathrm{Var}}(\widehat{S}^{\prime}_{J,t})$.

3.3 Unbiased Estimation of the Covariance ${\mathop{\mathrm{Cov}}(\widehat{n},w)}$

Following the notation in Section 3.2, for each edge $j$, the weight $w_{j,t}$ is a random quantity incremented by $1$ for each subgraph $J\in\Delta_{j,t}$ completed by the new edge arrival at time $t$. Thus, $w_{j,t}$ can be written as a sum of random counts, i.e., un-normalized indicator functions analogous to how $\widehat{n}_{j,t}$ is written as a sum of inverse probability estimators. Let $I_{J,t}=I(J\subset\widehat{K}_{t})$ be the indicator of subgraph $J$, and recall that $J^{(0)}$ is the subgraph $J$ without the last arriving edge $k_{\tau_{J}}$. Define $I^{\prime}_{J,t}=I_{J_{0},\tau_{J}-1}$, i.e., the indicator that all edges but the final edge are present in the sample $\widehat{K}_{t-1}$ immediately before the arrival of the final edge ($k_{\tau_{J}}$ of $J$). When the new edge $k_{\tau_{J}}$ arrives at time $t=\tau_{J}$, each edge in $J^{(0)}$ has its weight incremented; see line 14 of Algorithm 1. Thus, we can write $w_{j,t}=\sum_{J\in H_{j,t}}I^{\prime}_{J,t}$, analogous to Corollary 1, and decompose $w_{j,t}=w_{j,t-1}+\sum_{J\in\Delta_{j,t}}I^{\prime}_{J,t}$.

Computing the optimal skrinkage $\lambda$ estimator in Equation 4 requires estimates of the covariance ${\mathop{\mathrm{Cov}}(\widehat{n}_{j,t},w_{j,t})}$ for each edge $j\in\widehat{K}_{t}$, which is estimated in turn and follow by linearity from the estimates of the covariance ${\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{J,t},I^{\prime}_{J,t})}$. Theorem 4 establishes an unbiased estimator for the general case of $\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{J_{1},t_{1}},I^{\prime}_{J_{2},t_{2}})$, when $t_{1}\geq t_{2}$. Lemma 1 is central to both the proof of Theorem 4 and the computation of covariance estimates222The computational details and proofs for shrinkage estimation are discussed with examples in the supplementary materials.

Lemma 1**.**

*For $J_{1}\cap J_{2}=\emptyset$ and $t_{1}\geq t_{2}$, then ${\mathbb{E}}[\widehat{S}^{\prime}_{J_{1},t_{1}}I^{\prime}_{J_{2},t_{2}}]={\mathbb{E}}[I^{\prime}_{J_{2},t_{2}}]$ and hence $\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{J_{1},t_{1}},I^{\prime}_{J_{2},t_{2}})=0$. *

Theorem 4** (Unbiased Subgraph Covariance Estimation).**

1. (I)

When $t_{1}\geq t_{2}$, $\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{J_{1},t_{1}},I^{\prime}_{J_{2},t_{2}})$ has unbiased estimator $D_{J_{1},t_{1};J_{2},t_{2}}=\widehat{S}^{\prime}_{J_{1},t_{1}}I^{\prime}_{J_{2},t_{2}}-\widehat{S}^{\prime}_{J_{1}\setminus J_{2},t_{1}}\widehat{S}^{\prime}_{J_{1}\cap J_{2},t_{1}\vee t_{2}}P_{J_{1}\cap J_{2},t_{2}}I^{\prime}_{J_{2}\setminus J_{1},t_{2}}$. 2. (II)

$D_{J_{1},t_{1};J_{2},t_{2}}>0$* iff $\widehat{S}^{\prime}_{J_{1},t_{1}}>0$ and $I^{\prime}_{J_{2},t_{2}}>0$. Hence $D_{J_{1},t_{1};J_{2},t_{2}}$ can be computed from samples that have been taken.* 3. (III)

For the special case $J_{1}=J_{2}=J$ and $t_{1}=t_{2}=t$ then $D_{J,t;J,t}=\widehat{S}^{\prime}_{J,t}{\overline{P}}_{J,t}=I^{\prime}_{J,t}(P_{J,t}^{-1}-1)$.

4 Experiments & Discussion

Experimental Setup.

We test on graphs from different domains and with different characteristics; see [40] for data downloads. Table 1 provides a summary of dataset characteristics, where $|V|$ is the number of vertcies, $|K|$ is the number of edges, $T$ is the number of triangles, and $T_{\textrm{max}}$ is the maximum triangle count per edge. For all graph datasets, we consider an undirected, unweighted, simplified graph without self loops.

Edge streams are obtained by randomly permuting the edges in each graph, and the same edge order is used for all the methods. We repeat the experiment ten different times with sample fractions $f=\{0.10,0.20,0.40,0.50\}$. All experiments were performed using a server with two Intel Xeon E5-2687W 3.1GHz CPUs, 256GB of memory. The experiments are executed independently for each sample fraction. Additional results and ablation studies are discussed in the supplementary materials. Our experimental setup is summarized as follows:

• •

For each sample fraction, we use Algorithm 1 to collect a sample $\widehat{K}$, from edge stream $K$.

• •

The experiments in this section use triangles as an example of the motif pattern $M$. However, the approach itself is general and applicable to any motif patterns.

• •

During stream processing, we compute unbiased estimators and James-Stein shrinkage estimators of the local triangle counts for the sampled edges, as discussed in Sections 2–3.

• •

Given a sample $\widehat{K}\subset K$, we compute the mean squared error (MSE), and the relative spectral norm [1], $\lVert\mathrm{A}-\widehat{\mathrm{A}}\rVert_{2}/\lVert\mathrm{A}\rVert_{2}$, where $\mathrm{A}$ is the exact triangle-weighted adjacency matrix of the input graph, $\widehat{\mathrm{A}}$ is the average estimated triangle-weighted adjacency matrix of the sampled graph, and $\lVert\mathrm{A}\rVert_{2}$ is the spectral norm of $\mathrm{A}$.

• •

We compare the results of Algorithm 1 with uniform sampling (i.e., reservoir sampling [53]) using the Horvitz-Thompson estimator, and we also compare with Triest sampling [48]. All baseline methods use the same experimental setup as the proposed method.

4.1 Comparison to Baseline Methods

We collect a sample of edges $\widehat{K}\subset K$ from the edge stream $K$ in a single pass, which we use to construct the motif-weighted graph, where $M$ is the triangle motif and $\mathrm{A}$ is adjacency matrix of the triangle-weighted graph. We use $\widehat{\mathrm{A}}$ to denote the estimator of $\mathrm{A}$ obtained by sampling. We compute the shrinkage estimator as discussed in Section 3. And, we report the MSE at sample fraction $f=0.20$ in Table 2, which demonstrates the following insight: the shrinkage estimator applied to adaptive priority (APS) sampling significantly improves the performance of the vanilla APS which uses Horvitz-Thompson estimator for all graphs. This is particularly clear for soc-flickr and soc-orkut for which the APS shrinkage (APS JS) significantly outperforms all the other methods.

We also consider the spectral norm as another measure of approximation quality in addition to MSE. The spectral norm $\lVert\mathrm{A}-\widehat{\mathrm{A}}\rVert_{2}$ was previously used for matrix approximation [1]. $\lVert\mathrm{A}-\widehat{\mathrm{A}}\rVert_{2}$ measures the strongest linear trend of $A$ that is not captured by the estimator $\widehat{\mathrm{A}}$. This is different from the mean squared error which focused on the magnitude of the estimates.

We report the relative spectral norm (i.e., $\lVert\mathrm{A}-\widehat{\mathrm{A}}\rVert_{2}/\lVert\mathrm{A}\rVert_{2}$) at sample fraction $f=0.20$ for various graphs in Table 2. The experiments demonstrate that for all of the example graphs, both APS and APS with shrinkage significantly outperform uniform reservoir sampling and Triest sampling. One observed exception is the soc-flickr graph, where the estimates using APS is significantly high due to the high variance of Horvitz-Thompson estimation for edges with small counts. Under such scenarios, the APS with shrinkage significantly helps and improves the original APS estimates. We also notice the difference between how the MSE ranks the best methods versus the relative spectral norm. A good example of this is the soc-orkut graph, for which APS performs worse than the baselines. However, APS is superior to uniform sampling and Triest sampling for the relative spectral norm. Thus, despite of the large mean squared error, APS (even without shrinkage) captures the linear trend and structure of the data better than uniform reservoir sampling and Triest sampling. Finally, Figure 2 shows the convergence performance of relative spectral norm as a function of the sampling fraction. Notably, APS and APS with shrinkage converge faster than uniform and Triest sampling, and we observe that shrinkage estimation significantly improves the vanilla APS.

4.2 Analysis of the Estimated Distribution

We take the top-k non-zero edge weights of the exact triangle-weighted adjacency matrix $\mathrm{A}$, and we compare them against their corresponding estimates obtained by sampling. Figures 3 shows the top-1M weights for APS with shrinkage estimation. Similar figures for uniform sampling and Triest sampling are reported in Section D of the supplementary materials (Fig 8 and Fig 9 respectively). The results demonstrate the more accurate performance of APS with shrinkage estimation compared to the baseline methods; more specifically, APS with shrinkage estimation preserves the distribution and ranks of the top-k edge weights compared to uniform and Triest sampling. We report the analysis for two sampling fractions $f=\{0.20,0.40\}$.

In Figure 4, we compare APS against APS with shrinkage estimation for the soc-livejournal graph. The results show how the shrinkage estimator reduces the variance of APS, in particular for small local counts with high variance (i.e., as observed in the tail of the edge weight distribution). In Section C in the supplementary materials, we discuss an ablation study of Algorithm 1.

5 Related Work

Here, we categorize the related work in three research areas: (1) Higher-order Network Analysis, (2) Graph Approximation, and (3) IID Stream Sampling.

Higher-order Network Analysis. There has been an increasing interest in higher-order network analysis and modeling in particular to generalize pairwise links to many-body relationships with arbitrary node sets and motifs; see [34, 9, 52, 56, 4, 54, 43, 20, 41, 16]. The majority of these methods focus on small static networks that fit in memory.

Graph Approximation. Randomization in the context of graph approximation is a well-studied topic; see [13, 22, 29, 49] and [33, 3] for a survey. Much work was devoted for triangle count approximation and other motifs for static graphs (see [10, 51, 47, 50, 17, 38]) and for streaming graphs (see [8, 48, 5, 25, 32, 2]). In the streaming setting, most work focused on estimating point statistics using fixed probabilities, e.g., the global triangle or motif count using reservoir based sampling approaches; see [53]. In this paper, we focus instead on estimating the motif-weighted graph from a stream of unweighted edges, and propose a general novel methodology for adaptive priority sampling with shrinkage estimation. We compare against the state-of-the-art approach, Triest sampling [48] and we obtain significant improvement over their method. Triest sampling maintains a sample of edges from the stream using reservoir sampling [53] and random pairing [18] to exploit the available memory as much as possible. However, our approach provides a sampling framework in which edges are included in the reservoir sample based on their importance and topological relevance in the formation of local motifs and subgraphs of interest, and edge weights are allowed to adapt to the changing topology of the reservoir sample.

IID Stream Sampling. Prior work focused on IID streams (e.g., IP networks, DB transactions, etc), e.g., single-pass reservoir sampling ([27, 36, 53]), order and threshold sampling ([14, 39, 15]), and probability proportional to size sampling (IPPS). These methods were designed for sampling IID data streams (e.g., IP networks, DB transactions, etc). Here, we focus instead on streaming graphs (non-iid data). Thus, the prior work on IID streams cannot be directly applied in this setting where the focus is on higher-order subgraphs, and extending these methods to non-IID streams is subject to further research.

Broader Impact

There is a burgeoning recent literature of statistical estimation and adaptive data analysis of the higher-order structural properties of graphs in both the streaming and non streaming context that reflect the importance and interest of this topic for the graph algorithms and relational learning research community. On the other hand, shrinkage estimators are an established technique from more general statistics. This paper is the first to apply shrinkage based methods in the context of graph approximation. The expected broader impact is as a proof of concept that shows the way for other researchers in this area to improve estimation quality. Moreover, this work fits under statistical inference for temporal relational/network data, which would enable statistical analysis and learning for network data that appear in streaming settings, in particular when exact solutions are not feasible (similar to the important literature on randomization algorithms for data matrices [1]).

Furthermore, there are many applications where the data has a pronounced temporal, relational, and spatial structure (e.g., relational data). Examples of Non-IID streams include (i) non-independence due to temporal clustering in communication graphs on internet, online social networks, physical contact networks, and social media such as flash crowds and coordinated botnet activity; (ii) non-identical distributions in activity on these networks due to diurnal and other seasonal variations, synchronization of user network activity e.g., searches stimulated by hourly news reports. The proposed framework is suitable for these applications, because it makes no statistical assumptions concerning the arrival stream and the order of the arriving edges.

Appendix A Theorem Proofs

Proof of Theorem 1.

Proof.

Any subgraph $J$ can be defined as a subset of edges from the set of all edges $K$. Suppose $J_{t}\subset\widehat{K}^{\prime}_{t}$, then $J_{t}$ survives the sampling at time $t$ (i.e., $J_{t}\subset\widehat{K}_{t}$), if and only if another edge $j\in\widehat{K}^{\prime}_{t}\setminus J_{t}$ has minimum rank $z_{J,t}=\min_{j\in\widehat{K}^{\prime}_{t}\setminus J_{t}}r_{j,t}$, i.e., if $r_{i,t}>z_{J,t}$, or equivalently, $u_{i}<w_{i,t}/z_{J,t}$ for all $i\in J_{t}$. Denote $A_{i,J,s}=\{u_{i}<w_{i,s}/z_{J,s}\}$ as the event when $i\in J_{s}\cap\widehat{K}_{s}$. Then for $t_{J}\leq\tau_{J}\leq t$, the event $\{J\subset\widehat{K}_{t}\}$ decomposes as $\bigcap_{t_{J}\leq s\leq t}B_{J,s}$ where $B_{J,s}=\bigcap_{i\in J_{s}}A_{i,J,s}$.

1. (I)

The proof is by induction on $t$. For $t<t_{J}$ the conditioning is trivial and $u_{i}$ are IID on $(0,1]=(0,p_{i,J,t}]$. The same property holds at general $t$ for all $i\in J$ which have not yet arrived, i.e., for $i\in J\setminus J_{t}$. Consider now $t\geq t_{J}$ and assume that the result holds for $t-1$. The weights $w_{i,t}$ for $i\in J_{t}\cap\widehat{K}^{\prime}_{t}$ are fixed by the conditioning on the event $\{J_{t-1}\subset\widehat{K}_{t-1}\}$. Further conditioning on $z_{J,t}$ and $J_{t}\subset\widehat{K}_{t}$ requires $u_{i}<w_{i,t}/z_{J,t}$ for all $i\in J_{t}\subset\widehat{K}_{t}$. Imposing this condition on the assumed independent uniform distributions of $u_{i}$ on $(0,p_{i,J,t-1}]$ results in independent uniform distributions of $u_{i}$ on $(0,\min\{p_{i,J,t-1},w_{i,t}/z_{J,t}\}]=(0,p_{i,J,t}]$. 2. (II)

The conditional expectation of the indicator $I(J_{t}\subset\widehat{K}_{t})$ is,

[TABLE]

where in the last step we have used the statement of part (I) for the distribution of $u_{i}$ conditioning on ${\mathcal{Z}}_{J,t}$ and $\{J_{t-1}\subset\widehat{K}_{t-1}\}$, since $w_{i,t}$ is assumed determined given $\widehat{K}_{t-1}$. 3. (III)

By using (II), we find that the conditional expectation of $\widetilde{S}_{J,t}$ is:

[TABLE]

which is independent of the conditioning on $z_{J,t}$ and hence,

[TABLE]

The initial value (for the first edge arrival at time $t_{J}$) is $\widetilde{S}_{J,t_{J}}=I(t_{J}\in\widehat{K}_{t_{J}})/p_{t_{J},J,t_{J}}=I(u_{t_{J}}<w_{t_{J},t_{J}}/z_{J,t_{J}})/p_{t_{J},J,t_{J}}$. Clearly ${\mathbb{E}}[\widetilde{S}_{J,t_{J}}|z_{J,t_{J}}]=1$ and hence ${\mathbb{E}}[\widetilde{S}_{J,t_{J}}]=1$. Finally ${\mathbb{E}}[\widetilde{S}_{J,t}]=1$ for all $t\geq t_{J}$ by chaining the conditional expectations. 4. (IV)

Trivially $\widehat{S}_{J,t}=S_{J,t}=0$ for $t<\tau_{J}$. Since $z_{J,t}=z_{t}$ when $J\subset\widehat{K}_{t}$, $P_{J,t}=\widetilde{P}_{J,t}$ and hence $\widehat{S}_{J,t}=\widetilde{S}_{J,t}$ for $t\geq\tau_{J}$ and ${\mathbb{E}}[\widehat{S}_{J,t}]=1$ by (III).

∎

Proof of Theorem 2.

Proof.

1. (I)

If $d_{t}\neq t$, $t$ is admitted to the sample and hence

[TABLE]

for all $s\in[d_{t},t]$. Since edge $d_{t}$ is discarded at time $t$, and $d_{t}\neq t$, then the minimum rank $z_{t}=r_{d_{t},t}=w_{d_{t},t}/u_{d_{t}}$.

The first inequality follows from the non-decreasing property of $w_{d_{t},t}$. The second inequality follows since edge $d_{t}$ survives the sampling from time $d_{t}$ until $t$ and hence its rank cannot be lower than the threshold $z_{s}$ for any $s$ in that interval. But since the edge $d_{t}$ was admitted to the sample at time, we have $d_{d_{t}}\neq d_{t}$, where $d_{d_{t}}$ is the discarded edge at time $d_{t}$. Hence, we apply the argument back recursively to the first sampling time. Hence, $z^{*}_{t}=\max\{z^{*}_{t-1},z_{t}\}=z_{t}$. 2. (II)

By assumption if an edge $i$ is admitted to $\widehat{K}_{i}$, then $i\neq d_{i}$ and so by (I) and Equation 1, $p_{i,i}=\min\{1,w_{i,i}/z_{i}\}=\min\{1,w_{i,i}/z^{*}_{i}\}=p^{*}_{i,i}$. The general case is by induction. Assume $p_{i,s}=p^{*}_{i,s}$ for all times $s\in[i,t]$, and $z_{t+1}>z^{*}_{t}$, then $z^{*}_{t+1}=z_{t+1}$ hence $p^{*}_{i,t+1}=p_{i,t+1}$.

If $z_{t+1}\leq z^{*}_{t}$, then $z^{*}_{t+1}=z^{*}_{t}$ and hence

[TABLE]

Thus we can replace $z_{t+1}$ by $z^{*}_{t+1}$ in (1) but use of either leaves the iterated value unchanged, since by the induction hypothesis, both are greater than $p_{i,t}\leq w_{i,t}/z^{*}_{t}$

∎

Proof of Theorem 3.

Proof.

[TABLE]

From Theorem 1, and since $J\setminus L$, $L\setminus J$, and $J\cap L$ are disjoint subsets, we have,

[TABLE]

Thus, $\mathbb{E}[C_{J,t_{1};L,t_{2}}]=\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{J,t_{1}},\widehat{S}^{\prime}_{L,t_{2}})=\mathbb{E}[\widehat{S}^{\prime}_{J,t_{1}}\widehat{S}^{\prime}_{L,t_{2}}]-1$.

A special case of Theorem 3 happens when $J=L$ and $t_{1}=t_{2}=t$, which leads to $V(\widehat{S}^{\prime}_{J,t})=\widehat{S}^{\prime}_{J,t}(\widehat{S}^{\prime}_{J,t}-1)$, where $V(\widehat{S}^{\prime}_{J,t})$ is an unbiased estimator of $\mathop{\mathrm{Var}}(\widehat{S}^{\prime}_{J,t})$. ∎

Proof of Lemma 1.

Proof.

Let $J=J_{1}\cup J_{2}$. Chaining conditional expectations from Theorem 1(III)

[TABLE]

using Theorem 1(I). Hence ${\mathbb{E}}[\widehat{S}^{\prime}_{J_{1},t_{1}}I^{\prime}_{J_{2},t_{2}}]={\mathbb{E}}[I^{\prime}_{J_{2},t_{2}}]$ and since ${\mathbb{E}}[\widehat{S}^{\prime}_{J_{1},t_{1}}]=1$, then the $\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{J_{1},t_{1}},I^{\prime}_{J_{2},t_{2}})={\mathbb{E}}[\widehat{S}^{\prime}_{J_{1},t_{1}}I^{\prime}_{J_{2},t_{2}}]-{\mathbb{E}}[\widehat{S}^{\prime}_{J_{1},t_{1}}]{\mathbb{E}}[I^{\prime}_{J_{2},t_{2}}]=0$. ∎

Proof of Theorem 4.

Proof.

1. (I)

Since ${\mathbb{E}}[\widehat{S}^{\prime}_{J_{1},t_{1}}]=1$ it suffices to show that (the negative of) the second term in the definition of $D_{J_{1},t_{1};J_{2},t_{2}}$ in Theorem 4(i) has expectation ${\mathbb{E}}[I^{\prime}_{J_{2},t_{2}}]$. When $t_{1}\geq t_{2}$ then repeating the conditioning argument of Lemma 1, this term has conditional expectation

[TABLE]

and hence the stated property holds. 2. (II)

Holds since $\widehat{S}^{\prime}_{J,t}>0$ implies $I^{\prime}_{J,t^{\prime}}>0$ for $t\geq t^{\prime}\geq t_{J}$ 3. (III)

is a special case of (I).

∎

Appendix B Example: Estimators for Local Triangle Counts

Assume the motif $M$ is a triangle in the form $J=(i,j,k)$, where the edges in the triangle are ordered by their arrival times, i.e., $i<j<k$. Let $k$ denote the new edge arriving at time $t$, and $\widehat{\Delta}_{t}=\{J=(i,j,k)\subset\widehat{K}^{\prime}_{t}\}$ be the set of new triangles completed by $k$ at time $t$. We now show how the estimators can be incremented for each triangle. Note that edges $i,j\in\widehat{K}^{\prime}_{t}$ can participate in only one triangle at time $t$.

Unbiased estimator for ${\widehat{n}}$.

By applying Theorem 1, each triangle $J=(i,j,k)\in\widehat{\Delta}_{t}$ results in an increment of $\widehat{S}^{\prime}_{J,t}=1/({p_{i,t}p_{j,t}})$ in the count estimator for each edge in the triangle as follows:

[TABLE]

Unbiased estimator for ${\mathop{\mathrm{Var}}(\widehat{n})}$.

By applying Theorem 3, each triangle $J=(i,j,k)\in\widehat{\Delta}_{t}$ results in an increment of \mathop{\mathrm{Var}}(S^{\prime}_{J,t})=\big{(}1/(p_{i,t}p_{j,t})-1\big{)}/(p_{i,t}p_{j,t}) in the variance estimator of the count for each edge in the triangle as follows:

[TABLE]

${\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{J,t},\widehat{S}^{\prime}_{L,s})}$.

By applying Theorem 3, we detail all the possible cases for the computation of the covariance ${\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{J,t},\widehat{S}^{\prime}_{L,s})}$, where $L=(i^{\prime},j^{\prime},k^{\prime})$ is another triangle, and $L\neq J$:

$J\cap L=\emptyset$: if the two triangles are disjoint, then $\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{J,t},\widehat{S}^{\prime}_{L,s})=0$, see Figure 5 for an example. 2. 2.

$s=t$: assume $L=(i^{\prime},j^{\prime},k)\in\widehat{\Delta}_{t}$ is another triangle completed by $k$, and $L\neq J$. This means that $J\cap L=\{k\}$, (see Figure 5), and $\widehat{S}^{\prime}_{J\cap L,t\vee S}=1$. Then, the estimator of the covariance $\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{J,t},\widehat{S}^{\prime}_{L,s})=0$. 3. 3.

$s<t$: assume $L=(i^{\prime},j^{\prime},k_{s})\in\widehat{\Delta}_{s}$ is another triangle completed by edge $k_{s}$ at time $s$, for any $s<t$.

1. (a)

If $i=i^{\prime}$ and $L=(i,j^{\prime},k_{s})$, then the two triangles overlap in the edge $i$, and $\widehat{S}^{\prime}_{J\cap L,t\vee S}=1/p_{i,t}$. Thus, the estimator of the covariance is,

[TABLE]

Thus, for all triangles $L=(i,j^{\prime},k_{s})$, for $s<t$

[TABLE]

where U_{i,t}=\sum_{s<t}\big{(}p_{i,s}^{-1}-1\big{)}p_{j^{\prime},s}^{-1} 2. (b)

If $j=j^{\prime}$ and $L=(i^{\prime},j,k_{s})$, then similar to the previous case, then the estimator of the covariance is,

[TABLE]

Thus, for all triangles $L=(i^{\prime},j,k_{s})$, for any $s<t$

[TABLE]

where U_{j,t}=\sum_{s<t}\big{(}p_{j,s}^{-1}-1\big{)}p_{i^{\prime},s}^{-1}. 3. (c)

if $k_{s}=i$ or $k_{s}=j$, then the estimator of the covariance is zero,

[TABLE]

To facilitate incremental covariance computations for streaming data, we define $U_{i,t}$ and $U_{j,t}$ as the cumulative sum variables for edges $i$ and $j$ respectively, to keep track of previously sampled triangle estimators that contain $i$ and $j$ respectively, at any time $s<t$. Note that for the new arriving edge $k$, we have $U_{k,t}=0$. Now, we add the covariance increments to each edge as follows,

[TABLE]

Then, to update the cumulative variables for edges $i,j\in J=(i,j,k)$.

[TABLE]

Unbiased Estimator for ${\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{J,t},\widehat{I}^{\prime}_{L,s})}$.

By applying Theorem 4, we detail the computation of the covariance ${\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{J,t},\widehat{I}^{\prime}_{L,s})}$:

If $J\cap L=\emptyset$, then from Lemma 1, the ${\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{J,t},\widehat{I}^{\prime}_{L,s})}=0$. 2. 2.

If $s=t$ and $J=L$, then the ${\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{J,t},\widehat{I}^{\prime}_{J,t})}=(p_{i,t}p_{j,t})^{-1}-1$. 3. 3.

If $s=t$ and $J\neq L$, then $J\cap L=\{k\}$. And from Theorem 4 (I), the ${\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{J,t},\widehat{I}^{\prime}_{L,t})}=0$ 4. 4.

If $s<t$, and $L=(i^{\prime},j^{\prime},k_{s})$ is a triangle completed by edge $k_{s}$ at time $s$ then,

1. (a)

If $i=i^{\prime}$ and $L=(i,j^{\prime},k_{s})$, then $J\cap L=\{i\}$, and the covariance estimator is,

[TABLE]

And, for all triangles $L=(i,j^{\prime},k_{s})$, for any $s<t$,

[TABLE]

where $D_{i,t}=\sum_{s<t}(1-p_{i,s})$. 2. (b)

if $j=j^{\prime}$ and $L=(i^{\prime},j,k_{s})$, then $J\cap L=\{j\}$, the covariance estimator is, ${\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{J,t},\widehat{I}^{\prime}_{L,s})}=(p_{i,t}p_{j,t})^{-1}(1-p_{j,s})$.

[TABLE]

And, for all triangles $L=(i^{\prime},j,k_{s})$, for any $s<t$,

[TABLE]

where $D_{j,t}=\sum_{s<t}(1-p_{j,s})$. 3. (c)

If $k_{s}=i$ or $k_{s}=j$, then the ${\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{J,t},\widehat{I}^{\prime}_{L,s})}=0$.

We define $D_{i,t}$ and $D_{j,t}$ as the cumulative sum variables for edges $i$ and $j$ respectively, to keep track of previously sampled triangle indicators, that contain $i$ and $j$ respectively, at any time $s<t$. Note that for the new arriving edge $k$, we have $D_{k,t}=0$.

Estimating the ${\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{L,s},\widehat{I}^{\prime}_{J,t})}$.

For $s<t$, the estimate of the ${\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{L,s},\widehat{I}^{\prime}_{J,t})}$ is similar to the cases discussed previously. Thus, we adopt the same form in Theorem 4 (I). Note that while Theorem 4 (I) does not treat this case, it is straightforward to show that the estimator is also unbiased for the ${\mathop{\mathrm{Cov}}(\widehat{S}^{\prime}_{L,s},\widehat{I}^{\prime}_{J,t})}$. Hence, if $J\cap L=\{i\}$, the covariance estimator is,

[TABLE]

Thus, for all triangles $L=(i,j^{\prime},k_{s})$ and $s<t$,

[TABLE]

Similarly, if $J\cap L=\{j\}$, the covariance estimator is,

[TABLE]

And, for all triangles $L=(i^{\prime},j,k_{s})$ and $s<t$,

[TABLE]

Now, we add all the covariance increments for each edge as follows,

[TABLE]

Then, to update the cumulative variables for edges $i,j\in J=(i,j,k)$.

[TABLE]

We summarize all the variance and covariance computations in Algorithm 2, which is a supplementary to Algorithm 1 (in the case of triangle motifs).

Appendix C Ablation Study

To understand the effects of the various design choices in the proposed framework APS with shrinkage estimation, we conduct a thorough set of ablation study experiments. The proposed APS method provides a sampling framework that consists of two major parts: (1) Adaptive sampling with importance weights, and (2) James-Stein shrinkage estimation. Hence, there are several design choices to make, e.g., we could choose to use adaptive sampling without shrinkage estimation.

Results in Table 2 clearly show that shrinkage estimation significantly improves the performance of APS sampling. Another design choice is to use non-adaptive priority sampling where the edge weights/ranks are computed once at the time of sampling, and fixed during the rest of the streaming process. We conducted this experiment on the same datasets by using only the sampling weights assigned at arrival time (Line 14 in Alg 1) and fix it for the rest of the stream. We summarize the results in Table 4. For some graphs (e.g., soc-flickr), we observed that using non-adaptive weights in APS might perform better than using adaptive weights.

We conjecture this is due to the excessive variance of APS in the estimated count of the edges with small triangle counts, and can be observed in the tail of the distribution (see Figure 7). However, among all the design choices, the combination of (APS sampling + adaptive weights + shrinkage estimation) has the strongest regularization effect on the performance of graph sampling. We also observe that applying the shrinkage estimator to the non-adaptive sampling significantly improve the performance. These effects are demonstrated in Figures 6 and  7 which show the distribution of non-adaptive APS and adaptive APS respectively (with and without shrinkage estimation).

In summary, the results suggest that APS with shrinkage performs significantly better than related methods in previous work, and each of the design choices contributes to the final performance.

Appendix D Additional Plots

Appendix E Dataset Details

• •

soc-flickr: Crawl of the Flickr photo-sharing social network from May 2006. Nodes are users and edges represent that a user added another user to their list of contacts [19].

• •

soc-livejournal: LiveJournal is an online social community publishing platform, Nodes are users and edges are user-to-user links [35].

• •

soc-youtube: Youtube social network. Nodes are users and edges are user-to-user friendship links [35].

• •

wiki-Talk: Wikipedia network of user discussions from the inception of Wikipedia till January 2008. Nodes are Wikipedia users and edges are user-to-user edits of talk pages [31].

• •

web-BerkStan-dir: Web network where nodes represent webpages from Berkely and Stanford and edges represent hyperlinks among them [30].

• •

cit-Patents: The citation graph of US Patents includes all citations made by patents granted between 1975 and 1999 [29].

• •

soc-orkut-dir: Orkut online social network, where nodes represent users and edges represent user-to-user friendship links [35].