Fractality of Massive Graphs: Scalable Analysis
with Sketch-Based Box-Covering Algorithm∗
Takuya Akiba† , Kenko Nakamura‡ , Taro Takaguchi§
†

§

Preferred Networks, Inc. ‡ Recruit Communications Co., Ltd.
National Institute of Information and Communications Technology

arXiv:1609.07994v1 [] 26 Sep 2016

akiba@preferred.jp, kenkoooo@r.recruit.co.jp, taro.takaguchi@nict.go.jp

Abstract—Analysis and modeling of networked objects are
fundamental pieces of modern data mining. Most real-world
networks, from biological to social ones, are known to have
common structural properties. These properties allow us to
model the growth processes of networks and to develop useful
algorithms. One remarkable example is the fractality of networks,
which suggests the self-similar organization of global network
structure. To determine the fractality of a network, we need
to solve the so-called box-covering problem, where preceding
algorithms are not feasible for large-scale networks. The lack
of an efficient algorithm prevents us from investigating the
fractal nature of large-scale networks. To overcome this issue,
we propose a new box-covering algorithm based on recently
emerging sketching techniques. We theoretically show that it
works in near-linear time with a guarantee of solution accuracy.
In experiments, we have confirmed that the algorithm enables us
to study the fractality of million-scale networks for the first time.
We have observed that its outputs are sufficiently accurate and
that its time and space requirements are orders of magnitude
smaller than those of previous algorithms.

I. I NTRODUCTION
Graph representation of real-world systems, such as social relationship, biological reactions, and hyperlink structure,
gives us a strong tool to analyze and control these complex
objects [21]. For the last two decades, we have witnessed
the spark of network science that unveils common structural
properties across a variety of real networks. We can exploit
these frequently observed properties to model the generation
processes of real networked systems [20] and to develop
graph algorithms that are applicable to various objects [19].
A notable example of such properties is the scale-free property [3], [7], which manifests a power-law scaling in the vertex
degree distribution and existence of well-connected vertices
(often called hubs). The scale-free property, existence of hubs
especially, underlies efficient performance of practical graph
algorithms on realistic networks [1], [17].
Although the scale-free property inspires us to design better
network models and algorithms, it is purely based on the
local property of networks, i.e., the vertex degree. Real-world
networks should possess other common properties beyond the
local level. As a remarkable example of such non-local properties, the fractality of complex networks was found in network
science [27], [14]. The fractality of a network suggests that
∗ This work is done while all authors were at National Institute of
Informatics. A shorter version of this paper appeared in the proceedings of
ICDM 2016 [2].

the network shows a self-similar structure; if we replace
groups of adjacent vertices with supervertices, the resultant
network holds a similar structure to the original network
(see Section II-C for its formal definition). The fractality of
networks gives us unique insights into modeling of growth
processes of real-world networks [28]. In addition, fractal and
non-fractal networks, even with the same degree distribution,
indicate striking differences in facility of spreading [25] and
vulnerability against failure [16]. Aside from theoretical studies, the fractality provides us with useful information about
network topology. Examples include the backbone structure of
networks [15] and the hierarchical organization of functional
modules in the Internet [8], metabolic [28] and brain [13]
networks, to name a few.
Determination of the fractality of a network is based on the
so-called box-covering problem [27] (also see Section II-C).
We locally cover a group of adjacent vertices with a box such
that all vertices in a box are within a given distance from
each other, and then we count the number of boxes we use to
cover the whole network. In principle, we have to minimize
the number of boxes that cover the network, which is known
to be an NP-hard problem (see [26] and references therein).
Although different heuristic algorithms are proposed in the
previous work (e.g., [26], [24]), they are still not so efficient
as to be able to process networks with millions of vertices. This
limitation leaves the fractal nature of large-scale networks far
from our understanding.
Contributions: The main contribution of the present study
is to propose a new type of box-covering algorithm that is
much more scalable than previous algorithms. In general,
previous algorithms first explicitly instantiate all boxes and
then reduce the box cover problem to the famous set cover
problem. This approach requires quadratic Θ(n2 ) space for
representing neighbor sets and is obviously infeasible for
large-scale networks with millions of vertices. In contrast, the
central idea underlying the proposed method is to solve the
problem in the sketch space. That is, we do not explicitly
instantiate neighbor sets; instead, we construct and use the
bottom-k min-hash sketch representation [9], [10] of boxes.
Technically, we introduce several new concepts and algorithms. First, to make the sketch-based approach feasible,
we introduce a slightly relaxed problem called (1 − )B OX C OVER. We also define a key subproblem called the

(1 − )-S ET C OVER problem. The proposed box-cover algorithm consists of two parts. First, we generate min-hash
sketches of all boxes to reduce the (1−)-B OX C OVER problem
to the (1 − )-S ET C OVER problem. Our sketch generation
algorithm does not require explicit instantiation of actual boxes
and is efficient in terms of both time and space. Second, we
apply our efficient sketch-space set-cover algorithm to obtain
the final result. Our sketch-space set-cover algorithm is based
on a greedy approach, but is carefully designed with eventdriven data structure operations to achieve near-linear time
complexity.
We theoretically guarantee both the scalability and the
solution quality of the proposed box-cover algorithm. Specifically, for a given trade-off parameter k and radius parameter
`, it works in O((n + m)k log k min {`, log n}) time and
O(nk +m) space. The produced result is a solution of (1−)B OX C OVER within a factor
p 1 + 2 ln n of the optimum for
B OX C OVER for  ≥ 2 5(ln n)/k, with a high probability
that asymptotically approaches 1.
In experiments, we have confirmed the practicability of
the proposed method. First, we observed that its outputs are
quite close to those of previous algorithms and are sufficiently
accurate to recognize networks with ground-truth fractality.
Second, the time and space requirements are orders of magnitude smaller than those of previous algorithms, resulting in
the capability of handling large-scale networks with tens of
millions of vertices and edges. Finally, we applied our algorithm to a real-world million-scale network and accomplish its
fractality analysis for the first time.

TABLE I
F REQUENTLY USED NOTATIONS .

Notation

Description

(In
G = (V, E)
n, m
Nδ (v)
{Sp }p∈P
n
k
ri
Se
e S)
e
C(

the context of the box cover problem)
The graph.
The numbers of vertices and edges in G.
The vertices with distance at most δ from v.

(In the context of the set cover problem)
The set family.
The number of elements and collections.
(Bottom-k min-hash sketch)
The trade-off parameter of min-hash sketches.
The rank of an item i.
The min-hash sketch of set S.
The estimated cardinality of set S.

set of items. We first assign a random rank value ri ∼ U (0, 1)
to each item i ∈ X, where U (0, 1) is the uniform distribution
on (0, 1). Let S be a subset of X. For an integer k ≥ 1,
e where
the bottom-k min-hash sketch of S is defined as S,
e
i ∈ S ⇐⇒ ri ≤ k-th {rj | j ∈ S} . In other words, Se is
the set of vertices with the k smallest rank values. We define
Se = S if |S| < k.
For a set S ⊆ X, the threshold rank τ (S) is defined as
follows. If |S| ≥ k, τ (S) = k-th {ri | i ∈ S}. Otherwise,
e Using sketch S,
e
τ (S) = (k − 1)/ |S|. Note that τ (S) = τ (S).
e
e
we estimate the cardinality |S| as C(S) = (k − 1)/τ (S). Its
relative error is theoretically bounded as follows.

