# General expression for the component size distribution in infinite
configuration networks

Ivan Kryven<EMAIL_ADDRESS>University of Amsterdam, PO box 94214, 1090 GE,
Amsterdam, The Netherlands

###### Abstract

In the infinite configuration network the links between nodes are assigned
randomly with the only restriction that the degree distribution has to match a
predefined function. This work presents a simple equation that gives for an
arbitrary degree distribution the corresponding size distribution of connected
components. This equation is suitable for fast and stable numerical
computations up to the machine precision. The analytical analysis reveals that
the asymptote of the component size distribution is completely defined by only
a few parameters of the degree distribution: the first three moments, scale
and exponent (if applicable). When the degree distribution features a heavy
tail, multiple asymptotic modes are observed in the component size
distribution that, in turn, may or may not feature a heavy tail.

networks; degree distribution; connected components; heavy tail; power law;
giant component

###### pacs:

64.60.aq, 02.10.Ox, 89.75.Da

##  I Introduction

Random graphs provide models for complex networks, and in many cases, real-
world networks has been accurately described by such modelsBarabási and Albert
(1999); Newman  _et al._ (2002); Vázquez (2003); Kryven  _et al._ (2016).
Within the scope of random graph models one finds: Erdős-Rényi model,
Barabási-Albert modelBarabási and Albert (1999), node copying modelBhat  _et
al._ (2016), small world networkNewman and Watts (1999), configuration
networkMolloy and Reed (1998); Newman (2010) and many others. In the
configuration network $N$ nodes are assigned pre-defined degrees. The edges
connecting these nodes are then considered to be random, and every distinct
configuration of edges that satisfies the given degree sequence is treated as
a new instance of the network in the sense of random graphs. Interesting
properties emerge when the number of nodes, $N,$ approaches infinity, or at
the so-called thermodynamic limitChung  _et al._ (2003); Kryven (2016). In
this case, the infinite degree sequence, which provides the only input
information for the model, is equivalent to the frequency distribution of
degrees, $u(k),\;k=1,2,\dots,$ i.e. the probability that a randomly chosen
node has degree $k$.

Component-size distribution, $w(n),$ denotes probability that a randomly
chosen node is part of a connected component of finite size $n.$ Connected
components in the infinite configuration network can be of finite or infinite
size. Molloy and ReedMolloy and Reed (1995) showed that if an infinite
component exists then it is the only infinite component with probability 1.
Hence, the infinite component is referred to as the giant component.

Depending upon a specific context behind the network, the component size
distribution may summarise an important feature of the modelled system. In
polymer chemistry, for example, the infinite configuration network is used as
a toy model for hyper-branched and cross-linked polymers. In this context, the
component size distribution predicts viscoelastic properties of the material
while the emergence of the giant component is interpreted as a phase
transition from liquid to solid state of the soft matterKryven  _et al._
(2016); Kryven and Iedema (2014). Since connected components are closely
related to clusters in bond percolation processes, the distribution of
component sizes can be used to model outbreaks for SIR epidemiological
processesNewman  _et al._ (2002). In linguistics, component size distribution
of the sentence similarity graph is an important tool when studying structure
of natural languages Biemann and van den Bosch (2011). This brief list of
application cases is far from being exhaustive. Despite the vast applications,
the empirical component size distribution is hard to measure precisely unless
the whole topology of the network is known. On another hand, empirical
observations on the degree distribution, $u(k),$ are much easier to perform.

Ref.Molloy and Reed (1995) provides an elegant criterion that connects moments
of $u(k)$ to the fact that the network contains the giant component. A
somewhat deeper question further in this direction reads: providing $u(k)$ is
given, what is the component size distribution, $w(n)?$ In Ref.Newman  _et
al._ (2001) Newman et al. showed that the component size distribution can be
recovered by a numerical algorithm that involves solving a fixed point problem
followed by an inversion of a generation function. Such algorithm demonstrates
that indeed $u(k)$ and $w(n)$ can be put into a correspondence, however, it
becomes computationally infeasible for large values of $n$. This numerical
issues aries due to ill-posedness of the numerical generating function
inversion.

On another hand, the component size distribution has been analytically
resolved only for a limited number of partial cases of $u(k)$Newman (2007).
Within the scope of analytically solvable cases, only the Yule-Simon degree
distribution features a heavy tail, that is to say it decays proportionally to
an algebraic function, $n^{-\beta},\;\beta>0$ at large $n$. At the same time,
the heavy-tailed (or scale-free) distributions are commonly observed in the
empirical data collected from many real-world networksYook  _et al._ (2002);
Ravasz and Barabási (2003); Eguiluz  _et al._ (2005); Fu  _et al._ (2008).
Empirically observed exponents vary in a broad range. Some studies report
degree exponents that are as small as $\beta=0.81$ in the case of the Internet
topologyAdamic and Huberman (2002) and $\beta=1$ in social networks Timár  _et
al._ (2016). On the opposite side of this spectrum, one finds exponent
$\beta=5$ in the generalisation of preferential attachment model Dorogovtsev
_et al._ (2000).

The only asymptotic analysis available for component size distribution in the
configuration network states that for large $n,$ $w(n)$ is either proportional
to $n^{-3/2}$ or to $e^{-Cn},$ where $C>0$ is a constant Newman  _et al._
(2001). The current paper uncovers new asymptotic modes for $w(n)$ that emerge
only when the degree distribution features a heavy tail. The paper shows that
for an arbitrary $u(k),$ $w(n)$ can be expressed as a finite sum. In practice,
this sum can be stably computed up to the machine precision in the cost of
$O(n^{2})$ multiplicative operations. Finally, the paper discusses how a
finite cutoff introduced in the degree distribution reflects on the
distribution of component sizes.

##  II Component size distribution by Lagrange inversion

It has been noticed that all components in the infinite configuration model
are locally tree-like. Using this fact as a departure point, Newman et
al.Newman  _et al._ (2001) showed that the degree distribution can be put into
a correspondence to the component size distribution by applying the
generation-function (GF) formalism. Here, by a GF of
$u(k),\;\sum_{k=0}^{\infty}u(k)=1,$ we refer to the series,

$U(x)=\sum_{k=0}^{\infty}u(k)x^{k},\;x\in\mathbb{C},\;|x|\leq 1.$ (1)

According to the approach presented in Ref. Newman  _et al._ (2001), the
generating function for the component size distribution, $W(x)$ is found as a
solution of the following system of functional equations,

$\displaystyle W(x)=xU[W_{1}(x)],$ (2) $\displaystyle
W_{1}(x)=xU_{1}[W_{1}(x)].$ (2′)

where $U(x)$ is the GF of $u(k)$, and $U_{1}(x)$ is the GF for the excess
degree distribution

$u_{1}(k)=\frac{k+1}{\mu_{1}}u(k+1),$ (3)

where $\mu_{1}=\sum_{k=1}^{\infty}ku(k).$ Similarly to combinatorial tree-
counting problems, Eq. (2) can be solved by applying the Lagrange inversion
formulaBergeron  _et al._ (1998). The original formulation of the Lagrange
inversion principle is as follows. Suppose, $A(X),R(x),$ are such formal power
series that $A(x)=xR[A(x)]$ then for an arbitrary formal power series $F(x),$
the coefficient of power series $F[A(x)]$ at $x^{n}$ reads as,

$[x^{n}]F[A(x)]=\frac{1}{n}[t^{n-1}]F^{\prime}(t)R^{n}(t),\;n>0.$ (4)

