# Fine-grained Population Mapping from Coarse Census Counts and Open Geodata

Nando Metzger  ETH Zurich, Switzerland  John E. Vargas-Muñoz  École
Polytechnique Fédérale de Lausanne, Switzerland  Rodrigo C. Daudt  ETH Zurich,
Switzerland  Benjamin Kellenberger  Yale University, Untited States  École
Polytechnique Fédérale de Lausanne, Switzerland  Thao Ton-That Whelan
International Committee of the Red Cross, Switzerland  Ferda Ofli  Hamad Bin
Khalifa University, Qatar  Muhammad Imran  Hamad Bin Khalifa University, Qatar
Konrad Schindler  ETH Zurich, Switzerland  Devis Tuia  École Polytechnique
Fédérale de Lausanne, Switzerland

###### Abstract

Fine-grained population maps are needed in several domains, like urban
planning, environmental monitoring, public health, and humanitarian
operations. Unfortunately, in many countries only aggregate census counts over
large spatial units are collected, moreover, these are not always up-to-date.
We present Pomelo, a deep learning model that employs coarse census counts and
open geodata to estimate fine-grained population maps with 100 m ground
sampling distance. Moreover, the model can also estimate population numbers
when no census counts at all are available, by generalizing across countries.
In a series of experiments for several countries in sub-Saharan Africa, the
maps produced with Pomelo are in good agreement with the most detailed
available reference counts: disaggregation of coarse census counts reaches
$R^{2}$ values of 85–89%; unconstrained prediction in the absence of any
counts reaches 48–69%.

## Introduction

High-resolution population maps are crucial for many planning tasks, from
urban planning [1] to preparing humanitarian actions [2] and effective
disaster response [3]. Given the rapid population growth in many regions of
the world [4] and the increasing rate at which populations shift in response
to environmental and social changes, it is important to maintain accurate, up-
to-date maps. Unfortunately, census data are often only available at very
coarse spatial resolution (e.g., one aggregate number for a district with
hundreds or even thousands of km2) and therefore not suitable as a basis for
local planning: whether for sustainable land use and infrastructure management
or for targeted disaster relief, planners need to know in more detail where
the people are. The problem is especially prevalent in developing countries in
the global south, where humanitarian actions are more often needed yet census
data availability and quality are limited.

Remote sensing products and other openly available geographical datasets like
OpenStreetMap (OSM) can serve as auxiliary, high-resolution evidence to create
fine-grained population density maps [5]. Yet, the design of effective
population density models [6] that combine such data sources with low-
resolution census counts remains a challenge. Generally speaking, two
different approaches have been employed for population mapping [7]: bottom-up
and top-down. Bottom-up methods [8, 7] start from local surveys of population
density, collected at a number of sample locations, and attempt to generalize
from detailed but sparse samples to the unobserved regions to cover larger
areas. Researchers have proposed different ways to locally measure population
density, such as counting the (average) number of people per rooftop area [8,
7] or, if more resources are available for the local survey, specific average
densities for different types of residential zones (urban-, rural-, and non-
residential) [7, 9]. A main drawback of bottom-up methods is that local
surveys will necessarily remain extremely sparse and can hardly provide enough
data points to scale population mapping up to the country level. On the
contrary, top-down approaches [6, 10] rely on census data, which ensures
complete coverage at the expense of much lower spatial resolution, in some
cases down to a single head count per large district. The task then becomes to
disaggregate that data to a much finer resolution, often a regular grid.

Top-down approaches [6, 11] commonly use dasymetric disaggregation to
redistribute the known, spatially coarse population counts for census areas on
the order of many km2 across smaller spatial units [12] – for instance square
blocks of size 100$\times$100 meters – with the help of auxiliary variables
that covary with population density. Examples of such covariates for
population disaggregation are the presence of buildings [13], building counts
[10], building volumes [14], and cumulative road lengths [15]. Dobson et al.
[16] propose a weighted combination of several covariates, such as land cover
and proximity to roads, to compute what portion of the population to assign to
each target unit. Stevens et al. [6] resort to machine learning and train a
random forest model to predict population density from a set of covariates
like building maps and night light images. A limitation of the method proposed
by Stevens et al. [6] is that administrative regions with known census counts
are used as training units, which means that the input features (covariates)
must be aggregated across all pixels in such a region; whereas at prediction
(inference) time, density values are predicted for each individual pixel,
respectively feature vector. Consequently, the method’s efficacy is reduced in
countries with coarse census units, as aggregation over large spatial regions
induces distribution shifts between the training and test data.

Recent work has employed deep learning, in particular convolutional neural
networks (CNNs [17]) for population estimation, in an attempt to better
account for spatial context. Gervasoni et al. [18] used a CNN to map data
extracted from OSM to disaggregation weights per pixel with 200$\times$200 m
ground sampling distance (GSD). This allows one to super-resolve census counts
to a finer granularity; however, since the predicted weights are relative, the
method cannot be used to predict population numbers for regions without census
data. Jacobs et al. [19] trained a CNN to predict fine-grained population maps
from very high resolution (VHR) optical satellite imagery and census data.
This is viable for individual cities, but difficult to scale up to larger
geographic contexts, due to the limited availability and high cost of VHR
images. Both CNN-based methods [18, 19] depend on fairly fine-grained initial
counts with census units of few km2, e.g., urban census blocks in France,
respectively the USA. For comparison, the administrative units for which
census data are available in much of Africa have sizes on the order of several
hundred to several thousand km2, which precludes the use of such methods.

Here, we propose a methodology to estimate fine-grained population maps from
such very coarse census data by fusing them with covariate maps of higher
resolution. As inputs, we use publicly available products derived from remote
sensing images and other open data sources, e.g., building footprint maps,
night light images, and OSM road layers. They are used to train a model that
predicts population on a spatial grid with 100$\times$100 m GSD. Our method,
which we call Pomelo (short for “population mapping by estimation of local
occupancy rates”), is inspired by work on _guided super-resolution_ of low-
resolution images [20, 21, 22]. In an experimental evaluation with data from
three different countries in Sub-Saharan Africa (Tanzania, Zambia, Mozambique)
Pomelo delivers significantly more accurate population maps compared to
several baseline methods, including the pioneering work of Stevens et al. [6]
with which the widely used WorldPop maps [23] are created. Moreover, Pomelo
can not only disaggregate existing census data, but can also predict
population maps in the absence of census counts. Consequently, it can be
deployed to regions or countries where no suitable census information is
available – of course these estimates, inevitably, have higher uncertainties,
as the total population is no longer constrained by a known aggregate number.
In summary, Pomelo provides knowledge about where people are at the hectometer
scale, even if census data are not available.

## Results

We compare Pomelo with other state-of-the-art baseline models using three
different evaluation strategies, and three performance metrics (see Methods
section), namely $R^{2}$, mean absolute error (MAE), and mean absolute
percentage error (MAPE):

Coarse supervision: In this strategy, used by several works in the literature
[6, 24, 25], the census data at the finest available level are reserved
exclusively for performance evaluation. First, the available census data are
artificially coarsened by aggregating them to the next-higher (“coarse”) level
of administrative regions (e.g., counts at district level are aggregated into
a single count per province). Then, the coarse level regions are divided into
five random folds, of which three serve as training set, one as validation set
for hyper-parameter tuning and checkpointing the model on the best MAPE, and
one as test set to measure the model’s predictive skill. Finally, the model is
trained to predict gridded population numbers at 100 m GSD, such that they add
up to the coarse level counts. To measure performance, the gridded population
estimates are aggregated to obtain population numbers for each _fine_ level
administrative region in the test set and are then compared to the
corresponding census data.

This scheme corresponds to an application scenario where fine-level census
data are not available: all model fitting and mapping is only based on the
coarse level. The fine level counts are employed exclusively to evaluate
prediction performance.

The entire training and evaluation procedure is run as a five-fold cross-
validation, rotating the validation and test folds. Figure 1 depicts the
coarse census data used for training, the fine level census data used for
evaluation, and the fine-grained population maps obtained by our model, for
the surroundings of Zanzibar, Tanzania.

(a) Coarse (district level)

(b) Fine (ward level)

(c) Pixel (prediction level);
100 m resolution

Figure 1: Overview of the different level of population data. Shown location:
Zanzibar City, Tanzania. The map visualizations were created in QGIS 3.14
[26].

Fine supervision: In the second evaluation scheme, census data at the finest
available level are used for supervision. The data are _not_ coarsened, rather
the original fine-level data are divided into five random folds. These folds
are split along the same coarse-level administrative boundaries as above, such
that the training, validation and test portions during cross-validation are
identical to the coarse supervision scenario. The model is trained to predict
gridded population numbers at 100 m GSD, such that they add up to the fine-
level counts. Just like above, the gridded population estimates are aggregated
to obtain population numbers for each fine-level administrative region in the
held-out test fold (e.g., Figure 1(b)), which are then compared to the
corresponding census data to compute performance metrics.

This scheme corresponds to the actual, realistic population mapping task with
the available census granularity. The finest available resolution is used to
learn the best possible model, geographically separating training and testing
regions to ensure an unbiased evaluation. Note that the combined test set,
after holding out each fold once for testing, is the same for coarse
supervision and fine supervision, such that the error metrics are directly
comparable.

Transfer task: Besides census disaggregation, Pomelo can also be used to
estimate population density for a given country in the absence of census
numbers. To evaluate the performance on such more challenging task, we use
data from seven different countries. One of the countries is held out as test
set, and the data from the remaining six countries are randomly split into 80%
for training and 20% for validation (i.e., hyper-parameter tuning). Then, the
model is trained to predict gridded population numbers at 100 m GSD, using the
_fine-level_ strategy. Finally, the model is deployed for the held-out
country, and evaluated in the same way as above.

For each test country we run five models with different random parameter
initialization and train/validation splits and report the average and standard
deviation of each performance metric. We compare our proposed method with
other four baseline methods, described in the Methods section: Building count
disaggregation, the random forest model used by WorldPop [6], a Markov Random
Field based method for population disaggregation, and a Convolutional Neural
Network model.

For evaluation, we use data from Tanzania and Zambia with the three
aforementioned evaluation strategies. Mozambique is used only in the transfer
task evaluation. The data available for that country is not suitable for the
other two evaluation schemes: in Mozambique there are only 413 fine-level
regions grouped into 156 coarse-level regions. The low aggregation factor of
only $2.6\times$ causes the performance metrics to saturate at a similar value
of 83% $R^{2}$ for all methods, while giving little indication about the
correctness of the pixel-level maps.

