We study the time evolution of the mutual information between the
mass distributions in spatially separated but casually connected
regions in an expanding universe. The evolution of the mutual
information is primarily determined by the configuration entropy rate
which depends on the dynamics of the expansion and the growth of the
density perturbations. The joint entropy between the distributions
from the two regions plays a negligible role in such evolution. The
mutual information decreases with time in a matter dominated Universe
whereas it stays constant in a $\Lambda$-dominated Universe. The
$\Lambda$CDM model and some other models of dark energy predict a
minimum in the mutual information beyond which the dark energy
dominates the dynamics of the Universe. The mutual information may
have deeper connections to the dark energy and the accelerated
expansion of the universe.

§ INTRODUCTION

Entropy plays a key role in understanding a wide range of phenomena in
science. It plays an important role in deciding the evolution of the
universe. The Universe is regarded as a dynamical system in
Cosmology. The dynamics of the expansion of the universe is described
by the Friedmann equation which is based on the Einstein's equation of
General relativity and the cosmological principle. The cosmological
principle assumes that the universe is statistically homogeneous and
isotropic on sufficiently large scales. The validity of this
assumption is crucial to our understanding of the modern cosmology. A
large number of studies have been carried out to verify the
cosmological principle. Studies based on various cosmological
observations find that the universe is statistically homogeneous and
isotropic on scales somewhere beyond $70-150$ Mpc
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13].
The universe is highly inhomogeneous and anisotropic on smaller scales
due to the presence of a clear hierarchy of structures starting from
planets, stars, galaxies, groups, clusters to superclusters. All these
structures in the present universe are believed to emerge from the
growth of the primordial density fluctuations seeded in the early
universe. The observed CMBR temperature fluctuations of $\frac{\Delta
T}{T} \sim 10^{-5}$ at redshift $z\sim 1100$ suggests that the
universe was highly homogeneous and isotropic in the past. These tiny
fluctuations are amplified by the gravitational instability to produce
structures over a wide range of length scales.

Our universe is known to be expanding. Recent observations suggest
that the universe is currently undergoing an accelerated expansion
[14, 15].  Understanding the present accelerated
expansion of the universe is a major unsolved problem in
cosmology. The dynamics of the expansion affect the growth of
inhomogeneities in the universe. Conversely, the inhomogeneities may
also play an important role in the observed acceleration through their
effect on the large-scale dynamics of the universe [16, 17, 18, 19, 20]. [21] suggests that the
observed acceleration of the universe is consistent with the second
law of thermodynamics and the entropy of the universe in the
$\Lambda$CDM model tends to a finite maximum. Interestingly, the
alternative models such as non-singular bouncing Universes, modified
gravity theories and phantom fields do not lead to a state of maximum
entropy [22, 23, 24]. The continuous
dissipation of the information entropy of the matter distribution due
to the gravitational instability may also drive the accelerated
expansion of the universe [25, 26]. These studies
suggest that the large-scale inhomogeneity and the maximum entropy
production principle (MEPP) [27, 28] may
have important roles in the observed acceleration of the universe.

