# On the estimation of spatial density
from mobile network operator data

Fabio Ricciato and Angelo Coluccia F. Ricciato is with European Commission, DG
Eurostat, Unit B1 Methodology; Innovation in Official Statistics, Luxembourg.
Email<EMAIL_ADDRESS>Coluccia is with University of Salento,
Lecce, Italy. Email<EMAIL_ADDRESS>views expressed in this
paper are those of the authors and do not necessarily represent the official
position of the European Commission.

###### Abstract

Mobile Network Operator (MNO) data are increasingly used to infer mobility and
presence patterns of human population. Whether based on Cell Detail Records
(CDR) or signalling events, MNO data can be leveraged to estimate the spatial
distribution of mobile devices at a given time, and from there extrapolate the
distribution of humans. The process of transforming MNO data to a density map
involves a chain of multiple processing blocks. A key block relates to _geo-
location_ of individual radio cells or groups thereof, i.e., to methods for
determining the spatial footprint of radio cells. Traditionally, researchers
have resorted to geo-location methods based on Voronoi tessellations or
variants thereof: with this approach, cell locations are mutually disjoint and
the density estimation task reduces to a simple area-proportional computation.
More recently, some pioneering work have started to consider more elaborated
geo-location methods with partially overlapping (non-disjont) cell footprints.
Estimating the spatial density from a set of overlapping cell locations is
currently an open research problem, and it is the focus of this work. In this
contribution we start by reviewing three different estimation methods proposed
in the literature, for which we provide novel analytic results, based on
formal proofs, unveiling some key aspects of their mutual relationships.
Furthermore, we develop a novel estimator for which a closed-form solution can
be given. Numerical results based on synthetic data are presented to assess
the relative accuracy of each method.

##  I Introduction

Most people nowadays carry a mobile phone, and mobile phones interact several
times a day with the mobile network infrastructure. Every interaction between
the mobile device and the network reveals the approximate location of the
device, hence of the associated mobile user, at least at radio cell level.
Such interactions are recorded by the Mobile Network Operator (MNO) for
purposes related to the delivery of mobile communication services (e.g.
billing, network optimization and troubleshooting). Since more than two
decades, researchers have (re)used MNO data to study patterns of human
mobility and presence, effectively exploiting the mobile network
infrastructure as a large-scale “sensor of opportunity” for human mobility
flows (see e.g. [1] and references therein). More recently, MNO data have been
used to analyse human mobility in the context of Covid-19 pandemic and
associated containment measures [2, 3] notwithstanding a number of open
challenges in accessing and making use of such data [4]. Beyond the scope of
academic research, commercial companies and non-governmental organization are
offering analytics products based on MNO data in various application domains,
from humanitarian support to tourism, and statistical organizations are
looking with increasing interest at MNO data as a potential source for
compiling new official statistics [5, 6, 7]. However, despite the great volume
of research literature on the topic, several methodological aspects remain
open along the journey from raw MNO data to reliable summary information.

The focus of this contribution is on the problem of measuring the spatial
density of mobile users, at a given reference time, based on the data from a
single MNO. With proper re-scaling, the density of mobile users provides a
rough indication of the population present in the territory around that time,
also called “present population”, “hourly population” or _de facto_
population. The data processing flow, from raw MNO data to final density
estimates, involves a pipeline of different processing blocks, and current
efforts within the European Statistical System aim at elaborating a unifying
view of such processing pipeline, that is, a common methodological framework
(see [8] and also [7]). One key block along the pipeline relates to the geo-
location of radio cells, i.e., the mapping of radio cells (or groups thereof)
to the geographical territory. The geo-location block can be instantiated with
different methods with varying levels of sophistication. Choosing a more
elaborated method accounts to attempting a more detailed (and complex)
modeling of the mobile network infrastructure in order to potentially (but not
necessarily) achieve better accuracy.

All possible geo-location methods can be grouped into two classes,
discriminated by whether the set of cell locations delivered in output are
spatially disjoint or overlapping. The methods involving disjoint (non-
overlapping) cell locations, or “tessellations”, are simpler to implement and
by far more popular in the past research literature. Only recently a few
pioneering papers started to consider alternative, more elaborated schemes
based on overlapping cell locations [9, 10, 11].

A separate block along the pipeline, logically subsequent to the geo-location
block, performs the task of estimating (or inferring) the underlying spatial
density from the set of geo-localized data. There is a fundamental inter-
dependency between these two blocks: if the geo-localization method of choice
belongs to the specific class of tessellations, then the density estimation
reduces to a simple area-proportional solution. Conversely, if the geo-
localization method belongs to the more general class of overlapping
locations, determining the “best” estimation method is an open research
problem: this is indeed the focus of our paper.

Given this background, we make here the following contributions:

  * •

We provide a coherent representation of the general processing workflow, from
raw MNO data to spatial density estimates, with a clear definition of the geo-
location and estimation blocks;

  * •

Focusing on estimation methods for overlapping cell locations, we review three
distinct solutions from previous literature and provide novel analytic insight
about their mutual relationships;

  * •

We derive analytically the Maximum A Posteriori (MAP) estimator for the
problem at hand and, from there, elaborate a completely novel “data first”
estimator with closed-form solution;

  * •

We provide quantitative insight by comparing the accuracy of the different
solutions on simple synthetic scenarios.

Along the way, we identify directions for additional research.

The rest of this paper is organized as follows. We start by presenting a
general terminology and a comprehensive description of the methodological
framework underlying our work in Sec. II. Based on that, we define the system
model and formulate the estimation problem analytically in Sec. III. In Sec.
IV we review the existing estimation solutions that were proposed
independently in the recent literature, and in the following Sec. V we present
novel analytical results that re-interpret these schemes and provide insight
about their mutual relationships. In Sec. VI we derive the MAP estimator for
the problem at hand and, starting from its structure, in the following Sec.
VII we propose a completely novel estimator. Illustrative numerical results
based on synthetic data are presented in section VIII. Finally, in Sec. IX we
conclude.

##  II A general introduction to MNO data

In this section we provide a general view of “MNO data” and the key processing
steps involved in the extraction of spatial density estimates from such data
(the interested reader is referred to [7, 8, 12] for a more detailed
elaboration of the associated methodological issues).

###  II-A Event data and cell data

It is convenient to distinguish two kinds of input data that must be
necessarily combined in any study of human presence or mobility based on MNO
data:

  * •

Event data. These data result from transaction events between the mobile
device and the network infrastructure111The term “transactional data” was used
in an earlier study [12] as synonymous of “event data”.. A mobile device $z$
that at time $t$ interacts with the mobile network through the radio cell $i$
generates an event record $\langle z,t,i\rangle$. In practice $z$ represents
the device pseudonym (typically based on non-reversible hash of the
International Mobile Subscriber Identifier, IMSI), $t$ the event timestamp and
$i$ the radio cell identifier. Depending on how the data collection system is
configured, only a subset of all events are captured into the data set. Call
Detail Records (CDR) and the more informative signalling data fall in this
category222Under certain conditions (e.g. for signaling data captured on the
Radio Access Network links) additional variables may be available in addition
to the cell identifier $i$ that could help to identify the mobile user
location at sub-cell level, e.g. the so-called Timing Advance (TA) that allows
to determine upper and lower bounds on the distance between the device and the
base station antenna, or radio measurements from neighboring cells (as found
e.g. in Location Based Systems (LBS)). These additional variable(s) should be
considered absorbed by the the data element $i$ that, in such cases, should be
reinterpreted as the union of the cell identifier and any other available
variable carrying spatial information (including TA or LBS). The use of
additional information from the radio interface is still relatively infrequent
(see [13] for a prominent example) but could become more important in future
5G deployments. However, to ease the presentation, we do not elaborate further
this scenario and maintain the variable $i$ to denote solely the radio cell
identifier..

  * •

Cell data. Such data binds the network infrastructure, i.e., the collection of
“radio cells”, to the geographical territory333The term “auxiliary data” was
used in [12] as synonymous of “cell data”.. Considering a discretization of
the geographical territory (e.g. in a regular square grid), cell data might be
thought as the collection of records in the form $\langle
t,i,\bm{s}_{i}(t)\rangle$ associating to each cell $i$ at time444The temporal
dimension cannot be ignored for cell data since network coverage is not static
but changes in time (also) due to the continuous addition and deletion of
radio cell, to the manual and/or automatic (re)optimization of radio
configuration parameters, and to the physiological change of traffic load
(cell breathing in 3G).  $t$ a geographical profile, or “cell footprint”,
represented by the vector $\bm{s}_{i}$ (also called “event location” in [7]).
The $j$th element of such vector, denoted by $s_{ij}\geq 0$, encodes the
degree by which cell $i$ is expected to “cover” the grid point $j$. Depending
on the assumed geo-location model, the variables $s_{ij}$’s can be restricted
to take binary values, thus leading to on/off flat coverage patterns, or
alternatively they could let take continuous values aimed at representing the
strength of the cell signal in that area (an approach first proposed in [11,
14]).

Following [8] we call _geo-location_ the task of determining the spatial
footprint of each radio cell at a given time. Radio cells that have identical
footprints are grouped together, with index $i$ now referring to the whole
cell group. Clearly, using the pair of time and cell (group) identifier
$\langle t,i\rangle$ as matching key to combine event records $\langle
z,t,i\rangle$ and cell records $\langle t,i,\bm{s}_{i}(t)\rangle$, we can geo-
locate individual events, hence mobile users. Therefore, for a given reference
time or time window, we can count the number $c_{i}$ of devices that were
located in cell (or cell group) $i$ having footprint $\bm{s}_{i}$. Then the
set of counters and footprints over all cell (or cell groups)
$\left\\{c_{i},\bm{s}_{i}\right\\}_{i}$ represents the input for the density
estimation problem that is the focus of our contribution.

There are numerous geo-location methods possible, each making use of more or
less information and with different degrees of sophistication. Every geo-
location method underlies, explicitly or implicitly, an assumption about the
process by which the generic mobile device selects the radio cell to connect.
Generally speaking, we can classify all possible approaches into two large
families:

  * •

Tessellations: methods that partition the territory into a set of non-
overlapping (disjoint) footprints associated to different cells (or groups
thereof).

  * •

Overlapping cell locations: for which cell footprints in general overlap.

The family of tessellation approaches may be seen as a particular (degenerate)
case of the more general overlapping cells approach.

###  II-B Geo-location methods based on tessellations

The tessellations methods that have been considered in the literature on MNO
data typically resort to the principle of Voronoi partition. We recall that,
for a discrete set of $K$ points (Voronoi seeds) in a bounded planar region,
the Voronoi partition divides the region into a set of $K$ non-overlapping
sub-regions (Voronoi polygons) whereby each point is space is associated to
the closest seed. The Voronoi partition is well defined, in the sense that a
single partition (Voronoi diagram) exists for a given set of seeds in the
bounded region. However, there are multiple ways of projecting a set of mobile
radio cells into a set of Voronoi seeds (and multiple ways of bounding a
region of interest) and each method results in a different Voronoi diagram for
the same input set of cells.

The simplest — and by far the most popular — method takes antenna tower
positions as seeds. As in real networks multiple radio cells have their
antennas co-located on a single antenna tower, the number of unique seeds,
hence polygons, is lower than the number of radio cells. This method does not
require any information about the radio cell configuration other than the
antenna tower position, and therefore is very simple to implement. The
implicit modeling assumption underlying this method is that mobile devices
always connect deterministically to the closest antenna tower — an assumption
that ignores several fundamental aspects of mobile communications, including
antenna directionality, power control, load balancing, multi-layer radio
deployments, etc.

Another variant, first considered in [15, 9, 10] and more recently also in
[16], places seeds at the barycenter of radio cell coverage area. This
requires additional knowledge about the radio cell configuration parameters
(e.g., azimuth orientation, beam width, coverage range) in order to determine,
at least roughly, the cell coverage area.

The method presented by the Belgian MNO Proximus in [6] is more articulated.
First, they distinguish between large cells and small cells (femto-cells and
pico-cells) and apply Voronoi partition only to the large cells. Second, in
order to take into account the directionality of cell sectors, they places $N$
Voronoi seeds in the vicinity of each $N$-sector antenna (typically $N=3$ for
$120^{\circ}$ sectors) with a small offset in their respective azimuth
directions.

Generally speaking, all variants of the Voronoi approach have the disadvantage
that even a single cell change (i.e., addition, removal or shift of its
associated seed) produces a change in the neighboring cell footprints, and in
principle requires a new computation of the entire Voronoi diagram. Departing
from the Voronoi approach, other forms of tessellations may be obtained by
mapping each cell to one (and only one) particular unit of a predefined
partition, e.g. a regular square grid (as in [17]) or a variable resolution
quadtree [18].

Whether based on Voronoi or an independent fixed grid, the more sophisticated
variants of tessellation methods require additional information about the
radio cell configuration, beyond the antenna location. However, as more
detailed cell data are made available for the geo-location process, the
limitation of considering non-overlapping footprints seems to be less
justified, motivating the interest for methods accounting for the intrinsic
overlapping nature of real-world radio cells.

###  II-C Geo-location methods based on overlapping cell locations