### Results – coarse supervision

Table 1 presents the performance metrics for the coarse supervision strategy
in Tanzania and Zambia. The first two rows show learning-free methods that are
only capable of pure disaggregation, i.e., they cannot estimate population
numbers for regions without known (aggregate) census counts. The MRF
formulation brings a marked improvement compared to the simple disaggregation
scheme based only on building counts. For the three learning-based methods, we
first predict population numbers without using the target regions’ census
counts, then normalize those counts to relative fractions and use those as
weights for dasymetric disaggregation. Pomelo achieves the best performance in
all three metrics, with the closest competitor being the MRF. An interesting
observation is that the region-based RF method [6], when trained with such
coarse census data, does not even reach the performance of simple building
count disaggregation. The evaluation results for Zambia largely mirror those
for Tanzania, where the MRF is the closest to Pomelo in terms of $R^{2}$ and
MAE. Next follows the CNN, which handles Zambia much better than Tanzania.

Table 1: Performance with coarse supervision for Tanzania and Zambia. Evaluation Set | Method |  $R^{2}\uparrow$ [%] |  MAE $\downarrow$ |  MAPE $\downarrow$ [%]
---|---|---|---|---
Tanzania | Building disaggregation | 65.2 | 3’700 | 23.1
MRF | 78.6 | 3’300 | 21.7
RF per region [6] |  49.0 $\pm$ 3.0 |  4’380 $\pm$ 050 |  26.9 $\pm$ 0.2
CNN |  69.5 $\pm$ 2.0 |  4’200 $\pm$ 100 |  28.8 $\pm$ 0.5
Pomelo |  85.7 $\pm$ 0.9 |  3’100 $\pm$ 040 |  21.6 $\pm$ 0.2
Zambia | Building disaggregation | 82.1 | 4’300 | 53.9
MRF | 87.0 | 3’900 | 51.3
RF per region [6] |  76.4 $\pm$ 0.0 |  4’750 $\pm$ 090 |  60.9 $\pm$ 1.2
CNN |  83.9 $\pm$ 0.7 |  4’250 $\pm$ 040 |  51.4 $\pm$ 1.7
Pomelo |  87.9 $\pm$ 0.1 |  3’730 $\pm$ 030 |  50.4 $\pm$ 0.4

### Results – fine supervision

Table 2 shows quantitative performance when the models are trained and
evaluated on the finer census level. Note that we cannot use the learning-free
building disaggregation and MRF baselines with this evaluation strategy, since
they are only designed for disaggregation of regions with known overall
population numbers (which makes no sense in the fine supervision setting, as
the target quantity would have to be known in advance). For Tanzania, we again
observe the best performance with Pomelo, although the advantage is smaller
than in the case of coarse supervision (Table 1). The performance of the
region-based RF improves a lot when trained with fine supervision, showing
that the coarse supervision does not offer enough supervision signal for
training at the region level. The second section of Table 2 shows results for
Zambia. Pomelo maintains a slight edge in terms of $R^{2}$ and absolute error.
It does not fare as well in terms of MAPE, which is due to large relative
errors in few, sparsely populated regions (that nonetheless translate to small
errors of the absolute population count): just by removing Zambia’s smallest
census region with only 28 inhabitants from the evaluation, the MAPE of Pomelo
drops to 45%.

Figure 2 shows scatter plots (census counts versus our predicted counts) for
the fine-level administrative regions of both Tanzania and Zambia. In both
cases, the data are, with few exceptions, close to the (red, dashed) diagonal
that corresponds to the ideal result, where predictions and ground truth
coincide.

Table 2: Performance with fine supervision for Tanzania and Zambia. Evaluation Set | Method |  $R^{2}\uparrow$ [%] |  MAE $\downarrow$ |  MAPE $\downarrow$ [%]
---|---|---|---|---
Tanzania |  RF per region [6] |  79.1 $\pm$ 2.6 |  3320 $\pm$ 70 |  22.1 $\pm$ 0.3
CNN |  80.8 $\pm$ 2.1 |  3290 $\pm$ 80 |  21.7 $\pm$ 0.4
Pomelo |  87.6 $\pm$ 0.2 |  2890 $\pm$ 20 |  20.4 $\pm$ 0.3
Zambia |  RF per region [6] |  87.6 $\pm$ 0.0 |  3740 $\pm$ 40 |  46.5 $\pm$ 0.4
CNN |  88.2 $\pm$ 0.4 |  3680 $\pm$ 50 |  44.9 $\pm$ 0.4
Pomelo |  88.7 $\pm$ 0.3 |  3650 $\pm$ 20 |  48.0 $\pm$ 0.6

(a) Tanzania

(b) Zambia

Figure 2: Comparison between the estimates predicted by Pomelo and the
respective ground truth census in the fine level training task.

### Transfer Task

We also evaluate Pomelo in a scenario where no census data are available for
the target country. Here, the models must instead be trained on other
countries and should be able to generalize to a potentially different unknown
geographical context. From the seven countries considered (Tanzania, Zambia,
Mozambique, Rwanda, Uganda, Democratic Republic of Congo, and Nigeria), we
select six for training and validation and report the evaluation metrics on
the respective held out country (Tanzania, Zambia, and Mozambique) in Table 3.

The first rows of each section in the table (“average occupancy rate”) show
the results of a naïve baseline that computes the average occupancy rate in
the six training countries and multiplies it with the building count map of
the target country to obtain the population map. One can see that this simple
approach does not perform well in any of the countries, barely reaching a
positive $R^{2}$ score only in Zambia. The region-based RF model does not
generalize well either, to the point that the $R^{2}$ metric becomes negative
and the MAPE exceeds 100% for Zambia and Mozambique. The performance of the
CNN is still rather good across all evaluation settings, despite the more
challenging scenario. Pomelo exhibits even better overall performance, with
the exception of a slightly higher MAE in Mozambique.

As a proxy for a “general” model that is valid for an entire geographic area,
we also create a training set with data from all seven countries (including
the respective target country), and apply the resulting model to the (unseen)
test set of the target country, using five-fold cross-validation to cover the
whole country. For comparison, we also run models trained specifically for
each country, but without dasymetric rescaling, to keep the comparison fair.
That scenario would normally not be relevant in practice, since the training
implies that aggregate counts are available except in the special situation
that no current census counts are available and a previously trained model is
reused. However, by removing the influence of the postprocessing it shows the
ability of the country-specific model to predict unconstrained, absolute
numbers. As expected, the country-specific models perform better than the one
learned only with data from four other countries, but not as well as a model
trained on all seven countries, presumably due to the larger amount and
variability of training samples.

Table 3: Performance on the transfer task for all countries. Evaluation Set | Method |  $R^{2}\uparrow$ [%] |  MAE $\downarrow$ |  MAPE $\downarrow$ [%]
---|---|---|---|---
Tanzania | Average occupancy rate |  0-52.9 $\pm$ 00.0 |  13’200 $\pm$ 0’000 |  098.0 $\pm$ 00.0
RF per region [6] |  -022.1 $\pm$ 010.1 |  08’000 $\pm$ 0’700 |  066.5 $\pm$ 006.2
CNN |  -034.6 $\pm$ 036.2 |  06’800 $\pm$ 0’300 |  049.6 $\pm$ 003.0
Pomelo |  -068.1 $\pm$ 030.9 |  05’200 $\pm$ 0’200 |  035.6 $\pm$ 001.0
Pomelo – 7 countries |  -081.7 $\pm$ 030.6 |  03’170 $\pm$ 0’020 |  021.7 $\pm$ 000.6
|  Pomelo – Tanzania only |  -082.0 $\pm$ 031.5 |  03’800 $\pm$ 0’070 |  026.4 $\pm$ 000.2
Zambia | Average occupancy rate |  -014.4 |  12’300 $\pm$ 0’000 |  150.2 $\pm$ 000.0
RF per region [6] |  -193.4 $\pm$ 361.3 |  19’000 $\pm$ 8’000 |  271.9 $\pm$ 139.0
CNN |  -059.0 $\pm$ 001.6 |  06’600 $\pm$ 0’300 |  082.7 $\pm$ 005.6
Pomelo |  -068.7 $\pm$ 003.1 |  06’100 $\pm$ 0’200 |  075.6 $\pm$ 001.6
Pomelo – 7 countries |  -085.4 $\pm$ 000.5 |  03’810 $\pm$ 0’010 |  049.6 $\pm$ 000.2
|  Pomelo – Zambia only |  -074.0 $\pm$ 001.4 |  04’820 $\pm$ 0’070 |  047.6 $\pm$ 000.1
Mozambique | Average occupancy rate |  0-88.1 $\pm$ 000.0 |  82’500 $\pm$ 0’000 |  153.0 $\pm$ 000.0
RF per region [6] |  0-49.5 $\pm$ 021.1 |  59’000 $\pm$ 3’000 |  102.0 $\pm$ 003.6
CNN |  -043.6 $\pm$ 035.2 |  35’200 $\pm$ 1’300 |  063.0 $\pm$ 001.9
Pomelo |  -048.3 $\pm$ 003.2 |  36’900 $\pm$ 0’900 |  047.5 $\pm$ 000.6
Pomelo – 7 countries |  -083.0 $\pm$ 000.5 |  16’800 $\pm$ 0’200 |  034.1 $\pm$ 000.2
|  Pomelo – Mozambique only |  -061.0 $\pm$ 003.0 |  30’600 $\pm$ 0’900 |  045.0 $\pm$ 001.0

### Visual results

The finest census level available for validation are administrative regions
with areas of a few km2 to several hundred km2, and ranging from very low to
very high population densities. To visually examine the predicted population
maps at 100 m resolution, we select regions with known low/high densities and
resort to OSM for verification. Figure 3 presents example population
estimates, shown as heat maps overlaid on OSM for regions around the city of
Dar-es-Salaam in Tanzania. Each example covers an area of 400$\times$400 m,
corresponding to 16 cells in our gridded map. Figure 3a shows an example of a
very low population estimate in a rural area. The estimates in low-density
residential areas (Figure 3b) are somewhat higher and more variable, due to
varying numbers of buildings. In high-density residential areas (top part of
Figure 3c), the estimates are even higher. For non-residential areas with
mostly commercial buildings (bottom part of Figure 3c), the model predicts a
relatively low and uniform density. These variations can occur over a small
spatial distance: as an example, Figure 3c show a sudden drop in population
density between immediately adjacent city areas.

