{'fraction_non_alphanumeric': 0.03146250509228793, 'fraction_numerical': 0.01024725016452007, 'mean_word_length': 4.478454935622318, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Structurally Smoothed Graphlet Kernel

A commonly used paradigm for representing graphs is to use a vector that contains normalized frequencies of occurrence of certain motifs or sub-graphs. This vector representation can be used in a variety of applications, such as, for computing similarity between graphs. The graphlet kernel of Shervashidze et al. [32] uses induced sub-graphs of k nodes (christened as graphlets by Przulj [28]) as motifs in the vector representation, and computes the kernel via a dot product between these vectors. One can easily show that this is a valid kernel between graphs. However, such a vector representation suffers from a few drawbacks. As k becomes larger we encounter the sparsity problem; most higher order graphlets will not occur in a given graph. This leads to diagonal dominance, that is, a given graph is similar to itself but not to any other graph in the dataset. On the other hand, since lower order graphlets tend to be more numerous, using lower values of k does not provide enough discrimination ability. We propose a smoothing technique to tackle the above problems. Our method is based on a novel extension of Kneser-Ney and Pitman-Yor smoothing techniques from natural language processing to graphs. We use the relationships between lower order and higher order graphlets in order to derive our method. Consequently, our smoothing algorithm not only respects the dependency between sub-graphs but also tackles the diagonal dominance problem by distributing the probability mass across graphlets. In our experiments, the smoothed graphlet kernel outperforms graph kernels based on raw frequency counts.


INTRODUCTION
In this paper, we are interested in comparing graphs by computing a kernel between graphs [39]. Graph kernels are popular because many datasets from diverse domains such as bio-informatics [3,30], chemo-informatics [2], and web data mining [41] naturally can be represented as graphs. Almost all graph kernels (implicitly or explicitly) represent a graph as a (normalized or un-normalized) vector which contains the frequency of occurrence of motifs or subgraphs 1 . The key idea here is that well chosen motifs can capture the semantics of the graph structure while being computationally tractable. For instance, counting walks in a graph leads to the random walk graph kernel of Borgwardt et al. [3] (see Vishwanathan et al. [39] for an efficient algorithm for computing this kernel). Other popular motifs include subtrees [31], shortest paths [5], and cycles [16]. Of particular interest to us are the graphlet kernels of Shervashidze et al. [32]. The motif used in this kernel is the set of unique sub-graphs of size k, which were christened as graphlets by Przulj [28].
Observation 1: Computing meaningful graphlet kernels that have high discriminative ability requires a careful selection of k.
If k is small, then the number of unique graphlets is small (i.e., the length of the feature vector is small). See Figure 1.
Consequently the feature vector does not provide meaningful discrimination between two graphs. On the other hand, if k is large then a) the set of unique graphlets grows exponentially (i.e., the feature vector is very high dimensional) but b) only a small number of unique graphlets will be observed in a given graph (i.e., the feature vector is very sparse). Moreover, the probability that two graphs will contain a given large sub-graph is very small. Consequently, a graph is similar to itself but not to any other graph in the training data. This is well known as the diagonal dominance problem in the machine learning community [17], and the resulting kernel matrix is close to the identity matrix. In other words, the graphs are orthogonal to each other in the feature space. However, it is desirable to use large values of k in order to gain better discriminative ability. One way to circumvent the diagonal dominance problem is to view the normalized graphlet-frequency vector as estimating a multinomial distribution, and use smoothing.