Here $[t^{n-1}],$ as being the inverse operation to the GF transform (1),
refers to the coefficient at $t^{n-1}$ of the corresponding power series. By
substituting $A(x)=W_{1}(x)$, $F(x)=U(x)$ and applying Eq. (2) one transforms
the left hand side of (4),
$[x^{n}]F[A(x)]=[x^{n}]U[W_{1}(x)]=[x^{n-1}]W(x)=w(n-1).$ Further on, the
right hand side of (4) is transformed by substituting $R(x)=U_{1}(x)$ and
realising that, according to the definition (3),
$U^{{}^{\prime}}(x)=\mu_{1}U_{1}(x).$ Now, we are ready to write an expression
for $w(n),$ even though we have no explicit expression for generating function
$W(x)$ itself,

$w(n)=\frac{1}{n-1}[t^{n-2}]U^{{}^{\prime}}(x)U_{1}(x)^{n-1}=\\\
\frac{\mu_{1}}{n-1}[t^{n-2}]U_{1}(x)^{n},\;n>1.$ (5)

A similar equation was also derived in Ref.Newman (2007) by means of different
reasoning. In principle, Eq. (5) provides enough information to analytically
recover the component size distribution for a few special cases of the degree
distributionNewman (2007). In practice, however, the main difficulty when
applying Eq. (5) is that the equation employs the inverse GF transform,
$[t^{n}]$, which limits the choices one has when searching for an exact
solution or performing numerical computations. With this in mind, one may
rewrite (5) so that the new expression does not involve the GF concept at all.
It turns out that the only reason why Eq. (5) utilises the GF formalism is
that it provides means for _convolution power_.

The convolution of two distributions, $f(k)*g(k),\;k>0$ is defined as a binary
multiplicative operation,

$f(k)*g(k)=\sum\limits_{i+j=k}f(i)g(j),$

where the summation is performed over all non-negative ordered couples $i,j$
that sum up to $k$. This sum contains exactly $k+1$ of such couples. In this
paper, the order of operations is chosen in such a way that the point-wise
multiplication precedes convolution, for instance, $f(k)*kg(k)=f(k)*[kg(k)]$.
The convolution can be inductively extended to the $n$-fold convolution, or
_the convolution power_ ,

$f(k)^{*n}=f(k)^{*n-1}*f(k),$ (6)

where $f(k)^{*0}\equiv 1$ by the definition. It can be shown that the
convolution power can be expanded into a sum of products,

$f(k)^{*n}=\sum\limits_{\stackrel{{\scriptstyle
k_{1}+\dots+k_{n}=k}}{{k_{i}\geq 0}}}\prod\limits_{i=1}^{n}f(k_{i}).$ (7)

The convolution has a peculiar property in respect to the GF transform. If
$F(x),G(x),U(x)$ are GFs for $f(k),g(k),$ and $u(k)=f(k)*g(k)$ then
$U(x)=F(x)G(x)$. Furthermore, if $U(x)$ is GF for $u(k)$ then $U(x)^{n}$
generates $u(k)^{*n}.$ By exploiting this relation one immediately reduces Eq.
(5) to,

$w(n)=\begin{cases}\frac{\mu_{1}}{n-1}u_{1}^{*n}(n-2),&n>1,\\\
u(0)&n=1.\end{cases}$ (8)

Here, the value of $w(0)$ is derived directly from the formulation of the
problem: nodes with degree zero are also components of size one. This simple
equation is ready to be used: by combining (8) and the definitions (3),(7) one
may directly expresses the values of the component size distribution in terms
of $u(k):$ for $n>1,$

$w(n)=\frac{[ku(k)]^{*n}(2n-2)}{(n-1)\mu_{1}^{n-1}},\;n>1.$ (9)

For example, first five values of $w(n)$ read as,

$\displaystyle w(1)=$ $\displaystyle u(0),$ $\displaystyle w(2)=$
$\displaystyle\frac{1}{\mu_{1}}u(1)^{2},$ $\displaystyle w(3)=$
$\displaystyle\frac{3}{\mu_{1}^{2}}u(1)^{2}u(2),$ $\displaystyle w(4)=$
$\displaystyle\frac{4}{\mu_{1}^{3}}u(1)^{2}[2u(2)^{2}+u(1)u(3)],$
$\displaystyle w(5)=$
$\displaystyle\frac{5}{\mu_{1}^{4}}u(1)^{2}[4u(2)^{3}+6u(1)u(2)u(3)+u(1)^{2}u(4)].$

The number of terms in this expansion increases rapidly with $n$. That said,
the formula (8) can be easily readjusted for numerical computations. Namely,
one can use Fast Fourier Transform (FFT) to compute the convolution powers,
$u_{1}^{*n}(k)=\mathcal{F}^{-1}[\mathcal{F}[u_{1}(k)]^{n}].$ In this case,
$O(N^{2})$ multiplicative operations is sufficient to compute all values of
$w(n),n\leq N$. Alternatively, if $w(n-1)$ is known $w(n)$ can be found in the
cost of $O(n\log n).$ Besides FFT, there are algorithms that are specifically
designed for fast approximation of convolution powers, such as projection onto
basis functions that are invariant under convolutionKryven and Iedema (2013).

Analytic formulas for convolution powers (sometimes also referred to as
_compositas_ Kruchinin (2015)), were covered by literature for many elementary
functionsMa and Liu (2004); Ma and King (2002). Convolution powers of
$u_{1}(k)$ can also be found analytically by applying discrete functional
transforms, for instance Z-transform and discrete Fourier transform. A few
examples of such results are given in Table 1. Focusing on one of them, the
first curve in Figure 1 demonstrates that both analytical and numerical
results for the exponential degree distribution coincide.

Degree Distribution, | Component Size Distribution
---|---
Exponential distribution
$Ce^{-\lambda k}$ | $\frac{(1-e^{-\lambda})^{2n-1}}{e^{\lambda(n-1)}}\frac{\Gamma(3n-2)}{\Gamma(n)\Gamma(2n)}$
Geometric distribution
$(1-p)^{k-1}p$ | $(1-p)^{n-2}p^{2n-1}\frac{\Gamma(3n-2)}{\Gamma(n)\Gamma(2n)}$
Binomial distribution
$\binom{k_{\text{max}}}{k}(1-c)^{k_{\text{max}}-k}c^{k}$ | $\frac{1}{n-1}\binom{nk_{\text{max}}-n}{n-2}\times$
$\;\;\;\;\;\;\;(1-c)^{nk_{\text{max}}-2n-2}c^{n-2}$
Table 1: Exact expressions for component size distributions in configuration
network as evaluated with Eq. (8) via Z-Transform. Figure 1:  Examples of
component size distributions (solid lines) that feature fast (_1_) and slow
(_2_) convergence to their asymptotes (dashed lines). Both asymptotes are
covered by Case A, Table 2. _1_) $u(k)=Ce^{-1.05k},$ all three: the analytical
expression (see Table 1), numerical values (according to Eq. (8)) and the
asymptote practically coincide. _2_) $u(k)$ is non-zero in three points
$u(1)=0.97,\;u(2)=0.015,\;u(10)=0.015;$ the component size distribution
features oscillations before it converges to the asymptote.

##  III Asymptotic analysis

The format of Eq. (8), naturally suggests a straightforward way to perform an
asymptotic analysis for $n\to\infty.$ One may view $u(k)$ as a probability
mass function PMF (or alternatively discrete probability density function) of
some discrete random variables $k_{i}.$ Recall the following property of
convolution powers: if i.i.d. random variables $k_{i}$ have PMF $u_{1}(k)$
then $u_{1}^{*n}(k)$ gives the PMF for the sum $k_{1}+k_{2}+\dots+k_{n}$. The
central limit theorem (CLT) gives an estimate for this sum as $n\to\infty,$
and the idea is now to obtain the asymptotes of $w(n)$ by applying CLT to the
definition (8).

###  III.1 Light-tailed degree distributions

First, let us assume that distribution $u(k)$ decays faster than
algebraically, that is

$u(k)=o(k^{-\beta}),\;\beta>2,\;k\to\infty,$ (10)