This family of methods takes explicitly into account the fact that real-world
radio cells overlap by design. To the best of our knowledge, the first work to
consider an overlapping cells approach for the problem of density estimation
was proposed by Ricciato et al. in [9, 10] (see also the earlier work [19] for
a different application). Therein, the authors consider cell footprints that
are “flat”: for a given radio cell $i$, a generic point $j$ is either included
or excluded from its coverage area. This model implicitly assumes a on/off
reception model, therefore the coverage parameters
$s_{ij}\in\left\\{0,1\right\\}$ are binary. This method will be denoted
hereafter as “OF” for “Overlapping cells with Flat footprint”.

More recently, a more sophisticated approach was developed by Tennekes et al.
in [11] (see also the earlier presentation [14]) and implemented in the mobloc
R package [20]. This approach considers a non-uniform profile whereby the
parameters $s_{ij}$’s are continuous and vary within the cell footprint,
having a physical interpretation connected to signal strength (i.e., received
signal power). This method will be denoted hereafter as “OV” for “Overlapping
cells with Variable footprint”.

The geo-location methods OF and OV introduced above represent two points in a
larger space of possible approaches. For all possible proposals, the
determination of the cell footprints $\bm{s}_{i}$’s is a critical modeling
tasks: the set of possible implementation solutions is wide, ranging from
elementary geometric models (as done in [9, 10]) to more articulated but still
parsimonious parametric propagation models (as done in [11]), and even more
sophisticated (and less parsimonious) empirical numerical models that take
into account more detailed data about the territory (elevation maps, type of
building clutter) and about the radio network dynamics (power control, inter-
cell interference etc.), possibly reusing data and tools that are already
available and used for radio network planning and optimization tasks.

Regardless of how the cell footprints $\bm{s}$’s are empirically determined,
the next step in the modelling process is to encode such information into a
probabilistic data generation model that will then serve as basis for the
development of an estimation (inference) method. In the next section we fist
introduce the probabilistic model, and then explain how this model is linked
to the geo-location method of choice in the specific context of MNO data
analysis.

##  III Probabilistic model and problem formulation

For the rest of this paper we assume the geographic territory is discretized
into a regular square grid. We focus on the estimation problem at a given
reference time, therefore the temporal dimension can be dropped from the
formalism. To avoid excessive overload of terms, we resort to the term “tile”
to refer to the individual square grid unit, reserving the term “cell” to
denote the “radio cell” constituting the mobile (cellular) radio network.
Also, we shall use the terms “user”, “mobile user” and “mobile device”
interchangeably.

Let the $j$th element $u_{j}$ of the column vector
$\bm{u}\stackrel{{\scriptstyle\textrm{\tiny def}}}{{=}}[u_{1}\ u_{2}\ \cdots\
u_{J}]^{\mathsf{T}}$ denote the unknown number of mobile users in tile
$j=1,\ldots,J$. Let the $i$th element $c_{i}$ of the column vector
$\bm{c}\stackrel{{\scriptstyle\textrm{\tiny def}}}{{=}}[c_{1}\ c_{2}\ \cdots\
c_{I}]^{\mathsf{T}}$ denote the observed number of mobile users detected (or
counted) in cell $i=1,\ldots,I$ (cell count vector). We denote the total
number of users across all cells by $C\stackrel{{\scriptstyle\textrm{\tiny
def}}}{{=}}\sum_{i=1}^{I}{c_{i}}=\|\bm{c}\|_{1}=\bm{1}_{I}^{\mathsf{T}}\bm{c}$
where the symbol $\bm{1}_{k}$ denotes a column vector of size $k$ with all
elements equal to 1. Denote by $p_{ij}\in[0,1]$ the probability that a generic
device placed in tile $j$ will be detected (counted) in cell $i$. In other
words, $p_{ij}$ represents the conditional probability:

$p_{ij}\stackrel{{\scriptstyle\textrm{\tiny
def}}}{{=}}\text{Prob}\left\\{\text{detected in cell $i$}\;|\;\text{placed in
tile $j$}\right\\}.$ (1)

For the sake of a more compact notation we gather the individual probabilities
$p_{ij}$’s into a matrix $\bm{P}_{[I\times J]}$. The individual elements of
$\bm{P}$, which correspond to the _emission probabilities_ in the field of
emission tomography [21], will be referred to as _assignment probabilities_
throughout the paper. The elements of $\bm{P}$ sum to one along columns,
formally ${\bm{1}}_{I}^{\mathsf{T}}\bm{P}={\bm{1}}_{J}$, meaning that the
matrix $\bm{P}$ is column-stochastic555In this work we are interested to
determine the spatial distribution of all and only the mobile devices that
have connected to the network, not individual people. In other words, missing
observations (i.e., people that are not counted in any radio cell) and other
population coverage errors (i.e., users that are counted in multiple cells, or
spurious observations that do not correspond to a human) are not considered in
the model at this stage and are left for consideration in further extensions..
We consider for this work the ideal case that the matrix $\bm{P}$ is known
exactly. In other words, do not consider model mismatching errors due to
inaccurate knowledge of $\bm{P}$ at this stage. While a certain degree of
model mismatching is unavoidable in practical applications, we will show that
even the ideal case of perfectly known generative model is challenging (and
interesting) to address. This motivates leaving the issue of (sensitivity to)
model mismatching errors as a point for further study, outside the scope of
the present contribution.

The (measured) cell count vector $\bm{c}$ can be interpreted as the single
realization of a random vector $\tilde{\bm{c}}$ whose expected value is given
by:

$\overline{\bm{c}}\stackrel{{\scriptstyle\textrm{\tiny
def}}}{{=}}\mathrm{E}[\tilde{\bm{c}}]=\bm{P}\bm{u}.$ (2)

In the estimation problem we must solve for estimand $\bm{u}$ given the vector
of measurement data $\bm{c}$, representing the single available observation of
$\tilde{\bm{c}}$, and the model matrix $\bm{P}$. Being a sort of inversion
problem, the estimate $\hat{\bm{u}}$ can be written in general as:

$\hat{\bm{u}}=g(\bm{P},\bm{c})$ (3)

where $g(\cdot)$ denotes the estimator of choice. It is evident that in the
estimator design we must constrain the estimand variables to be non-negative,
i.e., $u_{i}\geq 0,\;\forall i$. However, we relax the constraint that such
variables should be integer since the rounding error can be safely neglected
vis-à-vis other sources of uncertainty.

######  Remark 1.

In cases of practical interest the number of tiles is (much) larger than the
number of cells, i.e., $J\gg I$. Therefore, even if the average value
$\overline{\bm{c}}$ were perfectly known, the direct inversion of (2) would
constitute an under-determined problem. For that reason, it should not come as
a surprise that, as we show below, the estimation problem in this case is
affected by issues of _structural non-identifiability_ (by the definition
given in [22]). Any additional external information that is available to help
the estimation process (e.g., prior distributions, spatial constraints derived
from geographical maps, known structural properties of the desired solution)
can be embedded in the estimator $g(\cdot)$ to resolve, or at least reduce,
the ambiguity among multiple solutions, as elaborated later.

If two generic tiles $j_{1}$ and $j_{2}$ yield equal assignment probabilities
for all cells, i.e., $p_{ij_{1}}=p_{ij_{2}}\;\forall i$, and are associated to
the same prior values (in case that prior information is used in the estimate)
then they are indistinguishable from each other. In such case, their
respective estimand variables $u_{j_{1}}$ and $u_{j_{2}}$ are perfectly
_collinear_ in the estimation problem and there is no way to resolve
differences between them. Therefore, it makes sense to merge both tiles into a
single super-tile, and then compute a single estimate for the latter666The
concept of super-tile that is introduced here algebraically corresponds to the
concept of _section_ introduced geometrically in the earlier study [10]. Here
we use the terms “section” and “super-tile” interchangeably.. Algebraically,
this corresponds to merging together identical columns of matrix $\bm{P}$. We
refer to this operation by the term consolidation. Note that the consolidation
process does not imply that the resulting (consolidated) matrix is necessarily
full rank, i.e., it does not guarantee the resulting consolidated instance of
the problem is fully identifiable.

(a) Multinomial

(b) Poisson

Figure 1: Hierarchical generative models

For the sake of completeness, we briefly illustrate hereafter how the
assignment matrix $\bm{P}$ can be derived from the cell footprints
$\bm{s}_{i}$’s. In the more general OV method [11], the variable $s_{ij}$ are
called “signal dominance” and are derived from the predicted strength of the
radio signal (power received for cell $i$ at tile $j$) through a logistic
transformation. The model assumes that a cell is selected randomly with a
probability that is proportional to the _relative_ signal dominance over the
other concurrent cells in the same tile, leading to

$p_{ij}=\frac{s_{ij}}{\sum_{n}{s_{nj}}}.$ (4)

The OF method can be considered as a particularization of the more general OV
approach. With OF, an on/off coverage model is assumed: the variables
$s_{ij}$’s are binary and all cells covering tile $j$ have the same
probability of being selected, therefore the non-zero elements of (4) take
fractional values:

$s_{ij}\in\\{0,1\\}\Rightarrow p_{ij}\in\left\\{0,\frac{1}{k_{j}}\right\\}$

with $k_{j}\stackrel{{\scriptstyle\textrm{\tiny def}}}{{=}}\sum_{n}{s_{nj}}$
denoting the (integer) number of concurrent cells at point $j$.

Tessellations fall into an even more special case of the above: as each point
in space is covered by exactly one and only one cell ($k_{j}=1$) the
assignment probabilities reduce to binary values
$p_{ij}\in\left\\{0,1\right\\}$. In other words, when cells are mutually non-
overlapping, the data generating process becomes deterministic: the binary
model matrix $\bm{P}$ after consolidation reduces to the identity matrix, and
the estimation problem becomes trivial777For this reason, the estimation
problem defined in the present contribution was not considered in past
literature dealing with MNO data that, aside of a few exceptions, ended up
considering implicitly or explicitly some variant of Voronoi tessellation for
the geo-location stage. as far as no external information is taken into
account (such as, e.g., prior information in Bayesian settings as considered
in [23]). In this contribution we focus on the more general (and non trivial)
case of overlapping cells, with the understanding that the proposed solution
will be applicable also to the special case of tessellations.

##  IV Estimation methods from previous literature

Whatever approach is chosen to populate the model matrix $\bm{P}$ during the
_modeling stage_ , a suitable estimation procedure (represented by the
function $g(\cdot)$ in (3)) is needed for the subsequent _resolution stage_.
This is the focus of the remaining part of this contribution. In this section
we present three different solutions proposed in the recent literature. To the
best of our knowledge, no other solution was previously considered for the
problem at hand.

###  IV-A MLE-Multinomial

The method elaborated in [10, Section 5.3] (previously appearing in [9])
derives a Maximum Likelihood Estimator (MLE) based on the hierarchical
generative model sketched in Fig. 1(a). In the first layer, starting from a
single root pool, the $C$ devices are randomly and independently allocated to
the $J$ tiles according to a vector of probabilities
$\bm{\mu}\stackrel{{\scriptstyle\textrm{\tiny def}}}{{=}}[\mu_{1}\ \cdots\
\mu_{J}]^{\mathsf{T}}$. Clearly, $\bm{\mu}$ is a vector with non-negative
elements summing to unity, i.e., $\bm{\mu}\geq\bm{0}$ and
$\|\bm{\mu}\|_{1}=\bm{1}_{J}^{\mathsf{T}}\bm{\mu}=1$. The resulting number of
units in the tiles is distributed as a multivariate random vector
$\tilde{\bm{u}}$ with Multinomial distribution (denoted by $\mathcal{M}$) of
parameters $C$ and $\bm{\mu}$, i.e.

$\tilde{\bm{u}}\sim\mathcal{M}\left(C,\bm{\mu}\right)\stackrel{{\scriptstyle\textrm{\tiny
def}}}{{=}}\frac{C!}{\tilde{u}_{1}!\tilde{u}_{2}!\cdots\tilde{u}_{J}!}\prod_{j=1}^{J}{\left(\mu_{j}\right)^{\tilde{u}_{j}}}$
(5)

with $\tilde{\bm{u}}\geq\bm{0}$ and $\|\tilde{\bm{u}}\|_{1}=C$. Therefore by
construction the mean value writes
$\overline{\bm{u}}\stackrel{{\scriptstyle\textrm{\tiny
def}}}{{=}}\mathrm{E}[\tilde{\bm{u}}\,|\,C,\bm{\mu}]=C\bm{\mu}$. In the second
layer (ref. Fig. 1(a)) the units are assigned randomly and independently from
tiles to cells according to the assignment matrix $\bm{P}$. Due to
independence of the random assignments at the two layers — from the root pool
to tiles, and from tiles to cells — the random vector $\tilde{\bm{c}}$ has
also a Multinomial distribution with parameters $C$ and $\bm{P}\bm{\mu}$:

$\tilde{\bm{c}}\sim\mathcal{M}\left(C,\bm{P}\bm{\mu}\right)\stackrel{{\scriptstyle\textrm{\tiny
def}}}{{=}}\frac{C!}{\tilde{c}_{1}!\tilde{c}_{2}!\cdots\tilde{c}_{I}!}\prod_{i=1}^{I}{\left(\sum_{j=1}^{J}{\mu_{j}\,p_{ij}}\right)^{\tilde{c}_{i}}}.$
(6)

Thus, the log-likelihood function for an observed value $\bm{c}$ is derived as
(neglecting an irrelevant constant factor):

$\mathcal{L}(\bm{\mu};\bm{c},\bm{P})=\sum_{i=1}^{I}{c_{i}\log{\sum_{j=1}^{J}{\mu_{j}p_{ij}}}}=\bm{c}^{\mathsf{T}}\,\log{\bm{P}\bm{\mu}}$
(7)

