# Approaching Globally Optimal Energy Efficiency in Interference Networks via
Machine Learning

Bile Peng, , Karl-Ludwig Besser, , Ramprasad Raghunath, , and Eduard A.
Jorswieck B. Peng, K.-L. Besser, R. Raghunath and E. A. Jorswieck are with the
Department of Information Theory and Communication Systems, TU Braunschweig,
Schleinitzstr. 22, 38112 Braunschweig, Germany (e-mail: {peng, raghunath,
besser, jorswieck}@ifn.ing.tu-bs.de).The work is supported by the Federal
Ministry of Education and Research Germany (BMBF) as part of the 6G Research
and Innovation Cluster 6G-RIC under Grant 16KISK020K.

###### Abstract

This work presents a machine learning approach to optimize the energy
efficiency (EE) in a multi-cell wireless network. This optimization problem is
non-convex and its global optimum is difficult to find. In the literature,
either simple but suboptimal approaches or optimal methods with high
complexity and poor scalability are proposed. In contrast, we propose a
machine learning framework to approach the global optimum. While the neural
network (NN) training takes moderate time, application with the trained model
requires very low computational complexity. In particular, we introduce a
novel objective function based on stochastic actions to solve the non-convex
optimization problem. Besides, we design a dedicated NN architecture for the
multi-cell network optimization problems that is permutation-equivariant. It
classifies channels according to their roles in the EE computation. In this
way, we encode our domain knowledge into the NN design and shed light into the
black box of machine learning. Training and testing results show that the
proposed method without supervision and with reasonable computational effort
achieves an EE close to the global optimum found by the branch-and-bound
algorithm. Hence, the proposed approach balances between computational
complexity and performance.

###### Index Terms:

Energy efficiency, non-convex optimization, machine learning, permutation-
equivariance, reparameterization trick.

##  I Introduction

Transmit power control belongs to the most fundamental problems in research on
cellular mobile networks [1] and it is an important measure to optimize
objectives such as weighted sum-rate [2], energy efficiency (EE) [3] and
fairness [4]. These problems share some significant similarities, among which
the non-convexity due to the fractional nature of the signal-to-interference-
noise ratio (SINR) expression is a major difficulty, because multiple local
optima could exist and convex optimization techniques are not applicable.

This paper focuses on the EE optimization since the energy consumption is a
key performance indicator for the future sustainable wireless communication
systems. Among different techniques, transmit power control is an important
approach to improve the EE in the air interface. On the other hand, given the
similarities of the above-mentioned power control problems in cellular
networks and the fact that the EE maximization problem with multiple users is
particularly difficult [5], the proposed solution is expected to have good
generalizability to other power control problems, if it can solve the EE
maximization problem effectively.

In the literature, fractional programming is introduced to optimize the EE in
[6], the dual problem of the EE maximization in a multi-cell orthogonal
frequency-division multiple access (OFDMA) system is solved in [7], the EE in
a multiple-input-multiple-output (MIMO) broadcast channel is optimized with an
iterative water filling scheme in [8], the sum of EEs of multiple users is
optimized with fractional-programming and alternating optimization method in
[9], the weighted sum EE is optimized with the two-layer optimization in [10],
the max-min fairness of the EE is optimized with the dual Dinkelbach-type
algorithm in [11], the EE in OFDMA is optimized in three steps in [12], the
global EE and the minimum EE are maximized with fractional programming with
sequential optimization theory in [13], the global optimum in the non-convex
problem is found with the successive incumbent transcending [14] and the
branch-and-bound algorithm [3]. These works either do not promise the global
optimum or have a large computation complexity. In order to solve this
problem, [3, 15] propose to train an artificial neural network (ANN) with
targets obtained by the branch-and-bound algorithm, i.e., supervised learning,
as a quicker online approach compared to the original solution. However, we
still need to perform the offline branch-and-bound algorithm in order to
obtain the target and the overall computational effort is still high.

Considering the limitations of the above-mentioned works, we propose an
unsupervised, model-based machine learning approach to solve the EE
maximization problem. The unsupervised learning does not require “correct”
answers as targets but lets the algorithm find the optimal transmit power by
itself. Therefore, the requirement on training data and computational power is
considerably lower than supervised learning. In model-based machine learning,
the optimizer knows the objective function and computes its gradient to
optimize the neural network (NN) (unlike the model-free machine learning,
which optimizes the NN based on feedback from unknown mechanics in a trial-
and-error manner [16, 17]). Since the model-based learning has more
information about the objective than the model-free learning, it usually has a
higher sample efficiency [18]. In the past years, machine learning has been
successfully applied to cellular network optimizations [19, 20] and in
particular, to EE optimization problems in cellular networks [21, 22].
Compared to these existing machine learning methods, our proposed approach has
the following two innovations:

  * •

We propose an innovative objective function, which replaces the original EE
expression in the gradient ascent method to overcome the non-convexity problem
(Section III). The non-convexity is a long-existing problem in machine
learning with possible counter-measures such as alternating optimization,
expectation-maximization (EM) algorithm, dual problem and stochastic
optimization [23]. Our proposed solution is based on the stochastic
optimization and the reparameterization trick, which approximates the
expectation of the objective with the mean of random variables. The
reparameterization trick makes the random variable differentiable and
therefore enables stochastic gradient ascent (SGA). It is shown that this
approach has a smaller variance compared to other methods [24]. In recent
years, it has been successfully applied to Bayesian optimization [25],
variational auto-encoder [26] and reinforcement learning [27], but not yet in
non-convex unsupervised learning, in particular, not for the cellular network
optimization problem.

  * •

We design a permutation-equivariant NN architecture, i.e., if the order of the
input pairs of user and base station (BS) are permuted, the order of the
output of the NN is permuted automatically in the same way [28, 29] (Section
IV). This increases the generality of the NN significantly, since the
permutation-equivariance is an inherent property of the power control problems
in cellular networks. However, classical NN architectures, such as NNs with
fully connected layers, do not have this property. In the literature, a
dedicated permutation-equivariant NN architecture has been applied to signal
detection [30] and the permutation-equivariant graphical neural network (GNN)
has been applied to capacity maximization [31, 32, 33]. However, these works
are not particularly optimized for the non-convexity of the considered problem
and the performance is not compared with the global optimum.

Note that this proposed approach is an _unsupervised_ learning method.
Therefore, we do not need to provide any target values for a given input. In
contrast, supervised learning methods require a target for each input. The
optimizer tries to minimize the difference between the output of the neural
network and the target (or label). The loss function is usually defined as
mean squared error (MSE) (for continuous output) or cross entropy (for
categorical output). By minimizing the loss function, the neural network
learns to reproduce the given target. However, in our problem, the target is
unknown and can therefore not be provided to the optimizer. By maximizing the
objective function with respect to massive data samples, the neural network
learns an optimized mapping from the channel matrix $\mathbf{G}$ to the
transmit powers $\mathbf{p}$ that maximizes the EE. If the training data are
massive and general enough, we can expect that the learnt mapping achieves a
good performance as well for data not in the training data set. We believe
that the unsupervised approach has a wider range of applications since the
target is unavailable for many difficult problems. Furthermore, it is also
very exciting to see that the neural network finds a good solution by itself
without any supervision.

_Notation:_ A normal small letter (e.g., $s$) denotes a scalar, a bold small
letter (e.g., $\mathbf{v}$) denotes a vector, a bold capitalized letter (e.g.,
$\mathbf{M}$) denotes a matrix or a tensor of dimension higher than 2,
$\mathbf{M}^{T}$ denotes the transpose of matrix $\mathbf{M}$, $|\mathcal{S}|$
denotes the cardinality of set $\mathcal{S}$, $\mathrm{dim}(\mathbf{v})$
denotes the number of dimensions of vector $\mathbf{v}$,
$\text{diag}(\mathbf{M})$ denotes the diagonal elements of matrix
$\mathbf{M}$, $\mathbf{I}$ denotes the identity matrix, $\mathbf{0}$ denotes
the matrix of all zeros and $\mathbf{1}$ denotes the matrix of all ones, the
extraction matrix $\mathbf{E}$ is defined as
$\mathbf{E}=\mathbf{1}-\mathbf{I}$, $\mathcal{U}(\mathbf{a},\mathbf{b})$
denotes a uniform distribution between $\mathbf{a}$ and $\mathbf{b}$ where the
$i$-th component is uniformly distributed between $a_{i}$ and $b_{i}$ and all
components in the distribution are independent from each other, $\odot$
denotes the element-wise product of two vectors,
$\prod\limits^{\circ}\mathbf{v}$ denotes the product of elements of
$\mathbf{v}$, $\sum\limits^{\circ}\mathbf{v}$ denotes the sum of elements of
$\mathbf{v}$, $\mathbf{T}[i,j,k]$ denotes element in position $(i,j,k)$ in
tensor $\mathbf{T}$, $\mathbf{T}[\cdot,j,k]$ denotes the vector of the
elements with position $(j,k)$ in the second and third dimensions in tensor
$\mathbf{T}$, $\mathbf{T}[i,\cdot,\cdot]$ denotes the matrix of the elements
with position $i$ in the first dimension in tensor $\mathbf{T}$,
$\mathbf{A}\bullet\mathbf{B}$ denotes multiplication between matrix
$\mathbf{A}$ and tensor $\mathbf{B}$, where $\mathbf{A}$ is multiplied with
all vectors $\mathbf{B}[\cdot,i,j]$ for all $i$ and $j$,
$\mathbf{A}\circ\mathbf{B}$ denotes multiplication between matrix $\mathbf{A}$
and tensor $\mathbf{B}$, where $\mathbf{A}$ is multiplied with all matrices
$\mathbf{B}[i,\cdot,\cdot]$ for all $i$, or multiplication between tensor
$\mathbf{A}$ and matrix $\mathbf{B}$, where all matrices
$\mathbf{A}[i,\cdot,\cdot]$ (for all $i$) is multiplied with $\mathbf{B}$.