which is also equivalent to $u_{1}(k)=o(k^{-\beta+1})$. Then according to CLT,
$u_{1}^{*n}(k)$ approaches the normal distribution,
$u_{1}^{*n}(k)\stackrel{{\scriptstyle
d}}{{\to}}(\sqrt{n}\sigma)^{-1}\mathcal{N}(\frac{k-nM}{\sqrt{n}\sigma},0,1),$
when $n\to\infty,$ where $M=\sum_{k=1}^{\infty}ku_{1}(k)$ and
$\sigma^{2}=\sum_{k=1}^{\infty}k(k-M)^{2}u_{1}(k)<\infty$ denote the mean
value and variance of $u_{1}(k)$. The normal distribution can now replace
$u_{1}^{*n}(k)$ in (8), which yields the asymptote for the component size
distribution,

$w(n)\sim\frac{\mu_{1}e^{-\frac{(n(1-M)-2)^{2}}{2n\sigma^{2}}}}{(n-1)\sqrt{2\pi
n\sigma^{2}}},\text{ as }n\to\infty.$ (11)

Quantities $M,\sigma^{2}$ are directly expressible in terms of moments of
degree distribution $u(k),$

$\displaystyle M=$
$\displaystyle\sum_{k=1}^{\infty}ku_{1}(k)=\frac{1}{\mu_{1}}\sum_{k=1}^{\infty}(k^{2}-k)u(k)=\frac{\mu_{2}-\mu_{1}}{\mu_{1}},$
(12) $\displaystyle\sigma^{2}=$
$\displaystyle\sum_{k=1}^{\infty}k(k-M)^{2}u_{1}(k)=$
$\displaystyle\frac{1}{\mu_{1}}\sum_{k=0}^{\infty}k(k-M-1)^{2}u(k)=\frac{\mu_{3}\mu_{1}-\mu_{2}^{2}}{\mu_{1}^{2}},$

where

$\mu_{i}=\mathbb{E}[k^{i}]=\sum_{k=1}^{\infty}k^{i}u(k),i=1,2,\dots$

Finally, substituting the expressions (12) into (11) and gives the final
version of the asymptote,

$w(n)\sim\frac{\mu_{1}^{2}n^{-3/2}e^{-\frac{(\mu_{2}-2\mu_{1})^{2}}{2(\mu_{1}\mu_{3}-\mu_{2}^{2})}n}}{\sqrt{2\pi(\mu_{1}\mu_{3}-\mu_{2}^{2})}},\text{
as }n\to\infty.$ (13)

Two examples of component size distributions that converge with various rates
to their asymptotes are given in Figure 1. Peculiarly, the only information on
$u(k)$ that is contained in the asymptote definition (13) is the first three
moments $\mu_{1},\mu_{2},\mu_{3}$. Furthermore, depending upon the value of
$\theta=\mu_{2}-2\mu_{1},$ the asymptotic expression (13) switches between the
_two modes_ : it either decays exponentially as $O(e^{-An}),$ when $\theta\neq
0$, or it decays as an algebraic function, $O(n^{-3/2}),$ when $\theta=0$ (see
also Table 2, Case A). The last equality is the well-known giant component
criterion,

$\mu_{2}-2\mu_{1}=0.$ (14)

The criterion (14) was obtained by Molloy and ReedMolloy and Reed (1995) by
means of a different reasoning. In Ref Molloy and Reed (1995), the authors
prove that $\theta>0$ implies existence of the giant component in the
configuration network, whereas $\theta<0$ implies non-existence of this
component. In Ref. Newman  _et al._ (2001), it was hypothesised that the
$-3/2$ exponent is universal and must hold for all degree distributions at the
critical point $\theta=0$. We will see now that when the condition (10) fails
to hold, distinct from $-3/2$ exponents may also appear in the asymptotic of
$w(n)$.

###  III.2 Heavy-tailed degree distributions

Suppose that, on the contrary to the condition (10), degree distribution
$u(k)$ features a heavy tail,

$u(k)\sim sk^{-\beta},\;\beta>2,\;k\to\infty,$ (15)

which is equivalent to $u_{1}(k)\sim sk^{-\alpha-1},\;\alpha={\beta-2}>0$,
$k\to\infty.$ It turns out that exponent $\alpha$ and the scale $s,$ together
with the moments $\mu_{1},\mu_{2},\mu_{3}$ provide enough information to
generalise the asymptote (13) for the case of heavy-tailed degree
distributions. Suppose $0<\alpha\leq 2.$ In terms of $u(k)$ moments this
condition casts out as $\mu_{3}=\infty.$ As follows from Gnedenko and
Kolmogorov’s generalisation of CLTGnedenko and Kolmogorov (1954) the mass
density distribution for $u_{1}^{*n}(k)$ approaches the stable law,

$u_{1}^{*n}(k)\stackrel{{\scriptstyle
d}}{{\to}}\frac{1}{\gamma(n)}G^{A}\left(\frac{k-\mu(n)}{\gamma(n)},\alpha,1\right),\;n\to\infty.$
(16)

Here, we use the notation of Uchaikin & Zolotarev Uchaikin and Zolotarev
(1999) which includes: _exponent parameter_ $\alpha,$ the _location parameter_

$\mu(n)=\begin{cases}n\frac{\mu_{2}-\mu_{1}}{\mu_{1}},&\alpha>1,\\\ sn\ln
n,&\alpha=1,\\\ 0,&0<\alpha<1,\end{cases}$ (17)

and the _scale parameter_

$\gamma(n)=\begin{cases}\sqrt{sn\ln n},&\alpha=2,\\\ \sqrt{\pi
s}[2\Gamma(\alpha)\sin\frac{\alpha\pi}{2}]^{-1/\alpha}n^{1/\alpha},&\alpha\in(0,1)\cup(1,2)\\\
\frac{\pi ns}{2},&\alpha=1.\end{cases}$ (18)

No general analytical expression is known for $G^{A}(x,\alpha,1),$ and the
stable law is defined via its Fourier transform

$\mathcal{F}[{G^{A}(x,\alpha,1)}]=\begin{cases}e^{-x^{\alpha}-i\tan\frac{\pi\alpha}{2}},&\alpha\in(0,1)\cup(1,2],\\\
e^{-x^{\alpha}+i\frac{2\alpha}{\pi}},&\alpha=1.\end{cases}$ (19)

Consider the case when $1<\alpha<2.$ According to (16) the point in which the
stable law is evaluated, $x(n)=\frac{n-\mu(n)}{\gamma(n)}$, approaches
positive or negative infinities depending upon the sign of
$\theta=\mu_{2}-2\mu_{1}.$ Indeed, as $n\to\infty,$

$x(n)\to\begin{cases}+\infty,&\theta<0,\\\ 0,&\theta=0,\\\
-\infty,&\theta>0.\end{cases}$ (20)

For these values of $\alpha,$ function $G^{A}(x,\alpha,1)$ is non-zero on
$(-\infty,+\infty).$ If $x(n)\to\infty,$ the function features an algebraic
decay, whereas if $x(n)\to-\infty$ the decay is exponential. Therefore, the
limiting value switching that takes place in (20) may reflect on the
asymptotic behaviour of $u_{1}^{*n}(n).$ To give a precise answer one has to
consider series expansions of $G^{A}(x,\alpha,1)$ around the points of
interest, $x\in\\{-\infty,0,+\infty\\}$. We use here the leading terms of
these seriesUchaikin and Zolotarev (1999),