wherein we have used the logarithm of vector notation
$\log{\bm{x}}\stackrel{{\scriptstyle\textrm{\tiny
def}}}{{=}}\left[\log{x_{1}}\ \log{x_{2}}\ \cdots\right]^{\mathsf{T}}$ to
refer compactly to the vector of element-wise logarithms. Let
$\hat{\bm{\mu}}_{\text{ML}}$ denote the value of the probability vector
$\bm{\mu}$ that maximizes the log-likelihood (7) subject to the constraints
$\bm{\mu}\geq\bm{0}$ and $\|\bm{\mu}\|_{1}=1$, formally:

$\hat{\bm{\mu}}_{\text{ML}}=\operatorname*{arg\,max}_{\begin{subarray}{c}\|\bm{\mu}\|_{1}=1\\\
\bm{\mu}\geq\bm{0}\end{subarray}}{\,\bm{c}^{\mathsf{T}}\log{\bm{P}\bm{\mu}}}.$
(8)

By definition $\hat{\bm{\mu}}_{\text{ML}}$ represents the Maximum Likelihood
(ML) estimate for $\bm{\mu}$ given the observed data $\bm{c}$. Given
$\hat{\bm{\mu}}_{\text{ML}}$, we still have to provide a point estimate for
$\bm{u}$. A quite natural choice (also taken by Shepp and Vardi in their
seminal paper [21] for the method presented in the next subsection) is to take
the corresponding _expected mean value_ as the final estimate, formally:

$\hat{\bm{u}}_{\text{ML}}\stackrel{{\scriptstyle\textrm{\tiny
def}}}{{=}}\mathrm{E}[\tilde{\bm{u}}\,|\,C,\bm{\hat{}}{\bm{\mu}}_{\text{ML}}]=C\,\bm{\hat{}}{\bm{\mu}}_{\text{ML}}$
(9)

where we retain the label “ML” for simplicity. This equation represents a
simple rescaling of $\hat{\bm{\mu}}_{\text{ML}}$ by the factor $C$. Therefore,
with a simple variable substitution, the minimization (8) can be rewritten to
deliver directly the final estimate $\hat{\bm{u}}_{\text{ML}}$, leading to the
following constrained optimization problem:

$\hat{\bm{u}}_{\text{ML}}=\operatorname*{arg\,max}_{\begin{subarray}{c}\|\bm{u}\|_{1}=C\\\
\bm{u}\geq\bm{0}\end{subarray}}{\,\bm{c}^{\mathsf{T}}\log{\bm{P}\bm{u}}}$ (10)

In the original papers [9, 10] the minimization (10) is conducted via standard
numerical solvers. No remark is made therein about the uniqueness (or lack
thereof) of the solution.

###  IV-B MLE-Poisson

Recently, the authors of [24] have considered to apply to this problem the ML
estimator developed in the field of emission tomography by Shepp and Vardi in
[21]. Like the previous approach, also this method is based on a hierarchical
generative model with two layers, as sketched in Fig. 1(b), but now the
elements of $\tilde{\bm{u}}$ are modelled as independent Poisson (instead of
Multinomial) random variables with parameters
$\bm{\rho}\stackrel{{\scriptstyle\textrm{\tiny def}}}{{=}}[\rho_{1}\ \cdots\
\rho_{J}]^{\mathsf{T}}$, i.e., $\tilde{u}_{j}\sim\mathcal{P}(\rho_{j})$
($\mathcal{P}$ denoting the Poisson distribution). Following the same
reasoning that led to (9), the ML estimate of $\bm{\rho}$ is taken as the
estimate for $\bm{u}$ (see [21, p. 113]). The log-likelihood function is not
given explicitly in [21] but it is proved that the ML estimate can be computed
iteratively through an Expectation Maximization (EM) procedure: at the generic
iteration $m$ the new estimate $\hat{u}_{j}^{m+1}$ is computed from the
previous estimate $\hat{u}_{j}^{m}$ according to the following formula (see
[21, eq. (2.13)] or, equivalently, [24, eq. (2)]):

$\hat{u}_{j}^{m+1}=\hat{u}_{j}^{m}\cdot\sum_{i=1}^{I}{c_{i}\frac{p_{ij}}{\sum_{k=1}^{J}{p_{ik}\hat{u}_{k}^{m}}}}.$
(11)

The authors of the original paper warn that the initialization point should
not contain zero elements, and by default they assume a flat (uniform) initial
solution $\hat{u}_{j}^{0}=\frac{C}{J}$, $\forall j$.

###  IV-C A simple estimator based on Bayes rule (SB)

The simple procedure presented hereafter was adopted in the mobloc R package
developed by Tennekes et al. [20] and elaborated in [11] (see also [14]). We
shall refer to this method as the “Simple Bayes-rule estimator” (SB for
short). Let $q_{ji}$ denote the conditional probability:

$q_{ji}\stackrel{{\scriptstyle\textrm{\tiny
def}}}{{=}}\text{Prob}\left\\{\text{placed in tile $j$}\;|\;\text{detected in
cell $i$}\right\\}.$ (12)

Note the inversion of the conditioning direction between $q_{ji}$ and $p_{ij}$
defined earlier in (1). Let

$\alpha_{j}\stackrel{{\scriptstyle\textrm{\tiny
def}}}{{=}}\text{Prob}\left\\{\text{placed in tile $j$}\right\\}$

denote the prior probability that a single generic unit falls in tile $j$,
before observing the measurement data. Recalling the Bayes rule

$\text{Prob}\left\\{j\,|\,i\right\\}=\frac{\text{Prob}\left\\{i\,|\,j\right\\}\cdot\text{Prob}\left\\{j\right\\}}{\text{Prob}\left\\{i\right\\}}$

it follows that

$q_{ji}=\frac{p_{ij}\,\alpha_{j}}{\sum_{k=1}^{J}{p_{ik}\,\alpha_{k}}}.$

Therefore, the estimate in each tile $j$ is computed directly as

$\hat{u}_{j}=\sum_{i=1}^{I}{q_{ji}\,c_{i}}=\alpha_{j}\,\sum_{i=1}^{I}{c_{i}\frac{p_{ij}}{\sum_{k=1}^{J}{p_{ik}\,\alpha_{k}}}}.$
(13)

######  Proposition 1.

The simple Bayes-rule estimator (13) belongs to the linear-type subclass of
(3); in fact, it can be rewritten in vector form as

$\hat{\bm{u}}_{\text{SB}}=\bm{Q}\bm{c},\quad\bm{Q}_{[J\times
I]}\stackrel{{\scriptstyle\textrm{\tiny
def}}}{{=}}\mathrm{diag}(\bm{\alpha})\bm{P}^{\mathsf{T}}\mathrm{diag}^{-1}(\bm{P}\bm{\alpha})$
(14)

where $\mathrm{diag}(\bm{v})$ is a diagonal matrix containing the entries of
vector $\bm{v}$ (and $\mathrm{diag}^{-1}(\bm{v})$ their reciprocals).

###### Proof.

We first notice that $\sum_{k=1}^{J}p_{ik}\alpha_{k}$ and
$\sum_{i=1}^{I}\xi_{i}p_{ij}$ are the $i$th and $j$th element, respectively,
of the vectors $\bm{P}\bm{\alpha}$ and $\bm{P}^{\mathsf{T}}\bm{\xi}$,
respectively, with $\xi_{i}=c_{i}/\sum_{k=1}^{J}p_{ik}\alpha_{k}$. Then the
thesis follows straight by rewriting in vector form
$\bm{\xi}=\mathrm{diag}^{-1}(\bm{P}\bm{\alpha})\bm{c}$ and exploiting the
associative property of the matrix product. ∎

Proposition 1 highlights that, compared to the MLE approach, the SB is
appealing due to its simplicity, ease of implementation, and low computational
cost: through direct computation it provides immediately a single non-
ambiguous solution for the given input model $\bm{P}$ and data $\bm{c}$.

###  IV-D Remarks on initial guess

All the three methods reviewed above provide a point estimate $\hat{\bm{u}}$
in output, based on the model matrix $\bm{P}$ and observed data $\bm{c}$ in
input. To do so, they all require implicitly or explicitly the provision of an
_initial point_ as input to the computation. The initial point takes either
the form of a (stochastic) vector of $J$ non-negative elements summing to
unity, hereafter denoted by $\bm{\alpha}\stackrel{{\scriptstyle\textrm{\tiny
def}}}{{=}}[\alpha_{1}\ \cdots\ \alpha_{J}]^{\mathsf{T}}$ (with
$\|\bm{\alpha}\|_{1}=1$ and $\bm{\alpha}\geq\bm{0}$) or, equivalently, of its
rescaled version $\bm{a}\stackrel{{\scriptstyle\textrm{\tiny
def}}}{{=}}C\bm{\alpha}$ (with $\|\bm{a}\|_{1}=C$ and $\bm{a}\geq\bm{0}$).
Such initial vector represents our ex-ante “best guess” about how the mobile
units may be spatially distributed, _before seeing the data_ ; in this sense,
as already noted by Shepp and Vardi888Quoting from [21]: “As a point of
philosophical interest, the choice of the initial estimate is somewhat akin to
the choice of a Bayes prior but there is actually no prior measure.”, it
represents a sort of prior information, and for this reason we shall refer to
$\bm{\alpha}$ (or equivalently to $\bm{a}$) as the _prior vector_. Note
however the difference between the notions of _prior vector_ and _prior
distribution_ : the former can serve as parameter for the latter.

The role of the prior vector $\bm{\alpha}$ (non rescaled) is explicit in the
SB estimator (14) where its elements represent prior probabilities (per
tiles). As for MLE-Poisson, the prior vector $\bm{a}$ (rescaled) serves as the
initial point for the iterative procedure (11) by setting
$\hat{\bm{u}}^{0}=\bm{a}$. Similarly, the prior vector $\bm{a}$ (rescaled) is
the natural starting point for the numerical minimization process (10) to
compute the MLE-Multinomial estimate. In summary, all three methods rely on a
prior vector (rescaled or not) as initial point for the computation.

A point of caution is needed in determining the value of the prior vector,
particularly for what concerns its zero elements. It can be immediately
recognized that setting the initial value $\alpha_{j}=0$ (or equivalently
$a_{j}=0$) for a generic tile $j$ will force the final solution to maintain a
zero value in the same tile in the EM iterative procedure (11) (i.e., zeros
are stable points of the iteration) as well as in the direct computation of SB
based on (14). In other words, with both these methods (namely SB and MLE-
Poisson with EM), zeroing the $j$th element in the prior vector overrides any
possible contribution carried by the measurement data for the corresponding
tile $j$. This is equivalent to completely excluding _a priori_ the
corresponding tile $j$ from the computation. But if we intend to do so, it is
certainly more practical to drop upfront the variable associated to such tile,
rather than instantiating a variable whose final value is fixed ex-ante to
zero. By eliminating upfront the zero tiles, we can practically assume that
all remaining elements in the prior vector take non-zero values.

##  V New insights on previous estimators

In this section we present our new results and reinterpret the previous
solutions. More specifically, the analytic insight presented in this section
will serve as basis for the development of a novel estimation method in a
later section.

###  V-A MLE-Multinomial and MLE-Poisson are equivalent

We show that the two ML procedures presented earlier in section IV-A and IV-B
are equivalent in the sense that they yield exactly the same set of solutions.
Beforehand, we need a technical lemma summarized in the following Proposition.

######  Proposition 2.

The negative log-likelihood (7) is convex, but not strictly convex and the
constraints are linear, so the overall problem is convex but not strictly
convex. As a consequence, all stationary points are equivalent global maxima.

###### Proof.

It is sufficient to show the positive semidefiniteness of the Hessian of the
log-likelihood (7). The generic $(k,h)$ element of the Hessian is given by the
second derivatives, i.e.

$H_{k,h}=-\sum_{i=1}^{I}c_{i}\frac{p_{ik}p_{ih}}{(\bm{p}_{i}^{\mathsf{T}}\bm{u})^{2}}$

wherein $\bm{p}_{i}^{\mathsf{T}}$ denotes the $i$th row of $\bm{P}$. Thus, for
any vector $\bm{x}=[x_{1}\cdots x_{I}]^{\mathsf{T}}$ and assuming $c_{i}>0$
(possible values $c_{i}=0$ do not count for the sake of the argument), by
recalling that the elements of $\bm{P}$ are non-negative (of which at least
one strictly positive element for each column), we can write

$\displaystyle\sum_{k,h=1}^{J}x_{k}x_{h}H_{k,h}$
$\displaystyle=-\sum_{k,h=1}^{J}x_{k}x_{h}\sum_{i=1}^{I}c_{i}\frac{p_{ik}p_{ih}}{(\bm{p}_{i}^{\mathsf{T}}\bm{u})^{2}}$
$\displaystyle=-\sum_{i=1}^{I}\frac{c_{i}}{(\bm{p}_{i}^{\mathsf{T}}\bm{u})^{2}}\left(\sum_{k=1}^{J}p_{ik}x_{k}\right)\left(\sum_{h=1}^{J}p_{ih}x_{h}\right)$
$\displaystyle=-\sum_{i=1}^{I}y_{i}^{2}$

where