Figure 3: Visual examples of population estimates on a 100 m grid, represented
as heat maps on OSM background [27]. The color scale for the heat maps is
shown in the top right in absolute population counts.

Figure 4 depicts the predicted building occupancy rates as heat maps for two
neighborhoods in the city of Dar-es-Salaam, that are located close to each
other but feature very different occupancy rates. The densely populated area
of the Mtoni region in Figure 4a has an estimated (mean) building occupancy
rate of 4.2, whereas the area in the Kijichi region in Figure 4b has a much
lower estimated occupancy rate of 2.8.

Figure 4: Estimated average building occupancy rate, depicted as heat map on OSM background, for the Mtoni and Kijichi regions of Dar-es-Salaam. Note the significantly different occupancy rates in the marked areas. Figure 5: Population and building occupancy rate per ward around Dar-es-Salaam, Tanzania: a) estimated population counts, b) true population counts, c) estimated building occupancy rate, d) true occupancy rate, e) signed relative errors of population estimates. The map visualizations were created in QGIS 3.14 [26] |  |  |  |
---|---|---|---|---
a) Sentinel-2 RGB bands | b) HRPDM | c) RF per region |  d) Pomelo |
Figure 6: Visual comparison of the Sentinel-2 imagery [28] High Resolution
Population Density Maps (HRPDM) [29] per-region RF disaggregation [6] and
Pomelo. Shown location: Zanzibar City, Tanzania. The Sentinel-2 mosaic was
created with Google Earth Engine [30] and visualized in QGIS [26].

Figure 5 visualizes population numbers and building occupancy rates per
administrative region (“ward”) in the surroundings of Dar-es-Salaam. One can
see that wards have rather heterogeneous sizes and populations. Figure 5a
depicts our model’s population estimates, obtained by aggregating the grid
cells within each ward. Figure 5b depicts the reference values from the
national census. Figure 5c are our model’s estimated building occupancy rates,
computed by dividing the estimated population by the number of buildings – or,
equivalently, by averaging the predicted densities, with weights proportional
to the building counts in their corresponding grid cells. Figure 5d are the
true occupancy rates according to the census (i.e., the population count
divided by the number of buildings). The relative error of the per-ward
population estimates, as a percentage of the true numbers, is shown in Figure
5e. In general, the predictions are in good agreement with the actual numbers.
We do observe a trend that Pomelo overestimates the populations in built-up
low-density regions, while it tends to underestimate the extremely high
occupancies in the centers of large cities. Finally, Figure 6 shows a visual
comparison of recent population maps of Zanzibar City, at 100m resolution;
Including the High Resolution Population Density Maps (HRPDM) project [29],
Random Forest disaggregation per census region following [6] (but trained with
the same data as our method), and Pomelo. It is apparent that Pomelo recovers
the population distribution in more detail (i.e., with a higher effective
resolution), whereas the HRPDM maps as well as the region-wise disaggregation
result appear overly smooth and do not recover high-frequency variations of
the population density.

## Discussion

We evaluate the performance of our method using three different scenario:
_coarse supervision_ , _fine supervision_ , _transfer task_. In the past few
years, the first method has been the _de facto_ standard for top-down
population estimation [6, 25] and connects our evaluations to related
literature. However, we argue that the _fine supervision_ scenario is actually
more representative of the procedure one would use in practice: On the one
hand, it would seem unnatural to not use the finest available spatial units if
the goal is to produce the best possible population maps from a given census
dataset. On the other hand, in the _coarse supervision_ setting one predicts
population maps at a level that is already known from the census. The coarse
scenario seems to have arisen largely from the desire to evaluate pure
disaggregation methods, which requires access to ground truth counts at two
different levels of the spatial hierarchy.

It is worth noting that the region-level RF method of Stevens et al. [6] does
not work well if the census is only available for spatially coarse units,
which is a rather frequent situation in developing countries. We see several
possible reasons for this: first, there simply are fewer regions at coarser
hierarchy levels and one may be left with too few training examples to learn a
good model. Second, performing feature aggregation leads to a domain shift
between coarse-level training data and fine-level test data – in particular at
the pixel level.

We found it advantageous to first estimate the spatial distribution of the
building occupancy rate, and then compute population numbers by multiplying
the occupancy rate and the building count at a given location. This is in
contrast to methods that directly estimate absolute population numbers [6] or
relative population fractions (which are the weights for dasymetric
disaggregation [18]). The finding that factoring into building counts and
occupancy reduces the estimation errors lends support to two assumptions
implicit in our approach: First, the available maps of building counts
(respectively, building footprints) are apparently rather accurate, so that
using them directly rather than as one of many “soft” covariates reduces the
estimation error. Second, the occupancy rate appears to have lower spatial
variability than the population density, making it easier to estimate from the
same covariates.

The important role of country-scale building counts (respectively, building
footprint maps) also merits some discussion. For our work we obtain them from
free sources, namely the Open Buildings dataset [31] and the Grid Maps of
Building Patterns [32]. Both these datasets are created with computer vision-
based building detectors on the basis of high-resolution satellite imagery
(GSD <1 m). Although the datasets are of high quality, such large-scale maps
inevitably contain errors and data gaps, especially in rural areas [33]. We
have tried to maximize completeness by fusing the two building datasets, but
note that missing buildings (and to a lesser degree perhaps also spurious
buildings) may cause errors in our population maps, mainly in scarcely
populated areas.

Although the mentioned building footprint datasets are available for free, the
underlying high-resolution images used to produce them are not. The processing
of high-resolution data at country scale also requires considerable
computational resources. Our method thus critically depends on data that are
oftentimes primarily produced for a different purpose or for philanthropic
reasons, and over whose production we have no control. If at some point in the
future no up-to-date building dataset is at hand, one could try to resort to
autoregressive settlement growth models [34, 35]. Alternatively, one may keep
the building detection implicit and also feed high-resolution images to the
population estimator. Such approaches have been developed [7, 36], but since
they rely on (commercial) high-resolution satellite imagery they can hardly be
scaled up to entire countries, except by entities who could equally supply the
building counts. Crowdsourced data (e.g., OSM roads) have their own challenges
due to incompleteness and temporal inaccuracies. However, in several regions
of developing countries, commercial products and official data are not
available, making crowdsourced data the only source of information with
sufficient coverage in those regions [37].

In addition to the building counts, Pomelo is driven by several other
geospatial data layers that are publicly available and correlate with
population. We have analyzed feature importance with the permutation method
[38], which essentially measures the performance drop caused by randomly
shuffling a single input layer, so as to render that covariate uninformative.
Exemplary feature importance scores for Tanzania (fine supervision) are shown
in Figure 7. Although the ordering of the covariates fluctuates depending on
the chosen performance measure, we noted one commonality, namely the
nightlight and settlement layers are among the most predictive inputs.
Moreover, also features with low scores appear to carry some relevant
information, as explicit feature selection based on the importance scores
tends to harm model performance. For more details please refer to the Methods
section.

Figure 7: Feature importance analysis for the fine supervision setting in
Tanzania.

Apart from the covariates used in the present work, there are others that
could potentially be useful, and that one could obtain from open geospatial
data sources. For instance, online land use maps as available in OSM could
help to identify non-residential buildings such as schools or shopping
centers, and building height estimates from Esch et al. [39], once publicly
available, could help to predict more precise building occupancy rates.
Furthermore, data extracted from social media could possibly also support
population mapping in certain areas. For example, the Facebook Marketing API
[40] allows one to access (anonymized and aggregated) information about the
platform’s users at a certain location [41, 42], such as the number of active
users, the type of internet connection used, etc. It is quite possible that
such information is to some degree predictive of population density. However,
social media and OSM data are typically incomplete and biased, which is a
challenge and would require computational methods that can exploit the
important features without relying on their completeness. Finally, it would be
important to develop techniques to handle the case where covariates are
missing for certain regions or time periods in a given country.

## Methods

### Data

To validate our proposed methodology, we use covariates that are related to
population, and that can be derived from remote sensing imagery, open geo-
spatial data (e.g., OSM) or governmental sources. For the present study we
rely on data that have been preprocessed by the WorldPop [23] project. We
collect the covariates listed in Table 4, with a resolution of 100 $\times$
100 m. A selection of variables is visualized in the left part of Figure 8.
Moreover, we obtain census data for the countries of Tanzania (ward-level,
$n=3654$) Zambia (ward-level, $n=1421$) and Mozambique (postos
administrativos-level, $n=413$) to evaluate the performance of our proposed
Pomelo method.

Table 4: Summary of the used covariates. Type |  |  Description
---|---|---
Buildings | 1 |  Building counts from Google Open Buildings [31]
2 |  Building counts from Gridded Maps of Building Patterns [32]
3 |  Mean building areas from Google Open Buildings [31]
4 |  Mean building areas from Gridded Maps of Building Patterns [32]
Accessibility | 5 |  Travel time to city with more than 50k inhabitants [23]
Nightlight | 6 |  Cloud-free DMSP nightlight composite [43]
7 |  Cloud-free VIIRS nightlight composite [43]
Settlement | 8 |  Distance to built-up area from Global Human Settlement layer (GHSL) [44]
9 |  Distance to built-up area from Built-Settlement Growth Model (BSGM) [23]
Topography | 10 |  Elevation model from the Shuttle Radar Topography Mission (SRTM) [45]
11 |  Slope from the Shuttle Radar Topography Mission (SRTM) [45]
Water | 12 |  Distance to Waterbody from ESA CCI Water [46]
13 |  Distance to Waterway from OpenStreetMap [27]
14 |  Distance to Coastline [23]
Land cover | 15 |  Distance to Nature Reserves from World Database on Protected Areas (WDPA) [47]
Roads | 16 |  Distance to Road from OpenStreetMap [27]
17 |  Distance to Road Intersection derived from OpenStreetMap [27]

We also aggregate the census counts to the second-finest administrative level
for the countries of Tanzania (district-level), Zambia (constituency-level),
and Mozambique (district-level), using administrative boundaries from the
Humanitarian Data Exchange [48], to simulate scenarios where the census counts
are coarser. We will later refer to these aggregated data as coarse census
counts. Finally, to assess the generalization performance of Pomelo across
countries, we collect the same source data for four further countries, Uganda
(UGA), Rwanda (RWA), Nigeria (NGA), and Democratic Republic of Congo (COD), to
be used as additional training data. The dataset is summarized in Table 5.