Numerous works in the literature point to the existence of very
large-scale structures in the universe. The `Sloan Great Wall' in the
nearby universe extends to scales greater than 400 Mpc [29].
[30] find a large quasar group that extends to
$500\,h^{-1}$ Mpc at $z \sim 1.3$.  More recently, [31]
report the discovery of an enormously large giant arc at $z\sim 0.8$
that spans $\sim 1$ Gpc. [32] find correlated orientations
of the axes of large quasar groups on Gpc scales. A supervoid of
diameter $\sim 600 \, h^{-1}$ Mpc detected by [33],
indicates the possible existence of very large underdense regions in
the universe. Some other studies report the evidence of bulk flow from
the analysis of Type Ia supernovae [34] and the existence
of anomalously large dipole in the distribution of quasars
[35]. These findings indicate the presence large-scale
inhomogeneity and anisotropy that are in apparent contradiction with
the cosmological principle. However, there can always be homogeneity
and isotropy on some larger scales. So it is difficult to falsify the
cosmological principle solely based on these observations. Further,
the statistical significance of these structures are questionable
[11, 36, 37]. In any case, these observations
are interesting in their own right and require further scrutiny to
arrive at a conclusion.

A wide variety of statistical tools are used to characterise the
inhomogeneities in the universe. The n-point hierarchy of the
correlation functions and their Fourier transforms, the polyspectra
[39] are widely used to study the inhomogeneous matter
distribution in cosmology. The Minkowski Functionals can measure the
morphology of the large-scale structures in the universe
[40]. The Kullback–Leibler Relative Information Entropy can
distinguish the local inhomogeneous mass density field from its
spatial average [41, 42]. [43] show
that a non-negligible dynamical entanglement may arise due to the the
mutual information between spatially separated but causally connected
regions. [44] use the Tsallis relative entropy to
calculate the mutual information between spatially separated but
causally connected regions of the universe. [46] study the
Renyi mutual information between distant spatial regions in the vacuum
state of a free scalar field.  [45] show that the mutual
information between two spatial regions may get enhanced due to
inflation. The mutual information of disjoint regions in higher
dimension is discussed in [47].

In the present work, we want to calculate the Shannon mutual
information between the disjoint but causally connected regions in an
expanding universe. The spatial distributions of matter in any two
distant regions may have finite mutual information due to the presence
of large-scale structures and long-range correlations. We do not
consider an inhomogeneous universe for our current analysis. We
consider a homogeneous and isotropic universe and study the time
evolution of the mutual information between any two distant
regions. It would be interesting to investigate the role of different
constituents of the universe on the time evolution of the mutual
information between distant regions.

The plan of our work is as follows. We describe the time evolution of
the mutual information in an expanding universe in section 2, describe
the results in section 3 and present our conclusions in section 4.

§ MUTUAL INFORMATION AND ITS TIME EVOLUTION

The configuration entropy associated with the matter distribution over
a significantly large volume $V$ of the universe is defined as

\begin{eqnarray}
S(t)=-\int \rho(\vec{r},t)\,\log \rho(\vec{r}, t)\,dV
\label{eq:sc}
\end{eqnarray}

The volume $V$ is subdivided into a number of subvolumes $dV$ and the
density $\rho(\vec{r},t)$ is measured within each of them. Here $r$
describes the comoving coordinate associated with the centre of the
subvolumes and $\rho(\vec{r},t)$ refers to the matter density in the
subvolumes. The density $\rho(\vec{r},t)$ is directly related with the
probability of finding a mass element within a subvolume.

Let us consider two large identical volumes $V$ separated by a large
distance. We label these two volumes as $A$ and $B$. The two regions
$A$ and $B$ are causally connected. The configuration entropy of the
two regions $A$ and $B$ can be written as,

\begin{eqnarray}
\rho_{{\scaleto{A}{3.5pt}}}(\vec{r}, t)\,\log
\rho_{{\scaleto{A}{3.5pt}}}(\vec{r}, t)\,dV_{{\scaleto{A}{3.5pt}}}
\label{eq:sca}
\end{eqnarray}

\begin{eqnarray}
\rho_{{\scaleto{B}{3.5pt}}}(\vec{r}, t)\,\log
\rho_{{\scaleto{B}{3.5pt}}}(\vec{r}, t)\,dV{{\scaleto{B}{3.5pt}}}
\label{eq:scb}
\end{eqnarray}

We consider $A$ and $B$ to be significantly large volumes so that the
universe can be treated as statistically homogeneous and isotropic on
those scales. We are only interested on scales where one can safely
use the linear perturbation theory to describe the evolution of the
configuration entropy in these regions.

This figure shows two large identical volumes $A$ and $B$
divided into equal number of subvolumes. Here
$\rho_{{\scaleto{A}{3.5pt}}}$ and $\rho_{{\scaleto{B}{3.5pt}}}$
refers to the density within any two subvolumes at a given instant
$t$ and $r_{{\scaleto{AB}{3.5pt}}}$ is the radial separation between
the two subvolumes under consideration. We consider $A$ and $B$ to
be causally connected.

One can define the mutual information between the mass distributions
within the two regions $A$ and $B$ as,

\begin{eqnarray}
I_{{\scaleto{AB}{3.5pt}}}(t)=\int \int
\rho_{{\scaleto{AB}{3.5pt}}}(\vec{r_{1}}, t\,;\vec{r_{2}}, t)\,\log
\frac{\rho_{{\scaleto{AB}{3.5pt}}}(\vec{r_{1}}, t\,;\vec{r_{2}},
\label{eq:miab}
\end{eqnarray}

where $\rho_{{\scaleto{AB}{3.5pt}}}(\vec{r_{1}}, t\,;\vec{r_{2}}, t)$
is the joint density distributions in the two volumes. It provides the
joint probability of finding a mass element within each of the two
subvolumes. One can simplify <ref> to write the mutual
information between the two regions as,

\begin{eqnarray}
I_{{\scaleto{AB}{3.5pt}}}(t) = S_{{\scaleto{A}{3.5pt}}}(t) +
S_{{\scaleto{B}{3.5pt}}}(t) - S_{{\scaleto{AB}{3.5pt}}}(t)
\label{eq:miabs}
\end{eqnarray}

where the joint entropy $S_{{\scaleto{AB}{3.5pt}}}(t)$ can be
expressed as,

\begin{eqnarray}
\begin{split}
S_{{\scaleto{AB}{3.5pt}}}(t)&=-\int \int
\rho_{{\scaleto{AB}{3.5pt}}}(\vec{r_{1}}, t\,;\vec{r_{2}}, t)\,\log
\rho_{{\scaleto{AB}{3.5pt}}}(\vec{r_{1}}, t\,;\vec{r_{2}},
&=-\int \int
\bar{\rho}_{{\scaleto{A}{3.5pt}}}(t)\,\bar{\rho}_{{\scaleto{B}{3.5pt}}}(t)\,[1+\xi(\vec{r}_{{\scaleto{AB}{3.5pt}}},t)]\,\log\Big[\bar{\rho}_{{\scaleto{A}{3.5pt}}}(t)\,\bar{\rho}_{{\scaleto{B}{3.5pt}}}(t)\,\left(1+\xi(\vec{r}_{{\scaleto{AB}{3.5pt}}},t)\right)\Big]\,dV_{{\scaleto{A}{3.5pt}}}\,dV_{{\scaleto{B}{3.5pt}}}
\end{split}
\label{eq:scab}
\end{eqnarray}

where $\xi(\vec{r}_{{\scaleto{AB}{3.5pt}}},t)$ is the two-point
correlation function and $\bar{\rho}_{{\scaleto{A}{3.5pt}}}(t)$,
$\bar{\rho}_{{\scaleto{B}{3.5pt}}}(t)$ are the mean density in the two
regions $A$ and $B$. The separation vector between the centre of the
two subvolumes is
$\vec{r}_{{\scaleto{AB}{3.5pt}}}=\vec{r_{2}}-\vec{r_{1}}$. The
two-point correlation function
$\xi(\vec{r}_{{\scaleto{AB}{3.5pt}}},t)$ provides the excess
probability of finding two mass element separated by
$\vec{r}_{{\scaleto{AB}{3.5pt}}}$. The mutual information
$I_{{\scaleto{AB}{3.5pt}}}(t)$ quantifies the reduction of uncertainty
in the mass distribution within one volume given that we have the
complete knowledge of the mass distribution in the other. In other
words, it quantifies the expected gain in information about the mass
distribution in one volume when the other volume is observed.

We assume that the matter distribution in the universe is homogeneous
and isotropic on a scale of $V^{\frac{1}{3}}$. This allows us to write
= \bar{\rho}(t)$, $S_{{\scaleto{A}{3.5pt}}}(t) \approx
S_{{\scaleto{B}{3.5pt}}}(t) = S(t)$,
and $ \frac{dS_{{\scaleto{A}{3.5pt}}}(t)}{dt} \approx
\frac{dS_{{\scaleto{B}{3.5pt}}}(t)}{dt} = \frac{dS(t)}{dt}$.

The mean density $\bar{\rho}(t)$ and the two-point correlation
function $\xi(r_{{\scaleto{AB}{3.5pt}}},t)$ would evolve differently
in different cosmological models. The mutual information between the
two regions and its time evolution would thus depend on the
cosmological model.

§ RESULTS

§.§ Mutual information in a matter dominated Universe

We would like to calculate the mutual information between the mass
distributions in the two regions $A$ and $B$ in a matter dominated
Universe ($\Omega_m=1$). The average density in a matter dominated
Universe is $\bar{\rho}(t)=\frac{1}{6\,\pi\,G\,t^2}$ and the growing
mode of density perturbations is $D(t) \propto t^{\frac{2}{3}}$. In
the linear regime, the time evolution of the two-point correlation
function can be described as,
$\xi(r_{{\scaleto{AB}{3.5pt}}},t)=D^2(t)\,\xi(r_{{\scaleto{AB}{3.5pt}}})$. Here
the proportionality constant is absorbed in

[25] shows that the configuration entropy rate in a matter
dominated Universe is always negative i.e. $ \frac{dS(t)}{dt}<0$. Let
us write $ \frac{dS(t)}{dt}=-f(t)$.

The time evolution of the mutual information between the regions $A$
and $B$ can be written from <ref> as,

\begin{eqnarray}
\frac{d I_{{\scaleto{AB}{3.5pt}}}(t)}{dt} = -2 \frac{dS(t)}{dt}+\frac{d}{dt}
\int \int
\Big[\bar{\rho}^2(t)\,\left(1+D^{2}(t)\,\xi(r_{{\scaleto{AB}{3.5pt}}})\right)\Big]
\,\log
\Big[\bar{\rho}^2(t)\,\left(1+D^{2}(t)\,\xi(r_{{\scaleto{AB}{3.5pt}}})\right)\Big]
\,dV_{{\scaleto{A}{3.5pt}}}\,dV_{{\scaleto{B}{3.5pt}}}
\label{eq:miedu}
\end{eqnarray}

Simplifying <ref>, we get,

\begin{eqnarray}
\frac{d I_{{\scaleto{AB}{3.5pt}}}(t)}{dt} =-2 f(t)+I_{1}+I_{2}+I_{3}+I_{4}+I_{5}
\label{eq:miedu1}
\end{eqnarray}

\begin{eqnarray}
I_{1} =\frac{2}{9\,\pi^{2}\,G^{2}}\,\log(6\pi G) \int \int
\left(t^{-5}+\frac{2}{3}\,t^{-\frac{11}{3}}
\,\xi(r_{{\scaleto{AB}{3.5pt}}})\right)\,dV_{{\scaleto{A}{3.5pt}}}\,dV_{{\scaleto{B}{3.5pt}}}
\label{eq:i1}
\end{eqnarray}

\begin{eqnarray}
I_{2} =-\frac{1}{9\pi^{2}G^{2}}\, \int \int t^{-5} \left(1-4\log
\label{eq:i2}
\end{eqnarray}

\begin{eqnarray}
I_{3} = -\frac{1}{9\pi^{2}G^{2}}\, \int \int
\xi(r_{{\scaleto{AB}{3.5pt}}})\,t^{-\frac{11}{3}}
\left(1-\frac{8}{3}\,\log
\label{eq:i3}
\end{eqnarray}

\begin{eqnarray}
I_{4} =-\frac{1}{9\pi^{2}G^{2}}\, \int \int
\left(t^{-5}+\frac{2}{3}\,t^{-\frac{11}{3}}\,\xi(r_{{\scaleto{AB}{3.5pt}}})\right)\,\log
\left(1+t^{\frac{4}{3}}\,\xi(r_{{\scaleto{AB}{3.5pt}}})\right)\,dV_{{\scaleto{A}{3.5pt}}}\,dV_{{\scaleto{B}{3.5pt}}}
\label{eq:i4}
\end{eqnarray}

\begin{eqnarray}
I_{5} = \frac{1}{27\pi^{2}G^{2}}\, \int \int
\left(t^{-4}+t^{-\frac{8}{3}}\,\xi(r_{{\scaleto{AB}{3.5pt}}})\right)\,\frac{t^{\frac{1}{3}}\,\xi(r_{{\scaleto{AB}{3.5pt}}})}{1+t^\frac{4}{3}\,\xi(r_{{\scaleto{AB}{3.5pt}}})}\,dV_{{\scaleto{A}{3.5pt}}}\,dV_{{\scaleto{B}{3.5pt}}}
\label{eq:i5}
\end{eqnarray}

The integrals in the expressions of $I_{1}, I_{2}, I_{3}, I_{4}$ and
$I_{5}$ can not be simplified further due to the lack of
symmetry. However, one can easily analyze the time dependence of these
expressions. Observations show that the galaxy two-point correlation
function has a nearly universal dependence on pair separation $r$ as
$\xi(r) \sim r^{-1.8}$. The terms involving
$\xi(r_{{\scaleto{AB}{3.5pt}}})$ in these expressions will have a
smaller magnitude. The <ref>, <ref>,
<ref>, <ref> and <ref> have strong time
dependence. $I_{2}$ and $I_{3}$ become positive for larger values of
time. Only $I_{4}$ remains negative at all time. The sum
$I=(I_{1}+I_{2}+I_{3}+I_{4}+I_{5})$ is positive but decays towards
zero with increasing time. The term $f(t)$ in <ref> has
a much weaker time dependence as compared to $I$
[25]. This leads to $\frac{d
I_{{\scaleto{AB}{3.5pt}}}(t)}{dt}<0$ which implies that the mutual
information between two independent regions decreases with time in a
matter dominated Universe.

It may be noted that <ref> depends on the size of the
regions $A$ and $B$ and the separation between them.  The integrals in
<ref> would be carried out over different volumes when
there is a change in the size of the two regions. The separations
between the different pairs of subvolumes would change with the
distance between the two regions. The integrals in <ref>
will have different values since the two-point correlation function
changes with the separation.  If the two regions are separated by a
very large distance compared to the dimensions of the two regions then
the integrals in <ref> would loss their physical

§.§ Mutual information in a $\Lambda$-dominated Universe

Now we would like to calculate $\frac{d
I_{{\scaleto{AB}{3.5pt}}}(t)}{dt}$ in a $\Lambda$-dominated
Universe. We have $\bar{\rho}(t)=\bar{\rho}=$constant and the growing
mode of density perturbations is $D(t)=k=$ constant in
$\Omega_{\Lambda}=1$ Universe. The time evolution of the mutual
information in such a universe can be expressed as,

\begin{eqnarray}
\frac{d I_{{\scaleto{AB}{3.5pt}}}(t)}{dt} = -2
\frac{dS(t)}{dt}+\frac{d}{dt} \int \int
\Big[\bar{\rho}^2\,\left(1+k^{2}\,\xi(r_{{\scaleto{AB}{3.5pt}}})\right)\Big]
\,\log
\Big[\bar{\rho}^2\,\left(1+k^{2}\,\xi(r_{{\scaleto{AB}{3.5pt}}})\right)\Big]
\,dV_{{\scaleto{A}{3.5pt}}}\,dV_{{\scaleto{B}{3.5pt}}}
\label{eq:milam}
\end{eqnarray}

[25] show that $\frac{dS}{dt}=0$ in a $\Lambda$-dominated
Universe. So we have $ \frac{d I_{{\scaleto{AB}{3.5pt}}}(t)}{dt}=0$.
Clearly, $I_{{\scaleto{AB}{3.5pt}}}=$ constant in a
$\Lambda$-dominated Universe. There would be a constant mutual
information between the regions $A$ and $B$ at all times in such a

§.§ Mutual information in the $\Lambda$CDM model, the dynamical dark energy models and the holographic dark energy models

It is clear that the joint entropy between the two regions $A$ and $B$
plays a negligible role in the time evolution of the mutual
information. The configuration entropy rate $\frac{dS}{dt}$ determines
the time evolution of the mutual information between the two
regions. The configuration entropy rates have been calculated for the
$\Lambda$CDM model, different dynamical dark energy models and the
holographic dark energy models in the literature [48, 49, 50, 51]. The configuration entropy rate
decreases to reach a minimum and then increases with time in all these
models. However, the location and the amplitude of the minimum depend
on the models. The location of the minimum precisely indicates the
epoch of dark energy domination predicted by the relevant model. The
average density will fall faster in such models as compared to the
matter dominated Universe. So the joint entropy term would contribute
negligibly to the evolution of the mutual information in all these
models. The time evolution of the mutual information in a given model
will be thus entirely determined by the behaviour of the configuration
entropy rate in that model. The mutual information between $A$ and $B$
in all these models would initially decrease with time and eventually
hit a minimum. The mutual information would increase after this
minimum once the dark energy starts to dominate the dynamics of the

§ CONCLUSIONS

We analyze the time evolution of the mutual information between
disjoint regions of the universe in different cosmological models. The
mutual information here quantifies the reduction of uncertainty in the
knowledge of the matter distribution in one region given that we have
complete knowledge of it in the other region. In another words, the
mutual information provides some knowledge about the matter
distribution in one region provided we have the complete knowledge of
the matter distribution in the another region. A zero mutual
information indicates that the mass distribution in the two regions
are statistically independent. We do not separately calculate the
mutual information between the two disjoint regions but obtain an
expression for its time evolution from the definition. It may be noted
that the mutual information is positive or zero by definition. It can
not assume negative values. However, the rate of change of the mutual
information can be negative, positive or zero.

The time evolution of the mutual information between disjoint regions
of the universe is primarily determined by the dynamics of the
expansion and the growth rate of the density perturbations. We find
that the mutual information decreases continuously in a purely matter
dominated Universe whereas it stays constant in a purely
$\Lambda$-dominated Universe. So the disjoint regions become
statistically independent in a matter-dominated Universe whereas they
remain entangled forever in a $\Lambda$-dominated Universe. The mutual
information decreases to reach a minimum and then increases with time
in the $\Lambda$CDM model, dynamical dark energy models and the
holographic dark energy models.  Clearly, the time evolution of the
mutual information is governed by the changes in the configuration
entropy of the matter distribution in the universe. The change in the
joint entropy between the mass distributions in the two regions does
not contribute significantly to this evolution.

The two regions $A$ and $B$ are causally connected. However, they may
be separated by large distance. In reality, one can not measure the
mass distributions in the two volumes simultaneously. One can infer
some information about the mass distribution in one volume while
observing the other. These information corresponds to the same cosmic
time. It is worthwhile to mention here that a mutual information
between the two causally connected regions may also introduce a
non-negligible dynamical entanglement [43]. The effect of
such dynamical entanglement is not considered in the present work.
Further, one can also consider the contributions from the higher order
correlations. A non-zero three-point correlation function would modify
the joint probabilities in <ref>. We plan to address these
issues in future works. It would be also interesting to study the
evolution of mutual information in the inhomogeneous cosmological
models. The presence of long range correlations in the mass
distributions can significantly modify the mutual information and its

The continuous dissipation of the configuration entropy during the
matter dominated era demands enormous entropy production that can
counterbalance this loss and maximize the entropy [25].
The accelerated expansion of the universe provides an avenue for
maximum entropy production in accordance with the second law of
thermodynamics. It is interesting to note that the evolution of the
mutual information is strongly sensitive to the cosmological constant
or dark energy. This implies that the mutual information may have
deeper connections to the dark energy and the accelerated expansion of
the universe.

§ ACKNOWLEDGEMENTS

The author thanks the anonymous reviewers for some valuable comments and suggestions on the dtaft. This research was funded by the SERB, DST, Government of India through the project CRG/2019/001110. The author would also like to acknowledge IUCAA, Pune for providing support through associateship programme.

[1] Martinez, V. J.; Coles, P. Correlations and Scaling in the QDOT Redshift Survey. Astrophys. J. 1994, 437, 550-555

[2] Borgani, S. Scaling in the Universe. Phys. Rep. 1995, 251, 1-152

[3]  Guzzo, L. Is the universe homogeneous? (On large scales) New Astron. 1997, 2, 517-532

[4] Cappi, A.; Benoist, C.; da Costa, L. N.; Maurogordato, S. Is the Universe a fractal? Results from the Southern Sky Redshift Survey 2. Astron. Astrophys. 1998, 335, 779-788

[5] Bharadwaj, S.; Gupta, A. K.; Seshadri, T. R. Nature of clustering in the Las Campanas redshift survey. Astron. Astrophys. 1999, 351, 405-412

[6] Pan, J.; Coles, P.  Large-scale cosmic homogeneity from a multifractal analysis of the PSCz catalogue. Mon. Not. R. Astron. Soc. 2000, 318, L51-L54

[7] Yadav, J.; Bharadwaj, S.; Pandey, B.; Seshadri, T. R. Testing homogeneity on large scales in the Sloan Digital Sky Survey Data Release One. Mon. Not. R. Astron. Soc. 2005, 364, 601-606

[8]  Hogg, D. W.; Eisenstein, D. J.; Blanton, M. R.; Bahcall, N. A.; Brinkmann, J.; Gunn, J. E.; Schneider, D. P. Cosmic Homogeneity Demonstrated with Luminous Red Galaxies. Astrophys. J. 2005, 624, 54-58

[9] Sarkar, P.; Yadav, J.; Pandey, B.; Bharadwaj, S. The scale of homogeneity of the galaxy distribution in SDSS DR6. Mon. Not. R. Astron. Soc. 2009, 399, L128-L131

[10] Scrimgeour, M. I.; Davis, T.; Blake, C. et al. The WiggleZ Dark Energy Survey: the transition to large-scale cosmic homogeneity. Mon. Not. R. Astron. Soc. 2012, 425, 116-134

[11] Nadathur, S. Seeing patterns in noise: gigaparsec-scale ‘structures’ that do not violate homogeneity. Mon. Not. R. Astron. Soc. 2013, 434, 398-406

[12] Pandey, B.; Sarkar, S. Testing homogeneity in the Sloan Digital Sky Survey Data Release Twelve with Shannon entropy. Mon. Not. R. Astron. Soc. 2015, 454, 2647-2656

[13] Pandey, B.; Sarkar, S. Probing large scale homogeneity and periodicity in the LRG distribution using Shannon entropy. Mon. Not. R. Astron. Soc. 2016, 460, 1519-1528)

[14] Riess A. G.; Filippenko A. V.; Challis P.; Clocchiatti A.; Diercks A.; Garnavich P. M.; Gilliland R. L., et al. Observational Evidence from Supernovae for an Accelerating Universe and a Cosmological Constant. Astron. J. 1998, 116, 1009-1038

[15] Perlmutter S.; Turner M. S.; White M. Constraining Dark Energy with Type Ia Supernovae and Large-Scale Structure. Phys. Rev. Lett. 1999 83, 670-673

[16] Buchert, T.; Ehlers, J. Averaging inhomogeneous Newtonian cosmologies. Astron. Astrophys. 1997, 320, 1-7

[17] Schwarz, D. J. Accelerated expansion without dark energy. 2002 arXiv:astro-ph/0209584

[18]  Kolb, E. W.; Matarrese, S.; Riotto, A. On cosmic acceleration without dark energy. New J. Phys. 2006, 8, 322-346

[19] Buchert, T.  Dark Energy from structure: a status report. Gen. Relativ. Gravit. 2008, 40, 467-527

[20]  Ellis, G. F. R. Inhomogeneity effects in cosmology. Class. Quantum Gravity 2011, 28, 164001-164030

[21] Pavon, D. Why the expansion of the Universe is accelerated? Int. J. Geom. Methods Mod. Phys. 2014 11, 1460007

[22] Radicella N.; Pavón D. A thermodynamic motivation for dark energy. Gen. Relativ. Gravit. 2012, 44, 685-702

[23] Mimoso J. P.; Pavón D. Entropy evolution of universes with initial and final de Sitter eras. Phys. Rev. D 2013, 87, 047302

[24] Ferreira P. C.; Pavón D. Thermodynamics of nonsingular bouncing universes. Eur. Phys. J. C 2016, 76, 37-41

[25] Pandey B. Does information entropy play a role in the expansion and acceleration of the Universe? Mon. Not. R. Astron. Soc. 2017, 471, L77-L81

[26] Pandey B. Configuration entropy of the Cosmic Web: Can voids mimic the dark energy? Mon. Not. R. Astron. Soc. 2019, 485, L73-L77

[27] Martyushev, L.; Seleznev, V.D. Maximum entropy production principle in physics, chemistry and biology. Phys. Rep. 2006, 426, 1-45

[28] Martyushev, L. Entropy and Entropy Production: Old Misconceptions and New Breakthroughs. Entropy 2013, 15, 1152-1170

[29] Gott, J. R., III; Jurić, M.; Schlegel, D., et al. A Map of the Universe. Astrophys. J. 2005, 624, 463-484

[30] Clowes R. G.; Harris K. A.; Raghunathan S.; Campusano L. E.; Söchting I. K.; Graham M. J. A structure in the early Universe at $z \sim 1.3$ that exceeds the homogeneity scale of the R-W concordance cosmology. Mon. Not. R. Astron. Soc. 2013, 429, 2910-2916

[31] Lopez A. M.; Clowes R. G.; Williger G. M. A Giant Arc on the Sky. Mon. Not. R. Astron. Soc. 2022, 516, 1557-1572

[32] Friday T.; Clowes R. G.; Williger G. M. Correlated orientations of the axes of large quasar groups on Gpc scales. Mon. Not. R. Astron. Soc. 2022, 511, 4159-4178

[33] Keenan R. C.; Barger A. J.; Cowie L. L. Evidence for a  300 Megaparsec Scale Under-density in the Local Galaxy Distribution. Astrophys. J. 2013, 775, 62

[34] Colin J.; Mohayaee R.; Rameez M.; Sarkar S. Evidence for anisotropy of cosmic acceleration. Astron. Astrophys. 2019, 631, L13

[35] Secrest N. J.; von Hausegger S.; Rameez M.; Mohayaee R.; Sarkar S.; Colin J. A Test of the Cosmological Principle with Quasars. Astrophys. J. Lett. 2021, 908, L51

[36] Sheth R. K.; Diaferio A. How unusual are the Shapley supercluster and the Sloan Great Wall? Mon. Not. R. Astron. Soc. 2011, 417, 2938-2949

[37] Park C.; Choi Y.-Y.; Kim J.; Gott J. R.; Kim S. S.; Kim K.-S. The Challenge of the Largest Structures in the Universe to Cosmology. Astrophys. J. Lett. 2012, 759, L7-L13

[39] Peebles, P. J. E. The Large-Scale Structure of the Universe. Princeton University Press, 1980, Princeton, New Jersey

[40] Mecke K. R.; Buchert T.; Wagner H. Robust morphological measures for large-scale structure in the Universe. Astron. Astrophys. 1994, 288, 697-704

[41] Hosoya A.; Buchert T.; Morita M. Information Entropy in Cosmology. Phys. Rev. Lett. 2004, 92, 141302

[42] Akerblom N.; Cornelissen G. Relative entropy as a measure of inhomogeneity in general relativity. J. Math. Phys. 2012, 53, 012502

[43] Wiegand A.; Buchert T. Multiscale cosmology and structure-emerging dark energy: A plausibility analysis. Phys. Rev. D 2010, 82, 023523

[44] Czinner V. G.; Mena F. C. Relative information entropy in cosmology: The problem of information entanglement. Phys. Lett. B 2016, 758, 9-13

[45] Espinosa-Portalés L.; García-Bellido J. Long-range enhanced mutual information from inflation. Phys. Rev. D 2021, 103, 043537

[46] Shiba N. Direct Calculation of Mutual Information of Distant Regions. J. High Energy Phys. 2020, 2020, 182

[47] Cardy J. Some Results on Mutual Information of Disjoint Regions in Higher Dimensions. J. Phys. A 2013, 46, 285402

[48] Das B.; Pandey B. Configuration entropy in the $\Lambda$CDM and the dynamical dark energy models: Can we distinguish one from the other?  Mon. Not. R. Astron. Soc. 2019, 482, 3219-3226

[49] Pandey B.; Das B. A new method to probe the mass density and the cosmological constant using configuration entropy. Mon. Not. R. Astron. Soc. Lett. 2019, 485, L43-L47

[50] Bhattacharjee S. Growth Rate and Configurational Entropy in Tsallis Holographic Dark Energy. Eur. Phys. J. C 2021 81, 217

[51] Das B.; Pandey B. A study of holographic dark energy models using configuration entropy. Res. Astron. Astrophys. 2023, 23, 065003


Guillaume Derval, Frédéric Docquier and Pierre Schaus

11institutetext: ICTEAM, UCLouvain, Louvain-la-Neuve, Belgium 22institutetext:
IRES, UCLouvain, Louvain-la-Neuve, Belgium
22email<EMAIL_ADDRESS>
# An aggregate learning approach for interpretable semi-supervised population
prediction and disaggregation using ancillary data

Guillaume Derval ✉ 11 0000-0002-6700-3519 Frédéric Docquier 22 Pierre Schaus
11 0000-0002-3153-8941

###### Abstract

Census data provide detailed information about population characteristics at a
coarse resolution. Nevertheless, fine-grained, high-resolution mappings of
population counts are increasingly needed to characterize population dynamics
and to assess the consequences of climate shocks, natural disasters,
investments in infrastructure, development policies, etc. Dissagregating these
census is a complex machine learning, and multiple solutions have been
proposed in past research. We propose in this paper to view the problem in the
context of the aggregate learning paradigm, where the output value for all
training points is not known, but where it is only known for aggregates of the
points (i.e. in this context, for regions of pixels where a census is
available). We demonstrate with a very simple and interpretable model that
this method is on par, and even outperforms on some metrics, the state-of-the-
art, despite its simplicity.

###### Keywords:

Disaggregation Aggregate Learning GIS

##  1 Introduction

Most countries periodically organize rounds of censuses of their population at
a granularity that differs from country to country. The level of
disaggregation is often governed by the administrative division of the
country. Although census data are usually considered as accurate in terms of
population counts and characteristics, the spatial granularity, that is
sometimes in the order of hundreds of square kilometers, is too coarse for
evaluating local policy reforms or for making informed decisions about health
and well-being of people, economic and environmental interventions, security,
etc. (see [23]). For example, fine-grained, high-resolution mappings of the
distribution of the population are required to assess the number of people
living at or near sea level, near hospitals, in the vicinity or airports and
highways, in conflict areas, etc. They are also needed to understand how
population movements react to various types of shocks such as natural
disasters, conflicts, plant creation and closures, etc.

Multiple methods can be used to produce gridded data sets (also called
rasters), with pixels of a relatively small scale compared to the
administrative units of the countries. Notably, the Gridded Population of the
World (GPW) [3] provides a gridded dataset of the whole world, with a
resolution of 30 arcseconds (i.e. pixels of approximately 1 km2 at the
equator). The method used in GPW assumes that population is uniformly
distributed within each administrative unit. With minor adjustments for
boundaries and water bodies, the density of the population in a pixel is
identical to the density of the population of the administrative unit in the
underlying census.

More advanced and successful models rely on ancillary and remotely sensed
data, and are trained using machine learning techniques (see [6, 20, 18]).
These data can include information sensed by satellite (nighttime light
intensities, temperature, etc.) or data provided by NGOs and governments (on
health facilities, road infrastructure, etc.). These models are used to
disaggregate coarse census data into a grid of small pixels, using ancillary
data to predict the population distribution within the census unit. The level
of resolution is either imposed by the method or decided by the user. As an
example, such models predict that people concentrate near the road
infrastructure rather than in cultivation fields; or that the density of the
population in highly elevated regions is smaller than near the sea level.

These models show various degrees of usability and accuracy. Among the methods
cited above, RF [20] (for Random Forest, the family of algorithms on which the
method is based) gives the best results using remotely sensed and ancillary
data. However, the results of RF are seemingly artificial. The right panel of
Fig 1 shows a fine-grained population mapping of Cambodia generated with the
RF model from [20]. Geometrical shapes clearly emerge from the population
disaggregation. On the contrary, as depicted on the left panel of Fig 1, our
PCD (for Pure Census Disagregation) method generates maps with a natural
distribution aspect.

Figure 1: Cambodian population maps obtained by PCD-LinExp-PPE (a
specialization of PCD using the PCD-LinExp model type, introduced in a later
section) and RF. The bottom maps are a zoom of a specific, moderately
populated region of Cambodia. The red box highlights a region where RF
produces seemingly artificial results: it creates a circle around a (non-
displayed) hospital, and saturates near the road network.

All methods in the literature are converted into standard supervised
regression learning tasks. Supervised regression learning aims to predict an
output value associated with a particular input vector. In its standard form,
the training set contains an individual output value for each input vector.
Unfortunately, the disaggregation problem does not directly fit into this
supervised regression learning framework since the prediction function is not
directly available for input vectors. In a disaggregation problem, the input
consists of a partition of the training set (the pixels of each unit) and for
each partition, the sum of the output values is constrained by the census
count. This framework is exactly the one introduced as the aggregated output
learning problem by [16]. Our PCD method conceives the formulation of the
disaggregation problem as an aggregated output learning problem. PCD is able
to train the model based on a much larger training set composed of pixels.
This approach is parameterized by the error function that the user seeks to
minimize.

PCD can be used with a large variety of model types. In this paper, we
consider a specific model, arguably simplistic, that is called PCD-LinExp. The
PCD-LinExp model can be written as:

$f_{\bm{\theta}=(\bm{a},b,c)}(\bm{D}_{p_{x},p_{y}})=\max(0,\exp(\bm{a}^{T}\bm{D}_{p_{x},p_{y}}+b)+c),$

where $f$ is the predicted population count, the vector $\bm{a}$ and scalars
$b$ and $c$ are parameters to estimate, and $\bm{D}_{p_{x},p_{y}}$ is a set of
sensed/ancillary covariates of a given pixel located at $p_{x},p_{y}$.

PCD-LinExp is arguably simplistic. This is by design, for two reasons. First,
it is interpretable: each available remotely sensed covariate is associated to
a single parameter in the vector $\bm{a}$. Second, the aim of this paper is to
show that the aggregate learning approach is particularly suitable and gives
good results on this particular problem. Future research will focus on more
complex models. Despite its simplicity, we actually show in the result section
that, due to its comparatively small unadjusted error compared to other
methods, PCD-LinExp can be used to predict population counts when census data
are not available, something that is not possible with existing methods.

As a case study, we experiment it on Cambodia using various sets of remotely
sensed/ancillary data. While this is not the main focus of this paper, we
shortly interpret the result of our method on Cambodia. Our main result, the
disaggregated map of Cambodia, is depicted on the left panel of Fig 1. This
paper also discusses methodological issues raised by existing approaches. In
particular, we demonstrate that the previously used error metrics are biased
when available census data involves administrative units with highly
heterogeneous surfaces and population densities. We propose alternative
metrics that better reflect the accuracy properties that should be fulfilled
by a sound disaggregation approach. We then present the results for Cambodia
and compare methods using various error metrics, providing statistical
evidence that PCD-LinExp generates the most accurate results.

###  1.1 Notations and definitions

A country is divided into administrative units (abbreviated as AUs), that are
part of various administrative levels and form a hierarchy. For example, in
Cambodia, the AUs at the finest level of the AU hierarchy are called communes,
which are grouped into districts. In most countries, the census is conducted
at the finest possible level of the AU hierarchy. For simplicity, throughout
this paper, we refer to an AU at the finest administrative level as units, and
use the term superunit for an AU at the coarser level.

######  Definition 1 (Units, Superunits)

Let $\mathcal{U}$ be the set of units and $\mathcal{S}$ the set of superunits.
Let us denote by $\operatorname{units}(s)=\\{u\mid u\in s\\}$ the set of units
in a superunit $s$. The set of pixels in a unit is
$\operatorname{pixels}(u)=\\{p\mid p\in u\\}$ and by extension the set of
pixels in a superunit is
$\operatorname{pixels}(s)=\bigcup_{u\in\operatorname{units}(s)}\operatorname{pixels}(u)$.
A pixel $p$ is a triple $(p_{x},p_{y},p_{w})$, where $(p_{x},p_{y})$ is the
position of the pixel and $p_{w}$ the weight of the pixel in the unit, which
is the area of the pixel covered by the unit;
$\operatorname{surface}(u)=\sum_{p\in\operatorname{pixels}(u)}p_{w}$ is the
total surface of the unit u; $\operatorname{surface}(s)$ for superunits is
defined similarly. We further define $\operatorname{census}(u)$ as the
population count of unit $u$ and
$\operatorname{census}(s)=\sum_{u\in\operatorname{units}(s)}\operatorname{census}(u)$
as the population count of superunit $s$.

An example illustrating the definitions is provided in Fig 2.

Figure 2: Notation example.

###  1.2 Existing methods for disaggregation

None of the existing approaches use the aggregate learning framework [16].
They use standard supervised learning approaches where the training set
contains an individual output value for each input vector. The input data
available for the disaggregation problem does not directly fit into this
framework since the prediction value for each pixel is not available; an
aggregate form is available, however (the census count for a unit). Existing
disaggregation methods first transform the data to make them fit into the
standard regression problem. Existing approaches are differentiated by their
transformation method, by the error metric used to evaluate the quality of the
learned function, and by the family of machine learning algorithm used to
train the model.

##### Areal weighting

[8, 10] is a straightforward method to redistribute the census data to a grid
of pixels. The population predicted for each pixel is the sum of the
contribution of each unit that contains the pixel. The contribution of a unit
to a pixel $p$ is defined as
$\operatorname{prediction}_{u}(p):=\frac{p_{w}}{\sum_{p^{\prime}\in\operatorname{pixels}(u)}p^{\prime}_{w}}\cdot\operatorname{census}(u).$
The Gridded Population of the World (GPW) dataset [3] provides a raster (grid)
of the population density/count at a fine resolution for the whole world.
After census data is gathered, it is converted into a grid using an additional
dataset on water bodies only. It is then a direct application of areal
weighting, where each pixel weight $p_{w}$ is adjusted for the ratio of water
bodies inside the pixel.

##### Dasymetric modeling

is an extension of areal weighting using ancillary data [7, 24, 14]. A weight
is inferred from this ancillary data, which redistribute the population
according to these weights. Land use/cover are often used as ancillary data
[5, 15, 21]. Methods used to select the weights vary from purely heuristic
choice of weights [5] to use of Expectation-Maximization or Logit regression
[9], and other linear regression models [2]. The methods below can all be seen
as members of the Dasymetric modeling family, with the specificity that they
all use more complex machine learning algorithm to predicts the weights
accorded to each ancillary data type. They are also (mostly) agnostic on the
type of the data and can handle multiple types at the same time, in a generic
way, which is not the case for the methods presented above.

##### Random Forest (RF) method: pixel aggregation

The methodology followed by [20] solves a regression problem in order to
predict, for any pixel on a grid, the logarithm of the population count in
this pixel. It is based on nonparametric machine learning techniques (RF) and
uses observable covariates such as presence of hospitals, road networks,
nightlight intensities, etc. Each covariate is mapped/projected on a grid of
pixels (also denoted raster) spreading over the selected country. The
resolution of this grid is decided by the user. In [20], it is set at
$8.33e^{-4}$ degrees, which approximately corresponds to 100 m sided square
pixels at the equator. A key contribution of [20] is the methodology to
extract a training set used by the RF [1] machine learning algorithm.

One data point is produced for each unit in the original training set. Each
data point is a vector, containing a value (binary or floating point) for each
covariate. These values are resulting from the pixel aggregation step, which
varies with the type of covariate: for continuous covariates, the result of
the aggregation of pixels is the mean of the covariate values of these pixels.
For binary covariates, the mode (class majority) of the pixel’s covariate
values is taken as a summary for the unit. We thus obtain $|U|$ data points
(remember that $|U|<<$ the number of pixels), having each a fixed size,
suitable for standard supervised machine learning techniques, such as random
forests, neural networks, or SVMs. The training set is then used to learn a RF
regression model.

Let us denote the training data as $\bm{D}$, a tensor in three dimensions,
such that $\bm{D}_{x,y,f}$ gives the value of the covariate $f$ at position
$x,y$. $\bm{D}_{x,y}$ gives the vector of all the covariate values at position
$x,y$. The regression $RF(\bm{D}_{x,y})$ is thus able to predict a population
density for a pixel with covariates values $\bm{D}_{x,y}$.

##### Corrected output of linear regression model (dasymetric redistribution)

The method used by [20] does not directly consider the output of
$RF(\bm{D}_{x,y})$ as the predictor for the population density. Instead, it
considers the output of the $RF$ for each pixel as a weighting layer used for
subsequent dasymetric redistribution (also called pycnophylactic
redistribution [22]).

For each unit $u$, a correcting factor $w(u)$ is computed as

$w(u)=\frac{\operatorname{census}(u)}{\sum_{p\in\operatorname{pixels}(u)}p_{w}\cdot
RF(\bm{D}_{p_{x},p_{y}})}.$

The predicted density estimator for the pixel $p$ lying in unit $u$ is then

$\operatorname{prediction}(p)=RF(\bm{D}_{p_{x},p_{y}})\cdot w(u).$

During the training of the model, this would, of course, lead to an error of
zero. Therefore, the dasymetric redistribution is only applied on superunits.
This part of the method obviously limits the usage of the method to situations
where the census data is available, at least at a coarse resolution. We call
the error _adjusted_ if it is dasymetrically redistributed, and _unadjusted_
otherwise.

##### Learning from heuristic, using satellite images and deep neural networks

[18] relax the problem in a different way; they use an external census grid,
produced in a heuristic way [4], as ground truth. More precisely, they train
their model based on a grid made from the 2000 US census, and validate it by
measuring the prediction accuracy on the 2010 US census. This “relaxation” of
the initial problem allows using standard supervised machine learning
techniques. As these census grids are extrapolated from the official census,
they suffer from higher errors, making the dataset noisier. However, compared
to the Pixel Aggregation relaxation, it creates a much larger dataset, with
one data point per pixel rather than one point per unit. In their paper,
Robinson et al. do not use ancillary data. They use satellite images in the
visible and invisible spectrum. They use a Convolutionnal Neural Network model
(abbreviated CNN, see [13]), which is a modified VGG-A neural net [19]. The
neural net is given a window of pixels of size $74\times 74$ (each pixel being
approximately $15m\times 15m$), and predicts a (single) class for this set of
pixels. It is a classification task, each class representing a given
population range, obtained after binning the US censuses. The method uses a
different type of data, among which some are not available in most countries
(census grids). Hence, the interest in comparing this method with others is
limited, and is thus not included here. Another method, LL [6] use the same
preprocessing method (using a nearest neighbor method to produce the heuristic
grid). LL is a prediction method that converts all covariates into a grid with
pixels of size 250m$\times$250m (approx. 0.0625 km2), and then assigns a
ground truth for the population count of each pixel. The latter is proxied
with the mean population density of the unit they belong to. As for the
previous method, we do not compare our method to LL. The data types are
different (they only rely on satellite images) and the results are too coarse
(both RF and our method uses pixels that are 6.25 times smaller in surface).
In the original paper [6], the LL method produces good results but only after
applying the dasymetric redistribution (as with RF).

##  2 Data preparation

In order to compare our PCD-LinExp method with existing ones, various data
sets have been collected for Cambodia. These include both data sensed by
satellites , and data provided by the Cambodian administration and by NGOs .
The description of the dataset, along with the preprocessing done, is
available in the additional resources of this paper. Both the data and the
preprocessing are similar to the ones in the article describing the RF method
[20].

In our experiments, we select a resolution of $8.33e^{-4}$ degrees, which
approximately corresponds to squares of 100m sides at the equator. We then
obtain a stack of rasters, which can be represented as a tensor in three
dimensions. We call this tensor $\bm{D}$, with $\bm{D}_{x,y,f}$ being the
pixel of the feature $f$ at position $x,y$. In the case of Cambodia, the
dimension of this tensor is $6617\times 5468\times 42$. We use $\bm{D}_{x,y}$
as notation for the vector of all the covariates at position $x,y$. One
difference between our method to prepare data and the one in [20] is that our
data are standardized. In order to reduce numerical errors and make covariates
equally important at the beginning of the computation, each feature is
centered (resulting covariates have their mean equal to zero) and transformed
such that they have a unit variance.

##  3 Evaluation

###  3.1 Error metrics

Assuming the finest ground truth we have is the administrative unit of the
census, all errors metrics are in the form
$e(\operatorname{census},\operatorname{prediction},\operatorname{surface},\mathcal{U})$,
where $\operatorname{census}$, $\operatorname{prediction}$ and
$\operatorname{surface}$ are functions $\mathcal{U}\rightarrow\mathbb{R}^{+}$
that give the effective and predicted population of a unit and its surface,
respectively.

Existing methods[20, 6, 18] mostly use the RMSE criterion and its variants to
estimate the prediction error. In our context, the RMSE is defined as:

$e_{RMSE}=\sqrt{\frac{\sum_{u\in\mathcal{U}}\left(\operatorname{prediction}(u)-\operatorname{census}(u)\right)^{2}}{|\mathcal{U}|}}.$

RMSE is biased towards highly populated, low-surface administrative units. An
analysis of the data on Cambodia shows that there are smaller-surface units
than high-surface ones, and that the small ones tend to be more populated. In
most countries, the division of the country in administrative units is such
that the number of inhabitants in each unit is not too large.

Since there are more small regions than large ones, methods that better
predict these regions have better scores. Said differently, $e_{RMSE}$
overemphasizes errors on small, highly populated regions. We are aware that
this is not always a desirable property. Users might be more interested in
accurately disaggregating high-surface administrative units than smaller units
with already more precise population count.

Fig 3 provides an example of two sets of units, with the same population, but
with different predictions. Both of them have $e_{RMSE}=700$. As can be seen,
the example on the left panel makes very important errors on units U3 (error =
496) and U2 (272), where actual population counts are relatively small.
However, these errors do not contribute a lot to the RMSE and are instead
absorbed by the error of unit U1 (1000), which is relatively small compared to
its population, but large in absolute value. On the contrary, the example on
the right panel of Fig 3 gives a more balanced error distribution. In absolute
value the error is greater in U1. However, compared to the actual population
counts, they are smaller in each administrative unit.

Figure 3: Example of administrative units with equal RMSE. In this diagram,
for simplicity, pop and pred are respectively abbreviations for
$\operatorname{census}$ and $\operatorname{prediction}$.

We postulate that a sound error metric for the disaggregation problem should
be expressed in per-pixel units, therefore being independent of the surface of
each unit and be based on relative population deviations rather than on
absolute population counts.

We introduce the Per-Pixel-Error (PPE):

$e_{PPE}=\frac{1}{\sum_{u\in\mathcal{U}}\operatorname{surface}(u)}\sum_{u\in\mathcal{U}}\operatorname{surface}(u)\cdot\frac{|\operatorname{prediction}(u)-\operatorname{census}(u)|}{\operatorname{census}(u)},$

which is a surface-weighted average of relative deviations from the effective
population (i.e. a Weighted Mean Absolute Percentage error, with surfaces as
weights). PPE is thus independent from the size of the unit. Moreover, the use
of relative deviations allows to effectively compare less-populated units to
more-populated ones. In the examples from Fig 3, considering that all units
have the same surface, the left example has a PPE of $86.5\%$ while the right
one has a PPE of $10\%$. Remember that they have exactly the same RMSE (700).

A squared variant of PPE, the Per-Pixel Squared Error (PPSE), is defined as:

$e_{PPSE}=\frac{1}{\sum_{u\in\mathcal{U}}\operatorname{surface}(u)}\sum_{u\in\mathcal{U}}\operatorname{surface}(u)\cdot\frac{(\operatorname{prediction}(u)-\operatorname{census}(u))^{2}}{\operatorname{census}(u)},$

(which is a Weighted Root Mean Squared Error metric weighed by surface) which
overweights greater error ratios.

In the quantitative analysis below, we conduct our experiments using PPE and
PPSE as error metrics. We also use RMSE so as to compare our results with
previous works. MAE was used in previous works, but it suffers from the same
shortcomings as RMSE and does not provide additional relevant information, and
as such is omitted in our experiments.

###  3.2 Comment on dasymetric redistribution

As discussed previously, dasymetric redistribution is a commonly used
processing method to rescale (ex-post) the population in a given surface using
a weight map. This method is used by RF [20]. The main argument justifying the
use of dasymetric redistribution is the full exploitation of the available
information included in the population census. If a disaggregation method
predicts that there are 3,500 people in a given AU, and that the census
reports an actual population of 4,500, dasymetric redistribution consists of
multiplying all pixels within the AU by $\frac{4500}{3500}$, thus reducing the
visible error.

While the method seems intuitive and duly justified, it actually increases the
prediction error in some cases. This is not only due to the fact that the
population census is, like all types of databases, subject to errors and
various noises. It is also due to the fact that large local errors in a given
region are transferred to the whole region.

As an example, let us consider the case where we have two units composed each
of one pixel. These units are grouped into a single superunit. Table 1 shows
(on the left) the actual and predicted population counts of these units, along
with the adjusted prediction after dasymetric redistribution on the superunit.
The total population of the superunit is $200$, but the predicted population
is $160$. The correcting factor to be applied is $\frac{200}{160}=1.25$. The
right part of Table 1 shows the error metrics applied before and after the
dasymetric redistribution. In our example, both metrics increase after the
redistribution. The reason is the error made originally on unit B. The low
prediction forces the redistribution to increase the global count,
dramatically increasing the error on unit A.

Table 1: Example for dasymetric redistribution

Unit | Real pop. | Pred. pop. | Adj. pop.
---|---|---|---
A | 100 | 120 | 150
B | 100 | 40 | 50

Metric | Unadjusted | Adjusted
---|---|---
RMSE | $\sqrt{2000}$ | $\sqrt{2500}$
PPE | $40\%$ | $50\%$

Another drawback of dasymetric redistribution is that it breaks the dichotomy
between inputs and outputs in the validation stage. Indeed, the dasymetric
redistribution is applied to superunits (it would lead to an error of zero if
applied to units). This implies that part of the prediction truth values are
used to correct the output, before prediction errors are computed. We leave as
an open question the relevance of using dasymetric redistribution. In
practice, we notice that it sometimes increases the error in our experiments.
Dasymetric redistribution also tends to make comparisons between methods
difficult, as it squeezes the original errors. To allow comparison with
previous works, we produce both the unadjusted and adjusted (by dasymetric
redistribution) results in our experiments.

##  4 Pure Census Disaggregation (PCD)

Existing methods require preprocessing the data before estimating the model.
We propose to formalize the disaggregation problem in the context of
aggregated output learning [16]. We call this formulation PCD. Aggregated
output learning is a family of methods where the response information (the
value to be predicted) about each data point is unknown. However, for specific
subsets of data points, the sum of their responses is known. The
disaggregation problem can be converted into an aggregated output learning
problem: each pixel is a data point, and the subsets of data points for which
we have the sum of the responses are AUs’ pixels (the aggregated responses are
the census counts).

In our context, solving the disaggregation problem consists of learning a
function $f$ that, when applied and aggregated over all pixels of an AU,
predicts the census count of this AU with accuracy. More precisely, we attempt
to find the parameters $\bm{\theta}$ of a model $f_{\bm{\theta}}$ such that a
given error metric ($e_{PPE}$, $e_{PPSE}$, $e_{RMSE}$, etc.) is minimized:

$\operatorname*{arg\,min}_{\bm{\theta}}error(\operatorname{census},\operatorname{prediction}_{\bm{\theta},\bm{D}},\operatorname{surface},\mathcal{U}),$

where $\operatorname{prediction}_{\bm{\theta},\bm{D}}$ is a function
$\mathcal{U}\rightarrow\Re^{+}$ in the following form:

$\operatorname{prediction}_{\bm{\theta},\bm{D}}(u)=\sum_{p\in\operatorname{pixels}(u)}p_{w}\cdot
f_{\bm{\theta}}(\bm{D}_{p_{x},p_{y}}).$

Unlike the alternative techniques discussed previously, the error function is
plugged into the PCD problem being globally minimized by gradient descent, in
order to find $\bm{\theta}$. Other techniques are applicable to solve
aggregated output learning problems [16] and can be applied to PCD. The
prediction function for a unit is not based on a summary; it is composed of
individual and unsupervised predictions for all pixels.

The method is semi-supervised: while information is given to the method about
administrative units of the census, no indication about the population is
available at the individual pixel level. We use the name PCD-LinExp to
indicate the usage of LinExp in the PCD framework.

####  4.0.1 An interpretable model

We decide to focus on a regression model that we call LinExp. By design, it is
easily interpretable. Remember the model can be written as:

$f_{(\bm{a},b,c)}(\bm{D}_{p_{x},p_{y}})=\max(0,\exp(\bm{a}^{T}\bm{D}_{p_{x},p_{y}}+b)+c),$

where $\bm{a}$ is a vector containing one parameter for each covariate, whose
scalar product with the vector of covariate at position ($p_{x},p_{y}$) is
given to the $\exp$ function; $b$ and $c$ are two scalars allowing correcting
the output. The output is then passed into a function $\max(0,\bullet)$ (a
ReLU [11], here used outside its traditional scope) to ensure the output is
positive. The presence of the exp function is motivated by the distribution of
the census.

As the covariates in the tensor $\bm{D}$ have been standardized (they are
centered and have unit variance), the vector $\bm{a}$ is directly
interpretable as the importance of each covariate:

  * •

a high absolute value indicates that the covariate contributes highly to the
prediction, while a low absolute value indicates its contribution is small;

  * •

positive values indicate that each increase of the covariate increases the
predicted population, while negative values do the opposite.

The model is convex; however, when plugged into an error function such as
$e_{RMSE}$ or $e_{PPE}$, it becomes non-convex, mainly due to the sum of
absolute (or squared) values of exponentials. Using a (sub-)gradient descent
algorithm minimizing the error function on the training set is then a simple
yet effective choice to find good values for parameters $\bm{a},b,c$.

The function to be optimized for PCD-LinExp, given that $e_{PPE}$ is used, is

$\operatorname*{arg\,min}_{\bm{a},b,c}\frac{1}{\sum_{u\in\mathcal{U}}\operatorname{surface}(u)}\sum_{u\in\mathcal{U}}\\\
\left(\operatorname{surface}(u)\cdot\frac{|(\sum_{p\in\operatorname{pixels}(u)}p_{w}\cdot\max(0,\exp(\bm{a}^{T}\bm{D}_{p_{x},p_{y}}+b)+c))-\operatorname{census}(u)|}{\operatorname{census}(u)}\right).$

####  4.0.2 Optimization

Due to the presence of the $\exp$ function in the model, it is quite sensitive
to sudden erroneous increases in the vector $\bm{a}$, which can make the
output explode and even overflow. Hence, relatively small learning rates are
used, and we choose to avoid the use of stochastic gradient descent (SGD): the
whole dataset is computed at each iteration. This is allowed by the relatively
small size of our dataset; bigger datasets (not fitting entirely in memory)
would require using SGD.

Optimization is done using the Adam algorithm [12]. Adam, for Adaptive moment
estimation, is a variation of the standard Stochastic Gradient Descent
algorithm; it uses on first-order gradients and estimations of the first and
second moments of these gradients to regularize learning steps.

We use learning rates of $0.01$ for vectors $a$ and $b$, and of $0.001$ for
vector $c$. These values result from manual tuning. We run $1000$ iterations
of Adam before finishing and returning the results. Experiments show that
convergence is most of the time reached before $\sim 500$ iterations.

####  4.0.3 Initialization

Concerning initialization, we attempt to have at the first iterations a small
value for $\bm{a}^{T}\bm{x}+b$, as a large value may produce very important
gradients, which is not desired. Experiments have shown that having
$\bm{a}=\bm{0}$ and a negative value for $b$ (we chose $-4$) works well in
practice. Indeed, it greatly underestimates the population during the first
few iterations of the gradient descent.

####  4.0.4 Experimental Setup

In order to evaluate correctly the different methods (PCD-LinExp, RF)
presented earlier, we implemented a common test framework in Python. The
PyTorch library [17] was chosen to create the PCD-LinExp models. The RF
implementation relies on the library provided in the first publication of the
method (i.e. the randomforests library of the R language). The new RF
implementation has been tested in similar conditions than in the original
paper and gives similar results (in fact, we obtain better results, which is
probably due to small differences in the data).

We used standard 10-fold cross-validation to test methods. Cambodia is divided
into 193 districts, which are themselves divided into 1621 communes. As we
have to test the dasymetric redistribution, we need to keep communes inside
the same district into the same fold. We thus obtain folds with 19 districts
(three folds have 20 districts). As far as the cross-validation works is
concerned, a fold is used as test set, and the other nine folds are used as
training set. The fold used to test the models is then changed. We then obtain
10 measures of the generalization of the methods, from which we take the mean.
The measures chosen are the Per-Pixel-Error (PPE), the Per-Pixel-Squared-Error
(PPSE), and the Root Mean Square Error (RMSE).

####  4.0.5 RF-adj

In order to produce more accurate comparisons between the different methods,
we introduce a small variation of the RF method, that we call RF-adjusted
(abbreviated as RF-adj). The motivation behind RF-Adj is that RF tends to
over-predict population by a fair margin, as it is not constructed to run on
unadjusted errors (it was originally described with dasymetric redistribution
included in the method [20]).

RF-Adj is a simple adjustment that minimizes the error on the training set.
Given a trained $RF$ function, the $RF_{ajd}$ function is defined as
$RF_{adj}(\bm{D}_{p_{x},p_{y}})=c\cdot RF(\bm{D}_{p_{x},p_{y}}),$ where $c$ is
the mean factor of error for each unit in the training set:

$c=\frac{\sum_{u\in\mathcal{U}_{\text{training}}}\frac{\operatorname{prediction}_{\text{RF}}(u)}{\operatorname{census}(u)}}{|\mathcal{U}_{\text{training}}|}.$

Experiments show that RF-Adj generates smaller unadjusted errors than RF.

##  5 Results and discussion

Fig 1 (at the beginning of the paper) compares the population maps generated
by the PCD-LinExp-PPE and RF methods, before dasymetric redistribution.
Results obtained by PCD-LinExp-RMSE and RF-Adj are not included as there are
visually similar to PCD-LinExp-PPE and RF, respectively. The two methods
identify similar regions of high population density, notably near the capital
Phnom Penh (at the south-east of the country) and around the lake Tonlé Sap,
in the center of the country. Visually, PCD-LinExp-PPE produces a seemingly
more detailed map than RF, with less “saturation”. On the one hand, PCD-
LinExp-PPE highlights the role of the road infrastructure. On the other hand,
RF produces circles of high population density around hospitals and roads.

Table 2 shows the error metrics of the methods described in our experimental
setup. Box plots containing the data of the individual folds are depicted in
Fig 4. Running a paired t-test (null hypothesis mean(A-B) = 0, and a p-value
limit of $0.05$) on the fold errors reveals that LinExp-PPE and LinExp-RMSE
significantly less PPE/PPSE than RF and RF-Adj on unadjusted metrics. LinExp-
RMSE is significantly better than all the other methods on unadjusted RMSE. On
adjusted metrics, no method is significantly better than another, but LinExp-
PPE that is dominated by all method on the adjusted RMSE metric.

PCD-LinExp-PPE and PCD-LinExp-RMSE better perform when using the unadjusted
PPE/PPSE metric; PCD-LinExp-RMSE generates better results for the unadjusted
RMSE metric. The PCD-LinExp unadjusted outcomes are comparable to the adjusted
outcomes obtained after dasymetric redistribution, which strengthens our claim
that PCD-LinExp can be used to predict the population where census data are
not available.

Table 2: Error metrics Method | Unadjusted  | Adjusted
---|---|---
| PPE | PPSE | RMSE  | PPE | PPSE | RMSE
Areal weighting  | N/A | N/A | N/A  | $87.90\ \%$ | $215.46\ \%$ | $6777.71$
PCD-LinExp-PPE  | $45.43\ \%$ | $45.46\ \%$ |  $6202.23$ | $39.75\ \%$ | $34.19\ \%$ | $4407.64$
PCD-LinExp-RMSE  | $54.20\ \%$ | $78.35\ \%$ |  $3889.35$ | $39.84\ \%$ | $34.45\ \%$ | $3335.68$
RF  | $178.07\ \%$ | $914.97\ \%$ |  $6584.09$ | $39.74\ \%$ | $32.44\ \%$ | $3175.68$
RF-Adj  | $127.50\ \%$ | $512.72\ \%$ |  $4742.39$ | $39.68\ \%$ | $32.24\ \%$ | $3199.69$
Figure 4: Scatter plot of the error metrics for each model and for each fold.
While it is difficult to compare results for the redistributed errors, PCD-
LinExp-PPE and PCD-LinExp-RMSE give the best result on unadjusted errors. The
methods PCD-LinExp-PPER, PCD-LinExp-RMSER, RFR and RF-AdjR are redistributed
counterparts of the original methods.

Fig 5 shows the unadjusted PPE error obtained for each unit as a function of
the surface of the unit. PCD-LinExp-PPE and PCD-LinExp-RMSE generate smaller
errors. The errors obtained under PCD-LinExp-PPE are globally independent of
the unit surface. This is because the PPE metric reduces the influence of the
surface. As expected, other methods that do not account for the surface
generate greater errors for less-populated, high-surface units. Smaller errors
are obtained for highly populated, small-surface units as this would greatly
impact the RMSE.

Figure 5: Error for each unit for all models.

Overall, all results presented in this section indicate that PCD-LinExp-PPE
and PCD-LinExp-RMSE generate better results than RF and RF-Adj, at least for
predicting the population (i.e. on unadjusted errors). The comparison for
population disaggregation (using adjusted errors) is more difficult, as all
experiments indicate ties. We explained previously that using dasymetric
redistribution breaks the training set/test set dichotomy, and thus should be
avoided. However, this demonstrates that using the aggregate learning
framework with a very simple model in this problem outperforms/is on par
(depending on the metric) with more complex models not using aggregate
learning.

####  5.0.1 Acknowledgments

Computational resources have been provided by the supercomputing facilities of
the Université catholique de Louvain (CISM/UCL) and the Consortium des
Équipements de Calcul Intensif en Fédération Wallonie Bruxelles (CÉCI) funded
by the Fond de la Recherche Scientifique de Belgique (F.R.S.-FNRS) under
convention 2.5020.11. We would like to thank Pavel Demin and the CP3 group
that shared with us part of their reserved resources. The second and third
authors acknowledge financial support from the ARC convention on “New
approaches to understanding and modeling global migration trends” (convention
18/23-091).

## References

  * [1] Breiman, L.: Random forests. Machine Learning 45(1), 5–32 (Oct 2001). https://doi.org/10.1023/A:1010933404324
  * [2] Briggs, D.J., Gulliver, J., Fecht, D., Vienneau, D.M.: Dasymetric modelling of small-area population distribution using land cover and light emissions data. Remote Sensing of Environment 108(4), 451 – 466 (2007). https://doi.org/https://doi.org/10.1016/j.rse.2006.11.020
  * [3] Center for International Earth Science Information Network - CIESIN - Columbia University: Gridded population of the world, version 4 (gpwv4): Population density, revision 10 (20180711 2017), https://doi.org/10.7927/H4DZ068D
  * [4] Center for International Earth Science Information Network - CIESIN - Columbia University: U.s. census grids 2010 (summary file 1) (20180719 2017), https://doi.org/10.7927/H40Z716C
  * [5] Dmowska, A., Stepinski, T.F.: High resolution dasymetric model of u.s demographics with application to spatial distribution of racial diversity. Applied Geography 53, 417 – 426 (2014). https://doi.org/https://doi.org/10.1016/j.apgeog.2014.07.003
  * [6] Doupe, P., Bruzelius, E., Faghmous, J., Ruchman, S.G.: Equitable development through deep learning: The case of sub-national population density estimation. In: Proceedings of the 7th Annual Symposium on Computing for Development. pp. 6:1–6:10. ACM DEV ’16, ACM, New York, NY, USA (2016). https://doi.org/10.1145/3001913.3001921
  * [7] Eicher, C.L., Brewer, C.A.: Dasymetric mapping and areal interpolation: Implementation and evaluation. Cartography and Geographic Information Science 28(2), 125–138 (2001)
  * [8] Flowerdew, R., Green, M.: Developments in areal interpolation methods and gis. In: Geographic Information Systems, Spatial Modelling and Policy Evaluation, pp. 73–84. Springer (1993)
  * [9] Gallego, F.J.: A population density grid of the european union. Population and Environment 31(6), 460–473 (Jul 2010). https://doi.org/10.1007/s11111-010-0108-y
  * [10] Goodchild, M.F., Anselin, L., Deichmann, U.: A framework for the areal interpolation of socioeconomic data. Environment and planning A 25(3), 383–397 (1993)
  * [11] Hahnloser, R.H., Sarpeshkar, R., Mahowald, M.A., Douglas, R.J., Seung, H.S.: Digital selection and analogue amplification coexist in a cortex-inspired silicon circuit. Nature 405(6789), 947 (2000)
  * [12] Kingma, D.P., Ba, J.: Adam: A method for stochastic optimization. CoRR abs/1412.6980 (2014), http://arxiv.org/abs/1412.6980
  * [13] LeCun, Y., Bottou, L., Bengio, Y., Haffner, P.: Gradient-based learning applied to document recognition. Proceedings of the IEEE 86(11), 2278–2324 (1998)
  * [14] Mennis, J.: Generating surface models of population using dasymetric mapping. The Professional Geographer 55(1), 31–42 (2003)
  * [15] Monmonier, M.S., Schnell, G.A.: Land use and land cover data and the mapping of population density. The International Yearbook of Cartography 24(115), e121 (1984)
  * [16] Musicant, D.R., Christensen, J.M., Olson, J.F.: Supervised learning by training on aggregate outputs. In: Seventh IEEE International Conference on Data Mining (ICDM 2007). pp. 252–261. IEEE (2007)
  * [17] Paszke, A., Gross, S., Chintala, S., Chanan, G., Yang, E., DeVito, Z., Lin, Z., Desmaison, A., Antiga, L., Lerer, A.: Automatic differentiation in pytorch (2017)
  * [18] Robinson, C., Hohman, F., Dilkina, B.: A deep learning approach for population estimation from satellite imagery. In: Proceedings of the 1st ACM SIGSPATIAL Workshop on Geospatial Humanities. pp. 47–54. ACM (2017)
  * [19] Simonyan, K., Zisserman, A.: Very deep convolutional networks for large-scale image recognition. arXiv preprint arXiv:1409.1556 (2014)
  * [20] Stevens, F.R., Gaughan, A.E., Linard, C., Tatem, A.J.: Disaggregating census data for population mapping using random forests with remotely-sensed and ancillary data. PLOS ONE 10(2), 1–22 (02 2015). https://doi.org/10.1371/journal.pone.0107042, https://doi.org/10.1371/journal.pone.0107042
  * [21] Tian, Y., Yue, T., Zhu, L., Clinton, N.: Modeling population density using land cover data. Ecological Modelling 189(1-2), 72–88 (2005)
  * [22] Tobler, W.R.: Smooth pycnophylactic interpolation for geographical regions. Journal of the American Statistical Association 74(367), 519–530 (1979)
  * [23] UN Economic and Social Council: Resolution adopted by the economic and social council on 10 june 2015 (2020 world population and housing census programme) (August 2015), http://www.un.org/ga/search/view_doc.asp?symbol=E/RES/2015/10
  * [24] Wright, J.K.: A method of mapping densities of population: With cape cod as an example. Geographical Review 26(1), 103–110 (1936)