Organization: The remainder of this paper is organized as
follows. We describe the definitions and notations in Section II.
In Section III, we present our algorithm for sketch-space
S ET C OVER. We explain our sketch construction algorithm to
complete the proposed method for B OX C OVER in Section IV.
In Section V, we present a few empirical techniques to further
improve the proposed method. We explain the experimental
evaluation of the proposed method in Section VI. We conclude
in Section VII.

Lemma 1 (Bottom-k cardinality estimator [9], [10]): The
e
cardinality estimation C(S)
is an unbiased estimator√of |S|,
and has a coefficient of variation (CV)1 of at most 1/ k − 2.

II. P RELIMINARIES

In addition, our algorithms heavily rely on the mergeability
of min-hash sketches. Suppose S1 , S2 ⊆ X and S3 = S1 ∪ S2 .
Then, since Se3 ⊆ Se1 ∪ Se2 , Se3 can be obtained only from Se1
and Se2 . We denote this procedure as Merge-and-Purify (e.g.,
Se3 = Merge-and-Purify(Se1 , Se2 )).
For simplicity, we assume that ri is unique for i ∈ X,
and sometimes identify i with ri . In particular, we use the
comparison between elements such as i < j for i, j ∈ X,
where we actually compare ri and rj . We also define k-th(S)
as the element with the k-th smallest rank in S ⊆ X.

A. Notations
We focus on networks that are modeled as undirected
unweighted graphs. Let G = (V, E) be a graph, where V
and E are the vertex set and edge set, respectively. We use
n and m to denote |V | and |E|, respectively. For d ≥ 0 and
v ∈ V , we define Nd (v) as the set of vertices with distance at
most d from v. We call Nd (v) the d-neighbor. When d = 1,
we sometimes omit the subscript, i.e., N (v) =SN1 (v). We also
define Nd (S) for a set S ⊆ V as Nd (S) = v∈S Nd (v). In
other words, Nd (S) represents the set of vertices with distance
at most d from at least one vertex in S. The notations we will
frequently use hereafter are summarized in Table I.
B. Bottom-k Min-Hash Sketch
In this subsection, we review the bottom-k min-hash sketch
and its cardinality estimator [9], [10]. Let X denote the ground

The following corollary can be obtained by applying Chernoff bounds [11].
Corollary 2: For  > 0 and c > 1, by setting k ≥ (2 +
c)−2 ln |X|, the probability of the estimation having a relative
c
error larger than  is at most 1/ |X| .

C. Problem Definition
1) Graph Fractality: The fractality of a network [27] is a
generalization of the fractality of a geometric object in Euclidean space [12]. A standard way to determine the fractality
of a geometric object is to use the so-called box-counting
1 The

CV is the ratio of the standard deviation to the mean.

106
10

∝
∝

5

b ( `)

104

2) Box Cover: As we described in the previous section, the
fractality of graphs is analyzed by solving the box covering
problem. The problem has two slightly different versions: the
diameter version [27] and the radius version [26]. It has been
empirically shown that these two versions yield negligible
difference in the results. In this study, we focus on the radius
version, which is defined as follows.

`−1.6
e−0.3 `

103
102
101
100 0
10

101

102

`

103

104

Fig. 1. Fitting of b(`) to the power-law and exponential functions for a
synthetic fractal network.

method; we tile the object with cubes of a fixed length and
count the number of cubes needed. If the number of cubes
follows a power-law function of the cube length, the object is
said to be fractal. A fractal object holds a self-similar property
so that we observe similar structure in it when we zoom in
and out to it.
The idea of the box counting method is generalized to
analyze the fractality of networks [27]. The box-covering
method for a network works by covering the network with
boxes of finite length `, which refers to a subset of vertices in
which all vertices are within distance `. For example, a box
with ` = 1 is a set of nodes all adjacent to each other. If
the number of boxes of length ` needed to cover the whole
network, denoted by b(`), follows a power-law function of `:
b(`) ∝ `−d , the network is said to be fractal. The exponent d is
called the fractal dimension. As can be noticed, b(`) crucially
depend on how we put the boxes. Theoretically, we have to
put boxes such that b(`) is minimized to assess its precise
scaling. However, this box-covering problem is NP-hard and
that is why we propose our new approximation algorithm of
this problem in the rest of this paper.
After computing b(`) for a network, we want to decide
whether the network is fractal or not. A typical indicator of
a non-fractal network is an exponential form: b(`) ∝ exp(c`),
where c is a constant factor [29]. Therefore, comparison of
the fitting of the obtained b(`) to power-law and exponential functions enables us to determine the fractality of the
network. Figure 1 illustrates the comparison for (3, 3, 7)flower, a network model with ground-truth fractality [23] (see
Section VI-A). Since b(`) is closer to the power-law than to the
exponential function in this case, the fitting procedure correctly
indicates the fractality of this network model.
It should be noted that the fractality of a network suggests
its self-similarity. Let us aggregate the vertices in a box into a
supervertex and then aggregate the edges spanning two boxes
into a superedge. Then we obtain a coarse-grained version
of the original network. If the original network is fractal,
the vertex degree distributions of the original and coarsegrained networks are (statistically) the same [27]. Note that
the fractality and self-similarity of a scale-free network are
not equivalent, and a non-fractal scale-free network can be
self-similar under certain conditions [18].