Observation 2:
The normalized graphlet-frequency vector exhibits power-law behavior, especially for large values of k. In other words, a few popular graphlets occur very frequently while a vast majority of graphlets will occur very rarely. Put another way, a few graphlets dominate the distribution. To see this, we randomly sampled a graph from six benchmark datasets (details of the datasets can be found in Section 6), and exhaustively computed occurrences of all graphlets of size k = 8 and plotted the resulting histogram on a log-log scale in Figure 2. As can be seen, the frequencies are approximately linear in the log-log scale which indicates power law behavior. Therefore, any smoothing technique that we use on the normalized graphlet-frequency vector must respect this power-law behaviour.   to each other. While many such relationships can be derived, we will work with perhaps the simplest one which is depicted in Figure 4. Here, we construct a directed acyclic graph (DAG) with the following property: a node at depth k denotes a graphlet of size k. Given a graphlet g of size k and other graphlet g of size k + 1 we add an edge from g to g if, and only if, g can be obtained from g by deleting a node of g . This shows that graphlets of size k have a strong relationship to graphlets of size k + 1 and one must respect this relationship when deriving a smoothing technique.
Our contributions in this paper are as follows. First, we propose a new smoothing technique for graphlets which is inspired by Kneser-Ney smoothing [20] used for language models in natural language processing. Our model satisfies the desiderata that we outlined above, that is, it respects the power law behavior of the counts and yet takes into account the structure of the space of graphlets. Second, we provide a novel Bayesian version of our model that is extended from the Hierarchical Pitman-Yor process of Teh [35]. Unlike the traditional Hierarchical Pitman-Yor Process (HPYP) where the base distribution is given by another Pitman-Yor Process (PYP), in our case it is given by a transformation of a PYP that is guided by the structure of the space. Third, we perform experiments to validate and understand how smoothing affects the performance of graphlet kernels.
The structure of the paper is as follows. In Section 2, we discuss background on graphlet kernels and smoothing techniques. In Section 3, we introduce our Kneser-Ney-inspired smoothing technique. In Section 4, we propose an alternate Bayesian version of our model. Related work is discussed in Section 5. In Section 6, we perform experiments and discuss our findings, and we conclude the paper with Section 7.

Notation
Graphlets are small, connected, non-isomorphic sub-graphs of a large network. They were introduced by Przulj [28] to design a new measure of local structural similarity between biological networks. Graphlets up to size five are shown in Figure 3.

The Graphlet kernel
Let G k = {g1, g2, . . . , gn k } be the set of size-k graphlets where n k denotes the number of unique graphlets of size k. Given a graph G, we define fG as a normalized vector of length n k whose i-th component corresponds to the frequency of occurrence of gi in G: Here ci denotes number of times gi occurs as a sub-graph of G. Given two graphs G and G , the graphlet kernel kg is defined as: which is simply the dot product between the normalized graphletfrequency vectors.

Smoothing multinomial distributions
In this section we will briefly review smoothing techniques for multinomial distributions and show that graphlet kernels are indeed based on estimating a multinomial. Suppose we observe a sequence e1, e2, . . . , en containing n discrete events drawn from a ground set of size M and we would like to estimate the probability P (ei) of observing each event ei. Maximum likelihood estimation based on sequence counts obtained from the observations provides a way to compute P (ei): where ci denotes the number of times the event ei appears in the observed sequence and j cj denotes the total number of observed events. Therefore, one can easily see that the representation used in graphlet kernels in Section 2.2 is actually an MLE estimate on the observed sequences of graphlets. Figure 3: Connected, non-isomorphic induced sub-graphs of size k ≤ 5. Plots are generated with NetworkX library [14].
However, MLE estimates of the multinomial distribution are spiky, that is, they assign zero probability to events that did not occur in the observed sequence. What this means is that an event with low probability is often estimated to have zero probability mass. This issue occurs in a number of different domains and therefore, unsurprisingly, has received significant research attention [42]; smoothing methods are typically used to address this problem. The general idea behind smoothing is to discount the probabilities of the observed events and to assign extra probability mass to unobserved events.
Laplace smoothing is the simplest and one of the oldest smoothing methods, where only a fixed count of 1 is added to every event. This results in the estimate or equivalently, where λ is a normalization factor which ensures that the distributions sum to one. The intuition behind Laplace smoothing is basically to interpolate a uniform distribution with the MLE distribution. Although Laplace smoothing resolves the zero-count problem, it does not produce a power law distribution which is a desirable feature in real-life models. Therefore, researchers have worked on finding smoothing techniques that respect power law behavior. The key idea behind these methods is to redistribute the probability mass using a so-called fallback model, where the fallback model is also recursively estimated. Kneser-Ney smoothing is a fallback based smoothing method which has been identified as the state-of-the-art smoothing in natural language processing by several studies [7]. Kneser-Ney smoothing computes the probability of an event by using the raw counts that are discounted by using a fixed mass. Then, the discounted mass is re-added equally to all event probabilities by using a base distribution: where d ≥ 0 is the discounting parameter, P0(·) is the base distribution and P0(ei) denotes the probability mass the base distribution assigned to event ei. The quantity n j=1 |{ej : cj > d}| is a normalization factor to ensure that the distribution sums to 1 and simply denotes the number of events the discount is applied. When discount parameter d = 0, we recover MLE estimation since no mass is taken away from any event. When d is very large, then we recover the base distribution on the events since we discount all the available mass. One should interpolate between these two extremes in order to get a reasonable smoothed estimate. In order to propose a new Kneser-Ney-based smoothing framework, one needs to specify the discount parameter d and the base distribution P0(·). In the next section we will show how one can derive a meaningful base distribution for graphlets.