$G^{A}(x,\alpha,1)=\\\
\begin{cases}\frac{\Gamma\left(1+\frac{1}{\alpha}\right)\sin\frac{\pi}{\alpha}}{\pi}+O(x),&x\to
0\\\
\frac{\Gamma(\alpha+1)x^{-\alpha-1}}{\Gamma(2-\alpha)\Gamma(\alpha-1)}+O(x^{-2\alpha-1}),&x\to\infty,\\\
\frac{e^{-(\alpha-1)\left(\frac{x}{\alpha}\right)^{\frac{\alpha}{\alpha-1}}}\left(\frac{x}{\alpha}\right)^{\frac{1}{2}\left(\frac{1}{\alpha-1}-1\right)}}{\sqrt{2\pi\alpha(\alpha-1)}}[1+O(x^{-\frac{\alpha-1}{\alpha}})],&x\to-\infty,\end{cases}$
(21)

By replacing the expression for the limiting distribution (16) with the
leading terms given in (21) one obtains the asymptotes for (8). This time, the
asymptote has _three modes_ : depending upon the value of $\theta,$ it either
features a heavy tail with exponent $-\alpha-1,$ a heavy tail with exponent
$-\frac{1}{\alpha}-1,$ or an exponential decay, as shown in Table 2, Case D. A
few examples of such asymptotic modes for a heavy-tailed degree distribution

$u(k)=\begin{cases}C&k=1,\\\ s(\beta-2)k^{-\beta}&k>1\end{cases}$ (22)

are given in Fig. 2. The degree distribution (22) is defined by two
parameters: exponent $\beta$ and scale $s$; whereas the constant $C$ is such
that the total probability is normalised, $\sum_{k}u(k)=1$.

Figure 2:  Component size distributions (_solid lines_) and their asymptotes
(_dashed lines_) as obtained for degree distributions with exponent
$\beta=3.5$ ($\alpha=1.5$) and various values of scale parameter. Three
distinct asymptotic modes are illustrated: _1_) $s=0.066:$ $\theta<0,$
$\eta=\frac{5}{2}$; _2_) $s=0.644:$ $\theta\approx 0$, $\eta=\frac{5}{3}$;
_3_) $s=0.8:$ $\theta>0$;

When $\alpha=2$, the behaviour of $\frac{n-\mu(n)}{\gamma(n)},$ $n\to\infty$
is identical to (20), but the expression for $\gamma(n)$ is different and the
series expansions (21) lead to somewhat different asymptotes, see Table 2,
Case C.

According to the definition (17), the location parameter vanishes,
$\mu(n)\equiv 0,$ when $\alpha<1$. In this case, $x(n)=\frac{n}{\gamma(n)}\to
0$ as $n\to\infty,$ and only _one_ asymptotic mode is possible for $w(n)$.
Stable law $G^{A}(x,\alpha,1)$ is supported on $(0,\infty),$ and we make use
of the series expansion around $x\to 0^{+},$

$G^{A}(x,\alpha,1)=\\\
\frac{e^{-(1-\alpha)\left(\frac{\alpha}{x}\right)^{\frac{\alpha}{1-\alpha}}}\left(\frac{\alpha}{x}\right)^{\frac{1}{2}\left(1+\frac{1}{1-\alpha}\right)}}{\sqrt{2\pi\alpha(1-\alpha)}}[1+O(x^{\frac{1-\alpha}{\alpha}})],\;x\to
0^{+},$ (23)

which when plugged in (16) yields faster then algebraic decay of the component
size distribution, see Table 2, Case F. Due to the parametrisation scheme for
the stable law, the point $\alpha=1$ needs to be considered separately. In
this case, $\frac{n-\mu(n)}{\gamma(n)}\to\infty$ when $n\to\infty,$ and we
utilise the leading term of the series expansion,

$G^{A}(x,\alpha,1)=\frac{1}{\sqrt{2\pi}}e^{\frac{x-1}{2}-e^{x-1}}[1+O(e^{1-x})],\;x\to\infty,$
(24)

which admits one sub-algebraic asymptotic mode for $w(n)$ as shown in Table 2,
Case E. This case is special in that the stable law $G^{A}(x,\alpha,1)$ is
supported on $x\in(-\infty,\infty),$ but asymptotically,
$\frac{1}{\gamma(n)}G^{A}\left(\frac{n-\mu(n)}{\gamma(n)},\alpha,1\right)$
always tends to $-\infty$ for large $n.$ At the same time, if for small $n$
the point $x(n)=\frac{n-\mu(n)}{\gamma(n)}$ stays on the positive half-axis
where (24) does not provide correct description for $G^{A}(x,\alpha,1),$ the
convergence to the asymptote will be slow. In other words, there is an
intermediate asymptote that the component size distribution can be
approximated with, before it eventually switches to Eq. (24). This switching
point is given by such $n_{0}$ that $x(n)$ changes the sign from
${}^{\prime}+^{\prime}$ to ${}^{\prime}-^{\prime},$ i.e. when $n$ becomes
greater then $n_{0}.$ By solving $x(n_{0})=0,$ one obtains
$n_{0}=e^{\frac{1}{s}}$, which means that in principle, the switching between
the intermediate and the final asymptotes may be indefinitely postponed if $s$
is small enough. The intermediate asymptote itself is deduced from the leading
term of the stable law expansion at $\infty,$ that is
$G^{A}(x,1,1)=\frac{2}{\pi}x^{-2}+O(x^{-3}).$ After the substitutions one
obtains,

$w(n)\simeq\frac{\mu_{1}s}{(s\log
n-1)^{2}}n^{-2},\;\frac{1}{s}>>0,\;n<e^{\frac{1}{s}},\;\alpha=1.$

As illustrated in Fig. 3, similar considerations are also valid for the case
$0<\alpha<1,$ where

$w(n)\simeq\mu_{1}s\frac{\Gamma(\alpha+1)}{\Gamma(\alpha)}n^{-\alpha-1},\;\frac{1}{s}>>0,\;0<\alpha<1.$

When occurs, such switching has a practical importance when dealing with
empirically observed component size data. Indeed, it may happen that one
observes only the intermediate asymptote and not the final one due a small
number of samples at the tail of the component size distribution. For
instance, the second curve in Fig. 3 does feature an exponential decay at
infinity, but if one limits the data points to $n<10^{6},$ the component size
distribution will seem to be a heavy-tailed one.

Figure 3:  Component-size distributions (_solid line_) corresponding to degree
distributions with exponent $\beta=2.6$ ($\alpha=0.6$). In this case, the
component size distributions can not feature a heavy tail, however, depending
upon the scale parameter $s$ a transient asymptote with exponent $-1.6$
(_dashed line_) emerges: _1_) $s=8.3\cdot 10^{-2},$ fast convergence to the
exponential asymptote. _2_) $s=8.3\cdot 10^{-5},$ the distribution transiently
follows what seems to be a heavy tail for $n<10^{6},$ whereas for larger $n$
the theory predicts no heavy tail.  Figure 4: Examples of component size
distributions (_solid lines_) that are associated with heavy-tailed degree
distributions with $\beta=6$ ($\alpha=4$). The dashed lines represent the
asymptotes in accordance with Case B in Table 2. Depending on the sign of
$\theta,$ three asymptotic modes are distinguished: _1)_ $s=1.93:$
$\theta=-0.8$, _1a)_ $s=9.42:$ $\theta=-0.027,$ _2)_ $s=9.69:$
$\theta=1.6\cdot 10^{-7}$, _3)_ $s=10.05:$ $\theta=0.038.$ When $\theta$ is a
small negative number (_curve 1a_), $w(n)$ first decays as $n^{-3/2}$ but
eventually switches to asymptote $n^{-\alpha-1}=n^{-5}.$

