\section{Introduction}

Deep convolution neural networks have shown impressive results in many domains such as in computer vision problems. The fundamental improvement of current deep learning methods is that, unlike earlier more shallow network layers, deep learning automatically identifies appropriate stratification of relevant attributes to construct a predictive model [2]. But despite improvements on quality and accuracy of predictions, these models are hard to interpret, especially when they become deeper and training weight adjustments propagate to fully connected layers. The problem of interpreting (or debugging) becomes extremely important in more sensitive domain such as dealing with medical data, i.e. any wrong decision might directly  affect patient's well-being. We therefore need to be able to build models which are interpretable and can support the explanation of their decisions. While neural networks can achieve high accuracy, they are not easily interpretable. Improved interpret-ability will help inform tradeoffs, e.g., methods like Bayesian nets are more easily interpretable but less accurate.

Recently, many methods have been proposed to address this trade-off  of interpret-ability [8-10]. One approach is based on computing the gradient of the output layer with respect to the feature maps [3], which is a relaxation generalization of [5]. The gradient is then back-propagated to a feature map layer, to try and distinguish the discriminative pixels for any input image. The gradient is determined based on network activation; to focus attention to where network activation provides insight into classifier output, we here consider two popular methods.

In the following, we explain our approach by describing a simple method to interpret deep neural network prediction. Our approach produces a map which highlights the salient region for the input image that represents the evidence related to the activations of the output layer of the network, i.e. why this network arrived at this decision?

The reminder of this paper is organized as follows. Section 2 presents our proposed approach. In Section 3 we describe the experiments with our initial framework, and  finally, Section 4 concludes with our preliminary results, followed by a discussion of extensive future work.

\section{Approach}
In the case of image classification, we focus on the loss of spatial information in fully connected layers of deep convolutional networks, and identify an approximate visual explanation using the feature maps $X$ at a particular convolution layer $l$. To be able to explain a prediction of the black box function $f$, we rely on computing the Kullback-Leibler (KL) divergence gradient of the raw class scores $y^{'}$ with respect to ground truth $y$, to estimate the discriminative localization map. We use this gradient as evidence for the network prediction. To achieve this, we  first compute the joint probabilities of  $y^{'}$ as follows:

\begin{equation}
\begin{split}
 p_{ij} =   \frac{(1+ ||k_{i}-k_{j}||^2)^{-1} }{\sum_{u \neq v} ( 1+||k_{u}-k_{v}||^2)^{-1}}
 \end{split}
\end{equation}