DEFINING A BASE DISTRIBUTION
The space of graphlets has an inherent structure. One can construct a directed acyclic graph (DAG) in order to show how different graphlets are related with each other. A node at depth k denotes a graphlet of size k. When it is clear from context, we will use the node of the DAG and a graphlet interchangeably. Given a graphlet gi of size k and another graphlet gj of size k + 1 we add an edge from gi to gj if, and only if, gi can be obtained from gj by deleting a node of gj. We first discuss how to construct this DAG and then discuss how we can use this DAG to define a base distribution.
The first step towards constructing our DAG is to obtain all unique graphlet types of size k. Therefore, we first exhaustively generate all possible graphs of size k (this involves a one time O(2 k ) effort), and use Nauty [25] to obtain their canonically-labelled isomorphic representations. In order to obtain the edges of the DAG, we take a node at depth k + 1, which denotes a canonically-labeled isomorphic graphlet and delete a node to obtain a size k graphlet. We use Nauty to canonically-label the size k graph, which in turn allows us to link to a node at depth k. By deleting each node of the k+1 sized graphlet we can therefore obtain k + 1 possible links. We repeat this for all nodes at level k + 1 before proceeding to level k. Figure 4 shows the constructed DAG for size k = 5 graphlets. Since all descendants of a given graphlet at level k are at level k + 1, a topological ordering of the vertices is possible, and hence it is easy to see that the resulting graph is a DAG. Now, let us define the edge weight between an arbitrary graphlet gj of size k + 1 and its parent gi of size k. Let sij denote the number of times gi occurs as a sub-graph of gj and Cg i denote all the children of graphlet gi in the DAG. Then, we define the edge Figure 4: Topologically sorted graphlet DAG where each node are within one delete-distance away. Nodes are colored by average degree. Image is generated by Gephi [1]. Figure 5: Graphlet g15 get the probability mass to its parents g7, g6, g5 according to the weights w1, w2, w3 respectively.
weight between graphlet gi and gj as Next we show how the DAG can be used to define a base distribution. Suppose we have a distribution over graphlets of size k. Then we can transform it into a distribution over size k+1 graphlets in a recursive way by exploiting edge connections in the DAG as follows: Here gj denotes a graphlet of size k + 1 and P a(gj) denotes the parents of graphlet gj in the DAG. LEMMA 1. Given the set of child nodes of a graphlet gi, if the edge weights on the DAG are all non-negative and satisfy then (8) defines a valid probability distribution.
PROOF. For clarity, we introduce a few notations to facilitate the proof. Assume that there are in total J child graphlets, which we denote by g1, g2, · · · , gJ . Further assume that there are in total I parent graphlets, which we denote by g1, · · · , gI . Let Ij denote the number of parents graphlet gj has, i.e. Ij = |P a(gj)|. Clearly, we have that J j=1 Ij = I. Thus, we have which completes the proof.
The base distribution we defined above respects the structural space of the graphlets. Pretend for a moment that we are given only the frequencies of occurrences of size k graphlets and are asked to infer the probability of occurrences of size k + 1 graphlets. Without any additional information, one can infer the distribution as follows: each parent graphlet casts a vote for every child graphlet based on how many times the parent occurs in the child. The votes of all parents are accumulated and this provides a distribution over size k + 1 graphlets. In other words, a natural way to infer the distribution at level k + 1 is to use how likely we are to see its subgraphs. Figure 5 illustrates the relationship between a graphlet of size k + 1 and its parent graphlets of size k. Here, edge weights denote how many times each parent occurs as a sub-graph of g15. In the case that we do not observe graphlet g15, it still gets probability mass proportional to the edge weight from its parents g7, g6, g5, thus overcoming the sparsity problem of unseen data. Our model combines this base distribution with the observed real data to generates the final distribution.
The way how the discounted mass is distributed is controlled by the edge weights between two graphlets. In Equation 7 we defined edge weights according to the number of times a parent node occurs in its children. However, one can explore different weighting schemes between the nodes on the DAG based on domain knowl-edge. For example, in the case of structured graphs such as social networks, one might benefit from weighting the edges according to the PageRank [26] score of the nodes. Similarly, other link analysis algorithms such as Hubs or Authority given by HITS algorithm [19] can be used in order to exploit the domain knowledge.

