Stochastic Geometry Meets Optimization: Data-Driven Cellular Planning

1

Stochastic Geometry Meets Optimization: Data-Driven Cellular Planning

arXiv:2610.31793v1 [cs.IT] 22 Oct 2026

Abhishek K. Shukla, Abhishek K. Gupta, and Ketan Rajawat

Abstract—Deployment optimization of base stations (BSs) is critical for improving cellular network performance metrics such as coverage and capacity, impacting user experience, scalability, and infrastructure cost-effectiveness. Traditional methods that optimize precise BS locations do not generally provide broader insights into spatial deployment priorities or alternative viable deployments. Stochastic geometry approaches, on the other hand, offer a more interpretable intensity-based framework, but rarely utilize real-world user location data or investigate optimal BS deployment strategies. In this work, we propose a data-driven stochastic geometry-based approach that optimizes BS deployment intensity to maximize a network performance metric using samples from an unknown user-location distribution. We specialize the framework to an inhomogeneous Poisson point process (PPP) for cellular coverage analysis. The resulting infinite-dimensional problem is discretized over the deployment area and solved using sampled user locations. We characterize the discretization and sampling errors, apply stochastic mirror ascent to solve the empirical problem, and analyze the selection of concrete BS deployments using independent validation samples. Unlike traditional methods, our approach yields an optimized spatial intensity that reveals deployment priorities across regions and accommodates practical deployment constraints, multi-tier networks, and incremental deployment scenarios. Numerical experiments with synthetic and real-world spatial demand data demonstrate the proposed approach under different network configurations. Index Terms—Stochastic geometry, base station deployment, optimization

I. I NTRODUCTION In cellular networks, deployment optimization refers to the problem of determining where to place the BSs or access points in a network, so as to maximize performance metrics such as user coverage and network capacity. As mobile networks become increasingly dense, heterogeneous, and complex, strategic BS placement plays a critical role in network planning, directly impacting user experience, scalability, and the cost-effectiveness of infrastructure investments [1], [2]. The concept of deployment optimization encompasses a range of scenarios, including selectively deactivating underutilized BSs for energy savings [3], [4] and temporarily deploying mobile or UAV-based BSs to address localized surges in demand [5], [6], [7], underscoring its broad applicability in modern network designs. The authors are with the Department of Electrical Engineering, Indian Institute of Technology Kanpur, Kanpur 208016, India (emails: abhisheks20@iitk.ac.in; gkrabhi@iitk.ac.in; ketan@iitk.ac.in). Part of this work was presented at the 2026 Allerton Conference. Supplementary material includes the MATLAB codes and additional proofs, available at https://home.iitk.ac.in/∼ ketan/Supplementary

Classical approaches have focused on finding the exact locations of BSs by directly optimizing a relevant performance metric, using a range of mathematical programming approaches [8], heuristics [9], and game-theoretic formulations [10]. These methods seek a single deployment for a given network state and feasible-site set. Such a solution is appropriate when the exact sites are known a priori, but provides limited guidance when site availability, deployment constraints, or network conditions change as the network evolves. In particular, it does not characterize alternative viable deployments or indicate which regions remain attractive for deployment under related network configurations. Recently, stochastic geometry has emerged as a tractable approach to model and analyze wireless networks by accounting for various practical effects including BS locations, fading coefficients, shadowing, blockage, and antenna directionality, under reasonable assumptions [11]. In particular, homogeneous Poisson point process models for users and BSs often yield simple expressions for network performance metrics such as coverage probability and rate coverage, and are widely used for analysis [12], [13] as well as BS intensity optimization [14]. More recent stochastic-geometry models also account for the spatial coupling between user hotspots and BS locations arising in user-centric small-cell deployments [15]. However, these works primarily analyze prescribed spatial models rather than optimize the BS intensity using observed user locations. A key challenge in computing the optimal BS deployment is the lack of knowledge of the underlying mobile user distribution. In theoretical studies, it is common to make simplifying assumptions that scatter users according to a stationary uniform distribution [11]. Such simplistic models, however, do not capture the behavior of real-world users, which may congregate at certain spots, avoid specific zones, or exhibit multi-scale behavior, motivating the need for datadriven modeling. In related contexts, recent studies have increasingly utilized actual user location data, collected from cellular traces [16] and Wi-Fi logs [17] over time, in order to build realistic user density maps using machine learning approaches such as kernel density estimation. We remark that an additional challenge here is that at any given time, only a subset of users is generally active and observable. In this work, we put forth a novel data-driven stochastic geometry-based approach for determining an optimal BS deployment to provide connectivity to users spread over a specified region. Despite the huge potential and multiple applications in 5G/6G network planning, such a formulation has never been considered in the literature. In particular, we seek to design an optimal BS deployment that provides coverage to

2

users whose underlying location distribution is unknown but observed via a few sample locations. To this end, we model the BSs as following an inhomogeneous PPP with intensity function λ(x). Our objective is to compute the optimal λ(x) that maximizes a specific network performance metric R(λ), such as coverage or throughput. Unlike traditional deployment approaches that yield precise BS locations, our approach yields an intensity function that reflects the optimal BS intensity across regions, offering insights into the deployment needs (e.g. densification vs. interference management), while also adapting to actual user patterns. The intensity function can also be interpreted as a probabilistic priority function for BS soft-deployment to guide an operator. Given the BS intensities, exact BS locations may subsequently be extracted while taking into account other geographical, economic, and social constraints. The proposed approach is also amenable to varied practical deployment constraints, multi-tier networks, and incremental deployment scenarios. Remarkably, it is known that under sufficiently strong shadowing, even deterministic BS deployments may exhibit propagation characteristics similar to those of Poisson networks [18], providing further motivation for studying intensity-based soft deployment. At the outset, determining the optimal BS intensity using real-world user data is challenging on two fronts: the intensity function λ(x) is infinite-dimensional and the underlying user distribution is unknown. To address the first challenge, we discretize the area of interest into D smaller sub-regions, so that λ(x) can be approximated with a piecewise constant function, rendering the optimization finite-dimensional. Second, to bypass the need for the unknown user distribution, we approximate R(λ) with the empirical average over the M sampled user locations. While the solution of the approximate problem is suboptimal with respect to the original problem, we establish that the expected optimality gap is no more than ǫ provided that D = Õ(ǫ−2 ) and M = Õ(ǫ−4 ). Interestingly, these bounds are generally applicable to a broader class of variational optimization problems, and may be of independent interest. Next, we provide a stochastic mirror ascent algorithm for solving the finite-dimensional empirical problem efficiently. We then derive the coverage probability expressions for an arbitrary BS intensity under instantaneous SINR-based and average SINR-based cell association rules, which can serve as the objectives of our maximization. We also propose a procedure for generating concrete BS deployments from the optimized intensity and selecting a suitable realization using independent validation samples, while quantifying the associated placement error. We highlight various practical applications of the proposed approach in 5G/6G networks, including multi-tier deployments, incremental deployment in the presence of existing BSs, and deployment with exclusion zones. Numerical investigations demonstrate the usefulness of the approach and provide insights into how the optimized BS intensity depends on user distribution, existing infrastructure, and deployment constraints. In particular, the optimized intensity reveals not only preferred deployment regions but also alternative viable locations and their relative priorities, which

are not evident from a single deterministic deployment. To the best of our knowledge, this is the first framework that learns a fully spatially varying BS point-process intensity directly from user-location samples while jointly quantifying discretization, statistical estimation, and realization-selection errors. Notation: We denote scalars by regular font letters, vectors by bold font letters, and sets by calligraphic capital letters. The notation k·k is used for both L2 and ℓ2 norms, depending on whether the argument is a function or a vector, respectively. II. P ROBLEM F ORMULATION We consider the design of a cellular network where the BSs follow a point process Φλ parameterized by an intensity function λ : A → R+ where A ⊂ R2 is a compact region with area A = |A| and diameter Ad . The users are also distributed according to an independent point process Ψu over A with user locations following distribution z ∼ Pz . Let r(λ, z) denote the performance metric of interest for a user located at z, such as coverage probability (CP) pc = P[SINR > τ ] or per-user rate coverage rc = P[Rate > τ ]. The average user performance is given by R(λ) = Ez∼Pz [r(λ, z)]. We seek an intensity λ that maximizes R(λ). Examples of R include the average SINR (or rate) coverage which denotes the fraction of users that are in coverage, i.e., their SINR (or rate) is above the target threshold. The optimized intensity λ provides a spatially adaptive guideline for BS deployment, where higher intensity values indicate regions requiring denser BS placement while lower values suggest sparser deployment. The problem of optimally deploying an average of N BSs in the region A can be written as λ⋆ = arg max E[r(λ, z)] λ∈L1+ (A) Z λ(x)dx = N s. t. A

(P) (1)

where L1+ (A) denotes the set of nonnegative Lebesgueintegrable functions on A. Denoting the Rset of admissible intensity functions by S := {λ ∈ L1+ (A) | A λ(x)dx = N }, we can write (P) compactly as λ⋆ = arg maxλ∈S R(λ). To avoid trivial cases, we will assume that maxλ∈S R(λ) is finite and achieved by some λ⋆ ∈ S. As discussed earlier, without an explicit user distribution, we cannot directly solve (P), as the expectation in the objective function cannot be evaluated. Instead, given M independent and identically distributed (i.i.d.) user-location samples Z := {zm }M replace R(λ) by m=1 , where zm ∼ Pz , we can 1 PM its empirical approximation R̂(λ) = M m=1 r(λ, zm ) ≈ E[r(λ, z)]. Such an empirical estimate accommodates arbitrary and spatially heterogeneous user distributions, and improves in accuracy as M increases. To obtain a finite-dimensional problem suitable for numerical optimization, we also need to discretize the intensity function. To this end, we partition the region A into D disjoint equal-area sub-regions denoted by A1 , . . ., AD such that ∪k Ak = A. The piecewise PD constant approximation of λ(x) is denoted by λd (x) := k=1 λdk 1 x∈Ak where Z 1 λ(x)dx (2) λdk = |Ak | Ak

3

and 1 denotes the indicator function which takes the value 1 if the condition in its subscript is true and zero otherwise. Observe further that since |Ak | = A/D, the left-hand side of the constraint in (1) simplifies as Z D D Z X AX λdk . (3) λ(x)dx = λ(x)dx = D A k=1 k=1 Ak P Therefore, the constraint (1) can be written as k λdk = N D/A. Collecting λd1 , . . ., λdD into the vector λd ∈ RD , we see that λd ∈ S is equivalent to λd ∈ Λ := {λd  0 | 1T λd = NAD }. Hence, the modified problem becomes M

1 X λ̂d = arg max r(λd , zm ) λd ∈Λ M m=1

(P̂)

