## 3 Proposed Optimization Algorithm
### BS Transmit Power Allocation
#### 3.1.2 Two-Layer Approach for Solving Problem (P1)  
According to Theorem 1, the problem (P1) can be solved via a two-layer iterative approach. We tackle the parametric problem \((\mathrm{P1}^{\prime})\) for a given pair \((\mathbf{\mu},\mathbf{\rho})\) in the inner layer. In the outer layer, we optimize \((\mathbf{\mu},\mathbf{\rho})\) to ensure that \((\mathbf{\mu},\mathbf{\rho})\), together with the corresponding inner-layer solution \(\mathbf{P}\), satisfy the conditions (28) and (29). It is straightforward to verify that, when \((\mathbf{\mu},\mathbf{\rho})\succeq 0\), the inner-layer problem \((\mathrm{P1}^{\prime})\) is a linear program (LP) due to the linear nature of the objective function and constraints. This analytical structure enables us to derive the optimal value in a straightforward manner. To proceed, we initially define \(v_{1}\), \(v_{2}\) and \(v_{3}\) as follows:  
\[v_{1}=\sum_{n\in\mathcal{N}}\mu_{n}a_{n}\sum_{m\in\mathcal{M}} \gamma_{m,n}d_{m,n}, \tag{30}\] \[v_{2}=\sum_{n\in\mathcal{N}}\mu_{n}\rho_{n}\left(\sum_{n\in \mathcal{N}}\frac{\sum\limits_{m\in\mathcal{M}}\gamma_{m,n}u_{m,n}d_{m,n}}{ \sum\limits_{m\in\mathcal{M}}\gamma_{m,n}d_{m,n}}-\frac{\sum\limits_{m\in \mathcal{M}}\gamma_{m,n}u_{m,n}d_{m,n}}{\sum\limits_{m\in\mathcal{M}}\gamma_{ m,n}d_{m,n}}\right),\] (31) \[v_{3}=\sum_{n\in\mathcal{N}}\mu_{n}\rho_{n}\sigma^{2}. \tag{32}\]  
Then, (26) can be expressed as  
\[F=(v_{1}-v_{2})P_{m,n}-v_{3}. \tag{33}\]  
From (33), we can see that \(P_{m,n}=P_{m}^{c}\) if \(F=F^{*}\) (\(F^{*}\) denotes the optimal value of \((\mathrm{P1}^{\prime})\)), otherwise \(P_{m,n}=0\). After solving the problem (P1\({}^{\prime}\)) with given \((\mathbf{\mu},\mathbf{\rho})\), we employ the modified Newton method [40] to update \((\mathbf{\mu},\mathbf{\rho})\) in the outer layer, which guarantees the convergence to \((\mathbf{\mu}^{*},\mathbf{\rho}^{*})\) (i.e., meets the conditions (28) and (29)). To begin with, we define some functions for conciseness as follows:  
\[\psi_{n}^{1}(\mu_{n})=\mu_{n}D_{n}-1, \tag{34}\] \[\psi_{n}^{2}(\rho_{n})=\rho_{n}D_{n}-U_{n}. \tag{35}\]  
According to ([40], Theorem 3.1), the optimal \((\mathbf{\mu}^{*},\mathbf{\rho}^{*})\) is unique and can be obtained by solving the non-linear equations: \([\psi_{1}^{1},\cdots\psi_{N}^{1},\psi_{1}^{2},\cdots,\psi_{N}^{2}]=\mathbf{0}\). The modified Newton method can be applied to update \((\mathbf{\mu},\mathbf{\rho})\) to meet these equations. Specifically, the point-wise updating equations of \(\mathbf{\mu}\) and \(\mathbf{\rho}\) are  
\[\mu_{n}^{i+1}=(1-\delta^{i})\mu_{n}^{i}+\delta^{i}\frac{1}{D_{n} ^{i}}, \tag{36}\] \[\rho_{n}^{i+1}=(1-\delta^{i})\rho_{n}^{i}+\delta^{i}\frac{U_{n} ^{i}}{D_{n}^{i}}, \tag{37}\]  
where \(i\) is the iteration index and \(\delta^{i}\in(0,1)\) is the step size at iteration \(i\). A discussion on choosing the step size can be found in [40; 41] and is thus ```
1:Input the BS selections \(\mathbf{\Gamma}\) and the update accuracy control factor \(\epsilon>0\). 2:Initialize \((\mathbf{\mu},\mathbf{\rho})\succeq 0\). 3:loop
4: Given \((\mathbf{\mu},\mathbf{\rho})\), obtain the optimal \(\mathbf{\mathsf{P}}\) to problem (P1\({}^{\prime}\)). 5: Update \((\mathbf{\mu},\mathbf{\rho})\) using (36) and (37). 6:if\(\sum_{n\in\mathcal{N}}(|\psi_{n}^{1}|^{2}+|\psi_{n}^{2}|^{2})<\epsilon\)then
7:break
8:endif
9:endloop
10:Output:\(\mathbf{\Gamma}\). ```

**Algorithm 1** Two-Layer Iterative Algorithm for Problem (P1)

### BS Selection Design

In this subsection, we develop an approach based on the reformulation-linearization technique to tackle the BS selection problem (P2). The objective function of problem (P2) also follows the sum-of-ratios format. Therefore, similar to the algorithm design for (P1), problem (P2) can be effectively solved using a two-layer iterative approach.