Kneser-Ney Smoothing with a structural distribution
We now have all the components needed to explain our Structural Kneser-Ney (SKN) framework. Given an arbitrary graphlet gj of size k + 1, we estimate the probability of observing that graphlet as follows: As can be seen from the equation, we first discount the count of all graphlets by d, and then redistribute this mass to all other graphlets. The amount of mass a graphlet receives is controlled by the base distribution. In order to automatically tune the discount parameter d, we use the Pitman-Yor process (a Bayesian approximation of Kneser-Ney) in the next section.

PITMAN-YOR PROCESS
We will only give a very high level overview of a Pitman-Yor process and refer the reader to the excellent papers by Teh [35] and [12] for more details. A Pitman-Yor process P on a ground set G k+1 of size-(k + 1) graphlets is defined via where d k+1 is a discount parameter 0 ≤ d k+1 < 1, θ > −d k+1 is a strength parameter, and P k is a base distribution. The most intuitive way to understand draws from the Pitman-Yor process is via the Chinese restaurant process (also see Figure 6). Consider a restaurant with an infinite number of tables. Customers enter the restaurant one by one. The first customer sits at the first Intuitively the reason this process generates power-law behavior is because popular graphlets which are served on tables with a large number of customers have a higher probability of attracting new customers and hence being generated again. This self reinforcing property produces power law behavior. In a hierarchical Pitman-Yor process, the base distribution P k is recursively defined via a Pitman-Yor process P k ∼ P Y (d k , θ k , P k−1 ). In order to label a table, we need a draw from P k , which is obtained by inserting a customer into the corresponding restaurant.
In our case P k+1 is defined over G k+1 of size n k+1 while P k is defined over G k of size n k ≤ n k+1 . Therefore, like we did in the case of Kneser-Ney smoothing we will use the DAG and Equation (8) to define a base distribution. This changes the Chinese Restaurant process as follows: When we need to label a table, we will first draw a size-k graphlet gi ∼ P k by inserting a customer into the corresponding restaurant. Given gi, we will draw a size-(k + 1) graphlet gj proportional to wij, where wij is obtained from the DAG. Deletion of a customer is handled similarly. Detailed pseudo-code can be found in Algorithms 1 and 2.