Finally, we consider the case when the condition (15) holds for $\beta>4:$
even though $u_{1}(k)$ has finite mean $M$ and variance $\sigma^{2}$ it also
features a heavy tail. Again, as $n\to\infty,$
$x(n)=\frac{n-\mu(n)}{\sigma(n)}$ features the limiting values that are
defined by the sign of $\theta$, see Eq. (20). One would expect that since
$\sigma^{2}$ is finite, this case should be also well approximated with (13).
This is indeed the case for $x(n)\approx 0.$ However, large deviations from
zero $x(n)>>0$ do not follow Gaussian statistics Ramsay (2006); Nacher and
Akutsu (2011), and we approximate $u_{1}^{*n}(k)$ with the Pareto stable law
$u_{1}^{*n}(k)\to\frac{1}{\sigma(n)}G^{P}(\frac{k-\mu(n)}{\sigma(n)},\alpha),\;n\to\infty.$
It turns out that $G^{P}(x,\alpha)$ behaves as the normal distribution for
$x<C,$ where $C$ is a finite positive constant, but features a heavy tail with
the same exponent as $u_{1}(k)$ when $x\to\infty,$ see Ref. Ramsay (2006)
Thus, when $\theta\geq 0$ the component size distribution features asymptotic
modes as in (13), while when $\theta<0$ it features a heavy tail with exponent
$-\alpha-1,$ see Table 2, Case B. Interestingly, when $\theta$ is a small
negative number, $w(n)$ transiently follows one asymptote and then switches to
the other as demonstrated in Figure 4. If there is a process that continuously
changes the degree distribution so that $\theta$ progresses from being
negative to positive, the exponent of the associated component size
distribution will jump from the sub-critical branch, at $\theta<0$, to the
critical one at $\theta=0.$ An example of such transition between two power-
law modes is given in Figure 4, where a component size distribution switches
between power laws with exponents $\eta=3$ and $\eta=1.5$.

##  IV Discussion and Conclusions

The broad generality of the results obtained in the previous section is
achieved due to the fact that the configuration networks are locally tree-like
and have vanishing probability of clustering in the thermodynamic limit, which
allows one to benefit from the available in analytic combinatorics tools. Eq.
(8), that was analysed in the previous section, connects the degree
distribution in a configuration network to the distribution of sizes for
connected components. The main conclusion one may draw from this equation is
that the convolution power provides a smoothing effect. This means that all
points of $u(k),\;k=1,\dots,\infty$ have a significant contribution to the
definition of $w(n)$, but as $n$ increases, the system ‘forgets’ the exact
shape of the degree distribution and the component size distribution tends to
the asymptote, that is defined by only a few parameters. The only information
that is still preserved at the limit $n\to\infty$ is the first three moments
of the degree distribution if such does not feature a heavy tail, see for
example Fig. 1. If $u(k)$ does feature a heavy tail then the information that
characterises the tail becomes also important: that is the scale parameter $s$
and the exponent $\beta$. Depending upon the values of these parameters, many
asymptotical modes exist.

The expression for the asymptote is framed in terms of small deviation
statistics for a sum of random variables and in some cases can be used as a
good approximation for the component size distribution. Table 2 contains the
analytical expressions for the asymptotes. Additionally, supporting code
computing the component size distribution and the corresponding asymptotes is
providedgit . When using the asymptotical expressions to approximate $w(n)$,
one should pay attention to two factors that follow from central limits:
firstly $n$ should be large, secondly the approximation is best for $\theta$
close to zero. Finally, small deviations or a cutoff in a heavy-tailed degree
distribution can trigger considerable and non-trivial changes in $w(n)$, for
instance, the change of the asymptotical mode of the latter.

###  IV.1 Degree distributions with a cutoff

In practice, no empirical degree distribution is a heavy-tailed one. Most of
the ‘real-world’ degree distributions feature a cutoff,
$u(k)=0,\;k>k_{\text{cut}},$ and therefore fail to be heavy-tailed in the
strict sense of the definition (15). It turns out that if a cutoff is featured
at large enough $k_{\text{cut}},$ the above-provided asymptotic analysis still
has a relevant meaning. This situation can be compared to how we commonly
attribute the fractal dimension to real-world geometric objects that fail to
be fractals on infinitesimal scales.

Suppose one applies a cutoff at $k_{\text{cut}}$ to a degree distribution,
$u(k),$ that features a heavy tail. Since $u(k)$ has a finite support, the
asymptote of associated $w(n)$ is covered by Case A (Table 2), however, if
$k_{cut}$ is large, $w(n)$ may also transiently follow the original asymptote.
Instead of an analytical investigation, we demonstrate the influence of the
cutoff with numerical examples obtained by computing (8). This influence
strongly depends on how the sign of $\theta$ is affected by the introduction
of the cutoff. For example, if $\theta>0$ even after the cutoff, the cutoff
will cause more nodes to appear in finite-size components, and thus the
component size distribution will shift towards larger sizes. The opposite case
is valid when $\theta\leq 0$ before (and after) the cutoff, then the cutoff
causes the component size distribution to shift towards smaller sizes. The
third option is when the cutoff changes the sign of $\theta$ form ‘$+$’ to
‘$-$’. In this case, both shifts are possible. Fig. 5 shows how a component
size distribution that corresponds to degree exponent $\beta=3.3$ is affected
by a cutoff with various vales of $k_{\text{cut}}$.

Figure 5: The effect of a cutoff imposed on a heavy-tailed degree distribution
with $\beta=3.3$ ($\alpha=1.3$) and $s=7.73.$ The solid curves correspond to
component size distributions with: _1)_ no cutoff, $\theta>0$; _2)_ cutoff at
$k=1000,$ $\theta=0$; _3)_ cutoff at $k=100,$ $\theta<0.$ The asymptote for
_1_ is covered by Case D, Table 2; while due to the cutoffs the asymptotes for
_2_ and _3_ are covered by Case A. Figure 6: The effect of a cutoff imposed on
a heavy-tailed degree distribution with exponent $\beta=1$ ($\alpha=-1$) and
$s=-2\cdot 10^{-4}$. The following values of the cutoff are considered: _1_)
$k_{\text{cut}}=80$, corresponds to $\theta<0$; _2_) $k_{\text{cut}}=100,$
corresponds to $\theta\approx 0;\emph{3})$ $k_{\text{cut}}=150$ corresponds to
$\theta>0,$ _4_) $k_{\text{cut}}=10^{3};$ _5_) $k_{\text{cut}}=10^{5}.$

###  IV.2 Excess degree distribution with no mean value

In principle, the excess degree distributions that do not have a mean value,
i.e. $\beta<2$, do not fall within any of the above categories. However, if
one introduces a cutoff, $u(k)$ will feature finite moments including,
$\mu_{3}<\infty,$ hence this case should be treated according to Case A of
Table 2. Fig. 6 shows how cutoffs at $k=k_{\text{cut}}$ influence an instance
of component size distribution with $\beta=1$. Unlike as in the previous
example, in which $u(k)$ with no cutoff generates a valid $w(n)$, here the
increase of $k_{\text{cut}}$ results in vanishing probability of finding a
finite-size component at all: for any $n,$ $w(n)\to 0$ when
$k_{\text{cut}}\to\infty$. This illustrates the fact that finite-size
components do not exist for $\beta\leq 1,$ and the whole configuration network
is connected almost surely. Non-existence of finite components for $\beta\leq
1$ also follows from the fact that in this case $\mu_{1}$ diverges and the
point values of $w(n),$ as given below the the definition (9), tend to zero.

Suppose the cutoff in the empirical, heavy-tailed degree distribution is due
to the fact that the network sample has a finite size,
$k_{cut}=N,\;N\neq\infty$, then one may approximate the expected number of
edges in this sample as

$n_{e}=\frac{N\mu_{1}}{2}=\frac{N}{2}\sum_{k=1}^{N}ku(k)\simeq\frac{N}{2}\sum_{k=1}^{N}k^{-\beta+1},\;N>>1,$

so that

$n_{e}\simeq\begin{cases}N(1-N^{2-\beta}),&\beta\neq 2,\\\ N\log
N,&\beta=2.\end{cases}$