$y_{i}=\frac{\sqrt{c_{i}}}{\bm{p}_{i}^{\mathsf{T}}\bm{u}}\sum_{h=1}^{J}p_{ih}x_{h}=\sqrt{c_{i}}\frac{\bm{p}_{i}^{\mathsf{T}}\bm{x}}{\bm{p}_{i}^{\mathsf{T}}\bm{u}}.$

It follows that $\bm{x}^{\mathsf{T}}\bm{H}\bm{x}\leq 0$, with equality only if
$\bm{y}=[y_{1}\cdots y_{I}]^{\mathsf{T}}=\bm{0}$. We can rewrite
$\bm{y}=\bm{D}\bm{P}\bm{x}$, where
$\bm{D}=\mathrm{diag}([\frac{\sqrt{c_{1}}}{\bm{p}_{1}^{\mathsf{T}}\bm{u}}\
\cdots\ \frac{\sqrt{c_{I}}}{\bm{p}_{I}^{\mathsf{T}}\bm{u}}])$. Since $\bm{D}$
has full rank $I$, $\bm{D}\bm{P}$ has the same rank of $\bm{P}$; the latter
has rank that cannot exceed $I<J$, hence the rank of $\bm{D}\bm{P}$ is smaller
than $J$. We conclude that there exist some $\bm{x}\neq\bm{0}$ such that
$\bm{y}=\bm{D}\bm{P}\bm{x}=\bm{0}$, hence $\bm{H}$ is negative semidefinite. ∎

The equivalence between MLE-Multinomial and MLE-Poisson (i.e., EM algorithm
(11)) is established as follows.

######  Proposition 3.

The optimality condition of MLE-Multinomial is identical to the convergence
condition of the EM procedure based on MLE-Poisson, i.e., the two methods
share the same solution space.

###### Proof.

We start by considering the iterative EM procedure derived from MLE-Poisson:
at convergence, the stable points of (11) must fulfill the condition
$\hat{\bm{u}}^{m+1}=\hat{\bm{u}}^{m}$ ($m$ denoting the iteration index)
therefore for the solution vector $\hat{\bm{u}}$ it holds:

$\sum_{i=1}^{I}{c_{i}\frac{p_{ij}}{\sum_{k=1}^{J}{p_{ik}\hat{u}_{k}}}}=\sum_{i=1}^{I}c_{i}\frac{p_{ij}}{\bm{p}_{i}^{\mathsf{T}}\hat{\bm{u}}}=1,\quad\forall
j$ (15)

wherein $\bm{p}_{i}^{\mathsf{T}}$ denotes the $i$th row of $\bm{P}$.

For the MLE-Multinomial formula, from Proposition 2 we know that all
stationary points are equivalent global maxima. To identify the stationary
points, we derive the Lagrangian function, i.e., the derivative with respect
to $u_{j}$ of the objective function from (10) augmented by the scaled version
of the equality constraint. By exploiting the fact that the $i$th element of
vector $\bm{P}\bm{u}$ is $\bm{p}_{i}^{\mathsf{T}}\bm{u}$, we obtain, for each
$j$:

$\frac{\partial(\bm{c}^{\mathsf{T}}\log{\bm{P}\bm{u}}-\lambda(\|\bm{u}\|_{1}-C))}{\partial
u_{j}}=\sum_{i=1}^{I}c_{i}\frac{p_{ij}}{\bm{p}_{i}^{\mathsf{T}}\bm{u}}-\lambda=0.$

Multiplying by $u_{j}$ and summing over $j$ we obtain

$\sum_{i=1}^{I}c_{i}\frac{\sum_{j}p_{ij}u_{j}}{\bm{p}_{i}^{\mathsf{T}}\bm{u}}=\lambda
C$

and since $\sum_{j}p_{ij}u_{j}=\bm{p}_{i}^{\mathsf{T}}\bm{u}$ it turns out
that $\lambda=1$. By substituting back in the first-order optimality condition
it follows that, for each $j$,

$\sum_{i=1}^{I}c_{i}\frac{p_{ij}}{\bm{p}_{i}^{\mathsf{T}}\bm{u}}=1.$ (16)

It is evident that the optimality condition (16) derived for MLE-Multinomial
is identical to the convergence condition (15) of the EM procedure derived for
MLE-Poisson. That means, the two methods share the same solution
space999Strictly speaking, this holds true as long as the initialization
vector in the iterative procedure (11) is free from zero elements, as
discussed in Sec. IV-D.. ∎

Looking retrospectively, such equivalence should have been expected in the
light of the close relationship between the Poisson and Multinomial
distributions [25, Ch. V]: if $K$ independent random variables follow a
Poisson distribution $x_{k}\sim\mathcal{P}(\rho_{k})$, $k=1,\ldots,K$, then
the distribution of the vector $\bm{x}\stackrel{{\scriptstyle\textrm{\tiny
def}}}{{=}}[x_{1}\ \cdots\ x_{K}]^{\mathsf{T}}$ _conditioned to the sum_
$\sum_{k=1}^{K}{x_{k}}=\|\bm{x}\|_{1}=X$ is Multinomial
$\mathcal{M}(X,\bm{\mu})$ with
$\mu_{k}=\frac{\rho_{k}}{\sum_{\iota=1}^{K}{\rho_{\iota}}}$.

As an additional result, we can easily recognize that the multiplicative
factor appearing in (11) equals the partial derivative of the objective
function in (10), i.e.

$\frac{\partial}{\partial
u_{j}}\bm{c}^{\mathsf{T}}\log{\bm{P}\bm{u}}=\sum_{i=1}^{I}{c_{i}\frac{p_{ij}}{\sum_{k=1}^{J}{p_{ik}{u}_{k}}}}.$
(17)

In other words, we can reinterpret the EM procedure in (11) (derived for the
Poisson generative model) as a purely multiplicative method to solve
iteratively the optimization (10) (derived for the Multinomial generative
model) starting from an initial guess101010Such purely multiplicative
procedure in the form $u_{j}^{m+1}=u_{j}^{m}\cdot\psi_{j}^{m}$ is distinct
from the standard gradient descent approaches that, in general, involves an
additive update in the form $u_{j}^{m+1}=u_{j}^{m}+\phi_{j}^{m}$..

###  V-B Maximum Likelihood Bounded Subspace (MLBS)

Recall that the vector $\bm{c}$ represents a _single realization_ of the
random vector $\tilde{\bm{c}}$ with unknown mean $\overline{\bm{c}}$. Since no
other measurements are available, in the absence of external information the
vector $\bm{c}$ represents a natural _estimate of the mean value
$\overline{\bm{c}}$_. Replacing the unknown term $\overline{\bm{c}}$ with its
estimate $\bm{c}$ in the generative relation (2) leads to the constraint
$\bm{P}\bm{u}=\bm{c}$. It is rather intuitive that, among all possible (non-
negative) values of $\bm{u}$, those respecting this constraint are the ones
that best conform with the available measurement $\bm{c}$ and with the model
$\bm{P}$. This argument leads us to restrict the search for a “good” estimate
$\hat{\bm{u}}$ within the bounded subspace $\mathcal{U}$ defined formally by:

$\mathcal{U}\stackrel{{\scriptstyle\textrm{\tiny
def}}}{{=}}\left\\{\bm{u}\;:\;\bm{P}\bm{u}=\bm{c},\,\bm{u}\geq\bm{0}\right\\}.$
(18)

Note that the constraint $\bm{P}\bm{u}=\bm{c}$ absorbs the condition
$\|\bm{u}\|_{1}=C$ on the total count. Since in our application $J\gg I$,
matrix $\bm{P}$ is rank deficient and the constraint $\bm{P}\bm{u}=\bm{c}$
admits multiple solutions in the variable $\bm{c}$. Therefore, while we have
restricted somewhat the search space, we still need to provide a criterion for
selecting unambiguously a single solution within that space. This aspect is
elaborated later in Sec. VII.

While this way of reasoning might seem heuristic, it turns out that the
bounded subspace defined in (18) is intimately connected to the MLE-
Multinomial and MLE-Poisson procedures described earlier.

######  Proposition 4.

The bounded subspace $\mathcal{U}$ represents a set of solutions for both MLE
procedures, i.e., any point in $\mathcal{U}$ is also a solution of MLE-
Multinomial and MLE-Poisson.

###### Proof.

The condition $\bm{P}\bm{u}=\bm{c}$ can be rewritten as
$\bm{p}_{i}^{\mathsf{T}}\bm{u}=c_{i}$, $\forall i$, and that implies that the
optimality condition (16), or equivalently (15), is always verified: with the
ansatz $\bm{p}_{i}^{\mathsf{T}}\bm{u}=c_{i}$, $\forall i$, (16) reduces to the
identity $\sum_{i=1}^{I}p_{ij}=1$ that is always satisfied by the column-
stochasticity of $\bm{P}$, i.e., $\bm{1}_{I}^{\mathsf{T}}\bm{P}=\bm{1}_{J}$.
Therefore, we have established that $\mathcal{U}$ represents a set of
solutions for both MLE procedures, i.e., any point in $\mathcal{U}$ is also a
solution of MLE-Multinomial and MLE-Poisson. ∎

Following Proposition 4, we shall refer to $\mathcal{U}$ as the Maximum
Likelihood Bounded Subspace (MLBS).

######  Remark 2.

The SB in (14) may fall outside the MLBS (a numerical example is given below
in the next subsection) therefore does not qualify in general as a ML
solution. Moreover, the SB solution coincides with the first point of the EM
sequence after the first iterative step when the starting point is set to the
prior vector $\bm{\alpha}$, as can be immediately verified by comparing (14)
and (11), the only difference being that $\hat{\bm{u}}^{m}_{j}$ in (11) is
replaced by $\bm{\alpha}_{j}$ in (14). In other words, the SB point lies
between the prior vector $\bm{\alpha}$ and the MLBS.

###  V-C Graphical illustration for a toy scenario

To illustrate graphically the above results, we present numerical results
referred to a simple toy scenario with $J=3$ tiles and $I=2$ cells in the
configuration sketched in Fig. 2, with observation count vector
$\bm{c}=[40\;70]^{\mathsf{T}}$, $C=110$ and model matrix
$\bm{P}=\begin{bmatrix}1&.25&0\\\ 0&.75&1\\\ \end{bmatrix}$. In this toy
instance there are only two degrees of freedom, as the $J=3$ variables are
constrained by the total sum $\sum_{j=1}^{J}{u_{j}}=C$, or equivalently
$\|\bm{u}\|_{1}=\|\bm{c}\|_{1}$. Therefore, we can map all admissible
solutions $\bm{u}$ to an horizontal plane. In the 3D plot of Fig. 3 the
horizontal dimensions represent $u_{1}$ and $u_{3}$, while the vertical
dimension reports the value of the likelihood function (7). The convexity of
the likelihood function (proven in Proposition 2) is graphically evident.
Notably, the locus of all global maxima of the maximum function lie on a
segment, that is the 1D MLBS of the whole 2D solution space. In the planar
plot of Fig. 4 we mark with blue dots the MLBS. The flat initial guess (for
which $u_{1}=u_{2}=u_{3}=36.7$), the corresponding SB estimate from (14) and
the EM solution from (11) are indicated by markers.

Figure 2: Toy example with $J=3$ tiles and $I=2$ cells. Figure 3: Log-
likelihood surface for the toy example of Fig. 2. Figure 4: Solution points
for the toy example of Fig. 2.

##  VI The Maximum A Posteriori (MAP) estimator

In this section we derive the exact and approximated forms of the Maximum A
Posteriori (MAP) estimator for the problem at hand. To the best of our
knowledge, this is a novel contribution not presented earlier in previous
literature. The approximated MAP form serves as basis for the formulation of a
novel type of estimators elaborated in the following section.

###  VI-A Derivation

We derive the MAP estimator based on the hierarchical generative model with
Multinomial distribution shown in Fig. 1(a). The result is summarized in the
following Proposition.

######  Proposition 5.

The MAP estimator of $\bm{u}$ can be obtained as

$\begin{split}\hat{\bm{u}}_{\text{MAP}}&=\operatorname*{arg\,max}_{\begin{subarray}{c}\|\bm{u}\|_{1}=C\\\
\bm{u}\geq\bm{0}\end{subarray}}{\big{\\{}\bm{c}^{\mathsf{T}}\log{\bm{P}\bm{u}}+\bm{u}^{\mathsf{T}}\log{\bm{\alpha}}-\sum_{j=1}^{J}{\log{u_{j}!}}\big{\\}}}\end{split}$
(19)

###### Proof.

Following the stochastic model shown in Fig. 1(a), $\bm{u}$ represents a
random vector with Multinomial distribution of unknown parameter $\bm{\mu}$.
Taking the _prior vector_ $\bm{\alpha}$ as the prior value (before seeing the
data) of $\bm{\mu}$, we obtain that the _prior probability distribution_ of
$\bm{u}$, denoted by $\mathbb{P}(\bm{u})$, is Multinomial with parameters $C$
and $\bm{\alpha}$, formally

$\mathbb{P}(\bm{u})=\mathcal{M}\left(C,\bm{\alpha}\right)=\frac{C!}{u_{1}!u_{2}!\cdots
u_{J}!}\prod_{j=1}^{J}{\left(\bm{\alpha}_{j}\right)^{u_{j}}}$ (20)

from where the logarithm

$\log{\mathbb{P}(\bm{u})}=\bm{u}^{\mathsf{T}}\log{\bm{\alpha}}-\sum_{j=1}^{J}{\log{u_{j}!}}$