Table 5: Overview of census data used and data extent in # of pixels. Each pixel covers a region of 100 $\times$ 100 m. The Category Others refers to the additional countries of Uganda, Rwanda, Nigeria and the Democratic Republic of Congo. | # census regions | # pixels | Avg. fine census region
---|---|---|---
| fine | coarse | total | built-up | area | population
Tanzania | 3654 | 170 | 111M | 7.5% |  0’333 km2 |  015’449
Zambia | 1421 | 150 |  090M | 3.4% |  0’633 km2 |  013’035
Mozambique |  0413 | 156 |  096M | 5.6% |  2’324 km2 |  075’095
Others | 2779 | 295 | 412M | 5.5% |  9’975 km2 | 139’171

### Proposed method: Population Mapping by Estimating Local Occupancy Rates

Our goal is to produce population maps with a finer resolution than the
underlying census counts. In machine learning terms we are thus faced with an
instance of _weakly supervised_ learning [49]: we do not have access to ground
truth values for the individual grid cells, rather we only have one target
value (the census count for a region) as supervision signal for a whole set of
outputs (all grid cells within the region).

This cumulative sum per region is the supervision signal used by our model.
More specifically, we employ a neural network, c.f. Figure 8, that maps the
covariates to population density estimates and that can be trained with lower-
resolution ground truth: rather than per-pixel supervision, which is not
available, the training procedure compares _aggregated_ estimates per
administrative region to the available region-level census data.

Figure 8: Overview of the proposed method: We use a neural network to create a
mapping between covariates (blue circles) and population density. The direct
output of the neural network are building occupancy rates (green circle) which
we multiply with the building count map to obtain population counts. We
receive the supervision signal for training by aggregating the population maps
to the respective census regions through summation, and comparing our
estimates with the available census data.

We denote the covariates at grid location $l_{i}$ as $X(l_{i})$, the
population estimate at the same location as $\hat{p}\left(l_{i}\right)$, and
$c_{j}$ as the true census data for region $A_{j}$. We observe that the
correlation between the number of buildings at a given location and its
population turns out to be so strong that even dasymetric disaggregation based
_only_ on building counts produces reasonable (albeit less precise) population
density maps, as shown in the Results section. We exploit this strong
correlation for the design of the neural network’s architecture as follows: To
obtain even more accurate estimates, we predict a spatially varying _building
occupancy rate_ $f_{\theta}(X(l_{i}))$, rather than the population count
directly. The population estimate $\hat{p}\left(l_{i}\right)$ is then computed
from it as

$\hat{p}\left(l_{i}\right)=\mathit{f_{\theta}}\big{(}X(l_{i})\big{)}\cdot
b\left(l_{i}\right),$ (1)

where $b\left(l_{i}\right)$ is the number of buildings at location $l_{i}$
(obtained from existing building products [31, 32]) and $\theta$ are the
parameters of the neural network. We use a neural network composed of a
sequence of convolutional layers with kernel sizes $1\times 1$, i.e., its
architecture is equivalent to a per-pixel multi-layer perceptron (MLP).
Empirically, we observe that adding more spatial context (e.g., via
convolutional layers with kernel sizes >1) does not significantly improve the
model’s prediction performance (see results in Table 7). At first this may
seem surprising, but at the low GSD of our maps most covariates are smooth.
Consequently, the disadvantages of a large receptive field outweigh the
potential benefits: feature values of neighboring pixels tend to be very
similar and hence do not add much information. Moreover, the spatial averaging
effect still causes some loss of high-frequency detail (a similar effect has
been observed by other authors, e.g., de Lutio et al. [22]). Also, the
increased number of parameters is more prone to overfitting. A practical
benefit of using 1$\times$1 kernels is that, in the absence of spatial
interactions, there is no spatial diffusion of the input information. All
computations can be restricted to the small fraction of grid cells with at
least one building. This makes the network memory-efficient, an important
feature when processing large countries with coarse census regions, as GPU
memory is often the computational bottleneck for modern neural networks.

The fundamental insight behind Pomelo (and any other top-down disaggregation
scheme) is that, given enough data, the weak supervision is sufficient: the
sum over many grid cells is only a weak constraint per census region, but it
accumulates over different regions because similar covariate features at
different locations must yield similar population densities; whereas the input
covariates have (at least) the target resolution and inject the missing high-
frequency information. Formally, measuring the loss $\mathcal{L}$ over all
grid cells $l_{i}$ in each region $A_{j}$ gives rise to the following
optimization problem:

$\begin{split}&\hat{\theta}=\underset{\theta}{\operatorname{argmin}}\sum_{j}\mathcal{L}\left(\hat{c}_{j},c_{j}\right)\;,\\\
&\text{with }\hat{c}_{j}=\sum_{i\in
A_{j}}\hat{p}\left(l_{i}\right).\end{split}$ (2)

Empirically, the loss function that achieves the best results is the $L_{1}$
distance between the log-transformed predictions and targets:

$\mathcal{L}\left(\hat{c}_{j},c_{j}\right)=|\log\hat{c}_{j}-\log c_{j}\,|\;.$
(3)

If census data are available at inference time, we can leverage the
associated, additional constraints to perform disaggregation (dasymetric
mapping) as a form of post-processing. In other words, we treat the final
outputs of our model $\hat{p}\left(l_{i}\right)$ not as absolute quantities
but as relative proportions. By linearly rescaling the predictions within a
census region $j$ such that they add up to 1, we obtain weights that indicate
what fraction of the total census count $c_{j}$ should be assigned to each
location $l_{i}$. By multiplying those weights with the $c_{j}$, we then
obtain the disaggregated population estimates $\hat{p}_{adj}$, adjusted to
exactly match the coarse census counts:

$\hat{p}_{adj}\left(l_{i}\right)=\frac{\hat{p}\left(l_{i}\right)}{\sum_{k\in
A_{j}}\hat{p}\left(l_{k}\right)}\cdot c_{j}.$ (4)

As long as the census data are correct, that adjustment can be expected to
improve the population density map compared to the raw, absolute estimates. If
no census counts are available at all, we just keep the raw estimates
$\hat{p}\left(l_{i}\right)$.

### Model Setup

Our neural network is composed of four convolutional layers with kernel sizes
$1\times 1$ or, equivalently, of a four-layer MLP that is applied per pixel.
The first three layers each have 128 filters, and the last layer has one
filter to output a scalar density value per location. We use dropout [50]
(with probability 0.4), and rectified linear unit (ReLU) layers after each
convolutional layer. We apply a softplus function to the output of the last
convolutional layer, constraining the occupancy rates to positive numbers. The
model is trained using the Adam optimizer [51] with a base learning rate of
$0.0001$. The weight decay parameter for regularization is optimized via grid-
search on the validation set. To mitigate the low number of training samples,
we propose a data augmentation strategy specific to our task, namely we create
artificial "pseudo-regions" from two real administrative regions, by merging
their pixels and summing their population counts.

The importance of the main components of our model design is ablated in Table
6. It shows quantitative results obtained with different model setups for the
scenario with coarse supervision in Tanzania. We start from a baseline that
has the same model architecture, but

  1. 1.

employs the standard $L_{1}$ distance between predictions and target values as
loss function (no $\log L_{1}$ loss);

  2. 2.

directly predicts population counts per grid cell, and

  3. 3.

is trained only with the actual administrative regions (no region-based
augmentation).

We then gradually add the components used in our proposed setting. First, we
find that computing the loss with log-transformed outputs significantly
improves over the standard $L_{1}$ loss. Second, predicting the building
occupancy rate instead of the population and transforming it to a population
count in postprocessing also clearly reduces the error, especially the $R^{2}$
metric. Third, data augmentation by synthetically merging census regions also
has a positive effect, confirming the beneficial effect of problem-specific
augmentation when training data is scarce. Finally, combining all three
measures, as in our proposed model, yields greatly improved predictive skill,
reaching $\approx$25% lower mean absolute error, respectively >14 percent
points higher $R^{2}$.

Table 6: Effect of modeling choices (Tanzania data, coarse supervision). Method | LogL1 | Occ. | Aug. |  $R^{2}$ $\uparrow$ [%] |  MAE $\downarrow$ |  MAPE $\downarrow$ [%]
---|---|---|---|---|---|---
Baseline |  |  |  |  71.4 $\pm$ 0.3 |  4’130 $\pm$ 180 |  27.5 $\pm$ 0.9
\+ Log$L_{1}$ loss | ✓ |  |  |  73.2 $\pm$ 0.3 |  3’580 $\pm$ 130 |  23.0 $\pm$ 0.2
\+ Occupancy rate |  | ✓ |  |  81.4 $\pm$ 1.2 |  3’820 $\pm$ 180 |  26.8 $\pm$ 1.2
\+ Augmentation |  |  | ✓ |  73.0 $\pm$ 1.2 |  3’820 $\pm$ 130 |  25.0 $\pm$ 0.4
Pomelo | ✓ | ✓ | ✓ |  85.8 $\pm$ 0.7 |  3’100 $\pm$ 130 |  21.7 $\pm$ 0.3

Table 7 show the effect of larger kernel widths on Pomelo’s performance for
Tanzania, in both the coarse and fine supervision settings. It can be seen
that larger kernel sizes are consistently detrimental, it appears that the
they increase model complexity beyond the level that can be learned with the
available supervision. On the one hand the bigger kernels inflate the number
of weights that must be trained ($9\times$ more, respectively $25\times$
more). These additional degrees of freedom increase the risk of overfitting
and make learning harder – particularly in our weakly supervised setting,
which becomes apparent with coarse supervision from only 170 regions. While on
the other hand the covariate maps in general exhibit a high degree of spatial
smoothness and do thus do not contain high-frequency context that larger
kernels could capture.