##  II Problem Statement

###  II-A Problem Formulation

We consider an uplink transmit power control problem in a multi-cell wireless
network, as shown in Figure 1, where $I$ single-antenna users are served by
$M$ BSs, which are equipped with $n_{R}$ antennas. We assume that every user
is assigned to one BS and denote the assigned BS of user $i$ as BS $i$. We
further assume that the BSs use the matched filter to process the received
signal. User $i$ chooses a transmit power $p_{i}\in[0,p_{\max}]$, where
$p_{\max}$ is the maximum transmit power111Information available at the user
is usually limited in practice, which might be insufficient to obtain the
global optimum. The BS can choose a $p_{i}$ based on the complete information
available at the BS and inform the user.. Our objective is to maximize the
sum-EE, which is defined as the sum of ratios between data rate and power
consumption of all user-BS pairs. The problem is formulated as

$\displaystyle\max_{\mathbf{p}}$ $\displaystyle
J=\sum_{i=1}^{I}\frac{\log\left(1+\frac{g_{ii}p_{i}}{1+\sum_{j\neq
i}g_{ji}p_{j}}\right)}{\mu p_{i}+P_{c}}$ (1) s.t. $\displaystyle 0\leq
p_{i}\leq p_{\max}\quad\text{ for }i=1,2,\dots,I,$

where $\mathbf{p}=(p_{1},p_{2},\dots,p_{I})$ is the vector of transmit powers
of all user-BS pairs, $g_{ji}$ is the equivalent channel gain from user $j$ to
BS $i$ assuming the matched filter receiver. Therefore, $g_{ii}$ is the
channel gain of the useful signal and $g_{ji}$, $j\neq i$, is the channel gain
of the interference from user $i$ to BS $j$, $\mu$ is the inefficiency of the
power amplifier and $P_{c}$ is the static power assumption.

BS $i$User $i$ Figure 1: Illustration of the considered scenario. The channel
of the signal of interest for the user-BS pair $i$ is drawn as the solid line
and the interference channels are drawn as the dotted lines.

###  II-B Properties of the Problem

####  II-B1 Non-convexity

A non-convex function is a function whose second-order derivative is not non-
negative on its entire domain. If a function is non-convex, it can have
multiple local optima. From (1), we notice that the objective $J$ is a non-
convex function in power $\mathbf{p}$. The conventional convex optimization
methods might converge to a poor local optimum, if we use them to optimize a
non-convex objective. In the literature, [3] proposes to solve the problem
with a branch-and-bound algorithm, which guarantees to converge to the global
optimum at the cost of high computational effort and poor scalability.
Alternatively, in standard NN optimization, we use gradient ascent/descent to
optimize the NN, which cannot not be prevented from converging to a local
optimum.

####  II-B2 Permutation-Equivariance

Denote the equivalent channel matrix as $\mathbf{G}$, where the element in row
$i$ and column $j$ is $g_{ij}$ as defined in (1). We expect that the power
allocation $\mathbf{p}$ is permutation-equivariant to $\mathbf{G}$, i.e., if
the user-BS pairs are permuted, the elements in $\mathbf{p}$ should be
permuted in the same way. This is because the order of the user-BS pairs is
inherently arbitrary and should have no impact on the optimal transmit power.
Mathematically, we use the permutation matrix $\mathbf{P}$ to describe a
permutation, whose elements are either 0 or 1 and the sum of every row and
every column is 1. For example, the matrix

$\mathbf{P}=\begin{pmatrix}1&0&0\\\ 0&0&1\\\ 0&1&0\end{pmatrix}$

is a valid permutation matrix. The product $\mathbf{P}\mathbf{G}$ permutes the
rows of matrix $\mathbf{G}$ whereas $\mathbf{G}\mathbf{P}^{T}$ permutes the
columns of $\mathbf{G}$. With the example $\mathbf{P}$ defined above, it
permutes the second and third rows (columns) of $\mathbf{G}$ while the first
row (column) is unchanged. If we apply both row and column permutation, the
resulting matrix

$\mathbf{G}^{\prime}=\mathbf{P}\mathbf{G}\mathbf{P}^{T}$ (2)

is permuted in both rows and columns in the same way simultaneously compared
to $\mathbf{G}$. Consider a mapping $\Phi:\mathbb{R}^{I\times
I}\to\mathbb{R}_{+}^{I}$ with $\Phi(\mathbf{G})=\mathbf{p}$, $\Phi$ is called
permutation-equivariant if

$\mathbf{p}\mathbf{P}^{T}=\Phi(\mathbf{P}\mathbf{G}\mathbf{P}^{T})$ (3)

holds for any permutation matrix $\mathbf{P}$. Obviously, a conventional NN
with fully connected layers is not permutation-equivariant.

##  III Application of Reparametrization Trick to Non-Convex Optimization
Objectives

In this section, we propose an innovative objective function with stochastic
actions and reparameterization trick that can allow a NN to converge to the
global optimum even in non-convex problems.

###  III-A Defining the Objective

First, we propose an innovative objective function for non-convex optimization
problems. Given a non-convex objective function $J(\mathbf{p})$, we can
maximize $J$ with gradient ascent, where $\mathbf{p}$ is updated as
$\mathbf{p}\leftarrow\mathbf{p}+r\nabla_{\mathbf{p}}J(\mathbf{p})$, with
learning rate $r$. Since $J$ is non-convex, there may exist multiple local
optima.

In Figure 2, it is illustrated that if $p$ is initialized as $p^{(0)}$ near a
poor local optimum $O_{L}$, the optimization with gradient ascent would
converge to $O_{L}$ rather than the global optimum $O_{G}$. Note that we
reduced the vector $\mathbf{p}$ to a one-dimensional scalar $p$ for
illustration simplicity.

$J$$p$$p^{(0)}$$O_{L}$$O_{G}$ Figure 2: Illustration of a non-convex objective
function $J$. When $p$ is initialized as $p^{(0)}$ near the local optimum
$O_{L}$, the classical gradient ascent optimization converges to the local
optimum $O_{L}$, which is not the global optimum $O_{G}$.

In order to solve this problem, we do not use a deterministic $\mathbf{p}$.
Instead, we assume that $\mathbf{p}$ is a random variable with
$\mathbf{p}\sim\mathcal{U}(\mathbf{a},\mathbf{b})$, where $\mathbf{a}$ and
$\mathbf{b}$ are lower and upper bound of the support of the distribution,
respectively. The expectation of the objective $J$ is

$K=\mathbb{E}_{\mathbf{p}\sim\mathcal{U}(\mathbf{a},\mathbf{b})}(J\left(\mathbf{p})\right).$
(4)

We can use the mean of multiple independent identically distributed (i.i.d.)
random variables to approximate the expectation in (4). However, we cannot
compute the gradient of a random variable and gradient ascent cannot be
applied to optimize $K$ as a result. Therefore, we use the reparameterization
trick. Instead of sampling $\mathbf{p}$ from the distribution
$\mathcal{U}(\mathbf{a},\mathbf{b})$ directly, we sample $\mathbf{p}$ from the
standard uniform distribution $\mathcal{U}(\mathbf{0},\mathbf{1})$ and rewrite
(4) as

$K=\mathbb{E}_{\mathbf{p}\sim\mathcal{U}(\mathbf{0},\mathbf{1})}(J\left(\mathbf{a}+(\mathbf{b}-\mathbf{a})\odot\mathbf{p})\right),$
(5)

where we use the obvious fact
$\mathbf{a}+(\mathbf{b}-\mathbf{a})\odot\mathbf{p}\sim\mathcal{U}(\mathbf{a},\mathbf{b})$
if $\mathbf{p}\sim\mathcal{U}(\mathbf{0},\mathbf{1})$. In this way, we can
sample $\mathbf{p}$ from a fixed distribution to approximate the expectation,
compute $\nabla_{\mathbf{a}}K$ and $\nabla_{\mathbf{b}}K$ and perform a
gradient ascent step to improve $K$.

It would be beneficial for the optimization to converge to the global optimum
if the global optimum lies within $[\mathbf{a},\mathbf{b}]$. Therefore, it is
desirable that $\mathbf{a}$ and $\mathbf{b}$ are initialized outside the
feasible region. It is then necessary to introduce penalty terms

$P(\mathbf{a})=-\min\left(\sum\limits^{\circ}\left(\mathbf{a}-\mathbf{p}_{\min}\right),0\right)$
(6)

and

$Q(\mathbf{b})=\max\left(\sum\limits^{\circ}\left(\mathbf{b}-\mathbf{p}_{\max}\right),0\right),$
(7)

where $\mathbf{p}_{\min}$ and $\mathbf{p}_{\max}$ define the feasible region
of the power allocation. The penalty terms (6) and (7) are 0 if $\mathbf{a}$
or $\mathbf{b}$ are inside the feasible region, otherwise they are
proportional to the distances from the feasible region.