Algorithm 1 Insert a Customer
Input: d k+1 , θ k+1 , P k t ← 0 // Occupied tables c ← () // Counts of customers l ← () // Labels of tables if t = 0 then t ← 1 append 1 to c draw graphlet gi ∼ P k // Insert customer in parent draw gj ∼ wij append gj to l return gj else with probability ∝ max(0, cj − d) cj ← cj + 1 return lj with probability proportional to θ + dt t ← t + 1 append 1 to c draw graphlet gi ∼ P k // Insert customer in parent draw gj ∼ wij append gj to l return gj end if Algorithm 2 Delete a Customer Input: d, θ, P0, C, L, t with probability ∝ c l c l ← c l − 1 gj ← lj if c l = 0 then P k ∝ 1/wij delete c l from c delete lj from l t ← t − 1 end if return g

RELATED WORK
The problem of estimating multinomial distributions is a classic problem. In natural language processing they occur in the following context: suppose we are given a sequence of words w1, . . . , w k and one is interested in asking what is the probability of observing word w next. Estimating this probability lies at the heart of language models, and many sophisticated smoothing techniques have been proposed. This is a classic multinomial estimation problem that suffers from sparsity since the event space is unbounded. Moreover, natural language exhibits power law behavior since the distribution tends to be dominated by a small number of frequently occurring words. In extensive empirical evaluation it has been found the Kneser-Ney smoothing is very effective for language Figure 6: An illustration of table assignment, adapted from Goldwater et al. [13]. In this example, labels at the tables are given by (l1, . . . , l4) = (G44, G30, G32, G44). Black dots indicate the number of occurrences of each label in 10 draws from the Pitman-Yor process. models [7], [24]. Here, the base distribution is constructed using smaller context of k−1 words which naturally leads to a denser distribution. Even though language models and graphlets have some similarities, there is a significant fundamental difference between the two. In language models, one can derive the base distribution using a smaller context. However, in the case of graphlets there is no equivalent concept of a fallback model. Therefore, we need to derive the base distribution by using smaller size graphlets. However, this leads to a problem since the distribution is now defined on a smaller space. Therefore, we need to apply a transformation by using the DAG in order to convert the distribution back into to the original space. Goldwater et al. [13] and Teh [35] independently showed that Kneser-Ney can be explained in a Bayesian setting by using the Pitman-Yor Process (PYP) [27]. In the Bayesian interpretation, a hierarchical PYP where the Pitman-Yor prior comes from another PYP is used. Similar to Kneser-Ney, this interpretation is not directly applicable to our model since the previous PYP has a different space, thus we need to apply a transformation.
Graph kernels can be considered as special cases of convolutional kernels proposed by Haussler [15]. In general, graph kernels can be categorized into three classes: graph kernels based on walks and paths [11], [18], [5], graph kernels based on limited-size subgraphs [16], [31], [34] and graph kernels based on subtree patterns [29]. Shervashidze et al. [33] performs a relaxation on the vertices and exploit labeling information embedded in the graphs to derive their so-called Weisfeiler-Lehman kernels. However, their kernel is applicable only to labeled graphs. The sparsity problem of graphlet kernels has been addressed before. Hash kernels proposed by Shi et al. [34] addresses the sparsity problem by applying a sparse projection into a lower dimensional space. The idea here is that many higher order graphlets will "collide" and therefore be mapped to the same lower dimensional representation, thus avoiding the diagonal dominance problem. Unfortunately, we find that in our experiments the hash kernel is very sensitive to the hash value used for embedding and rarely performed well as compared to the MLE estimate.