is obtained, up to an irrelevant (constant) additive term.

The posterior distribution is obtained by the product of likelihood and prior
distribution $\mathbb{P}(\bm{u})$. Accordingly, the logarithm of the posterior
distribution equals the sum of the log-likelihood (7) and
$\log\mathbb{P}(\bm{u})$. Therefore, maximizing the (logarithm of) posterior
distribution leads to maximizing the following function of $\bm{\mu}$ and
$\bm{u}$

$\operatorname*{arg\,max}_{\begin{subarray}{c}\|\bm{\mu}\|_{1}=1,\;\|\bm{u}\|_{1}=C\\\
\bm{\mu}\geq\bm{0},\;\bm{u}\geq\bm{0}\end{subarray}}{\big{\\{}\bm{c}^{\mathsf{T}}\log{\bm{P}\bm{\mu}}+\log\mathbb{P}(\bm{u})\big{\\}}}.$
(21)

Recalling (9) we bind the estimate $\hat{\bm{u}}$ to equal the rescaled value
of $\hat{\bm{\mu}}$, i.e., we impose $\hat{\bm{u}}=C\hat{\bm{\mu}}$. With this
simple variable substitution, from (21) we obtain the following MAP
estimator111111In the derivation we have used the following invariance
property: a constant factor $\kappa$ in the argument of the logarithm terms
does not influence the solution to the maximization. Formally:
$\operatorname*{arg\,max}_{\bm{x},\bm{y}}\\{\log{\frac{\bm{x}}{\kappa}}+\bm{y}\\}=\operatorname*{arg\,max}_{\bm{x},\bm{y}}\\{\log{\bm{x}}-\log{\kappa}+\bm{y}\\}=\operatorname*{arg\,max}_{\bm{x},\bm{y}}\\{\log{\bm{x}}+\bm{y}\\}.$:

$\begin{split}\hat{\bm{u}}_{\text{MAP}}=\operatorname*{arg\,max}_{\begin{subarray}{c}\|\bm{u}\|_{1}=C\\\
\bm{u}\geq\bm{0}\end{subarray}}{\big{\\{}\bm{c}^{\mathsf{T}}\log{\bm{P}\bm{u}}+\log\mathbb{P}(\bm{u})\big{\\}}}\end{split}$
(22)

which can be rewritten to finally yield the thesis (note that the term
$\bm{\alpha}$ can be replaced by $\bm{a}=\bm{\alpha}/C$ for the same
invariance property recalled in the previous footnote). ∎

In eq. (19) the term $\bm{c}^{\mathsf{T}}\log{\bm{P}\bm{u}}$ carries the data
information (from the measurements $\bm{c}$) while the other two terms carry
information from the _prior distribution_ $\mathbb{P}(\bm{u})$. More
specifically, the term $\bm{u}^{\mathsf{T}}\log{\bm{a}}$ carries information
from the _prior vector_ 121212Recall that in our model the prior vector
represents the parameter of the prior distribution. $\bm{a}$, or equivalently
$\bm{\alpha}$, while the term $\sum_{j=1}^{J}{\log{u_{j}!}}$ captures the
_combinatorial diffusion_ effect that is intrinsic to the adoption of a
Multinomial distribution, and induces a preference for solutions with higher
entropy (it can be easily proved that this terms is maximized when all
elements of $\bm{u}$ are equal).

######  Remark 3.

The discrete factorial term appearing in (19) could be replaced by its
analytical continuation, i.e., by the Gamma function $\Gamma(n+1)=n!$, leading
to an equivalent continuous function that, in principle, can be minimized
numerically via standard numerical methods (e.g., gradient descent). However,
for very large problem instances, with $I$ and $J$ in the range of tens of
thousands, numerical resolution with general purpose solvers might still be
too impractical, motivating the derivation of a computationally simpler
alternative.

###  VI-B Approximated MAP

In Sec. VII we will present a novel alternative estimator, labelled DF for
“‘data first”, that exploits the particular structure of the problem at hand
and yields a closed-form analytic solution. In order to provide additional
insight, we will also elaborate on the relation and conceptual difference
between the new DF estimator and the classical MAP estimator derived above. To
this aim, it is preliminarily convenient to derive an approximated version of
the MAP estimator, as follows.

######  Proposition 6.

The MAP estimator (19) can be approximated by solving the following simpler
optimization problem:

$\hat{\bm{u}}_{\text{A.M.}}=\operatorname*{arg\,max}_{\begin{subarray}{c}\|\bm{u}\|_{1}=C\\\
\bm{u}\geq\bm{0}\end{subarray}}{\big{\\{}\bm{c}^{\mathsf{T}}\log{\bm{P}\bm{u}}-\frac{1}{2}\|\bm{u}-\bm{a}\|^{2}_{\bm{a}}\big{\\}}}$
(23)

where $\|\bm{u}-\bm{a}\|^{2}_{\bm{a}}\stackrel{{\scriptstyle\textrm{\tiny
def}}}{{=}}(\bm{u}-\bm{a})^{\mathsf{T}}\bm{A}^{-1}(\bm{u}-\bm{a})=\sum_{j=1}^{J}\frac{(u_{j}-a_{j})^{2}}{a_{j}}$,
with $\bm{A}=\mathrm{diag}(\bm{a})$, is the _weighted $\ell^{2}$-norm_ of the
difference vector $(\bm{u}-\bm{a})$.

###### Proof.

We start by approximating the (prior) Multinomial distribution with a
multivariate normal

$\mathbb{P}(\bm{u})=\mathcal{M}(C,\bm{\alpha})\simeq\mathcal{N}(\bm{a},\bm{\Sigma}_{\bm{a}})$
(24)

with mean value $\overline{\bm{u}}=\bm{a}=C\bm{\alpha}$ and covariance matrix
$\Sigma_{\bm{a}}$ composed of the following elements:

$\sigma_{mn}^{2}=\left\\{\begin{array}[]{rcl}C\,a_{m}\left(1-a_{m}\right)&\mbox{for}&m=n\\\
C\,a_{m}\,a_{n}&\mbox{for}&m\neq n\\\ \end{array}\right..$ (25)

Plugging the normal approximation (24) into (22) leads to the approximate MAP
estimator

$\hat{\bm{u}}_{\text{A.M.}}=\operatorname*{arg\,max}_{\begin{subarray}{c}\|\bm{u}\|_{1}=C\\\
\bm{u}\geq\bm{0}\end{subarray}}{\big{\\{}\bm{c}^{\mathsf{T}}\log{\bm{P}\bm{u}}-\frac{1}{2}(\bm{u}-\bm{a})^{\mathsf{T}}\bm{\Sigma}_{\bm{a}}^{-1}(\bm{u}-\bm{a})\big{\\}}}.$
(26)

In our application we are considering very large number of tiles ($J\gg 1$),
which implies that the individual probabilities are very small in absolute
terms ($\alpha_{j}\ll 1$). This justifies the further simplification of the
covariance matrix: neglecting the product terms in the off-diagonal elements,
$C\alpha_{m}\,\alpha_{n}\simeq 0$, and approximating the diagonal elements
with $C\,\alpha_{m}\left(1-\alpha_{m}\right)\simeq C\,\alpha_{m}=a_{m}$, the
covariance matrix reduces to the diagonal matrix
$\bm{A}\stackrel{{\scriptstyle\textrm{\tiny def}}}{{=}}\mathrm{diag}([a_{1}\
\cdots\ a_{J}])$ and its inverse to
$\bm{\Sigma}_{\bm{a}}^{-1}\simeq\bm{A}^{-1}=\mathrm{diag}([\frac{1}{a_{1}}\
\ldots\ \frac{1}{a_{J}}])$. In this way the approximate MAP estimator rewrites

$\hat{\bm{u}}_{\text{A.M.}}\approx\operatorname*{arg\,max}_{\begin{subarray}{c}\|\bm{u}\|_{1}=C\\\
\bm{u}\geq\bm{0}\end{subarray}}{\big{\\{}\bm{c}^{\mathsf{T}}\log{\bm{P}\bm{u}}-\frac{1}{2}(\bm{u}-\bm{a})^{\mathsf{T}}\bm{A}^{-1}(\bm{u}-\bm{a})\big{\\}}}$
(27)

which finally yields the thesis. ∎

##  VII A novel class of prior-friendly MLE

###  VII-A Motivations

In Sec. V we have established that all points within MLBS attain the same
maximum value of the likelihood function. Therefore, the estimation problem is
structurally non-identifiable, following the definition given in [22], in the
sense that it is not possible to pick a unique solution based solely on the
information contained _in the data_ $\bm{c}$. As all points within the MLBS
conform to the data equally well. in order to select a particular point within
the MLBS, or at least restrict to a subset of preferred solutions, we must
necessarily resort to some additional assumption and/or to external
information.

One possibility is to demand that, in addition to conforming to the measured
data, the solution shall conform also to some particular structural property
characterizing the physical process we aim to measure, e.g., in the form of
some smoothness criterion, minimum gradient between adjacent tiles, etc. A
second approach is to resort to some kind of _prior information_ to aid the
solution selection process. Both approaches amount to adding an additional
component to the objective function to be minimized, and therefore may be
interpreted as different forms of regularization. For our specific
application, dealing with human distribution in space, there are no obvious
structural constraints tied to the “physics” of the underlying process, and it
appears more natural to encode external information in the form of prior
information (as done, e.g., in [23]).

In principle, a possible way to account for prior information (however
determined) is through the MAP estimator derived in the previous section.
Generally speaking, the MAP solution combines information _from the data_ ,
captured by the likelihood term embedding the measurement vector $\bm{c}$,
with information _before the data_ , captured by the prior distribution
embedding the _prior vector_ $\bm{\alpha}$ (or equivalently $\bm{a}$). In
fact, the posterior distribution is composed by the product of the two
components: $\text{posterior}=\text{likelihood}\times\text{prior}$. The
resulting solution is a trade-off, i.e., a compromise between the two sources
of information. Likewise for other Bayesian estimators, as more data are
available (more samples) the solution component driven by the data information
increases its relative importance and eventually dominates over the prior
information. On the opposite direction, when the data are scarce, the prior
information may be dominating. In the extreme situation that _only a single
sample measurement is available_ , as is specifically the case at hand in our
application, adopting a MAP estimator involves a certain risk of diminishing
the relative weight of the measurement information to the point that it almost
vanishes. In other words, as with MAP (exact or approximated) we cannot
control explicitly how much weight to put on the data vis-à-vis the prior,
there is a certain risk of ending up with a solution that reflects mostly the
prior information, only slightly perturbed by the measurements.

To avoid this undesirable effect, in the following we develop an alternative
estimation procedure for the problem at hand based on the “data first”
principle, where in-data information (likelihood) is given priority over other
off-data knowledge (prior distribution and/or any other structural property).

###  VII-B A “data first” approach

The goal expressed above is achieved though an estimator built according to
the following structure:

$\hat{\bm{u}}=\operatorname*{arg\,min}_{\bm{u}\in\mathcal{U}}{f\big{(}\bm{u},\bm{a}\big{)}}$
(28)

with $f(\bm{u},\bm{a})$ denoting a distance function between vectors $\bm{u}$
and $\bm{a}$. With this structure, _the ML property is imposed as hard
constraint_. In fact, the condition

$\bm{u}\in\mathcal{U}\Leftrightarrow\big{\\{}\bm{P}\bm{u}=\bm{c},\;\bm{u}\geq\bm{0}\big{\\}}$
(29)

forces the solution to lie within the MLBS that is by definition the locus of
all ML points, i.e., the points that _best conform to the measured data_
$\bm{c}$. Among these points, then we select the point that additionally _best
conforms with the off-data knowledge_ encoded in the prior vector $\bm{a}$,
through the minimization of some distance function
$f\big{(}\bm{u},\bm{a}\big{)}$. Among many possible choices for the distance
function, opting for the weighted $\ell^{2}$-norm defined in Proposition 6:

$f\big{(}\bm{u},\bm{a}\big{)}=\|\bm{u}-\bm{a}\|^{2}_{\bm{a}}\stackrel{{\scriptstyle\textrm{\tiny
def}}}{{=}}(\bm{u}-\bm{a})^{\mathsf{T}}\bm{A}^{-1}(\bm{u}-\bm{a})$ (30)

with $\bm{A}\stackrel{{\scriptstyle\textrm{\tiny
def}}}{{=}}\mathrm{diag}(\bm{a})$, allows us to establish a direct connection
between the newly proposed estimator and the approximate MAP form derived in
(23), as elaborated below. Plugging (29) and (30) into the structure leads to
the following novel estimator (labelled “DF” for “Data First”):

$\hat{\bm{u}}_{\text{DF}}=\operatorname*{arg\,min}_{\begin{subarray}{c}\bm{P}\bm{u}=\bm{c}\\\
\bm{u}\geq\bm{0}\end{subarray}}{(\bm{u}\\!-\\!\bm{a})^{\mathsf{T}}\bm{A}^{-1}(\bm{u}\\!-\\!\bm{a})}.$
(31)

Remarkably, the DF solution (31) can be given in closed-form, as summarized in
the following result.

######  Proposition 7.

The DF estimator (31) is given by

$\hat{\bm{u}}_{\text{DF}}=\max(\bm{A}\bm{P}^{\mathsf{T}}(\bm{P}\bm{A}\bm{P}^{\mathsf{T}})^{-1}(\bm{c}-\bm{P}\bm{a})+\bm{a},\bm{0})$
(32)

where the maximum function is intended element-wise.

###### Proof.

The constrained minimization problem can be solved by considering the Karush-
Kuhn-Tucker conditions:

$\left\\{\begin{array}[]{l}\nabla_{\bm{u}}\left[\frac{1}{2}(\bm{u}\\!-\\!\bm{a})^{\mathsf{T}}\bm{A}^{-1}(\bm{u}\\!-\\!\bm{a})+\bm{\lambda}^{\mathsf{T}}(\bm{P}\bm{u}-\bm{c})-\bm{\nu}^{\mathsf{T}}\bm{u}\right]=\bm{0}\\\
\bm{P}\bm{u}=\bm{c}\\\ \nu_{j}\geq 0,\quad\nu_{j}u_{j}\geq 0,\quad u_{j}\geq
0,\quad j=1,\ldots,J\end{array}\right.$

The calculation of the gradient in the first line yields

$\bm{A}^{-1}(\bm{u}\\!-\\!\bm{a})+\bm{P}^{\mathsf{T}}\bm{\lambda}-\bm{\nu}=\bm{0}$

from which

$\bm{u}=\bm{A}\bm{\nu}-\bm{A}\bm{P}^{\mathsf{T}}\bm{\lambda}+\bm{a}.$

The vector of the Lagrange multipliers $\bm{\lambda}$ corresponding to the
equality constraint can be explicitly obtained by exploiting the constraint
itself, that is

$\bm{P}\bm{u}=\bm{c}\Leftrightarrow\bm{P}\bm{A}\bm{\nu}-\bm{P}\bm{A}\bm{P}^{\mathsf{T}}\bm{\lambda}+\bm{P}\bm{a}=\bm{c}$

which leads to

$\bm{\lambda}=(\bm{P}\bm{A}\bm{P}^{\mathsf{T}})^{-1}(\bm{P}\bm{A}\bm{\nu}+\bm{P}\bm{a}-\bm{c}).$

Substituting back one obtains

$\displaystyle\hat{\bm{u}}$
$\displaystyle=\bm{A}\bm{\nu}-\bm{A}\bm{P}^{\mathsf{T}}(\bm{P}\bm{A}\bm{P}^{\mathsf{T}})^{-1}(\bm{P}\bm{A}\bm{\nu}+\bm{P}\bm{a}-\bm{c})+\bm{a}$
$\displaystyle=\bm{A}^{\frac{1}{2}}\left[\bm{\mathfrak{P}}_{\bm{A}^{\frac{1}{2}}\bm{P}^{\mathsf{T}}}(\bm{c}-\bm{P}\bm{a})-\bm{\mathfrak{P}}_{\bm{A}^{\frac{1}{2}}\bm{P}^{\mathsf{T}}}^{\perp}\bm{A}^{-\frac{1}{2}}\bm{\nu}\right]+\bm{a}$

where
$\bm{\mathfrak{P}}_{\bm{M}}=\bm{M}(\bm{M}^{\mathsf{T}}\bm{M})^{-1}\bm{M}^{\mathsf{T}}$
is the projection matrix onto the subspace spanned by the columns of the
matrix $\bm{M}$, and
$\bm{\mathfrak{P}}_{\bm{M}}^{\perp}=\bm{I}-\bm{\mathfrak{P}}_{\bm{M}}$ the
projection matrix onto its orthogonal complement. The solution is thus a
shifted version (by $\bm{a}$) of the sum of two vectors, obtained by
projecting the vectors $(\bm{c}-\bm{P}\bm{a})$ and
$\bm{A}^{-\frac{1}{2}}\bm{\nu}$ onto the column-space spanned by
$\bm{A}^{\frac{1}{2}}\bm{P}^{\mathsf{T}}$ and its orthogonal complement,
respectively. Due to the constraints $u_{j}\geq 0$, $\nu_{j}\geq 0$ and
$\nu_{j}u_{j}=0$, where $\nu_{j}=0$ the solution is given by the corresponding
component in the first term only; conversely, where $\nu_{j}>0$ then
necessarily $u_{j}=0$. The overall solution can be finally recast into eq.
(32). ∎

######  Remark 4.

Comparing (31) with (23), we can immediately recognize the key difference
between the DF estimator and the approximate MAP estimator. The latter
minimizes the sum of the data and prior components in the objective function,
encoded respectively in the term $\bm{c}^{\mathsf{T}}\log{\bm{P}\bm{u}}$ and
$(\bm{u}\\!-\\!\bm{a})^{\mathsf{T}}\bm{A}^{-1}(\bm{u}\\!-\\!\bm{a})$, and in
this way it strikes a compromise between in-data and off-data information seen
as peer conflicting goals. Instead, the “data first” DF estimator (31), easily
computable via eq. (32) from Proposition 7, gives absolute priority to the in-
data information, encoded in the constraint $\bm{P}\bm{u}=\bm{c}$ that ensures
maximal adherence to the data (ML condition), and secondarily seeks to conform
to the off-data knowledge encoded in the objective function.

###  VII-C Advantages of the closed-form solution

The availability of a closed-form solution for DF is not necessarily an
advantage in terms of computational efficiency vis-à-vis the EM/ML procedure,
since the analytic form (32) still involves the computation of a matrix
inverse that, for very large instances, might be computationally heavier than
running a handful of iterations of the form (11). However, a careful look at
(32) reveals that the heaviest computation part depends only on the model
matrix $\bm{P}$ and on the prior $\bm{a}$, and is independent from the data
$\bm{c}$. In fact, setting $\bm{F}\stackrel{{\scriptstyle\textrm{\tiny
def}}}{{=}}\bm{A}\bm{P}^{\mathsf{T}}(\bm{P}\bm{A}\bm{P}^{\mathsf{T}})^{-1}$
and $\bm{g}\stackrel{{\scriptstyle\textrm{\tiny
def}}}{{=}}(\bm{F}\bm{P}-\bm{I}_{J})\bm{a}$, eq. (32) can be rewritten as

$\hat{\bm{u}}_{\text{DF}}=\max(\bm{F}\bm{c}-\bm{g},\bm{0})$ (33)

where the maximum function is intended element-wise. While matrix $\bm{F}$ and
vector $\bm{g}$ are heavy to compute (especially because $\bm{F}$ embeds a
matrix inversion) they are both independent from the measured data $\bm{c}$.
Therefore, in scenarios where multiple estimates must be computed for the same
model matrix $\bm{P}$ but with different measurement vectors, $\bm{F}$ and
$\bm{g}$ can be pre-computed only once. Along the same reasoning, the closed-
form (33) is also appealing for designing linear filters where the spatial
density estimate is continuously updated as new measurements become available
from the same network. Finally, the form (33) is also convenient for
investigating, through differentiation, local sensitivity to perturbation of
the measurement input vector $\bm{c}$.

##  VIII Numerical results

In this section we present numerical results comparing the performance of the
different estimators for a sample synthetic scenario. Testing on synthetic
data has two important benefits: _(i)_ it allows one to control explicitly the
data generation parameters, thus enabling full control over the degree of
(mis-)matching between the data and estimation model, and _(ii)_ it allows one
to quantify the absolute estimation error against a known “ground truth”. Both
aspects are important for our study that, as stated in the introduction, is
focused at this stage on investigating alternative solutions for a specific
block in a larger methodological chain. For the sake of completeness we report
also estimates based on simple Voronoi tessellation based on cell towers, that
is the most popular approach in the literature.

###  VIII-A Scenario

We briefly describe the data generation scenario. The reference area consists
of a regular square grid of $100\times 100$ tiles (the same size considered
earlier in [9, 10]). The radio network topology consists of two layers, with a
total of 115 tower locations placed as shown in Fig. 5(a). The first radio
layer consists of 42 antenna towers (red triangles in Fig. 5(a)) with each
tower hosting three $120^{\circ}$ sectors with a maximum range (normalized to
the tile size) equal to 15 and azimuth equal to $\pi/5$. The second layer
consists of 73 antenna towers (blue circles in Fig. 5(a)) with each tower
hosting three $120^{\circ}$ sectors with a maximum range equal to 6 and
$0^{\circ}$ azimuth. The cell coverage is approximated to the tile level: a
tile is considered covered by a radio cell if the tile center falls within the
radio cell coverage perimeter.

Taking the 115 tower locations as seeds for a Voronoi tessellation results in
the pattern shown in Fig. 5(b), where each colour marks a distinct Voronoi
polygon.

For the cell overlapping pattern (with flat footprints), the consolidation
process results in the set of _sections_ (or equivalently _super-tiles_) shown
in Fig. 5(c). This is also a partition (by definition, each section gathers
all tiles covered by the same set of radio cells) but much finer than Voronoi,
with a total of 1654 sections, i.e., $\times 14$ more than Voronoi polygons in
this specific scenario.

The actual “ground truth population” (GTP) of mobile devices is shown in Fig.
6(a): it consists of four squared clusters representing densely populated
areas, and a clutter of moderately populated tiles sparse over the region of
interest.

For a given partition (Voronoi polygons or sections) averaging the GTP within
each (disjoint) partition area represents the best possible “oracle estimate”
for that specific partition. The results for the Voronoi and section
partitions are shown in Figs. 6(b) and 6(c), respectively. As expected, the
finer texture of the latter renders a spatial distribution that is closer to
the original GTP, and such qualitative visual impression will be confirmed
quantitatively by the numerical results presented later.

As for the assignment of population units to radio cell, we follow the same
probabilistic generative model considered earlier in [9, 10]: a generic mobile
device placed selects randomly the radio cell among those in range, with equal
probabilities and independently from other devices. Therefore, if tile $j$ is
covered by radio cell $i$ plus other $k_{j}-1$ radio cells, then a generic
mobile device in this tile will select cell $i$ with probability
$p_{ij}=\frac{1}{k_{j}}$. The resulting matrix $\bm{P}$ has non-negative
fractional values. The same matrix $\bm{P}$ is used in the data generation
process and in the estimation models, i.e., we test the different methods in
model matching condition.

(a) Tower locations

(b) Voronoi polygons

(c) Sections (aka super-tiles)

Figure 5: Network topology.

(a) Ground truth population (GTP)

(b) GTP mapped to Voronoi

(c) GTP mapped to sections (super-tiles)

Figure 6: Spatial distribution of actual population per grid unit (a), per
Voronoi polygon (b) and per section (c).

(a) Voronoi estimate

(b) SB estimate

(c) DF estimate

Figure 7: Estimated density maps (with flat prior). The ML/EM estimate is
omitted as it is visually indistinguishable from DF.  Figure 8: Kantorovich-
Wasserstein Distances (KWD) between GTP and different spatial distributions.

###  VIII-B Similarity metric

In order to evaluate the goodness of a generic estimation method X we must
resort to some similarity metric
$w(\hat{\bm{u}}_{\text{X}},\bm{u}_{\text{GTP}})$ between the estimated density
$\hat{\bm{u}}_{\text{X}}$ and the ground truth population
$\bm{u}_{\text{GTP}}$. To this aim, we resort to the Kantorovich-Wasserstein
Distance (KWD for short) with Euclidean ground distance. KWD is known in the
literature with different names, e.g., Earth Mover Distance, Mallow distance,
etc. (see [26] and references therein).

For two input distributions $\bm{u}$ and $\bm{v}$ defined over the same grid
$\mathcal{B}$ and with same total mass $\|\bm{u}\|_{1}=\|\bm{v}\|_{1}=C$, KWD
may be interpreted as the minimum cost of transporting the mass from
configuration $\bm{u}$ to $\bm{v}$ (or vice-versa) when the cost of
transporting a unit of mass between two generic tiles $j$ and $k$ equals the
ground distance $d_{kj}$ between the tile centres. KWD can be expressed as the
solution $w$ of the following Linear Programming (LP) optimization problem:

$\begin{array}[]{ll@{\quad}ll}\text{minimize}&w(\bm{u},\bm{v})\stackrel{{\scriptstyle\textrm{\tiny
def}}}{{=}}\frac{1}{C}\sum_{k}{\sum_{j}{x_{kj}d_{kj}}}&&\\\ \text{subject
to}&\sum_{j}{x_{kj}}=u_{k}&\forall\,k\in\mathcal{B}\\\
&\sum_{k}{x_{kj}}=v_{j}&\forall\,j\in\mathcal{B}\\\ &x_{kj}\geq
0&\forall\,k,j\in\mathcal{B}\end{array}$ (34)

The direct resolution of such LP problem is computationally expensive,
preventing the computation of exact KWD values on very large grids. In a
recent work Bassetti, Gualandi and Veneroni [26] have shown that a close
approximation to the exact KWD value, within a provable deterministic bound,
can be obtained by solving a transportation problem over a regular lattice. In
the proposed solution, the integer parameter $L$ determines the density of the
lattice and can be used as a tuning knob to trade-off computation time (and
memory) with approximation error: they show that $L=4$ is sufficient to
achieve a KWD approximation within 1% of the exact KWD value. They also
release an open-source implementation of their code131313Available at
https://github.com/stegua which we used to compute the KWD (with Euclidean
ground distance and lattice parameter $L=4$) between each estimate and GTP.

A physical interpretation of KWD
$w(\hat{\bm{u}}_{\text{X}},\bm{u}_{\text{GTP}})$ in our context is the
following. Imagine that, starting from the estimated density
$\hat{\bm{u}}_{\text{X}}$, we move each population unit (i.e., a mobile
device) across tiles so as to arrive to the final GTP configuration
$\bm{u}_{\text{GTP}}$ (transportation plan), and we do so in the way that
minimizes the total travelled distance (_optimal_ transportation plan). Some
population units will not travel, while others will travel over shorter or
longer distances. _On average, every population unit will travel a distance
equal to the KWD value $w$_. In other words, since KWD is normalized to the
total mass (note the factor $\frac{1}{C}$ in the defining LP formulation (34))
its value can be interpreted as the _average distance_ in the optimal
transportation plan between GTP and the estimated map (or vice-versa).

###  VIII-C Estimation accuracy assessment

In Fig. 7 we report the density maps obtained with the different estimation
methods. When needed, the prior vector or initial point was set to the (non-
informative) uniform vector with equal elements. For the sake of completeness,
we include also the map obtained with the simple Voronoi method with seeds on
tower locations. The map obtained by the ML/EM estimator (11) is visually
indistinguishable from the DF map obtained with (32) and is therefore omitted
(note however that the ML/EM solution is not exactly equal to the DF
solution). For each estimated map, the KWD values (from GTP) are plotted in
Fig. 8 along with additional reference values. The KWD value between GTP and
the flat distribution serves as a naif reference. The KWD values between GTP
and the oracle estimates (refer to Fig. 6(b) for Voronoi and Fig. 6(c) for
sections) represent best-case references for their corresponding methods.

Based on the numerical results reported in Fig. 8 we can summarize the
following findings141414Similar plots and numerical results were obtained for
other choices of the population and/or network deployment scenario obtained
with different values for the generative parameters, and in this sense the
reported qualitative findings are representative of a wider exploration of
several test scenarios. . First, as expected, all estimation methods improve
over the naif reference: the average travel distance decreases from 12 in the
flat map down to 4 in the Voronoi estimate and below 1 in the ML/EM and DF
methods.

Second, as expected, all methods based on overlapping cells outperform the
Voronoi based solution. Even the simple SB method (14) improves by a factor of
$\times 2$ over Voronoi (from 4 to 1.9) and with ML/EM the gain factor
increases to $\times 4$ (from 4 to 1).

Third, we note that the absolute gap between the two oracle estimates
(labelled in Fig. 8 as “Or. Vor.” for Voronoi, and “Or. Sec.” for the section
partition resulting from the overlapping cells method) is relatively small:
that means the poorer performance of Voronoi can be accounted only partially
to the coarser spatial resolution (Voronoi polygons are fewer and larger than
sections, compare Figs. 5(b) and 5(c)) and should be rather seen as a
consequence of assuming that mobile devices always connect to the nearest cell
tower. This assumption is clearly violated in the considered simulation
scenario, but it is also highly problematic in real-world scenarios when one
considers the typical features of cellular radio planning (multi-layer radio
coverage design, mixing of small/large cells, power control, etc.).

Fourth, within the scope of estimation methods for overlapping cells, we find
that the ML/EM method outperforms the SB method by a factor of $\times 2$ (see
inset in Fig. 8). We highlight that such improvement comes with a relatively
small additional computation cost, since the EM procedure (11) is iterative
and converges quickly in a few iterations. We remark that the EM solution is
analytically equivalent to the MLE solution developed independently in [10],
but the resolution procedure given in (11) is much simpler to implement and
faster to compute than resorting to general-purpose solvers as done in [10].

Fifth, we find that the proposed DF method improves further over ML/EM, even
though the gain is small (in the order of 5%). This indicates that the EM
iteration tends to reach a final point in the MLBS that is oriented in the
direction of the starting point $\bm{a}$, and therefore closer than most other
points in the MLBS, though not necessarily the closest one. In other words,
the “implicit” tie to the starting point $\bm{a}$ that is intrinsic to the EM
procedure (11) can be improved, though marginally, by the “explicit”
optimization of the weighted distance in (31).

Finally, we find that the accuracy of DF is relatively close to corresponding
oracle solution (labelled “Or. Sec.” in Fig. 8) indicating that little room is
left for further improvement with alternative estimators derived for the same
model.

##  IX Conclusions and future work

The analytic and numerical results presented in this study indicate that, when
the model matrix $\bm{P}$ is known accurately, the best estimates are obtained
with EM/ML and DF. Compared to EM/ML, the DF solution brings a (slight)
accuracy gain and, importantly, yields a closed-form analytic solution whose
particular structure is advantageous for recurrent estimation tasks (e.g.,
with same model matrix but different measurement data, as hinted in Sec.
VII-C).

Our results contribute to advance our understanding of the relative advantages
and costs of the (recently proposed) geo-location methods based on overlapping
cells vis-à-vis the (more traditional) tessellation methods. At the end of the
day, in all practical applications the choice between one approach or the
other is a matter of balancing benefits and costs. On the benefit side, we
have established that EM/ML and DF estimation methods bring a potential
accuracy gain to the overlapping cells approach under matched model condition,
i.e., in the ideal case that the assignment probabilities $p_{ij}$’s, or
equivalently the radio cell footprints $s_{ij}$’s, can be determined
precisely. Compared with the simpler (but less accurate) SB method, the
additional _computation cost_ is contained (especially with the iterative
EM/ML procedure) and therefore should not play a discriminant role in the
choice of the estimator.

The natural next step in the research is to investigate the _sensitivity_ of
the final solution to model mismatching error, i.e., the impact of uncertainty
in determining the model parameters $p_{ij}$’s (or equivalently $s_{ij}$’s) on
the final spatial accuracy.

## References

  * [1] D. Naboulsi, M. Fiore, S. Ribot, and R. Stanica.  Large-scale mobile traffic analysis: A survey.  IEEE Communications Surveys and Tutorials, 18(1), 2016.  DOI 10.1109/COMST.2015.2491361.
  * [2] C. Santamaria Serna, F. Sermi, S. Spyratos, S. Iacus, A. Annunziato, D. Tarchi, and M. Vespe.  Measuring the impact of covid-19 confinement measures on human mobility using mobile positioning data, 2020.  JRC Report EUR-30290-EN. DOI 10.2760/913067.
  * [3] S. Iacus, C. Santamaria Serna, F. Sermi, S. Spyratos, D. Tarchi, and M. Vespe.  Mapping mobility functional areas (mfa) using mobile positioning data to inform covid-19 policies, 2020.  JRC Report EUR-30291-EN. DOI 10.2760/076318.
  * [4] N. Oliver and al.  Mobile phone data and covid-19: Missing an opportunity?, May 2020.  https://arxiv.org/pdf/2003.12347.pdf.
  * [5] M. Tiru.  Overview of the sources and challenges of mobile positioning data for statistics, 2014.  https://unstats.un.org/unsd/trade/events/2014/Beijing/Margus%20Tiru%20-%20Mobile%20Positioning%20Data%20Paper.pdf.
  * [6] F. De Meersman et al.  Assessing the quality of mobile phone data as a source of statistics.  In European Conference on Quality in Official Statistics, June 2016\.  https://tinyurl.com/y7pbracn.
  * [7] D. Salgado and al.  A proposed production framework with mobile network data, May 2020.  ESSnet Big Data II - Deliverable I.3 (draft version). Available online from https://webgate.ec.europa.eu/fpfis/mwikis/essnetbigdata/index.php/WPI_Milestones_and_deliverables.
  * [8] F. Ricciato, G. Lanzieri, and A. Wirthmann.  Towards a methodological framework for estimating present population density from mobile network operator data.  Working Paper (submitted), 2020.  https://europa.eu/!Hm66Ng.
  * [9] F. Ricciato, P. Widhalm, M. Craglia, and F. Pantisano.  Estimating population density distribution from network-based mobile phone data.  JRC Technical Report, 2015.  https://tinyurl.com/ydz4mgaw.
  * [10] F. Ricciato, P. Widhalm, M. Craglia, and F. Pantisano.  Beyond the “single-operator, cdr-only” paradigm: An interoperable framework for mobile phone network data analyses and population density estimation.  Pervasive and Mobile Computing, May 2016.
  * [11] M. Tennekes, Y. Gootzen, and S. H. Shah.  A bayesian approach to location estimation of mobile devices from mobile network operator data.  CBDS Working paper 06-20 https://www.cbs.nl/en-gb/background/2020/22/geographic-location-estimation-of-mobile-devices, 2020\.
  * [12] F. Ricciato.  Towards a reference methodological framework for processing mno data for official statistics.  In 15th Global Forum on Tourism Statistics, Cusco, Peru, November 2018.  https://tinyurl.com/ycgvx4m6.
  * [13] F. Calabrese, M. Colonna, P. Lovisolo, D. Parata, and C. Ratti.  Real-time urban monitoring using cell phones: A case study in rome.  IEEE Transactions on Intelligent Transportation Systems, 12(1), 2011\.
  * [14] M. Tennekes.  Statistical inference on mobile phone network data.  Presentation at European Forum for Geography and Statistics (EFGS 2018) https://www.efgs.info/wp-content/uploads/conferences/efgs/2018/presentations/TENNEKES_V1_EFGS2018-Statistical-interference-on-mobile-phone-network-data.pdf, 2018\.
  * [15] Center for Spatial Information Science — Univ. of Tokyo.  A study on urban mobility and dynamic population estimation by using aggregate mobile phone sources.  http://www.csis.u-tokyo.ac.jp/dp/115.pdf.
  * [16] D. Bachir, G. Khodabandelou, V. Gauthier, M. El Yacoubi, and J. Puchinger.  Inferring dynamic origin-destination flows by transport mode using mobile phone data.  Transportation Research Part C: Emerging Technologies, 101, 2019\.  https://doi.org/10.1016/j.trc.2019.02.013.
  * [17] G. Khodabandelou, V. Gauthier, M. El-Yacoubi, and M. Fiore.  Population estimation from mobile network traffic metadata.  In IEEE WoWMoM, 2016.  DOI 10.1109/WoWMoM.2016.7523554.
  * [18] R Lagonigro, R. Oller, and Martori J.C.  A quadtree approach based on european geographic grids: reconciling data privacy and accuracy.  SORT Statistics and Operations Research Transactions, 41(1):139–158, 2017.
  * [19] A. Janecek, D. Valerio, K. Hummel, F. Ricciato, and H. Hlavacs.  The cellular network as a sensor: From mobile phone data to real-time road traffic monitoring.  IEEE Transactions on Intelligent Transportation Systems, 16(5), October 2015.
  * [20] M. Tennekes.  R package for mobile location algorithms and tools, April 2017.  https://github.com/MobilePhoneESSnetBigData/mobloc.
  * [21] L. Shepp and Y. Vardi.  Maximum likelihood reconstruction for emission tomography.  IEEE Transactions on Medical Imaging, 1982.
  * [22] A. Raue, C. Kreutz, F. J. Theis, and J. Timmer.  Joining forces of bayesian and frequentist methodology: a study for inference in the presence of non-identifiability.  Phil. Trans. R. Soc., 371, 2013.  http://jeti.uni-freiburg.de/papers/20110544.full.pdf.
  * [23] B. Sakarovitch, M.P. de Bellefon, P. Givord, and M. Vanhoof.  Estimating the residential population from mobile phone data, an initial exploration.  Economie et Statistique, 505-506, 2018.  https://www.persee.fr/docAsPDF/estat_0336-1454_2018_num_505_1_10870.pdf.
  * [24] J. van der Laan and E. de Jongey.  Maximum likelihood reconstruction of population densities from mobile signalling data.  In NetMob’19, 2019.
  * [25] H.M. Taylor and S. Karlin.  An introduction to Stochastic Modeling, 3rd edition.  Academic Press, 1998.
  * [26] F. Bassetti, S. Gualandi, and M. Veneroni.  On the computation of kantorovich-wasserstein distances between 2d-histograms by uncapacitated minimum cost flows, July 2019.  arXiv:1804.00445v3.


# Finite groups with many isolated subgroups

Marius Tărnăuceanu

(November 28, 2021)

###### Abstract

Given a finite group $G$, we denote by $L(G)$ the subgroup lattice of $G$ and
by ${\rm Isolated}(G)$ the set of isolated subgroups of $G$. In this note, we
describe finite groups $G$ such that $|{\rm Isolated}(G)|=|L(G)|-k$, where
$k=0,1,2$.

MSC 2020 : Primary 20D30; Secondary 20D60.

Key words : finite group, subgroup lattice, isolated subgroup, ${\rm
CP}_{1}$-group.

##  1 Introduction

Let $G$ be a finite group and $L(G)$ be the subgroup lattice of $G$. We say
that a subgroup $H\in L(G)$ is isolated in $G$ if for every $x\in G$ we have
either $x\in H$ or $\langle x\rangle\cap H=1$. Groups with isolated subgroups
were studied in [1, 2]. However, this concept appears much earlier (see for
instance Section 66 of [5] and the entry ”isolated subgroup” in the
Encyclopedia of Mathematics, c.f. [11]).

Recently, the set ${\rm Isolated}(G)$ of isolated subgroups of $G$ has been
investigated for certain classes of non-abelian $p$-groups (see [4]) and for
abelian groups (see [10]). Also, the structure of finite $p$-groups with
isolated minimal non-abelian subgroups and of finite $p$-groups with an
isolated metacyclic subgroup has been determined in [8].

The current note deals with the following natural question:

How large can the set ${\rm Isolated}(G)$ be?

Our main result is stated as follows. We recall that a group is called a ${\rm
CP}_{1}$-group if all its elements have prime order; the structure of such a
group is given by [3] (see also [6]).

######  Theorem 1.1.

Let $G$ be a finite group. Then

  * a)