Subsequently, three scenarios are possible here:
i) _sparse network,_ $n_{e}=CN,\;C>0,\;\beta>2:$ the asymptotic modes are
given in Table 2;
ii) _semi-dense network,_ either $n_{e}=CN\log N,\;\beta=2$ or
$n_{e}=CN^{3-\beta},\;1<\beta<2:$ the mean value of excess distribution
diverges; there are finite components but no power law in the distribution of
component sizes;
iii) _dense network,_ $n_{e}=CN^{3-\beta},\;\beta\leq 1:$ the mean value of
degree distribution $\mu_{1}\to\infty$, and finite components vanish as
$N\to\infty$.

Figure 7: The correspondence between the exponent $\beta$ in a heavy-tailed
degree distribution, and the exponent $\eta$ in the associated heavy-tailed
component size distribution. The critical brunch corresponds to $\theta=0,$
subcritical branch to $\theta<0.$ Positive $\theta$ is not associated with
heavy-tailed component size distributions.

###  IV.3 The role of the giant component

All the cases presented in Table 2 depend in some way on the value of
$\theta$. This is not a coincidence as the sign of $\theta$ is the indicator
for the giant component existence. If the degree distribution features a heavy
tail with exponent $\beta\geq 3$, depending upon the value of $\theta,$ there
are two possible heavy-tail exponents for the component size distribution:
subcritical branch $\eta=\beta-1$ when $\theta<0$, and critical branch
$\eta=\min\\{\frac{3}{2},\frac{\beta-1}{\beta-2}\\}$ when $\theta=0.$ This
relation is illustrated in Fig. 7, where the component size distribution
exponent $\eta$ is plotted versus the degree-distribution exponent $\beta$. We
can see that if the giant component exists, $\theta>0,$ then irrespectively of
what is the degree distribution, the component size distribution always decays
faster then the power law. Therefore it can be concluded that the giant
component is not compatible with a heavy-tailed component size distribution.
Any degree distribution with $\beta<3,$ leads to a giant component since
$\theta$ can only be positive in this case. Furthermore, if $\beta\leq 1$ the
giant component is also the only component: with probability 1 the
configuration network is fully connected.

Finite moments of $u(k)\;\;$ |  $u(k)$, $k\to\infty$ | $\theta=\mu_{2}-2\mu_{1}$ | Asymptote of $w(n)\;\;$
---|---|---|---
$\mu_{3}<\infty$ | A. $o(k^{-\beta}),\;\beta>4$ | $\theta\neq 0$ | $C_{1}e^{-C_{2}n}n^{-3/2}$
$\theta=0$ | $C_{1}n^{-3/2}$
B. $O(k^{-\beta}),\;\beta>4$ | $\theta<0$ | $C_{3}n^{-\alpha-1}$
$\theta=0$ | $C_{1}n^{-3/2}$
$\theta>0$ | $C_{1}e^{-C_{2}n}n^{-3/2}$
$\mu_{3}=\infty,$ $\mu_{2}<\infty$ | C. $O(k^{-\beta}),\;\beta=4$ | $\theta<0$ | $C_{3}n^{-\alpha-1}$
$\theta=0$ | $C_{1}^{\prime}\frac{n^{-3/2}}{\sqrt{\log n}}$
$\theta>0$ | $C_{1}^{\prime}\frac{n^{-3/2}}{\sqrt{\log n}}e^{-C_{2}^{\prime}\frac{n}{\log n}}$
D. $O(k^{-\beta}),\;3<\beta<4$ | $\theta<0$ | $C_{3}n^{-\alpha-1}$
$\theta=0$ | $C_{4}n^{-\frac{1}{\alpha}-1}$
$\theta>0$ | $C_{5}e^{-C_{6}n}n^{-3/2}$
$\mu_{2}=\infty$ | E. $O(k^{-\beta}),\;\beta=3$ | $\theta>0$ | $C_{7}e^{-C_{8}-C_{9}n^{\frac{2}{\pi}}}n^{\frac{1}{\pi}-2}$
F. $O(k^{-\beta}),\;2<\beta<3$ | $\theta>0$ | $C_{10}e^{-C_{11}n}n^{-3/2}$

$C_{1}=\frac{\mu_{1}^{2}}{\sqrt{2\pi(\mu_{1}\mu_{3}-\mu_{2}^{2})}}$, $\;\;C_{1}^{\prime}=\frac{\mu_{1}}{\sqrt{2\pi s}},$ |  $C_{7}=\frac{\sqrt{2}\mu_{1}}{\pi^{3/2}s}$,
---|---
$C_{2}=\frac{(\mu_{2}-2\mu_{1})^{2}}{2(\mu_{1}\mu_{3}-\mu_{2}^{2})}$, $\;\;C_{2}^{\prime}=\frac{(\mu_{2}-2\mu_{1})^{2}}{2s\mu_{1}^{2}},$ |  $C_{8}=\frac{1}{\pi s}+\frac{1}{2}$,
$C_{3}=\frac{s\mu_{1}^{\alpha+2}\Gamma(\alpha+1)}{(2\mu_{1}-\mu_{2})^{\alpha+1}\Gamma(\alpha)},$ |  $C_{9}=e^{-1-\frac{2}{\pi s}}$,
$C_{4}=\mu_{1}\Gamma\left(1+\frac{1}{\alpha}\right)\sin\frac{\pi}{\alpha}\left(\frac{2\Gamma(\alpha)\sin\frac{\pi\alpha}{2}}{\pi^{\alpha+1}s}\right)^{1/\alpha}\\!,$ | $C_{10}=\frac{\mu_{1}}{\sqrt{2-2\alpha}}\left(\frac{\sqrt{2}\Gamma(\alpha)\sin\frac{\pi\alpha}{2}}{\alpha\pi^{\alpha}s}\right)^{\frac{1}{2\alpha-2}},$
$C_{5}=\frac{\mu_{1}}{\sqrt{\alpha-1}}\left(\frac{2^{2-\alpha}\left(\frac{\mu_{2}}{\mu_{1}}-2\right)^{2-\alpha}\Gamma(\alpha)\sin\frac{\pi\alpha}{2}}{\alpha\pi^{\alpha}s}\right)^{\frac{1}{2\alpha-2}}\\!\\!\\!\\!,\;\;$ | $C_{11}=(1-\alpha)\left(\frac{\sqrt{2}\Gamma(\alpha)\sin\frac{\pi\alpha}{2}}{\pi\alpha^{\alpha}s}\right)^{\frac{1}{\alpha-1}},\;\;\;\;$
$C_{6}=(1-\alpha)\left(\frac{2\left(\frac{\mu_{2}}{\mu_{1}}-2\right)^{\alpha}\Gamma(\alpha)\sin\frac{\pi\alpha}{2}}{\alpha^{\alpha}\pi s}\right)^{\frac{1}{\alpha-1}},$ | $\alpha=\beta-2.$

Table 2: Asymptotic behaviour of component sizes $w(n),$ in terms of degree
distribution parameters: the first three moments $\mu_{1},\mu_{2},\mu_{3}$,
scale parameter $s$, and exponent $\beta$. Supporting source code available in
Ref. git .

###### Acknowledgements.

This work is part of the research program Veni with project number
639.071.511, which is financed by the Netherlands Organisation for Scientific
Research (NWO).

