Moment based estimation of stochastic Kronecker graph parameters

Introduction

Stochastic Kronecker graphs were introduced by as a method for simulating very large random graphs. Random synthetic graphs are used to test graph algorithms and to understand observed properties of graphs. By using simulated graphs, instead of real measured ones, it is possible to test algorithms on graphs larger or denser than presently observed ones. Simulated graphs also allow one to judge which features of a real graph are likely to hold in other similar graphs and which are idiosyncratic to the given data.

Stochastic Kronecker graphs are able to serve these purposes through a model that has only three or four parameters. Parameter estimation poses unique challenges for those graphs. The main problem is that for a graph with $N$ nodes, the likelihood has contributions from $N!$ permutations of the nodes . In practice, many thousands or millions of randomly sampled permutations are used to estimate the likelihood. Even then it takes more than $O(N^{2})$ work to evaluate the likelihood contribution from one of the permutations.

In this paper we present a method of moments strategy for parameter estimation. While moment methods can be inefficient compared to maximum likelihood, statistical efficiency is of reduced importance for enormous samples and in settings where the dominant error is lack of fit. The method equates expected to observed counts for edges, triangles, hairpins ($2$-stars or wedges) and tripins ($3$-stars). The Kronecker model gives quite tractable formulas for these moments.

The outline of this paper is as follows. Section 2 defines Kronecker graphs and introduces some notation. Section 3 derives the expected feature counts. Section 4 describes how to solve method of moment equations for the parameters of the Kronecker graph model. Section 5 presents some examples on fitting Kronecker models to some real world graphs. We compare several moment based ways to estimate Kronecker graph parameters and find the most reliable results come from a criterion that sums squared relative errors between observed and expected features. We find that the fitted Kronecker models usually underestimate the number of triangles compared to the real graphs. While our parameter estimates underestimate triangle counts and some other features, we find that they provide much closer matches than some previously published parameters fit by KronFit. Section 6 fits parameters to graphs that were randomly generated from the Kronecker model. We find that the estimated parameters closely track their generating values, with some small bias when a parameter is at the extreme range of valid values. Section 7 has our conclusions.

along with the code used to estimate Kronecker parameters.

The Kronecker model

Given a node set $\mathcal{N}$ of cardinality $N\geq 1$, and a matrix $P_{ij}\in$ defined over $i,j\in\mathcal{N}$, a random graph $G^{*}(P)$ is one where the edge $[ij]$ exists with probability $P_{ij}$ and all $N^{2}$ edges exist or don’t independently. The graph $G^{*}$ includes loops $[ii]$ and may possibly include both $[ij]$ and $[ji]$. We snip these out by defining the random graph $G(P)$ with edges $[ij]$ only when $i\neq j$ and $[\max(i,j),\min(i,j)]\in G^{*}$, using any non-random ordering of $\mathcal{N}$. Both $G^{*}$ and $G$ are in fact probability weighted ensembles of graphs, but for simplicity we describe them as single random graphs. We assume that $P$ is a symmetric matrix and so the ordering of nodes does not affect the distribution.

An extremely parsimonious stochastic Kronecker graph takes $P$ to be the $r$-fold Kronecker product of $\Theta=\begin{pmatrix}a&b\\ b&c\end{pmatrix},$ for $a,b,c\in$. That is

If the power $r$ is known, then only three numbers need to be specified, and with them we can then simulate other graphs that are like the original. Perhaps surprisingly, stochastic Kronecker graphs imitate many, but of course not all, of the important features seen in large real world graphs. See for example .

We would like to pick parameters $a,b,c\in$ to match the properties seen in a real and large graph. Parameter matrices $\Theta=\begin{pmatrix}a&b\\ b&c\end{pmatrix}$ and $\Theta^{*}=\begin{pmatrix}c&b\\ b&a\end{pmatrix}$ give rise to the same graph distribution. To force identifiability, we may assume that $a\geq c$.

Moment formulas

The Kronecker structure in $P$ makes certain aspects of $G$ very tractable. For example, the number $E$ of edges in $G$ can be shown to have expectation

This section derives equation (1) and similar formulas for the expected number of features of various types. The expected feature counts require sums over various sets of nodes. Section 3.1 records some summation formulas that simplify that task. Then Section 3.2 turns expected feature counts into sums and Section 3.3 shows how those sums simplify for stochastic Kronecker matrices.