Table 7: Performance with varying kernel sizes on the dataset of Tanzania. Training Setting | Kernel Size | #params |  $R^{2}\uparrow$ [%] |  MAE $\downarrow$ |  MAPE $\downarrow$ [%]
---|---|---|---|---|---
coarse |  $1\times 1$ (ours) |  034k |  85.7 $\pm$ 0.9 |  3’100 $\pm$ 040 |  21.6 $\pm$ 0.2
$3\times 3$ | 313k |  83.5 $\pm$ 1.3 |  3’700 $\pm$ 110 |  26.6 $\pm$ 0.7
$5\times 5$ | 870k |  80.9 $\pm$ 1.0 |  3’900 $\pm$ 100 |  28.4 $\pm$ 0.8
fine |  $1\times 1$ (ours) |  034k |  87.6 $\pm$ 0.2 |  2’890 $\pm$ 020 |  20.4 $\pm$ 0.3
$3\times 3$ | 313k |  86.0 $\pm$ 2.1 |  3’010 $\pm$ 100 |  20.9 $\pm$ 0.5
$5\times 5$ | 870k |  86.8 $\pm$ 1.0 |  3’030 $\pm$ 050 |  21.5 $\pm$ 0.4

We also empirically study the choice of input covariates. To that end, we only
keep the five input layers with the highest individual importance according to
the permutation method [38] and discard the remaining features. The results
with this reduced input set for Tanzania (fine supervision setting) are
reported in Table 8. Empirically, explicit variable selection noticeably harms
performance. Apparently features with lower importance scores still carry
valuable information that complements the most important features. In this
context, note that in the proposed scheme the number of input layers is not
critical: with the Pomelo architecture, an additional covariate only
introduces 128 learnable parameters. The computational savings achievable with
feature selection are negligible.

Table 8: Performance on Tanzania data with reduced covariate set (fine supervision setting). Method |  $R^{2}\uparrow$ [%] |  MAE $\downarrow$ |  MAPE $\downarrow$ [%]
---|---|---|---
all features (ours) |  87.6 $\pm$ 0.2 |  2’890 $\pm$ 20 |  20.4 $\pm$ 0.3
top 5 features |  85.7 $\pm$ 0.4 |  2’980 $\pm$ 10 |  20.7 $\pm$ 0.0

### Methods used for numerical comparisons

We compare Pomelo with four other methods:

Building count disaggregation: Dasymetric disaggregation using only the
available building counts per pixel as weights.

Random Forest (RF) at region level: Our own re-implementation of the random
forest model used by WorldPop [6]. The model is fed the same covariates
(features) as the Pomelo network. The training units for that scheme are not
grid cells but administrative regions, with features aggregated over all
pixels within a region. The map visualizations were created in QGIS 3.14 [26].

Markov Random Field (MRF): Here, disaggregation is based directly on the
assumption that locations with similar features have similar population
densities. To enforce such prior knowledge explicitly, we resort to an MRF
[52] that compares pixel-wise predictions to each other and to their
encompassing region. This involves finding the combination of population
values that minimizes the following (negative log-likelihood) energy function:

$\displaystyle\mathbf{\hat{p}}_{MRF}=\operatorname*{arg\,min}_{\mathbf{p}\in\mathbb{R}^{N}}\sum_{i}^{N}\sum_{k\in
Q_{i}}|p_{i}-p_{k}|+\lambda\sum_{j}|c_{j}-\sum_{k\in A_{j}}p_{k}|\;,$ (5)

where $N$ is the number of grid cells and $Q_{i}$ the set of nearest neighbors
to cell $l_{i}$ in the (normalized) feature space. The first term encourages
locations with similar features (i.e., the $k$ nearest neighbors to pixel $i$
in feature space) to have similar population counts; the second term pushes
the aggregate count over an administrative region towards the region’s total
census. Parameter $\lambda$ determines the balance between the two
constraints; for all the experiments we set $\lambda=1$. We create the
(approximate) nearest-neighbor graph with the fast ANN method [53], initialize
the population counts with the building count disaggregation described above,
and find an optimal configuration by minimizing the energy function with the
Iterated Conditional Modes algorithm [52], using update steps of $\pm 1\%$.
Empirically, we found that using too many covariates as features harms the
performance of the MRF model, likely due to the higher dimension of the
associated feature space. The best performance, reported in the Results
Section, are obtained with only three covariates: building count, average
building size, and night lights.

Convolutional Neural Network (CNN): It has been proposed to estimate
population counts with a CNN with spatial context [19], trained with the same
form of aggregated supervision as Pomelo. We instead advocate for a shared
per-pixel architecture. Hence, we further include a baseline that has the same
network architecture as Pomelo, with two exceptions: the kernel size is set to
3$\times$3 for all convolutional layers, and the network is trained to
directly output population as suggested by [19] (rather than densities, as in
our system). We point out that, while our CNN baseline is inspired by that
study, the two are not directly comparable: Jacobs et al. target local
population mapping in urban areas and base their estimates on high-resolution
Planet imagery with 3 m GSD, and on city block-level population and housing
counts from the US census. Consequently, they can also afford to train a much
larger (U-net [54]) model.

### Evaluation metrics

For all three evaluation strategies, shown in the Results section, the
estimated maps are evaluated by aggregating the per-pixel counts back to a
list of $n_{c}$ population numbers $\mathbf{\hat{c}}$ at finest available
census level and comparing them to the actual census counts $\mathbf{c}$ in
terms of $R^{2}$, mean absolute error (MAE), and mean absolute percentage
error (MAPE):

$\begin{split}&R^{2}(\mathbf{\hat{c}},\mathbf{c})=1-\frac{\sum_{j=1}^{n_{c}}(c_{j}-\hat{c}_{j})^{2}}{\sum_{j=1}^{n_{c}}(c_{j}-\bar{c})^{2}}\\\
&\text{MAE}(\mathbf{\hat{c}},\mathbf{c})=\frac{1}{n_{c}}\sum_{j=1}^{n_{c}}|c_{j}-\hat{c}_{j}|\\\
&\text{MAPE}(\mathbf{\hat{c}},\mathbf{c})=100\cdot\frac{1}{n_{c}}\sum_{j=1}^{n_{c}}\bigg{|}\frac{c_{j}-\hat{c}_{j}}{c_{j}}\bigg{|}\end{split}$
(6)

## Data Availability

All data used in this study is publicly available. The population maps can be
accessed via https://doi.org/10.6084/m9.figshare.21444282.v1, while the code
is available in the following repository:
https://github.com/jvargasmu/population_estimation.

## Acknowledgments

This research is funded by the Science and Technology for Humanitarian Action
Challenges (HAC) project.

## Author contributions

D.T., K.S., T.W., and F.O. conceived the study, N.M., J.V., R.D., K.S., and
D.T. designed the research, N.M., and J.V. implemented the methods. All the
authors discussed the results and wrote the manuscript.

## Competing interests

The authors declare no competing interests.