where, with some abuse of notation, we allow functionals r, R, and R̂ to also take the vectors such as λd and λ̂d (in addition to the functions λd and λ̂d ) as inputs, i.e., r(λd , z) := r(λd , z) and R(λd ) := R(λd ). In a similar vein, the optimal value of (P̂) is given byPR̂(λ̂d ) := R̂(λ̂d ) where D we define the function λ̂d (x) := k=1 λ̂dk 1 x∈Ak . Hence, (P̂) can be seen as a discretized and sampled version of (P). As (P̂) is a finite-dimensional optimization problem, it can be solved using standard numerical methods. In the next section, we characterize the excess error associated with solving (P̂) instead of (P). In practice, the optimization algorithm itself may incur additional error, e.g., if (P̂) cannot be solved exactly due to nonconvexity. We do not analyze that component of the error in this section. III. B OUND ON THE EXCESS ERROR ⋆

We define the excess error as E = E[|R(λ ) − R(λ̂d )|] = E[|R(λ⋆ ) − R(λ̂d )|] where the expectation is with respect to the random samples {zm }M m=1 . As we describe below, this error consists of two factors arising from discretizing the intensity function and sampling the user distribution. This section characterizes E in terms of the average BS count N , discretization level D, and number of user-location samples M , thereby quantifying how these parameters affect the quality of the learned BS intensity. To separate the effects of discretization and sampling, we introduce the oracle discretized solution λ̃d = arg max R(λd ), λd ∈Λ

(4)

which optimizes the exact average metric R(λd ) = E[r(λd , z)] over the same discretized feasible set Λ, while the expectation is with respect to the true user distribution. Let λ̃d denote the piecewise constant intensity function corresponding to λ̃d . Then, E = E[|R(λ⋆ ) − R(λ̃d ) + R(λ̃d ) − R(λ̂d )|] ⋆

≤ |R(λ ) − R(λ̃d )| + E[|R(λ̃d ) − R(λ̂d )|]

(5) (6)

where we have used the triangle inequality and the fact that R(λ̃d ) = R(λ̃d ) as explained earlier. Here, the first term captures the loss due to restricting the original infinitedimensional problem to piecewise constant intensities, while

the second term captures the additional loss due to replacing the exact objective by its empirical approximation. Accordingly, we define the discretization and sampling errors as Ed := |R(λ⋆ ) − R(λ̃d )|,

(7)

Es := E[|R(λ̃d ) − R(λ̂d )|].

(8)

The subsequent subsections provide bounds on these two error terms, respectively. A. Bound on the discretization error Let λ⋆d denote the vector obtained by discretizing the ⋆ optimal λ⋆ as described earlier in (2) P so that λd (x) = P ⋆ k λ̃dk 1 x∈Ak . Since k λdk 1 x∈Ak . Similarly, let λ̃d (x) := λ̃d ∈ S, it follows from the definition of λ⋆ that R(λ⋆ ) ≥ R(λ̃d ). Likewise, since λ⋆d ∈ Λ, it follows from the definition of λ̃d that R(λ̃d ) ≥ R(λ⋆d ). Combining these inequalities, we obtain 0 ≤ R(λ⋆ ) − R(λ̃d ) ≤ R(λ⋆ ) − R(λ⋆d ) ⋆

⇒ Ed = |R(λ ) − R(λ̃d )| ≤ |R(λ ) − R(λ⋆d )|. ⋆

(9) (10)

It remains to bound |R(λ⋆ ) − R(λ⋆d )|, which captures the loss incurred by replacing λ⋆ with its piecewise constant approximation. To obtain such a bound, we need to impose continuity of R with respect to an appropriate metric. Since λ describes the spatial distribution of a fixed BS mass over A, we can use the Wasserstein-1 distance metric, also known as the earth mover’s distance. Specifically, let L1 be the set of all 1Lipschitz functions on A so that L1 = {ϕ : A → R | |ϕ(x) − ϕ(y)| ≤ kx − yk , ∀x, y ∈ A} where k·k is the Euclidean norm. Then, the Wasserstein-1 distance between two intensities λ and λ′ is given by Z  ′ ′ W1 (λ, λ ) = sup ϕ(x)(λ(x) − λ (x))dx (11) ϕ∈L1

A

where A λ(x)dx = A λ (x)dx = N . Unlike the L2 metric, W1 accounts for the distance over which spatial mass must be transported and is a more natural way of comparing such equal-mass measures. We now state the key regularity condition on r(·, z). R

R

′

Assumption 1. For every z ∈ A, the function r(·, z) is LLipschitz continuous on S with respect to the W1 metric, i.e., |r(λ, z) − r(λ′ , z)| ≤ LW1 (λ, λ′ ),

λ, λ′ ∈ S.

(12)

The continuity assumption is natural for many performance metrics such as coverage and network capacity. Assumption 1 will later be verified for the coverage probability of a downlink cellular network in which the BSs follow an IPPP. Assumption 1, together with Jensen’s inequality for the absolute value, implies that R is also L-Lipschitz continuous with respect to W1 since |R(λ) − R(λ′ )| = |E[r(λ, z) − r(λ′ , z)]| ≤ E[|r(λ, z) − r(λ′ , z)|] ≤ LW1 (λ, λ′ ).

(13) (14) (15)

4

We also have the following lemma. Lemma 1. Under Assumption 1, and if the partitions satisfy r A , (16) max sup kx − xk k ≤ cg 1≤k≤D x∈Ak D where xk ∈ Ak is the centroid of Ak and cg is a geometric constant independent of D, we have the bound r A Ed ≤ 2cg N L . (17) D Proof: Using (13)-(15) in (10), we obtain Ed ≤ LW1 (λ⋆ , λ⋆d ).

(18)

and therefore it suffices to bound W1 (λ⋆ , λ⋆d ). Let us denote Z D Z X ⋆ ϕ(x)δλ⋆ (x)dx (19) ϕ(x)δλ (x)dx = w(ϕ) := A

⋆

k=1

⋆

where δλ := λ |w(ϕ)| = | ≤

− λ⋆d . Since

D Z X

k=1 Ak D Z X Ak

k=1

(16)

≤ cg

q

A D

Ak

⋆

δλ (x)dx = 0, it holds that (20)

⋆

|ϕ(x) − ϕ(xk )||δλ (x)|dx

k=1

D Z X

k=1

R

(ϕ(x) − ϕ(xk ))δλ⋆ (x)dx|

D Z ϕ ∈ L1 X

≤

Ak

Ak

Ak

(21)

kx − xk k |δλ⋆ (x)|dx

|λ⋆ (x) − λ⋆d (x)|dx ≤ 2cg N

(22) q

A D,

(23)

where R (23) follows from P R the triangle inequality and the fact that A λ⋆ (x)dx = k Ak λ⋆ (x)dx = N and likewise for λ⋆d . Since this bound holds for every ϕ ∈ L1 , it follows that q A W1 (λ⋆ , λ⋆d ) = sup w(ϕ) ≤ sup |w(ϕ)| ≤ 2cg N D . (24) ϕ∈L1

ϕ∈L1

Combining this with (18) gives the required bound. The geometric condition in Lemma 1 is satisfied by many standard partitions. For instance, if the partition is a regular polygon with at least 3 sides, the condition holds with cg = 1. More generally, the condition holds for any shape-regular familyof partitions whose maximum centroidal radius scales as  p A/D . It also accommodates nonuniform and adaptive O meshes, such as those considered in [19], provided that their maximum centroidal radius satisfies the stated bound. Intuition for the discretization error bound comes from viewing the replacement of λ⋆ by λ⋆d as a mass transport problem. The piecewise constant approximation only redistributes BS mass within each sub-region, and the maximum displacement frompany point within a sub-region to its centroid is bounded by cg A/D. Since the total intensity mass is N , the resulting W1 distance, and hence the discretization error, scales linearly with N and with p the characteristic cell length p A/D, so that Ed = O(N L A/D). The shape of the subregions affects the geometric constant cg , but not the scaling of the bound with D or N .

Finally, the dependence on D can be improved by making additional but unverifiable assumption on λ⋆ as stated in the following remark. Remark 1. If the optimal normalized intensity µ⋆ = N1 λ⋆ R ⋆ has bounded variation, i.e., A k∇µ (x)k dx ≤ B⋆ , then it holds that Ed ≤ LNDAB⋆ , as shown in Appendix A in the supplementary material. However, since λ⋆ is itself unknown, it is not possible to verify if it has bounded variation or to obtain B⋆ . Therefore, this improved rate cannot be used as an a priori bound without additional information about λ⋆ . B. Bound on the sampling error Since λ̃d maximizes R over Λ and λ̂d ∈ Λ, it follows that R(λ̃d ) − R(λ̂d ) ≥ 0 so that Es = E[R(λ̃d ) − R(λ̂d )].

(25)

The term within the expectation can be decomposed as:     R(λ̃d ) − R(λ̂d ) = R(λ̃d ) − R̂(λ̃d ) + R̂(λ̃d ) − R̂(λ̂d )   + R̂(λ̂d ) − R(λ̂d ) . (26)

Since λ̂d = arg maxλd ∈Λ R̂(λd ) and λ̃d ∈ Λ, it follows that R̂(λ̂d ) ≥ R̂(λ̃d ) and we can drop the second term on the right. Therefore, the difference can be bounded as   R(λ̃d ) − R(λ̂d ) ≤ max R(λd ) − R̂(λd ) λd ∈Λ   + max R̂(λd ) − R(λd ) . (27) λd ∈Λ

Observe that the two terms in (27) measure, in either direction, the worst-case difference between the exact objective R and its empirical estimate R̂ over all λd ∈ Λ. To bound these terms, we use the Rademacher complexity of the class of performance metrics {r(λd , ·) | λd ∈ Λ}, defined as M

1 X εm r(λd , zm )] λd ∈Λ M m=1