Equation (4) is more complicated than the others. It can be proved by defining $g_{ijk}=\sum_{l}f_{ijkl}-f_{ijki}-f_{ijkj}-f_{ijkk}$, writing $\sum\!^{{}^{*}}_{ijkl}\,f_{ijkl}=\sum\!^{{}^{*}}_{ijk}g_{ijk}$ and then applying (3).

In some of our formulas below, the first index is singled out but the others are exchangeable. By this we mean that $f_{ijk}=f_{ikj}$, when there are three indices, while $f_{ijkl}=f_{ijlk}=f_{ikjl}=f_{iklj}=f_{iljk}=f_{ilkj}$ is the version for four indices.

When indices after the first are exchangeable, then equation (3) simplifies to

When all indices $ijk$ are exchangeable, so that $f_{ijk}=f_{ikj}=f_{jik}=f_{jki}=f_{kij}=f_{kji}$, then equation (5) simplifies to

2 Expected feature counts for independent edges

The graph features we describe are shown in Figure 1. In addition to edges, there are hairpins ($2$-stars) where two edges share a common node, tripins ($3$-stars) where three edges share a node, and triangles. The Kronecker model has independent edges. Here we find the expected feature counts for any random graph where edge $[ij]$ appears with probability $P_{ij}$ and edges are independent.

Recall that $G^{*}$ is a random graph with $\Pr([ij]\in G^{*})=P_{ij}$ (independently). Let it have incidence matrix $A^{*}$. There may be loops $A_{ii}^{*}=1$, and for $i\neq j$, $A_{ij}^{*}$ and $A_{ji}^{*}$ are independently generated. The graph $G$ is formed by deleting loops from $G^{*}$ and symmetrizing the incidence matrix via

The number of edges in $G$ is $E=(1/2)\sum\!^{{}^{*}}_{ij}A_{ij}$. The expected number of edges satisfies

The number of hairpins in $G$ is $H=(1/2)\sum\!^{{}^{*}}_{ijk}A_{ij}A_{ik}$. Dividing by two adjusts the sum for counting $\{[ij],[ik]\}$ twice. The expected value of $H$ satisfies

by letting $f_{ijk}=P_{ij}P_{ik}$, for which $f_{ijk}=f_{ikj}$, and applying equation (5).

The number of triangles in $G$ is $\Delta=(1/6)\sum\!^{{}^{*}}_{ijk}A_{ij}A_{ik}A_{jk}$, because the sum counts each triangle $3!=6$ times. The expected value of each term is $f_{ijk}=P_{ij}P_{ik}P_{jk}$ which is symmetric in its three arguments and so we may apply equation (7) to get

The number of tripins in $G$ is $T=(1/6)\sum\!^{{}^{*}}_{ijkl}A_{ij}A_{ik}A_{il}$. The final three indices in $f_{ijkl}=P_{ij}P_{ik}P_{il}$ are exchangeable, and so equation (6) applies. Thus

3 Simplifying the sums

The sums in the expected counts simplify, because of the properties of the Kronecker graph. Let the node set be $\mathcal{N}=\mathcal{N}_{r}=\{0,1,\dots,2^{r}-1\}$. For $i\in\mathcal{N}$ write $i=\sum_{s=1}^{r}2^{s-1}i_{s}$ for $i_{s}\in\{0,1\}$. Similarly let $j$, $k$, and $l$ be described in terms of $j_{s},k_{s},l_{s}\in\{0,1\}$ for $s=1,\dots,r$.

The matrix entry $P_{ij}=P^{(r)}_{ij}$ may be written

For $r\geq 2$, we simplify the expression by induction using a smaller version of the problem defined via $P^{(r-1)}$. Specifically,

where indices $i_{s}$, $j_{s}$ and $k_{s}$ are summed over their full ranges, and the indices $i$, $j$, $k$ for $P^{(r-1)}_{ij}P^{(r-1)}_{ik}$ are summed over the node set $\mathcal{N}_{r-1}=\{0,\dots,2^{r-1}-1\}$.

All of the sums of products of elements of $P^{(r)}_{ij}$ listed in the previous section, with summation over all levels of each index, also reduce this way to $r$’th powers of their value for the case $r=1$.

The first expression (9) follows by summing over the $8$ rows of Table 1. As a result

In the rest of this section, we record the other sums we need. First, the sums over one index variable take the form

where cases $m=1,2,3$ are used in our expected feature counts. The sums over two index variables are

The cases we need are for $(m,n)\in\{(0,1),(0,2),(0,3),(1,1),(1,2)\}$.

Four sums over three indices are used. They are:

Finally, one sum over four indices is used:

4 Expected feature counts