In one iteration, if $J(a)<\mathbb{E}_{\mathbf{p}\sim\mathcal{U}(a,b)}(J)$ and
$J(b)<\mathbb{E}_{\mathbf{p}\sim\mathcal{U}(a,b)}(J)$, as shown in Figure
3(a), a gradient ascent step will increase $a$ and decrease $b$. If
$J(a)<\mathbb{E}_{\mathbf{p}\sim\mathcal{U}(a,b)}(J)$ and
$J(b)=\mathbb{E}_{\mathbf{p}\sim\mathcal{U}(a,b)}(J)$, as shown in Figure
3(b), $\mathbf{a}$ increases and $\mathbf{b}$ is unchanged in a gradient
ascent step. If $J(a)=J(b)=\mathbb{E}_{\mathbf{p}\sim\mathcal{U}(a,b)}(J)$, as
shown in Figure 3(c), neither increasing $\mathbf{a}$ nor decreasing
$\mathbf{b}$ improves
$\mathbb{E}_{\mathbf{p}\sim\mathcal{U}(\mathbf{a},\mathbf{b})}(J)$. In this
case, $\mathbf{a}$ and $\mathbf{b}$ cannot converge if we use (5) as the
objective function. Therefore, we introduce an entropy regularization term in
the objective function as an additional incentive to reduce the entropy. For
multivariate uniform distribution with independent components, its entropy is
given as

$H(\mathbf{a},\mathbf{b})=\ln\left(\prod\limits^{\circ}(\mathbf{b}-\mathbf{a})\right).$
(8)

The total objective function $L$, including the penalties and the entropy
regularization, is

$\displaystyle L=$
$\displaystyle\mathbb{E}_{\mathbf{p}\sim\mathcal{U}(\mathbf{0},\mathbf{1})}J(\mathbf{a}+(\mathbf{b}-\mathbf{a})\odot\mathbf{p})$
(9) $\displaystyle-\epsilon P(\mathbf{a})-\epsilon Q(\mathbf{b})-\kappa
H(\mathbf{a},\mathbf{b}),$

where $\epsilon$ is the Lagrange multiplier for support region constraint and
$\kappa$ is the coefficient for the entropy regularization term, whose value
needs to be tuned during the optimization process, which is discussed in
Section III-B. If we compute the gradients of (9) with respect to $\mathbf{a}$
and $\mathbf{b}$, i.e., $\nabla_{\mathbf{a}}L$ and $\nabla_{\mathbf{b}}L$,
perform gradient ascents
$\mathbf{a}\leftarrow\mathbf{a}+r\cdot\nabla_{\mathbf{a}}L$ and
$\mathbf{b}\leftarrow\mathbf{b}+r\cdot\nabla_{\mathbf{b}}L$, where $r$ is the
learning rate, and repeat the process iteratively until the entropy of
$\mathcal{U}(\mathbf{a},\mathbf{b})$ is smaller than a predefined threshold
$H_{0}$, $\mathbf{a}$ and $\mathbf{b}$ have a high possibility to converge to
the global optimum.

$J$$p$$a$$b$$\mathbb{E}_{\mathbf{p}\sim\mathcal{U}(a,b)}(J)$ (a) In the
beginning of optimization

$J$$p$$a$$b$$\mathbb{E}_{\mathbf{p}\sim\mathcal{U}(a,b)}(J)$ (b) No entropy
regularization required

$J$$p$$a$$b$$\mathbb{E}_{\mathbf{p}\sim\mathcal{U}(a,b)}(J)$ (c) Entropy
regularization required

Figure 3: A non-convex objective function and its mean in the interval
$[a,b]$. (a): Both $J(a)$ and $J(b)$ are less than
$\mathbb{E}_{\mathbf{p}\sim\mathcal{U}(a,b)}(J)$, performing a gradient ascent
step of $\mathbb{E}_{\mathbf{p}\sim\mathcal{U}(a,b)}(J)$ increases $a$ and
decreases $b$. (b): After a few iterations,
$J(a)<\mathbb{E}_{\mathbf{p}\sim\mathcal{U}(a,b)}(J)$ and
$J(b)=\mathbb{E}_{\mathbf{p}\sim\mathcal{U}(a,b)}(J)$. A gradient ascent step
increases $a$ and does not change $b$. (c): Now we have
$J(a)=J(b)=\mathbb{E}_{\mathbf{p}\sim\mathcal{U}(a,b)}(J)$. Performing
gradient ascent changes neither $a$ nor $b$. We need the entropy
regularization term to help $a$ to “cross over” the local optimum.

######  Remark 1.

Note that the gradient of the first term in (9) (the integral of $J$) is $J$
itself. Therefore, whether $\mathcal{U}(\mathbf{a},\mathbf{b})$ reduces in a
gradient ascent step depends on the values $J(\mathbf{G},\mathbf{a})$ and
$J(\mathbf{G},\mathbf{b})$ themselves, not on their gradients. A flat global
optimum is therefore more advantageous than a sharp local optimum. In the case
of Figure 4,
$\nabla_{a}\mathbb{E}_{\mathbf{p}\sim\mathcal{U}(a,b)}(J)<\nabla_{b}\mathbb{E}_{\mathbf{p}\sim\mathcal{U}(a,b)}(J)$
because $J(a)<J(b)$. With a properly chosen $\kappa$, $\mathbf{a}$ can “jump”
over the sharp and poor local optimum and $\mathbf{b}$ stays in front of the
global optimum.

$J$$p$$a$$b$ Figure 4: An objective function with a flat global optimum and a
sharp local optimum.

###  III-B Tuning the Lagrange Multiplier for Entropy Term

From Remark 1, we realize that tuning $\kappa$ in (9) is crucial to the
success of the training because it has to be large enough such that the
training does not get stuck at local optima and it has to be small enough such
that the global optimum stays in the support of
$\mathcal{U}(\mathbf{a},\mathbf{b})$.

In the literature, it is suggested to optimize the factor of the entropy term
via the dual problem for a similar original problem [27]. However, since it is
very difficult to obtain the infimum of the Lagrangian of our complicated
objective function (unlike the canonical definition of the dual problem),
there is no guarantee of convexity of the dual problem, which would be a fatal
disadvantage in the considered problem. We propose a simple heuristic approach
to tune $\kappa$ in this section. Let $\kappa_{i}$ be the $\kappa$ in
iteration $i$ of the gradient ascent optimization. $\kappa_{i}=0$ for $i\leq
h$, where $h$ is a hyperparameter for constant $\kappa$. For $i>h$, $\kappa$
is computed as

$\kappa_{i+1}=\begin{cases}\kappa_{i}+\Delta\kappa&\textbf{if
}\sum_{j=1}^{h}H_{i-j}/h\leq H_{i}\\\
\max(0,\kappa_{t}-\Delta\kappa/2)&\text{otherwise,}\\\ \end{cases}$ (10)

where $\Delta\kappa$ is a constant small learning rate. The intuition behind
(10) is that we carefully increase $\kappa$ if the entropy does not reduce
(case 1) and decrease $\kappa$ otherwise such that $\kappa$ is not too big to
make $\mathbf{a}$ or $\mathbf{b}$ cross the global optimum (case 2).

######  Example.

In order to demonstrate the ability of the proposed objective function, we
apply the widely-used Rastrigin function for performance testing [34], which
is the overlapping of a parabolic function and cosine functions:

$f(\mathbf{x})=An+\sum_{i=1}^{n}(x_{i}^{2}-A\cos(2\pi x_{i})),$ (11)

where $A=10$ and $n$ is the number of dimensions. $f(\mathbf{x})$ is non-
convex with many local minima and its global minimum is $f(0,0,\dots,0)=0$. We
assume $n=10$ and use the proposed method (9) and conventional gradient
descent to minimize (11). $\kappa$ is initialized as 0 and penalty terms (6)
and (7) are omitted because the function is defined for all real numbers. The
results are shown in Figure 5, while the conventional gradient descent gets
stuck in a poor local minimum, the proposed method converges to the global
optimum successfully. This result demonstrates the ability of the proposed
method to find the global optimum of a simple non-convex objective function.

In principle, we can apply (9) to the power control problem (1) to obtain the
globally optimal transmit powers $\mathbf{p}$. However, its high computation
effort makes it difficult to be applied in real time (this fact can be
recognized by looking at the number of required iterations in Figure 5).
Therefore, we do not optimize $\mathbf{a}$ and $\mathbf{b}$ directly. Instead,
we propose to compute the parameters with neural networks and train the neural
networks with massive data. If the training data is general enough, the
trained neural network can return a good output given a new input that is not
part of the training data. Since the computation of the output given the input
in the neural network (the forward propagation) is simple and fast, we can
expect an efficient solution that can be applied in real time.

$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$1$$1.1$$1.2$$\cdot
10^{5}$$0$$100$$200$$300$$400$IterationFunction valueR-trickGradient
descent$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$1$$1.1$$1.2$$\cdot
10^{5}$$0$$2$$4$$6$$8$$\displaystyle\kappa$ Figure 5: Optimization of the
Rastrigin function with the proposed method and gradient descent. The
parameters of the Ackley function are set to $a=20$, $b=0.2$, $c=2\pi$, and
$d=6$.

Define $\mathbf{a}=\alpha_{\delta}(\mathbf{G})$ where $\alpha_{\delta}$ is the
mapping implemented by a neural network parameterized by $\delta$ and define
$\mathbf{\ell}=\mathbf{b}-\mathbf{a}=\beta_{o}(\mathbf{G})$ where $\beta_{o}$
is a neural network parameterized by $o$ (omicron). The objective is defined
as