RM := E[ max

(28)

where {εm } are i.i.d. random variables drawn from the Rademacher distribution, i.e., P(εm = 1) = P(εm = −1) = 1/2. The expectation is with respect to both {εm } and {zm }. In the context of learning, the Rademacher complexity measures how well a particular hypothesis class can fit to random labels. In the present setting, the relevance of RM is made precise by the following lemma. Lemma 2. It holds that   E[ max R(λd ) − R̂(λd ) ] ≤ 2RM λd ∈Λ   E[ max R̂(λd ) − R(λd ) ] ≤ 2RM λd ∈Λ

(29) (30)

where the expectations are with respect to the random user locations {zm } appearing in the definition of R̂.

The proof of Lemma 2 is provided in Appendix B of the supplementary material and follows along the lines of similar results in statistical learning theory; see e.g. [20, Sec. 4.5]. A small value of RM implies a small expected uniform deviation

5

between R̂(λd ) and R(λd ) over λd ∈ Λ. Using Lemma 2 in (25) and (27), we obtain   Es ≤ E[ max R(λd ) − R̂(λd ) ] λd ∈Λ   Lemma 2 + E[ max R̂(λd ) − R(λd ) ] ≤ 4RM . (31) λd ∈Λ

It therefore suffices to bound RM . Since Assumption 1 characterizes the sensitivity of r in terms of the W1 distance, whereas RM is indexed by the finite-dimensional vector λd , we first relate changes in r to the Euclidean distance between discretized intensity vectors. Lemma 3. It follows from Assumption 1 that LAd A kλd − λ′d k2 |r(λd , z) − r(λ′d , z)| ≤ √ D for all λd , λ′d ∈ Λ and z ∈ A. Proof: From (11), we have: Z W1 (λd , λ′d ) = sup ϕ(x)(λd (x) − λ′d (x))dx ϕ∈L1 A Z = sup (ϕ(x) − ϕ(x′ ))(λd (x) − λ′d (x))dx ϕ∈L1

(32)

(34)

R R since A λd (x)dx = A λ′d (x)dx for any fixed x′ ∈ A. Since ϕ ∈ L1 , we obtain Z ′ kx − x′ k |λd (x) − λ′d (x)|dx (35) W1 (λd , λd ) ≤ ≤ diam(A)

D Z X

k=1

Ak

|λdk − λ′dk |dx

A kλd − λ′d k1 ≤ Ad √AD kλd − λ′d k2 = Ad D

Substituting into (31), we obtain the desired result. In summary, the bound on the sampling error is   p Õ D/M where the Õ notation hides the logarithmic dependencies  √  as well as other constants. In this bound, the O 1/ M dependence reflects the improvement obtained from averaging over additional user samples. Interestingly, increasing D has two opposing effects:  the Euclidean Lipschitz √ constant Ld decreases as O 1/ D while the Gaussian complexity term associated with Λ increases as Õ (D). Thus, finer discretization reduces the discretization error but requires more user samples to bound the sampling error. C. Bound on excess error

(33)

A

A

PD P ND Since j ηj λdj ≤ (maxj ηj ) j=1 λdj = A maxj ηj , it follows from [22, Lemma 2.3] that q q p πD log(D) π ND . 2M log(D) = N LA RM ≤ LAMd A 2D d A M

(36) (37)

which yields the required Lipschitz √ constant. Accordingly, let Ld := LAd A/ D denote the resulting Euclidean Lipschitz constant. A key ingredient for bounding RM is the so-called vector contraction bound, first established in [21], which is re-stated below using the current notation:   r M X D X π Ld  (38) gij λdj  E max RM ≤ λd ∈Λ 2M i=1 j=1

where gij are i.i.d. standard normal random variables. The vector contraction inequality is stated in√[21] using Rademacher random variables, with the constant 2. The same argument applied to standard normal random variables yields (38), with p the constant π/2, as shown in the proof of [21, Proposition 1]. The following lemma gives the resulting sampling-error bound. Lemma 4. Under Assumption 1, we have the bound r πD log(D) Es ≤ 4N LAd (39) M Proof: Substituting Ld from Lemma 3 and using the fact P M that i=1 gij = ηj ∼ N (0, M ) in (38), we obtain   r D X π  LAd A (40) ηj λdj  E max RM ≤ λd ∈Λ M 2D j=1

Combining the bounds in Lemmas 1 and 4, we obtain the total bound on the excess error for square cells as  q q  q πD log(D) 1 D √ . + 4N LA = Õ + E ≤ N L 2A d D M M D where the Õ notation hides the constants including N , A, and Ad . As with most worst-case generalization bounds, the expression is informative when its right-hand side is small relative to the range of the performance metric. In particular, the bound reveals how the excess error scales with D, M , and N , and clarifies the tradeoff between its discretization and sampling components. Roughly, the error bound would be minimum if the two terms match, i.e., when D2 ≈ M . Hence, for fixed N , A, and Ad , if we want to ensure that the
excess error is at
most ǫ, then we should choose D = Õ ǫ12 and M = Õ ǫ14 , ignoring the logarithmic terms. We refer to these bounds as discretization and sampling complexities of the proposed scheme, respectively. If the underlying intensity has bounded variation
(cf. Remark 1), then the complexities
improve to D = Õ 1ǫ and M = Õ ǫ13 . Thus, the proposed formulation learns a spatial BS intensity directly from user-location samples without requiring prior knowledge or parametric estimation of Pz , while the bound quantifies the number of samples needed to support a chosen spatial resolution. IV. O PTIMAL D EPLOYMENT OF BS S We now present the optimization approach to solve (P̂) and discuss how the solution provides soft guidelines along with their connection to the real-world BS-deployment. A. Soft Deployment: Solving (P̂)

Since (P̂) is high-dimensional, we obtain the softdeployment using only first-order methods. Projected gradient methods are less attractive because maintaining feasibility at every iteration entails solving a Euclidean projection subproblem over the scaled simplex Λ, typically using a dedicated

6

sorting- or active-set-based procedure [23]. Instead, we propose to use the mirror ascent algorithm that yields a closedform update. Starting at an arbitrary λ(0) in the relative interior (0) of Λ (such as the point λdi = N/A for all i), the mirror ascent updates at the t-th iteration are evaluated as 1 (t) (t+1) (t) λd = arg maxh∇R̂(λd ), λi − (t) DKL (λ, λd ) (41) λ∈Λ η (t)

η (t) is the step-size and DKL (λ, λd ) := (t) i=1 λi log(λi /λdi ) is the KL divergence. These updates can be written in closed form as in [24]. The empirical form of R̂ further motivates the use of the stochastic mirror ascent (SMA) algorithm:   (t) (t) (t) λ exp η g i di ND (t+1)   λdi = (42) A PD λ(t) exp η (t) g (t)

where PD

j=1

dj

j

(t) for 1 ≤ i ≤ D, where {gi } are entries of

g(t) :=

X 1 (t) ∇r(λd , zm ) |Mt |

(43)

m∈Mt

with Mt being a random mini-batch of users at iteration t. After running for T iterations, the optimization algorithm may yield a point λ̄Td . Without additional global structure, such as concavity of R̂, it is not possible to find a bound on the optimization error Eo = E[R̂(λ̂d ) − R̂(λ̄Td )] that decreases to zero with T . If, however, we assume that r(·, z) is relatively smooth, which follows if k∇r(λd , z) − ∇r(λ′d , z)k∞ ≤ L∇ kλd − λ′d k1 , then existing results for nonconvex stochastic mirror ascent guarantee convergence to a stationary point of (P̂) [24]. Further, Appendix F of the supplementary material shows that under the noise-limited IPPP approximation, the coverage probability is concave and yields the bound Eo ≤ q D from [25]. N 2 log T Finally, since the empirical problem is generally nonconcave, the stationary point reached by mirror ascent may depend on the initialization. We initialize the BS intensity using the empirical user-density map, since it provides a data-driven starting point that already allocates greater deployment mass to regions of high user concentration. B. From Guidelines to BS Deployment the optimized intensity function λ̄T (x) := T k=1 λ̄dk 1 x∈Ak serves as a guideline for BS deployment, we use it to generate BS deployment by sampling the underlying spatial point process with optimized intensity. Here, each realization serves as a potential deployment candidate with different realized performance. To further optimize the deployment, we generate several candidates and select the one with the highest estimated performance using an independent set of user-location samples. For a realization φ of the spatial point process Φλ , let r(φ, z) denote, with a slight abuse of notation, the performance of the corresponding BS deployment for a user at z, and define R(φ) := Ez [r(φ, z)]. We then have

PAs D

R(λ) = Eφ [R(φ)] = Eφ Ez [r(φ, z)].

(44)

where the outer expectation is over realizations φ of Φλ . We draw K independent candidate deployments {φk }K k=1 from the point process Φλ̄T . Note that the original userlocation samples Z that were used to obtain λ̄Td should not also be used to select the final realization. Since λ̄Td was chosen specifically to perform well at these locations, reusing them may favor a deployment that performs well on Z but performs poorly for the underlying user population. We therefore draw an independent hold-out validation set of user v locations Zv := {zvm }M m=1 , which is used only to select the final deployment φ̂K = φk̂ , where Mv 1 X k̂ = arg max R̂(φk , Zv ) := r(φk , zvm ). 1≤k≤K Mv m=1

(45)

In the next subsection, we show that for sufficiently large values of K and Mv , the selected deployment cannot be much worse than the average performance R(λ̄T ) of the optimized intensity function. 1) Placement complexity: Define the placement error as i h (46) Ep = E (R(λ̄T ) − R(φ̂K ))+

where the positive part (a)+ = max{0, a} is used to exclude the realizations that perform better than the average. The following lemma quantifies the additional error incurred in obtaining a concrete deployment. Lemma 5. If r(φ, z) ∈ [0, Br ] for all φ and z, then the placement error satisfies q Ep ≤ BKr + Br log(K) (47) 2Mv Proof: The placement error can be split as

Ep ≤E[(R(λ̄T ) − max R(φk ))+ ] + E[max R(φk ) − R(φ̂K )] k

k

where the two terms capture, respectively, the errors due to generating only K candidate deployments and selecting φ̂K in (45) using the validation estimate R̂ instead of R. For the first term, observe that conditioned on λ̄T , the random variables {R(φk )}K k=1 are i.i.d. in [0, Br ] with mean Eφk R(φk ) = R(λ̄T ). Hence, using Markov’s inequality on the random variable Br − R(φk ) we obtain the required bound:   T (48) E R(λ̄ ) − max R(φk ) ≤ BKr k

+

which is proved in Appendix C of the supplementary material. In the second term, conditioned on the candidate deployments {φk }K k=1 , each R̂(φk , Zv ) is the average of Mv i.i.d. random variables taking values in [0, Br ]. Therefore, Hoeffding’s lemma implies that ∆R̂k := R̂(φk , Zv ) − EZv [R̂(φk , Zv )] is sub-Gaussian with parameter Br2 /(4Mv ) [20, Sec. 1.2.1]. Hence, an expectation-of-the-maximum bound for sub-Gaussian random variables yields [20, Prop. 1.5]: q K E[max ∆R̂k ] ≤ Br log (49) 2Mv . k

Let k ⋆ = arg maxk R(φk ) so that maxk R(φk ) = R(φk⋆ ). Since φ̂K maximizes R̂, we have

R(φk⋆ ) − R(φ̂K ) ≤ R(φk⋆ ) − R̂(φk⋆ , Zv ) + R̂(φ̂K , Zv ) − R(φ̂K ).

7

