{'fraction_non_alphanumeric': 0.04658260559615042, 'fraction_numerical': 0.013968098378185707, 'mean_word_length': 4.0946544092611505, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 6, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
PRESTO: Simple and Scalable Sampling Techniques for the Rigorous Approximation of Temporal Motif Counts

The identification and counting of small graph patterns, called network motifs, is a fundamental primitive in the analysis of networks, with application in various domains, from social networks to neuroscience. Several techniques have been designed to count the occurrences of motifs in static networks, with recent work focusing on the computational challenges provided by large networks. Modern networked datasets contain rich information, such as the time at which the events modeled by the networks edges happened, which can provide useful insights into the process modeled by the network. The analysis of motifs in temporal networks, called temporal motifs, is becoming an important component in the analysis of modern networked datasets. Several methods have been recently designed to count the number of instances of temporal motifs in temporal networks, which is even more challenging than its counterpart for static networks. Such methods are either exact, and not applicable to large networks, or approximate, but provide only weak guarantees on the estimates they produce and do not scale to very large networks. In this work we present an efficient and scalable algorithm to obtain rigorous approximations of the count of temporal motifs. Our algorithm is based on a simple but effective sampling approach, which renders our algorithm practical for very large datasets. Our extensive experimental evaluation shows that our algorithm provides estimates of temporal motif counts which are more accurate than the state-of-the-art sampling algorithms, with significantly lower running time than exact approaches, enabling the study of temporal motifs, of size larger than the ones considered in previous works, on billion edges networks.


Introduction
The identification of patterns is a ubiquitous problem in data mining [8] and is extremely important for networked data, where the identification of small, connected subgraphs, usually called network motifs [18] or graphlets [23], have been used to study and characterize networks from various domains, including biology [16], neuroscience [3], social networks [28], and the study of complex systems in general [17]. Network motifs have been used as building blocks for various tasks in the analyses of networks across such domains, including anomaly detection [25] and clustering [4].
A fundamental problem in the analysis of network motifs is the counting problem [5,2], which requires to output the number of instances of the given topology defining the motif. This challenging computational problem has been extensively studied, with several techniques designed to count the number of occurrences of simple motifs, such as triangles [26,21,24] or sparse motifs [7].
Most recent work has focused on providing techniques for the analysis of large networks, which have become the standard in most applications. However, in addition to a significant increase in size, modern networks also feature a richer structure, in terms of the type of information that is available for their vertices and edges [6]. A type of information that has drawn significant attention in recent years is provided by the temporal dimension [11,12]. In several applications edges are supplemented with timestamps describing the time at which an event, modeled by an edge, occurred: for example, in the analysis of spreading processes in epidemics, nodes are individuals, an edge represents a physical interaction between two individuals, and the timestamp represents the time at which the interaction was recorded [22].
When studying motifs in temporal networks, one is usually interested in occurrences of a given topology whose edge timestamps all appear in a small time span [11,20]. Discarding the temporal information of the network, i.e. ignoring the timestamps, may lead to incorrect characterization of the system of interest, while the analysis of temporal networks can provide insights that are not revealed when the temporal information is not accounted • We provide tight and efficiently computable bounds to the number of samples required by our algorithms to achieve (multiplicative) approximation error ≤ ε with probability ≥ 1−η, for given ε > 0 and η ∈ (0, 1). Our bounds are obtained through the application of advanced concentration results (i.e., Bennett's inequality) for the sum of independent random variables.
• We perform an extensive experimental evaluation on real datasets, including a dataset with more than 2.3 billion edges, never examined before. The results show that on large datasets our algorithm PRESTO significantly improves over the state-of-the-art sampling algorithms in terms of quality of the estimates while requiring a small amount of memory.

Preliminaries
In this section we introduce the basic definitions used throughout this work. We start by formally defining temporal networks. Given (x, y, t) ∈ E, we say that t is the timestamp of the edge (x, y). For simplicity in our presentation we assume the timestamps to be unique, which is without loss of generality since in practice our algorithms also handle non-unique timestamps. We also assume the edges to be sorted by increasing timestamps, that is t 1 < · · · < t m . Given an interval or window [t B , t E ] ⊆ R we will denote |t E − t B | as its length.  [13]. (1c): sequences of edges of T that map topologically on the Bi-Fan motif, i.e., in terms of static graph isomorphism. For δ = 10 only the green sequence is a δ-instance of the Bi-Fan, since the timestamps respect σ and t − t 1 = 20 − 10 ≤ δ. The red sequences are not δ-instances, since they do not respect such constraint or do not respect the order in σ.
We are interested in temporal motifs 1 , which are small, connected subgraphs whose edge timestamps satisfy some constraints. In particular, we consider the following definition introduced in [20]. Note that a k-node -edge temporal motif M = (K, σ) is also identified by the sequence (x 1 , y 1 ), . . . , (x , y ) of edges ordered according to σ. Given a k-node -edge temporal motif M , k and are determined by V K and E K . We will therefore use the term temporal motif, or simply motif, when k and are clear from context. Given a temporal motif M , we are interested in counting how many times it appears within a time duration of δ, as captured by the following definition. Definition 2.3. Given a temporal network T = (V, E) and δ ∈ R + , a time ordered sequence S = (x 1 , y 1 , t 1 ), . . . , (x , y , t ) of unique temporal edges from T is a δ-instance of the temporal motif M = (x 1 , y 1 ), . . . , (x , y ) if: 1. there exists a bijection f on the vertices such that f (x i ) = x i and f (y i ) = y i , i = 1, . . . , ; and 2. the edges of S occur within δ time, i.e., t − t 1 ≤ δ.
Note that in a δ-instance of the temporal motif M = (K, σ) the edge timestamps must be sorted according to the ordering σ. See Figure 1 for an example.
Let U(M, δ) = {u : u is a δ-instance of M } be the set of (all) δ-instances of the motif M in T , denoted only with U when M and δ are clear from the context. Given a δ-instance u ∈ U(M, δ), we denote the timestamps of its first and last edge with t u 1 and t u , respectively. The count of M is C M (δ) = |U(M, δ)|, denoted with C M when δ is clear from the context. The notation used throughout this work is summarized in Table 1. We are interested in solving the following problem. Probability of selecting a sample that contains u ∈ U in PRESTO-A and PRESTO-E respectively.
Motif counting problem. Given a temporal network T , a temporal motif M = (K, σ), and δ ∈ R + , compute the count C M (δ) of δ-instances of motif M in the temporal network T .
Solving the motif counting problem exactly may be infeasible for large networks, since even determining whether a temporal network contains a simple star motif is NP-hard [13]. State-of-the-art exact techniques [15,20] require exponential time and memory in the number of edges of the temporal network, which renders them impractical for large temporal networks. We are therefore interested in obtaining efficiently computable approximations of motif counts, as follows.
We call an algorithm that provides such guarantees an (ε, η)-approximation algorithm.

Related Work
Various definitions of temporal networks and motifs have been proposed in the literature; we refer the interested reader to [10,12,14]. Here we focus on those works that adopted the same definitions used in this work. The definition of temporal motif we adopt was first proposed by Paranjape et al. [20], which provided efficient exact algorithms to solve the motif counting problem for specific motifs. Such algorithms are efficient only for specific motif topologies and do not scale to very large datasets. An algorithm for the counting problem on general motifs has been introduced by Mackey et al. [15]. Their algorithm is the first exact technique allowing the user to enumerate all δ-instances u ∈ U without any constraint on the motif's topology. The major back-draw of such algorithm is that it may be impractical even for moderately-sized networks, due to its exponential time complexity and memory requirements.
Liu et al. [13] proposed the first sampling algorithm for the temporal motif counting problem. The main strategy of [13] is to partition the time interval containing all the edges of the network into non overlapping and contiguous windows of length cδ (i.e., a grid-like partition), for some c > 1. The partition is then randomly shifted (i.e., the starting point of the first window may not coincide with the smallest timestamp of the network). The edges in each partition constitute the candidate samples to be analyzed using an exact algorithm. An importance sampling scheme is used to sample (approximately) r windows among the candidates, with a window being selected with probability proportional to the fraction of edges it contains. The estimate for each sampled window is obtained by weighting each δ-instance in the window, and the estimates are averaged across windows.
This procedure is repeated b times to reduce the variance of the estimate. While interesting, this partition-based approach prevents such algorithm to provide (ε, η)-guarantees (see Section 3.1).
Recently Wang et al. [29] proposed an (ε, η)-approximation algorithm for the motif counting problem. Their approach selects each edge in T with a user-provided probability p. Then for each selected edge e = (x, y, t), the algorithm collects the edges with timestamps in the edge-centered window [t − δ, t + δ], of length 2δ, computes on these edges all the δ-instances u ∈ U containing e, weights the instances, and combines the weights to obtain the final estimate. From the theoretical point of view, the main drawback of this approach is that in order to achieve the desired guarantees one has to set p ≥ 1/(1 + ηε 2 ), resulting in high values of p (i.e., almost all edges are selected) for reasonable values of η and ε (e.g., p > 0.97 for η = 0.1 and ε = 0.5). In addition, such approach is impractical on very large datasets, mainly due its huge memory requirements, and does not provide accurate estimates (see Section 4).

PRESTO: Approximating Temporal Motif Counts with Uniform Sampling
We now describe and analyze our algorithm PRESTO (apPRoximating tEmporal motifS counTs with unifOrm sampling) for the motif counting problem. We start by describing, in Sec. 3.1, the common strategy underlying our algorithms. We then briefly highlight the differences between PRESTO and the existing sampling algorithms for the counting problem. In Sec. 3.2 and Sec. 3.3 we present and analyze two variants of the common strategy introduced in Sec. 3.1. We conclude with the analysis of PRESTO's complexity in Sec. 3.4.

PRESTO: General Approach
The general strategy of PRESTO is presented in Algorithm 1. Given a network T , PRESTO collects s samples, where each sample is obtained by gathering all edges e ∈ E with timestamp in a small random window [t r , t r + cδ] of length cδ, using SampleWindowStart(E) (Lines 2-5) to select t r (i.e., the starting time point of the window) through a uniform random sampling approach. For each sample, an exact algorithm is then used to extract all the δ-instances of motif M in such sample (Line 6). The weight of each δ-instance extracted is computed with the call ComputeWeight(u) (Line 8), and the sum of all such weights constitutes the estimate provided by the sample (Line 9). The weights account for the probability of sampling the δ-instances in the sample, making the final estimate unbiased. The final estimate produced in output is the average of the samples' estimates (Lines 10-11). Note that the for cycle of Line 1 is trivially parallelizable. The two variants of PRESTO we will present differ in the way i) SampleWindowStart(E) and ii) ComputeWeight(u) are defined.
Differently from [13], our algorithm PRESTO does not partition the edges of T in non overlapping windows, and relies instead on uniform sampling. We recall (see Section 2.1) that in [13], after computing all the non overlapping intervals defining the candidate samples, there may be several δ-instances u ∈ U that cannot be sampled (i.e., all δ-instances having t u 1 in window j and t u in window j + 1). This significantly differs from PRESTO, which instead samples at each iteration a small random window from [t 1 , t m ] without restricting the candidate windows, allowing to sample any δ-instance u ∈ U(M, δ) at each iteration. This enables us to provide stronger guarantees on the quality of the output, since each δ-instance has a non-zero probability of being sampled at each step, leading to (ε, η)-approximation guarantees.
Differently from the work of Wang et al. [29] PRESTO samples temporal windows of length cδ and does not follow the edge-centric approach (i.e., sampling temporal edges with a user-provided probability p) of [29]. In addition, the approach of [29] collects edges in temporal-windows of length 2δ (see Section 2.1), while in PRESTO the window size is controlled by the parameter c, which, when fixed to be < 2, leads to much more scalability than [29] while requiring less memory (see Section 4).

PRESTO-A:
Sampling among All Windows In this section we present and analyze PRESTO-A, our first (ε, η)-approximation algorithm obtained by specifying i) how the starting point t r of the temporal window defining a sample T i is chosen (function SampleWindowStart(E) in Line 2) and ii) how the weight w(u) of a δ-instance u in a sample is computed (ComputeWeight(u), Line 8).
The starting point t r of sample T i is sampled uniformly at random in the interval [t − cδ, t m− ] ⊆ R, where we recall = |E K | and m = |E|. Regarding the choice of the weight w(u) for each instances u ∈ S i , PRESTO-A considers w(u) = 1/p u , with p u being the probability of u to be in S i , that is is the length of the interval in which t r must be chosen for u to be in We now present the theoretical guarantees of PRESTO-A and give efficiently computable bounds for the sample size s needed for the (ε, η)-approximation to hold. Recall the definition of U(M, δ), which is the set of δinstances of M in T . Let u be an arbitrary δ-instance of motif M , and let T i be an arbitrary sample obtained by PRESTO-A at its i-th iteration. We define the following set of indicator random variables, for u ∈ U and for i = 1, . . . , s: We therefore have the following result.
Proof. Considering the linearity of expectation and the definition of the variables X i and X i u , i = 1, . . . , s, u ∈ U, we have: Combining with the linearity of expectation to C M the statement follows.
The following result provides a bound on the variance of the estimate of C M provided by C M .
Proof. We start by observing that, Var 1 Var(X i ) by the mutual independence of the variables X 1 , . . . , X s . We observe that Var( and (2.) from the definition of p u2 . Based on the above, we can bound the variance of the variables X i , i = 1, . . . , s, as follows, Var( Note that such bound does not depend on the index i = 1, . . . , s, thus substituting in the original summation we obtain, We now present a first efficiently computable bound on the number of samples s required to have that C M is a relative ε-approximation of C M with probability ≥ 1 − η. Such bound is based on the application of Hoeffding's inequality (see [19]), an advanced but commonly used technique in the analysis of probabilistic algorithms. Let us first introduce Hoeffding's inequality, then the bound on the sample size we derive based on the above is as follows, Theorem 3.2. Given ε ∈ R + , η ∈ (0, 1) let X 1 , . . . , X s be the random variables associated with the counts computed at iterations 1, . . . , s, respectively, of PRESTO-A. If s ≥ , we want to show that for s as in statement it holds P r ≤ η. Observe that the variables X 1 , . . . , X s , are mutually independent. And, by Lemma Then for X i , i = 1, . . . , s, the following holds: where (1.) follows from the definition of p u , by considering each of the variables X i u = 1, and noting that for each (c−1) 2 δ 2 and applying Theorem 3.1 to bound P r we obtain, by the choice of s as in statement.
We now show that by using Bennett's inequality (see [33]), a more advanced result on the concentration of the sum of independent random variables, we can derive a bound that is much tighter than the above, while still being efficiently computable. We first state Bennett's inequality, . . , s and for any t ≥ 0, the following holds One of the difficulties in applying such result to our scenario is that we know only an upper bound to the variance of the random variables X i , i = 1, . . . , s. We discuss why Bennett's inequality can be applied in such situation in Supplement A. We now state our main result. Theorem 3.4. Given ε ∈ R + , η ∈ (0, 1) let X 1 , . . . , X s be the random variables associated with the counts computed at iterations 1, . . . , s, respectively, of PRESTO-A. If s ≥ we want to show that for s as in the statement it holds P r ≤ η. In order to apply Bennett's bound, note that 1 applying Theorem 3.3, with the upper boundv to v we obtain, by substituting the quantities above and by the choice of s as in statement, which concludes the proof.
Note that the bound in Theo. 3.4 is significantly better than the one in Theo. 3.2, since the former has a quadratic dependence from ∆ T,1 while the latter enjoys a linear dependence from ∆ T,1 . Furthermore, differently from the bounds in [29], our bounds depend on characteristic quantities of the datasets (∆ T,1 ) and of the algorithm's input (c, δ). Therefore we expect our bounds to be more informative, and also possibly tighter than the bounds of [29] (with the improvement due, at least in part, to the use of Bennett's inequality, while [29] leverages on Chebyshev's inequality, which usually provides looser bounds [19]). We conclude by giving an intuition on the bound in Theorem 3.4, observe in fact that the bound is directly proportional to the quantity ∆ T,1 /((c−1)δ), such quantity (approximatively) corresponds to the fraction between the timespan of the temporal network T and the length of each window we draw, thus such bound suggests (in a very intuitive way) that as long as the windows we draw have a larger length (i.e., δ or c become larger) then less samples are needed for PRESTO-A to provide concentrate estimates within the desired theoretical guarantees. Note that windows with a larger duration will require a higher running time to be processed, as we will discuss in Section 3.4.

PRESTO-E: Sampling the Start from Edges
In this section we present and analyse our second (ε, η) approximation algorithm, PRESTO-E, obtained with a different variant of the general strategy of Algorithm 1.
PRESTO-E selects the starting point t r (Line 2 in Alg. 1) of sample T i only from the timestamps of edges of T . In particular, t r is chosen uniformly at random in {t 1 , . . . , t last } ⊆ {t 1 , . . . , t m }, where t last = min{t : (x, y, t) ∈ E ∧ t ≥ t m − cδ}. When the edges of T are far from each other or have a skewed distribution in time, or occur mostly in some well-spaced subsets {[t a , t b ] ⊆ [t 1 , t m ] : t a ≤ t b }, we expect PRESTO-E to collect samples with more δ-instances than PRESTO-A (which samples t r uniformly at random almost from the entire time interval [t 1 , t m ] but without restricting to existing edges). Similarly to PRESTO-A, the weight w(u) (Line 8 of Alg. 1) is computed by w(u) = 1/p u , wherep u is the probability that u ∈ U is in the set S i . Due to the (random) choice of t r ,p u =r u /∆ T,2 , where ∆ T,2 = |{t : (x, y, t) ∈ E ∧ t ∈ [t 1 , t last ]}| is the number of possible choices for t r (if t last = t m then ∆ T,2 = m) andr u = |{t : (x, y, t) ∈ E ∧ t ∈ [max{t 1 , t u − cδ}, min{t last , t u 1 }]}| is the number of choices for t r such that u ∈ U is in S i , i = 1, . . . , s.
Similarly as done for PRESTO-A we prove that PRESTO-E provides an unbiased estimate of C M , with bounded variance (the proof is in Supplement B). By applying Bennett's inequality, we derive the following. Theorem 3.5. Given ε ∈ R + , η ∈ (0, 1) let X 1 , . . . , X s be the random variables associated with the counts computed at iterations 1, . . . , s, respectively, of PRESTO-E. If s ≥ (∆ T ,2 −1) (1+ε) ln(1+ε)−ε ln 2 η , then it holds As for PRESTO-A, by using Bennett's inequality we obtain a linear dependence between the sample size s and ∆ T,2 , while commonly used techniques (e.g., Hoeffding's inequality) lead to a quadratic dependence.

PRESTO: Complexity Analysis
In this section we analyse the time complexity of PRESTO's versions introduced in the previous sections.
Note that our algorithms employ an exact enumerator as subroutine. For the sake of the analysis we consider the complexity when the subroutine is the algorithm by Mackey et al. [15], which we used in our implementation. Let us first start with a definition, given a temporal network T = (V, E), we denote with dim(T ) the set {t i |(x, y, t i ) ∈ E, i = 1, . . . , m}. The algorithm in [15] has a worst case complexity O(mκ ( −1) ) when executed on a motif with edges and a temporal network with m temporal edges, whereκ is the maximum number of edges within a window of length δ, i.e.,κ = max{|S(t)| : Observe that such complexity is exponential in the number of edges of the temporal motif. Obviously our algorithms benefit from any improvement to the state-of-the-art exact algorithms.
Complexity of PRESTO-A. The worst case complexity is O(smκ ( −1) + sC * M ) when executed sequentially, wherem is the maximum number of edges in a window of length cδ, i.e.,m = max{|S(t)| :

Experimental evaluation
In this section we present the results of our extensive experimental evaluation on large scale datasets. To the best of our knowledge we consider the largest dataset, in terms of number of temporal edges, ever used for the motif counting problem with more than 2.3 billion temporal edges. After describing the experimental setup and implementation (Sec. 4.1), we first compare the quality of the estimated counts provided by PRESTO with the estimates from state-of-the-art sampling algorithms from Liu et al. [13], and Wang et al. [29] (Sec. 4.2). Then we compare the memory requirements of the algorithms (Sec. 4.3), which may be a limiting factor for the analysis of very large datasets. Additional experiments on PRESTO's running time comparison with the exact algorithm by Mackey et al. [15], and with the current state of-the art sampling algorithms are in Supplement G. While results showing the linear scalability of PRESTO's parallel implementation with the available number of processors can be found in Supplement H.

Experimental Setup and Implementation
The characteristics of the datasets we considered are in Table 2. EquinixChicago [1] is a bipartite temporal network that we built, its description can be found in Supplement C. Descriptions of the other datasets are in [20,9,13].
We implemented our algorithms in C++14, using the algorithm by Mackey et al. [15] as subroutine for the exact solution of the counting problem on samples. We ran all experiments on a 64-core Intel Xeon E5-  2698 with 512GB of RAM machine, with Ubuntu 14.04. The code we developed, and the dataset we build are entirely available at https://github.com/VandinLab/PRESTO/, links to additional resources (datasets and other implementations) are in Supplement D. We tested all the algorithms on the motifs from Figure 2, that are a temporal version of graphlets from [23], often used in the analysis of static networks. We considered motifs with at most = 4 edges since the implementation by Wang et al. [29] does not allow for motifs with a higher in input. Since the EC dataset is a bipartite network, it can contain all but the motifs marked with a rectangle in Fig. 2. Therefore, for such dataset we did not consider the motifs in the rectangles of Fig. 2. We compared our algorithms PRESTO-A and PRESTO-E, collectively denoted as PRESTO, with 2 baselines: the sampling algorithm by Liu et al. [13], denoted by LS and the sampling algorithm by Wang et al. [29], denoted by ES.

Quality of Approximation
We start by comparing the approximation provided by PRESTO with the current state-of-the-art sampling algorithms. To allow for a fair comparison, since the algorithms we consider depend on parameters that are not directly comparable and that influence both the algorithm's running time and approximation, we run all the algorithms sequentially and measure the approximation they achieve by fixing their running time (i.e., we fix the parameters to obtain the same running time). While PRESTO can be trivially modified in order to work within a user-specified time limit, the same does not apply to LS or ES. Thus, to fix the same running time for all algorithms we devised the following procedure: i) for a fixed dataset and δ, we found the parameters of LS such its runtime is at least one order of magnitude 2 smaller than the exact algorithm by Mackey et al. [15] on motif G1-1 ii) we run LS on all the motifs on the fixed dataset with the parameters set as from i) and measure the time required for each motif; iii) we run ES such that it requires the same time 3 as LS; iv) we ran PRESTO-A and PRESTO-E with c = 1.25 and a time limit that is the minimum between the running time of LS and the one of ES. The parameters resulting from such procedure for all algorithms are in Supplement E. Further discussion on how to set c when running PRESTO is in Supplement F.
For each configuration, we ran each algorithm ten times, with a limit of 400 GB of (RAM) memory. We  computed the so called Mean Average Percentage Error (MAPE) after removing the best and worst results, to remove potential outliers. MAPE is defined as follows: let C M be the estimate of the count C M produced by one run of an algorithm, then the relative error of such estimate is |C M − C M |/C M ; the MAPE of the estimates is the average of the relative errors in percentage. In addition, we also computed the standard deviations of the estimates in output to the algorithms for each configuration, which will be reported to complement the mean values, small standard deviations correspond to more stable estimates of the algorithms around their MAPE values. Table 3 shows the MAPE values of the sampling algorithms on datasets SO, and BI. On the SO dataset we observe that both versions of PRESTO provide more accurate estimates over current state-of-the-art sampling algorithms. PRESTO-A provides the best approximation on 6 out of 12 motifs, with PRESTO-E providing the best results in all the other cases except for motif G2-1. Furthermore, both variants of PRESTO improve on all motifs w.r.t. LS and on 11 out of 12 motifs w.r.t. ES, with the error from ES being more than three times the error of PRESTO on such motifs. We note PRESTO-A and PRESTO-E achieve similar results, supporting the idea that similar performances are obtained when the network's timestamps are distributed evenly, that is the case for such dataset (see Supplement I).
For the BI dataset the results are similar to the ones for the SO dataset. PRESTO-A achieves the lowest approximation error on 7 out of 12 motifs with PRESTO-E achieving the lowest approximation error on all the other motifs but G2-2. PRESTO-A achieves better estimates than LS on all motifs and on 11 out of 12 motifs w.r.t. ES, PRESTO-E behaves similarly, except that it improves over LS on 10 out of 12 motifs.  Figure 2, "" denotes out of memory.  Table 4 instead shows the results on the RE and EC datasets. On the former, PRESTO-E achieves the lowest approximation error on 7 out of 12 motifs and it improves over PRESTO-A for most of the motifs. This is not surprising, since [29] showed that RE has a more skewed edge distribution (see also the discussion in Supplement I), a scenario for which we expect PRESTO-E to improve over PRESTO-A (see Section 3.3). Nonetheless, PRESTO-A achieves lower approximation error over 6 motifs when comparing to LS and on 10 motifs w.r.t. ES, while PRESTO-E achieves the best approximation error on most of the motifs, improving over LS on 8 motifs and over ES on 10 motifs.
Finally, we discuss the results on the EC dataset (again in Table 4), which is a 2.3 billion edges bipartite temporal network and thus it does not contain the motifs marked with a rectangle in Figure 2. Note that the results for ES are missing, since ES did not complete any run with 400GB of memory on the motifs tested. (We discuss the high memory usage of ES in Section 4.3). On such dataset both variants of PRESTO perform better than LS on 7 out 8 motifs and achieve the lowest approximation error on 4 out 8 motifs each. Interestingly, the standard deviations of PRESTO are often similar or lower than the one of LS, which is another of the advantages of our algorithms. We observed similar trends, and even smaller variances for PRESTO across all the experiments on the other datasets (see Tables 3, 4).
These results, coupled with the theoretical guarantees of Section 3, show that PRESTO outperforms the stateof-the-art sampling algorithms LS and ES for the estimation of most of the motif counts on temporal networks. Based on the results we discussed, PRESTO-E seems also to usually provide similar estimates to PRESTO-A, while significantly improving over PRESTO-A when the network edges are not uniformly distributed such the EC dataset.

Memory usage
In this section we discuss the RAM usage of the various sampling algorithms. The results are in Table 5, which shows for each algorithm the minimum and maximum amount of (RAM) memory used over all motifs in Figure 2 (collected on a single run). Both PRESTO's versions have lower memory requirements than ES, and equal or lower memory requirements than LS, on all datasets. ES ranges from requiring 6.5× up to 8.5× the memory used by PRESTO, making ES not practical for very large datasets such as EC (where ES did not completed any of the runs since it exceeded the 400GB of memory allowed for its execution). While LS and PRESTO have similar memory requirements, LS displays greater variance on BI and RC, which may be due to the fact that choosing windows with more edges, as done by LS, may require more memory to run the exact algorithm as subroutine.
PRESTO is, thus, a memory efficient algorithm and hence a very practical tool to tackle the motif counting problem on temporal networks with billions edges.

Conclusions
In this work we introduced PRESTO, a simple yet practical algorithm for the rigorous approximation of temporal motif counts. Our extensive experimental evaluation shows that PRESTO provides more accurate results and is much more scalable than the state-of-the-art sampling approaches. There are several interesting directions for future research, including developing rigorous techniques for mining several motifs simultaneously and developing efficient algorithms to maintain estimates of motif counts in adversarial streaming scenarios.

Supplementary material.
A Note on Bennett's Inequality In this section we will discuss how to apply Bennett's inequality when the exact variance term, i.e. v, from Theorem 3.3 is unknown and we only have an upper boundv to such term (which is our case in practice).
Let us consider the function f (v) = αv h(γ/v) with α = s/B 2 , γ = tB where h(·) is the same function h(·) appearing in Bennett's inequality (Theorem 3.3). With a somehow tedious analysis one can prove that the function f (v) is continuous and monotonic non increasing for every v, γ > 0 thus by definition, . Observe that we can write the right-hand side of the Bennett's bound as 2 exp(−f (v)), now is clear that given v 1 ≥ v 2 > 0 it holds 2 exp(−f (v 2 )) ≤ 2 exp(−f (v 1 )) by the property of negative exponential functions. This observation allows us to apply Bennett's inequality using an upper boundv to v.
By the above observation we also get the intuition that a smaller gap betweenv, the upper bound to the variance, and v the actual variance will result in a tighter bound to the sample size s (i.e., Bennett's inequality will provide more strict bounds). In our work we focused on a tradeoff between the sharpness of the boundv and its efficient computability in practice, since sharper bounds onv (than the one we provide) could be proved but they are much harder to compute and thus may not be of practical interest. We still leave for future works the space of improving the bounds in this work.

B PRESTO-E -Theoretical Guarantees
In this section we will prove that, for an appropriate choice of the number s of samples, PRESTO-E is an (ε, η)approximation algorithm. The analysis is similar to the one of PRESTO-A.
We define the indicator random variables, for each u ∈ U(M, δ) and for i = 1, . . . , s:X i u = 1 if u ∈ S i , 0 otherwise. Each variableX i u , i = 1, . . . , s , u ∈ U is a Bernoulli random variable for which it holds With a similar analysis to the one for PRESTO-A, the following results hold.
By applying Bennett's inequality (Theorem 3.3), we derive the following bound on the sample size s which is also reported in the main text.

C Description of the Equinix-Chicago Dataset
In this section, we briefly describe the Equinix-Chicago used in our experimental evaluation. We built the Equinix-Chicago dataset 4 starting from the internet packets collected by two monitors situated at Equinix data center, which captured the packets from Chicago to Seattle and vice versa on the 17 February 2011 for an observation period of 62 minutes. Each edge (x, y, t) represents a packet sent at time t, with microseconds precision, from the IPv4 address x to the IPv4 destination y 5 . Such network is thus a bipartite temporal network since edges, i.e., packet exchanges, occur only between nodes belonging to different cities, with no interactions captured among nodes of the same city. The final temporal network has more than 2.325 billion temporal edges and 11.157 million nodes. This dataset is publicly available at https://github.com/VandinLab/PRESTO/.

D Reproducibility -Links to External Resources
For reproducibility purposes, in this section we report the links to the datasets and implementations of LS and ES we used.

E Reproducibility -Parameters used in the experimental evaluation
We now discuss discuss the choice of the parameters resulting from our procedure for comparing the different sampling algorithms in Section 4.2.
Note that our algorithms PRESTO-A and PRESTO-E have only 1 tuning parameter in the experimental setting of Section 4.2, since δ is provided in input (i.e., c, for which we also discuss how to set it to obtain the maximum efficiency in Supplement F). The same holds for ES which has the paramer p, that controls the fraction of sampled edges (see Section 2.1). While LS has instead non trivial dependencies among the parameters b and r which are not entirely easy to set when it comes to different values of c and δ also due to the lack of a qualitative theoretical analysis relating the parameters to the approximation quality of the results. Table 6 reports all the parameters used on the different datasets. Given the procedure for running PRESTO's versions (described in Section 4.2), their sample size s is not deterministic, since only the execution time is fixed. In Table 6 we report the average sample sizes s PRESTO-A and s PRESTO-E , over all the runs on all the motifs on the given network, of PRESTO-A and PRESTO-E respectively. Recall that the values of p for ES where adapted to each motif on each configuration following the approach described in Section 4.2 selecting such values from a grid of values that is, we selected the best p in order to have a final running time of ES to be close to the one of LS choosing p from the following grid [10 −6 , 5 · 10 −6 , 10 −5 , 5 · 10 −5 , 10 −4 , 5 · 10 −4 , 10 −3 , 2 · 10 −3 , 5 · 10 −3 , 10 −2 ]. Thus, we report in Tab. 6 the range from which p was chosen in practice as result of our procedure for fixing the parameters.  The major advantages of PRESTO, come from the efficiency in time, small memory usage and from the scalability it achieves, which coupled with the theoretical guarantees provided and the precise estimates achieved in practice, make our algorithm a fundamental tool when analysing large temporal networks, it is thus important to set PRESTO's parameter (i.e., c) correctly to exploit such features at their best. Our versions of PRESTO have only the parameter c to be setted, since δ is part of the problem's input. We now discuss the importance of keeping c small (e.g., c < 2). To do this, we consider a concrete example to show that a small c leads to more efficiency and scalability. Suppose we want to approximate the count of motif G5-1 from Figure 2 on EC dataset 6 and we want to cover with samples a fixed length of the interval [t 1 , t m ] over the s iterations, i.e.,  Table 7 presents the running times for various values of c with fixed δ = 10 5 , V = 1.25 · 10 8 . When c decreases there is a substantial decrease in running time, in fact with c = 10 the running time is more than 3.9× the running time with c = 1.25. This should not surprise the reader, since such results is in accordance with our complexity analysis (see Section 3.4) which shows that larger samples in terms of edges (as obtained by sampling windows with a higher c) result in a higher running time.
Thus the best choice in order to set c is to keep c < 2 based on all the experiments we performed, to exploit efficiency and scalability at the maximum of their level.

G Running time comparison
In this section we discuss the running times under which we obtained the results in Section 4.2.
We measured the running time of PRESTO-A, PRESTO-E, LS, and ES over all experiments (i.e., datasets and motifs) described in Section 4.2, we also measured the running time of BT over three runs on the same configurations (i.e., dataset ad δ). For all algorithms we then computed the geometric average of their running times for each pair (dataset, motif). The results are shown in Figure 3. As expected, PRESTO runs with at least an order of magnitude less than BT. The reader may appreciate how the results in Section 4.2 were obtained with PRESTO's running time bounded by the minimum running time among ES and LS (fixed the motif and the dataset). Furthermore we note that on several motifs (G10-1, G9-1) on almost all datasets, even if we set p small, i.e., 10 −5 , ES runs within a higher time than LS and thus than PRESTO, we believe that this may be due to ES's procedure to explore the search space which may not suite such motifs. To this end we find that PRESTO is not affected by the motif topology. We present the results on PRESTO's running time against BT on the EC dataset in Figure (3d) which corresponds to the times under which we obtained the results in Table ??. Note that PRESTO runs orders of magnitude faster than BT while still providing accurate estimates. For example on G3-2 PRESTO achieves a 3% relative error (see Table 4) while being 53× faster than BT.
Coupled with the results of Sections 4.2 and 4.3, these aspects show that PRESTO is an efficient algorithm to obtain accurate estimates of motif counts from large temporal networks, especially when the running time to obtain an estimate is very limited, this may be the case on very large datasets.  Running Time (sec) Running Time (sec) Temporal Motif H Scalability of PRESTO's parallel implementation As discussed in Section 3.1, PRESTO can be easily parallelized by executing the s iterations of Algorithm 1 in parallel, in this section we discuss such implementation. We implemented such strategy through a thread pooling design pattern. We present the results of the multithread version of PRESTO-A. Results for PRESTO-E are similar. We present the results on EC for motif G5-1 (see Figure 2), which has the highest running time over all the experiments in Section 4.2 (see Figure (3d)), similar results are observed on other motifs and datasets. Figure 4 shows the speed-up T 1 /T j where T j , j ∈ {1, 2, 4, 8, 16} is the geometric average of the running time over five runs (with fixed parameters) executing PRESTO-A with j threads. In Figure 4 (left) we fixed c and δ as in Supplement E (Table 6), and we show how the speed-up changes by varying the sample size s. We observe a nearly linear speed-up for up to 8 threads, and a speed-up for 16 threads that is slightly less than linear. In Figure 4 (right) we fix c = 1.25 and s = 1000 and we show how the speed-up changes across different values of δ ∈ {75, 100, 125} · 10 3 . As we can see there is a linear and slightly superlinear behaviour up to 8 threads, with the performances becoming slightly less than linear when using 16 threads, similarly to the results obtained for the tests with fixed δ. This behaviour may be due to design pattern or to the time needed to process each sample, intuitively if the time to process each sample is not large, then there may be much time spent by the algorithm in synchronization when using a large number (≥ 16) of threads and this may be the case on EC dataset.
All together, the results of our experimental evaluation show that PRESTO obtains up to 500× of speed-up over exact procedure on the tests we performed while providing accurate estimates through its parallel implementation (53× faster through its sequential implementation, while its parallel implementation with 16 threads executes with a speedup of more than 10× over the sequential one), enabling the efficient analysis of billion edges networks.

I Edge Timestamps distribution -Skewed (Dataset RE) vs Uniform (Dataset SO)
In this section we compare the edge distributions, with respect to the timestamps, of the datasets SO and RE used in our experimental evaluation in Section 4. To obtain such distribution we computed for each timestamp t i , i = 1, . . . , m s.t. ∃(u, v, t i ) ∈ E, the number of edges in the forward window of length cδ starting from t i , i.e., we computed |{(x, v, t) ∈ E|t ∈ [t i , t i + cδ]}| for each t i , i = 1, . . . , m, where we kept the parameters c and δ as in Table 6. This corresponds, for example, to the number of edges in each sample candidate of PRESTO-E. The results are shown in Figure 5. We first discuss the distribution of the timestamps of the edges in the RE dataset ( Figure 5 (left)), we observe immediately that such distribution is very skewed. In particular, the number of edges in the windows of length cδ starting from edges with timestamps in the first quarter of [t 1 , t m ], i.e., t < 1.25 · 10 9 , do not exceed 10 3 edges. While most of the temporal windows starting at timestamps in the middle interval of [t 1 , t m ], i.e., 1.25 · 10 9 < t < 1.4 · 10 9 , have much more edges ranging from 10 3 to 10 5 (several orders of magnitude edges more than windows from the first quarter of the network's timespan). In the last quarter of [t 1 , t m ] instead (t > 1.4·10 9 ) the dataset presents very sparse edges with most of the windows (which are sample candidates in PRESTO) not exceeding 10 2 edges, note that the windows starting in the middle section of [t 1 , t m ] have at least one order of magnitude additional edges (i.e., [10 3 , 10 5 ] edges). This shows how the timestamps of the edges in RE have a very skewed distribution. Figure 5 (right) instead shows the distribution of the timestamps of the edges on the dataset SO. As we can see, except for a very brief transient state, the timestamps are almost uniformly distributed on [t 1 , t m ], with the windows of length cδ centered at each timestamp t i , i = 1, . . . , m having almost the same number of edges,