Now we can specialize the results of Section 3.2 to the Kronecker graph setting. Gathering together the previous developments, we find

The first term will dominate for large $r$ unless $b\ll a+c$. The relative magnitude of the second term is

5 Illustrations

If $a=b=c=1$, then $G^{*}$ has every possible edge and loop with probability $1$. As a result $G$ is the complete graph on $N=2^{r}$ nodes. Then it has $N(N-1)/2$ edges, $N(N-1)(N-2)/2$ hairpins, $N(N-1)(N-2)/2$ triangles, and it has $N(N-1)(N-2)(N-3)/6$ tripins.

Solving for aa, bb, and cc

There are four equations in Section 3.4. To estimate $a$, $b$, and $c$ will require at least three of them. Because they are high order polynomials it is possible that there are multiple solutions or even none at all. The latter circumstance would provide some evidence of lack of fit of the stochastic Kronecker model to a given graph. Regardless, each of the equations involves the count of a feature in the graph.

Three of the features we use are easily obtainable from the degrees of the nodes. Let $d_{i}=\sum_{j\in\mathcal{N}}A_{ij}$ be the degree of node $i$ in graph $G$. Then

give the number of edges, hairpins (or wedges), and tripins in terms of the degrees $d_{i}$.

The number of triangles $\Delta$ is not a simple function of $d_{i}$. Algorithms to count triangles are considered in . The time complexity can be as low as $O(E^{3/2})$, and sometimes even lower for approximate counting .

2 Objective functions

A pragmatic way to choose $a$, $b$, and $c$ is to solve

where the sum is over three or four of the features $F\in\{E,H,T,\Delta\}$ from Section 3.4 and the minimization is taken over $0\leq c\leq a\leq 1$ and $0\leq b\leq 1$. The terms in (11) are scaled by an approximate variance. A sharper expression would account for correlations among the features used. That should increase statistical efficiency, but in large problems lack of fit to the Kronecker model is likely to be more important than inefficiency of estimates within it.

Many real world networks may not have good fits in terms of these three Kronecker parameters. This is the case for most of the forthcoming experiments. The following more general objective can be more robust to these deviances:

Here $D$ is either of the two distance functions:

We will find in Section 5 below that robust results arise from the combination $D_{\text{sq}}$ and $N_{F^{2}}$, for which (12) reduces to

Because there are only three parameters, the criterion (12) can simply be evaluated over a grid inside $\{(a,b,c)\in^{3}\mid a\geq c\}$. To be sure of taking a point within $\varepsilon$ of the minimizer takes work $O(\varepsilon^{-3})$. An alternative is to employ a general nonlinear minimization procedure. The remainder of this section looks at a method to reduce that effort.

3 Matching leading terms

In a synthetic graph $N=2^{r}$ is known. When fitting to a real world graph a pragmatic choice is $r=\lceil\log_{2}(N)\rceil$. The interpretation is that the random graph $G^{*}$ may have had isolated nodes that were then dropped when forming $G$, but we suppose that fewer than half of the nodes in $G^{*}$ have been dropped.

If we consider just the lead terms, then we could get estimates $\hat{a}$, $\hat{b}$, and $\hat{c}$ by solving three of the equations:

The equations for $e$ and $h$ together can be solved to get

where we have assumed that $a\geq c$. The transformed tripin count $t$ matches $\hat{x}^{3}+\hat{y}^{3}$ and so it is redundant given $e$ and $h$, if we are just using lead terms. We must either count triangles, or use higher order terms.

Equation (14) may fail to have a meaningful solution. At a minimum we require $e^{2}\leq 2h$ and $e\geq\sqrt{2h-e^{2}}$. These translate into

The left hand inequality in (15) holds for any graph, but the right hand side need not. It holds when $N^{-1}\sum_{i}(d_{i}-\bar{d})^{2}\geq\bar{d}=N^{-1}\sum_{i}d_{i}$. If the variance of the node degrees $d_{i}$ is smaller than their mean, then equation (14) does not have real valued solutions. The degree distribution of a stochastic Kronecker graph has heavy tails . Therefore in applications where that model is suitable equation (14) will give a reasonable solution.

When $d_{i}$ have a variance larger than their mean, then we can do a univariate grid search for $b\in$ using equation (14) to get $a=x-b\equiv a(b)$ and $c=y-b\equiv c(b)$. The choice of $b$ can then be made as the minimizer of $|a(b)^{3}+c(b)^{3}+3b^{2}(a(b)+c(b))-\delta|$.

Examples