Since k ⋆ depends only on {φk }, the first term on the right has zero mean. The second term is bounded by E[maxk ∆R̂k ]. Therefore from (49), we obtain r   log K E max R(φk ) − R(φ̂K ) ≤ Br (50) 1≤k≤K 2Mv which, together with (48), yields the desired bound. Lemma 5 shows that Ep ≤ ǫ can be achieved
with K = O (1/ǫ) candidate deployments and Mv = Õ 1/ǫ2 validation samples. Since sampling candidate deployments from Φλ̄T is usually inexpensive, the validation sample complexity is generally the more demanding requirement. Algorithm 1 End-to-End Base Station Deployment 1: Partition A into D subregions and construct Λ 2: Obtain M user-location samples Z 3: Construct the empirical objective R̂ 4: Run updates in (42) for T iterations to obtain λ̄T d 5: Generate K candidate deployments {φk }K k=1 from Φλ̄T 6: Obtain Mv validation user-location samples Zv 7: Solve (45) and return φ̂K 2) End-to-end performance guarantee: Algorithm 1 summarizes the complete procedure for obtaining a concrete deployment from user-location samples. The following theorem combines the earlier bounds on discretization, sampling, optimization, and placement. Theorem 1. Under the conditions of Lemmas 1, 4, and 5, the deployment returned by Algorithm 1 satisfies    ⋆ E R(λ ) − R(φ̂K ) + q q q  A D 1 1 . (51) + + + ≤ Eo + Õ D M K Mv The result follows from the decomposition   R(λ⋆ ) − R(φ̂K ) ≤ R(λ⋆ ) − R(λ̃d ) + R(λ̃d ) − R(λ̄Td ) +   + R(λ̄T ) − R(φ̂K ) (52) +

where the first two differences are nonnegative by the optimality of λ⋆ and λ̃d , respectively. Applying the decomposition in (26) with λ̄Td in place of λ̂d , and using R̂(λ̃d ) ≤ R̂(λ̂d ), gives E[R(λ̃d ) − R(λ̄Td )] ≤ 4RM + Eo . Taking expectations and applying Lemmas 1, 4, and 5, together with the definition of Eo , proves (51). Except for the numerical optimization error Eo , the other terms can small by setting
be made arbitrarily
D = Mv = Õ 1/ǫ2 , M = Õ 1/ǫ4 , and K = Õ (1/ǫ). As stated earlier, Eo can be made small in certain cases, such as when the deployment process is IPPP and the system is noise-limited. Beyond BS deployment, Theorem 1 applies to intensitybased spatial design whenever the point-process model and performance metric satisfy the stated regularity conditions. To the best of our knowledge, it provides the first endto-end nonasymptotic error decomposition that incorporates discretization, finite-sample learning, and realization selection components, up to the numerical optimization error Eo .

V. O PTIMAL D EPLOYMENT UNDER S PECIFIC P ERFORMANCE M ETRICS We now specialize the framework to an inhomogeneous Poisson point process (PPP) Φλ on A with intensity λ. In particular, we consider performance of a downlink cellular network at conventional sub-6 GHz frequencies under two cell association rules: instantaneous signal-to-interferenceplus-noise ratio (SINR)-based association and average SINRbased association. We consider an orthogonal multiple-access downlink system model, where the BSs transmit constantly with a fixed power P to a single desired mobile user on any particular time-frequency resource. As in [26], we consider a generic path-loss function R ∞ ℓ(r) that satisfies (a) 0 < ℓ(r) ≤ 1, (b) ℓ′ (r) ≤ 0, and (c) 0 ℓ(r)rdr < ∞, e.g., ℓ(r) = 1+r1α /c where α ∈ (2, 5] is the path-loss exponent. For our analysis, we additionally require that the path-loss function does not decay faster than an exponential, i.e., there exists Lℓ < ∞ such that ℓ(r)eLℓ r becomes an increasing function of r. In practice, path-loss is typically modeled using a power law rather than exponential or super-exponential functions, so Lℓ tends to be small across this class of models. The average received power at the user is given by P ℓ(kx − zk) and the corresponding average SNR is P ℓ(kx − zk)/σ 2 where σ 2 is the noise power. The interference power received by a user located at z is denoted by I(z). We consider Rayleigh fading, where the fading gain for each link is i.i.d. exponential with unit mean, in which case, the Laplace transform of I(z) is known to be [27, Sec. III-C]:  Z  λ(y)dy LI(z) (s) = exp − (53) 1 R 1+ sP ℓ(ky−zk)

where R is the (sub-)region in A where the interfering BSs lie. As we shall see later, the expression allows us to calculate the expectation of various functions of I(z). In particular, the signal-to-interference-plus-noise ratio (SINR) [28] at the user located at z from a BS located at x is given by ℓ(kx−zk) . SINRx = GPI(z)+σ 2

(54)

The mobile user located at z ∈ A is said to be in coverage if the BS (located at x) selected under the association rule provides an SINR at z exceeding the threshold τ . We are interested in the performance metric coverage probability which captures the probability that a typical user is covered, i.e., pc (λ) = Ez [pc (λ, z)] = P[SINR(z) ≥ τ ]. A. Instantaneous SINR-Based Association In this association model, each mobile user connects to the BS offering the highest instantaneous SINR, considering both fading and path-loss. Interestingly, a user may not always be associated with the closest BS, since a higher instantaneous SINR may be received from a different BS. Let us denote the SINR for the user at z when it connects to a BS at x by SINRx (z), or simply by SINRx . Hence, the probability that a user at z is covered is given by   pc (λ, z) = P max SINRx ≥ τ (55) x∈Φλ

8

where the probability is calculated with respect to the random fading gains and the random BS point process. In order to simplify (55), we assume that τ > 1 so that for a given user, the instantaneous SINR can exceed τ for at most a single BS, as in [27]. Subsequent applications of the Campbell-Mecke as well as Slivnyak’s theorems [29] yield Z pc (λ, z) = P [SINRx > τ ] λ(x)dx. (56) A

Since G ∼ Exp(1), we have that

P [SINRx > τ ] = EI(z) [P [SINRx > τ | I(z)]]    σ 2 + I(z) = EI(z) exp −τ P ℓ(kx − zk)     τ τ σ2 LI(z) = exp − P ℓ(kx − zk) P ℓ(kx − zk)

(57) (58) (59)

where the expression for LI(z) (s) is given in (53) with R = A. Substituting (59) into (56), we obtain the required coverage probability expression Z R ℓ(kx−zk) τ σ2 λ(x)e− P ℓ(kx−zk) − A λ(y)/(1+ τ ℓ(ky−zk) )dy dx pc (λ, z) = ZA λ(x)h(λ, x)dx (60) = A

where we have suppressed the dependence on z and defined   Z h(λ, x) := exp −a(x) − λ(y)b(x, y)dy (61) A

κ a(x) := ℓ(dx )

b(x, y) :=

τ ℓ(dy ) . τ ℓ(dy ) + ℓ(dx )

(62)

2

with dx := kx − zk, dy := ky − zk, and κ := τ σ /P . The following lemma, whose proof is provided in Appendix D in the supplementary, verifies Assumption 1 for (60). Lemma 6. The coverage probability in (60) under instantaneous SINR association rule is 23 Lℓ -Lipschitz continuous with respect to the W1 metric. As an example, consider the path-loss ℓ(r) = 1+r1α /c for some c > 0. For this case, it can be seen that 1/α  ′ α−1 (r)| sup |ℓℓ(r) ≤ (α−1) =: Lℓ (63) c r≥0

so that pc (λ, z) is 23 Lℓ -Lipschitz. For c = 1 and α ∈ (2, 5], this quantity is upper bounded by approximately 3. The expression in (60) can be discretized as demonstrated earlier. 1) Extension to multi-tier networks: The current formulation also applies to a heterogeneous cellular network with Q overlaid independent IPPP tiers of BSs, each of which may differ in terms of transmit power, deployment density, target SINR, etc. Denoting the intensity of the BSs in the q-th tier by λ[q] (x), the coverage probability, with instantaneous SINRbased association, may similarly be derived as [27]: Q Z τi σ 2 X − P ℓ(kx−zk) −wi (x,z) [i] [q] i λ (x)e pc ({λ }, z) = dx wi (x, z) =

i=1 A Q Z X j=1

λ[j] (y)dy

Pi ℓ(kx−zk) A 1 + τ P ℓ(ky−zk) i

j

where Pi and τi > 1 are the powers and SINR thresholds for the i-th tier. Hence, the goal now is to maximize E[pc ({λ[q] }, z)] under the constraints under the prescribed mean BS counts for each tier. Discretizing these densities yields a QD-dimensional optimization problem. B. Average SINR-Based Association In this association model, the mobile user connects to the BS that offers the highest average SINR. Under the equalpower and monotone path-loss model considered here, this rule translates to associating each user with the nearest BS. For a user located at z, let xz be the nearest BS and Rz := kxz − zk be the serving distance. In case the IPPP realization does not contain any BS, we set Rz = ∞ and the corresponding event contributes zero coverage. Then, the coverage probability is the probability that the SINR at xz is greater than the threshold τ , i.e., pac (λ, z) = P [SINRxz > τ ] Z ∞ = P [SINRxz > τ | Rz = r] fRz (r)dr

(64) (65)

0

where fRz (r) is the density of Rz for finite r. To calculate fRz , recall that for a PPP, the probability that no BS lies within the set B(z, r) = {x ∈ A | kx − zk ≤ r} is given by ! Z P [Rz > r] = exp − λ(y)dy . (66) B(z,r)

Therefore, for finite r, the density of Rz is given by d P [Rz > r], which evaluates to − dr ! "Z # Z d fRz (r) = exp − λ(y)dy λ(y)dy . (67) dr B(z,r) B(z,r) The derivative in (67) can be evaluated analytically for simple regions or numerically for a general deployment region A. GP ℓ(r) Next, given Rz = r, we have that SINRxz = I(z)+σ 2. Hence the conditional probability becomes P [SINRxz > τ | Rz = r] = EI(z) [P [SINRxz > τ | Rz = r, I(z)]] i h  2 (68) = EI(z) exp −τ σP+I(z) ℓ(r)     2 τ σ LI(z) P ℓ(r) (69) = exp − Pτℓ(r)

where the expression for LI(z) (s) is given in (53) with R = A \ B(z, r), since only BSs farther than the serving BS contribute interference. Hence, (69) becomes P [SINRxz > τ | Rz = r] Z   τ σ2 λ(y)dy = exp − − ℓ(r) P ℓ(r) A\B(z,r) 1+ τ ℓ(ky−zk)

(70)

which can be substituted into (65) to obtain the coverage probability. We remark that for other fading distributions as well, the distance-based association rule remains unchanged, and only the conditional probability in (69) changes. Moreover, pac (λ, z) is Lℓ /e-Lipschitz with respect to W1 , as proved in Appendix D of the supplementary material.

9