${\rm Isolated}(G)=L(G)$ if and only if $G$ is a ${\rm CP}_{1}$-group;

  * b)

$|{\rm Isolated}(G)|=|L(G)|-1$ if and only if $G\cong\mathbb{Z}_{p^{2}}$ with
$p$ a prime;

  * c)

$|{\rm Isolated}(G)|=|L(G)|-2$ if and only if $G\cong\mathbb{Z}_{p^{3}}$ with
$p$ a prime or $G\cong\mathbb{Z}_{pq}$ with $p$ and $q$ distinct primes.

The classification of finite groups $G$ such that

$|{\rm Isolated}(G)|=|L(G)|-k$

can be continued for $k\geq 3$. We remark that cyclic groups of order
$p^{k+1}$, where $p$ is a prime, or of order $p_{1}p_{2}\cdots p_{k}$, where
$p_{1},p_{2},\dots,p_{k}$ are distinct primes, satisfy this property. In other
words, for every positive integer $k$ there is a finite group with exactly $k$
non-isolated subgroups.

For the proof of the above theorem, we need the following well-known results:

  * $\bullet$

A finite group with one or two maximal subgroups is cyclic.

  * $\bullet$

A finite group has a unique minimal subgroup if and only if it is either a
cyclic $p$-group or a generalized quaternion $2$-group.

  * $\bullet$