Here $p_{ij}$ denotes the joint probabilities, $k$ is the raw class scores before softmax , $i$ indexes a neuron value and  $\sum_{u \neq v}$ combines all the values. For the ground truth we estimate the pairwise affinities with perplexity.  We then compute the KL-divergence gradient i.e. $\frac{\delta y^{'}}{\delta y} \Rightarrow z$ derived here [6]. We also normalize the gradient to a zero mean and unit variance as follows:

\begin{equation}
\begin{split}
 \alpha =   \frac{z - \mu }{\sigma z}
 \end{split}
\end{equation}

The obtained weights $\alpha$ capture the relevant information in the feature maps acquired by the network. These weights are applied to every feature map $x_{i} \in X$ as to identify the discriminative pixels which influence the final prediction output as follows:

\begin{equation}
\begin{split}
 E_{KL-divergence} =  \sum_{i}\sum_{j} x_{i} * |\alpha_{j}|
 \end{split}
\end{equation}

The $\sum_{i}$ runs over all the features $x_{i} \in X$ and $j$ runs over all the obtained weights $\alpha$. The output of (3) is the localization map which can be used as the evidence to explain the prediction of the black box function. Finally, we normalize the heat-map $E_{KL-divergence}$ for visualization. Algorithm 1 summarizes the implementation method and the overall methodology is depicted in Figure \ref{fig:main_d}.

\begin{algorithm}
\caption{Proposed approach}\label{algorithm1}
\begin{algorithmic}
\INPUT image, ground truth $y$
\OUTPUT Discriminative localization map  $\Rightarrow E_{KL-divergence}$\\
Apply a single forward-pass to estimate $\Rightarrow y^{'}$\\
Compute the joint probabilities for both $y^{'}$ and $y$\\
Compute the gradient and normalize using {\color{red} {$(2)$}} $\Rightarrow \alpha$ \\
initialize  $E_{KL-divergence}$ to zero
\For{$i = 1$ \textbf{to} ${nFeatureMaps}$}
\State Initialize temp to zero
\For{$j = 1$ \textbf{to} ${sizeof\alpha}$}
\State $temp \gets  temp +(x_{i}*|\alpha_{j}|)$
\EndFor
\State\hskip-\ALG@thistlm {\textbf{end for}}
\State $E_{KL-divergence} \gets E_{KL-divergence} + temp$
\EndFor

\State\hskip-\ALG@thistlm {\textbf{end for}}

\end{algorithmic}
\end{algorithm}

\begin{figure}[H]

\centering
\includegraphics[width=13cm,height=5.5cm]
{framework.png}
\caption{Our proposed methodology.}
\label{fig:main_d}
\end{figure}

\section{Experiments}
To evaluate the proposed method in the context of visual explanation, we select images from the common objects in context (COCO) challenge set [1]. For the learned models, we employed the publicly available pre-trained model VGG-16 [4] . The initial results of our visualization algorithm are illustrated in Figure \ref{fig:vggexamp}. \par
\begin{figure}[H]
\centering
\includegraphics[width=13cm,height=5cm]
{vgg-example.png}
\caption{Explaining the decisions made by the VGG-16 (park bench, street sign, racket, cockatoo, traffic light and chihuahua), our approach highlights the most discriminative region in the image.}
\label{fig:vggexamp}
\end{figure}
We also compared our explanation algorithm with other similar methods i.e., guided back- propagation [9] and gradient weighted class activation [3]. The comparison results are depicted in Figure \ref{fig:compmethods}.
\begin{figure}[H]

\centering
\includegraphics[width=13cm,height=5cm]
{comp1.png}
\caption{Results on randomly chosen images from COCO visualizing VGG-16 predictions (street sign, racket and racket).}
\label{fig:compmethods}
\end{figure}
\subsection{Comparing two deep neural networks}
It is also important to understand the influence of the network architectures on the visualization algorithms. For the purpose of this experiment, we tested our algorithm on VGG-16 and AlexNet [7] networks. The comparison results are shown in Figure \ref{fig:compnetworks}. From these results, we can see that the quality of the visualization depends on the overall accuracy of the network. Since VGG-16 is more reliable, we can observe that the network pays attention to a smaller salient area, compared to the AlexNet network. Therefore, the better the accuracy the better the visualization interpretation.
\begin{figure}[H]

\centering
\includegraphics[width=13cm,height=5cm]
{comp.png}
\caption{Comparing visualization results of two networks.}
\label{fig:compnetworks}
\end{figure}
\section{Conclusion}
We introduce a new approach to explaining deep neural networks prediction. This particular solution focuses identification of salient regions by computing the KL-divergence gradient from the ground truth label and the output layer. Our preliminary experiments on image classification problems compare the performance of two networks. We believe that this work can further be extended to evaluate decisions in a variety of applications of deep neural networks, in image classification, speech, and text.  We are especially keen to consider text domains, as interpretability will expose patterns of semantic relatedness beyond that of images.

\setcounter{secnumdepth}{0}
\section{Acknowledgment}

We thank colleagues from the Alberta Machine Intelligence Institute for their discussions and advice. This work was supported by NSERC and AMII.

\section*{References}
\medskip

\small
[1] Lin T.Y., Maire M., Belongie S., Hays J. Perona, P. Ramanan D., Dollár P. and Zitnick C.L. Microsoft coco: Common objects in context. In \textit{ECCV}, 2014.

[2] LeCun Y., Bengio Y., and Hinton G. Deep learning. \textit{Nature} pages 436-444, 2015.

[3] Selvaraju R. R., Das A., Vedantam R., Cogswell M. Parikh, D. and Batra D. Grad-CAM: Why did you say that?. In \textit{Workshop on Interpretable Machine Learning in Complex Systems, NIPS}, 2016.

[4] Simonyan K., and Zisserman A. Very deep convolutional networks for large-scale image recognition. \textit{arXiv preprint arXiv:1409.1556}, 2014.

[5] Zhou B., Khosla A., Lapedriza, A. Oliva, A. and Torralba A. (2016). Learning deep features for discriminative localization. In \textit{Proc. CVPR}, pages 2921-2929, 2016.

[6] Maaten L. V. D., and Hinton G. Visualizing data using t-SNE. \textit{Journal of Machine Learning Research}, pages 2579-2605, 2008.

[7]Krizhevsky A., Sutskever I., and Hinton G. E. Imagenet classification with deep convolutional neural networks. In \textit{Proc. NIPS}, pages 1097-1105, 2012.

[8] Zintgraf L. M., Cohen T. S., Adel T., and Welling, M. Visualizing deep neural network decisions: Prediction difference analysis. In \textit{Proc. ICLR}, 2017.

[9] Springenberg J. T., Dosovitskiy A., Brox T., and Riedmiller M. (2014). Striving for simplicity: The all convolutional net. \textit{arXiv preprint arXiv:1412.6806 }, 2014.

[10] Mahendran A., and Vedaldi A. Salient deconvolutional networks. In \textit{Proc. ECCV}, pages 120-135, 2016.

\end{document}


\section{Introduction}
A plethora of systems in engineering, finance, and physical sciences exhibit stochastic dynamics \cite{allen07,lande03,malliaris82,gardiner86,oksendal03}. These systems are mathematically characterized in terms of how their probability density functions (or characteristic functions) evolve with time. However, it is generally difficult to solve for these quantities other than a few special cases. One often resorts to stochastic simulations or other approximation schemes \cite{rif96, socha2008linearization, Socha05,hes05,hespanha2006modelling}. Here, we propose a method to estimate stationary characteristic functions for a class of stochastic systems. The state of this class evolves as per a stochastic differential equation driven by both white noise and L\'evy noise. For such systems, a partial differential equation (PDE) governing the time evolution of the characteristic function can be written \cite{gri04,gri13}. In stationary state, this PDE becomes an ordinary differential equation (ODE) whose order is determined by the degrees of polynomials in drift and diffusion terms. The solution to this ODE requires as many unknown coefficients as its order. Our method relies upon relating these coefficients to the stationary moments, and determining the bounds on them by exploiting the positive definiteness of the characteristic function.

Since the proposed method involves estimating some moments of a stochastic system, it is closely related with moment closure methods \cite{KuntzStatMoments16,soltani2015conditional,SinghHespanhaDM}. Specifically, recent works have shown that the positive semidefiniteness of the moments can be used to find exact lower and upper bounds on moments of a stochastic system \cite{lgs16, lgs17, gvl17, kts17}. Here, we show that such bounds can also be obtained by using the fact that the characteristic function is a positive definite function. The proposed method is superior in the sense that it links the moment bounds to approximate the stationary characteristic function and thereby the stationary probability density function. This work is also related with another category of works wherein the characteristic function is used to describe stochastic dynamical systems, including those driven by Poisson/L\'evy noise \cite{epp93, gri96, wjb99, cgk02, gri04,  dgs07, sag07, dhh09, cot11, gri13}. However, these works are typically restricted to cases where the characteristic function can be exactly solved for. Our method approximates the characteristic function beyond these.

The rest of the paper is organized as follows. Section~\ref{sec:background} provides background results on characteristic functions and their relevant properties. Section~\ref{sec:EstimatingCF} describes generators of stochastic dynamical system and proposes our method to estimate moments and characteristic function in stationary state. Section~\ref{sec:examples} illustrates the approach via examples, and finally section~\ref{sec:discussion} concludes the paper.
\paragraph*{Notation}
The set of real numbers, and natural numbers are respectively denoted by $\mathbb{R}$ and $\mathbb{N}$.
The state of a stochastic process is denoted by $X \in \mathbb{R}$, with a specific value taken by it denoted by $x$.
The squared root of $-1$ is denoted by the letter $j$. The expectation operator is denoted by $\mathbb{E}$. Characteristic function is denoted by $\varphi(\omega)$.

\section{Mathematical Preliminaries}\label{sec:background}
In this section, we provide background results pertaining characteristic functions and their properties. The reader is referred to \cite{Fel682, ush99, biz00,sas13} for proofs and more details. For simplicity, we only consider one dimensional systems.

For a univariate random variable $X$, with distribution function $F$, the characteristic function $\varphi: \mathbb{R} \to \mathbb{C}$ is defined by
\begin{equation}
\varphi(\omega):=\mathbb{E}\left(e^{j\omega X}\right)=\int_{-\infty}^{\infty}e^{j\omega x} dF(x).
\end{equation}
For any random variable, the characteristic function always exists and it uniquely determines the distribution. If $F(x)$ has a density $h(x)$ then $\varphi(\omega)$ can be written in terms of $h$
\begin{equation}
\varphi(\omega) = \int_{-\infty}^{\infty}e^{j\omega x} h(x) dx.
\end{equation}
Given $\varphi(\omega)$, the distribution function $F(x)$ and/or the density $h(x)$ can be obtained via inversion.

As an immediate consequence of its definition, the characteristic function of a random variable has the following properties
\begin{enumerate}
\item $\varphi(0)=1$.
\item $\left \lvert \varphi(\omega) \right \rvert \leq 1$.
\item $\varphi(-\omega)=\overline{\varphi(\omega)}$.
\item If the random variable has a finite $m^{th}$ order moment, then it is given by
\begin{equation}\label{eq:cfMom}
\mathbb{E}\left(X^m\right)=(-j)^m \frac{d^m \varphi}{d \omega^m}\lvert_{\omega=0}.
\end{equation}
\end{enumerate}
Another important property of a characteristic function that we particularly use in this work is that it is a positive definite function.
The following theorem of Bochner forms basis of our analysis.
\begin{theorem}[Bochner]
A continuous complex-valued function $\varphi$ on $\mathbb{R}$ is a characteristic function if and only if $\varphi(0)=1$ and $\varphi$ is positive definite \cite{ush99, biz00}.
\end{theorem}

The positive definiteness of a function required by this theorem is defined as follows.
A complex valued function $\varphi:\mathbb{R} \to \mathbb{C}$ is said to be positive definite if the inequality
\begin{equation}
\sum_{i,j=1}^{q}\varphi(\omega_i-\omega_j)c_i\overline{c_j} \geq 0
\end{equation}
holds for every positive integer $q$, for all $\omega_1, \ldots, \omega_q \in \mathbb{R}$, and for all $c_1, \ldots, c_q \in \mathbb{C}$.
In other words, $\varphi$ is positive definite if and only if the matrix
\small{
\begin{equation}\label{eq:M}
M=
\begin{bmatrix}
\varphi(0) & \varphi(\omega_1-\omega_2) & \ldots &  \varphi(\omega_1-\omega_q) \\
\varphi(\omega_2-\omega_1)  & \varphi(0) & \ldots &  \varphi(\omega_2-\omega_q) \\
\vdots & \vdots & \ddots & \vdots \\
\varphi(\omega_q-\omega_1)  & \varphi(\omega_q-\omega_2)  & \ldots &\varphi(0) \\
\end{bmatrix}
\succeq 0.
\end{equation}
}\normalsize
for an arbitrary choice of $q \in \mathbb{N}$ and $\omega_1, \ldots, \omega_q \in \mathbb{R}$. Here $\succeq 0$ denotes the positive semidefiniteness.

As a consequence of this definition, a variety of properties of the characteristic function (including properties 2 and 3 above) can be established by choosing some test points and enforcing $M \succeq 0$. For example, consider the case $q=1$. Then, we must have that
\begin{equation}
\varphi(0) \geq 0.
\end{equation}
Likewise, for $q=2$, we should have
\begin{equation}
\begin{bmatrix}
\varphi(0) & \varphi(\omega_1 - \omega_2) \\
\varphi(\omega_2-\omega_1) & \varphi(0)
\end{bmatrix}
\succeq 0.
\end{equation}
Without loss of generality, we can assume $\omega_1=0$ and $\omega_2=\omega \in \mathbb{R}$. Then, we should have $\varphi(-\omega) = \overline{\varphi(\omega)}$ and $|\varphi(\omega)|^2 \leq \varphi(0)^2$.

\section{Stationary Characteristic Function of a Stochastic Process} \label{sec:EstimatingCF}
In this section, we describe the generator of a stochastic process and get a PDE for time evolution of the characteristic function. Then, we use this PDE to obtain an ODE for the stationary characteristic function and discuss its solution.

\subsection{Stochastic dynamical system and its generator}
We consider stochastic differential equations driven by L\'evy
noise, known as \emph{It\^o-L\'evy} processes. In
this paper, we will restrict ourselves to It\^o-L\'evy processes of the form:
\begin{multline}\label{eq:JD}
  dX(t)=f(X(t))dt+g(X(t))dW(t) \\
  + \int_R y \tilde{N}(dt,dy) .
\end{multline}
Here $X \in \mathbb{R}$ denotes the state, $W$ is the Wiener process,
$\tilde{N}$ is called a \emph{compensated Poisson measure}, and $f:\mathbb{R} \to \mathbb{R}$, $g:\mathbb{R} \to \mathbb{R}$ characterize the dynamics.
In order for \eqref{eq:JD} to be well-defined we must assume that the
driving L\'evy process has finite variance.
See \cite{oksendalapplied2009,applebaum2009levy} for more details on the formalism.

The jump process generalizes the jumps of a Poisson process. In a Poisson
process, all of the jumps have value $1$, and $N(t)$ simply counts the
number of jumps. In a L\'evy process, the jumps can take on arbitrary
values. The \emph{Poisson random measure}, $N(t,\cdot)$, is a
measure-valued stochastic process such that for a Borel set $S$, $N(t,S)$ is a Poisson
process that counts the number of jumps that took values in the set
$S$. The intensity of the Poisson random measure given by the \emph{L\'evy
  measure}, $\nu$:
\begin{equation*}
  \mathbb{E}[N(t,S)] = t \nu(S)
\end{equation*}
The compensated Poisson measure from \eqref{eq:JD} is the
measure-valued stochastic process defined by
$\tilde{N}(t,S) = N(t,S) - t \nu(S)$.

Let $\psi$ be a  twice continuously differentiable function.  A
standard result in stochastic differential equations \cite{oksendalapplied2009} shows that the
dynamics of $\mathbb{E}[\psi(X(t))]$ are given by:

\begin{equation}
\frac{d}{dt} \mathbb{E}[\psi(X(t))] = \mathbb{E}[L\psi(X(t))]
\end{equation}
where $L\psi$ is called the \emph{generator} of the process. The
generator is defined by

\begin{multline}
\label{eq:generator}
L\psi(X)
=
\frac{\partial \psi(X)}{\partial X}
f(X)
+\frac{1}{2}
\frac{\partial^2 \psi(X)}{\partial X^2}
g(X)^2
 + \\
\int_\mathbb{R}\left(\psi(X+y) - \psi(X)-
  \frac{\partial \psi(X)}{\partial X} y
  \right) \nu(dy).
\end{multline}

In the following, we use the generator to find a PDE that governs the
evolution of the characteristic function. We will see that for the
commonly-studied L\'evy processes, formulas exist to enable explicit
calculation of the required integral.

\subsection{Characteristic function of the process}
Consider the stochastic dynamics defined in \eqref{eq:JD}. We restrict
ourselves to the cases for which the functions $f$, and $g$ are
polynomials in the state $X$. Our goal is to compute the time
evolution of
\begin{equation}
\psi(X(t))=e^{j \omega X(t)}.
\end{equation}

The following theorem shows how a partial differential equation governing the evolution of the characteristic function can be obtained for \eqref{eq:JD}. We wish to point out that it is presented here for a formal statement, and it has been used in some form or other in several works, e.g., see \cite{gri13}.

\begin{theorem}\label{thm:ode}
Consider a one dimensional stochastic process defined in \eqref{eq:JD}. Let $f$ and $g^2$ be finite polynomials of degrees $d_f \in \mathbb{N}$ and $d_{g} \in \mathbb{N}$ respectively. Assuming that a stationary distribution exists, the characteristic function of the stationary distribution satisfies the following ordinary differential equation of order $n=\max \{d_f, d_g\}$:
\begin{equation}\label{eq:StODE}
j \omega  \sum_{l=0}^{d_f}a_{f_l} j^{-l} \frac{\partial^l \varphi}{\partial \omega^l} +\frac{1}{2} (j \omega)^2  \sum_{l=0}^{d_g}a_{g_l}  j^{-l} \frac{\partial^l \varphi}{\partial \omega^l}+\eta(\omega)\varphi=0,
\end{equation}
where $\eta(\omega)= \int_\mathbb{R}\left(e^{j\omega y} - 1 -
  j \omega y
  \right) \nu(dy)$.
\end{theorem}
\begin{proof}
Since we assume that $f$, and $g^2$ are polynomials, without loss of generality we can take their forms to be
\begin{equation}
f(X)=\sum_{l=0}^{d_f}a_{f_l} X^{l}, \quad g^2=\sum_{l=0}^{d_g}a_{g_l} X^{l},
\end{equation}
where $a_{f_l} \in \mathbb{R}$, and $a_{g_l} \in \mathbb{R}$ are coefficients. Taking $\psi(X)=e^{j\omega X},\; \omega \in \mathbb{R}$, we have that
\begin{align}\label{eq:derivativesPsi}
\frac{\partial^m \psi}{\partial X^m}=(j\omega)^m\psi, \quad \frac{\partial^m\psi}{\partial \omega^m}=(j X)^m\psi \quad \forall m \in \mathbb{N}
\end{align}

Using these, we can write $L\psi(X)$ as
\begin{multline}
L\psi(X)=j \omega   \sum_{l=0}^{d_f}a_{f_l} X^{l} +\frac{1}{2} (j \omega)^2  \sum_{l=0}^{d_g}a_{g_l} X^{l}
+ \\
\psi(X) \int_\mathbb{R}\left(e^{j\omega y} - 1 -
  j \omega y
  \right) \nu(dy).
\end{multline}
Using \eqref{eq:derivativesPsi}, the terms $X^l \psi$ can be replaced by $(j)^{-l}\frac{\partial^l \psi}{\partial \omega^l}$. This yields
\begin{multline}
L\psi(X)=j \omega  \sum_{l=0}^{d_f}a_{f_l} j^{-l} \frac{\partial^l \psi}{\partial \omega^l} +\frac{1}{2} (j \omega)^2  \sum_{l=0}^{d_g}a_{g_l}  j^{-l} \frac{\partial^l \psi}{\partial \omega^l}
+ \\
\psi(X) \int_\mathbb{R}\left(e^{j\omega y} - 1 -
  j \omega y
  \right) \nu(dy).
\end{multline}

Taking expectation, we get a partial differential equation in characteristic function $\varphi=\mathbb{E}(\psi(X))$
\begin{multline}\label{eq:pde}
\frac{\partial \varphi}{\partial t}=j \omega  \sum_{l=0}^{d_f}a_{f_l} j^{-l} \frac{\partial^l \varphi}{\partial \omega^l}
+\\
\frac{1}{2} (j \omega)^2  \sum_{l=0}^{d_g}a_{g_l}  j^{-l} \frac{\partial^l \varphi}{\partial \omega^l}
+
\eta(\omega) \varphi,
\end{multline}
where we used $\eta(\omega)$ to denote the integral $\int_\mathbb{R}\left(e^{j\omega y} - 1 - j \omega y \right) \nu(dy)$.
If the stationary distribution exists (see \cite{met12} for details), then we must have that $\frac{\partial \varphi}{\partial t}=0$. This results in the ordinary differential equation \eqref{eq:StODE}.
The degree of this ODE is $n=\max \{d_f, d_g \}$.
\end{proof}
\begin{remark}
  The function $\eta(\omega)$ is known in as the \emph{characteristic
    exponent}, \cite{applebaum2009levy}. A special case of the
  L\'evy-Khintchine formula shows that the compensated Poisson process
  has characteristic function given by:
\begin{equation}
  \mathbb{E}\left[
    e^{j\omega \int_0^ty \tilde{N}(dt,dy)}
    \right] = e^{t\eta(\omega)}.
  \end{equation}
  Thus, the statistics of the jump process $\int_0^t y
  \tilde{N}(dt,dy)$ are entirely determined by $\eta(\omega)$.

  For commonly-studied L\'evy processes, the characteristic exponent
  has an explicit formula. For example, the gamma process $\tau(t)$ has L\'evy
  measure
  \begin{equation}
    \nu(dy) = ay^{-1} e^{-by} dy,
  \end{equation}
  where $a$ and $b$ are positive parameters. It can be shown that
  \begin{equation}
    \mathbb{E}[e^{j\omega \tau(t)}] =
    e^{-ta\log\left(1-j\frac{\omega}{b}\right)}
    \quad \textrm{ and } \quad \mathbb{E}[\tau(t)] = t\frac{a}{b}
  \end{equation}
  It follows that the compensated Gamma process, $\tau(t)-\mathbb{E}[\tau(t)]$,
  has characteristic exponent given by
  $-a\log\left(1-j\frac{\omega}{b}\right)-j\omega a/b$.

  A more complex example, which will be studied in an example below is
  the variance-gamma process: $\sigma W(\tau(t))$. This process is
  formed by composing a standard Brownian motion $W(t)$ with a gamma
  process. In this case, it can be shown that
  \begin{equation}
    \mathbb{E}[e^{j\omega W(\tau(t))}] = e^{-ta\log\left(
        1+\frac{(\sigma \omega)^2}{2b}
        \right)} \quad \textrm{and} \quad \mathbb{E}[W(\tau(t))] = 0.
  \end{equation}
  It follows that compensated variance-gamma process is simply the
  variance-gamma process. Furthermore, the characteristic exponent is
  given by
  \begin{equation}
    \label{eq:vgExp}
    \eta(\omega) = - a \log\left(
1+\frac{(\sigma\omega)^2}{2b}
\right),
  \end{equation}

\end{remark}
Next, we discuss the solution of the ODE for the stationary characteristic function.

\subsection{Solving for stationary characteristic function}
The ODE in \eqref{eq:StODE} cannot be solved analytically except for a  handful of cases. However, it can be solved via numerical techniques. Either way it would require $n$ initial/intermediate/boundary values to find the solution. Other than the usual $\varphi(0)=1$, previous works have either utilized prior knowledge about the system (e.g., the distribution is symmetric), or used $\lim_{|\omega|\to \infty}\varphi(\omega) = 0$ and $\lim_{|\omega|\to \infty}\frac{\partial^l \varphi(\omega)}{\partial \omega^l} = 0$ for some $l$ \cite{gri13}. In practice these are hard to incorporate in a solution. Furthermore, if one is interested only in stationary moments, then solution of $\varphi(\omega)$ only in neighborhood of zero is sufficient.

We propose a different approach to compute both the moments and the characteristic function. This approach relies on two ideas. First being the fact that the characteristic function is related with the moments as
\begin{equation}\label{eq:momMu}
 \frac{\partial^l \varphi}{\partial \omega^l}\left |_{\omega=0} \right. = j^{l} \mu_l, \; l=\{1, \ldots, n-1\},
 \end{equation}
where $\mu_l \in \mathbb{R}$ represents the $l^{th}$ order moment. Thus, the moments are natural quantities to be used in computing the characteristic function. Second idea is to utilize the Bochner's theorem to estimate the moments $\mu_l$. In particular, we can use $\varphi(0)=1$ and positive semidefinite property of the matrix $M$ in \eqref{eq:M}. Using these, a semidefinite program can be formulated that gives lower and upper bounds on $\mu_l$ as stated in the theorem below.

\begin{theorem}
 Consider the one dimensional system defined by \eqref{eq:JD}. Assuming that $f$ and $g^2$ are polynomials of degree $d_f$ and $d_g$, a lower bound on a moment $\mu_l$ can be obtained via the semidefinite program
\begin{subequations}\label{eq:SDPsetup}
\begin{align}
& \min \;  \mu_k  \label{eq:SDPsetupOBJ}\\
& j \omega  \sum_{l=0}^{d_f}a_{f_l} j^{-l} \frac{\partial^l \varphi}{\partial \omega^l} +\frac{1}{2} (j \omega)^2  \sum_{l=0}^{d_g}a_{g_l}  j^{-l} \frac{\partial^l \varphi}{\partial \omega^l} + \eta(\omega) \varphi=0, \label{eq:SDPsetupODE}\\
& \frac{\partial^l \varphi}{\partial \omega^l}\left |_{\omega=0} \right. = j^{l} \mu_l, \; \; l={1, 2, \ldots n}, \label{eq:SDPsetupMom}\\
& \varphi(0)=1, \label{eq:SDPsetupNorm}\\
& M \succeq 0.\label{eq:SDPsetupM}
\end{align}
\end{subequations}
Here $k=\{1, \ldots, n\}$ with $n=\max \{d_f, d_g\}$ and $M$ is defined as in \eqref{eq:M}. Further, the minimum value obtained by the program increases as size of $M$ is increased by including more test points.
\end{theorem}
\begin{proof}
Since $f$ and $g^2$ are assumed to be polynomials, Theorem~\ref{thm:ode} implies that the characteristic function $\varphi(\omega)$ satisfies the ODE of order $n$ given by \eqref{eq:SDPsetupODE}. The moments are related with the derivatives of the characteristic function by the linear constraints in \eqref{eq:SDPsetupMom}. The constraint in \eqref{eq:SDPsetupNorm} and \eqref{eq:SDPsetupM} are a consequence of the Bochner's theorem. Since the objective function is linear in decision variables $\mu_k$ and the constraints are either equality or semidefinite constraints, the optimization problem is a semidefinite program \cite{boyd2004convex}.

Now suppose that the size of $M$ is increased by including more test points $\omega_1, \ldots, \omega_q$. This corresponds to adding more constraints in the program, and the solution cannot get worse by doing so.
\end{proof}

The upper bound on $\mu_k$ can be found by minimizing $-\mu_k$. Note that we can choose any test points $\omega_1, \omega_2, \ldots, \omega_q \in \mathbb{R}$ in order to generate the matrix $M$. For sake of simplicity, we will choose uniformly spaced values on the real-line.

The above semidefinite program can be used to compute lower and upper bounds on each of the moments $\mu_k$. These values can be then used to determine the solution to the ODE for characteristic function and thereby finding an approximation of the characteristic function. If the lower and upper bounds on each of the moments are reasonably close, then the approximate characteristic function would be quite close to the true characteristic function. It can be further used to compute the stationary probability density via inversion.

\begin{remark}
  \label{rem:OptControl}
The formulation in \eqref{eq:SDPsetup} can be interpreted as an optimal control problem for a linear time varying system. Specifically, consider a state vector $z=\begin{bmatrix} \varphi & \frac{\partial \varphi}{\partial \omega} & \cdots & \frac{\partial^{n-1} \varphi}{\partial \omega^{n-1}}\end{bmatrix}^\top$. Then, the differential equation describing the characteristic function becomes a linear time (in $\omega$-space) varying system
\begin{equation}
\frac{dz}{d\omega}=\mathcal{P}(\omega) z(\omega),
\end{equation}
for an appropriately defined matrix $\mathcal{P}$. In this setup, the objective would be to optimize the elements of $z(0)$ subject to the linear matrix inequality \eqref{eq:SDPsetupM}. Note that the first element of $z(0)$ is given by $1$.

\end{remark}

\begin{remark}
The semidefinite program in \eqref{eq:SDPsetup} can be used to find bounds on first $n$ moments where $n$ is the degree of the ODE \eqref{eq:SDPsetupODE}. If one is interested in computing the higher order moments, the approximate characteristic function can be differentiated and computed at $\omega=0$. By doing so, bounds on the higher order moments can also be computed. Alternatively, for systems with finite moments, one can compute the bounds on first $n$ moments via the proposed method, and then use the fact that stationary moments are related via a linear system of equations given by \cite{lgs17}
\begin{equation}\label{eq:XXbar}
A \mathcal{X} + B \overline{\mathcal{X}}=0.
\end{equation}
Here $\mathcal{X}$ is collection of moments up to some order and $\overline{\mathcal{X}}$ contains moments of order higher than those in $\mathcal{X}$. The number of elements of $\overline{\mathcal{X}}$ is as many as the degree of nonlinearity in the system given by $n$. While the usual moment closure methods estimate elements of $\overline{\mathcal{X}}$ in terms of those of $\mathcal{X}$, we simply supplement \eqref{eq:XXbar} with lower and upper bounds on $n$ moments and thereby compute bounds on all other moments in \eqref{eq:XXbar}.
\end{remark}
In the next section, we illustrate the proposed method on some simple examples and verify its performance.

\section{Examples}\label{sec:examples}
\begin{example}[Stochastic Logistic Model]
\begin{figure}
\centering
\includegraphics[width=0.95\linewidth]{fig1.pdf}\\
\caption{Comparison of estimated mean and characteristic function with exact values. {\it Top}: The mean $\mu_1$ is estimated via semidefinite program in \eqref{eq:SDPsetup}. The matrix $M$ is constructed by taking frequencies spaced by $one$ (e.g., $\omega_1=1, \omega_2=2, \ldots$). Increasing the number of test frequencies leads to a better estimation of the moment as both lower and upper bounds converge.Estimation of mean via SDP. {\it Bottom}: The characteristic function is estimated by using the middle value of lower and upper bounds for $\mu_1$ obtained from $30$ test frequencies in \eqref{eq:phi}--\eqref{eqn:cfCoeff}. The estimated probability densify function is obtained by numerically inverting the characteristic function. Notably, this particular system is exactly solvable, and the exact probability densify function is provided for comparison purpose, showing an excellent match with its estimated counterpart.}
\label{fig:mu1Logistic}
\end{figure}

Consider the following modified stochastic logistic growth model
\begin{equation}
dX = (1+X-0.1 X^2)dt + \sqrt{2} X dW.
\end{equation}
Without the constant term in the drift, this model is widely-used in modeling population growth \cite{allen07}. We added the constant term so that the trivial solution $X=0$ is ruled out.

Using \eqref{eq:pde}, the characteristic function evolves as per
\begin{align}
\frac{\partial \varphi}{\partial t}
=
j \omega \varphi + \omega \frac{\partial \varphi}{\partial \omega}+j 0.1 \omega \frac{\partial^2 \varphi}{\partial \omega^2} + \omega^2 \frac{\partial^2 \varphi}{\partial \omega^2}.
\end{align}
Therefore, the stationary characteristic function is the solution to the following differential equation
\begin{align}
j \omega \varphi + \omega \frac{\partial \varphi}{\partial \omega}+\left(j 0.1 \omega  + \omega^2\right) \frac{\partial^2 \varphi}{\partial \omega^2} =0.
\end{align}

It can be shown that the above ODE has the following generalized solution
\begin{multline}\label{eq:phi}
\varphi (\omega )
=
\frac{c_1}{\sqrt{5}} I_{0}\left(2 \sqrt{0.1-j \omega} \right)+\frac{c_2}{\sqrt{5}} K_{0}\left(2 \sqrt{0.1-j \omega}\right),
\end{multline}
where $I$ and $K$ denote the modified Bessel functions of first and second kinds, and $c_1$, $c_2$ are unknown coefficients. As expected,  the number of unknown coefficients is same as the order of nonlinearity in the dynamics.
\end{example}

To determine the coefficients, we can use the fact that $\varphi(0)=1$ and $\frac{\partial \varphi}{\partial \omega}\left |_{\omega=0} \right.=j \mu_1$, where $\mu_1 \in \mathbb{R}$ is the mean that is to be determined. This results in
\begin{subequations}\label{eqn:cfCoeff}
\begin{align}
& \frac{c_1}{\sqrt{5}} I_{0}\left(2 \sqrt{0.1} \right)+\frac{c_2}{\sqrt{5}} K_{0}\left(2 \sqrt{0.1}\right)=1 \\
& -j \sqrt{2} c_1 I_1\left(2\sqrt{0.1}\right)+j \sqrt{2} c_2 K_1\left(2\sqrt{0.1}\right)=j \mu_1
\end{align}
\end{subequations}
Using these, the stationary characteristic function can be written in terms of only one unknown $\mu_1$. Now, we can compute bounds on $\mu_1$ using the semidefinite program as in \eqref{eq:SDPsetup}. By choosing uniformly spaced values of $\omega_1, \omega_2, \ldots, \omega_n$, we computed the maximum and minimum allowable values $\mu_1$. We also find that increasing the size of the program by choosing more test points improves both lower and upper bounds. Taking $30$ test points at $\omega_1=1, \ldots \omega_{30}=30$, we get $5.2024 \leq \mu_1 \leq 5.2025$ (see Fig.~\ref{fig:mu1Logistic}, Top). This result is in excellent agreement with Monte Carlo simulations which yield an value of $5.2$ for $10000$ simulations.

Using the value of $\mu_1$ obtained here, we can use the characteristic function to reconstruct the probability density function of the stationary distribution (see Fig.~\ref{fig:mu1Logistic}, Bottom). As mentioned earlier, the bounds on $\mu_1$ can be used to estimate bounds on higher order moments as well (results not shown here).

\begin{example}[Variance Gamma Process]
\begin{figure}
\includegraphics[width=\linewidth]{Levy.pdf}
\caption{Comparison of probability density functions obtained from the proposed method and simulations. The histogram is obtained via $1000$ simulations of the process, whereas the estimated probability density is inverted from the estimated characteristic function by solving the semidefinite program in \eqref{eq:SDPsetup}. The optimization was performed by casting the problem as an optimal control problem described in Remark 2, with evenly spaced $80$ test frequencies.}
\label{fig:levy}
\end{figure}

Consider the following process
\begin{equation}
  dX=(-5X^3+10X)dt+dW(\tau(t)),
\end{equation}
where $W(\tau(t))$ is the variance-gamma process described
above. Recall, in particular, $W(\tau(t))$ has zero mean, and so the
compensated variance-gamma process is the same as the original
variance-gamma process. Furthermore, recall the expression for the
characteristic exponent, $\eta(\omega)$, from \eqref{eq:vgExp}.

Using \eqref{eq:pde}, the characteristic function of this process is governed by
\begin{equation}
\frac{\partial \varphi}{\partial t}=5 \omega^3 \frac{\partial^3 \varphi}{\partial \omega^3}+10 \omega \frac{\partial \varphi}{\partial \omega}+\eta(\omega)\varphi
\end{equation}
As before, this ODE solution will have three unknown coefficients, out
of which two could be readily computed by using the facts that
$\varphi(0)=1$ and the stationary mean value of $X$ is zero due to
symmetry of the process. Numerically solving the optimization problem
that minimizes $\mu_2$ (the second moment), and inverting the solution
yields quite accurate estimate of the stationary probability density
function (see Fig.~\ref{fig:levy}).

The optimization was performed by discretizing the optimal control
problem described in Remark~\ref{rem:OptControl} via the trapezoidal
rule, and then directly solving the corresponding SDP. Here, the
decision variables were the state variables $z(\omega_i)$, where
$\omega_i$ are the discretization frequencies. One subtlety is that in
order to express the matrix $M$ in terms of decision variables, the
frequencies, $\omega_i$, must be evenly spaced.
\end{example}

\section{Conclusion}\label{sec:discussion}
In this paper, we proposed a method to estimate stationary characteristic function for a class of stochastic dynamical system driven by both white and L\'evy noise.  The characteristic function for these stochastic systems is governed by an ODE. The method relies upon restricting the solutions of the ODE to positive definite functions and casts the problem as a semidefinite program. In future work, we would extend this framework to multidimensional cases. It would also be interesting to explore whether transient solutions of the characteristic function (described via PDE) can be obtained via a similar method.

\section*{Acknowledgment}
AS is supported by the National Science Foundation Grant ECCS-1711548.


\section{Introduction}

What are cities? In {\it The Death and Life of the Great American Cities},
Jacobs argues that human relations can be seen as a proxy for places within
cities \cite{jacobs1961}. A modern view of cities establishes that they can be
defined by the interactions among several types of networks \cite{batty2013,
west2017}, from infrastructure networks to social networks. In recent years, an
increasing number of studies have been proposed to define cities through
consistent mathematical models \cite{makse1995, makse1998, rozenfeld2008,
rozenfeld2011, rybski2013, murcio2013, hernando2014, frasco2014, arcaute2014,
masucci2015, arcaute2016, fluschnik2016} and to investigate urban indicators at
inter- and intra-city scales, in order to shed some light on problems faced by
decision makers \cite{bettencourt2007, arbesman2009, bettencourt2010a,
bettencourt2010b, gomez-lievano2012, gallos2012a, bettencourt2013, fragkias2013,
oliveira2014, melo2014, louf2014, alves2015, li2015, hanley2016, caminha2017,
operti2017}. Despite the efforts of such studies, to properly defining the
boundaries of urban settlements remains an open problem in the Science of
Cities. A minimum criterion of acceptability for any model of cities seems to be
the one that retrieves a conspicuous scaling law found for United States (US),
United Kingdom (UK), and other countries, known as Zipf's law
\cite{rozenfeld2008, rozenfeld2011, gabaix1999a, gabaix1999b, ioannides2003,
cordoba2008, giesen2010, giesen2011, jiang2011, giesen2012, ioannides2013,
duraton2014, watanabe2015}.

In 1949, Zipf \cite{zipf1949} observed that the frequency of words used in the
English language obeys a natural and robust power law behavior, {\it i.e.} a few
words are used many times, while many words are used just a few times. Zipf's
law can be represented generically by the following relation between the size
$S$ of objects from a given set and its rank $R$:
\begin{equation}\label{eq1}
R \propto S^{-\zeta},
\end{equation}

\noindent
where $\zeta=1$ is Zipf's exponent. The size of objects is, in the original
context, the frequency of used words. On the other hand, if such objects are
cities, then the sizes stand for the population of each city, taking into
account Zipf's law and reflecting the fact that there are more small towns than
metropolises in the world. We emphasize that it is not straightforward that the
Zipf's law, despite its robustness, should hold independently of the city
definition, since other scaling relations are not, such as the allometric
exponents for CO$_2$ emissions and light pollution \cite{oliveira2014,
operti2017}. Many other man-made and natural phenomena also exhibit the same
persistent result, {\it e.g.} earthquakes and incomes \cite{sornette1996,
okuyama1999}.

Here, we propose a worldwide model to define urban settlements beyond their
usual administrative boundaries through a bottom-up approach that takes into
account cultural, political, and geographical biases naturally embedded in the
population distribution of continental areas. After all, it is not surprising
that two regions, {\it e.g.} one in western Europe and another one in eastern
Asia, spatially contiguous in population or in commuting level have different
cultural, political or geographical characteristics. Thus, it is also not
surprising that such issues yield different stages of the same mechanics of
growth. The main goal of our model is to be successful in defining cities even
in large regions. Our conjecture is straightforward: there are hierarchical
mechanisms, similar to those present in previous studies of cities in the UK
\cite{arcaute2016} and brain networks \cite{gallos2012b}, behind the growth and
innovation of urban settlements. These mechanisms are ruled by a combination of
general measures, such as the population and the area of each city, and
intrinsic factors which are specific to each region, {\it e.g.} topographical
heterogeneity, political and economic issues, and cultural customs and
traditions. In other words, if political turmoil or economic recession plagues a
metropolis for a long time, all of its satellites are affected too, {\it i.e.}
the entire region ruled by the metropolis will be negatively impacted.

\section{The models}

\subsection{City Clustering Algorithm (CCA)}

In 2008, Rozenfeld {\it et al.} \cite{rozenfeld2008} proposed a model to define
cities beyond their usual administrative boundaries using a notion of spatial
continuity of urban settlements, called the City Clustering Algorithm (CCA)
\cite{rozenfeld2008, rozenfeld2011, rybski2013, oliveira2014, frasco2014,
fluschnik2016, caminha2017, operti2017}. The CCA is defined for discrete or
continuous landscapes \cite{rozenfeld2011} by two parameters: a population
density threshold $D^*$ and a distance threshold $\ell$. These parameters
describe the populated areas and the commuting distance between areas,
respectively. Here, we adopt the following strategy to improve the discrete CCA
performance: (i) Supposing a regular rectangular lattice $L_x \times L_y$ of
sites where the population density of the $k$-th site is $D_k$, we perform an
initial agglomeration by $D^*$ to identify all clusters. If $D_k>D^*$, then the
$k$-th site is populated and we aggregate it with its populated nearest
neighbors. Otherwise, the $k$-th site is unpopulated. (ii) For each populated
cluster, we define its {\it shell sites}, {\it i.e.} sites in the interface
between populated and unpopulated areas. (iii) Lastly, we perform a final
agglomeration by $\ell$, taking into account only the shell sites. If
$d_{ij}<\ell$, where $d_{ij}$ is the distance between the $i$-th and $j$-th
shell sites, and if they belong to different clusters, then the $i$-th and
$j$-th sites belong to the same CCA cluster, even with spatial discontinuity.
Otherwise, they indeed belong to different CCA clusters. This simple strategy
improves the algorithm's computational performance because the number of shell
sites is proportional to $L$, where $L=L_x \approx L_y$ is a linear measure of
the lattice.

\subsection{City Local Clustering Algorithm (CLCA)}

We propose a worldwide model based on the CCA, called the City Local Clustering
Algorithm (CLCA), not only to define cities beyond their usual administrative
boundaries, but also to take into account the intrinsic cultural, political and
geographical biases associated with most societies and reflected in their
particular growing dynamics. The traditional CCA, with fixed $\ell$ and $D^*$,
when applied to a large population density map, can introduce biases defining a
lot of clusters in some regions, while in others just a few. We present the CLCA
with the aim of defining cities even in large regions in order to overcome such
CCA weakness. Hence, it is possible that other models, such as the models based
on street networks proposed by Masucci {\it et al.} \cite{masucci2015} and
Arcaute {\it et al.} \cite{arcaute2016}, carry the same CCA burden and that
local adaptations are necessary for their applications into large regions.

The main idea of our model is to analyze the change of the CCA clusters
through the variation of $D^*$ under the perspective of different regions.
First, we define a regular rectangular lattice $L_x \times L_y$ of sites, where
the population density of the $k$-th site is $D_k$. We sort all the sites in a
list according to the population density, in descending order. Therefore, the
site with the greatest population density is the first entry in this list, which
we call the first {\it reference site}. The reference site can be considered as
the current core of the analyzed region. Second, we apply the CCA to the
lattice, keeping a fixed value of $\ell$, for a range of $D^*$ decreasing from a
maximum value $D^{(max)}$ to a minimum value $D^{(min)}$ with a decrement
$\delta$. During the decreasing of $D^*$, clusters are formed and they spread
out to all regions of the lattice. Eventually, the cluster that contains the
reference site (from now on the {\it reference cluster}), together with one or
more of the other clusters, will merge from $D^{(i)}$ to $D^{(i+1)}$, where
$D^{(i+1)} = D^{(i)}-\delta$. In order to accept or deny the merging of these
clusters, we introduce three conditions:

\begin{enumerate}[(i)]
\item If the area $A_r(D^{(i)})$ of the reference cluster $r$, {\it i.e.} the
cluster that contains the $r$-th reference site at $D^{(i)}$, obeys
\begin{equation}\label{eq2}
A_r(D^{(i)}) < A^*,
\end{equation}

\noindent
then the reference cluster $r$ always merges with other clusters, because it is
still considered very small. In this context, the area $A^*$ can be understood
as the minimal area of a metropolis.

\item If the difference between the areas of the reference cluster $r$ at
$D^{(i+1)}$ and $D^{(i)}$ obeys
\begin{equation}\label{eq3}
A_r(D^{(i+1)})-A_r(D^{(i)}) > H^* A_r(D^{(i)}),
\end{equation}

\noindent
then the reference cluster $r$ has grown without merging (Fig.\ref{clca}a) or
there is a merging of at least two large clusters (Fig.\ref{clca}b). In the last
case, we emphasize that if there are more than two clusters involved in the
merging process, the reference cluster $r$ may not be one of the largest. As the
first case is not desirable, we can avoid it by reducing the value of $\delta$
and keeping the value of $H^*$ relatively high. The parameter $H^*$ can be
understood as the percentage of the area of the reference cluster $r$ at
$D^{(i)}$. If the second case happens, we consider the entire region inside of
the reference cluster $r$ at $D^{(i+1)}$, but the clusters of this region (which
we call the {\it usual clusters}) are defined by those at $D^{(i)}$. The usual
clusters are the CCA clusters at the imminence of the merging process between
$D^{(i)}$ and $D^{(i+1)}$. This includes the reference cluster $r$ itself and
one or more of the other clusters before the merging (Fig.\ref{clca}b).
Furthermore, all of the sites of the reference cluster $r$ at $D^{(i+1)}$ are
removed from the initial list of reference sites. This condition is necessary
because we should not merge two large metropolises.

\item In condition (ii), when a reference cluster $r$ is merging with another
cluster that covers one or more regions already defined by previous reference
clusters at different values of $D^*$, there is a strong likelihood of the
emergence of a {\it forbidden region} within that cluster. In this case, we
force the region already defined by the largest value of $D^*$ to grow to the
limits of the forbidden region (Fig.\ref{clca}c). The forbidden regions are the
complementary areas of the reference clusters already defined within the usual
clusters. As a consequence of this procedure, some CCA clusters that were hidden
after the analysis of the previous reference cluster arise in this forbidden
region. We justify this condition by the idea that a metropolis rules the growth
of its satellites, since it plays a fundamental role in their socioeconomic
relations.
\end{enumerate}

\noindent
We apply the same procedure to the second reference cluster, to the third
reference cluster, and so on. Finally, we also define the {\it isolated
clusters} with the minimum value of $D^*$ for all the cases accepted in
condition (ii). In order to make our model clearer, we chose the descending
order to sort the population density for one reason: To favor the merging
process of the high density clusters that rose from the decreasing of the $D^*$.
In practice, we run our revised discrete CCA just once for the entire range of
input parameters and store all of the outputs in order to improve the
performance of the model. The apparent simplicity of this task hides a RAM
memory management problem of storing all of the outputs in a medium-performance
computer. We overcome such a barrier through the {\it zram module}
\cite{zram2017}, available in the newest linux kernels. The zram module creates
blocks which compress and store information dynamically in the RAM memory
itself, at the cost of processing time.

\section{The dataset}

We use the {\it Global Rural-Urban Mapping Project} (GRUMPv1) \cite{grump2011},
available from the {\it Socioeconomic Data and Applications Center} (SEDAC) at
Columbia University, to apply the CLCA to a single global dataset. The GRUMPv1
dataset is composed of georeferenced rectangular population grids for 232
countries around the world in the year 2000 (Fig. \ref{grumpv1}). Such a dataset
is a compilation of gridded census and satellite data for the populations of
urban and rural areas. These data are provided at a high resolution of
$30\;\text{arc-seconds}$, equivalent to $30/3600\;\text{degrees}$ or a grid of
$0.926\;km \times 0.926\;km$ at the Equator. We note that despite of the
heterogeneous population distributions that built the GRUMPv1, its overall
resolution is tolerable to the CLCA, since we can identify well-defined clusters
around all continents in the raw data.

We calculate the area of each site by the composition of two spherical
triangles \cite{snyder1987}. The area of a spherical triangle with edges $a$,
$b$ and $c$ is given by
\begin{equation}\label{eq4}
A=4 R_e^2 \tan^{-1}\left[\tan\left(\frac{s}{2}\right)\tan\left(\frac{s_a}{2}\right)\tan\left(\frac{s_b}{2}\right)\tan\left(\frac{s_c}{2}\right)\right]^{1/2},
\end{equation}

\noindent
where $s=(a/R_e+b/R_e+c/R_e)/2$, $s_a=s-a/R_e$, $s_b=s-b/R_e$, and
$s_c=s-c/R_e$. In this formalism, $R_e=6,378.137\;km$ is the Earth's radius and
the edge lengths are calculated by the great circle (geodesic) distance between
two points $i$ and $j$ on the Earth's surface:
\begin{equation}\label{eq5}
d_{ij}=R_e\cos^{-1}[\sin(\phi_i) \sin(\phi_j)+\cos(\phi_i)\cos(\phi_j)\cos(\lambda_j-\lambda_i)].
\end{equation}

\noindent
The values of $\lambda_i$ ($\lambda_j$) and $\phi_i$ ($\phi_j)$, measured in
radians, are the longitude and latitude, respectively, of the point $i$ ($j$).
Thus, we are able to define the population density for each site of the lattice,
since its population and area are known.

We also pre-process the GRUMPv1 dataset, dividing all countries and continents
--- and even the entire world --- into large regions which we call {\it clusters
of regions}, to apply our model in a feasible computational time using
medium-performance computers. These regions are defined by the CCA with lower
and upper bound parameters $D^*=50\;\text{people}/km^2$ and $\ell=10\;km$,
respectively. We believe that such large clusters can hold the socioeconomic and
cultural relations among different urban settlements of a territory. Fig.
\ref{us_east}a shows the largest clusters of regions in the US; as we can see,
all of the Eastern US is considered a single cluster.

\section{Results}

To show the relevance of our model, we apply the CLCA to the GRUMPv1 dataset at
three different geographic levels: countries, continents and the entire world.
For each case, we consider only a single set of CLCA parameters. We justify our
choices with the following assumptions: (i) $D^{(min)}=100\;\text{people}/km^2$,
a value slightly greater than the lower bound CCA parameter
($D^*=50\;\text{people}/km^2$) used to define the regions of clusters; (ii)
$D^{(max)}=1000\;\text{people}/km^2$, a loosened value of $D^{(max)}=\infty$;
(iii) $\delta=10\;\text{people}/km^2$, a small enough value to avoid the
reference clusters growing without merging; (iv) $\ell=3\;km$, the critical
distance threshold, already extensively analyzed by previous CCA studies
\cite{rozenfeld2008, rozenfeld2011, oliveira2014}; (v) $A^*=50\;km^2$, the
minimum area of a metropolis, as it is required that $A^*$ be reasonably greater
than the minimum unit of area from the dataset and smaller than a metropolis'
area; and (vi) $H^*=0.05$, a large enough value to favor the merging of clusters
which are similar in size. The Fig. \ref{us_east}b shows the CLCA cities defined
by the single set of CLCA parameters. For other regions, see the Supplementary
Information (SI).

We study the population distribution using the Maximum Likelihood Estimator
(MLE) proposed by Clauset {\it et al.} \cite{clauset2009}. Their approach
combines maximum-likelihood fitting methods with goodness-of-fit tests based on
Kolmogorov-Smirnov (KS) statistic. The Fig. \ref{g_countries_panel} shows the
log-log behavior of the Cumulative Distribution Function (CDF) for the
population of the CLCA cities, considering only the countries with the highest
number of CLCA cities for each continent (for other countries, see the SI). The
$Pr(\mathcal{P} \geq P)$ represents the probability that a random population
$\mathcal{P}$ takes on a value greater or equal to the population $P$. In all
CDF plots, we also show the maximum likelihood power-law fit, as well as, the
value of the exponent $\zeta=\alpha-1$, where $\alpha$ is the MLE exponent, and
the value of $P_{min}$, the lower bound of the MLE.

In Fig. \ref{g_countries_hist}, we show a normalized histogram, with frequency
$F$, of the $\zeta$ exponents for all countries ( 145 out of 232) with at least
10 CLCA cities in the region covered by the maximum likelihood power-law fit.
The mean value of the $\zeta$ exponents is $\bar{\zeta}=0.98$, with variance
$\sigma^2=0.09$. The dashed red line stands for the normal distribution
$\mathcal{N}(\bar{\zeta},\sigma^2)$. In spite of the $\zeta$ exponent
heterogeneity illustrated by Fig. \ref{g_countries_hist}, Zipf's law holds for
most countries around the globe. We emphasize that such results corroborate with
previous studies performed for one country or a small number of countries
\cite{rozenfeld2008, rozenfeld2011, gabaix1999a, gabaix1999b, ioannides2003,
cordoba2008, giesen2010, giesen2011, jiang2011, giesen2012, ioannides2013,
duraton2014, watanabe2015}. In special, the Fig. \ref{g_countries_hist} also
endorses an astute meta-analysis performed by Cottineau \cite{cottineau2017}.
Cottineau provided a comparison among the Zipf's law exponents found in 86
studies. Our results strongly corroborate those presented in such study, except
that our exponents are ranged between 0 and 2.

Furthermore, we challenge the robustness of our model at higher geographic
levels: continents and the entire world. We performed the same analyses and find
that our results persist on both scales, {\it i.e.} the CLCA cities follows
Zipf's law for continents and the entire world, as illustrated in Figs.
\ref{g_continents_panel} and \ref{g_world}.

We summarize our results in a set of 7 tables: Tables \ref{t_countries_af},
\ref{t_countries_as}, \ref{t_countries_eu}, \ref{t_countries_na} ,
\ref{t_countries_oc}, and \ref{t_countries_sa}, for countries from Africa; Asia;
Europe; North America; Oceania; and South America, respectively. Table
\ref{t_continents_world} contains similar information for all continents and the
entire world. In all cases, we show the name of the considered region (country,
continent or globe), the ISO 3166-1 alpha-3 code associated (only for
countries), the number of cities obtained by the CLCA and those covered by the
MLE, the lower bound $P_{min}$, and the Zipf's exponent $\zeta$.

It is remarkable that the top CLCA city, with a population of 63,585,039 people,
is composed of three large urban settlements (Alexandria, Cairo, and Luxor)
connected by several small ones. Figs. \ref{cairo}a, \ref{cairo}b, and
\ref{cairo}c show the largest cluster of regions in Egypt for the GRUMPv1
dataset; CLCA cities; and night-time lights from the National Aeronautics and
Space Administration (NASA) \cite{nasa2016}, respectively. We believe the main
reason for this finding has been present in the Northeast of Africa since before
the beginning of ancient civilization --- namely, the Nile river. Actually, it
is well known that almost the entire Egypt population lives in a strip along the
Nile river, in the Nile delta, and in the Suez canal on 4\% of the total country
area ($~10^6\;km^2$), where there are arable lands to produce food
\cite{gehad2003}. The river and delta regions are composed by some large cities
and a lot of small villages, making them extremely dense. Therefore, our results
raise the hypothesis that the cities and villages across the Nile can be seen as
a kind of ``megacity'', despite spatially non-contiguous, due to the
socioeconomic relation, reflected on the high commuting levels, among close
subregions.

The Table \ref{t_world} shows the top 25 CLCA cities in the entire world by
population, and their associated areas. After the top CLCA city,
Alexandria-Cairo-Luxor, we emphasize that the 13 next-largest CLCA cities are in
Asia. Indeed, we can see that the shape of the tail end of the entire world
population distribution (in Fig. \ref{g_world}) is roughly ruled by the greater
CLCA city in Africa and several CLCA cities in Asia.

These facts are not in line with what was recently reported by the {\bf United Nations (UN)} \cite{un2016}, {\it e.g.} the largest CLCA city, Alexandria-Cairo-Luxor, is just the 9th largest city according to the UN, and the largest UN city, Tokyo, is just the 4th largest according to our analyses.

\section{Conclusions}

We propose a model to define urban settlements through a bottom-up approach
beyond their usual administrative boundaries, and moreover to account for the
intrinsic cultural, political, and geographical biases associated with most
societies and reflected in their particular growing dynamics. We claim that such
a property qualifies our model to be applied worldwide, without any regional
restrictions. We also propose an alternative strategy to improve the
computational performance of the discrete CCA. We emphasize that the CCA can
still be used to define cities; however, it depends upon a different tuning of
its parameters for each large region without direct socioeconomic and political
relations. Furthermore, we show that the definition of cities proposed by our
approach is robust and holds to one of the most famous results in Social
Science, Zipf's law, not only for previously studied countries, {\it e.g.} the
US, the UK, or China, but for all countries (145 from 232 provided by GRUMPv1)
around the world. We also find that Zipf's law emerges at different geographic
levels, such as continents and the entire world. Another highlight of our study
is the fact that our model is applied upon one single dataset to define all
cities. Furthermore, we find that the most populated cities are not the major
players in the global economy (such as New York City, London, or Tokyo). The
largest CLCA city, with a population of 63,585,039 people, is an agglomeration
of several small cities close to each other which connects three large cities:
Alexandria, Cairo, and Luxor. Finally, after the top CLCA city of
Alexandria-Cairo-Luxor, we find that the next-largest 13 CLCA cities are in
Asia. These facts are not in full agreement with a recent UN report
\cite{un2016}. According to our results, the largest CLCA city,
Alexandria-Cairo-Luxor, is just the 9th largest city according to the UN, while
the largest UN city, Tokyo, is just the 4th largest according to our analyses.

\section{Data Accessibility}
The dataset supporting this article are available at
\url{http://sedac.ciesin.columbia.edu/data/collection/grump-v1}. More
specifically, the reader can click on ``Data sets'' and, after that, on
``Population Count Grid, v1 (1990,1995,2000)''. We also provide the codes for
the proposed model that are available at
\url{https://doi.org/10.5061/dryad.968nq8n}.

\section{Competing Interests}
We have no competing interests.

\section{Author's Contributions}
EAO performed the data analysis, the algorithm of the proposed model, and the
statistical analysis. He also participated in the design of the study and
drafted the manuscript. VF carried out the funding acquisition and helped draft
the manuscript. JSA participated in the design of the study, carried out the
funding acquisition, and helped draft the manuscript. HAM conceived, designed,
and coordinated the research, as well as, carried out the funding acquisition
and helped draft the manuscript. All authors approved the manuscript.

\section{Acknowledgements}
We thank the Global Rural-Urban Mapping Project (GRUMPv1) team for the dataset
provided. Furthermore, we would like to thank X. Gabaix for helpful comments and
discussions.

\section{Funding}
We gratefully acknowledge funding by CNPq, CAPES, FUNCAP, NSF, and ARL
Cooperative Agreement Number W911NF-09-2-0053 (the ARL Network Science CTA).