## References

  * [1] Naess, P.  Urban planning and sustainable development.  _European Planning Studies_ 9, 503–524 (2001).
  * [2] Lang, S. _et al._ Earth observation tools and services to increase the effectiveness of humanitarian assistance.  _European Journal of Remote Sensing_ 53, 67–85 (2019).
  * [3] Maas, P. _et al._ Facebook disaster maps: Aggregate insights for crisis response & recovery.  _ACM SIGKDD Conference on Knowledge Discovery & Data Mining_ 19, 3173 (2019).
  * [4] Thomson, D. R., Rhoda, D. A., Tatem, A. J. & Castro, M. C.  Gridded population survey sampling: a systematic scoping review of the field and strategic research agenda.  _International Journal of Health Geographics_ 19, 1–16 (2020).
  * [5] Vargas-Muñoz, J. E., Srivastava, S., Tuia, D. & Falcão, A. X.  OpenStreetMap: Challenges and opportunities in machine learning and remote sensing.  _IEEE Geoscience and Remote Sensing Magazine_ 9, 184–199 (2020).
  * [6] 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, e0107042 (2015).
  * [7] Weber, E. M. _et al._ Census-independent population mapping in northern Nigeria.  _Remote Sensing of Environment_ 204, 786–798 (2018).
  * [8] Hillson, R. _et al._ Methods for determining the uncertainty of population estimates derived from satellite imagery and limited survey data: a case study of Bo City, Sierra Leone.  _PloS one_ 9, e112241 (2014).
  * [9] Leasure, D. R., Jochem, W. C., Weber, E. M., Seaman, V. & Tatem, A. J.  National population mapping from sparse survey data: A hierarchical Bayesian modeling framework to account for uncertainty.  _Proceedings of the National Academy of Sciences_ 117, 24173–24179 (2020).
  * [10] Huang, X., Wang, C., Li, Z. & Ning, H.  A 100 m population grid in the CONUS by disaggregating census data with open-source Microsoft building footprints.  _Big Earth Data_ 5, 112–133 (2021).
  * [11] Freire, S., MacManus, K., Pesaresi, M., Doxsey-Whitfield, E. & Mills, J.  Development of new open and free multi-temporal global population grids at 250 m resolution.  _Association of Geographic Information Laboratories in Europe (AGILE)_ (2016).
  * [12] Sapena, M. _et al._ Empiric recommendations for population disaggregation under different data scenarios.  _Plos one_ 17, e0274504 (2022).
  * [13] Tiecke, T. G. _et al._ Mapping the world population one building at a time.  _arXiv preprint arXiv:1712.05839_ (2017).
  * [14] Sridharan, H. & Qiu, F.  A spatially disaggregated areal interpolation model using light detection and ranging-derived building volumes.  _Geographical Analysis_ 45, 238–258 (2013).
  * [15] Reibel, M. & Bufalino, M. E.  Street-weighted interpolation techniques for demographic count estimation in incompatible zone systems.  _Environment and Planning A_ 37, 127–139 (2005).
  * [16] Dobson, J. E., Bright, E. A., Coleman, P. R., Durfee, R. C. & Worley, B. A.  LandScan: a global population database for estimating populations at risk.  _Photogrammetric Engineering and Remote Sensing_ 66, 849–857 (2000).
  * [17] LeCun, Y., Bengio, Y. & Hinton, G.  Deep learning.  _nature_ 521, 436–444 (2015).
  * [18] Gervasoni, L., Fenet, S., Perrier, R. & Sturm, P.  Convolutional neural networks for disaggregated population mapping using open data.  _IEEE Conference on Data Science and Advanced Analytics (DSAA)_ 594–603 (2018).
  * [19] Jacobs, N., Kraft, A., Rafique, M. U. & Sharma, R. D.  A weakly supervised approach for estimating spatial density functions from high-resolution satellite imagery.  _ACM SIGSPATIAL Conference on Advances in Geographic Information Systems_ 33–42 (2018).
  * [20] Diebel, J. & Thrun, S.  An application of Markov random fields to range sensing.  _Advances in Neural Information Processing Systems_ (2006).
  * [21] Hui, T.-W., Loy, C. C. & Tang, X.  Depth map super-resolution by deep multi-scale guidance.  _European Conference on Computer Vision (ECCV)_ 353–369 (2016).
  * [22] de Lutio, R., D’Aronco, S., Wegner, J. D. & Schindler, K.  Guided super-resolution as pixel-to-pixel transformation.  _IEEE/CVF International Conference on Computer Vision (ICCV)_ (2019).
  * [23] Worldpop.  www.worldpop.org (2022).  Accessed: 2022-06-20.
  * [24] Mennis, J.  Generating surface models of population using dasymetric mapping.  _The Professional Geographer_ 55, 31–42 (2003).
  * [25] Grippa, T. _et al._ Improving urban population distribution models with very-high resolution satellite information.  _Data_ 4, 13 (2019).
  * [26] QGIS Development Team.  _QGIS Geographic Information System, Version 3.14_.  Open Source Geospatial Foundation (2009).
  * [27] OpenStreetMap contributors.  Data retrieved from https://planet.osm.org .  https://www.openstreetmap.org (2022).
  * [28] European Space Agency.  _Sentinel-2_.
  * [29] Meta.  _Facebook data for good high resolution population density maps_ (2022).
  * [30] Gorelick, N. _et al._ Google earth engine: Planetary-scale geospatial analysis for everyone.  _Remote sensing of Environment_ 202, 18–27 (2017).
  * [31] Sirko, W. _et al._ Continental-scale building detection from high resolution satellite imagery.  _arXiv preprint arXiv:2107.12283_ (2021).
  * [32] Dooley, C., Leasure, D. & Tatem, A.  Gridded maps of building patterns throughout sub-Saharan Africa, version 1.1.  Tech. Rep., University of Southampton, Southampton, UK (2020).
  * [33] Vargas-Muñoz, J. E., Lobry, S., Falcão, A. X. & Tuia, D.  Correcting rural building annotations in OpenStreetMap using convolutional neural networks.  _ISPRS journal of photogrammetry and remote sensing_ 147, 283–293 (2019).
  * [34] Nieves, J. J. _et al._ Annually modelling built-settlements between remotely-sensed observations using relative changes in subnational populations and lights at night.  _Computers, Environment and Urban Systems_ 80, 101444 (2020).
  * [35] Nieves, J. J. _et al._ Predicting near-future built-settlement expansion using relative changes in small area populations.  _Remote Sensing_ 12, 1545 (2020).
  * [36] Neal, I., Seth, S., Watmough, G. & Diallo, M. S.  Census-independent population estimation using representation learning.  _Scientific Reports_ 12, 1–12 (2022).
  * [37] Barrington-Leigh, C. & Millard-Ball, A.  The world’s user-generated road map is more than 80% complete.  _PloS one_ 12, e0180698 (2017).
  * [38] Breiman, L.  Random forests.  _Machine Learning_ 45, 5–32 (2001).
  * [39] Esch, T. _et al._ World settlement footprint 3d-a first three-dimensional survey of the global building stock.  _Remote Sensing of Environment_ 270, 112877 (2022).
  * [40] Meta.  _Facebook Marketing API_ (2022).
  * [41] Fatehkia, M., Coles, B., Ofli, F. & Weber, I.  The relative value of Facebook advertising data for poverty mapping.  _AAAI Conference on Web and Social Media_ 934–938 (2020).
  * [42] Leasure, D. R. _et al._ Ukraine crisis: Monitoring population displacement through social media activity.  _SocArXiv, doi:10.31235/osf.io/6j9wq_ (2022).
  * [43] Elvidge, C. D., Baugh, K. E., Zhizhin, M. & Hsu, F.-C.  Why viirs data are superior to dmsp for mapping nighttime lights.  _Proceedings of the Asia-Pacific Advanced Network_ 35, 62 (2013).
  * [44] Corbane, C., Florczyk, A., Pesaresi, M., Politis, P. & Syrris, V.  An overview and comparison of smooth labeling methods for land-cover classification.  _IEEE transactions on geoscience and remote sensing_ DOI: 10.2905/jrc-ghsl-10007 (2018).
  * [45] Farr, T. G. _et al._ The shuttle radar topography mission.  _Reviews of geophysics_ 45 (2007).
  * [46] Lamarche, C. _et al._ Compilation and validation of SAR and optical data products for a complete and global map of inland/ocean water tailored to the climate modeling community.  _Remote Sensing_ 9, 36 (2017).
  * [47] UNEP-WCMC & IUCN.  _Protected Planet: The World Database on Protected Areas (WDPA)_ (2021).
  * [48] Humanitarian data exchange.  https://data.humdata.org/ (2022).  Accessed: 2022-06-20.
  * [49] Zhou, Z.-H.  A brief introduction to weakly supervised learning.  _National Science Review_ 5, 44–53 (2018).
  * [50] Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I. & Salakhutdinov, R.  Dropout: a simple way to prevent neural networks from overfitting.  _The journal of machine learning research_ 15, 1929–1958 (2014).
  * [51] Kingma, D. P. & Ba, J.  Adam: A method for stochastic optimization.  _International Conference on Learning Representations (ICLR)_ (2014).
  * [52] Besag, J.  On the statistical analysis of dirty pictures.  _Journal of the Royal Statistical Society_ 68, 259–302 (1986).
  * [53] Johnson, J., Douze, M. & Jégou, H.  Billion-scale similarity search with GPUs.  _IEEE Transactions on Big Data_ 7, 535–547 (2019).
  * [54] Ronneberger, O., Fischer, P. & Brox, T.  U-net: Convolutional networks for biomedical image segmentation.  _International Conference on Medical Image Computing and Computer-assisted Intervention_ 234–241 (2015).


# Moments of Subsets of General Equiangular Tight Frames

Marina Haikin
Amazon, Tel Aviv, Israel §
<EMAIL_ADDRESS>
Ram Zamir
Tel Aviv University, Israel
<EMAIL_ADDRESS>
(April 2021)

###### Abstract

This note outlines the steps for proving that the moments of a randomly-
selected subset of a general ETF (complex, with aspect ratio $0<\gamma<1$)
converge to the corresponding MANOVA moments. We bring here an extension for
the proof of the ’Kesten-Mckay’ moments (real ETF, $\gamma=1/2$) [1]. In
particular, we establish a recursive computation of the $r$th moment, for
$r=1,2,\ldots$, and verify, using a symbolic program, that the recursion
output coincides with the MANOVA moments.

§§footnotetext: Now at Amazon. Previously at Tel Aviv University where this
work was performed.

##  I Recursive computation of moments

Let $F$ be a unit-norm equiangular tight frame comprised of $n$ vectors in
${\cal C}^{m}$, $\gamma=\frac{m}{n}<1$ is the aspect ratio of the frame. We
define $s\stackrel{{\scriptstyle\Delta}}{{=}}(1-x)/2\sqrt{x}$, where
$x=\gamma^{-1}-1$ and denote a generalized conference matrix as

$S=\sqrt{\frac{n-1}{{x}}}\cdot\left(F^{\prime}F-\frac{1}{2\gamma}\cdot
I_{n}\right)$ (1)

which satisfies:

  * (i)

$S_{i,i}$ = constant =
$\sqrt{n-1}(1-x)/2\sqrt{x}\stackrel{{\scriptstyle\Delta}}{{=}}\sqrt{n-1}\cdot
s$,

  * (ii)

$S$ is conjugate-symmetric and $|S_{i,j}|=1$ for every $i,j\in[n]$ with $i\neq
j$, and

  * (iii)

$S^{\prime}S\propto I_{n}$, i.e., the off-diagonal elements of $S^{2}$ are
zero, while the diagonal elements are equal to $(n-1)(x+1)^{2}/4x$.

Note that $x/(n-1)$ is the Welch bound, i.e., the squared absolute cross
correlation between the frame vectors. For $\gamma=0.5$, we have $x=1,s=0$, so
$S=\sqrt{n-1}\cdot\left(F^{\prime}F-I_{n}\right)$, and $S^{\prime}S=(n-1)I$,
in agreement with the properties of the conference matrix in [1].

We extend the proof in [1] to derive the expression for moments of the random
submatrix $X_{S}=PSP$ defined as

$m^{S}_{k}=\frac{1}{n^{k/2+1}}\cdot E\left\\{{\rm
trace}(X_{S}^{k})\right\\}=\sum_{t=1}^{k}\left(\sum_{\pi\in\prod(k,t)}V_{n}(\pi)\right)\cdot
p^{t}$ (2)

where $P$ is a diagonal matrix such that the diagonal elements are independent
random variables Bernoulli($p$). For every partition $\pi$ of $[k]$ into $t$
blocks, $V_{n}(\pi)$ is sum of products of $k$ elements from the generalized
conference matrix $S$, as a function of the partition. The sum is over all
possible $t$ distinct values from $[n]$ according to the blocks of $\pi$
(value per block),

$V_{n}(\pi)\stackrel{{\scriptstyle\Delta}}{{=}}\frac{1}{n^{k/2+1}}\sum_{a\in
L_{n}(\pi)}S_{a(1),a(2)}S_{a(2),a(3)}\cdots S_{a(k),a(1)},$ (3)
$L_{n}(\pi)\stackrel{{\scriptstyle\Delta}}{{=}}\\{a:[k]\rightarrow[n]:\\{a^{-1}(a(i)):i\in[k]\\}=\pi\\}.$
(4)

The main difference from the KM case will be the computation of $V_{n}(\pi)$.
Lemma 16 in [1] will be replaced by

######  Lemma 1.

For every non-crossing partition $\pi\in\prod(k,t)$, the following holds: if
the edges of $G_{\pi}$ partition into m simple cycles of sizes
$l_{1},\dots,l_{m}$, then as $n\rightarrow\infty$,

$V_{n}(\pi)\rightarrow A_{l_{1}}\cdots A_{l_{m}}$ (5)

where

$A_{l+1}=-\sum_{i=1}^{l}A_{i}A_{l+1-i},\,\,\,A_{1}=s,A_{2}=1.$ (6)

Note that for $\gamma=0.5$, $s=0$ and we get that for odd cycle, i.e., odd
$l_{i}$, $A_{l_{i}}=0$, and for even $l_{i}$, $A_{l_{i}}=(-1)^{l_{i}/2-1}\cdot
C_{l_{i}/2-1}$. Thus if $G_{\pi}$ contains any odd cycle, then
$V_{n}(\pi)\rightarrow 0$ otherwise $V_{n}(\pi)$ converges to a product of
Catalan numbers, in agreement with Lemma 16 in [1].