$\displaystyle L_{\delta,o}(\mathbf{G})$
$\displaystyle=\mathbb{E}_{\mathbf{p}\sim\mathcal{U}(\mathbf{a},\mathbf{b})}\left(J(\mathbf{G},\mathbf{p})\right)$
(12) $\displaystyle\qquad-\epsilon P(\alpha_{\delta}(\mathbf{G}))-\epsilon
Q(\beta_{o}(\mathbf{G}))-\kappa
H(\alpha_{\delta}(\mathbf{G}),\alpha_{\delta}(\mathbf{G})+\beta_{o}(\mathbf{G}))$
$\displaystyle=\mathbb{E}_{\mathbf{p}\sim\mathcal{U}(\mathbf{0},\mathbf{1})}\left(J(\mathbf{G},\alpha_{\delta}(\mathbf{G})+\beta_{o}(\mathbf{G})\cdot\mathbf{p})\right)$
$\displaystyle\qquad-\epsilon P(\alpha_{\delta}(\mathbf{G}))-\epsilon
Q(\beta_{o}(\mathbf{G}))-\kappa
H(\alpha_{\delta}(\mathbf{G}),\alpha_{\delta}(\mathbf{G})+\beta_{o}(\mathbf{G})),$

where the second line is again the reparametrization trick. We maximize $L$ by
performing gradient ascent with many channel realizations, i.e.,

$\displaystyle\delta\leftarrow\delta+r\cdot\frac{1}{S}\sum_{s=1}^{S}\nabla_{\delta}L_{\delta,o}(\mathbf{G}_{i})$
(13)

and

$\displaystyle o\leftarrow
o+r\cdot\frac{1}{S}\sum_{s=1}^{S}\nabla_{o}L_{\delta,o}(\mathbf{G}_{i}),$ (14)

where $\mathbf{G}_{i}$ is one channel realization and $s=1,2,\dots,S$
represents a batch of data samples. When this optimization step is performed
sufficiently often, $\beta_{o}(\mathbf{G}_{i})$ will converge to a value close
to 0 due to the entropy regularization described in Section III-B and
$\alpha_{\delta}(\mathbf{G}_{i})$ will be a vector containing good transmit
powers for the channel matrix $\mathbf{G}_{i}$.

As introduced in Section II-B2, we expect that $\alpha_{\delta}$ and
$\beta_{o}$ are equivariant to permutation of user-BS pairs. In Section IV, we
introduce a permutation-equivariant neural network architecture
_InterferenceNet_ that achieves a high performance in the considered multi-
cell wireless network optimization problem.

##  IV InterferenceNet: a Permutation-Equivariant Neural Network Architecture
for Interference Network Optimization

In this section, we introduce the neural network architecture
_InterferenceNet_. Compared to the trivial neural networks with fully
connected layers, the design of InterferenceNet takes the domain knowledge of
cellular network optimization into account and is permutation-equivariant,
which reflects the physical property of the permutation-equivariance of the
cellular network and improves the generalizability of the model.

###  IV-A Network Architecture

The InterferenceNet has $L$ layers. The channel from user $j$ to BS that
serves user $i$ is represented by a feature vector $\mathbf{f}_{ij,l}$ as the
input of layer $l$. In particular, $\mathbf{f}_{ij,1}=\log_{10}(g_{ij})$,
where $g_{ij}$ is the channel gain as defined in (1).

For each channel from user $j$ to BS $i$, denoted as $c_{ij}$, we define the
following four categories of channels according to their roles in the system
performance: Category 1 contains channel $c_{ij}$ itself, i.e.,
$\mathcal{C}^{1}(i,j)=\\{(i,j)\\}$, where $\mathcal{C}^{1}(i,j)$ is the set of
channels in category 1 for channel from user $j$ to BS that serves user $i$.
Category 2 contains channels with user $j$ but different BSs, i.e.,
$\mathcal{C}^{2}(i,j)=\\{(k,j),k=1,\dots,N,k\neq i\\}$. These are channels on
whose receivers user $j$ contributes interference to. Category 3 contains
channels with BS $i$ but different users, i.e.,
$\mathcal{C}^{3}(i,j)=\\{(i,k),k=1,\dots,N,k\neq j\\}$. These are channels
that contribute to the interference of user $j$. Category 4 contains channels
with different users and different BSs, i.e.,
$\mathcal{C}^{4}(i,j)=\\{(k,p),k,p=1,\dots,N,k\neq i,p\neq j\\}$. These are
channels that do not contribute to the performance of link from user $j$ to BS
$i$. These categories are illustrated in Figure 6.

122334444$i$$j$User (transmitter)BS (receiver) Figure 6: Illustration of 4
categories of channels for the channel $c_{ij}$: Category 1 is $c_{ij}$
itself. Category 2 contains channels sharing the same user with $c_{ij}$ but
with different BSs. Category 3 contains channels sharing the same BS with
$c_{ij}$ but with different users. Category 4 are the remaining channels.

For $l<L$ (i.e., for layers before the last layer), the output feature vector
of channel $c_{ij}$ in layer $l$ (i.e., the input feature in layer $l+1$) in
category $k$ is computed as

$\mathbf{f}^{c}_{ij,l+1}=\sum_{(k,p)\in\mathcal{C}^{k}(i,j)}\text{ReLU}(\mathbf{W}_{l}^{k}\mathbf{f}_{kp,l}+\mathbf{b}_{l}^{k})/|\mathcal{C}^{k}(i,j)|,$
(15)

where $\mathbf{W}^{k}_{l}$ and $\mathbf{b}^{k}_{l}$ are the trainable weights
and bias of layer $l$ for category $k$, respectively. The reason we use the
mean instead of the sum of
$\text{ReLU}(\mathbf{W}_{l}^{k}\mathbf{f}_{kp,l}+\mathbf{b}_{l}^{k})$ for
$(k,p)\in\mathcal{C}^{k}(i,j)$ is that different categories have different
numbers of channels. For example, if we have 7 users, category 1 has 1 channel
and category 4 has 36 channels. If we use the sum instead of the mean, the
variance of features in category 4 would be much higher than the variance of
features in category 1. This is numerically difficult for the neural network
to learn. Therefore, we use the mean instead of the sum to make the features
of different categories roughly identically distributed.

The total output feature vector of channel $c_{ij}$ in layer $l$ is the
concatenation of channel gain and the feature vectors in all four categories,
i.e.,

$\mathbf{f}_{ij,l+1}=\left(\log_{10}(g_{ij}),\mathbf{f}^{1}_{ij,l+1},\mathbf{f}^{2}_{ij,l+1},\mathbf{f}^{3}_{ij,l+1},\mathbf{f}^{4}_{ij,l+1}\right).$
(16)

For $l<L$, $\mathbf{f}_{ij,l+1}$ is the input for the next layer.

For $l=L$ (i.e., for the last layer), the output of channel $c_{jj}$ is
computed as

$f_{jj,L+1}=\mathbf{w}_{L}\mathbf{f}_{jj,L}+b_{L},$ (17)

where the output $f_{jj,L+1}$ is a scalar instead of a vector (unlike the
previous layers), $\mathbf{w}_{L}$ and $b_{L}$ are trainable weights and bias
of layer $L$, respectively. Note that only the diagonal elements are processed
in the last layer (i.e., only $c_{jj}$ is processed and $c_{ij}$ is omitted
for $i\neq j$) since the output per user-BS pair is a scalar, therefore the
indices are $f_{jj,L+1}$ instead of $f_{ij,L+1}$, which is the action of user
$j$.

We can organize the input features of layer $l$ of all channels as a three-
dimensional tensor $\mathbf{F}_{l}$ of shape
$(\text{dim}(\mathbf{f}_{ij,l}),I,I)$, where
$\mathbf{F}_{l}[\cdot,i,j]=\mathbf{f}_{ij,l}$. For $l<L$, the output feature
tensors in categories $1$-$4$ are computed by

$\displaystyle\mathbf{F}_{l+1}^{1}$
$\displaystyle=\operatorname{ReLU}(\mathbf{W}_{l}^{1}\bullet\mathbf{F}_{l}+\mathbf{b}_{l}^{1}\bullet\mathbf{1}_{1\times
I\times I})$ (18) $\displaystyle\mathbf{F}_{l+1}^{2}$
$\displaystyle=\mathbf{E}\circ\operatorname{ReLU}(\mathbf{W}_{l}^{2}\bullet\mathbf{F}_{l}+\mathbf{b}_{l}^{2}\bullet\mathbf{1}_{1\times
I\times I})/(I-1)$ (19) $\displaystyle\mathbf{F}_{l+1}^{3}$
$\displaystyle=\operatorname{ReLU}(\mathbf{W}_{l}^{3}\bullet\mathbf{F}_{l}+\mathbf{b}_{l}^{3}\bullet\mathbf{1}_{1\times
I\times I})\circ\mathbf{E}/(I-1)$ (20) $\displaystyle\mathbf{F}_{l+1}^{4}$
$\displaystyle=\mathbf{E}\circ\operatorname{ReLU}(\mathbf{W}_{l}^{4}\bullet\mathbf{F}_{l}+\mathbf{b}_{l}^{4}\bullet\mathbf{1}_{1\times
I\times I})\circ\mathbf{E}/((I-1)^{2})\,,$ (21)

respectively. The final output $\mathbf{F}_{l+1}$ is the concatenation of
$\log_{10}(\mathbf{G})$, $\mathbf{F}_{l+1}^{1}$, $\mathbf{F}_{l+1}^{2}$,
$\mathbf{F}_{l+1}^{3}$, and $\mathbf{F}_{l+1}^{4}$ along the first dimension.
In this way, the computation in the neural network can be performed
efficiently with tensor products.

In the last layer, the raw output is computed as