EXPERIMENTS
To compare the efficacy of our approach, we compare our Kneser-Ney and Pitman-Yor smoothed kernels with state-of-the-art graph kernels namel the graphlet kernel [32], the hash kernel [34], the random walk kernel [11], [18], [38], and the shortest path kernel [5]. For random walk kernel, we uniformly set the decay factor λ = 10 −4 , for shortest path we used the delta kernel to compare the shortest-path distances, and for the hash kernel we used a prime number of 11291. We adopted Markov chain Monte Carlo sampling based inference scheme for the hierarchical Pitman-Yor language model from [35] and modified the open source implementation of HPYP from https://github.com/redpony/cpyp. Due to lack of space we will only present a subset of our experimental results. Full results including the source code and experimental scripts will be made available for download from http: //cs.purdue.edu/~ypinar/kdd.

Datasets
In order to test the efficacy of our model, we applied smoothing to real-world benchmark datasets, namely MUTAG, PTC, NCI1, NCI109, ENZYMES and DD. MUTAG [9] is a binary data set of 188 mutagenic aromatic and heteroaromatic nitro compounds, labeled whether they have mutagenicity in Salmonella typhimurium. The Predictive Toxicology Challenge (PTC) [36] dataset is a chemical compound dataset that reports the carcinogenicity for male and female rats. NCI1 and NCI109 [40], (http://pubchem.ncbi.nlm.nih.gov) datasets, made publicly available by the National Cancer Institute (NCI), are two subsets of balanced data sets of chemical compounds screened for ability to suppress or inhibit the growth of a panel of human tumor cell lines. Enzymes is a data set of protein tertiary structures obtained from [4]. DD [10] is a data set of protein structures where each protein is represented by a graph and nodes are amino acids that are connected by an edge if they are less than 6 Angstroms apart. Table 2 shows summary statistics for these datasets. Note that we did not use edge or node labels in our experiments.

Experimental Setting
All data sets we work with consist of sparse graphs. However, counting all graphlets of size k for a graph with n nodes requires O(n k ) effort which is intractable even for moderate values of k. Therefore, we use random sampling, as advocated by [32], in order to obtain an empirical distribution of graphlet counts that is close to the actual distribution of graphlets in the graph. For each value of for each k ∈ {2, . . . , 8} we randomly sampled 10,000 sub-graphs, and used Nauty [25] to get canonically-labeled isomorphic representations which are then used to construct the frequency representation.  We performed 5-fold cross-validation with C-Support Vector Machine Classification using LibSVM [6], using 4 folds for training and 1 for testing. We used a linear kernel. Since we are interested in understanding the difference in performance between smoothed and un-smoothed kernels we did not tune the value of C; it was simply set to 1. In order to tune the discount parameter for Kneser-Ney based smoothed kernel, we tried different parameters vary from 0.01 to 10,000 and report results for the best one. Effect of discounting parameter on performance First, we investigate the effect of the discounting parameter on the classification performance. Since the trends are similar across different datasets, we pick PTC as a representative dataset to report results. Figure 8 shows the classification accuracy on the PTC dataset with different discounting parameters for Kneser-Ney smoothing. As expected, applying very large discounts decreases the performance because the distribution (10) degrades to the base distribution. On the other hand, applying a very small discount also decreases the accuracy since the distribution degrades to the MLE estimate. From our experiments, we observe that the best performance in all datasets is achieved by using an intermediate discount value between these two extremes. However, the specific discount value is data dependent. Effect of graphlet size on performance Next we investigate how the value of k affects performance. Again, we show results for a representative dataset namely PTC. Figure 7 shows the classification accuracy on the PTC dataset as a function of graphlet size k for MLE (the graphlet kernel), Pitman-Yor smoothed kernel and Kneser-Ney smoothed kernel. From the figure, we can see that small graphlet sizes such as k = 2, 3 do not perform well and are not informative since the number of unique graphlets is very small (see Table 1). On the other hand, MLE does not perform well for large graphlet sizes such as k = 7, 8 because of the diagonal dominance problem. On the other hand, smoothed kernels in general obtain a balance between these two extreme situations and tend to yield better performance. Pitman-Yor smoothed kernel tends to achieve a better performance than MLE, but doesn't perform as good as Kneser-Ney. This is expected since we didn't tune the hyperparameters for Pitman-Yor process. Teh [35] shows that Pitman-Yor yields a better performance if one tune the hyperparameters. Therefore, our Pitman-Yor kernel is open to improvement.

Comparison with related work
We compare the proposed smoothed kernel with graphlet kernel and hash kernel on the benchmark data sets in Table 1. The results for the Shortest Path and Random Walk graph kernels are included mainly to show what is the state-of-the-art using other representations. We fixed the k = 5 which is observed to be the best value for the MLE based graphlet kernel on most datasets. We randomly sample 10,000 graphlets from each graph and feed the same frequency vectors to graphlet kernel (GK), hash kernel (HK), Kneser-Ney smoothed kernel (KN), and Pitman-Yor smoothed kernel (PYP). Therefore, the differences in performance that we observe are solely due to the transformation of the frequency vectors that these kernels perform. We performed an unpaired t-test and use bold numbers to indicate that the results were statistically significant at p < 0.0001.
We can see that KN kernel outperforms MLE and hash kernels on all of the benchmark data sets. The accuracy of the PYP kernel is usually lower than that of the KN kernel. We conjecture that this is because the Pitman-Yor process is sensitive to the hyperparameters and we do not carefully tune the hyper-parameters in our experiments. On PTC, DD and Enzymes datasets, the KN kernel reached the highest accuracy. On MUTAG, NCI1 and NCI109 datasets, KN kernels also yield good results and got comparable classification accuracies to shortest path and random walk kernels. For the DD dataset, shortest path and random walk kernels were not able to finish in 24 hours, due to the fact that this dataset has a large maximum degree.
To summarize, smoothed kernels turns out to be competitive in terms of classification accuracy on all datasets and are also applicable to very large graphs.

Effect of exhaustive sampling on performance
Next, we investigate whether any of the difference in performance can be attributed to sampling a small number of graphlets. In other words, we ask do the results summarily change if we performed exhaustive sampling instead of using 10,000 samples. We give MLE an unfair advantage by performing exhaustive sampling on MUTAG, PTC, NCI and NCI109 datasets for k = 5 by using a distributed memory implementation. Table 4 shows mean, median and standard deviations of number of samples in bruteforce sampled datasets for k = 5. Here, we can see that the original frequencies of the graphlets are quite high in most of the datasets. Even though our algorithm only uses 10,000 samples, it outperforms the graphlet kernel with exhaustive sampling on MUTAG, PTC and NCI1 and achieves competitive performance on NCI109 dataset. Results are summarized in Table 3.
Even though distribution with a small number of samples is close to the original distribution in the L1 sense, bruteforce sampling reveals that the true underlying distribution of the datasets contains a larger number of unique graphlets comparing to random sampling. Since the graphlet kernel uses a MLE estimate its performance degrades. On the other hand, our smoothing technique uses structural

DISCUSSION
We presented a novel framework for smoothing normalized graphletfrequency vectors inspired by smoothing techniques from natural language processing. Although our models are inspired by work done in language models, they are fundamentally different in the way they define a fallback base distribution. We believe that our framework has applicability beyond graph kernels, and can be used in any structural setting where one can naturally define a relationship such as the DAG that we defined in Figure 4. We are currently investigating the applicability of our framework to string kernels [22,23,37] and tree kernels [8]. It is also interesting to investigate if our method can be extended to other graph kernels such as random walk kernels. Our framework is also applicable to nodelabeled graphs since they also suffer from similar sparsity issues. We leave the application of our framework to labelled graphs to an extended version of this paper. We are also investigating better strategies for tuning the hyper-parameters of the Pitman-Yor kernels.