C. Practical Applications We now present several additional practical applications of the proposed soft deployment framework, where the derived guidelines can be adapted to more realistic environments. 1) Additional deployment constraints: An important feature of the current framework is that it is relatively straightforward to incorporate physical constraints into the formulation. We discuss some common deployment constraints. It is common to have certain areas to be off-limits for BS deployments, e.g., due to regulations or geographic limitations. Such exclusion zones may be enforced by setting λdk = 0 whenever Ak is part of an exclusion zone. More generally, if the average number of BSs within a region B ⊂ A consisting of complete grid cells is to be restricted to at most NB < N , we can include the constraint X D (71) λdk ≤ NB A . k:Ak ⊂B

In contrast, certain high-priority regions may necessitate minimum deployment intensity requirements, which can be imposed by requiring λdk ≥ λmin . Finally, in multi-tier deployments, we might specifically impose constraints of the PQ [q] form q=1 λdk ≤ D A so as to ensure that each sub-region gets, on average, at most one BS of any tier. 2) Incremental deployment: Incremental deployment refers to the scenario in which a deterministic set Φ0 of N0 BSs has already been installed and the objective is to optimize the intensity of an additional BS deployment having mean cardinality N . The new BS intensity λ can help identify regions where further infrastructure is most beneficial given the coverage and interference provided by Φ0 . The coverage expressions are modified to account for the superposition Φ0 ∪ Φλ so that the fixed BSs participate in association and appear explicitly as interferers in the SINR expressions, while their locations remain unchanged during the optimization. The corresponding IA and AA coverage expressions are provided in Appendix E of the supplementary material. 3) Adaptive deployment: The proposed framework separates the deployment problem into two stages: soft deployment and concrete deployment, allowing additional site-specific constraints to be accommodated during the second stage. For instance, detailed site availability or deployment restrictions may only become known after the intensity has been optimized. If a small portion of a selected grid cell is unavailable, the BS can still be placed within the feasible part of that cell. If an entire cell becomes unavailable, its intensity mass may be redistributed among nearby feasible cells, allowing the deployment to be adjusted without solving the original optimization problem again. 4) Sustainable deployment: Green 6G targets can be imposed within the framework by attaching a cost density c(x) to every candidate site. The cost density may model locationspecific carbon footprint, a proxy for electromagnetic exposure, or backhaulR rental costs. Defining the expected aggregate cost C(λ) := A c(x)λ(x)dx, we can (a) include a softpenalty and maximize R(λ) − βC(λ) to explore the coveragecost tradeoff, or (b) impose a budget on the expected cost of the form C(λ) ≤ Cmax , which becomes a linear constraint after discretization.

5) Changing user demand: The optimized intensity naturally assigns different deployment priorities to different regions. Examples include high-demand regions such as malls and stadiums. The framework can also be naturally extended to scenarios where the user distribution changes over time, by updating the user samples and re-optimizing the intensity. VI. N UMERICAL R ESULTS In this section, we evaluate the proposed approach with the instantaneous coverage probability derived in Sec. V-A. A comparison of instantaneuous and average association scenarios is provided in Appendix H of the supplementary material. We begin by visualizing the optimized intensity under different user distributions and deployment constraints (Sec. VI-A). Next, we compare the performance of the proposed and baseline deployment approaches, and validate the theoretically predicted effects of discretization, finite number of userlocation samples, and validation-based deployment selection (Sec. VI-B). Finally, we apply the proposed framework to find an optimal incremental BS deployment on top of an existing cellular deployment in the city of Chicago. TABLE I C OMMON NUMERICAL PARAMETERS . S CENARIO - SPECIFIC VALUES ARE STATED IN THE CORRESPONDING SUBSECTIONS . Parameter Value Deployment area A 10 × 10 km2 Path-loss exponent α 3.5 SINR threshold τ 2 Transmit power P 1 Path-loss scale c (0.2 km)α Noise power σ2 10−4

Parameter Value Number of cells D 1600 iterations T 100 Training samples M 1600 Validation samples Mv 400 Test samples Mt 10000 Candidates K 1000

All experiments use the default parameters in Table I unless otherwise noted. Along with the unconstrained (greenfield) deployments, we will also consider deployments under two types of deployment constraints: exclusion zones and incremental deployments with existing BSs at known locations. The code is provided in the supplementary material. The mirror-ascent updates in (41) are initialized at λ(0) = 0.99λu + 0.01λc where λu and λc refer to the demandproportional and uniform intensity functions, respectively. The small uniform component ensures a strictly positive initialization on the relative interior of the feasible simplex, while the demand-proportional component serves as a data-dependent warm start. For experiments involving concrete deployments, the training set is used to learn the intensity while an independent validation set is used to select a deployment, whose final performance is evaluated on a test set. A. Optimized Intensity Maps: Guidelines and Trends We begin with several examples with different scenario/setups to demonstrate the application of the proposed soft deployment approach. The solutions of these examples jointly illustrate how the optimized intensity changes with the user distribution and deployment constraints, providing important insights. In all figures, the color scale represents the spatial distribution of the BS intensity while the training user locations are represented by red dots.

10

Users

2

0.2

change with the installed network. Interestingly, the lowintensity ridge that connects the bottom two clusters suggests that they can also be served by the intermediate locations.

1 5

0.15

Users Existing BS

0

2 1.8

0.1

1.6 1.4

-1 0.05

1.2

0

1

-2

0.8

-2

-1

0

1

0

0.6

2.5

0.4

2

4 Users Existing BS

3

0.2

2

2 1

1.5

0

0

0

-5

0

5

Fig. 2. Optimized BS intensity for a scenario with mixed urban demand under two different planning contexts: (a) unconstrained and (b) incremental with an existing BS and an exclusion zone.

1

-1 -2

0.5

-3 -4 -4

-5 -5

0

-2

0

2

4

Fig. 1. Optimized BS intensities for representative user distributions. (top) Annular user-location distribution with mean N = 1. (bottom) Six-cluster user-location distribution with one already installed BS and mean N = 5 additional BSs.

We first consider a scenario with an annular user demand with N = 1 as shown in Figure 1(top). Observe that the optimized intensity does not simply follow the user density, but instead forms a concentric annulus of smaller radius. The inward displacement reflects a trade-off between remaining close to the individual users and providing comparable service around the entire ring. In additional runs, see Appendix G of the supplementary material, we observed that the optimized annulus initially contracts towards the center but again expands towards the edges, as noise level decreases. In contrast, a hard single-BS deterministic optimizer places the BS at the center in the low-noise regime and switches to an offcenter location near the user annulus as the noise increases, thereby masking the gradual radial transitions revealed by the optimized intensity. Further, due to radial symmetry of user distribution, for deterministic optimizer, any other location at the same distance from the center but at a different angle is also optimal. The optimizer picks one of these arbitrarily, again masking the equivalence of these points. Instead the soft deployment approach reveals the equivalence of points providing an additional flexibility to let us choose any BS location from these equivalent options. We next consider a scenario where the users form six spatial clusters, one of which is already served by an installed BS, over an 8 × 8 km2 area, as shown in Figure 1(bottom). As expected, for N = 5, the optimized incremental-BS intensity is only concentrated primarily at the remaining five clusters. Thus, the intensity map also reveals how deployment priorities

To further highlight that the optimized BS intensity represents deployment utility rather than just user density, we consider a mixed urban demand pattern under two different planning contexts. The user distribution is kept identical across the two cases and consists of a combination of localized demand regions, an annular venue, and a corridor. Fig. 2(a) considers a greenfield (unconstrained) deployment, whereas Fig. 2(b) introduces an incremental deployment over an existing BS as well as a circular exclusion zone. With N = 5, the resulting intensity maps differ substantially: the existing BS reduces the value of placing additional mass in the region it serves and the exclusion zone forces the intensity to concentrate along its outer edge. Thus, the optimized intensity captures not only where additional infrastructure is useful, but also how deployment opportunities shift in response to existing infrastructure and spatial constraints. B. Deployment performance and finite-sample behavior We next evaluate the quality of the concrete BS deployments generated from the optimized intensity and compare them with four baseline deployment strategies: two randomized and two deterministic. The demand-proportional baseline assigns the BS intensity λu in proportion to the empirical (training) user count in each feasible cell, while the uniform intensity baseline assigns constant BS intensity λc to each feasible cell. Both these randomized approaches subsequently follow the same validation-based deployment procedure as outlined in Sec. IV-B. The greedy baseline sequentially adds one BS at the feasible candidate location that maximizes training coverage probability, with all previously selected (and existing) BSs contributing to the SINR calculation. Finally, in K-means clustering, the training users are divided into N clusters and the BSs are placed at the resulting cluster centroids. In case the cluster centroids lie inside the exclusion zones, they are projected to the nearest point in the feasible region. We note that the optimal deterministic BS deployment can, in theory, be found by maximizing the training coverage probability over the set of feasible BS locations. However, the resulting problem has a non-convex non-smooth objective and

11

We consider two metrics for evaluating the randomized methods: (a) the conditional complementary cumulative distribution 1 PK 1 {p̂ function (CCDF) given by K c (φk , Zt ) ≥ ρ} which k=1 gives the fraction of N -BS deployments whose coverage probability is at least ρ; and (b) the best deployment test coverage probability p̂sel c = p̂c (φ̂K , Zt ) where the deployment is selected using validation set as in (45). A higher CCDF indicates a higher probability of generating a high-quality deployment, while a higher p̂sel indicates that the selected c deployment is indeed high-quality. As can be seen from Fig. 3, the proposed approach yields a markedly more favorable distribution of deployment quality, as compared to the other randomized schemes, for both N = 2 and N = 4 cases. We also observe that for the N = 2 case, the best validation-selected deployments from the demand-proportional, uniform, and optimized intensities have similar coverage probabilities. This is because in the scarceBS regime, only a few good-quality deployments possible, and sufficiently many candidate draws can yield a good deployment from any starting λ. However, as evident from the conditional CCDF curves, a large fraction of deployments resulting from the optimized λ are high-quality. The separation is clearer in the N = 4 case, where the validation-selected deployment of the optimized λ has very high coverage. As expected, the demand-proportional intensity often yields poor-

quality deployments, as it does not account for the coverage and interference of the installed network. The same is true for the deterministic K-means algorithm, which also achieves poor performance in such settings. Interestingly, the greedy placement overcomes both these issues, but since it selects the best BS locations in a sequential manner, the placement of its first BS is often suboptimal. Further, such deployments may overfit a specific scenario, resulting in a lower performance on the test set.

Conditional CCDF

0.8

0.6

0.4

0.2

0.8

Conditional CCDF