For general $\gamma$ ($s\neq 0$), we don’t have yet an a analogy to Lemma 8,
but $\sum_{\pi\in\prod(k,t)}V_{n}(\pi)$ can be computed by counting all
$N(k,t)$ (Narayana number) non-crossing partitions of $[k]$ into $t$ blocks
and for each induced cactus to compute $V_{n}(\pi)$ according to Lemma 1.

The desired moments of a (Bernoulli-$p$ selected) subset of ETF $F$ are then
derived by using (1) and the binom formula:

$\displaystyle m_{k}$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}E\left\\{\frac{1}{n}\cdot{\rm
trace}((PF^{\prime}FP)^{k})\right\\}$ $\displaystyle=\sum_{j=0}^{k}{k\choose
j}\cdot
x^{j/2}\cdot\left(\frac{x+1}{2}\right)^{k-j}p^{j==0}\cdot\frac{1}{n^{j/2+1}}\cdot
E\left\\{{\rm trace}(X_{S}^{j})\right\\}$
$\displaystyle=\left(\frac{x+1}{2}\right)^{k}\cdot p+\sum_{j=1}^{k}{k\choose
j}\cdot x^{j/2}\cdot\left(\frac{x+1}{2}\right)^{k-j}\cdot m^{S}_{j}$ (7)

where $m^{S}_{j}$ is defined in (2), and $j==0$ denotes $1$ if $j=0$ or $0$
otherwise.

Before we turn to the proof, we apply the recursion for several first moments:

$\displaystyle m^{S}_{0}=1\,\,\,\,\,\,$ $\displaystyle m_{0}=1$ $\displaystyle
A_{1}=s\,\,\,\,\,\,$ $\displaystyle m^{S}_{1}=sp\,\,\,\,\,\,$ $\displaystyle
m_{1}=p$ $\displaystyle A_{2}=1\,\,\,\,\,\,$ $\displaystyle
m^{S}_{2}=s^{2}p+p^{2}\,\,\,\,\,\,$ $\displaystyle m_{2}=p+p^{2}x$
$\displaystyle A_{3}=-2s\,\,\,\,\,\,$ $\displaystyle
m^{S}_{3}=s^{3}p+3sp^{2}+(-2s)p^{3}\,\,\,\,\,\,$ $\displaystyle\vdots$
$\displaystyle A_{4}=4s^{2}-1\,\,\,\,\,\,$ $\displaystyle
m^{S}_{4}=s^{4}p+6s^{2}p^{2}+(-8s^{2}+2)p^{3}+(4s^{2}-1)p^{4}$ $\displaystyle
A_{5}=-8s^{3}+6s\,\,\,\,\,\,$ $\displaystyle\vdots$ $\displaystyle
A_{6}=16s^{4}-24s^{2}+2$

###### Proof.

We follow the notations in [1]
$\Delta(a(1),a(2),\dots,a(k))=S_{a(1),a(2)}S_{a(2),a(3)}\cdots S_{a(k),a(1)}$.
If $a_{j}=a_{j+1}$ for any $j\in[k-1]$ or $a_{k}=a_{1}$,

$\displaystyle V_{n}(\pi)$ $\displaystyle=\frac{1}{n^{k/2+1}}\sum_{a\in
L_{n}(\pi)}\Delta(a(1),\dots a(j-1),a(j+1),a(j+1),\dots,a(k))$ (8)
$\displaystyle=\frac{\sqrt{n-1}s}{n^{k/2+1}}\sum_{a\in
L_{n}(\pi)}\Delta(a(1),\dots,a(j-1),a(j+1),\dots,a(k))$ (9)
$\displaystyle=\frac{s}{n^{(k-1)/2+1}}\sum_{a\in
L_{n}(\pi\setminus\\{j\\})}\Delta(a(1),\dots,a(j-1),a(j+1),\dots,a(k))+o(1)$
(10)

were the first equality follows from the fact that the diagonal entries of $S$
are $s\sqrt{n-1}$, thus $S_{a(j+1)a(j+1)}=s\sqrt{n-1}$. The restriction of
$\pi\setminus\\{j\\}$ to $[k]\setminus{j}$ results in a partition
$\pi^{\prime}$ of $[k]\setminus{j}$ into same $t$ blocks besides removing $j$
from the block $\pi(j)$. Thus the above expression of $V_{n}(\pi)$ implies

$V_{n}(\pi)\rightarrow s\cdot V_{n}(\pi^{\prime}),$ (11)

where throughout the proof, $\rightarrow$ means as $n\rightarrow\infty$. This
means that two consecutive indices which belong to same block and contribute a
loop to the graph $G_{\pi}$, contribute a constant factor $s$ to $V_{n}(\pi)$.
This claim, as well as the derivation in (8)-(10), extends for several loops.
Given partition $\pi$ with $r$ loops in total (in one or more blocks), we can
consider a squeezed partition $\pi^{\prime}$ with same number of blocks and
$k^{\prime}=k-r$. Each block $B^{\prime}_{i}$ for $i=1,\dots,t$, is a squeezed
block $B_{i}$, with only one representative index of sequence of several
consecutive equal indices.

$V_{n}(\pi)\rightarrow s^{r}\cdot V_{n}(\pi^{\prime})$ (12)

For example a partition 122234422333341 translated to 1234234 with a factor of
$s^{8}$.

Now we focus on estimation of $V_{n}(\pi^{\prime})$. $G_{\pi^{\prime}}$ is a
loop-free graph and thus Lemma 6 in [1] holds for the same reasons. Lemma
12(iv),12(v) also holds due to similar the properties of $S$, and thus Lemma
14 - which relies on all of the above - holds. This proves the claim of
vanishing crossing partitions, i.e., let $\pi^{\prime}\in\prod(k^{\prime},t)$
be a crossing partition, then $V_{n}(\pi^{\prime})\rightarrow 0$.

We follow the proof by induction of Lemma 16 in [1]. Lemma 12(i) holds again
by $S$ properties and thus for $t=2$,
$V_{n}(\pi^{\prime})=V_{n}(\\{1\\},\\{2\\})=1$. As for $G_{\pi}$, with assumed
$r$ cycles of sizes 1, proving the lemma for $V_{n}(\pi^{\prime})$, will
immediately imply the proof for $V_{n}(\pi)$, as $A_{1}^{r}=s^{r}$. Assuming
that the (loops-removed) partition $\pi^{\prime}$ is non-crossing, it follows
that the graph $G_{\pi^{\prime}}$ is a cactus, and by Lemma 15 it is
guaranteed to contain a singleton leaf (block in the partition). As we deal
with $\pi^{\prime}$, with $m^{\prime}=m-r$ cycles, we aim to prove that

$V_{n}(\pi^{\prime})\rightarrow A_{l_{1}}\cdots\ A_{l_{m^{\prime}}}.$ (13)

Case I of the proof refers to the case when the singleton block resides in a
cycle of length 2, i.e., if the singleton is $\\{k\\}\in\pi$, $\pi(1)=\pi(k)$.
A similar derivation as in [1] implies

$V_{n}(\pi^{\prime})=V_{n}(\pi^{\prime\prime})+o(1).$ (14)

Our induction hypothesis and (14), imply

$V_{n}(\pi^{\prime})\rightarrow A_{l_{1}}\cdots\ A_{l_{m^{\prime}-1}}.$ (15)

Since $l_{m}=2$ and $A_{2}=1$, this establishes (13).

Case II of the proof refers to the case when $\\{k\\}\in\pi$ (the singleton
block) resides in a cycle of length $l\geq 3$. With similar analysis to [1] we
get:

$V_{n}(\pi^{\prime})=-\sum_{i=2}^{l}V_{n}({\pi^{\prime}}^{i})+o(1).$ (16)

Assuming $l_{m}=l$ and applying the induction hypothesis, we have

$V_{n}(\pi^{\prime})\rightarrow-A_{l_{1}}\cdots\
A_{l_{m^{\prime}-1}}\cdot\sum_{i=2}^{l}A_{l-1}A_{l-i+1}=-A_{l_{l}}\cdots\
A_{l_{m^{\prime}-1}}\cdot\sum_{i=1}^{l-1}A_{l}A_{l-i}.$ (17)

Applying the identity (6), we get $V_{n}(\pi^{\prime})\rightarrow
A_{l_{1}}\cdots\ A_{l_{m^{\prime}}}$, thereby establishing (13). ∎

##  II IMPLEMENTATION OF THE DERIVED ALGORITHM

⬇

import numpy as np

import sympy

from sympy.utilities.iterables import multiset_partitions

import scipy.special

# recursuve calculation of A(l) for l=1..n, arr_A[i]=A(i+1)

def calc_A(n,s):

arr_A = np.zeros(n)

arr_A[0] = s

arr_A[1] = 1

for i in np.arange(n-2)+2:

arr_A[i] = -sum([arr_A[j]*arr_A[i-j-1] for j in np.arange(i)])

return arr_A

# generation of all non crossing partitions of k to t blocks

def calc_partitions_k_t(k,t):

partitions = list(multiset_partitions(np.arange(k),t))

non_crossing = []

non_crossing_len = []

for part in partitions:

cross = 0

n = len(part)

if n == 1:

non_crossing.append(part)

non_crossing_len.append([len(b) for b in part])

else:

for b1 in part:

for b2 in part:

if (b2[0]>b1[0])&(len(np.array(np.where((b1>b2[0])&(b1<b2[-1])))[0])>0):

cross = 1

break

if cross==1:

break

if cross==0:

non_crossing.append(part)

non_crossing_len.append([len(b) for b in part])

return non_crossing

’’’

Find␣length␣of␣cactus␣leafs␣which␣correspond␣to␣a␣non-crossing␣partition

’’’

# genrate lists off all intervals within blocks of partition

# list_pairs: a list of consecutive pairs in blocks [block[i],block[i+1]] for
every block in the partition

# list_pairs_intervals: a list of intervals between consecutive pairs in
blocks block[i+1]-block[i] for every

# block in the partition

def create_pairs_lists(part):

list_pairs = []

list_pairs_intervals = []

for block in part:

b_len = len(block)

if b_len > 1:

for i in range(b_len-1):

list_pairs.append([block[i],block[i+1]])

list_pairs_intervals.append(block[i+1]-block[i])