## References

  * Barabási and Albert (1999) A.-L. Barabási and R. Albert, science 286, 509 (1999).
  * Newman  _et al._ (2002) M. E. Newman, D. J. Watts, and S. H. Strogatz, PNAS 99, 2566 (2002).
  * Vázquez (2003) A. Vázquez, Physical Review E 67, 056104 (2003).
  * Kryven  _et al._ (2016) I. Kryven, J. Duivenvoorden, J. Hermans, and P. D. Iedema, Macromolecular Theory and Simulations 25, 449 (2016).
  * Bhat  _et al._ (2016) U. Bhat, P. L. Krapivsky, R. Lambiotte, and S. Redner, Physical Review E 94, 062302 (2016).
  * Newman and Watts (1999) M. E. J. Newman and D. J. Watts, Physical Review E 60, 7332 (1999).
  * Molloy and Reed (1998) M. Molloy and B. Reed, Combinatorics, probability and computing 7, 295 (1998).
  * Newman (2010) M. Newman,  _Networks: an introduction_ (Oxford university press, Oxford, 2010).
  * Chung  _et al._ (2003) F. Chung, L. Lu, and V. Vu, PNAS 100, 6313 (2003).
  * Kryven (2016) I. Kryven, Phys. Rev. E 94, 012315 (2016).
  * Molloy and Reed (1995) M. Molloy and B. Reed, Random structures & algorithms 6, 161 (1995).
  * Kryven and Iedema (2014) I. Kryven and P. D. Iedema, Chemical Engineering Science 126, 296 (2014).
  * Biemann and van den Bosch (2011) C. Biemann and A. van den Bosch,  _Structure discovery in natural language_ (Springer Science & Business Media, 2011).
  * Newman  _et al._ (2001) M. E. J. Newman, S. H. Strogatz, and D. J. Watts, Physical review E 64, 026118 (2001).
  * Newman (2007) M. E. J. Newman, Phys. Rev. E 76, 045101 (2007).
  * Yook  _et al._ (2002) S.-H. Yook, H. Jeong, and A.-L. Barabási, PNAS 99, 13382 (2002).
  * Ravasz and Barabási (2003) E. Ravasz and A.-L. Barabási, Physical Review E 67, 026112 (2003).
  * Eguiluz  _et al._ (2005) V. M. Eguiluz, D. R. Chialvo, G. A. Cecchi, M. Baliki, and A. V. Apkarian, Physical review letters 94, 018102 (2005).
  * Fu  _et al._ (2008) F. Fu, L. Liu, and L. Wang, Physica A: Statistical Mechanics and its Applications 387, 675 (2008).
  * Adamic and Huberman (2002) L. A. Adamic and B. A. Huberman, Glottometrics 3, 143 (2002).
  * Timár  _et al._ (2016) G. Timár, S. N. Dorogovtsev, and J. F. F. Mendes, Physical Review E 94, 022302 (2016).
  * Dorogovtsev  _et al._ (2000) S. N. Dorogovtsev, J. F. F. Mendes, and A. N. Samukhin, Physical review letters 85, 4633 (2000).
  * Bergeron  _et al._ (1998) F. Bergeron, G. Labelle, and P. Leroux,  _Combinatorial species and tree-like structures_ (Cambridge University Press, Cambridge, 1998).
  * Kryven and Iedema (2013) I. Kryven and P. Iedema, Macromolecular Theory and Simulations 22, 89 (2013).
  * Kruchinin (2015) D. V. Kruchinin, Advances in Difference Equations 2015, 1 (2015).
  * Ma and Liu (2004) N.-Y. Ma and F. Liu, Applied mathematics and computation 158, 225 (2004).
  * Ma and King (2002) N.-Y. Ma and R. King, Applied mathematics and computation 133, 83 (2002).
  * Gnedenko and Kolmogorov (1954) B. V. Gnedenko and A. Kolmogorov,  _Limit distributions for sums of independent variables_ (Addison-Wesley, Cambridge, 1954).
  * Uchaikin and Zolotarev (1999) V. V. Uchaikin and V. M. Zolotarev,  _Chance and stability: stable distributions and their applications_ (Walter de Gruyter, Urecht, 1999).
  * Ramsay (2006) C. M. Ramsay, Communications in Statistics—Theory and Methods 35, 395 (2006).
  * Nacher and Akutsu (2011) J. Nacher and T. Akutsu, Physica A: Statistical Mechanics and its Applications 390, 4636 (2011).
  * (32) “Matlab/GNU Octave source code for cacluating component size distribution in configuraiton network,” https://github.com/ikryven/PhysRevE_2017_GECS.


B. Yip et al.

# Image Pre-processing Using OpenCV Library on MORPH-II Face Database

B. Yip, R. Towner, T. Kling, C. Chen, and Y. Wang  Department of Mathematics
and Statistics, The University of North Carolina Wilmington, Wilmington, NC
28403, USA<EMAIL_ADDRESS>
###### Abstract.

This paper outlines the steps taken toward pre-processing the 55,134 images of
the MORPH-II non-commercial dataset. Following the introduction, section two
begins with an overview of each step in the pre-processing pipeline. Section
three expands upon each stage of the process and includes details on all
calculations made, by providing the OpenCV functionality paired with each
step. The last portion of this paper discusses the potential improvements to
this pre-processing pipeline that became apparent in retrospect.

###### Key words and phrases:

Image pre-processing; OpenCV Library; Image Analysis; MORPH-II.

###### 2010 Mathematics Subject Classification:

68U10; 97R50

##  1\. Introduction

The MORPH data is one of the largest publicly available longitudinal face
databases (Ricanek and Tesafaye, 2006). Since its first release in 2006, it
has been cited by over 500 publications. Multiple versions of MORPH have been
released, but for our face image analysis study, we use the 2008 MORPH-II non-
commercial release. The MORPH-II dataset includes 55,134 mugshots with
longitudinal spans taken between 2003 and late 2007. For each image, the
following metadata is included: subject ID number, picture number, date of
birth, date of arrest, race, gender, age, time since last arrest, and image
filename. Because of its size, longitudinal span, and inclusion of relevant
metadata, the MORPH-II dataset is widely utilized in the field of computer
vision and pattern recognition, including a variety of race, gender, and age
face imaging tasks.

(a) Variation in head tilt.

(b) Variation in camera distance.

(c) Variation in illumination.

(d) Variation in overall appearance.

Figure 1. Different types of variation present in the MORPH-II dataset.

However, despite the fairly standard format of police photography, many of the
images vary greatly in terms of head-tilt, camera distance, and illumination.
A great number of the images also contain large, empty backgrounds or excess
occlusion that add a corresponding amount of noise to the data. This
longitudinal dataset had an average of approximate 4 images per subject and
some appearances varied greatly from one image to the next. Preliminary
results showed women having an increased overall variation in their images due
to changes in makeup and hairstyle. Figure 1 showcases some variety found
during initial examination of the image dataset.

Consequently, the pre-processing step is crucial for the image analysis on the
MORPH-II dataset. For our purposes, we utilized the Open Source Computer
Vision 2 (OpenCV) library in Python to extract the face from each mugshot
using the image vectors (Bradski, 2000). The stages of this process are
outlined in section 2.

##  2\. Procedural Overview

This section provides a global description for the six stages of our pre-
processing algorithm. The premise is to minimize image noise by the use of
bounding boxes around necessary region of interest (ROI). Both Figures 2 and 3
are visual representations of each step and what is accomplished, from
different prospects.

Initial face and eyes detection.

Rotation.

Face and eye re-detection.

Cropping and scaling

Pre-processed image.

Figure 2. Stages of the pre-processing pipeline with successful face and eye
detection

Note: While the intermediate stages of the process in Figure 2 are shown with
color images, all computer vision tasks were done only on the grayscale
versions of each image (converted with OpenCV).

(a) Original image.

(b) Grayscale image.

(c) Initial face detection.

(d) Eye detection.

(e) Rotation & re-detection.

(f) Cropping and scaling.

Figure 3. Face preprocessing pipeline with successful face and eye detection.

###  2.1. Grayscale