$\mathbf{f}_{L+1}=\mathbf{w}_{L}\mathrm{diag}(\mathbf{F}_{L})+b_{L}\mathbf{1}_{1\times
N}.$ (22)

For $\alpha_{\delta}$, the final layer does not have an activation function.
For $\beta_{o}$, a special activation function is necessary since a uniform
distribution $\mathcal{U}(\mathbf{a},\mathbf{a}+\ell)$ can only be defined
when every dimension of $\ell$ is positive, output of $\beta_{o}$ is defined
as

$\ell=\max(\mathbf{f}_{L+1}^{\beta},\ell_{\min}),$ (23)

where $\mathbf{f}_{L+1}^{\beta}$ is the raw output of the last layer of
$\beta_{o}$ and $\ell_{\min}$ is a very small number (e.g., $10^{-6}$). For
all other layers in both $\alpha_{\beta}$ and $\beta_{o}$, the activation
function is ReLU.

The above-described information processing is illustrated in Figure 7.

ReceiverChannelTransmitterLayer $l$ for category 1Layer $l$ for category
2Layer $l$ for category 3Layer $l$ for category 4ReceiverChannelCategory
1Category 2Category 3Category 4Transmitter

(a) First layer

ReceiverChannelCategory 1Category 2Category 3Category 4TransmitterLayer $l$
for category 1Layer $l$ for category 2Layer $l$ for category 3Layer $l$ for
category 4ReceiverChannelCategory 1Category 2Category 3Category 4Transmitter

(b) Intermediate layers

ReceiverChannelCategory 1Category 2Category 3Category 4TransmitterLayer
$L$ReceiverTransmit powerTransmitter

(c) Final layer

Figure 7: Illustration of information processing from first over intermediate
to final output layers.

###  IV-B Proof of Permutation-Equivariance

The proposed InterferenceNet architecture is permutation-equivariant
intuitively because the above-defined 4 categories do not change by any
permutation and the summation within a category is permutation-invariant.
Mathematically, the permutation-equivariance is proved as follows.

######  Theorem 1.

The proposed InterferenceNet $\Phi$ is permutation-equivariant, i.e.,

$\mathbf{p}\mathbf{P}^{T}=\Phi(\mathbf{P}\mathbf{G}\mathbf{P}^{T})$ (24)

holds if $\mathbf{p}=\Phi(\mathbf{G})$ for any permutation matrix
$\mathbf{P}$.

###### Proof.

See Appendix A. ∎

###  IV-C Equivalence to Graphical Neural Network and Distributed
Implementation

The GNN is a novel neural network architecture that performs information
processing for multiple nodes connected by (directional) edges [35]. A simple
example of a GNN is shown in Figure 8.

$\mathbf{f}_{i}^{(l)}$$\mathbf{f}_{j}^{(l)}$$\mathbf{f}_{k}^{(l)}$$\mathbf{e}_{j,i}$
Figure 8: A simple graph on which a GNN is defined.

The node feature of node $i$ on layer $l$ is computed as [36]

$\mathbf{f}_{i}^{(l)}=\gamma^{(l)}\left(\mathbf{f}_{i}^{(l-1)},\square_{j\in\mathcal{N}(i)}\phi^{(l)}\left(\mathbf{f}_{i}^{(l-1)},\mathbf{f}_{j}^{(l-1)},\mathbf{e}_{j,i}\right)\right),$
(25)

where $\mathbf{f}_{i}^{(l)}$ is the node feature of node $i$ on layer $l$,
$\mathbf{e}_{j,i}$ is the edge feature of the directional edge from node $j$
to node $i$, $\phi^{(l)}$ is a neural network that is applied on each edge and
performs message passing (in the example in (25), the edge is from node $j$ to
node $i$, therefore its input is the concatenation of
$\mathbf{f}_{j}^{(l-1)}$, $\mathbf{e}_{j,i}$ and $\mathbf{f}_{j}^{(l-1)}$),
$\square$ is a symmetric, permutation-invariant function (i.e., the function
value is constant if the arguments are permuted, for example,
$\square(x_{1},x_{2},\dots,x_{n})=\square(x_{2},x_{1},\dots,x_{n})$), such as
sum or max, $\gamma^{(l)}$ is a neural network that is applied to each node on
layer $l$ and processes the information of the node itself and the messages
from its neighbors. In the example in (25), the node is node $i$, therefore
its input is the concatenation of $\mathbf{f}_{i}^{(l-1)}$ and
$\square_{j\in\mathcal{N}(i)}\phi^{(l)}\left(\mathbf{f}_{i}^{(l-1)},\mathbf{f}_{j}^{(l-1)},\mathbf{e}_{j,i}\right)$.
The distributed nature of the GNN and its ability to perform message passing
and information aggregation make it a competitive candidate to optimize
cellular networks [31, 32, 33].

Although the InterferenceNet architecture is introduced as a centralized
neural network, we can understand each row in Figure 6 as a node and easily
prove that it can be reshaped to fit into the form of (25). Based on (18) -
(21), we define $\phi$, $\square$ and $\gamma$ as

$\displaystyle\phi^{(l)}(\mathbf{F}_{l}[\cdot,i,\cdot])=\left(\begin{array}[]{l}\text{ReLU}(\mathbf{W}_{l}^{2}\bullet\mathbf{F}_{l}[\cdot,i,\cdot]+\mathbf{b}_{l}^{2}\bullet\mathbf{1}_{1\times
I\times I})\\\
\text{ReLU}(\mathbf{W}_{l}^{4}\bullet\mathbf{F}_{l}[\cdot,i,\cdot]+\mathbf{b}_{l}^{4}\bullet\mathbf{1}_{1\times
I\times I})\circ\mathbf{E}/(I-1)\end{array}\right),$ (26)

$\square_{j\in\mathcal{N}(i)}\phi^{(l)}(\mathbf{F}_{l}[\cdot,j,\cdot])=\sum_{j}\phi^{(l)}(\mathbf{F}_{l}[\cdot,j,\cdot])/(I-1),$
(27)

and

$\displaystyle\gamma^{(l)}(\mathbf{F}_{l}[\cdot,i,\cdot],\square_{j\in\mathcal{N}(i)}\phi^{(l)}(\mathbf{F}_{l}[\cdot,j,\cdot]))=$
$\displaystyle\left(\begin{array}[]{l}\mathbf{G}\\\
\text{ReLU}(\mathbf{W}_{l}^{1}\bullet\mathbf{F}_{l}[\cdot,i,\cdot]+\mathbf{b}_{l}^{1}\bullet\mathbf{1}_{1\times
I\times I})\\\
\mathbf{E}\circ\text{ReLU}(\mathbf{W}_{l}^{3}\bullet\mathbf{F}_{l}[\cdot,i,\cdot]+\mathbf{b}_{l}^{3}\bullet\mathbf{1}_{1\times
I\times I})/(I-1)\\\
\square_{j\in\mathcal{N}(i)}\phi^{(l)}(\mathbf{F}_{l}[\cdot,j,\cdot])\end{array}\right),$
(28)

respectively, where the concatenation in (26) and (28) is along the first
dimension, we conclude that the computation in one layer is presented in the
form of (25). This fact proves that the proposed InterferenceNet architecture
is equivalent to the GNN, thus inherits its advantages such as distributed
implementation in each BS and limited data amount of message passing in the
fronthaul.

###  IV-D Complexity

Besides the consideration of domain knowledge in the neural network
architecture design and the permutation-equivariance, another advantage of the
proposed method is the low complexity due to the small number of trainable
parameters. Compared to the conventional deep neural network (DNN) with fully
connected layers proposed in [3], which use different weights and biases to
process information of different channels, all channels use the same filters
((18) - (21)) to process the information in the proposed InterferenceNet
architecture. Therefore, the number of trainable neural network parameters is
independent from the number of user-BS pairs. In the proposed architecture,
there are 39844 trainable parameters. Compared to it, the neural network with
fully connected layers proposed in [3] has 8708189 trainable parameters, which
is approximately 218 times more than the InterferenceNet architecture. The
less number of trainable parameters results in a more compact neural network,
less computational effort in training and less risk of overfitting.

Summarizing the above sections, the proposed algorithm to solve the non-convex
power control problem is formulated as Algorithm 1.

Algorithm 1 Neural network training with a non-convex objective with
stochastic action and reparameterization trick

Formulate the objective function according to (12).

Initialize $\alpha_{\delta}$ and $\beta_{o}$ such that
$\alpha_{\delta}(\mathbf{G})<0$ and $\beta_{o}(\mathbf{G})>1$ for all possible
$\mathbf{G}$.

Initialize $\kappa=0$ and $s_{0}=p_{\max}$.

while $H(\alpha_{\delta}(\mathbf{G}),\beta_{o}(\mathbf{G}))>H_{0}$ do

Perform a gradient ascent step of (9) w.r.t. $\delta$ and $o$.

Update $\kappa$ according to (10).

end while

##  V Training and Testing Results

In this section, we will evaluate the introduced InterferenceNet architecture
on a multi-user communication scenario.

###  V-A Scenario and Channel Models

We consider a multi-cell mobile network scenario with four BSs, as shown in
Figure 1. Coordinates of the four BSs are $(0.5,0.5)$, $(0.5,1.5)$,
$(1.5,1.5)$ and $(0.5,0.5)$, with the unit of $\mathrm{km}$, and each BS is
equipped with $n_{R}=2$ antennas. We consider $I=7$ and $I=4$ single-antenna
users who share the same physical resource block (PRB) and thus have
interference on each other. Please note that the total number of users in this
cellular network can be more than $I$ because a cellular network can have more
than one PRB. Users allocated with different PRBs have negligible interference
on each other because their signals are well separated in time or frequency
domain. Furthermore, the transmit power control problems of different PRBs can
be considered decoupled from each other and can be solved separately. The
network can therefore serve many more users while we focus on the difficult
part: users sharing the same PRB and causing interference to each other. The
same power control strategy can be applied to every PRB individually.