Problem 1 (B OX C OVER): In the B OX C OVER problem, given
a graph G and a radius limit ` > 0, the objective is to find a
set S ⊆ V of the minimum size such that N` (S) = V .
The size of set S is equal to b(`) discussed in Section II-C1.
In this study, we consider a slightly relaxed variant of the
B OX C OVER problem, named (1−)-B OX C OVER. The (1−)B OX C OVER problem is defined as follows.
Problem 2 ((1 − )-B OX C OVER): In the (1 − )-B OX C OVER
problem, we are given a graph G, a radius limit ` > 0 and
an error tolerance parameter  > 0. The objective is to find a
set S ⊆ V of the minimum size such that |N` (S)| ≥ (1 − )n.
3) Set Cover: The B OX C OVER problem is a special case
of the S ET C OVER problem, which is defined as follows.
Problem 3 (S ET C OVER): In the S ET C OVER problem, we are
given a set family {Sp }p∈P . The objective
is to S
find a set
S
R ⊆ P of the minimum size such that p∈R Sp = p∈P Sp .
The proposed box-covering algorithm deals with a slightly
different version of S ET C OVER, named (1 − )-S ET C OVER
with sketched input as a key subproblem, which is defined as
follows.
Problem 4 ((1 − )-S ET C OVER with sketched input): In the
sketched input version of the (1 − )-S ET C OVER problem, we
are given the min-hash sketches of a set family {Sep }p∈P and
an error tolerance parameter  > 0. The objective
S is to find
a set R S⊆ P of the minimum size such that | p∈R Sp | ≥
(1 − )| p∈P Sp |.
We first design an efficient approximation algorithm for
(1 − )-S ET C OVER (Section III). We then propose a new boxcovering algorithm using it (Section IV).
III. S ET C OVER IN S KETCH S PACE
In this section, we design an efficient approximation algorithm for the sketched-input version of (1 − )-S ET C OVER
(Problem 4). We call each p ∈ P a collection and i ∈ Sp
an element. Because of the connection to the B OX C OVER
problems, we assume that the numbers of collections and
elements
are equal. We denote them by n, that is,
S
S |P | =
| p∈P Sp | = n. For R ⊆ P , we define SR = p∈R Sp .
e R ) by C(R),
e
Moreover, for simplicity, we denote C(S
which
can be calculated from merged min-hash sketch SeR .
We first explain the basic greedy algorithm that runs in
O(n2 k) time, and then present its theoretical solution guarantee. Finally, we propose an efficient greedy algorithm, which
runs in O(nk log n) time and produces the exact same solution
as the basic algorithm.

Algorithm 1 Select-Greedily-Naive({Sep }p∈P )

At the i-th iteration, there is some collection p such that


1
n − C(Ri ∪ {p}) ≤ 1 − ∗ (n − C(Ri )) ,
(1)
|R |

eR ← {}
1: R ← {} , S
e
2: while R 6= P and C(R)
< (1 − /2)n do
e ∪ {p}) | p ∈ P \ R}
3:
p ← argmax{C(R
4:
R ← R ∪ {p} , SeR ← Merge-and-Purify(SeR , Sep )
5: return R

A. Basic Greedy Algorithm
Our basic greedy algorithm Select-Greedily-Naive is described as Algorithm 1. We start with an empty set R = {}. In
e ∪ {p}) for every p ∈ P \ R,
each iteration, we calculate C(R
and select p that maximizes the estimated cardinality, and add
e
it to R. We repeat this until C(R)
gets at least (1 − /2)n,
and the resulting R is the solution.
e ∪ {p}), together with R, we manage the
To calculate C(R
merged min-hash sketch SeR , so that SeR always corresponds
to the min-hash sketch of SR . To this end, we use the merger
operation of min-hash sketch. Let us assume that the items in
a min-hash sketch are stored in the ascending order of their
ranks. Then, merging two min-hash sketches can be done in
O(k) time like in the merge sort algorithm; we just need to
pick the top-k distinct items with the lowest ranks in the two
min-hash sketches. The complexity analysis of the algorithm
is as follows.
Lemma 3: Algorithm Select-Greedily-Naive
O(n2 k) time and O(nk) space.

runs

in

Proof Sketch. This algorithm always terminates as, even in the
worst case, after the n-th iteration, R gets P . Each iteration
takes O(nk) time, and the number of iterations is at most n.
Therefore, the time complexity is O(n2 k) time.
B. Theoretical Solution Guarantee
We can guarantee the quality of the solution produced by
the above algorithm as follows.
p
Lemma 4: For  ≥ 2 5(ln n)/k, algorithm SelectGreedily-Naive produces a solution of (1 − )-S ET C OVER
within a factor 1 + 2 ln n of the optimum for S ET C OVER with
a probability of at least 1 − 1/n.
In other words, with a high probability that asymptotically
approaches 1, R is the solution of (1 − )-S ET C OVER and
|R| ≤ (1 + 2 ln n)|R∗ |, where R is the output of algorithm
Select-Greedily-Naive and R∗ is the optimum solution of
S ET C OVER (with the same set family as the input).
Proof. Let R be the output of algorithm Select-GreedilyNaive, and R∗ be the optimum solution to S ET C OVER. Let
R0 = ∅ and Ri ⊆ P be the currently selected sets after
the i-th iteration
of algorithm Select-Greedily-Naive. Let
S
C(R) = | p∈R Sp |. From Corollary 2, for set R ⊆ P , the
e
probability of C(R)
having a relative error larger than /2 is
3
at most 1/n .

since otherwise there would not be any solution R∗ of that
size to S ET C OVER. During the i-th iteration, there are at
most n new sets to be examined, and thus the union bound
e is at most /2
implies that the relative error between C and C
2
with a probability of at least 1 − 1/n . Therefore, with that
probability,


e i+1 ) ≥ C(R
e i ∪ {p}) ≥ 1 −  C(Ri ∪ {p}).
(2)
C(R
2
From Inequalities 1 and 2, with some calculation, we have

 



1

e
e i) ,
1−
n−C(Ri+1 ) ≤ 1 − ∗
n − C(R
1−
2
|R |
2
with a probability of at least 1 − 1/n2 . As the number of
iterations is at most n, by applying the union bound over all
iterations, we obtain

i

i

1
e
1−
n − C(Ri ) ≤ 1 − ∗
n < e− |R∗ | n,
2
|R |
with a probability of at least 1−1/n. If i is at least 2|R∗ | ln n,
it becomes strictly less than 1/n, which is smaller than the
e Therefore, the number of iterations is at most
resolution of C.
∗
d2|R | ln ne, and thus |R| ≤ d2|R∗ | ln ne ≤ (1 + 2 ln n)|R∗ |.
e
Moreover, C(R) ≥ (1 − /2)C(R)
≥ (1 − )n, and thus R is
the solution to the (1 − )-S ET C OVER problem.
C. Near-Linear Time Greedy Algorithm
Algorithm Select-Greedily-Naive takes quadratic time,
which is unacceptable for large-scale set families. Therefore,
we then design an efficient greedy algorithm Select-GreedilyFast, which produces the exact same output as algorithm
Select-Greedily-Naive but runs in O(nk log n) time. As the
input size is O(nk), this algorithm is near-linear time.
The behavior of Select-Greedily-Fast at a high level is
the same as that of Select-Greedily-Naive. That is, we start
with an empty set R = {}, and, at each iteration, it adds
e to R. The central
p ∈ P \ R with the maximum gain on C
idea underlying the speed-up is to classify the state of each
p ∈ P at each iteration into two types and manage differently
to reduce the reevaluation of the gain. To this end, we closely
look at the relation between sketches Sep and SeR .
Types and Variables: Let us assume that we are in the main
loop of the greedy algorithm. Here, we have a currently
incomplete solution R ⊂ P . Let p ∈ P \ R. We define
that p belongs to type A if the k-th element of SeR∪p is in
Sp , i.e., k-th(SeR∪p ) ∈ Sp . Otherwise, p is type B. Note that
SeR∪p = Merge-and-Purify(SeR , Sep ).
We define ap = |SeR∪p ∩ Sep |, bp = |SeR∪p ∩ SeR |, and
cp = |SeR ∩ Sep |. Please note that, if p is a type-A collection,
e
C(R
∪ {p}) is determined by k-th(SeR∪p ) = ap -th(Sep ).
e
Similarly, if p is a type-B collection, C(R∪{p})
is determined
e
e
by k-th(SR∪p ) = bp -th(SR ).

Algorithm 2 Select-Greedily-Fast({Sep }p∈P )
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:

// Initialization
R ← {} , SeR ← {}
QA ← an empty min-queue (key: ranks, value: collections)
QB ← an empty min-queue (key: integers, value: collections)
for j = 1, 2, . . . , k do
T (j) ← a binary search tree (key: ranks, value: collections)
Ii ← {p ∈ P | i ∈ Sep }
for all p ∈ P do
Insert (p, k-th{ri | i ∈ Sep }) to QA and T (1)
ap ← k, bp ← 0
// Main loop
e
while R 6= P and C(R)
< (1 − /2)n do
// Selection
e ∪ {p}) | p is at the top of QA or QB }
p ← argmax{C(R
R ← R ∪ {p}
Remove p from QA , QB and T
0
SeR
← Merge-and-Purify(SeR , Sep )
0
0
∆ ← SeR
\ SeR , SeR ← SeR
// Notifying events of type 3
for all i ∈ ∆ do
for all p ∈ Ii do
if p ∈ QB then
Move p from T (bp ) to T (bp + 1)
else if ap -th(Sep ) = i then
while ap -th(Sep ) = i do ap ← ap − 1
Update p’s key in T (bp + 1) to ap -th(Sep )
Remove p from QA and Insert (p, bp + 1) to QB
else
Move p from T (bp + 1) to T (bp + 2)
bp ← bp + 1
// Notifying events of types 1 and 2-1
for j = 1, 2, . . . , k do
P 0 ← Retrieve those with keys ≥ j-th(SeR ) from T (j)
for all p ∈ P 0 do
r ← j-th(SeR )
Remove p from QA , QB and T (j)
if p ∈ QA then
ap ← ap − 1, bp ← bp + 1
if ap -th(Sep ) ∈ SeR then
bp ← bp − 1, r ← (j − 1)-th(SeR )
while ap -th(Sep ) ∈ SeR do ap ← ap − 1
if ap -th(Sep ) > r then
Insert (p, ap -th(Sep )) to QA and T (bp + 1)
else
Insert (p, bp ) to QB
Insert (p, ap -th(Sep )) to T (bp )
return R

Events to be Captured: Suppose that we have decided to
adopt a new collection and R is about to be updated to R0
(i.e., R0 = R ∪ {p0 } for some collection p0 ∈ P ). Let us first
assume that a single element appeared in the merged sketch,
i.e., SeR0 \ SeR = {i}. Let p ∈ P \ R0 . In the following, we
examine and classify the events where the evaluation of p is
e ∪ {p}) 6= C(R
e 0 ∪ {p}) (types 1 and 2), or
updated, i.e., C(R
cp is updated (type 3).
Type 1: We assume that i 6∈ Sep and p is type A. From the
definition, τ (SeR∪p ) = ap -th(Sep ), and, bp -th(SeR ) ≤ ap -th(Sep )
< (bp + 1)-th(SeR ). Therefore, τ (SeR∪p ) 6= τ (SeR0 ∪p ) if and
only if (bp +1)-th(SeR0 ) 6= (bp +1)-th(SeR ) and (bp +1)-th(SeR0 )

< ap -th(Sep ). We define that a type-1 event happens to p when
this condition holds.
Type 2: Similarly, we assume that i 6∈ Sep and p is
type B. From the definition, τ (SeR∪p ) = bp -th(SeR ) and
ap -th(Sep ) < bp -th(SeR ). Thus τ (SeR∪p ) 6= τ (SeR0 ∪p ), if
and only if bp -th(SeR ) 6= bp -th(SeR0 ). There are two cases:
bp -th(SeR0 ) ≤ ap -th(Sep ) (type 2-1), after which p becomes
type A, or bp -th(SeR0 ) > ap -th(Sep ) (type 2-2), after which p
still belongs to type B.
Type 3: If i ∈ Sep , then cp will be incremented.
The following lemma is the key to the efficiency of our
algorithm.
Lemma 5: For each p ∈ P , throughout the algorithm execution, events of type 1, type 2-1, or type 3 occur at most 3k
times in total.
Proof Sketch. We use the progress indicator Φ = k − ap +
bp + cp . Initially, ap = k and bp = cp = 0; hence Φ = 0. For
each event occurrence, Φ increases by at least one. As ap ≥ 0
and bp , cp ≤ k, Φ ≤ 3k.
Please note that events of type 2-2 are not considered in
the above lemma, and, indeed, they happen Θ(n) times in
the worst case for each collection. Therefore, we design the
algorithm so that we do not need to capture type-2-2 events.
Finding the Maximum Gain: To adopt a new collection in
each iteration, we need to efficiently find the collection that
gives the maximum gain. We clarify the ordering relation in
each type.
e ∪
Type A: Let p, q ∈ P \ R be type-A collections, then C(R
e
e
e
{p}) ≥ C(R ∪ {q}) if and only if ap -th(Sp ) ≤ aq -th(Sq ).
e ∪
Type B: Let p, q ∈ P \ R be type-B collections, then C(R
e ∪ {q}) if and only if bp -th(SeR ) ≤ bq -th(SeR ),
{p}) ≥ C(R
which is equivalent to bp ≤ bq .
Data Structures: We use the following data structures to
notify the collections about an event occurrence.
Type 1: For each type-A collection p, as we observed above,
p wants to be notified about a type-1 event when (bp +
1)-th(SeR ) becomes smaller than ap -th(Sep ). Therefore, for
j = 1, 2, . . . , k, we prepare a binary search tree T (j), where
values are collections and keys are ranks (i.e., collections are
managed in the ascending order of ranks in each tree). For
each type-A collection p, we put p in T (bp + 1) with key
ap -th(Sep ). Then, when j-th(SeR ) is updated to a new value,
from T (j), we retrieve collections with keys larger than or
equal to the new value and notify them about an event.
Type 2-1: Similarly, for each type-B collection p, p wants to
be notified about a type-2-1 event when bp -th(SeR ) becomes
smaller than or equal to ap -th(Sep ). Thus, we store p in T (bp )
and set its key to ap -th(Sep ). Then, when j-th(SeR ) is updated
to a new value, we retrieve those in T (j) with keys larger than
or equal to the new value and notify them about an event.

Type 3: To capture type-3 events, the use of an inverted index
suffices. That is, for each i ∈ X, we precompute Ii = {p ∈
P | i ∈ Sep }. When i comes to SeR , we notify the collections
in Ii .
Moreover, we also need data structures to find the collections with the maximum gain as follows.
Type A: Type-A collections are managed in a minimumoriented priority queue, where the key of a collection p is
ap -th(Sep ).
Type B: Type-B collections are managed in another minimumoriented priority queue, where the key of a collection p is bp .
Overall Set-Cover Algorithm: The overall algorithm of
Select-Greedily-Fast is described as Algorithm 2. In each
iteration, we adopt the new collection with maximum gain,
which can be identified by comparing the top elements of
the two priority queues. Then, we process events to update
variables and data structures. At the beginning of Section III-C,
we assumed that a single element appears in the new sketch.
When more than one elements come to the new sketch, we
basically process each of them separately. See Algorithm 2 for
the details of the update procedure. The algorithm complexity
and solution quality are guaranteed as follows.
Lemma 6: Algorithm Select-Greedily-Fast
O(nk log n) time and O(nk) space.

runs

in

Proof Sketch. Each data structure operation takes O(log n)
time, which, from Lemma 5, happens at most 3k times for
each collection.
Lemma 7: Algorithm Select-Greedily-Fast produces the
same solution as algorithm Select-Greedily-Naive.

Algorithm 3 Build-Sketches(G, `)
1: Xv ← {v} for all v ∈ V .
2: for ` times do
3:
for all v ∈ V in the increasing order of rv do
4:
Av ← {u ∈ V | v is added to Xu in the last iteration}
5:
for all w ∈ N (Av ) do
6:
Xw ← Merge-and-Purify(Xw , {v})
7:
if Xv was not modified for any v ∈ V then
8:
break
9: return {Xv }v∈V

In each iteration, for each vertex, we essentially merge the
sketches of its neighbors into its sketch in a message-passinglike manner. Two speed-up techniques are employed here to
avoid an unnecessary insertion check. For v ∈ V , let Av
be the vertices in whose sketches v is added to in the last
iteration. First, for each v ∈ V , we try to insert v only into
the sketches of the vertices that are neighbors of Av , as v
cannot be inserted into other vertices. Second, we conduct the
procedure above in the increasing order of ranks, since this
decreases the unnecessary insertion. We prove its correctness
and complexity as follows.
Lemma 8: In algorithm Build-Sketches, after the i-th iterei (v) for all v ∈ V .
ation, Xv = N
Proof Sketch. We prove the lemma by mathematical induction
e0 (v), it is true for i = 0. Now we assume
on i. Since {v} = N
it holds for i and prove it also holds for i + S
1. Let B =
{u ∈ V | (v, u) ∈ E}. Since Ni+1 (v) = {v} ∪ u∈B Ni (u),
ei+1 (v) can be obtained by
and {v} ∈ Ni (v) ⊆ Ni+1 (v), N
ei (u) for all u ∈ B ∪ {v}.
merging N

Proof Sketch. Both algorithms choose the collection with the
maximum gain in each iteration.

e` (v) for
Corollary 9: Algorithm Build-Sketches computes N
all v ∈ V .

