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 F. Ricciato is with European Commission, DG Eurostat, Unit B1 Methodology; Innovation in Official Statistics, Luxembourg. Email: fabio.ricciato@ec.europa.eu.
A. Coluccia is with University of Salento, Lecce, Italy. Email: angelo.coluccia@unisalento.it.
The views expressed in this paper are those of the authors and do not necessarily represent the official position of the European Commission. [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 rescaling, 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 geolocation 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 closedform 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.
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).
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 infrastructure 1 . A mobile device z that at time t interacts with the mobile network through the radio cell i generates an event record z, t, i . 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. 1 The term "transactional data" was used in an earlier study [12] as synonymous of "event data".
Call Detail Records (CDR) and the more informative signalling data fall in this category 2 . • Cell data. Such data binds the network infrastructure, i.e., the collection of "radio cells", to the geographical territory 3 . 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 t, i, s i (t) associating to each cell i at time 4 t a geographical profile, or "cell footprint", represented by the vector s i (also called "event location" in [7] ). The jth element of such vector, denoted by s ij ≥ 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 t, i as matching key to combine event records z, t, i and cell records t, i, s i (t) , 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 s i . Then the set of counters and footprints over all cell (or cell groups) {c i , s i } 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). 2 Under 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 subcell 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. 3 The term "auxiliary data" was used in [12] as synonymous of "cell data". 4 The 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).
• 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.
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 • 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.
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 ∈ {0, 1} 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 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, intercell 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 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.
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 jth element u j of the column vector u def = [u 1 u 2 · · · u J ] T denote the unknown number of mobile users in tile j = 1, . . . , J. Let the ith element c i of the column vector c def = [c 1 c 2 · · · c I ] T denote the observed number of mobile users detected (or counted) in cell i = 1, . . . , I (cell count vector). We denote the total number of users across all cells by C
where the symbol 1 k denotes a column vector of size k with all elements equal to 1. Denote by p ij ∈ [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:
For the sake of a more compact notation we gather the individual probabilities p ij 's into a matrix P [I×J] . The individual elements of 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 P sum to one along columns, formally 1 T I P = 1 J , meaning that the matrix P is column-stochastic 5 . We consider for this work the ideal case that the matrix P is known exactly. In other words, do not consider model mismatching errors due to inaccurate knowledge of 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 c can be interpreted as the single realization of a random vectorc whose expected value is given by:
In the estimation problem we must solve for estimand u given the vector of measurement data c, representing the single available observation ofc, and the model matrix P . Being a sort of inversion problem, the estimateû can be written in general as:û = g(P , c)
where g(·) 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 ≥ 0, ∀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. 5 In 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. Remark 1. In cases of practical interest the number of tiles is (much) larger than the number of cells, i.e., J I. Therefore, even if the average value 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(·) 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 ij1 = p ij2 ∀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 j1 and u j2 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 latter 6 . Algebraically, this corresponds to merging together identical columns of matrix 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.
For the sake of completeness, we briefly illustrate hereafter how the assignment matrix P can be derived from the cell footprints 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
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:
= 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 ∈ {0, 1}. In other words, when cells are mutually non-overlapping, the data generating process becomes deterministic: the binary model matrix P after consolidation reduces to the identity matrix, and the estimation problem becomes trivial 7 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.
Whatever approach is chosen to populate the model matrix P during the modeling stage, a suitable estimation procedure (represented by the function g(·) 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.
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 µ
Clearly, µ is a vector with non-negative elements summing to unity, i.e., µ ≥ 0 and µ 1 = 1 T J µ = 1. The resulting number of units in the tiles is distributed as a multivariate random vectorũ with Multinomial distribution (denoted by M) of parameters C and µ, i.e.
withũ ≥ 0 and ũ 1 = C. Therefore by construction the mean value writes u def = E[ũ | C, µ] = Cµ. In the second layer (ref. Fig. 1(a) ) the units are assigned randomly and independently from tiles to cells according to the assignment matrix 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 vectorc has also a Multinomial distribution with parameters C and P µ:
Thus, the log-likelihood function for an observed value c is derived as (neglecting an irrelevant constant factor):
wherein we have used the logarithm of vector notation
T to refer compactly to the vector of element-wise logarithms. Letμ ML denote the value of the probability vector µ that maximizes the log-likelihood (7) subject to the constraints µ ≥ 0 and µ 1 = 1, formally:
By definitionμ ML represents the Maximum Likelihood (ML) estimate for µ given the observed data c. Givenμ ML , we still have to provide a point estimate for 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:û
where we retain the label "ML" for simplicity. This equation represents a simple rescaling ofμ ML by the factor C. Therefore, with a simple variable substitution, the minimization (8) can be rewritten to deliver directly the final estimateû ML , leading to the following constrained optimization problem:
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.
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ũ are modelled as independent Poisson (instead of Multinomial) random variables with parameters ρ def = [ρ 1 · · · ρ J ] T , i.e.,ũ j ∼ P(ρ j ) (P denoting the Poisson distribution). Following the same reasoning that led to (9) , the ML estimate of ρ is taken as the estimate for 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 estimatê u m+1 j is computed from the previous estimateû m j according to the following formula (see [21, eq. (2.13)] or, equivalently, [24, eq 
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û 0 j = C J , ∀j.
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:
Note the inversion of the conditioning direction between q ji and p ij defined earlier in (1). Let
denote the prior probability that a single generic unit falls in tile j, before observing the measurement data. Recalling the Bayes rule
Therefore, the estimate in each tile j is computed directly aŝ
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 aŝ
where diag(v) is a diagonal matrix containing the entries of vector v (and diag −1 (v) their reciprocals).
Proof. We first notice that J k=1 p ik α k and I i=1 ξ i p ij are the ith and jth element, respectively, of the vectors P α and P T ξ, respectively, with ξ i = c i / J k=1 p ik α k . Then the thesis follows straight by rewriting in vector form ξ = diag −1 (P α)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 P and data c.
All the three methods reviewed above provide a point estimateû in output, based on the model matrix P and observed data 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 α def = [α 1 · · · α J ] T (with α 1 = 1 and α ≥ 0) or, equivalently, of its rescaled version a def = Cα (with a 1 = C and a ≥ 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 Vardi 8 , it represents a sort of prior information, and for this reason we shall refer to α (or equivalently to 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 α (non rescaled) is explicit in the SB estimator (14) where its elements represent prior probabilities (per tiles). As for MLE-Poisson, the prior vector a (rescaled) serves as the initial point for the iterative procedure (11) by settingû 0 = a. Similarly, the prior vector 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 α 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 jth 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.
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.
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.
wherein p T i denotes the ith row of P . Thus, for any vector x = [x 1 · · · x I ] 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 P are non-negative (of which at least one strictly positive element for each column), we can write J k,h=1
It follows that x T Hx ≤ 0, with equality only if y = [y 1 · · · y I ] T = 0. We can rewrite
. Since D has full rank I, DP has the same rank of P ; the latter has rank that cannot exceed I < J, hence the rank of DP is smaller than J. We conclude that there exist some x = 0 such that y = DP x = 0, hence H is negative semidefinite.
The equivalence between MLE-Multinomial and MLE-Poisson (i.e., EM algorithm (11)) is established as follows. Proof. We start by considering the iterative EM procedure derived from MLE-Poisson: at convergence, the stable points of (11) must fulfill the conditionû m+1 =û m (m denoting the iteration index) therefore for the solution vectorû it holds:
wherein p T i denotes the ith row of 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 ith element of vector P u is p T i u, we obtain, for each j:
Multiplying by u j and summing over j we obtain
and since j p ij u j = p T i u it turns out that λ = 1. By substituting back in the first-order optimality condition it follows that, for each j,
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 space 9 .
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 ∼ P(ρ k ), k = 1, . . . , K, then the distribution of the vector
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.
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 guess 10 .
Recall that the vector c represents a single realization of the random vectorc with unknown mean c. Since no other measurements are available, in the absence of external information the vector c represents a natural estimate of the mean value c. Replacing the unknown term c with its estimate c in the generative relation (2) leads to the constraint P u = c. It is rather intuitive that, among all possible (non-negative) values of u, those respecting this constraint are the ones that best conform with the available measurement c and with the 9 Strictly 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. 10 Such purely multiplicative procedure in the form u m+1 j = u m j · ψ m j is distinct from the standard gradient descent approaches that, in general, involves an additive update in the form u m+1
model P . This argument leads us to restrict the search for a "good" estimateû within the bounded subspace U defined formally by:
Note that the constraint P u = c absorbs the condition u 1 = C on the total count. Since in our application J I, matrix P is rank deficient and the constraint P u = c admits multiple solutions in the variable 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) Proof. The condition P u = c can be rewritten as p T i u = c i , ∀i, and that implies that the optimality condition (16), or equivalently (15), is always verified: with the ansatz p T i u = c i , ∀i, (16) reduces to the identity I i=1 p ij = 1 that is always satisfied by the column-stochasticity of P , i.e., 1 T I P = 1 J . Therefore, we have established that U represents a set of solutions for both MLE procedures, i.e., any point in U is also a solution of MLE-Multinomial and MLE-Poisson.
Following Proposition 4, we shall refer to U as the Maximum Likelihood Bounded Subspace (MLBS). (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 α, as can be immediately verified by comparing (14) and (11), the only difference being thatû m j in (11) is replaced by α j in (14) . In other words, the SB point lies between the prior vector α and the MLBS.
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 u j = C, or equivalently u 1 = c 1 . Therefore, we can map all admissible solutions 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.
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.
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 u can be obtained aŝ
Proof. Following the stochastic model shown in Fig. 1(a) , u represents a random vector with Multinomial distribution of unknown parameter µ. Taking the prior vector α as the prior value (before seeing the data) of µ, we obtain that the prior probability distribution of u, denoted by P(u), is Multinomial with parameters C and α, formally (20) from where the logarithm log P(u) = u T log α − J j=1 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 P(u). Accordingly, the logarithm of the posterior distribution equals the sum of the log-likelihood (7) and log P(u). Therefore, maximizing the (logarithm of) posterior distribution leads to maximizing the following function of µ and u arg max
Recalling (9) we bind the estimateû to equal the rescaled value ofμ, i.e., we imposeû = Cμ. With this simple variable substitution, from (21) we obtain the following MAP estimator 11 :
which can be rewritten to finally yield the thesis (note that the term α can be replaced by a = α/C for the same invariance property recalled in the previous footnote). 11 In the derivation we have used the following invariance property: a constant factor κ in the argument of the logarithm terms does not influence the solution to the maximization. Formally:
In eq. (19) the term c T log P u carries the data information (from the measurements c) while the other two terms carry information from the prior distribution P(u). More specifically, the term u T log a carries information from the prior vector 12 a, or equivalently α, while the term J j=1 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 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 Γ(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.
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:
Proof. We start by approximating the (prior) Multinomial distribution with a multivariate normal
with mean value u = a = Cα and covariance matrix Σ a composed of the following elements:
Plugging the normal approximation (24) into (22) leads to the approximate MAP estimator
In our application we are considering very large number of tiles (J 1), which implies that the individual probabilities are very small in absolute terms (α j 1). This justifies the further simplification of the covariance matrix: neglecting the product terms in the off-diagonal elements, Cα m α n 0, and approximating the diagonal elements with C α m (1 − α m ) C α m = a m , the covariance matrix reduces to the diagonal matrix A def = diag([a 1 · · · a J ]) and its inverse to Σ −1
. In this way the approximate MAP estimator rewriteŝ
(27) which finally yields the thesis.
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 nonidentifiable, 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 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 c, with information before the data, captured by the prior distribution embedding the prior vector α (or equivalently a). In fact, the posterior distribution is composed by the product of the two components: posterior = likelihood × 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).
The goal expressed above is achieved though an estimator built according to the following structure:
with f (u, a) denoting a distance function between vectors u and a. With this structure, the ML property is imposed as hard constraint. In fact, the condition
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 c. Among these points, then we select the point that additionally best conforms with the off-data knowledge encoded in the prior vector a, through the minimization of some distance function f u, a . Among many possible choices for the distance function, opting for the weighted 2 -norm defined in Proposition 6:
with A def = diag(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"):
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 bŷ
where the maximum function is intended element-wise.
Proof. The constrained minimization problem can be solved by considering the Karush-Kuhn-Tucker conditions:
The calculation of the gradient in the first line yields
The vector of the Lagrange multipliers λ corresponding to the equality constraint can be explicitly obtained by exploiting the constraint itself, that is P u = c ⇔ P Aν − P AP T λ + P a = c which leads to λ = (P AP T ) −1 (P Aν + P a − c).
Substituting back one obtainŝ (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 c T log P u and (u−a) T A −1 (u−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 P u = 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.
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 P and on the prior a, and is independent from the data c. In fact, setting F def = AP T (P AP T ) −1 and g def = (F P − I J )a, eq. (32) can be rewritten aŝ
where the maximum function is intended element-wise. While matrix F and vector g are heavy to compute (especially because F embeds a matrix inversion) they are both independent from the measured data c. Therefore, in scenarios where multiple estimates must be computed for the same model matrix P but with different measurement vectors, F and 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 c.
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.
We briefly describe the data generation scenario. The reference area consists of a regular square grid of 100 × 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 • sectors with a maximum range (normalized to the tile size) equal to 15 and azimuth equal to π/5. The second layer consists of 73 antenna towers (blue circles in Fig. 5 (a)) with each tower hosting three 120 • sectors with a maximum range equal to 6 and 0 • 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., ×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 = 1 kj . The resulting matrix P has nonnegative fractional values. The same matrix P is used in the data generation process and in the estimation models, i.e., we test the different methods in model matching condition.
In order to evaluate the goodness of a generic estimation method X we must resort to some similarity metric w(û X , u GTP ) between the estimated densityû X and the ground truth population u 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 u and v defined over the same grid B and with same total mass u 1 = v 1 = C, KWD may be interpreted as the minimum cost of transporting the mass from configuration u to 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:
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 code 13 which we used to compute the KWD (with 13 Available at https://github.com/stegua Euclidean ground distance and lattice parameter L = 4) between each estimate and GTP.
A physical interpretation of KWD w(û X , u GTP ) in our context is the following. Imagine that, starting from the estimated densityû X , we move each population unit (i.e., a mobile device) across tiles so as to arrive to the final GTP configuration u 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 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).
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 findings 14 . 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 ×2 over Voronoi (from 4 to 1.9) and with ML/EM the gain factor increases to ×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 14 Similar 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. 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 ×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 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 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.
The analytic and numerical results presented in this study indicate that, when the model matrix 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) geolocation 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.