In this section, we experiment with different techniques for fitting the parameters of the Kronecker model. These experiments involve 8 real world networks whose statistical properties are listed in the rows of the forthcoming tables labeled “Source.”

The networks ca-GrQc, ca-HepTh, ca-HepPh are co-authorship networks from arXiv . The nodes of the network represent authors, and there is an edge between two nodes when the authors jointly wrote a paper. Likewise, the hollywood-2009 network is a collaboration graph between actors and actresses in IMDB . Nodes are actresses or actors, and edges are collaborations on a movie, as evidenced by jointly appearing on the cast. These networks are naturally undirected and all edges are unweighted.

Both as20000102 and as-Skitter are technological infrastructure networks . Each node represents a router on the internet and edges represent a physical or virtual connection between the routers. Again, these networks are undirected and unweighted.

The wikipedia-20051105 graph is a symmetrized link graph of the articles on Wikipedia generated from a data download on November 5th, 2005 . The underlying network is directed, but in these experiments, we have converted it into an undirected network by dropping the direction of the edges.

All of the previously described networks have distinctly skewed degree distributions. That is, there are a few authors, actors, routers, or articles with a large number of links, despite the overall network having a small average degree. The final network we study is usroads, a highway-level network from the National Highway Planning Network (http://www.fhwa.dot.gov/planning/nhpn/), which does not have a highly skewed distribution. We include it as an example of a nearly planar network. It is also naturally undirected.

In two of the experiments, we generate synthetic Kronecker networks. The algorithm to realize these networks is an explicit coin-flipping procedure instead of the more common ball-dropping method . For each cell $i,j$ in the $2^{r}-1\choose 2$ upper triangular portion, we first determine the log of the probability of a non-zero value in that cell, then generate a random coin flip with that probability as heads and record an edge when the coin comes up heads. This procedure is scalable because the full matrix of probabilities is never formed. It is also easily parallelizable. Our implementation uses pthreads to exploit multi-core parallelism. It takes somewhat more work than the ball-dropping procedure, scaling as $O(r2^{2r})$ instead of $O(rm)$, where $m$ is the number of balls dropped. Often $m\approx 2^{r+3}$, that is, $8$ balls per vertex . Each ball generates about one edge; see for a more thorough analysis. Coin-flipping preserves the exact Kronecker distribution whereas ball-dropping is an approximation.

The experiments with these networks investigate (i) the difference in results from the various choices of $D$ and $N$ in the objective (12); (ii) the fitted parameters to the 8 real world networks; and (iii) the difference in fitted parameters when only using three of the four graph features.

The first study regards the choice of objective function. Of eight possible combinations of distance and normalization, we considered two to be unreasonable a priori. Here we investigate the other six pairs.

Table 2 shows the different parameters $a,b,$ and $c$ chosen by each objective function, as well as the expected feature counts for those parameters for three graphs: a single realization of a Kronecker graph with $a=0.99,b=0.48,c=0.25$, the collaboration network ca-GrQc, and the infractucture network as20000102. The rows labeled “Source” contain the actual feature counts in each network. The optimization algorithm to pick $a,b,c$ uses the best objective value from three procedures. First, it tries 50 random starting points for the fmincon function in Matlab R2010b, an active set algorithm. Then, it performs a grid search procedure with 100 equally spaced points in each dimension. Finally, it tries the leading term matching algorithm from Section 4.3, and considers those parameters.

Based on these results, either of the objectives $D_{\text{sq}},N_{F^{2}}$ or $D_{\text{abs}},N_{F}$ appears to be a robust choice when the model does not fit exactly. Due to the continuity of the $D_{\text{sq}}$ function, the rest of our fits in this manuscript uses the $D_{\text{sq}},N_{F^{2}}$ variation.

2 Parameters for real-world networks

For the 8 networks previously described, we use the objective function (12) with $D_{\text{sq}},N_{F^{2}}$ to fit the parameters $a,b,c$. The results, along with the expected feature counts for the fitted parameters, are presented in Table 3. We show the minimizer for the three different strategies to optimize the objective described in the previous section: a direct minimization procedure, the grid search procedure, and the leading term matching approach (Section 4.3). For each approach, the table also shows the time required for that algorithm and the value of the objective function at the minimizer.

Leskovec et al. provide the fitted parameters $a,b,$ and $c$ from their KronFit algorithm for the networks ca-GrQc, ca-HepTh, ca-HepPh, and as20000102. We include them in Table 3 for comparison. In all cases but one, the expected feature count using KronFit is farther from the observed feature count than the expectation under our moment based fits. Sometimes it is much farther. There was one exception. For the graph as20000102, KronFit gave a better estimate of the number of edges than our moment method gave.

KronFit typically underestimates the feature counts. The effect is severe for triangles. Kronecker random graphs commonly have many fewer triangles than the real world graphs to which they are fit. Our moment based estimators find parameters leading to many more triangles than the KronFit parameters do.

In fairness, we point out that our method is designed to match expected to observed feature counts, while KronFit fits by maximum likelihood. Therefore the evaluation criterion is closer to the fitting criterion for us. But maximum likelihood ordinarily beats or matches the method of moments in large samples from parametric models; it’s mismatching criteria are more than compensated for by superior statistical efficiency. The explanation here may involve maximum likelihood being less robust to lack of fit of the Kronecker model, or it may be that KronFit is not finding the MLE.

The results in Table 3 show small differences in the fits between the direct and grid algorithms, although the direct algorithm is much faster. The leading term matching algorithm, when it succeeds, generates similar Kronecker parameters, although with a distinctly worse objective value. The results from the KronFit algorithm differ and likely match the graph in another aspect.

Lead term matching is tens of times faster than direct search and roughly $1000$ times faster than grid search. But even the grid search takes under a minute in our examples, so the speed savings from the lead term approach is of little benefit here. For the large graphs, the time to compute the network features dominates the time to fit the parameters, showing that this approach scales to large networks.

Overall, the results indicate that the Kronecker models tend not to be a good fit to the data. The model appears to have a considerable difference in at least once of the graph features. Usually, it’s the number of triangles, which differs by up to two orders of magnitude for many of the collaboration networks.

3 Fitting partial sets of features

The previous set of experiments illustrated that the Kronecker graphs may not simultaneously fit all four of the network features: edges, hairpins/wedges, tripins, and triangles. In Table 4, we examine the change in fits when only using three of the four network features in the summation in the objective (12). We take the set of parameters with the smallest objective among all the procedures investigated in the previous section. The results show small changes to the parameters and expected feature fits. Nonetheless, the minimizer remains mostly unchanged.

Table 4 provides a kind of cross-validated feature estimation, showing the accuracy of a feature’s estimate when it is not included in the fitting. Apart from the exception noted before (the edge counts for as20000102) our moment based estimates give closer matches to the source feature counts than KronFit provides, whether the moment being studied is part of the fitting process or not.

We see some examples where leaving out one feature seems to improve the fitting of another. For instance, in three of the four graphs, leaving out the tripin count improved the match for triangles and conversely.

Synthetic examples

The results from the previous section show that there can often by a large deviation in the expected moments of the best Kronecker fit. In this section, we investigate the accuracy of the fitting procedure when the graph is a realization of a stochastic Kronecker network.

we generate 50 realizations of each Kronecker graph. For each realization, we compute a fit using the objective (12) with the choices $D_{\text{sq}},N_{F^{2}}$ and using the combination of approaches from the previous section. Figure 2 shows distribution of fitted parameters to these 50 samples. For all four sets of parameters, the fitted results closely match the true values, with fairly small variation.

For these synthetic problems, we also study how the empirical and fitted features differ. Figure 3 shows the distribution of the relative difference between the expectation of the fitted Kronecker features and the actual feature of each realization. It also shows the difference between the original feature count and the feature count of a re-realization. In other words, generate a Kronecker graph, fit the parameters, and re-generate with the fitted parameters. The figures show that the fitted parameters closely match the realizations. A curious property is that the fitted triangle count is always smaller than the empirical count. The difference in the re-realization can be large, almost 20% in the case of tripins or triangles for the first set of Kronecker parameters.

Our final study is the distributions of the graph features given the Kronecker parameters, the expected features of the fitted parameters, and the graph features of a re-realized Kronecker graph. The plots in Figure 4 show that these distributions are all quite similar.

Conclusions

We have presented formuals for expected feature counts in Kronecker graphs and used them to generate a method of moments fitting strategy. We found that summing squared relative feature count errors was robust and easy to optimize. For graphs generated by the Kronecker model, our parameter and feature estimates closely match those of the fitted graph. For real world graphs we often find that the fitted Kronecker model implies smaller feature counts (apart from edges) than are seen in the real graph. The moment estimators typically come closer to the counts than those from KronFit.

Acknowledgments

We thank Tamara Kolda, C. Seshadhri, and Ali Pinar for helpful discussions. This work was supported by DMS-0906056 of the National Science Foundation.

References