Assumptions on channel models and parameters are in line with [3] in general.
In particular, the propagation path-loss is computed with the model presented
in [37]. The carrier frequency is assumed to be $1.8\text{\,}\mathrm{GHz}$
whereas the power decay factor is $4.5$. The fast fading is modeled as i.i.d.
complex random variables, where the amplitude is drawn from the standard
normal distribution and the complex phase is uniformly distributed between $0$
and $2\pi$. The static power consumption is assumed to be
$1\text{\,}\mathrm{W}$ and the power amplifier inefficiency $\mu_{i}=4$ for
all $i$. Noise power at the receiver is computed as $\sigma^{2}=FN_{0}B$,
where $F=$3\text{\,}\mathrm{dB}$$ is the noise figure at the receiver,
$N_{0}=$-174\text{\,}\mathrm{dBm}\text{/}\mathrm{Hz}$$ is the noise spectrum
density, $B=$180\text{\,}\mathrm{kHz}$$ is the bandwidth. All users have the
same maximum transmit power $p_{\max}$.

The InterferenceNet is implemented with PyTorch 1.11 [38], in which the
reparameterization trick is implemented as the _rsample_ method. There are
$L=5$ layers in the InterferenceNet and the feature dimension is 20 for each
category and each layer. The ADAM [39] optimizer is chosen to train the neural
network and the learning rate is set to be $2\cdot 10^{-7}$ for
$p_{\max}>$0\text{\,}\mathrm{dBm}$$ and $10^{-7}$ otherwise. The training set
contains 6000 samples whereas the testing set has 1000 samples. The batch size
is 512.

###  V-B Training and Testing Results

$0$$0.5$$1$$1.5$$2$$2.5$$3$$3.5$$4$$4.5$$5$$\cdot 10^{4}$$0$$1$$2$EpochEE
[Mbit/Joule]SGAR-trick Figure 9: Improvement of objective with training with 4
users using SGA and reparameterization trick.

As mentioned above, we consider 4 and 7 users. Figure 9 shows the improvement
of objective using SGA and reparameterization trick with 4 users. The maximum
transmit power is $5\text{\,}\mathrm{d}\mathrm{B}\mathrm{m}$. Although the EE
optimization problem is non-convex, SGA is able to obtain the same EE as the
reparameterization trick, indicating that both methods find the same solution.
In the following text, we will therefore focus on the results with 7 users.

Figure 10 shows the curves of achieved EE with 7 users during training. The
blue, orange and green curves are achieved EE with
$p_{\max}=$-20\text{\,}\mathrm{dBm}$$ and reparameterization trick,
$p_{\max}=$5\text{\,}\mathrm{dBm}$$ and reparameterization trick, and
$p_{\max}=$5\text{\,}\mathrm{dBm}$$ and conventional SGA, respectively. All
three setups use the same neural network architecture and hyperparameters
except the learning rate (see Section V-A). When the maximum transmit power is
low ($-20\text{\,}\mathrm{dBm}$), the optimal transmit power is the maximum
transmit power $p_{\max}$ and the achieved EE is limited by it. When the
maximum transmit power is high ($5\text{\,}\mathrm{dBm}$), the optimal
transmit power lies between 0 and $p_{\max}$ and the achieved EE is higher
than with a low maximum transmit power. Unlike the four-user scenario, if we
use the conventional SGA instead of the objective function (12), it converges
to a poor local optimum much earlier than training with the reparameterization
trick, as the green curve shows. This result shows the advantage of the
proposed solution compared to the conventional SGA for non-convex objectives.

$0$$0.2$$0.4$$0.6$$0.8$$1$$1.2$$1.4$$\cdot 10^{4}$$0$$1$$2$$3$EpochEE
[Mbit/Joule]$p_{\max}=$-20\text{\,}\mathrm{dBm}$$,
r-trick$p_{\max}=$5\text{\,}\mathrm{dBm}$$,
r-trick$p_{\max}=$5\text{\,}\mathrm{dBm}$$, SGA Figure 10: Improvement of
objective with training over 15000 thousand epochs for three cases with 7
users.

In order to gain more insights into the training process, Figure 11 shows the
penalty introduced by (6) and (7). We can observe that the penalty reduces to
0 quickly at the beginning of the training, suggesting that $\mathbf{a}$ and
$\mathbf{b}$ lies within $[0,p_{\max}]$ for all $\mathbf{G}$ after training of
some epoches. Please note that it is important to initialize $\mathbf{a}$ and
$\mathbf{b}$ outside the feasible region because the global optimum cannot be
considered in (12) if it lies outside $[\mathbf{a},\mathbf{b}]$ (an example
is, as shown above, the optimal transmit power is $p_{\max}$ for a low maximum
transmit power). Since it is very difficult to initialize $\mathbf{a}$ and
$\mathbf{b}$ exactly on the boundaries, it is easier to initialize them
outside the feasible region and use the penalty terms to reduce them to the
feasible region at the beginning of the training.

$0$$0.2$$0.4$$0.6$$0.8$$1$$1.2$$1.4$$\cdot
10^{4}$$0$$2$$4$EpochPenalty$p_{\max}=$-20\text{\,}\mathrm{dBm}$$$p_{\max}=$5\text{\,}\mathrm{dBm}$$
Figure 11: Value of penalty during training.

Figure 12 illustrates the mean entropy during training for the two
aforementioned maximum transmit powers. Since $\mathbf{a}$ and $\ell$ are
initialized in such a way that the initial transmit power distribution covers
the whole interval $[0,p_{\max}]$ for any $\mathbf{G}$. As a result, the
entropy is high at the beginning of the training. It begins to decrease
quickly first due to the penalty terms (6) and (7) because the penalty term is
independent from $\mathbf{G}$. Since $p_{\max}$ is assumed to be the same for
all users, the penalty term is the same for all samples. This results in the
high decreasing speed of entropy. After $\mathbf{a}$ and $\mathbf{a}+\ell$ lie
within $[0,p_{\max}]$, the entropy decreases more slowly because the objective
function is different given different $\mathbf{G}$. The entropy reduces to
roughly $-80$ at the end of training due to the lower bound of $\ell$ posed in
(23). Training with higher maximum transmit power converges more slowly due to
the larger and more complex feasible region.

$0$$0.2$$0.4$$0.6$$0.8$$1$$1.2$$1.4$$\cdot
10^{4}$$-80$$-60$$-40$$-20$$0$EpochEntropy$p_{\max}=$-20\text{\,}\mathrm{dBm}$$$p_{\max}=$5\text{\,}\mathrm{dBm}$$
Figure 12: Value of entropy of action distribution during training.

Figure 13 shows the mean $\kappa$ during training. As described in Section
III-B, $\kappa$ is initialized to be 0 and is increased if the entropy does
not decrease in the previous $h$ epochs. We can observe that the mean $\kappa$
is 0 at first and begins to increase after about $4000$ epochs, indicating
that the entropy of the distribution does not decrease for some samples (as
shown in Figure 3(c)). Comparing Figure 12 and Figure 13, we can see that the
mean entropy is far above the final entropy as $\kappa$ begins to increase
(after about 400 epochs), indicating that the increasing $\kappa$ is due to
multiple local optima rather than the lower bound of the distribution (23). An
increased $\kappa$ helps to overcome the suboptimal local optima and to
further decrease the entropy. Due to the lower bound on $\ell$ (23), the
entropy stops decreasing despite the increasing $\kappa$ in the final phase of
training. This happens after roughly $8000$ epochs for
$p_{\max}=$-20\text{\,}\mathrm{dBm}$$ and $15\,000$ epochs for
$p_{\max}=$5\text{\,}\mathrm{dBm}$$. The increasing speed of mean $\kappa$ is
slower for $p_{\max}=$5\text{\,}\mathrm{dBm}$$ because the convergence speed
with the high maximum transmit power is slower (see Figure 12).

$0$$0.2$$0.4$$0.6$$0.8$$1$$1.2$$1.4$$\cdot 10^{4}$$0$$0.5$$1$$\cdot
10^{-2}$EpochMean
$\displaystyle\kappa$$p_{\max}=$-20\text{\,}\mathrm{dBm}$$$p_{\max}=$5\text{\,}\mathrm{dBm}$$
Figure 13: Value of $\kappa$ during training.

Above text shows details and unveils insights of the training. In Figure 14,
we show the achieved EE with the testing data (i.e., data that have not been
used in the training phase) and compare it to the global optimum computed with
the branch-and-bound algorithm [3] for different maximum transmit powers and
different numbers of users. For both setups, the optimal transmit power is the
maximum transmit power when it is low (therefore the global optimum is easier
to find) and it lies between $[0,p_{\max}]$ for higher maximum transmit power
(therefore the EE stops increasing with further increasing $p_{\max}$), which
makes the problem more challenging. From the figures we observe that the
proposed method achieves good results close to the global optimum for all
setups.

$-30$$-25$$-20$$-15$$-10$$-5$$0$$5$$10$$15$$20$$0$$1$$2$$3$$4$$\displaystyle
p_{\max}$ [dBW]EE [Mbit/Joule]Global optimum, 4 usersTesting result, 4
usersGlobal optimum, 7 usersTesting result, 7 users Figure 14: Comparison
between global optimum and achieved performance.