non-convex constraints, making it difficult to solve. Further, such deterministic methods return a single deployment and do not directly characterize alternative viable configurations, making them less flexible under operational uncertainty. On the other hand, for randomized methods, alternative highquality deployments can be obtained from the same optimized intensity without re-solving the problem. For these experiments, we consider a synthetic user-location distribution, where 30% of the users are uniformly distributed over an annulus of outer radius 1.7 and thickness 0.7, centered at (−1, −1) while the remaining users follow a Gaussian distribution centered at (3, 3). An existing BS is located at the center of the Gaussian component. We consider two cases: scarce deployment with N = 2 and an adequate deployment with N = 4. Recall that the candidate deployments {φk }K k=1 generated by the randomized approaches, namely those involving uniform, demand-proportional, and optimized intensities, have N BSs on average. However, a randomly generated deployment φk could contain more or fewer than N BSs, making its coverage probability unsuitable for a direct comparison with the deterministic, namely greedy or K-means, methods that place exactly N BSs. For the purpose of comparing concrete deployments under the same physical BS budget, we therefore generate deployments {φk }K k=1 that have exactly N BSs. This conditioning is used only for evaluation in Fig. 3 and the proposed intensity is still obtained by optimizing the original IPPP objective. Let us denote the test set of user locations by t Zt := {ztm }M m=1 and the empirical coverage probability of a |Z| deployment φ over the user set Z = {zm }m=1 by 1 X|Z| pc (φ, zm ). (72) p̂c (φ, Z) = m=1 |Z|

0.6

0.4

0.2

0 0.4

0.5

0.6

0.7

0.8

0.9

1

Test coverage probability

Fig. 3. Empirical distribution of pc evaluated over the common test set, conditioned on the deployments having exactly N BSs. The legend provides test coverage probability of the validation-selected deployment.

We next examine empirically the discretization and sampling errors bounded in Sec. III and illustrate that their dependence on D and M follows the trends predicted by the bounds. All subsequent expressions are evaluated for the original mean-N PPP model. Since the optimal intensities appearing in Ed and Es cannot be computed exactly, we use the following approximations. 1) For each D, a solution of (P̂) obtained using a large M approximates λ̃d . 2) Solution of (P̂) obtained using large M and D approximates λ⋆ . 3) For a given λ, pc (λ) is estimated by averaging pc (λ, z) over a large independent test set. 4) For each M , the finite-sample experiment is repeated using independent training sets, so that the average of the resulting errors approximates Es . 5) Finally, in all cases, the output of the mirror-ascent algorithm is used to approximate the corresponding optimizer. Fig. 4 shows the approximate discretization and sampling p √ M/D, respectively. As errors as functions of D and √ expected, the discretization error decreases with D, since a finer piecewise-constant intensity can better represent a given continuous intensity. Likewise, as expected the sampling error

12

0.12

0.1

0.008

Sampling error

Discretization error

C. Incremental deployment in Chicago

0.01

0.08

0.06

0.04

0.006

0.004

0.002

0.02

0

0 5

10

20

0

2

4

6

Fig. 4. Discretization and sampling errors decrease inversely with p M/D, consistent with the theoretical bounds.

8

√

D and

0.95 0.9

Test coverage probability

0.85 0.8 0.75 0.7 0.65

Finally, we consider incremental deployment over Chicago. For the user distribution, we utilize Facebook high resolution population density maps, which estimate resident population on approximately 30 m grid cells from satellite imagery and census data [30]. The resulting population map is used as a stand-in for the user-location distribution and deployment region, with cells outside its convex hull or predominantly within large water, wetland, or protected open-space regions treated as exclusion zones. Further, for existing infrastructure, we utilize AT&T cell records obtained from OpenCellID [31]. Nearby records belonging to the same carrier are merged into spatial site proxies. After geographic and clustering filters, the experiment contains 56 existing AT&T site proxies, which are treated as existing BSs. The 14,646 user-location records are split into 70% training, 15% validation, and 15% test samples. The existing deployment is highly nonuniform and contains large spatial regions where a limited incremental deployment can be useful. Fig. 6 shows the user demand, existing site proxies, and the optimized BS intensity for mean N = 5. We observe that the optimized intensity is concentrated in regions having both substantial demand and comparatively sparse existing infrastructure. The experiment demonstrates that the proposed framework can be used to obtain insights regarding incremental deployment in existing urban centers.

0.6

VII. C ONCLUSION

0.55 0.5 100

101

102

103

104

Fig. 5. Test coverage probability of selected deployment for N = 4.

p also decreases with M/D, since a finer grid increases the size of the vector being estimated and hence requires more user-location samples. The trends are similar for both values of N . For the sampling error, error bars show 95% confidence intervals over 10 independent user realizations. We also examine the coverage probability of the validationselected deployment as a function of K and Mv . Fig. 5 shows the coverage probability for the selected deployment, estimated by averaging over a large number of independent test user-location samples. For comparison, we also show the coverage probability obtained by oracle selection using the test set, as well as the average coverage probability pc (λ̄T ) of the optimized point process. We observe that as K increases, the coverage improves rapidly at first and then more gradually, reflecting the benefit of generating a larger candidate set. Further, the gain from increasing K is larger for Mv = 200 than for Mv = 10, since a larger validation set more reliably identifies the better deployments. The curves corresponding to different values of Mv are close, suggesting that modest values of K and Mv may be sufficient to obtain a deployment that is better than pc (λ̄) with high probability. The worst-case dependence on K and Mv obtained in Lemma 5 can therefore be conservative in structured deployment problems.

We studied base station (BS) deployment optimization using stochastic geometric principles, where the deployment intensity is optimized using sampled user-location data. The proposed framework converts the infinite-dimensional problem into a finite-dimensional empirical problem that can be solved using standard optimization algorithms. We established bounds on the discretization and sampling errors arising from these approximations, as well as the error incurred in selecting a concrete BS deployment from the optimized intensity. The framework was applied to maximize coverage probability in downlink cellular networks under different association rules and practical deployment constraints. Numerical experiments demonstrate how the optimized intensity captures the effects of user distribution, existing infrastructure, and deployment constraints. More importantly, unlike a single deterministic deployment, the resulting intensity reveals broader spatial deployment priorities and alternative viable deployment regions, providing useful insights for cellular network planning. A PPENDIX A B OUND IN R EMARK 1 The bound follows directly from the application of a wellknown Poincaré–Wirtinger inequality [32, Sec. 5.8.1], which yields: Z Z diam(Ak ) |∇µ⋆ (x)|dx. kµ⋆ (x) − µ⋆d (x)k dx ≤ 2 Ak Ak (73)

13

of the training sample with the same distribution as {zm }M m=1 . Defining M

R̂′ (λd ) :=

1 X r(λd , z′m ) M m=1

(76)

and using E′ to denote the expectation with respect to {z′m }M m=1 , we can write     max R(λd ) − R̂(λd ) = max E′ [R̂′ (λd )] − R̂(λd ) λ ∈Λ λd ∈Λ  d  = max E′ [R̂′ (λd ) − R̂(λd )] λd ∈Λ   ≤ E′ [ max R̂′ (λd ) − R̂(λd ) ] (77) λd ∈Λ

where we have used the convexity of the max function. Next, we use the observation that a symmetric random variable X has the same distribution as εX, where ε is an independent Rademacher random variable. Applying it to X = r(λ, z′m ) − r(λ, zm ), substituting into (77), and taking expectation with respect to {zm }, we obtain   E[E′ [ max R̂′ (λd ) − R̂(λd ) ]] λd ∈Λ

M

1 X εm (r(λd , z′m ) − r(λd , zm ))] λd ∈Λ M m=1

= E[ max

M

1 X εm r(λd , z′m )] λd ∈Λ M m=1

≤ E[ max

M

1 X −εm r(λd , zm )]. λd ∈Λ M m=1

+ E[ max

Fig. 6. Optimized incremental BS intensity in Chicago with existing infrastructure and exclusion regions. The inset shows the zoomed-in version of the dotted rectangle.

where the expectations on the right are with respect to {zm , z′m , εm }. Again, since εm r(λd , z′m ) has the same distribution as −εm r(λd , zm ), the last two terms in (78) are the same. Hence, substituting into (77), we obtain   E[ max R(λd ) − R̂(λd ) ] λd ∈Λ

M

1 X εm r(λd , z′m )] = 2RM λd ∈Λ M m=1