Research in computer vision showed converting images to grayscale increased
the accuracy of locating the necessary facial features as it reduced the
effect illumination variance had on the images. Doing so was the first
introduction to the OpenCV library. For each image, we utilized the OpenCV
function $cv2.cvtColor(src,channel)$ where $src$ is the input image and
$channel$ represents the color channel for the output image. In our case we
used $Color_{B}GR2Gray$. This results in the new image pixel value, $Y$, where
$Y=0.299\cdot R+0.587\cdot G+0.114\cdot B$, and $R$, $G$, $B$ represent Red,
Green and Blue respectively. The values of $R$, $G$, $B$, are from the
original image pixel. Effects are shown in Figure 3(B).

###  2.2. Face Detection

The initial face detection step located and marked the position of a face
within an image, which can be seen in Figure 3(C). This eliminates backgrounds
and hairstyles, which are image properties that are not useful for computer
vision tasks. This procedure increased the accuracy of future detection steps.
If a face was not successfully detected, the image was stored in a face not
found (fnf) folder for manual detection later on. Both face and eye detection
steps were accomplished using Haar-Feature based cascade classifiers from
OpenCV. The function used was
$cv2.cascadeClassifier.detectMultiScale(src,sf,mn)$ where $src$ is the input
image, $sf$ is the scale factor at each image scale, and $mn$ is the minimum
amount of neighbors each candidate face rectangle should acquire.

Note that: Both face and eye detection were done using the the Haar-feature
based Cascade Classifiers from OpenCV (the .xml files can be obtained from the
OpenCV GitHub repository). For our purposes, we only had to adjust the
parameters of the OpenCV detection function to locate the face in each image.
Eye detection, however, required additional steps.

###  2.3. Eye Detection

We implemented a very similar algorithm for locating the eyes, illustrated in
Figure 3(D). After the face is detected, the eyes are located within the
region of interest (ROI) determined by the face (i.e. within the bounding box
for the face), and the eye centers are computed as the center of the bounding
box for each eye. The new domain and range of the image matrix come from the
bounding box around a face that was found successfully. We then marked
bounding boxes around each eye with the same Haar cascade function from OpenCV
in section 2.2 with different parameter values to account for the smaller
scope.

In many cases, wrinkles, shadows, and other facial blemishes were detected as
eyes. To account for this, we implemented two conditions to eliminate as many
of the incorrect eye detections as possible: if 1) the angle between eye
centers was greater than fifteen degrees, or 2) the interocular distance
(number of pixels between eye centers) was less than one fifth of the image
width, the located features were discarded. Following the test of these
conditions, a while loop conditioned on successful eye detection was used to
refine the parameters of the detect function when eyes were not found. If this
too proved unsuccessful, the image was stored for manual detection.

When both eyes were successfully found, figure 2(D), we captured the
coordinate location of the right $(x_{r},y_{r})$ and left $(x_{l},y_{l})$ eye
centers by calculating the center location of the new bounding boxes. These
eye centers were crucial for future steps.

###  2.4. Rotation

Given successful eye-detection, the image is rotated based on the angle
between the eye centers, as illustrated in Figure 3(E). Rotating the image
began with the eye centers $(x_{r},y_{r})$ and $(x_{l},y_{l})$. We added a
conditional to ensure the left eye stored in $(x_{l},y_{l})$ was actually the
left eye of the subject. $\theta$ was then calculated by subtracting the left
eye from the right eye, $(x_{r}-x_{l},y_{r}-y_{l})$, yielding the displacement
between the eyes. We had to convert $\theta$ to the complex plane to allow use
of the numpy angle function. $\theta$ was then a parameter used in the
getRotationMatrix2D OpenCV function to create the necessary transformation
matrix, M.

    ΨΨΨΨΨcv2.getRotationMatrix2D(center, theta, scale)
    ΨΨΨΨΨcenter = center of rotation source image
    ΨΨΨΨΨscale = scale factor
    ΨΨΨΨ

$\begin{bmatrix}\alpha&\beta&(1-\alpha)\cdot center.x-\beta\cdot center.y\\\
-\beta&\alpha&\beta\cdot center.x+(1-\alpha)\cdot center.y\\\ \end{bmatrix},$

where, $\alpha=scale\cdot\cos(angle),\beta=scale\cdot\sin(angle)$.

After the transformation matrix is calculated, it is applied to each pixel in
the source image to produce the necessary rotated image. The OpenCV function
warpAffine reconstructs the source image by producing the desired rotated
image, dst:

    ΨΨΨΨΨcv2.warpAffine(src, M, (column, row))
    ΨΨΨΨ

$dst(x,y)=src(M_{11}x+M_{12}y+M_{13},M_{21}x+M_{21}y+M_{21}).$

###  2.5. Face and Eye Re-detection

Following rotation, the face and eyes are re-detected as above, as illustrated
in Figures 2 and 3(E). Images with unsuccessfully detected faces were stored
in a ”fnf_r” (face not found, rotated image) folder for manual detection later
on. Images with undetected eyes were stored similarly, in a ”enf_r” (eyes not
found, rotated image) folder.

###  2.6. Cropping and Scaling

After the eye centers were successfully re-located in the rotated image, a new
bounding box for the face was determined based on the interocular distance.
This step ultimately found the white bounding box in Figure 3(F) for the
cropped image based on the interocular distance. Eye centers were relocated in
the rotated image using the Haar cascade in section 2.3. The new eye center
coordinates gave the interocular distance used for defining the height and
width of the new bounding box. Finding the upper left corner of the box began
with subtracting one interocular distance from the x-value midpoint of the
eyes. The y-coordinate was calculated by adding four fifths times the
interocular distance to the current eye height $(x=eyemidpoint-
interoc,y=eyeheight+0.8\cdot interoc)$. [Note: 0.8 was chosen as the best
option for capturing an appropriate region of the face]. The final bounding
box was a slice of the rotated image with proportions 2:2.35 times the
interocular distance.

The image was then cropped according to this frame and scaled down to 70
pixels tall by 60 pixels wide, using the following OpenCV function:

    ΨΨΨΨΨcv2.resize(cropped_img, (60, 70))
    ΨΨΨΨ

If an image could not be reduced to this size (i.e. something went wrong
earlier on), it was stored in another folder for manual detection.

###  2.7. Manual Pre-processing

The images in which the face or eyes were not successfully detected were
handled by manually clicking the eye centers of each subject. These new eye
centers were then used for rotating the image and the rest of the pipeline
followed suit. The images below in Figure 4 are examples of problem images
that included errant detection.

Figure 4. Examples of problem images and errant detection

###  2.8. In Retrospect

When performing the manual eye-detection, the images were rotated based on the
eye centers. However, we avoided re-clicking the eyes in the rotated image by
simply applying the rotation matrix to the coordinates of the eye centers in
the unrotated image. Had this been recognized in the original code, the face
and eye re-detection step could have been skipped entirely. This would likely
have drastically reduced run time and decreased the number of images
necessitating manual eye detection.

##  3\. Conclusion

When performing the manual eye-detection, the images were rotated based on the
eye centers (as above). However, we avoided re-clicking the eyes in the
rotated image by simply applying the rotation matrix to the coordinates of the
eye centers in the unrotated image. Had this been recognized in the original
code, the face and eye re-detection step could have been skipped entirely
(effectively merging stages 2.3 and 2.4). This would likely have drastically
reduced run time and decreased the number of images necessitating manual eye
detection.

##  4\. Acknowledgments

This material is based in part upon work supported by the National Science
Foundation under Grant Numbers DMS-1659288. Any opinions, findings,
conclusions, or recommendations expressed in this material are those of the
author(s) and do not necessarily reflect the views of the National Science
Foundation.

## References

  * Bradski (2000) Bradski, G. (2000).  The OpenCV Library.  Dr. Dobb’s Journal of Software Tools.
  * Ricanek and Tesafaye (2006) Ricanek, K. and Tesafaye, T. (2006).  Morph: A longitudinal image database of normal adult age-progression.  In Automatic Face and Gesture Recognition, 2006. FGR 2006. 7th International Conference on, pages 341–345. IEEE.