##  VI Conclusion

The EE is an important objective of the future wireless communication systems.
Due to the non-convexity of the objective, this problem has multiple local
optima and cannot be solved with conventional convex optimization tools.
Although the global optimum can be found with the branch-and-bound algorithm,
this solution does not scale well due to the high complexity. This paper
presents a machine learning solution to the problem. We first introduce a
reparameterization trick based objective to solve the non-convexity problem.
Instead of optimizing an operation point of transmit powers, we apply
stochastic action, let the transmit powers uniformly distributed in the
interval $[\mathbf{a},\mathbf{b}]$ and define the objective as the expected EE
given the transmit power distribution, where we approximate the expectation
with the mean of random variables and use the reparameterization trick to make
the random variable differentiable. The interval is initialized to cover the
whole feasible region. By optimizing the interval boundaries $\mathbf{a}$ and
$\mathbf{b}$ with SGA, the interval reduces gradually and converges to the
global optimum. In case that the interval stops shrinking due to multiple
strong local optima, we introduce the entropy regularization term to provide
additional incentive for convergence. The second contribution of this paper is
a dedicated neural network architecture _InterferenceNet_ , which encodes the
domain knowledge of interference channels and realizes permutation-
equivariance (i.e., when the BS-user pairs are permuted, the output of the
neural network is permuted in the same way automatically). The InterferenceNet
defines four categories of channels for each user-BS pair: channel of the
considered transmitter and receiver, channels of other transmitters and
considered receiver (interference sources), channels of the considered
transmitter and other receivers (interference on other users) and channels of
different transmitters and receivers (channels not directly related to the
considered user-BS pair). It is also proven that the proposed architecture is
equivalent to the GNN and therefore it allows for a distributed
implementation. Training and testing results show solid improvement behavior.
Compared to the state-of-the-art, the achieved EE is close to the global
optimum obtained with the brand-and-bound algorithm. Possible extensions of
this work are mainly in two directions: on one hand, we can apply the proposed
method to other objectives with additional constraints, such as sum-rate
subject to minimum rate requirement. On the other hand, we can extend the
proposed InterferenceNet architecture to more complicated systems, e.g., cell-
free massive MIMO and multihop communication systems. Source code and data
will be available for public once the paper is accepted.

## Acknowledgement

The authors would like to thank Mr. F. Siegismund-Poschmann for his suggestion
on the test function and Mr. R. Wang for helping the training and testing
process.

##  Appendix A Proof of Theorem 1

We first extend the permutation of matrix (2) to permutation of tensor as

$\mathbf{F}^{\prime}=\mathbf{P}\circ\mathbf{F}\circ\mathbf{P}^{T},$ (29)

i.e., the second and third dimensions (corresponding to the rows and columns
of a matrix) are permuted by $\mathbf{P}$ simultaneously. Then we prove that
the operation of layer $l$ is permutation-equivairant, i.e.,
$\mathbf{P}\circ\mathbf{F}_{l+1}\circ\mathbf{P}^{T}=\Phi_{l}(\mathbf{P}\circ\mathbf{F}_{l}\circ\mathbf{P}^{T})$
holds if $\mathbf{F}_{l+1}=\Phi_{l}(\mathbf{F}_{l})$ for $l<L$, where
$\Phi_{l}$ is the information processing in layer $l$, and
$\mathbf{p}\mathbf{P}^{T}=\Phi_{L}(\mathbf{P}\circ\mathbf{F}_{L}\circ\mathbf{P}^{T})$
holds if $\mathbf{p}=\Phi_{L}(\mathbf{F}_{L})$. Finally, we prove the whole
InterferenceNet is also equivariant.

Combining (18) and (29), the permuted output for category 1 is computed as

$\displaystyle(\mathbf{F}_{l+1}^{1})^{\prime}$
$\displaystyle=\text{ReLU}(\mathbf{W}_{l}^{1}\bullet\mathbf{P}\circ\mathbf{F}_{l}\circ\mathbf{P}^{T}+\mathbf{b}_{l}^{1}\bullet\mathbf{1}_{1\times
N\times N})$ (30)
$\displaystyle=\text{ReLU}(\mathbf{P}\circ\mathbf{W}_{l}^{1}\bullet\mathbf{F}_{l}\circ\mathbf{P}^{T}+\mathbf{P}\circ\mathbf{b}_{l}^{1}\bullet\mathbf{1}_{1\times
N\times N}\circ\mathbf{P}^{T})$
$\displaystyle=\mathbf{P}\circ\text{ReLU}(\mathbf{W}_{l}^{1}\bullet\mathbf{F}_{l}+\mathbf{b}_{l}^{1}\bullet\mathbf{1}_{1\times
N\times N})\circ\mathbf{P}^{T}$
$\displaystyle=\mathbf{P}\circ\mathbf{F}_{l+1}^{1}\circ\mathbf{P}^{T}.$

In the second line of (30), exchanging order of $\mathbf{P}$ and
$\mathbf{W}_{l}^{1}$ is allowed because the multiplication between
$\mathbf{W}_{l}^{1}$ and $\mathbf{F}_{l}$ is in the first dimension and does
not affect the permutation in the second and third dimension. Multiplying
$\mathbf{P}$ and $\mathbf{P}^{T}$ on the left and right hand sides of
$\mathbf{b}_{l}^{1}\bullet\mathbf{1}_{1\times I\times I})$ in the second and
third dimensions does not change $\mathbf{b}_{l}^{1}\bullet\mathbf{1}_{1\times
I\times I}$ because this action permutes its second and third dimension, which
are identical. In the third line of (30), we move the permutation out of the
summation and ReLU operation because the ReLU is an elementwise operation and
does not depend on the order of elements in the tensor. Therefore, the
operation of a layer on category 1 is permutation-equivariant.

In order to prove the permutation-equivariance of the other three categories,
we first prove that the product between $\mathbf{E}$ and any permutation
matrix $\mathbf{P}$ is commutative, which is straightforward:

$\displaystyle\mathbf{E}\mathbf{P}$
$\displaystyle=\mathbf{1}\mathbf{P}-\mathbf{I}\mathbf{P}$ (31)
$\displaystyle=\mathbf{P}\mathbf{1}-\mathbf{P}\mathbf{I}$
$\displaystyle=\mathbf{P}\mathbf{E},$

where we use the property
$\mathbf{1}\mathbf{P}=\mathbf{P}\mathbf{1}=\mathbf{1}$ because the sum of rows
and columns of $\mathbf{P}$ is 1.

Combining (19) and (29), the permuted output for category 2 is computed as

$\displaystyle(\mathbf{F}_{l+1}^{2})^{\prime}=$
$\displaystyle\mathbf{E}\circ\text{ReLU}(\mathbf{W}_{l}^{2}\bullet\mathbf{P}\circ\mathbf{F}_{l}\circ\mathbf{P}^{T}+\mathbf{b}_{l}^{2}\bullet\mathbf{1}_{1\times
N\times N})/(2(N-1))$ (32) $\displaystyle=$
$\displaystyle\mathbf{E}\circ\mathbf{P}\circ\text{ReLU}(\mathbf{W}_{l}^{2}\bullet\mathbf{F}_{l}+\mathbf{b}_{l}^{2}\bullet\mathbf{1}_{1\times
N\times N})\circ\mathbf{P}^{T}/(2(N-1))$ $\displaystyle=$
$\displaystyle\mathbf{P}\circ\mathbf{E}\circ\text{ReLU}(\mathbf{W}_{l}^{2}\bullet\mathbf{F}_{l}+\mathbf{b}_{l}^{2}\bullet\mathbf{1}_{1\times
N\times N})\circ\mathbf{P}^{T}/(2(N-1))$ $\displaystyle=$
$\displaystyle\mathbf{P}\circ\mathbf{F}_{l+1}^{2}\circ\mathbf{P}^{T}.$

In the second line of (32), we reuse the result of (30). In the third line of
(32), we use the fact proved above that the product between $\mathbf{E}$ and
$\mathbf{P}$ is commutative. Proofs of equivariance for categories 3 and 4 are
similar to proof for category 2 and are omitted. Since the operation of a
layer is equivariant for all four categories, the operation is equivariant for
the total feature as well, i.e.,
$\mathbf{P}\mathbf{F}_{l+1}\mathbf{P}^{T}=\Phi_{l}(\mathbf{P}\mathbf{F}_{l}\mathbf{P}^{T})$
holds if $\mathbf{F}_{l+1}=\Phi_{l}(\mathbf{F}_{l})$ for $l<L$.

For the last layer, the proof of equivariance is similar to the previous
layers because the additional diagonalization operation is also equivariant.

Having proved all layers in InterferenceNet is equivariant, it is easy to
prove that the InterferenceNet itself is also equivariant by recursively
pulling $\mathbf{P}$ out of the operators:

$\displaystyle\mathbf{p}^{\prime}=$
$\displaystyle\Phi_{L}(\dots\Phi_{2}(\Phi_{1}(\mathbf{P}\circ\mathbf{G}\circ\mathbf{P}^{T}))\dots)$
(33) $\displaystyle=$
$\displaystyle\Phi_{L}(\dots\Phi_{2}(\mathbf{P}\circ\Phi_{1}(\mathbf{G})\circ\mathbf{P}^{T})\dots)$
$\displaystyle=$ $\displaystyle\dots$ $\displaystyle=$
$\displaystyle\Phi_{L}(\mathbf{P}\circ\Phi_{L-1}(\dots\Phi_{1}(\mathbf{G})\dots)\circ\mathbf{P}^{T})$
$\displaystyle=$
$\displaystyle\Phi_{L}(\Phi_{L-1}(\dots\Phi_{1}(\mathbf{G})\dots))\mathbf{P}^{T}$
$\displaystyle=$ $\displaystyle\mathbf{p}\mathbf{P}^{T}.$