return [list_pairs,list_pairs_intervals]

# sort the lists by the length of the intervals.

# start computing the the length of cycle which corresponds to an interval,
from shorter to longer -

# for each interval substruct the cycles which correspond to smaller intervals
contained in the interval

def update_lists(list_pairs,list_intervals):

zipped_lists = zip(list_intervals, list_pairs)

sorted_pairs = sorted(zipped_lists)

tuples = zip(*sorted_pairs)

list_intervals, list_pairs = [list(tuple) for tuple in tuples]

for i in range(len(list_intervals)):

intersections_length = 0

for j in range(i):

if (list_pairs[j][0] > list_pairs[i][0]) & (list_pairs[j][1] <
list_pairs[i][1]):

intersections_length = intersections_length + list_intervals[j]

list_intervals[i] = list_intervals[i] - intersections_length

return list_intervals

# generated the updated lists and add the cycle which correspond to an outer
interval [k, 0],

# the length of which is k - sum of all cycles

def partition_to_cactus_sizes(part,k):

[list_pairs,list_intervals] = create_pairs_lists(part)

if len(list_pairs)>0:

list_intervals = update_lists(list_pairs,list_intervals)

cactus_sizes = list_intervals

cactus_sizes.append(k - sum(cactus_sizes))

else:

cactus_sizes = [k]

return cactus_sizes

’’’

calculation␣of␣the␣’centralized’␣moments

’’’

# calculation of asymptotic V_n(pi) for a given partition by eq. 5 and sum
over all possible (k,t) non-crossing

# partitions

def calc_V_pi_sum(k,t):

V_pi_sum = 0

cactuses = [partition_to_cactus_sizes(part,k) for part in
calc_partitions_k_t(k,t)]

for cactus in cactuses:

V_pi_sum = V_pi_sum + np.prod(arr_A[np.array(cactus)-1])

return V_pi_sum

# calculation of m_^s_k by eq. 2

def calc_m_s(k,p):

m_s = 0

for t in np.arange(k)+1:

m_s = m_s + calc_V_pi_sum(k,t)*p**t

return m_s

’’’

binom␣formula␣for␣’de-centralization’␣of␣the␣moments,␣eq.7

’’’

def calc_m(k,x,p):

m = ((x+1)/2)**k*p

for j in np.arange(k)+1:

m = m + scipy.special.binom(k,j)*x**(j/2)*((x+1)/2)**(k-j)*(m_s_vec[j-1])

return m

# run calculation of 10 first moments

gamma = np.sqrt(0.5)

p = 0.6

x = 1/gamma - 1

s = (1-x)/(2*np.sqrt(x))

d = 10

arr_A = calc_A(d,s)

m_s_vec = [calc_m_s(k,p) for k in np.arange(d)+1]

m_vec = [calc_m(k,x,p) for k in np.arange(d)+1]

And the output is:

⬇

\------------------------------------

gamma = 0.7071067811865476

p = 0.6

\------------------------------------

x = 0.4142135623730949

s = 0.45508986056222756

\------------------------------------

calc ETF moments until d = 10

\------------------------------------

A[l] for l = 1..10

[0.45508986056222756, 1.0, -0.9101797211244551, -0.1715728752538092,
1.9765215939999434, -2.2842712474619047,

-2.1862701245144653, 10.126983722080917, -8.324794719173925, -21.347181155833745]

A[l] for l = 1..6 by manual application of eq.6

[0.45508986056222756, 1.0, -0.9101797211244551, -0.17157287525380915,
1.9765215939999432, -2.2842712474619047]

\------------------------------------

example for all partitions (5,3) and the corresponding cactus lengths

[[0, 1, 2], [3], [4]] [1, 1, 3]

[[0, 1, 3], [2], [4]] [1, 2, 2]

[[0, 1], [2, 3], [4]] [1, 1, 3]

[[0, 1, 4], [2], [3]] [1, 3, 1]

[[0, 1], [2, 4], [3]] [1, 2, 2]

[[0, 1], [2], [3, 4]] [1, 1, 3]

[[0, 2, 3], [1], [4]] [1, 2, 2]

[[0, 2, 4], [1], [3]] [2, 2, 1]

[[0, 2], [1], [3, 4]] [1, 2, 2]

[[0, 3], [1, 2], [4]] [1, 2, 2]

[[0], [1, 2, 3], [4]] [1, 1, 3]

[[0, 4], [1, 2], [3]] [1, 3, 1]

[[0], [1, 2, 4], [3]] [1, 2, 2]

[[0], [1, 2], [3, 4]] [1, 1, 3]

[[0, 3, 4], [1], [2]] [1, 3, 1]

[[0, 4], [1, 3], [2]] [2, 2, 1]

[[0], [1, 3, 4], [2]] [1, 2, 2]

[[0, 4], [1], [2, 3]] [1, 3, 1]

[[0], [1, 4], [2, 3]] [1, 2, 2]

[[0], [1], [2, 3, 4]] [1, 1, 3]

example of cactus length for partition
[[0,2],[1],[3],[4,5,7,10],[6],[8,9],[11]]

[1, 1, 2, 2, 2, 4]

\------------------------------------

m_s[k] for k = 1..10

[0.2730539163373365, 0.4842640687119286, 0.35144954734733014,
0.5249702161277665, 0.44014600164558493,

0.6031675926696468, 0.5443213610399722, 0.7080038062958746,
0.6683603744969959, 0.8395068505692707]

m_s[k] for l = 1..4 by manual application of eq.2

[0.2730539163373365, 0.4842640687119286, 0.35144954734733014,
0.5249702161277665]

\------------------------------------

m[k] for k = 1..10

[0.6, 0.7491168824543142, 0.9949402589451768, 1.3553641553671059,
1.8704311672940084, 2.600699842184038,

3.633212856926222, 5.091635541233709, 7.151079693455817, 10.05923853882072]

m[k] for l = 1..10 by Dubbs&Edelman (Manova moments triangles in Fig.1)

[0.6, 0.7491168824543142, 0.994940258945177, 1.355364155367106,
1.8704311672940082, 2.6006998421840386,

3.633212856926222, 5.091635541233708, 7.1510796934558165, 10.059238538820706]

The output of a symbolic calculation is:

⬇

\------------------------------------

gamma = gamma

p = p

\------------------------------------

x = x

s = s

\------------------------------------

calc ETF moments until d = 10

\------------------------------------

A[l] for l = 1..10

[s, 1, -2*s, 4*s**2 - 1, -8*s**3 + 6*s, 16*s**4 - 24*s**2 + 2,

-32*s**5 + 80*s**3 - 20*s, 64*s**6 - 240*s**4 + 120*s**2 - 5, -128*s**7 + 672*s**5 - 560*s**3 + 70*s,

256*s**8 - 1792*s**6 + 2240*s**4 - 560*s**2 + 14]

A[l] for l = 1..6 by manual application of eq.6

[s, 1, -2*s, 4*s**2 - 1, -8*s**3 + 6*s, 16*s**4 - 24*s**2 + 2]

\------------------------------------

m_s[k] for k = 1..10 [to save space we show only first 6]

[p*s, p**2 + p*s**2, -2*p**3*s + 3*p**2*s + p*s**3,

p**4*(4*s**2 - 1) + p**3*(-8*s**2 + 2) + 6*p**2*s**2 + p*s**4,

p**5*(-8*s**3 + 6*s) + p**4*(5*s*(4*s**2 - 1) - 10*s) + p**3*(-20*s**3 + 10*s)
+ 10*p**2*s**3 + p*s**5,

p**6*(16*s**4 - 24*s**2 + 2) + p**5*(36*s**2 + 6*s*(-8*s**3 + 6*s) - 6) +

p**4*(15*s**2*(4*s**2 - 1) - 60*s**2 + 5) + p**3*(-40*s**4 + 30*s**2) +
15*p**2*s**4 + p*s**6, ...]

m_s[k] for l = 1..4 by manual application of eq.2

[p*s, p**2 + p*s**2, -2*p**3*s + 3*p**2*s + p*s**3,

p**4*(4*s**2 - 1) + p**3*(-8*s**2 + 2) + 6*p**2*s**2 + p*s**4]

\------------------------------------

gamma = gamma

p = p

\------------------------------------

x = x

s = (-x + 1)/(2*sqrt(x))

\------------------------------------

calc ETF moments until d = 10

\------------------------------------

m[k] for k = 1..10

[0.5*p*(-x + 1) + p*(x/2 + 1/2),

1.0*p*(-x + 1)*(x/2 + 1/2) + p*(x/2 + 1/2)**2 + 1.0*x**1.0*(p**2 + p*(-x +
1)**2/(4*x)),

1.5*p*(-x + 1)*(x/2 + 1/2)**2 + p*(x/2 + 1/2)**3 + 3.0*x**1.0*(p**2 + p*(-x +
1)**2/(4*x))*(x/2 + 1/2) +

1.0*x**1.5*(p**3*(x - 1)/sqrt(x) + 3*p**2*(-x + 1)/(2*sqrt(x)) + p*(-x +
1)**3/(8*x**(3/2))),

2.0*p*(-x + 1)*(x/2 + 1/2)**3 + p*(x/2 + 1/2)**4 + 6.0*x**1.0*(p**2 + p*(-x +
1)**2/(4*x))*(x/2 + 1/2)**2 +

4.0*x**1.5*(x/2 + 1/2)*(p**3*(x - 1)/sqrt(x) + 3*p**2*(-x + 1)/(2*sqrt(x)) +
p*(-x + 1)**3/(8*x**(3/2))) +

1.0*x**2.0*(p**4*(-x + (x - 1)**2)/x + p**3*(2 + 2*(-x + 1)*(x - 1)/x) +
3*p**2*(-x + 1)**2/(2*x) +

p*(-x + 1)**4/(16*x**2)), ...]

m[k] for l = 1..10 by Dubbs&Edelman (Manova moments triangles in Fig.1)

[p,

p**2*x + p,

p**3*(x**2 - x) + 3*p**2*x + p,

p**4*(x**3 - 3*x**2 + x) + p**3*(6*x**2 - 4*x) + 6*p**2*x + p, ...]

\------------------------------------

difference between the simplified symbolic expressions of first 10 moments

[0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

## References

  * [1] M. Magsino, D. G. Mixon, and H. Parshall, “Kesten–Mckay law for random subensembles of paley equiangular tight frames,” _Constructive Approximation_ , pp. 1–22, 2020.