≤ 2E[ max Multiplying both sides by N , summing over k = 1, . . . , D, and q 2A using the bound diam(Ak ) ≤ D for square-sub-regions, we obtain r D Z X A ⋆ ⋆ (74) kλ (x) − λd (x)k dx ≤ N B⋆ 2D Ak k=1

Substituting into (23) and proceeding as in the proof of Lemma 1, we obtain Ed ≤

LN AB⋆ D

(75)

which is the required result. A PPENDIX B P ROOF OF L EMMA 2 The proof follows along the lines of the symmetrization result in [20, Sec. 4.5]. Let {z′m }M m=1 be an independent copy

(78)

(79)

which is the required result. The proof for the second inequality also follows similarly. A PPENDIX C B OUND IN (48) Consider i.i.d. random variables {Xk }K k=1 in [0, Br ] with cumulative distribution F (x) and mean µ := EXk . For any 0 ≤ t ≤ µ, applying Markov’s inequality to the nonnegative random variable Br − Xk gives [20, (1.5)] F (µ − t) = P [Xk ≤ µ − t] = P [Br − Xk ≥ Br − µ + t] r −µ ≤ 1 − Btr . ≤ BBr −µ+t

(80)

Since Xmax := maxk Xk has distribution F (t)K , the tailintegral formula yields: Z µ E[(µ − Xmax )+ ] = P [Xmax < µ − t] dt (81) 0 K Z µ Z Br  t 1− = F (µ − t)K dt ≤ dt ≤ BKr . (82) Br 0 0

14

A PPENDIX D L IPSCHITZ C ONTINUITY OF COVERAGE PROBABILITIES IN (60) AND (65) In this section we show that the coverage probabilities in (60) and (65) are Lℓ (e−1 + 4−1 )-Lipschitz and Lℓ e−1 Lipschitz respectively, with respect to the W1 metric, hence verifying Assumption 1. The approach is similar for both the coverage probability expressions. The functional mean-value theorem gives [33] Z 1Z r(λ, z) − r(λ′ , z) = ϕt (y)(λ(y) − λ′ (y))dydt 0 A Z ′ = ϕ(y)(λ(y) − λ (y))dy (83) A

t ,z) where ϕt (y) := dr(λ (y), λt = (1 − t)λ′ + tλ, and ϕ(y) := dλ R1 1 0 ϕt (y)dt. If ϕ is Lϕ -Lipschitz, we can write Lϕ ϕ ∈ L1 and use (11) to obtain

|r(λ, z) − r(λ′ , z)| ≤ Lϕ W1 (λ, λ′ )

(84)

It remains to determine Lϕ for the two cases.

We begin by calculating its functional derivative [34, Sec. 2.2]: Z dpc (λ, z) λ(x)h(λ, x)b(x, y)dx. (85) (y) = h(λ, y) − dλ A We observe from the triangle inequality that Z 1 Z 1 ∇y ϕt (y)dt ≤ k∇y ϕt (y)k dt ≤ sup k∇y ϕt (y)k t∈[0,1]

0

≤ sup k∇y h(λt , y)k t∈[0,1]

+ sup t∈[0,1]

Z

λt (x)h(λt , x)∇y b(x, y)dx =: T1 + T2 (86)

A

For the first term, for any t ∈ [0, 1], the triangle inequality yields Z k∇y h(λt , y)k = h(λt , y) ∇y a(y) + λt (u)∇y b(y, u)du A   Z ≤ h(λt , y) k∇y a(y)k + λt (u)∇y b(y, u)du . A

Recall that since e

Lℓ r

ℓ(r) is increasing in r, we have that |ℓ′ (r)| ≤ Lℓ . r≥0 ℓ(r)

sup

For y 6= z, since k∇y dy k = 1, we can write κ κ k∇y a(y)k = 2 |ℓ′ (dy )| ≤ Lℓ = Lℓ a(y) ℓ (dy ) ℓ(dy ) τ ℓ(du ) ℓ′ (dy )(∇y dy ) ∇y b(y, u) = − (τ ℓ(du ) + ℓ(dy ))2 so using the fact that ℓ(du ) ≥ 0, we obtain Z Z λt (u)∇y b(y, u)du ≤ Lℓ λt (u)b(y, u)du A

A

∇y b(x, y) =

τ ℓ(dx )ℓ′ (dy ) (∇y dy ) (τ ℓ(dy ) + ℓ(dx ))2

(91)

so using k∇y dy k = 1 for y 6= z and the inequality 4ab ≤ (a + b)2 , we have the bound Z λt (x)h(λt , x)∇y b(x, y)dx A Z |ℓ′ (dy )| τ ℓ(dx )ℓ(dy ) ≤ dx (92) λt (x)h(λt , x) ℓ(dy ) A (τ ℓ(dy ) + ℓ(dx ))2 Z Lℓ λt (x)h(λt , x)dx = 14 Lℓ pc (λt , z) ≤ Lℓ /4 (93) ≤ 4 A for t ∈ [0, 1]. Therefore we have that T2 ≤ Lℓ /4. Adding the bounds on T1 and T2 , we can take Lϕ = Lℓ (e−1 + 4−1 ) < 2Lℓ /3 and hence pc (λ, z) is 32 Lℓ -Lipschitz in λ with respect to the W1 metric. B. Lipschitz continuity of pac (·, z) in (65)

A. Lipschitz continuity of pc (·, z) in (60)

0

+ R Combining these and substituting w := a(y) −w ≤ A λt (u)b(y, u)du, we obtain k∇y h(λt , y)k ≤ Lℓ we Lℓ Lℓ e . Hence the bound on T1 becomes T1 ≤ e . For the second term in (86), observe that

(87)

Defining p(r, z) :=

Z

λ(y)dy =

Z

λ(y)11{dy ≤r} dy

(94)

the coverage probability in (65) takes the form Z ∞ dp(r, z) pac (λ, z) = ha (λ, r)e−p(r,z) dr dr 0

(95)

B(z,r)

A

where dy := ky − zk and the expression for ha (λ, r) is provided in (70). Note that for a fixed y ∈ A,   dp(r, z) dp(r, z) d = δ(r − dy ), (y) = 1 {dy ≤r} dλ dλ(y) dr (96) where δ(·) denotes the Dirac delta function. Hence, the functional derivative of pac is given by Z ∞ 1 {dy >r} dp(r, z) dpac (λ, z) (y) = − ha (λ, r)e−p(r,z) dr ℓ(r) dλ dr 1 + τ ℓ(d 0 ) y Z ∞ dp(r, z) a −p(r,z) − h (λ, r)e 1 {dy ≤r} dr dr Z0 ∞ + ha (λ, r)e−p(r,z) δ(r − dy )dr (97) 0 Z ∞ Z dy a h (λ, r) ha (λ, r)fRz (r)dr fRz (r)dr − =− 1 + τ ℓ(r) d 0 y ℓ(dy ) + ha (λ, dy )e−p(dy ,z) =: ϕa (y)

(88) (89)

(90)

(98)

Since all the terms in ϕa (y) depend on y only through dy , let s := dy , so that it suffices to use the fact that |dy − dy′ | ≤ d a ϕt (s)| where ϕat (s) is defined ky − y′ k and only bound | ds in (98) with λ replaced by λt = (1 − t)λ′ + tλ. With a slight abuse of notation, define a(s) :=

κ ℓ(s)

b(r, s) :=

τ ℓ(s) , τ ℓ(s) + ℓ(r)

(99)

15

which are simply (62) when the serving and interfering BSs are at distances r and s, respectively. Then Z s ϕat (s) = ha (λt , s)e−pt (s,z) − ha (λt , r)b(r, s)fRz ,t (r)dr 0 Z ∞ − ha (λt , r)fRz ,t (r)dr (100) s

where pt is defined as in (94) and fRz ,t in (67) with λ replaced with λt . Hence, we obtain ha (λt , s)e−pt (s,z) = exp −a(s) − pt (s, z) −

Z

λt (u)b(s, du )du A\B(z,s)

!

.

(101)

Using the Leibniz rule and the fact that b(s, s) = τ /(1 + τ ), we obtain i d h a h (λt , s)e−pt (s,z) = −ha (λt , s)e−pt (s,z) × (102) ds # " Z 1 dpt (s, z) ∂b(s, du ) ′ . du + a (s) + λt (u) ∂s 1+τ ds A\B(z,s)

A PPENDIX E I NCREMENTAL D EPLOYMENT As earlier, let Φ0 denote the set of existing BSs. Under the instantaneous-SINR association rule, a user z may be covered either through one of the installed BSs u ∈ Φ0 with probability X Y (1 − b(u, v)) (107) p1c = h(λ, u) u∈Φ0

v∈Φ0 \{u}

or through one of the newly deployed BSs with probability Z Y 2 (1 − b(x, u))dx. (108) λ(x)h(λ, x) pc = A

u∈Φ0

Since τ > 1, at most one BS can satisfy the SINR threshold. Therefore, the total coverage probability is given by p1c + p2c . Coverage probability of a specified set Φ0 of BSs can likewise be obtained by substituting λ ≡ 0 in (107). Under average-SINR association, the serving BS is the nearest BS to the user. Let uz denote the existing BS that is closest to z and ru := kuz − zk. Recall that xz is the newly deployed BS closest to z and Rz = kxz − zk. There are two distinct possibilities: Rz ≥ ru so that the user is served by an existing BS, and Rz < ru so that the user is served by a new BS. The contribution to the coverage probability from the first case is Y (1 − b(uz , v)) (109) p1c =

Likewise, differentiating the two integrals in (100) gives Z s  d ha (λt , r)b(r, s)fRz ,t (r)dr ds 0 dpt (s, z) τ ha (λt , s)e−pt (s,z) = τ ds Z1 + s ∂b(r, s) a v∈Φ0 \{uz } + h (λt , r) fRz ,t (r)dr (103) ! ∂s Z Z 0 Z ∞  d λ(y)b(uz , y)dy λ(y)dy − × exp −a(uz ) − ha (λt , r)fRz ,t (r)dr A\B(z,ru ) B(z,ru ) ds s dpt (s, z) where a(·) is defined in (62). When a newly deployed BS . (104) = −ha (λt , s)e−pt (s,z) serves z, the coverage probability is given by ds   Substituting (102)–(104) into the derivative of (100), the terms Z Z ru 2 (s,z) τ σ λ(y)dy cancel, and we obtain containing dptds  exp − − p2c = ℓ(r) " # P ℓ(r) Z 0 A\B(z,r) 1 + τ ℓ(ky−zk) ∂b(s, du ) (ϕat )′ (s) = − ha (λt , s)e−pt (s,z) a′ (s) + λt (u) du Y ℓ(r) ∂s A\B(z,s) × (110) fR (r)dr. Z s ℓ(r) + τ ℓ(kv − zk) z ∂b(r, s) v∈Φ0 a fRz ,t (r)dr. (105) − h (λt , r) ∂s 0 Finally, the total coverage probability is given by p1c + p2c . whose terms can be bounded using the same arguments as those used for T1 and T2 in the instantaneous-association A PPENDIX F proof. Specifically, defining C ONVEXITY IN THE N OISE -L IMITED C ASE Z In this section, we show that pc (λ) for IPPP deployment w := a(s) + λt (u)b(s, du )du, (106) A\B(z,s) is concave when the network is noise-limited, i.e., when the we have ha (λt , s)e−pt (s,z) ≤ e−w , and therefore the mag- aggregate interference power I(z) is much smaller than the 2 nitude of the first term in (105) is bounded as in T1 by noise power σ . L Lℓ we−w ≤ eℓ . For the second term, we can likewise bound Lemma 7. In the noise-limited case, the coverage probability − ∂b(r,s) ≤ L4ℓ so that expressions in (60) and (65) are both concave in λ. ∂s Z s Lℓ Lℓ a ∂b(r, s) Proof: In the instantaneous-SINR association, the noisefRz ,t (r)dr ≤ pc (λt , z) ≤ . − ha (λt , r) ∂s 4 4 limited coverage probability is given by [27, Sec. II]: 0  Z  Different from the earlier proof, the first term in (105) −a(x) pc (λ, z) = 1 − exp − λ(x)e dx (111) is nonpositive, while the second is nonnegative, and hence A − Leℓ ≤ ϕ′t (s) ≤ L4ℓ which implies that we can take Lϕ = Leℓ . In conclusion, pac (λ, z) is Lℓ /e-Lipschitz with respect to the which is concave since the function 1 − e−x is concave in x W1 metric. and the integral is linear in λ.

16

In the average-SINR association, the noise-limited version of (65) becomes   Z ∞ τ σ2 (112) fRz (r)dr pac (λ, z) = exp − P ℓ(r) 0    Z ∞ τ σ2 = P [Rz ≤ r] d − exp − (113) P ℓ(r) 0   R where P [Rz ≤ r] = 1 − exp − B(z,r) λ(y)dy is concave h  i σ2 in λ and the measure d − exp − Pτℓ(r) is nonnegative. Hence, pac (λ, z) is a nonnegative weighted integral of concave functions and is therefore concave in λ. Finally, expectation over z and empirical averaging preserve concavity, and hence the coverage probability is concave in both the cases. In the noise-limited case therefore, mirror ascent converges to the global optimum and Eo can be made arbitrarily small. To obtain the bounds in closed-form however, we must first establish that the discretized gradients of pc (λd , z) and pac (λd , z) have bounded components. Lemma 8. In the noise-limited case, it holds that A k∇pc (λd , z)k∞ ≤ D

A k∇pac (λd , z)k∞ ≤ D

(114)

where the expressions for pc and pac are provided in (60) and (65), respectively. Proof: Observe that in (85), both h(λ, y) and b(x, y) lie in the interval [0, 1]. Therefore, the right-hand side is bounded (λ,z) (y)| ≤ 1. The required bound follows as as | dpcdλ k∇pc (λd , z)k∞ = max|[∇pc (λd , z)]i | (115) i Z dpc (λ, z) A (y)dy| ≤ |Ai | = D = max| . i dλ Ai For the average-SINR association case, we can  again use

these bounds and the fact that 0 ≤

1 + τ ℓ(r) ℓ(dy )

dpa (λ,z) | cdλ (y)|

−1

≤ 1 in

≤ 1. Therefore, (98), so as to conclude that proceeding as above, we obtain the required bound. Since the stochastic gradient g(t) in (43) is unbiased and, from the preceding lemma, satisfies g(t) ∞ ≤ A/D, we can directly apply the standard stochastic mirror-ascent bound for concave objectives [25]. In particular, using λ̄Td as the stepsize-weighted average of the iterates and starting from uniform initialization and constant step-size, we have r 2 log D . (116) Eo ≤ N T We can then add this bound to the other terms in Thm. 1 to obtain the full end-to-end bound. A PPENDIX G E FFECT OF NOISE POWER In this section, we again consider the experiment with users distributed along an annulus, as in Fig. 1(a), but with different values of noise power σ 2 . As expected, Fig. 7 shows that when the noise is high, the optimized intensity remains close to the users to reduce path loss. As noise decreases, a more

centrally concentrated intensity can provide coverage to users around the annulus. However, as noise is reduced further, interference becomes relatively more important, and deploying multiple BSs at the center results in interference. Spreading the intensity back towards the annulus increases the likelihood that a user is close to one BS while the others are farther away, thereby improving the SINR and coverage probability. A PPENDIX H E FFECT OF THE A SSOCIATION RULE In this section, we illustrate how the assumed association rule can affect the optimized intensity. We consider the region A = [−5, 5]2 with two existing BSs placed symmetrically on both sides of a corridor containing 80% of the 300 users. We also have a separate group of users at (3, 3) containing the remaining users. We take σ 2 = 2 × 10−5 , while the remaining propagation parameters are the same as specified in Table I. The noise parameter is chosen to make the effect of the association rule more visible. Fig. 8 compares the intensities optimized under instantaneous-SINR association (IA) and average-SINR association (AA). The difference is concentrated in the corridor, where the users lie between the two existing BSs. Under AA, each user is associated with the nearest BS, while the other BS contributes only interference; hence, some additional intensity is allocated within the corridor. Under IA, either existing BS may serve a user depending on the instantaneous fading, reducing the need for additional deployment in this region. The user group centered at (3, 3) is far from both existing BSs, and therefore the two association rules yield nearly identical intensities there. Finally, the optimized training coverage probability is pc = 0.71 under instantaneous-SINR association and pc = 0.61 under averageSINR association. Thus, the association rule can significantly affect both coverage and deployment priorities in regions containing multiple existing BSs. R EFERENCES [1] S. Lee, S. Lee, K. Kim, and Y. H. Kim, “Base station placement algorithm for large-scale LTE heterogeneous networks,” PloS one, vol. 10, no. 10, p. e0139190, 2015. [2] R. Jurdi, J. G. Andrews, D. Parsons, and R. W. Heath, “Identifying coverage holes: Where to densify?” in Proc. Asilomar Conference on Signals, Systems, and Computers, 2017, pp. 1417–1421. [3] L. Suárez, L. Nuaymi, and J.-M. Bonnin, “Energy-efficient bs switchingoff and cell topology management for macro/femto environments,” Computer Networks, vol. 78, pp. 182–201, 2015. [4] Q.-N. Le-The, T. Beitelmal, F. Lagum, S. S. Szyszkowicz, and H. Yanikomeroglu, “Cell switch-off algorithms for spatially irregular base station deployments,” IEEE Wireless Commun. Lett., vol. 6, no. 3, pp. 354–357, 2017. [5] H. Wang, H. Zhao, W. Wu, J. Xiong, D. Ma, and J. Wei, “Deployment algorithms of flying base stations: 5g and beyond with uavs,” IEEE Internet of Things Journal, vol. 6, no. 6, pp. 10 009–10 027, 2019. [6] J. Lyu, Y. Zeng, R. Zhang, and T. J. Lim, “Placement optimization of uav-mounted mobile base stations,” IEEE Commun. Lett., vol. 21, no. 3, pp. 604–607, 2016. [7] H. Zhao, H. Wang, W. Wu, and J. Wei, “Deployment algorithms for UAV airborne networks toward on-demand coverage,” IEEE J. Sel. Areas Commun., vol. 36, no. 9, pp. 2015–2031, 2018. [8] H. Huang and A. V. Savkin, “Deployment of heterogeneous UAV base stations for optimal quality of coverage,” IEEE Internet of Things Journal, vol. 9, no. 17, pp. 16 429–37, 2022.

17

3

3

3

0.7

0

0

0.6

0

0.5

-3 -3

0

3

3

-3 -3

0

3

-3 -3

0

3

0.4

3

3

0.3

0.2 0

0

0

0.1

-3 -3

0

3

-3 -3

0

3

-3 -3

0 0

3

Fig. 7. Intensity maps across a range of noise powers.

[9] T. Akram, M. Awais, R. Naqvi, A. Ahmed, and M. Naeem, “Multicriteria UAV base stations placement for disaster management,” IEEE Syst. Journal, vol. 14, no. 3, pp. 3475–3482, 2020. [10] S. Fu, X. Feng, A. Sultana, and L. Zhao, “Joint power allocation and 3D deployment for uav-bss: A game theory based deep reinforcement learning approach,” IEEE Trans. Wireless Commun., vol. 23, no. 1, pp. 736–748, 2023. [11] J. Andrews, A. Gupta, A. Alammouri, and H. Dhillon, An introduction to cellular network analysis Using stochastic geometry. Springer, 2023. [12] H. S. Dhillon, R. K. Ganti, F. Baccelli, and J. G. Andrews, “Modeling and analysis of K-tier downlink heterogeneous cellular networks,” IEEE J. Sel. Areas Commun., vol. 30, no. 3, pp. 550–60, 2012. [13] A. K. Gupta, J. G. Andrews, and R. W. Heath, “On the feasibility of sharing spectrum licenses in mmwave cellular systems,” IEEE Trans. Commun., vol. 64, no. 9, pp. 3981–95, 2016. [14] M. Di Renzo, A. Zappone, T. T. Lam, and M. Debbah, “Systemlevel modeling and optimization of the energy efficiency in cellular networks—a stochastic geometry framework,” IEEE Trans. Wireless Commun., vol. 17, no. 4, pp. 2539–56, 2018. [15] C. Saha, H. S. Dhillon, N. Miyoshi, and J. G. Andrews, “Unified analysis of hetnets using Poisson cluster processes under max-power association,” IEEE Trans. Wireless Commun., vol. 18, no. 8, pp. 3797–3812, 2019. [16] S. E. Hammami, H. Afifi, M. Marot, and V. Gauthier, “Network planning tool based on network classification and load prediction,” in Proc. IEEE WCNC, 2016, pp. 1–6. [17] Ł. Kułacz and A. Kliks, “Dynamic spectrum allocation using multisource context information in openRAN networks,” Sensors, vol. 22, no. 9, p. 3515, 2022. [18] B. Błaszczyszyn, M. K. Karray, and H. P. Keeler, “Wireless networks appear Poissonian due to strong shadowing,” IEEE Trans. Wireless Commun., vol. 14, no. 8, pp. 4379–90, 2015. [19] A. T. McRae, C. J. Cotter, and C. J. Budd, “Optimal-transport–based mesh adaptivity on the plane and sphere using finite elements,” SIAM Journal on Scientific Computing, vol. 40, no. 2, pp. A1121–A1148, 2018.

[20] F. Bach, Learning theory from first principles. MIT press, 2024. [21] A. Maurer, “A vector-contraction inequality for Rademacher complexities,” in Intl. Conf. on Algorithmic Learning Theory, 2016, pp. 3–17. [22] P. Massart, Concentration inequalities and model selection: Ecole d’Eté de Probabilités de Saint-Flour XXXIII-2003. Springer, 2007. [23] L. Condat, “Fast projection onto the simplex and the ℓ1 ball,” Mathematical Programming, vol. 158, no. 1, pp. 575–585, 2016. [24] I. Fatkhullin and N. He, “Taming nonconvex stochastic mirror descent with general bregman divergence,” in International Conference on Artificial Intelligence and Statistics. PMLR, 2024, pp. 3493–3501. [25] A. Nemirovski, A. Juditsky, G. Lan, and A. Shapiro, “Robust stochastic approximation approach to stochastic programming,” SIAM Journal on optimization, vol. 19, no. 4, pp. 1574–1609, 2009. [26] M. Haenggi, Stochastic geometry for wireless networks. Cambridge University Press, 2012. [27] J. G. Andrews, A. K. Gupta, and H. S. Dhillon, “A primer on cellular network analysis using stochastic geometry,” arXiv preprint arXiv:1604.03183, 2016. [28] M. Haenggi, J. G. Andrews, F. Baccelli, O. Dousse, and M. Franceschetti, “Stochastic geometry and random graphs for the analysis and design of wireless networks,” IEEE J. Sel. Areas Commun., vol. 27, no. 7, pp. 1029–46, 2009. [29] F. Baccelli and B. Błlaszczyszyn, “Stochastic geometry and wireless networks: Volume II applications,” Found. Trends Netw., vol. 4, no. 1–2, p. 1–312, Jan. 2009. [30] Facebook, “Data for Good,” 2025, [Online; accessed 31-March-2025]. [Online]. Available: https://dataforgood.facebook.com/dfg [31] OpenCelliD, “OpenCelliD: Open Database of Cell Towers,” https://opencellid.org/, accessed: Aug. 2026. [32] J. Muscat, Functional analysis: an introduction to metric spaces, Hilbert spaces, and Banach algebras. Springer Nature, 2024. [33] P. Sahoo and T. Riedel, Mean value theorems and functional equations. World Scientific, 1998.

18

Fig. 8. Optimized intensities for the instantaneous-SINR (left) and average-SINR (right) association rules.

[34] K. W. Cassel, Variational methods with applications in science and engineering. Cambridge University Press, 2013.

Related documents

• Stochastic Geometry Meets Optimization: Data-Driven Cellular Planning#1446838
• Optimal Location Configuration of Base Sytations and Frequency Assignment for Cellular Mobile Networks.#384996
• From Data to Deployment: A Scoping Review of Empirical Tools, Techniques, and Operational Readiness in Data-Driven 5G/6G Network Optimization#1401728
• Enabling Mobile Base Stations in 5G via Wireless Access Backhaul (WAB): A Multi-Band Experimental Study#798676
• Enabling Mobile Base Stations in 5G via Wireless Access Backhaul (WAB): A Multi-Band Experimental Study#797262