This concludes the proof.

## References

  * [1] Mung Chiang, Prashanth Hande, Tian Lan and Chee Wei Tan  “Power control in wireless cellular networks”  In _Foundations and Trends® in Networking_ 2.4 Now Publishers, Inc., 2008, pp. 381–533
  * [2] Honghai Zhang et al.  “Weighted sum-rate maximization in multi-cell networks via coordinated scheduling and discrete power control”  In _IEEE Journal on Selected Areas in Communications_ 29.6 IEEE, 2011, pp. 1214–1224
  * [3] Bho Matthiesen et al.  “A globally optimal energy-efficient power control framework and its efficient implementation in wireless interference networks”  In _IEEE Transactions on Signal Processing_ 68 IEEE, 2020, pp. 3887–3902  DOI: 10.1109/TSP.2020.3000328
  * [4] Amin Ghazanfari, Hei Victor Cheng, Emil Björnson and Erik G Larsson  “Enhanced fairness and scalability of power control schemes in multi-cell massive MIMO”  In _IEEE Transactions on Communications_ 68.5 IEEE, 2020, pp. 2878–2890
  * [5] Alessio Zappone and Eduard Jorswieck  “Energy efficiency in wireless networks via fractional programming theory”  In _Foundations and Trends in Communications and Information Theory_ 11.3-4 Now Publishers Inc. Hanover, MA, USA, 2015, pp. 185–396
  * [6] Christian Isheden, Zhijiat Chong, Eduard Jorswieck and Gerhard Fettweis  “Framework for link-level energy efficiency optimization with informed transmitter”  In _IEEE Transactions on Wireless Communications_ 11.8 IEEE, 2012, pp. 2946–2957
  * [7] Derrick Wing Kwan Ng, Ernest S Lo and Robert Schober  “Energy-efficient resource allocation in multi-cell OFDMA systems with limited backhaul capacity”  In _IEEE Transactions on Wireless Communications_ 11.10 IEEE, 2012, pp. 3618–3631
  * [8] Jie Xu and Ling Qiu  “Energy efficiency optimization for MIMO broadcast channels”  In _IEEE Transactions on Wireless Communications_ 12.2 IEEE, 2013, pp. 690–701
  * [9] Shiwen He, Yongming Huang, Shi Jin and Luxi Yang  “Coordinated beamforming for energy efficient transmission in multicell multiuser systems”  In _IEEE Transactions on Communications_ 61.12 IEEE, 2013, pp. 4961–4971
  * [10] Shiwen He, Yongming Huang, Luxi Yang and Björn Ottersten  “Coordinated multicell multiuser precoding for maximizing weighted sum energy efficiency”  In _IEEE Transactions on Signal Processing_ 62.3 IEEE, 2013, pp. 741–751
  * [11] Bo Du, Cunhua Pan, Wence Zhang and Ming Chen  “Distributed energy-efficient power optimization for CoMP systems with max-min fairness”  In _IEEE Communications Letters_ 18.6 IEEE, 2014, pp. 999–1002
  * [12] Quansheng Xu, Xi Li, Hong Ji and Xiaojiang Du  “Energy-efficient resource allocation for heterogeneous services in OFDMA downlink networks: Systematic perspective”  In _IEEE Transactions on Vehicular Technology_ 63.5 IEEE, 2014, pp. 2071–2082
  * [13] Alessio Zappone et al.  “Energy-efficient power control: A look at 5G wireless technologies”  In _IEEE Transactions on Signal Processing_ 64.7 IEEE, 2015, pp. 1668–1683
  * [14] Bho Matthiesen and Eduard Axel Jorswieck  “Efficient global optimal resource allocation in non-orthogonal interference networks”  In _IEEE Transactions on Signal Processing_ 67.21 IEEE, 2019, pp. 5612–5627
  * [15] Karl-Ludwig Besser, Bho Matthiesen, Alessio Zappone and Eduard A. Jorswieck  “Deep Learning Based Resource Allocation: How Much Training Data is Needed?”  In _2020 IEEE 21st International Workshop on Signal Processing Advances in Wireless Communications (SPAWC)_ IEEE, 2020  DOI: 10.1109/SPAWC48557.2020.9154298
  * [16] Pieter-Tjerk De Boer, Dirk P Kroese, Shie Mannor and Reuven Y Rubinstein  “A tutorial on the cross-entropy method”  In _Annals of operations research_ 134.1 Springer, 2005, pp. 19–67
  * [17] Benjamin Recht  “A tour of reinforcement learning: The view from continuous control”  In _Annual Review of Control, Robotics, and Autonomous Systems_ 2 Annual Reviews, 2019, pp. 253–279
  * [18] Anusha Nagabandi, Gregory Kahn, Ronald S Fearing and Sergey Levine  “Neural network dynamics for model-based deep reinforcement learning with model-free fine-tuning”  In _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , 2018, pp. 7559–7566  IEEE
  * [19] Dhananjay Kumar, NN Kanagaraj and R Srilakshmi  “Harmonized Q-learning for radio resource management in LTE based networks”  In _2013 Proceedings of ITU Kaleidoscope: Building Sustainable Communities_ , 2013, pp. 1–8  IEEE
  * [20] Kan Zheng et al.  “Big data-driven optimization for mobile networks toward 5G”  In _IEEE network_ 30.1 IEEE, 2016, pp. 44–51
  * [21] Amna Mughees, Mohammad Tahir, Muhammad Aman Sheikh and Abdul Ahad  “Towards energy efficient 5G networks using machine learning: Taxonomy, research challenges, and future research directions”  In _IEEE Access_ 8 IEEE, 2020, pp. 187498–187522
  * [22] Ramprasad Raghunath, Bile Peng, Karl-Ludwig Besser and Eduard A. Jorswieck  “Reinforcement Learning-Based Global Programming for Energy Efficiency in Multi-Cell Interference Networks”  In _ICC 2022 – IEEE International Conference on Communications_ IEEE, 2022, pp. 5499–5504  DOI: 10.1109/ICC45855.2022.9838408
  * [23] Prateek Jain and Purushottam Kar  “Non-convex optimization for machine learning”, 2017  arXiv:1712.07897 [stat.ML]
  * [24] Ming Xu, Matias Quiroz, Robert Kohn and Scott A Sisson  “Variance reduction properties of the reparameterization trick”  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , 2019, pp. 2711–2720  PMLR
  * [25] James T Wilson, Riccardo Moriconi, Frank Hutter and Marc Peter Deisenroth  “The reparameterization trick for acquisition functions”, 2017  arXiv:1712.00424 [stat.ML]
  * [26] Anthony L Caterini, Arnaud Doucet and Dino Sejdinovic  “Hamiltonian variational auto-encoder”  In _Advances in Neural Information Processing Systems_ 31, 2018
  * [27] Tuomas Haarnoja et al.  “Soft actor-critic algorithms and applications”, 2018  arXiv:1812.05905 [cs.LG]
  * [28] Jonathan Gordon, David Lopez-Paz, Marco Baroni and Diane Bouchacourt  “Permutation equivariant models for compositional generalization in language”  In _International Conference on Learning Representations_ , 2019
  * [29] Siamak Ravanbakhsh, Jeff Schneider and Barnabas Poczos  “Equivariance through parameter-sharing”  In _International Conference on Machine Learning_ , 2017, pp. 2892–2901  PMLR
  * [30] Kumar Pratik, Bhaskar D Rao and Max Welling  “RE-MIMO: Recurrent and permutation equivariant neural MIMO detection”  In _IEEE Transactions on Signal Processing_ 69 IEEE, 2020, pp. 459–473
  * [31] Mark Eisen and Alejandro Ribeiro  “Optimal wireless resource allocation with random edge graph neural networks”  In _IEEE Transactions on Signal Processing_ 68 IEEE, 2020, pp. 2977–2991
  * [32] Zhiyang Wang, Mark Eisen and Alejandro Ribeiro  “Learning decentralized wireless resource allocations with graph neural networks”  In _IEEE Transactions on Signal Processing_ 70 IEEE, 2022, pp. 1850–1863
  * [33] Arindam Chowdhury et al.  “Unfolding WMMSE using graph neural networks for efficient power allocation”  In _IEEE Transactions on Wireless Communications_ 20.9 IEEE, 2021, pp. 6004–6017
  * [34] Günter Rudolph  “Globale Optimierung mit parallelen Evolutionsstrategien”, 1990
  * [35] Matthew J Johnson et al.  “Composing graphical models with neural networks for structured representations and fast inference”  In _Advances in neural information processing systems_ 29, 2016
  * [36] Matthias Fey and Jan Eric Lenssen  “Fast graph representation learning with PyTorch Geometric”, 2019  arXiv:1903.02428 [cs.LG]
  * [37] George Calcev et al.  “A wideband spatial channel model for system-wide simulations”  In _IEEE Transactions on Vehicular Technology_ 56.2 IEEE, 2007, pp. 389–403
  * [38] Adam Paszke et al.  “PyTorch: An Imperative Style, High-Performance Deep Learning Library”  In _Advances in Neural Information Processing Systems 32_ Curran Associates, Inc., 2019, pp. 8024–8035  URL: http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf
  * [39] Diederik P Kingma and Jimmy Ba  “Adam: A method for stochastic optimization”, 2014  arXiv:1412.6980 [cs.LG]