A finite group has exactly two minimal subgroups if and only if it is either
cyclic of order $p^{m}q^{n}$ with $p,q$ distinct primes or a direct product of
a cyclic $p$-group of odd order and a generalized quaternion $2$-group.

By a generalized quaternion $2$-group we mean a group of order $2^{n}$ for
some positive integer $n\geq 3$, defined by

$Q_{2^{n}}=\langle a,b\mid
a^{2^{n-2}}=b^{2},a^{2^{n-1}}=1,b^{-1}ab=a^{-1}\rangle.$

Finally, we note that the dual problem of finding finite groups with few
isolated subgroups seems to be also very interesting. We say that a finite
group is isolated-simple if it has no proper isolated subgroups. Some obvious
examples of such groups are cyclic groups, generalized quaternion $2$-groups
or abelian $p$-groups of type
$\mathbb{Z}_{p^{\alpha_{1}}}\times\mathbb{Z}_{p^{\alpha_{2}}}\times\cdots\times\mathbb{Z}_{p^{\alpha_{k}}}$
with $2\leq\alpha_{1}\leq\alpha_{2}\leq...\leq\alpha_{k}$ (see Lemma 2.3 of
[10]).

Open problem. Determine the structure of isolated-simple groups.

Most of our notation is standard and will not be repeated here. Basic
definitions and results on group theory can be found in [5, 9]. For subgroup
lattice concepts we refer the reader to [7].

##  2 Proof of Theorem 1.1

a) Assume that ${\rm Isolated}(G)=L(G)$, but $G$ is not a ${\rm
CP}_{1}$-group. Then there exists $a\in G$ with $o(a)=pq$, where $p$ and $q$
are possibly equal primes. Since $a\notin\langle a^{p}\rangle$ and $\langle
a\rangle\cap\langle a^{p}\rangle=\langle a^{p}\rangle\neq 1$, it follows that
$\langle a^{p}\rangle$ is not isolated, a contradiction.

Conversely, let $H\in L(G)$ and $a\notin H$. Then $\langle a\rangle\cap
H\leq\langle a\rangle$ and $\langle a\rangle\cap H\neq\langle a\rangle$. Since
$\langle a\rangle$ is of prime order, we get $\langle a\rangle\cap H=1$. This
shows that $H\in{\rm Isolated}(G)$, as desired.

b) Assume that ${\rm Isolated}(G)=L(G)\setminus\\{H_{0}\\}$. Then all
conjugates of $H_{0}$ are also not isolated and therefore $H_{0}$ is normal in
$G$.

Let $a\in G\setminus H_{0}$ such that $\langle a\rangle\cap H_{0}\neq 1$. Then
$a\notin\langle a\rangle\cap H_{0}$ and

$\langle a\rangle\cap(\langle a\rangle\cap H_{0})=\langle a\rangle\cap
H_{0}\neq 1,$

which implies that $\langle a\rangle\cap H_{0}$ is not isolated in $G$. It
follows that $\langle a\rangle\cap H_{0}=H_{0}$, i.e. $H_{0}\subseteq\langle
a\rangle$. Clearly, $H_{0}$ cannot have proper subgroups because they would
not be isolated, contradicting our assumption. Thus, $|H_{0}|=p$ with $p$ a
prime. We also observe that any proper subgroup of $\langle a\rangle$ must
coincide with $H_{0}$ because it is not isolated. This leads to $|\langle
a\rangle|=p^{2}$.

Suppose that $G$ possesses a subgroup $H$ of order $p^{3}$ containing $\langle
a\rangle$. Then $H$ is of one of the following types: $\mathbb{Z}_{p^{3}}$,
$\mathbb{Z}_{p}\times\mathbb{Z}_{p^{2}}$, the non-abelian group of order
$p^{3}$ and exponent $p^{2}$, $D_{8}$ or $Q_{8}$. It is easy to see that if
$H\cong\mathbb{Z}_{p^{3}}$ then $\langle a\rangle$ is not isolated, while if
$H\cong D_{8}$ then the non-cyclic subgroups of order $4$ of $H$ are not
isolated, contrary to our assumption. In the other cases there exists $\langle
b\rangle\leq H$ with $|\langle b\rangle|=p^{2}$ and $\langle
b\rangle\cap\langle a\rangle=H_{0}$, proving that $\langle a\rangle$ is again
not isolated, a contradiction. Thus, $\langle a\rangle$ is a Sylow
$p$-subgroup of $G$.

Suppose now that $|G|$ has a prime divisor $q\neq p$ and let $K\leq G$ with
$|K|=q$. Then $H_{0}K\leq G$, $a\notin H_{0}K$ and

$\langle a\rangle\cap H_{0}K=H_{0}\neq 1,$

showing that $H_{0}K$ is not isolated in $G$, a contradiction.

Hence $G=\langle a\rangle\cong\mathbb{Z}_{p^{2}}$, as desired.

c) Assume that ${\rm Isolated}(G)=L(G)\setminus\\{H_{1},H_{2}\\}$. Similarly,
$H_{1}$ and $H_{2}$ are normal in $G$.

Let $a\in G\setminus H_{1}$ such that $\langle a\rangle\cap H_{1}\neq 1$. Then
$\langle a\rangle\cap H_{1}$ is not isolated, implying that $\langle
a\rangle\cap H_{1}\in\\{H_{1},H_{2}\\}$. We distinguish the following two
cases:

Case 1. $\langle a\rangle\cap H_{1}=H_{1}$

Then $H_{1}\subset\langle a\rangle$. Since proper subgroups of $\langle
a\rangle$ are not isolated, we infer that they are at least one, namely
$H_{1}$, and at most two, namely $H_{1}$ and $H_{2}$. So, we have two
possibilities:

Subcase 1.1. $H_{1}$ is the unique proper subgroup of $\langle a\rangle$

Then $|\langle a\rangle|=p^{2}$ and $|H_{1}|=p$, where $p$ is a prime.

Suppose first that $H_{1}$ is the unique minimal subgroup of $G$. This implies
that $G$ is either a cyclic $p$-group or a generalized quaternion $2$-group.
Then it is easy to see that $G$ has exactly two non-isolated subgroups if and
only if $G\cong\mathbb{Z}_{p^{3}}$.

Suppose now that $G$ possesses minimal subgroups different from $H_{1}$ and
let $K$ be such a subgroup. Then $a\notin H_{1}K$ and

$\langle a\rangle\cap H_{1}K=H_{1}\neq 1,$

that is $H_{1}K$ is not isolated in $G$. Since $H_{1}K\neq H_{1}$, we get
$H_{1}K=H_{2}$. Thus $H_{2}$ contains all minimal subgroups of $G$.

Observe also that $\langle a\rangle$ is normal in $G$. Indeed, for every $x\in
G$, if $\langle a\rangle^{x}\neq\langle a\rangle$ then $\langle
a\rangle\cap\langle a\rangle^{x}=1$ because $\langle a\rangle^{x}$ is isolated
in $G$. It follows that $\langle a^{p}\rangle^{x}$ is not isolated in $G$ and
since $\langle a^{p}\rangle^{x}\neq H_{1},H_{2}$ we get a contradiction.

Let $H\neq\langle a\rangle,H_{2}$ be a maximal subgroup of $\langle a\rangle
H_{2}$. Then $H$ is isolated in $G$ and therefore $\langle a\rangle\cap H=1$.
Since $\langle a\rangle H=\langle a\rangle H_{2}$, it follows that $H$ is a
minimal subgroup of $G$, which leads to $H\subset H_{2}$. Thus $\langle
a\rangle H_{2}$ has exactly two maximal subgroups, namely $\langle a\rangle$
and $H_{2}$, and consequently it is cyclic. This implies that $\langle
a\rangle$ is not isolated in $G$, a contradiction.

Subcase 1.2. $H_{1}$ and $H_{2}$ are the unique proper subgroups of $\langle
a\rangle$

Then we have either

$|\langle a\rangle|=p^{3}\mbox{ and }|H_{1}|=p,|H_{2}|=p^{2},\mbox{ where
}p\mbox{ is a prime}$

or

$|\langle a\rangle|=pq\mbox{ and }|H_{1}|=p,|H_{2}|=q,\mbox{ where }p,q\mbox{
are distinct primes}.$

In both cases, if $G$ would have a minimal subgroup $K$ different from $H_{1}$
and $H_{1},H_{2}$, respectively, then $a\notin H_{1}K$ and

$\langle a\rangle\cap H_{1}K=H_{1}\neq 1,$

i.e. $H_{1}K$ is non-isolated in $G$. Since obviously $H_{1}K\neq
H_{1},H_{2}$, this contradicts our hypothesis. Thus $H_{1}$ and $H_{1},H_{2}$,
respectively, are the unique minimal subgroups of $G$, proving that $G$ is one
of the following groups:

  * -

$\mathbb{Z}_{p^{m}}$, where $p$ is a prime and $m\in\mathbb{N}$;

  * -

$Q_{2^{n}}$, where $n\in\mathbb{N}$, $n\geq 3$;

  * -

$\mathbb{Z}_{p^{m}q^{n}}$, where $p,q$ are distinct primes and
$m,n\in\mathbb{N}$;

  * -

$\mathbb{Z}_{p^{m}}\times Q_{2^{n}}$, where $p$ is an odd prime and
$m,n\in\mathbb{N}$, $n\geq 3$.

In all cases one can easily check that $G$ has exactly two non-isolated
subgroups if and only if it coincides with $\langle a\rangle$, as desired.

Case 2. $\langle a\rangle\cap H_{1}=H_{2}$

Then $H_{2}\subseteq\langle a\rangle$ and $H_{2}\subset H_{1}$. If
$H_{2}=\langle a\rangle$, then we get $\langle a\rangle\subseteq H_{1}$, a
contradiction. Thus $H_{2}\neq\langle a\rangle$. Again, proper subgroups of
$\langle a\rangle$ are not isolated in $G$. If $H_{2}$ is the unique proper
subgroup of $\langle a\rangle$, then we proceed similarly with Subcase 1.1,
while if $H_{1}$ and $H_{2}$ are the unique proper subgroups of $\langle
a\rangle$, then we proceed similarly with Subcase 1.2.

This completes the proof.∎

## References

  * [1] V.M. Busarkin, The structure of isolated subgroups in finite groups, Algebra i Logika, Sem. 4 (1965), 33-50.
  * [2] V.M. Busarkin, Groups containing isolated subgroups, Sibirsk. Mat. Zh. 9 (1968), 752-756.
  * [3] M. Deaconescu, Classification of finite groups with all elements of prime order, Proc. Amer. Math. Soc. 106 (1989), 625–629.
  * [4] Z. Janko, Finite $p$-groups with some isolated subgroups, J. Algebra 465 (2016), 41-61.
  * [5] A.G. Kurosh, The theory of groups, 2nd English ed., 2 volumes, Chelsea Publishing Co., New York, 1960.
  * [6] K.N. Cheng, M. Deaconescu, M.L. Lang and W. Shi, Corrigendum and addendum to ”Classification of finite groups with all elements of prime order”, Proc. Amer. Math. Soc. 117 (1993), 1205-1207.
  * [7] R. Schmidt, Subgroup lattices of groups, de Gruyter Expositions in Mathematics 14, de Gruyter, Berlin, 1994.
  * [8] Q. Song and L. An, Isolated subgroups of finite $p$-groups, to appear in Bull. Austral. Math. Soc., 2021.
  * [9] M. Suzuki, Group theory, I, II, Springer Verlag, Berlin, 1982, 1986.
  * [10] M. Tărnăuceanu, Isolated subgroups of finite abelian groups, to appear in Czech. Math. J., arXiv:2102.04598.
  * [11] Isolated subgroup, Encyclopedia of Mathematics, Webpage url: https://encyclopediaofmath.org/wiki/Isolated subgroup.

Marius Tărnăuceanu

Faculty of Mathematics

“Al.I. Cuza” University

Iaşi, Romania

e-mail<EMAIL_ADDRESS>