IV. S KETCH -BASED B OX C OVERING

Lemma 10: Algorithm Build-Sketches runs in O((n + m)k
log k min {`, log n}) expected time and O(nk + m) space.

In this section, we complete our sketch-based box-covering
algorithm for the (1 − )-B OX C OVER problem (Problem 2).
We first propose an efficient algorithm to construct min-hash
sketches representing the `-neighbors, and then present and
analyze the overall box-covering algorithm,
A. Sketch Generation
For v ∈ V , we denote the min-hash sketch of N` (v) as
e` (v). Here, we construct N
e` (v) for all vertices v ∈ V
N
to reduce the (1 − )-B OX C OVER problem to the (1 − )S ET C OVER problem (Problem 4). Our sketch construction
algorithm Build-Sketches is described as Algorithm 3.
It receives a graph G and a radius parameter `. Each vertex
v manages a tentative min-hash sketch Xv . Initially, Xv only
includes the vertex itself, i.e., Xv = {v}, which corresponds to
e0 (v). Then, we repeat the following procedure for ` times so
N
ei (v). This algorithm has
that, after the i-th iteration, Xv = N
a similar flavor to algorithms for approximated neighborhood
functions and all-distances sketches [22], [5], [10].

Proof Sketch. In addition to the graph, the algorithm stores a
sketch of size k for each vertex, and hence it works in O(nk +
m) space. Each insertion trial takes O(log k) time (Line 6).
Therefore, it suffices to prove that the number of traversed
edges is O((n + m)k`) and O((n + m)k log n). The former
bound is easier, since, in each iteration, the number of last
inserted elements in each sketch is at most k, and thus we
traverse each edge at most k times.
For the latter bound, we count the expected number of
vertices that are inserted once into Xv for a vertex v ∈ V .
The vertex that is i-th to arrive at v is inserted into Xv with
a probability of min {1, k/i}, and thus it is at most


n
X
k
min 1,
= k + k(H(n) − H(k)) = O(k log n),
i
i=1
where H(i) is the i-th Harmonic number. Therefore, each edge
is traversed at most O(k log n) times in total.

B. Overall Box-Cover Algorithm
The overall box-covering algorithm Sketch-Box-Cover is
as follows. We first construct the min-hash sketches using algorithm Build-Sketches and then solve the set cover problem
in the sketch space using algorithm Select-Greedily-Fast.
The guarantees on performance and accuracy are immediate
from the previous lemmas and corollaries as below.
Theorem 11 (Scalability guarantee): Algorithm Sketch-BoxCover works in O((n + m)k log k min {`, log n}) time and
O(nk + m) space.
Theorem 12 (Solution accuracy guarantee):
p With a probability of at least 1 − 1/n, for  ≥ 2 5(ln n)/k, algorithm Sketch-Box-Cover produces a solution to the (1 − )B OX C OVER problem within a factor 1 + 2 ln n of the optimum
for the B OX C OVER problem.
Assuming k is a constant, the time and space complexities
are near-linear. Similarly, given a constant , the time and
space complexities are still near-linear, since it suffices to set
k = d20−2 ln ne. In practice, as seen in our experiments,
the algorithm produces solutions that are much closer to the
optimum than what is expected from the above approximation
ratio with much smaller k.
V. P RACTICAL I MPROVEMENT
In this section, we propose techniques to improve the
practicality of the proposed method.
Exact Coverage Management: For the termination condition
in the greedy selection algorithm (i.e., Line 12 in Algorithm 2),
when applied to the box cover problem, we propose to use
the exact coverage C(R) instead of the estimated coverage
e
C(R).
This technique makes the results more stable. We can
efficiently manage the exact coverage as follows.
First, we prepare an array δ, and initialize it as δ[v] = ∞
for all v ∈ V . After selecting a vertex v in each iteration, we
conduct a pruned breadth-first search (BFS) from v. Suppose
we are visiting vertex u with distance d in this BFS. If
δ[u] ≤ d, then we prune this BFS, i.e., we do not traverse the
edges from u. Otherwise, we set δ[u] = d and continue the
search. We do not visit vertices with a distance larger than `.
The number of covered vertices is the number of non-infinity
values in array δ. Since the value of δ[u] changes at most `+1
times, each vertex or edge is visited O(`) times. Therefore, the
total time consumption of this process throughout all iterations
is O((n + m)`).
Multi-Pass Execution: On the basis of the above exact coverage management technique, we sometimes detect that, even
e
e )),
while the estimated coverage is saturated (i.e., C(R)
= C(P
the actual coverage is below the specified threshold. In that
case, to choose more vertices, we propose to repeat the
algorithm from sketch construction until the actual coverage
becomes higher than the threshold.
In the i-th pass, we only care for vertices that are not
covered by the previous passes. This can be easily realized

by modifying the algorithm Build-Sketches so that, at Line
1, we set Xv = ∅ for already covered vertices. For accurate
results, node ranks should be reassigned for each pass.
Exact Neighborhood: To further improve the accuracy, we
propose to combine our sketch-based algorithm with a nonsketch-based algorithm. For a very small radius parameter
`, neighborhood N` (v) is sometimes much smaller than k.
Moreover, even for a larger `, when using the above multi-pass
execution technique, the remaining neighbors may become
small in later passes. In these cases, the sketching approach has
little advantage. Therefore, we detect such circumstances and
switch to a non-sketch-based greedy algorithm. Interestingly,
e` (v)| ≤ k, then
this switching can be done seamlessly. If |N
e` (v) = N` (v). Therefore, under such circumstances, the
N
output of algorithm Build-Sketches can be immediately given
to the non-sketch-based greedy algorithm.
The proposed overall procedure is as follows. We specify
a parameter α. We start by constructing the “sketches” with
algorithm Build-Sketches, but, at first, we apply the algorithm
as if k = ∞, i.e., we do not conduct purification on the minhash sketches. During the construction, if the total number of
elements in all min-hash sketches exceeds αnk at some point,
then we conduct purification on all the min-hash sketches,
continue the construction with the actual k value, and pass
the resulting sketches to the sketch-based greedy algorithm.
Otherwise, we apply the non-sketch-based greedy algorithm
to the resulting “sketches,” which are actually exact neighborhood sets. Assuming parameter α is a constant, the total time
and space complexity remain the same.
Exact Box Covering: Together with the preceding three techniques, to further make the results reliable, we propose to use
our algorithm for solving the original B OX C OVER problem
(Problem 1) rather than (1 − )-B OX C OVER problem (Problem 2). In other words, we recommend setting  = 0 to ensure
that all vertices are completely covered. As we will see in
the experiments, even with this seemingly extreme threshold,
thanks to the above techniques, both the running time and the
solution quality are reasonable.
VI. E XPERIMENTS
In this section, we present our experimental results to verify
the performance of the algorithm. Specifically, we compare our
algorithm with other preceding algorithms in terms of accuracy
and computation time.
We mainly focus on model networks, instead of empirical
ones, in order to validate the results of our algorithm with
ground-truth theoretical solutions and to investigate the scalability of the algorithm for various sizes of networks. However,
we also demonstrate the practicalness of our algorithm by
applying it to a real million-scale web graph. On the basis
of the result, we reveal the fractality of such large-scale real
graph for the first time.
A. Setup
Environment: Experiments were conducted on a Linux server
with Intel Xeon X5650 (2.67 GHz) and 96GB of main

approximation ratio

(3, 0, 6)-SHM
(3, 3, 7)-flower
(2, 0, 8)-SHM

24 25 26 27 28 29 210

k

1.8
1.6
1.4
1.2
1.0

approximation ratio

1.8
1.6
1.4
1.2
1.0

(3, 0, 6)-SHM
(3, 3, 7)-flower
(2, 0, 8)-SHM

2-3 2-2 2-1 20 21 22 23

Fig. 2. Average approximation ratio to the theoretical solutions for various
k and α.

memory. Algorithms were implemented in C++ and compiled
by gcc 4.8.4 with -O3 option.
Algorithms: For comparison, we used a naive algorithms
named greedy coloring (GC) and three advanced and popular algorithms, named maximum excluded mass burning
(MEMB), minimal value burning (MVB), and compact box
burning (CBB). GC, MEMB, and CBB were introduced in
[26] and MVB was in [24].
Network Models: We used two network models with groundtruth fractality: the (u, v)-flower [23] and the Song-HavlinMakse (SHM) [28] model. These models have power-law degree distributions, the representative characteristic of complex
networks. Both models can be either fractal or non-fractal,
depending on the structural parameter values. We refer to them
as the (u, v, g)-flower and (c, e, g)-SHM model to indicate the
parameter settings. The common parameter g (g = 1, 2, 3, . . . )
determines the network size n: n = (w − 2/w − 1)wg +
w/w − 1, where w ≡ u + v for the (u, v, g)-flower, and
n = (2c+1)g n0 for the (c, e, g)-SHM model. In addition to the
flower and SHM models, we considered the Barabási-Albert
(BA) network model [3] as one of the most famous models of
complex networks. The BA model is not fractal [27]. We refer
to this model as (c, t)-BA, where c is the number of edges that
a new node has and t sets the network size as n = 125 × 2t .
Fractality Decision Procedure: After the computation of the
box-covering algorithms, we determined whether the obtained
b(`) indicates the fractality or not. This task was done by
fitting the b(`) curve with a power-law function (i.e., fractal) and an exponential function (i.e., non-fractal) by using optimize.leastsq function in SciPy package of
Python. We used the parameters estimated by fitting the curves
to linearized models as the initial values for the nonlinear
fitting. The key quantity was the ratio between the residual
error of fitting to a power-law function and that to an exponential function, denoted by rfit . If − log10 rfit is postive (i.e.,
rfit < 1), the network was supposed to be fractal. Otherwise,
it was supposed to be non-fractal. This procedure of fitting
and comparison follows that used in [29].
B. Parameter Settings
First of all, we have to decide the parameter values of our
algorithm:  (error tolerance), k (sketch size), and α (exact
neighborhood switch threshold). In principle, the accuracy
of results as well as running time increases with k and α,
and it decreases with . As we discussed in Section V, we

fixed  = 0. To choose k and α, we plotted the average
approximation ratio of our results to the theoretical solutions
for several fractal network models as a function of k and α
in Figure 2. The average approximation ratio is defined by
ρ ≡ hbsketch (`)/btheory (`)i` , where h·i` is the average over
`. To compute bsketch (`), we executed the algorithm for ten
times and took the average of the resulting b(`) over the ten
runs.
In the left panel of Figure 2, we varied 24 ≤ k ≤ 210
while fixing α = 1. The ρ values were affected slightly by
k for the SHM models and tended to decrease with k for
the flower network. On the basis of the results, we decided
to use k = 27 throughout the following experiments. In the
right panel of Figure 2, we varied 2−3 ≤ α ≤ 23 while fixing
k = 27 . The ρ values were almost constant regardless of the
α values for all of the three networks considered. Therefore,
taking into account the running time, we decided to use α = 1
throughout the following experiments. It is worth noting that
Figure 2 also demonstrates the high accuracy and robustness
of our algorithm over a broad range of parameter values.
C. Accuracy and Scalability
Table II summarizes the main results of this paper and shows
the comparison of our algorithm (Sketch) with other preceding
algorithms for fractal and non-fractal network models with
various sizes. We evaluated the performance of algorithms by
two measures. The first was the accuracy given by − log10 rfit
(Section VI-A). If this measure took a positive (negative) value
for a fractal (non-fractal) network, the algorithm correctly
distinguished the fractality of the network. The second was
computation time in seconds.
Discrimination Ability: As we can see in Table II, the sketch
algorithm perfectly distinguishes between the fractal and nonfractal networks as the other algorithms do (except for CBB for
(3, 4, 5)-flower). The proposed algorithm shows its advantage
in computation time: the algorithm is generally faster than
other algorithms and is able to handle large networks that
other algorithms do not terminate. Although MEMB is faster
than Sketch for some relatively small network models, this
result is expected because actual neighborhood sets are not
significantly larger than sketch sizes in these networks. As
a summary, (i) the sketch algorithm correctly detected the
fractality of network models with around ten times smaller
computation time than the fastest previous algorithm. In addition, (ii) the algorithm was able to deal with networks with
millions of nodes with acceptable computation time (within
1 day), whereas other algorithms could not in our machine
environment.
Time and Memory Consumption: The proposed algorithm
is scalable for not only for computation time but also for
memory usage. In Figure 3, computation time (seconds) and
memory usage (KB) of the five algorithms were plotted as
a function of the number of vertices. We use (2, 2, g)-flower
(3 ≤ g ≤ 11) and (2, t)-BA (0 ≤ t ≤ 15) networks as the
example of a fractal and a non-fractal network, respectively.

TABLE II
RUNNING TIME IN SECONDS (Time) AND THE RELATIVE ERROR RATIO OF A POWER - LAW FUNCTION , − log10 rfit (Fit). DNF MEANS THAT IT DID NOT
FINISH IN ONE DAY OR RAN OUT OF MEMORY.

Graph
|V |

Model

|E|

Sketch
Time
Fit

MEMB [26]
Time
Fit

GC [26]
Time
Fit

MVB [24]
Time
Fit

CBB [26]
Time
Fit

5 Networks with ground-truth fractality (“Fit” values are expected to be positive.)
(2,
(2,
(2,
(2,
(2,
(2,
(2,
(2,
(2,
(3,
(3,
(3,
(3,
(2,
(2,
(2,
(3,

2,
2,
2,
2,
3,
3,
3,
4,
4,
3,
3,
4,
4,
0,
0,
0,
0,

4)-flower
7)-flower
10)-flower
11)-flower
6)-flower
7)-flower
8)-flower
6)-flower
7)-flower
6)-flower
7)-flower
5)-flower
7)-flower
6)-SHM
7)-SHM
8)-SHM
6)-SHM

172
10,924
699,052
2,796,204
11,720
58,595
292,970
37,326
223,950
37,326
223,950
14,007
686,287
12,501
62,501
312,501
67,229

256
16,384
1,048,576
4,194,304
15,625
78,125
390,625
46,656
279,936
46,656
279,936
16,807
823,543
12,500
62,500
312,500
67,228

0
15
8,628
62,138
26
286
2,913
138
1,526
148
1,779
34
8,873
33
224
2,728
207

0.8
2.5
3.5
4.0
1.2
1.1
1.0
1.1
1.0
1.1
1.2
0.7
0.8
1.2
1.2
1.1
1.0

0
10
DNF
DNF
14
377
DNF
121
DNF
101
DNF
16
DNF
8
206
DNF
190

1.0
3.4
—
—
1.1
1.0
—
1.0
—
1.3
—
0.9
—
1.1
1.1
—
0.9

0
228
DNF
DNF
146
8,538
DNF
2,422
DNF
10,751
DNF
560
DNF
872
48,116
DNF
21,726

28.7
27.7
—
—
0.1
0.1
—
1.7
—
2.1
—
0.1
—
1.1
1.1
—
0.9

199
DNF
DNF
DNF
DNF
DNF
DNF
DNF
DNF
DNF
DNF
DNF
DNF
32
1,126
DNF
628

1.0
—
—
—
—
—
—
—
—
—
—
—
—
1.1
1.1
—
0.9

0
122
DNF
DNF
5,593
DNF
DNF
2,559
DNF
1,284
61,562
3,380
DNF
325
6,579
DNF
4,623

28.0
27.4
—
—
0.5
—
—
0.6
—
1.4
1.5
-0.4
—
0.7
0.9
—
0.9

-5.4
-6.2
-7.0
-4.8
-6.0
—
-1.8
-1.8
-1.9
-0.4
-0.4
-0.3
-0.3
-0.6
-0.6
-0.6
-0.6
—
—

364
3,610
DNF
DNF
DNF
DNF
DNF
DNF
DNF
126
8,615
25
2,070
54
DNF
DNF
DNF
DNF
DNF

-2.2
-2.6
—
—
—
—
—
—
—
-3.5
-4.9
-2.6
-4.2
-0.5
—
—
—
—
—

21,833
DNF
DNF
1,862
61,953
DNF
3,781
DNF
DNF
7,129
DNF
1,224
DNF
0
404
DNF
DNF
DNF
DNF

-2.6
—
—
-3.3
-4.4
—
-1.9
—
—
-2.5
—
-2.9
—
-0.3
-0.1
—
—
—
—

5 Networks with ground-truth non-fractality (“Fit” values are expected to be negative.)

102

102

Sketch
MEMB
CBB
GC
MVB

100

10-4 2
10
108
107
106
105
104
103 2
10

time (sec)

104

104
103
number of vertices
Sketch
MEMB
CBB
GC
MVB

104
103
number of vertices

100
-2

105

108
466
1,774
20
123
4,195
23
223
1,678
31
390
12
210
0
1
17
377
6,474
36,125

-2.9
-3.8
-4.6
-3.0
-4.7
-6.0
-1.0
-0.8
-0.7
-3.5
-4.9
-2.6
-4.1
-0.9
-2.7
-1.3
-1.5
-1.4
-1.4

Sketch
MEMB
CBB
GC
MVB

10
102 103 104 105 106 107
number of vertices
108
Sketch
MEMB
107
CBB
106
GC
MVB
5
10
104
103 2
10 103 104 105 106 107
number of vertices

105

Fig. 3. Scalability of computation time (top) and memory usage (bottom) for
(2, 2, g)-flower (left) and (2, t)-BA (right) networks.

The symbols corresponding to an algorithm were not shown
if the algorithm did not stop within 24 hours or could not

197
1,641
DNF
16
280
DNF
20
548
DNF
32
1,397
8
580
0
0
76
3,535
DNF
DNF

-2.9
-3.8
—
-2.7
-3.2
—
-0.6
-0.7
—
-3.5
-4.9
-2.6
-4.2
-0.9
-2.0
-1.3
-1.5
—
—

105
104
103
102
101
100 0
10

286
2,999
38,278
44
826
DNF
53
1,598
67,866
433
17,703
97
9,504
0
2
154
12,457
DNF
DNF

theory

0.20
0.15
0.10

CV

59,049
177,147
531,441
16,384
65,536
1,048,576
15,625
78,125
390,625
31,104
186,624
16,384
131,072
497
3,997
31,997
255,997
2,047,997
8,191,997

104

10-2

memory (KB)

29,526
88,575
265,722
10,924
43,692
699,052
11,720
58,595
292,970
24,885
149,301
14,045
112,349
250
2,000
16,000
128,000
1,024,000
4,096,000

b ( `)

2, 10)-flower
2, 11)-flower
2, 12)-flower
3, 7)-flower
3, 8)-flower
3, 9)-flower
4, 6)-flower
4, 7)-flower
4, 8)-flower
1, 6)-SHM
1, 7)-SHM
1, 5)-SHM
1, 6)-SHM
1)-BA
4)-BA
7)-BA
10)-BA
13)-BA
15)-BA

memory (KB)

time (sec)

(1,
(1,
(1,
(1,
(1,
(1,
(1,
(1,
(1,
(2,
(2,
(3,
(3,
(2,
(2,
(2,
(2,
(2,
(2,

0.05
101

102

`

103

104

0.00
100

101

102

`

103

104

Fig. 4. Results of different runs for (3, 3, 7)-flower. (Left) b(`) and (Right)
CV of them as a function of `.

execute owing to memory shortage. The performance of the
proposed algorithm is comparable to or worse than some other
algorithms when the network is relatively small (i.e., n < 104 ).
However, the algorithm is orders of magnitude faster than
other algorithms for large networks. Also, it achieves such
high a high speed with incomparably smaller memory usage
than MEMB, the second fastest algorithm.
Robustness over Randomness: The sketch algorithm accurately recovers b(`) of theoretical prediction for fractal network

b ( `)

106
10

5

10

4

∝
∝

VII. C ONCLUSIONS

`−3.4
e−0.3 `

103
102
101
100 0
10

101

`

102

Fig. 5. Results for a real web graph.

models, and the results are robust over different execution runs.
The left panel of Figure 4 shows b(`) of ten different runs
of the proposed algorithm on (3, 3, 7)-flower. The b(`) values
follow well the theoretical solution, which is indicated by the
solid line. As we can clearly observe, the fluctuation in the b(`)
values due to the randomness is very small. The consistency
over different runs is captured by the CV of b(`) (i.e., the ratio
of the standard deviation of b(`) to its average over ten runs)
as a function of ` (right panel of Figure 4). The CV values
tend to increase with `. This tendency can be explained by
the following two factors. First, the b(`) value takes a positive
integer value and monotonically decreases with ` by definition.
Thus, even a change of ±1 in b(`) might cause a large CV
value if ` is large. Second, our algorithm intrinsically fluctuates
more when ` is larger. This could be because the sizes of the
solutions become smaller for larger `, and hence the algorithm
gets a little more sensitive to estimation errors. Nevertheless,
it should be noted that the variance of our algorithm was
considerably small even for large ` (i.e., CV ∼ 0.19 at most).
This magnitude of variance would have little impact on the
estimation of fractality.
D. Application to Real Large Network
In closing this section, we applied the sketch algorithm to a
large-scale real graph to show the scalability of the proposed
algorithm with an empirical instance. The results also gave
us some insight on the fractality of large-scale real-world
networks, which is beyond the reach of previous algorithms.
As a representative instance of a real-world large graph, we
considered the in-2004 network [6], [4], which is a crawled
web graph of 1, 382, 908 vertices and 16, 917, 053 edges. We
discarded the direction of the edges (i.e., hyperlinks) to make
the network undirected. The algorithm took 11.7 hours in total.
The resulting b(`) of the sketch algorithm and the fitting
curves are shown in Figure 5. We omitted the three points
with the smallest ` values from the fitting because empirical
networks would not show a perfect fractality, contrary to welldesigned network models. A large part of the points fall on the
line of the fitted power-law function, and indeed, our fractality
decision procedure yielded − log10 rfit = 0.79, which suggests
the fractality of the in-2004 network. It is worth mentioning
that the fractality of this network was unveiled for the first time
for the sake of our algorithm.

Fractality is an interesting property that appears in some
classes of real networks. In the present study, we designed a
new box-covering algorithm, which is useful for analyzing the
fractality of large-scale networks. In theory, we have shown
desirable guarantees on scalability and solution quality. In the
experiments, we confirmed that the algorithm’s outputs are
sufficiently accurate and that it can handle large networks with
millions of vertices and edges. We hope that our method enables further exploration of graph fractality and its applications
such as graph coarsening.
Repeatability: Our implementation of the proposed and previous box-cover algorithms is available at http://git.io/fractality.
It also contains the generators of the synthetic network models,
and thus the results in this paper can be perfectly replicated.
We hope that our public software will enable further exploration of graph fractality and its applications.
Acknowledgment: This work was supported by JSPS KAKENHI (No. 15H06828), JST, ERATO, Kawarabayashi Large
Graph Project, and JST, PRESTO. Web graph data was downloaded from http://law.di.unimi.it/datasets.php. T.T. thanks to
K. Takemoto for valuable discussions.
R EFERENCES
[1] T. Akiba, Y. Iwata, and Y. Yoshida. Fast exact shortest-path distance
queries on large networks by pruned landmark labeling. In SIGMOD,
page 349, 2013.
[2] T. Akiba, K. Nakamura, and T. Takaguchi. Fractality of massive graphs:
Scalable analysis with sketch-based box-covering algorithm. In ICDM,
2016. to appear.
[3] A.-L. Barabási and R. Albert. Emergence of Scaling in Random
Networks. Science, 286(5439):509–512, 1999.
[4] P. Boldi, M. Rosa, M. Santini, and S. Vigna. Layered label propagation: A multiresolution coordinate-free ordering for compressing social
networks. In WWW, pages 587–596, 2011.
[5] P. Boldi, M. Rosa, and S. Vigna. HyperANF: Approximating the
neighbourhood function of very large graphs on a budget. In WWW,
pages 625–634, 2011.
[6] P. Boldi and S. Vigna. The WebGraph framework I: Compression
techniques. In WWW, pages 595–601, 2004.
[7] G. Caldarelli. Scale-Free Networks. Oxford University Press, 2007.
[8] S. Carmi, S. Havlin, S. Kirkpatrick, Y. Shavitt, and E. Shir. A model
of Internet topology using k-shell decomposition. Proc. Natl. Acad. Sci.
USA, 104(27):11150–11154, 2007.
[9] E. Cohen. Size-estimation framework with applications to transitive
closure and reachability. J. Comput. Syst. Sci., 55(3):441–453, 1997.
[10] E. Cohen. All-distances sketches, revisited: HIP estimators for massive
graphs analysis. IEEE TKDE, 27(9):2320–2334, 2015.
[11] E. Cohen, D. Delling, T. Pajor, and R. F. Werneck. Sketch-based
influence maximization and computation: Scaling up with guarantees.
In CIKM, pages 629–638, 2014.
[12] K. Falconer. Fractal Geometry: Mathematical Foundations and Applications, Second Edition. Wiley-Blackwell, 2003.
[13] L. K. Gallos, H. A. Makse, and M. Sigman. A small world of weak ties
provides optimal global integration of self-similar modules in functional
brain networks. Proc. Natl. Acad. Sci. USA, 109(8):2825–2830, 2012.
[14] L. K. Gallos, C. Song, and H. A. Makse. A review of fractality and
self-similarity in complex networks. Physica A, 386:686–691, 2007.
[15] K.-I. Goh, G. Salvi, B. Kahng, and D. Kim. Skeleton and Fractal Scaling
in Complex Networks. Phys. Rev. Lett., 96(1):018701, 2006.
[16] T. Hasegawa and K. Nemoto. Hierarchical scale-free network is fragile
against random failure. Phys. Rev. E, 88(6):062807, 2013.
[17] D. Kempe, J. Kleinberg, and É. Tardos. Maximizing the spread of
influence through a social network. In KDD, pages 137–146, 2003.

[18] J. S. Kim, K.-I. Goh, B. Kahng, and D. Kim. Fractality and selfsimilarity in scale-free networks. New J. Phys., 9(6):177, 2007.
[19] J. Kleinberg. The small-world phenomenon. In STOC, pages 163–170,
2000.
[20] J. Leskovec, J. Kleinberg, and C. Faloutsos. Graphs over Time :
Densification Laws, Shrinking Diameters and Possible Explanations. In
KDD, pages 177–187, 2005.
[21] M. E. J. Newman. Networks: an Introduction. Oxford University Press,
2010.
[22] C. R. Palmer, P. B. Gibbons, and C. Faloutsos. ANF: A fast and scalable
tool for data mining in massive graphs. In KDD, pages 81–90, 2002.
[23] H. D. Rozenfeld, S. Havlin, and D. Ben-Avraham. Fractal and transfractal recursive scale-free nets. New J. Phys., 9(6):175, 2007.
[24] C. M. Schneider, T. A. Kesselring, J. S. Andrade, and H. J. Herrmann.
Box-covering algorithm for fractal dimension of complex networks.
Phys. Rev. E, 86(1):016707, 2012.
[25] M. Á. Serrano, D. Krioukov, and M. Boguñá. Percolation in Self-Similar
Networks. Phys. Rev. Lett., 106(4):048701, 2011.
[26] C. Song, L. K. Gallos, S. Havlin, and H. A. Makse. How to calculate
the fractal dimension of a complex network: the box covering algorithm.
J. Stat. Mech., 2007(03):P03006, 2007.
[27] C. Song, S. Havlin, and H. A. Makse. Self-similarity of complex
networks. Nature, 433(7024):392–395, 2005.
[28] C. Song, S. Havlin, and H. A. Makse. Origins of fractality in the growth
of complex networks. Nat. Phys., 2(4):275–281, 2006.
[29] K. Takemoto. Metabolic networks are almost nonfractal: A comprehensive evaluation. Phys. Rev. E, 90(2):022802, 2014.

