# DEBOSH: Deep Bayesian Shape Optimization

Nikita Durasov,1 Artem Lukoyanov,2 Jonathan Donier,2 Pascal Fua1

###### Abstract

Shape optimization is at the heart of many industrial applications, such as
aerodynamics, heat transfer, and structural analysis. It has recently been
shown that Graph Neural Networks (GNNs) can predict the performance of a shape
quickly and accurately and be used to optimize more effectively than
traditional techniques that rely on response-surfaces obtained by Kriging.
However, GNNs suffer from the fact that they do not evaluate their own
accuracy, which is something Bayesian Optimization methods require. Therefore,
estimating confidence in generated predictions is necessary to go beyond
straight deterministic optimization, which is less effective.

In this paper, we demonstrate that we can use Ensembles-based technique to
overcome this limitation and outperform the state-of-the-art. Our experiments
on diverse aerodynamics and structural analysis tasks prove that adding
uncertainty to shape optimization significantly improves the quality of
resulting shapes and reduces the time required for the optimization.

##  1 Introduction

Simulation is a powerful tool to refine the 3D shape of a design. For example,
Computational Fluid Dynamics (CFD) is central to maximizing the performance of
aircraft wings, windmill blades, hydrofoils, car bodies, ship hulls, and
propellers, while Structural Analysis is key to designing structures that are
both light and strong. In both cases, because the simulation must be re-run
each time the design is changed, a widespread engineering practice is to test
only a few designs, without searching the design space. The same can be said
in many other fields, such as heat transfer management or electronic component
design, that involve complicated simulations to optimize shapes.

Figure 1: Shape Optimization. Starting from 3D shapes in initial training
dataset, DEBOSH uses a GCNN to compute a response surface and Bayesian
Optimization to efficiently produce an optimized shape that it has never seen
before.

One way to explore the shape space more thoroughly is to use genetic
algorithms (Gosselin, Tye-Gingras, and Mathieu-Potvin 2009). However, in their
simplest form, they require many evaluations of a fitness function and
therefore many expensive simulations, which makes them inefficient. Relying on
Kriging, also known as Gaussian-process (GP) regression (Rasmussen and
Williams 2006), is one of the most popular ways to reduce the required number
of evaluations (Jeong, Murayama, and Yamamoto 2005; Laurenceau et al. 2010;
Toal and Keane 2011; Xu et al. 2017). Unfortunately, Kriging works best for
models that can be parameterized using relatively few parameters, which limits
its applicability.

Other methods rely on topology optimization (Saviers, Ranjan, and Mahmoudi
2019) or adjoint differentiation (Allaire 2015; Gao, Wu, and Xia 2017; Behrou,
Ranjan, and Guest 2019). The first is highly effective but only applicable in
very specific cases, while the second approximates the solution of a so-called
adjoint form of the partial differential equations to compute the gradients of
the fitness function with respect to the 3D simulation mesh. However, this
only applies to relatively small deformations with respect to an original
shape and requires careful re-meshing of the shape and re-running the
simulation once the changes exceed a threshold. Furthermore, this method is
even more computationally expensive than those discussed above (Alexandersen,
Sigmund, and Aage 2016).

A radically different approach was proposed in (Baqué et al. 2018): Given a
collection of 3D surface meshes and corresponding results of simulations, it
is possible to train a Geodesic Convolutional Neural Networks (GCNNs)
(Boscaini et al. 2016; Monti et al. 2017) to emulate a complex fluid-dynamics
simulator. This makes it possible to write a differentiable objective function
that returns aerodynamic and hydrodynamic properties—drag, lift,
stability—very fast and therefore to maximize performance with respect to the
3D shape of the target objects, no matter how many parameters are required to
define its shape. That gives it an edge over Kriging that can only deal with
relatively fewer parameters. However, unlike Kriging, it suffers from the fact
that it does not account for the potential imprecision of its prediction.

To remedy this, we introduce a new Bayesian shape optimization algorithm,
which we dub DEBOSH. We start with the GCNN of (Baqué et al. 2018) and give it
the ability to estimate its own uncertainly using an ensemble-based approach
similar to (Lakshminarayanan, Pritzel, and Blundell 2017; Gal and Ghahramani
2016a; Durasov et al. 2021; Wen, Tran, and Ba 2020). This enables us to rely
on Bayesian Optimization (BO) to explore the shape space and re-train the
network when necessary to produce shapes such as the one depicted by Fig. 1.
As a result DEBOSH outperforms both earlier deep learning and classical
approaches in terms of quality of the optimized shapes and required time to
obtain them, which we demonstrate on several benchmarks.

##  2 Related Work

In all engineering fields that rely on running computationally demanding
simulations on 3D shapes to estimate their performance, optimizing these
shapes to maximize this performance is a difficult problem because, in
practice, the simulator can only be run so many times. In this section, we
briefly review some of the dominant approaches to addressing this problem. We
also introduce standard Bayesian Optimization, which we will rely on in the
method section.

###  2.1 Shape Optimization

A popular way to automate the shape optimization processes is to use genetic
algorithms (Gosselin, Tye-Gingras, and Mathieu-Potvin 2009) to explore the
shape space. However, as genetic algorithms require many evaluations of the
fitness function, which involves running an expensive physical simulation, a
naive implementation would be inefficient. Kriging, or Gaussian-process (GP)
regression (Rasmussen and Williams 2006), is one of the most popular ways to
reduce the required number of evaluations (Jeong, Murayama, and Yamamoto 2005;
Laurenceau et al. 2010; Toal and Keane 2011; Xu et al. 2017; Umetani and
Bickel 2018). It estimates a response surface that is used to accelerate the
exploration of the shape space. Although effective when dealing with models
controlled by relatively few design parameters, it does not scale well to
high-dimensional inputs or to sparse GPs (Herbrich, Lawrence, and Seeger
2003), due to the curse of dimensionality (Bellman 1961). As a result, these
methods often depend on handcrafted low-dimensional shape parameterizations
that can be difficult to design and severely restrict how much of the shape
space can be explored.

To alleviate this problem, it has been proposed to apply deep learning methods
(Vaswani et al. 2017; Chen et al. 2017; Redmon et al. 2016; Durasov et al.
2019) and Convolutional Neural Networks (CNNs) in 3D voxelized grids. This has
been used to accelerate the search for solutions of the discrete Poisson
equation (Tompson et al. 2017), to directly regress the fluid velocity fields
given an implicit surface description (Guo, Li, and Iorio 2016), or to compute
fluid simulations velocities from a set of reduced parameters (Kim et al.
2019). However, the underlying 3D CNN architectures working with spatial
voxels tend to feature large memory footprints and computationally demanding
inference. To work around this problem, it has been proposed to replace the
GPs with Geodesic Convolutional Neural Networks (GCNNs) (Boscaini et al. 2016;
Monti et al. 2017) to emulate a simulator (Baqué et al. 2018; Remelli et al.
2020). Given a set of generic surfaces parametrized as meshes, a GCNN is
trained to predict their aerodynamic characteristics, as computed by standard
CFD packages (Drela 1989; Weller et al. 1998; Mountrakis et al. 2015). These
characteristics are then used to write an objective function that is
differentiable with respect to the shape parameters. To maximize performance,
the objective function is minimized with respect to these parameters. This
makes it possible to deal with more complex parameterization but fails to
estimate uncertainties, which hampers the use of so-called exploration-and-
exploitation algorithms such as those discussed below, which are very powerful
for these kinds of applications (Laurenceau et al. 2010).

###  2.2 Bayesian Optimization

Bayesian Optimization (BO) (Mockus 2012) is an exploration-and-exploitation
approach to finding global minima of a multivariate function
$g:\mathbf{{A}}\rightarrow\mathbb{R}$ that does not assume any specific
functional form for $g$ and is usually used for functions that are expensive
to evaluate. In our application $g$ is implemented by a physics-based
simulator.

BO typically starts with a surrogate model
$\widetilde{g}_{\Theta}:\textbf{A}\rightarrow\mathbb{R}$ whose output depends
on a set of parameters $\Theta$. $\widetilde{g}_{\Theta}$ is assumed to
approximate $g$, to be fast to compute, and to be able to evaluate the
reliability of its own predictions in terms of a variance. It is used to
explore $\mathbf{A}$ quickly in search of a solution of
$\textbf{z}^{*}=\operatornamewithlimits{argmin}_{\textbf{z}\in\mathbf{{A}}}g(\mathbf{z})$.
Given an initial training set $\\{(\textbf{z}_{i},r_{i})\\}_{i}$ of input
vectors $\textbf{z}_{i}$ and outputs $r_{i}=g(\textbf{z}_{i})$, it iterates
the following steps:

  1. 1.

Find $\Theta$ that yields the best possible prediction by
$\widetilde{g}_{\Theta}$.

  2. 2.

Generate new samples not present in the training set.

  3. 3.

Evaluate an acquisition function on these samples.

  4. 4.

Add the best ones to the training set and go to step 1.

An effective approach to generating the new samples is to apply a genetic
algorithm (Gosselin, Tye-Gingras, and Mathieu-Potvin 2009) to the surrogate
model $\widetilde{g}_{\Theta}$, which is much faster than applying it to the
original function $g$. For each new sample z, one must compute the acquisition
function. It is often taken to be the Expected Improvement (EI), which, in
essence, tries to find the point that minimizes the distance to the objective
evaluated at the maximum and under Gaussian assumptions can be written as

$\displaystyle EI(\textbf{z})$
$\displaystyle=\left(\widetilde{r}-g(\textbf{z}^{+})-\epsilon\right)\Phi(Z)+\widetilde{\sigma}\phi(Z)\;,$
(1) $\displaystyle Z$
$\displaystyle=\frac{\widetilde{r}-g(\textbf{z}^{+})-\epsilon}{\widetilde{\sigma}}\;,$

where $\widetilde{r}$ is the value predicted by $\widetilde{g}_{\Theta}$,
$\widetilde{\sigma}$ is the corresponding standard deviation, $\textbf{z}^{+}$
is the best exemplar found so far, $\Phi(\cdot)$ is cumulative distribution
function (CDF) and $\phi(\cdot)$ is probability density function (pdf) of
standard Gaussian distribution, $\epsilon\in\mathbb{R}$ is a hyperparameter.

##  3 Method

In this section, we introduce our DEBOSH approach to shape optimization. As
depicted by Fig. 2, the algorithm iterates the following steps:

  1. 1.

Given a set of shapes, run physical simulations on each one.

  2. 2.

Use the training shapes and the simulation results to train one GCNN, or more,
to emulate the physical simulations.

  3. 3.

Generate new shapes and use GCNN(s) to evaluate their acquisition function, as
described in 2.2.

  4. 4.

Add the best new shapes to the set of shapes and iterate.

These four steps are performed either a preset number of times or until
reaching given optimization objectives. The main obstacle to implementing this
scheme is that GCNNs do not naturally return the uncertainties that evaluation
of the acquisition function requires. We overcome it using either MC-Dropout
(Gal and Ghahramani 2016b) or Deep Ensembles (Lakshminarayanan, Pritzel, and
Blundell 2017), as described below.

Figure 2: DEBOSH pipeline: (1) Run physical simulations. (2) Train the GCNN.
(3) Perform Bayesian Optimization. (4) Select promising samples, add them to
the training set, and go back to step 1.

###  3.1 Training Set

Given a set of $N$ 3D shapes $\\{\textbf{x}_{i}\\}_{1\leq i\leq N}$
represented by triangulated meshes, we run a physics-based simulator numerical
solver to create a corresponding set $\\{\textbf{y}_{i}\\}_{1\leq i\leq N}$ of
physical values, such as pressure or mechanical stress at each vertex. Given a
field y, we derive a performance value $r=\textbf{R}(\textbf{y})$, such as
overall drag for a car or maximum stress for a mechanical part, using a
function R that is task-specific. Assuming that each mesh x is parameterized
by a low-dimensional latent vector z and that there is a differentiable
mapping $\textbf{P}:\textbf{z}\rightarrow\textbf{x}$, this gives us the
initial training set $\\{(\textbf{z}_{i},r_{i})\\}_{i}$ that we need to
initialize Bayesian Optimization, as described in Section 2.2.

The $\textbf{x}_{i}$ are represented by triangulated meshes with vertices and
edges $(V_{i},E_{i})$, along with node and edge attributes
$\\{\textbf{u}_{i}^{j}\in\mathbb{R}^{M}\\}_{j=1}^{|V_{i}|}$ and
$\\{\textbf{v}_{i}^{j}\in\mathbb{R}^{L}\\}_{j=1}^{|E_{i}|}$, where $M$ and $L$
are the dimensions of these attributes. In practice, we take
$\textbf{u}_{i}^{j}$ to be the concatenation of vertex’s xyz-coordinates with
the result of the sinusoidal transformation
$\\{(\sin{ax},\sin{ay},\sin{az})\\}_{a=1}^{A}$ of these coordinates to encode
high frequencies. Edge features $\textbf{v}_{i}^{j}$ are the difference of
xyz-coordinates of nodes belonging to a given edge.

###  3.2 Emulator Architecture

Given the training set of Section 3.1, we train a GCNN (Monti et al. 2017) to
predict the field of physical values y given a mesh x, as depicted by Fig. 3.
The forward pass through a single layer can be written as

$\widetilde{\textbf{u}}_{i}^{j}=\frac{1}{|\mathcal{N}(j)|}\sum_{n\in\mathcal{N}(j)}\frac{1}{K}\sum_{k=1}^{K}\textbf{w}_{k}(\textbf{v}_{i}^{j,n})\otimes\Theta_{k}\textbf{u}_{i}^{j}$
(2)

where $\mathcal{N}(j)$ is the neighborhood of the $j$-th node and

$\textbf{w}_{k}(\textbf{e})=\exp\left(-\frac{1}{2}(\textbf{e}-\mu_{k})^{T}\Sigma_{k}^{-1}(\textbf{e}-\mu_{k})\right)$
(3)

denotes a weighting function based on trainable mean vector $\mu_{k}$ and
diagonal covariance matrix $\Sigma_{k}^{-1}$. Such a layer involves three
types of trainable weights: $\Theta_{k}\in\mathbb{R}^{\widetilde{L}\times L}$
maps input $\textbf{u}_{i}^{j}$ to the new dimensionality, and $\mu_{k}$ and
$\Sigma_{k}$ define the Gaussian kernel used to convolve a node with its
neighbors. Our final model consists of 25 GNN layers with ELU nonlinearities
(Clevert, Unterthiner, and Hochreiter 2015) and skip connections (He et al.
2016).

Figure 3: Emulation model for pressure prediction. Given an input geometry
mesh our GCNN outputs a field of pressure values on the sample’s surface
without having to run a computationally expensive simulation.

###  3.3 Bayesian Optimization

Recall from Section 2.2 that Bayesian Optimization (BO) is an exploration-and-
exploitation approach to minimizing a function $g$ given a trainable estimator
$\widetilde{g}_{\Theta}$ that can evaluate its own uncertainty. In our case,
$g$ is the function that takes as input a latent vector z, produces a
triangulated mesh $\textbf{x}=\textbf{P}(\textbf{z})$, computes the associated
field of physical values y using a simulator, and finally outputs a
performance value $r=R(\textbf{y})$.

We now need to define $\widetilde{g}_{\Theta}$. Given a latent vector z and
the GCNN of Section 3.2 with weights $\Theta$, we can approximate the
performance by feeding the triangulated mesh $\textbf{P}(\textbf{z})$ to the
network to produce an estimated field of values $\tilde{\textbf{y}}$ and an
estimated performance value $\tilde{r}=R(\tilde{\textbf{y}})$. We take this to
be the estimate that $\widetilde{g}_{\Theta}$ but we also need an uncertainty
estimate in the form of a variance, which means assessing the reliability of
our network. Among all the methods that can be used to assess to this end, MC-
Dropout (Gal and Ghahramani 2016b) and Deep Ensembles (Lakshminarayanan,
Pritzel, and Blundell 2017) have emerged as two of the most popular ones. MC-
Dropout involves randomly zeroing out network weights and assessing the
effect, whereas Ensembles involve training multiple networks starting from
different initial conditions. In practice, the latter tends to perform better
but also to be more computationally demanding due to the requirement of
restarting the training procedure from scratch many times.

In practice, we implemented both approaches. For ensembles approaches, for
each training set, we train the $N$ different networks with weights
$\Theta_{1},\ldots\Theta_{n}$ and use the mean and variance of the estimates

$\displaystyle\mu(\textbf{z})$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\widetilde{\textbf{g}}_{\Theta_{i}}(\textbf{x})\;,$
(4) $\displaystyle\sigma^{2}(\textbf{z})$
$\displaystyle=\frac{1}{N-1}{\sum_{i=1}^{N}(\widetilde{\textbf{g}}_{\Theta_{i}}(\textbf{x})-\mu(x))^{2}}\;.$

For MC-Dropout, we proceed similarly but, instead of using different weights
$\Theta_{i}$, we run the inference several times to get multiple estimates,
along with their mean and variance.

##  4 Experiments

In this section, we compare DEBOSH against state-of-the-art alternatives in
several scenarios.

###  4.1 Baselines and Variants

We will compare four different approaches to performing shape optimization,
which we describe briefly below.

  * •

Kriging: Using a GP to estimate a response surface and corresponding
uncertainty from a number of shape exemplars and corresponding simulations
results is known as Kriging (Laurenceau et al. 2010). This response surface
can then be used to optimize shapes using Bayesian Optimization. As discussed
in Section 2.1, it is a widely used practice and is a standard approach to BO.

  * •

GCNN: An alternative to using GPs to estimate the response surface is to use a
standard GCNN (Baqué et al. 2018; Remelli et al. 2020), as discussed in
Section 2.1. Since a GCNN does not return uncertainties, the optimization is
performed using deterministic gradient-descent.

  * •

Ours-MCD: Our method when using MC-Dropout to estimate uncertainties.

  * •

Ours-Ens: Our method when using deep ensembles to estimate uncertainties.

We report the results obtained by all four methods on three different types of
optimization problems in Table. 1. Ours-Ens clearly outperforms the other
approaches with Ours-MCD being second best. In other words, using a GCNN to
estimate both the fitness of a shape and the corresponding uncertainty is
beneficial and using ensembles is even more so than using drop-out. We discuss
how we obtained these results in more details below.

| Kriging | GCNN | Ours-MCD | Ours-Ens
---|---|---|---|---
Foils: Lift-to-Drag | $1.62$ | $1.81$ | 1.94 | 2.08
Cars: Drag | 3.4 | 2.8 | 2.6 | 2.3
Knuckles: Stress | 105. | 100. | 95. | 91.
Table 1: Optimized fitness values for airfoils, cars, and knuckles. Highest
lift-to-drag and lowest drag and stress values for the three types of shapes
we tested all four methods on. Ours-Ens systematically delivers the best
values with Ours-MCD being second best.

###  4.2 2D Airfoil Optimization

Figure 4: Lift-to-drag maximization for airfoils. We ran each training
procedure 5 times with different random seeds. Each line represents the
average negative lift-to-drag-ratio of the best performing profile as a
function of the number of training iterations. The shaded areas depict the
corresponding variances. They are color coded according to the method used to
perform the minimization, Ours-Ens, Ours-MCD, Kriging, and GCNN. The Ours-Ens
curve is below the others.

2D airfoil profile optimization has become a de facto standard for
benchmarking shape optimization in the CFD community. In industrial practice,
profiles have long been parameterized using a three-dimensional NACA parameter
vector (Jacobs, Ward, and Pinkerton 1948) and optimized using conventional
Kriging-based methods (Jeong, Murayama, and Yamamoto 2005; Chiplunkar, Bosco,
and Morlier 2017). We generated an initial dataset of 200 shapes for airfoil
optimization by randomly selecting NACA parameters $\textbf{z}_{i}$ and then
producing corresponding 2D contours $\textbf{x}_{i}$ for each one.

We use the popular XFoil simulator (Drela 1989) to compute the pressure
distribution $\textbf{y}_{i}$ over the surface of each $\textbf{x}_{i}$. Even
though GCNNs models are primarily designed to handle 3D shapes, they can also
handle 2D ones by considering the 2D equivalent of a surface mesh, which is a
discretized 2D contour. As in (Baqué et al. 2018), train one to predict
pressure values $\tilde{\textbf{y}_{i}}$ at each vertex along with the overall
lift-to-drag ratio, which we take to be $R(\textbf{y}_{i})$, which could also
be obtained by integrating the pressure values over the surface of the shape.

Figure 5: Optimized 2D wing profiles. Best performing profile found by each
one of the four methods we tested.

As discussed in Section 3.3, for Ours-Ens we train five GCNNs at each
iteration while we train only one Ours-MCD. In both cases, after each training
iteration, we use the genetic algorithm NSGAII (Deb et al. 2002) coupled to
the GCNN(s) to sample 30 new airfoil shapes and to evaluate the corresponding
acquisition function. For Kriging, we proceed similarly but use GPs instead of
a GCNN to run the generic algorithm and evaluate the acquisition function. For
GCNN, we use the network to predict lift-to-drag values
$R(\widetilde{\textbf{y}}_{i})$ and pick sampled shapes for simulations
according to this metric instead of the acquisition function, which would
require an uncertainty estimate.

Fig. 4 depicts the lift-to-drag-ratios for all four approaches as function of
the number of training iterations. Fig. 5 depicts the best foils found by each
one of the four approaches. We report the corresponding lift-to-drag ratios in
Table 1.

###  4.3 Minimizing Car Drag

Figure 6: Drag minimization for cars. As in Fig. 4, the individual curves are
color coded. They represent drag values and their variances as a function of
the number of training iterations for Ours-Ens and Ours-MCD, Kriging, and
GCNN. The Ours-Ens curve is again below the others.

Even though airfoil optimization using the NACA parameterization is a standard
benchmark, it relies on a 3-Dimensional latent vector and is therefore very
low-dimensional. We now turn to a higher-dimensional problem in which the
latent vector is of dimension 256: car-shape optimization to minimize drag.

As in (Remelli et al. 2020), we use a cleaned-up and processed subset of the
ShapeNet dataset (Chang et al. 2015) that features $N=1400$ car meshes that
are suitable for CFD simulation. For each such mesh $\textbf{x}_{i}$, we run
OpenFOAM (Jasak et al. 2007) to estimate the pressure field $\textbf{y}_{i}$
created by air travelling at 15 meters per second towards the car. We also
train the MeshSDF (Remelli et al. 2020) mapping P on this training dataset so
that for each $\textbf{x}_{i}$ there is 256-dimensional latent vector
$\textbf{z}_{i}$ such that $\textbf{x}_{i}=\textbf{P}(\textbf{z}_{i})$. Though
many mesh and graph embedding methods were proposed recently in literature
(Park et al. 2019; Mescheder et al. 2019; Ivanov, Durasov, and Burnaev 2018;
Baqué et al. 2018), MeshSDF proved to be one of the most efficient and
expressive in CFD context.

Figure 7: Gradient descent-based drag minimization. Starting from the initial
shape on the left, $R(\widetilde{\textbf{z}})$ from Eq. 6 is minimized with
respect to z using gradient descent steps. The middle shapes depict the
optimization process and the relative improvements in terms of defined
$R(\widetilde{\textbf{z}})$. The final result is shown on the right. The
mapping P of (Remelli et al. 2020) enforces implicit semantic priors while
being able to effect topology changes that result in plausible cars shapes.

Given such a training set, we can train a GCNN such as the one of Section 3.2
to predict a pressure field $\widetilde{\textbf{y}}$ that approximates the
true one for a car shape $\textbf{x}=\textbf{P}(\textbf{z})$. Its drag can be
approximated by the surface integral

$D(\widetilde{\textbf{y}})=\iint_{S}\widetilde{y}\,\textbf{n}_{x}\,dS\;,$ (5)

where $\textbf{n}_{x}$ denotes the projection of the surface normals along the
airflow direction and $\widetilde{y}$ is the scalar value of the pressure
field at a given surface point.

Figure 8: Optimized 3D car shapes. Best car shapes found by each one of the
four methods we tested. (A) Kriging (B) GCNN (C) Ours-MCD (D) Ours-Ens. Note
the artifact produced by Kriging that the other methods do not produce. Note
also, that the topology of the Ours-MCD and Ours-Ens cars are different, which
is made possible by the fact that the latent representation we use allows
topology changes.

We could take $D$ to be the function we minimize but all the shapes would
reduce to a point. Instead, we write

$R(\widetilde{\textbf{z}})=D(\widetilde{\textbf{y}})+R_{aux}(\textbf{z})\;,$
(6)

where $R_{aux}$ is a penalty term that becomes large if the car shape strays
away too much from those in the training set. In practice, we take it to be
the average Euclidean distance to the five closest neighbors of z in the
training set. During minimization, keeping $R_{aux}$ forces z to remain close
from the parts of the shape space for which the parameterizer $\mathbf{P}$ has
been trained and to keep on generating plausible cars. Adding this penalty
term constrains the exploration of the latent space but guarantees more robust
and realistic optimization outcomes.

When performing BO to minimize $R$, we could use a genetic algorithm to create
new samples as we did for the air foils. However, we now have to explore a
256-dimensional shape space, and each evaluation of the surrogate model, while
much faster than an actual simulation, requires computing P. This involves
running the Marching Cube algorithms to compute the mesh vertices (Lorensen
and Cline 1987), which takes time and makes vanilla genetic algorithms very
slow.

Instead, we take advantage of the fact that, even though the Marching Cube
algorithms is non-differentiable, the mapping P is differentiable (Guillard et
al. 2021) as follows. At each training iteration, we start from the existing
elements of the training set and perform gradient descent using them as a
starting point. Fig. 7 depicts this process. This gives a new set of shapes
and we select the best ones according to the acquisition function, as in
standard BO.

Fig. 6 depicts the resulting drag coefficients for all four approaches as
function of the number of training iterations. Again, Ours-Ens outperforms the
other approaches, followed by Ours-MCD. We report the final numbers in Table 1
and the corresponding car shapes are shown in Fig. 8.

###  4.4 Mechanical Stress Minimization

Figure 9: Stress minimization problem. The plot represents a setup for stress
prediction and minimization for special detail called knuckle. Horizontal
force $\vec{\textbf{F}}$ applied to the tip of the handle induces tension
distributed throughout the entire volume of the detail. Proper design of the
joint is required to produce durable and long-lasting construction.

To show that we can handle a very different kind of problem, we now consider a
structural engineering problem. Our goal is to minimize the mechanical
stresses in the knuckle depicted by Fig. 9. It is made of a metal ring with a
long handle attached to it. Fixing the inner circle of the ring and applying a
horizontal force to the handle tip creates mechanical stress. Minimizing it is
desirable because it significantly extends the lifespan of the part and can be
achieved by shaping properly the most fragile part, which is the joint between
the handle and the ring.

Figure 10: Optimized 3D knuckle. Starting from 3D shapes such as the one at
the top, Ours-Ens searches for optimized RBF parameters to strengthen the
joint between the handle and the ring, the most fragile part of the assembly.
The final result is shown at the bottom. As expected, the method makes the
joint heavier and therefore better able to handle the external force
$\vec{\textbf{F}}$.

To test this, we have generated $N=200$ knuckle shapes with different joint
configurations. Variation in joint shapes is parametrized using the Radial
Basis Function (RBF) parametrizer (Forti and Rozza 2014) P that relies on
control points and RBF interpolation. For each shape $\textbf{x}_{i}$, we use
a finite element solver (Cimrman 2021) that computes the von Mises stresses
$\textbf{y}_{i}$ across the joint and we write

$R({\textbf{y}}_{i})=\max_{i}{\textbf{y}}_{i}\;.$ (7)

Since the latent vector is of dimension 9, we use the same approach to BO as
for the airfoil. At every training iteration, we sample 30 new shapes and use
the acquisition function to retain 10 shapes and add them to original training
dataset. Fig. 11 depicts the resulting stresses for all four approaches as
function of the number of training iterations. Again, Ours-Ens outperforms the
other approaches. We report the final numbers in Table 1 and show the best
final shape in Fig. 10.

Figure 11: Stress minimization for knuckles. As in Fig. 4, the individual
curves are color coded. They represent stress values and their variances as a
function of the number of training iterations for Ours-Ens and Ours-MCD,
Kriging, and GCNN. The Ours-Ens curve is again below the others.

##  5 Conclusion

In this work, we have shown that neural networks could be used as surrogate
models for Bayesian Optimization. We have demonstrated on several realistic
scenarios that using either MC-dropout or deep ensembles is an effective way
to estimate the uncertainty of the network’s prediction, which is a key
ingredient of BO. MC-Dropout tends to be faster but Ensembles deliver better
results. In future work, we will focus on developing approaches to estimating
uncertainty that combine the accuracy of Ensembles with the speed of MC-
dropout.

## References

  * Alexandersen, Sigmund, and Aage (2016) Alexandersen, J.; Sigmund, O.; and Aage, N. 2016.  Large Scale Three-Dimensional Topology Optimization of Heat Sinks Cooled by Natural Convection.  _International Journal of Heat and Mass Transfer_ , 100(Supplement C): 876–891.
  * Allaire (2015) Allaire, G. 2015.  A Review of Adjoint Methods for Sensitivity Analysis, Uncertainty Quantification and Optimization in Numerical Codes.  _Ingénieurs de l’Automobile_ , 836: 33–36.
  * Baqué et al. (2018) Baqué, P.; Remelli, E.; Fleuret, F.; and Fua, P. 2018.  Geodesic Convolutional Shape Optimization.  In _International Conference on Machine Learning_.
  * Behrou, Ranjan, and Guest (2019) Behrou, R.; Ranjan, R.; and Guest, J. K. 2019.  Adaptive Topology Optimization for Incompressible Laminar Flow Problems with Mass Flow Constraints.  _Computer Methods in Applied Mechanics and Engineering_ , 346: 612–641.
  * Bellman (1961) Bellman, R. 1961.  _Curse of Dimensionality_.  Princeton university press.
  * Boscaini et al. (2016) Boscaini, D.; Masci, J.; Rodolà, E.; and Bronstein, M. 2016.  Learning Shape Correspondence with Anisotropic Convolutional Neural Networks.  In _Advances in Neural Information Processing Systems_ , 3189–3197.
  * Chang et al. (2015) Chang, A.; Funkhouser, T.; G., L.; Hanrahan, P.; Huang, Q.; Li, Z.; Savarese, S.; Savva, M.; Song, S.; Su, H.; Xiao, J.; Yi, L.; and Yu, F. 2015.  Shapenet: An Information-Rich 3D Model Repository.  In _arXiv Preprint_.
  * Chen et al. (2017) Chen, L.-C.; Papandreou, G.; Kokkinos, I.; Murphy, K.; and Yuille, A. L. 2017.  Deeplab: Semantic image segmentation with deep convolutional nets, atrous convolution, and fully connected crfs.  _IEEE transactions on pattern analysis and machine intelligence_ , 40(4): 834–848.
  * Chiplunkar, Bosco, and Morlier (2017) Chiplunkar, A.; Bosco, E.; and Morlier, J. 2017.  Gaussian Process for Aerodynamic Pressures Prediction in Fast Fluid Structure Interaction Simulations.  In _World Congress of Structural and Multidisciplinary Optimisation_ , 221–233.
  * Cimrman (2021) Cimrman, R. 2021.  Fast Evaluation of Finite Element Weak Forms Using Python Tensor Contraction Packages.  _Advances in Engineering Software_ , 159.
  * Clevert, Unterthiner, and Hochreiter (2015) Clevert, D.-A.; Unterthiner, T.; and Hochreiter, S. 2015.  Fast and accurate deep network learning by exponential linear units (elus).  _arXiv preprint arXiv:1511.07289_.
  * Deb et al. (2002) Deb, K.; Pratap, A.; Agarwal, S.; and Meyarivan, T. 2002.  A fast and elitist multiobjective genetic algorithm: NSGA-II.  _IEEE transactions on evolutionary computation_ , 6(2): 182–197.
  * Drela (1989) Drela, M. 1989.  XFOIL: An Analysis and Design System for Low Reynolds Number Airfoils.  In _Conference on Low Reynolds Number Aerodynamics_ , 1–12.
  * Durasov et al. (2021) Durasov, N.; Bagautdinov, T.; Baque, P.; and Fua, P. 2021.  Masksembles for Uncertainty Estimation.  In _Conference on Computer Vision and Pattern Recognition_.
  * Durasov et al. (2019) Durasov, N.; Romanov, M.; Bubnova, V.; Bogomolov, P.; and Konushin, A. 2019.  Double Refinement Network for Efficient Monocular Depth Estimation.  In _2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 5889–5894.
  * Forti and Rozza (2014) Forti, D.; and Rozza, G. 2014.  Efficient geometrical parametrisation techniques of interfaces for reduced-order modelling: application to fluid–structure interaction coupling problems.  _International Journal of Computational Fluid Dynamics_ , 28(3-4): 158–169.
  * Gal and Ghahramani (2016a) Gal, Y.; and Ghahramani, Z. 2016a.  Dropout as a bayesian approximation: Representing model uncertainty in deep learning.  In _international conference on machine learning_ , 1050–1059. PMLR.
  * Gal and Ghahramani (2016b) Gal, Y.; and Ghahramani, Z. 2016b.  Dropout as a Bayesian Approximation: Representing Model Uncertainty in Deep Learning.  In _International Conference on Machine Learning_ , 1050–1059.
  * Gao, Wu, and Xia (2017) Gao, Y.; Wu, Y.; and Xia, J. 2017.  Automatic Differentiation Based Discrete Adjoint Method for Aerodynamic Design Optimization on Unstructured Meshes.  _Chinese Journal of Aeronautics_ , 30(2): 611–627.
  * Gosselin, Tye-Gingras, and Mathieu-Potvin (2009) Gosselin, L.; Tye-Gingras, M.; and Mathieu-Potvin, F. 2009.  Review of Utilization of Genetic Algorithms in Heat Transfer Problems.  _International Journal of Heat and Mass Transfer_ , 52(9): 2169–2188.
  * Guillard et al. (2021) Guillard, B.; Remelli, E.; Lukoianov, A.; Richter, S.; Bagautdinov, T.; Baque, P.; and Fua, P. 2021.  DeepMesh: Differentiable Iso-Surface Extraction.  _arXiv Preprint_.
  * Guo, Li, and Iorio (2016) Guo, X.; Li, W.; and Iorio, F. 2016.  Convolutional Neural Networks for Steady Flow Approximation.  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’16, 481–490. New York, NY, USA: Association for Computing Machinery.  ISBN 9781450342322.
  * He et al. (2016) He, K.; Zhang, X.; Ren, S.; and Sun, J. 2016.  Deep residual learning for image recognition.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 770–778.
  * Herbrich, Lawrence, and Seeger (2003) Herbrich, R.; Lawrence, N.; and Seeger, M. 2003.  Fast Sparse Gaussian Process Methods: The Informative Vector Machine.  In _Advances in Neural Information Processing Systems_ , 625–632.
  * Ivanov, Durasov, and Burnaev (2018) Ivanov, S.; Durasov, N.; and Burnaev, E. 2018.  Learning node embeddings for influence set completion.  In _2018 IEEE International Conference on Data Mining Workshops (ICDMW)_ , 1034–1037. IEEE.
  * Jacobs, Ward, and Pinkerton (1948) Jacobs, E.; Ward, K.; and Pinkerton, R. 1948.  The Characteristics of 78 Related Airfoil Sections from Tests in the Variable Density Wind Tunnel.  Technical report.
  * Jasak et al. (2007) Jasak, H.; Jemcov, A.; Tukovic, Z.; et al. 2007.  OpenFOAM: A C++ library for complex physics simulations.  In _International workshop on coupled methods in numerical dynamics_ , volume 1000, 1–20.
  * Jeong, Murayama, and Yamamoto (2005) Jeong, S.; Murayama, M.; and Yamamoto, K. 2005.  Efficient Optimization Design Method Using Kriging Model.  _Journal of Aircraft_ , 42(2): 413–420.
  * Kim et al. (2019) Kim, B.; Azevedo, V. C.; Thuerey, N.; Kim, T.; Gross, M.; and Solenthaler, B. 2019\.  Deep fluids: A generative network for parameterized fluid simulations.  In _Computer Graphics Forum_ , 59–70.
  * Lakshminarayanan, Pritzel, and Blundell (2017) Lakshminarayanan, B.; Pritzel, A.; and Blundell, C. 2017.  Simple and Scalable Predictive Uncertainty Estimation Using Deep Ensembles.  In _Advances in Neural Information Processing Systems_.
  * Laurenceau et al. (2010) Laurenceau, J.; Meaux, M.; Montagnac, M.; and Sagaut, P. 2010.  Comparison of Gradient-Based and Gradient-Enhanced Response-Surface-Based Optimizers.  _American Institute of Aeronautics and Astronautics Journal_ , 48(5): 981–994.
  * Liu, Durasov, and Fua (2021) Liu, W.; Durasov, N.; and Fua, P. 2021.  Leveraging Self-Supervision for Cross-Domain Crowd Counting.  _arXiv preprint arXiv:2103.16291_.
  * Lorensen and Cline (1987) Lorensen, W.; and Cline, H. 1987.  Marching Cubes: A High Resolution 3D Surface Construction Algorithm.  In _ACM SIGGRAPH_ , 163–169.
  * Mescheder et al. (2019) Mescheder, L.; Oechsle, M.; Niemeyer, M.; Nowozin, S.; and Geiger, A. 2019.  Occupancy networks: Learning 3d reconstruction in function space.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 4460–4470.
  * Mockus (2012) Mockus, J. 2012.  _Bayesian approach to global optimization: theory and applications_ , volume 37.  Springer Science & Business Media.
  * Monti et al. (2017) Monti, F.; Boscaini, D.; Masci, J.; Rodolà, E.; Svoboda, J.; and Bronstein, M. M. 2017.  Geometric Deep Learning on Graphs and Manifolds Using Mixture Model CNNs.  In _Conference on Computer Vision and Pattern Recognition_ , 5425–5434.
  * Mountrakis et al. (2015) Mountrakis, L.; Lorenz, E.; Malaspinas, O.; Alowayyed, S.; Chopard, B.; and Hoekstra, A. 2015.  Parallel Performance of an IB-LBM Suspension Simulation Framework.  _Journal of Computational Science_.
  * Park et al. (2019) Park, J. J.; Florence, P.; Straub, J.; Newcombe, R.; and Lovegrove, S. 2019.  Deepsdf: Learning continuous signed distance functions for shape representation.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 165–174.
  * Rasmussen and Williams (2006) Rasmussen, C. E.; and Williams, C. K. 2006.  _Gaussian Process for Machine Learning_.  MIT Press.
  * Redmon et al. (2016) Redmon, J.; Divvala, S.; Girshick, R.; and Farhadi, A. 2016.  You only look once: Unified, real-time object detection.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 779–788.
  * Remelli et al. (2020) Remelli, E.; Lukoianov, A.; Richter, S.; Guillard, B.; Bagautdinov, T.; Baque, P.; and Fua, P. 2020.  MeshSdf: Differentiable Iso-Surface Extraction.  In _Advances in Neural Information Processing Systems_.
  * Saviers, Ranjan, and Mahmoudi (2019) Saviers, K. R.; Ranjan, R.; and Mahmoudi, R. 2019.  Design and Validation of Topology Optimized Heat Exchangers.  In _AIAA Scitech 2019 forum_ , 1465.
  * Toal and Keane (2011) Toal, D.; and Keane, A. 2011.  Efficient Multipoint Aerodynamic Design Optimization via Cokriging.  _Journal of Aircraft_ , 48(5): 1685–1695.
  * Tompson et al. (2017) Tompson, J.; Schlachter, K.; Sprechmann, P.; and Perlin, K. 2017.  Accelerating eulerian fluid simulation with convolutional networks.  In _International Conference on Machine Learning_ , 3424–3433. PMLR.
  * Umetani and Bickel (2018) Umetani, N.; and Bickel, B. 2018.  Learning Three-Dimensional Flow for Interactive Aerodynamic Design.  _ACM Transactions on Graphics_ , 37(4): 89.
  * Vaswani et al. (2017) Vaswani, A.; Shazeer, N.; Parmar, N.; Uszkoreit, J.; Jones, L.; Gomez, A. N.; Kaiser, Ł.; and Polosukhin, I. 2017.  Attention is all you need.  In _Advances in neural information processing systems_ , 5998–6008.
  * Weller et al. (1998) Weller, H. G.; Tabor, G.; Jasak, H.; and Fureby, C. 1998.  A Tensorial Approach to Computational Continuum Mechanics Using Object-Oriented Techniques.  _Computational Physics_ , 12(6): 620–631.
  * Wen, Tran, and Ba (2020) Wen, Y.; Tran, D.; and Ba, J. 2020.  Batchensemble: an alternative approach to efficient ensemble and lifelong learning.  _arXiv preprint arXiv:2002.06715_.
  * Xu et al. (2017) Xu, G.; Liang, X.; Yao, S.; Chen, D.; and Li, Z. 2017.  Multi-Objective Aerodynamic Optimization of the Streamlined Shape of High-Speed Trains Based on the Kriging Model.  _PloS one_ , 12(1): 1–14.


# Distributed Combinatorial Optimization of Downlink User Assignment in mmWave
Cell-free Massive MIMO Using Graph Neural Networks

Bile Peng1, Bihan Guo1, Karl-Ludwig Besser 3, Luca Kunz1, Ramprasad
Raghunath1, Anke Schmeink2,
Eduard A. Jorswieck1, Giuseppe Caire4, H. Vincent Poor 3 Email: {b.peng,
bihan.guo, luca.kunz, r.raghunath<EMAIL_ADDRESS>{karl.besser<EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS>The work of B. Peng, R.
Raghunath, A. Schmeink, E. Jorswieck and G. Caire 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 16KISK031. The work of K.-L.
Besser is supported by the German Research Foundation (DFG) under grant BE
8098/1-1. The work of L. Kunz is supported by the BMBF under grant 16KISK074.
The work of H. V. Poor was supported in part by the U.S National Science
Foundation under Grants CNS-2128448 and ECCS-2335876. 1Institute for
Communications Technology, Technische Universität Braunschweig, Germany
3Department of Electrical and Computer Engineering, Princeton University, USA
2Chair of Information Theory and Data Analytics, RWTH Aachen University,
Germany  4Communications and Information Theory Group, Technische Universität
Berlin, Germany

###### Abstract

Millimeter wave (mmWave) cell-free massive MIMO (CF mMIMO) is a promising
solution for future wireless communications. However, its optimization is non-
trivial due to the challenging channel characteristics. We show that mmWave CF
mMIMO optimization is largely an assignment problem between access points
(APs) and users due to the high path loss of mmWave channels, the limited
output power of the amplifier, and the almost orthogonal channels between
users given a large number of AP antennas. The combinatorial nature of the
assignment problem, the requirement for scalability, and the distributed
implementation of CF mMIMO make this problem difficult. In this work, we
propose an unsupervised machine learning (ML) enabled solution. In particular,
a graph neural network (GNN) customized for scalability and distributed
implementation is introduced. Moreover, the customized GNN architecture is
hierarchically permutation-equivariant (HPE), i.e., if the APs or users of an
AP are permuted, the output assignment is automatically permuted in the same
way. To address the combinatorial problem, we relax it to a continuous
problem, and introduce an information entropy-inspired penalty term. The
training objective is then formulated using the augmented Lagrangian method
(ALM). The test results show that the realized sum-rate outperforms that of
the generalized serial dictatorship (GSD) algorithm and is very close to the
upper bound in a small network scenario, while the upper bound is impossible
to obtain in a large network scenario.

###### Index Terms:

Assignment, augmented Lagrangian method, cell-free massive MIMO, graph neural
network, unsupervised machine learning.

##  I Introduction

Cell-free massive MIMO (CF mMIMO) is a promising solution for next-generation
wireless communications. By deploying access points (APs) in a high density
and enabling cooperation between APs to serve the same user, CF mMIMO can
realize high received signal strength (RSS), low interference thanks to the
cooperation between APs, and high reliability due to the resulting macro
diversity [1]. These properties make CF mMIMO suitable for millimeter wave
(mmWave) communication, which enjoys a large bandwidth but suffers from low
RSS due to high propagation loss and low amplifier output power. In addition,
weak diffraction at high frequencies [2] makes blockage a serious problem. By
using the mmWave spectrum and cooperation of multiple APs in CF mMIMO, the RSS
is enhanced and the blockage problem is mitigated through macro-diversity.

Given the hardware and channel properties, optimizing mmWave CF mMIMO is
primarily an assignment problem between APs and users. Sophisticated precoding
techniques, such as zero-forcing (ZF) and minimum mean square error (MMSE),
are neither feasible (due to the immature hardware) nor indispensable (the
spectral efficiency is not crucial given a large bandwidth). Instead,
conjugate precoding is more feasible with simple hardware, realizes high
antenna gain given many antennas, and reduces interference significantly
because multi-user channels become asymptotically orthogonal as the number of
antennas increases [3]. Moreover, transmit power control is not necessary
because the optimal transmit power is always the maximum available power given
the high path loss for typical objectives.

In the literature, a mixed-integer programming problem for CF mMIMO fronthaul
resource allocation is considered in [4]. APs and pilots are assigned in a
rule-based, distributed way in [5]. An alternative selection method is
proposed for AP-user assignment in [6]. The block-sparsity norm $l_{2,1}$ is
applied to optimize AP-user assignment and power control simultaneously in
[7]. The pilot assignment problem is considered based on serving set in [8]
and graph coloring in [9].

It can be concluded from the above that the assignment problem in CF mMIMO
optimization is difficult to solve due to its non-differentiability. The
existing algorithms are either heuristic or do not scale well. In recent
years, machine learning (ML) has been proposed to solve combinatorial problems
in [10]. Such methods train a neural network (NN) on massive amounts of data
to find an optimized mapping from problem representation to solution. In the
literature, an encoder-decoder like NN is applied to optimize AP scheduling in
supervised manner in [11]. In [12], a “meta-algorithm” is proposed to generate
a greedy policy for combinatorial problems by joining reinforcement learning
(RL) and graph embedding. RL is applied in [13] to optimize pilot assignment
in a decentralized way. However, if we do not have “correct” labels and must
use unsupervised learning for optimization, trial-and-error based RL is not
necessarily the optimal choice, since it has a poor scalability due to the
“curse of dimensionality” [14]. On the other hand, if the objective is known,
model-based optimization is a better choice because it does not rely on
sampling solutions through trial-and-error. Therefore, the dimensionality has
only a moderate impact.

Among different contributions to CF mMIMO optimization with ML, the graph
neural network (GNN) architecture has gained particular attention due to its
distributed nature and high scalability [15, 16, 17, 18]. Despite its
significant potential, combinatorial assignment problems with connection
number constraints remain open to the best of our knowledge.

We propose a model-based ML approach to optimize the assignment between APs
and users using a customized GNN. Compared to existing works, the proposed
solution is unsupervised (i.e., it does not need “correct” labels) and model-
based. Our main contributions are summarized below.

  * •

We propose a customized hierarchically permutation-equivariant (HPE) [19] GNN,
i.e., any permutation of APs and users leads to the same permutation of the
solution. Moreover, data traffic in the fronthaul between APs (i.e., message
passing between nodes of the GNN) is reduced significantly compared to the
canonical GNN.

  * •

We relax the original combinatorial problem to a continuous one, design a
penalty inspired by information entropy to enforce discreteness, and apply the
augmented Lagrangian method (ALM) for training.

Beyond the AP-user assignment problem, this work is an early contribution to
unsupervised ML for general scalable combinatorial optimization problems.

##  II Problem Statement

###  II-A Problem Formulation

We consider a mmWave CF mMIMO network with $N$ APs and $K$ users, as
illustrated in Figure 1.

AP 1AP 2AP 3AP 4AP 5User 1User 2 Figure 1: The mmWave CF mMIMO network
scenario in the downlink, where dashed arrows indicate assignments between APs
and users.

In this CF mMIMO network, every AP is equipped with a large antenna array with
an analog beamformer, which performs conjugate precoding. Due to the narrow
beamwidth and height difference between APs and users, interference between
users is negligible even with similar azimuths of arrival111The elevations of
arrival are different for users with different distances to their common
serving AP.. Therefore, we make the pragmatic assumption that interference can
be omitted. Moreover, due to high propagation path loss at mmWave frequencies
and the limited output power of the amplifier, we assume that the optimal
transmit power is the maximum available power. This transforms the
conventional precoding and transmit power control problem into a combinatorial
assignment problem: The optimization variable $s_{kn}\in\\{0,1\\}$ indicates
whether AP $n$ serves user $k$ ($s_{kn}=1$) or not ($s_{kn}=0$). Since phase
synchronization is difficult for mmWave communications, we assume non-coherent
phases between APs. The capacity of user $k$ is therefore

$r_{k}=\log_{2}\left(1+\frac{\sum_{n=1}^{N}g_{kn}s_{kn}}{\sigma^{2}}\right),$
(1)

where $g_{kn}$ is the effective channel gain from AP $n$ to user $k$ including
transmit power, antenna gain and path loss, and $\sigma^{2}$ is the noise
power (we assume the noise is additive white and Gaussian).

We assume that each AP can serve at most $U$ users due to hardware
constraints, and each user must be served by at least $L$ APs to ensure
reliability via macro-diversity. Our objective is to maximize the sum rate.
The problem is formulated as

$\displaystyle\max_{\mathbf{S}}\quad$ $\displaystyle f=\sum_{k=1}^{K}r_{k}$
(2a) s.t. $\displaystyle s_{kn}\in\\{0,1\\},$ (2b)
$\displaystyle\sum_{k=1}^{K}s_{kn}\leq U,\quad\forall n$ (2c)
$\displaystyle\sum_{n=1}^{N}s_{kn}\geq L,\quad\forall k$ (2d)

where $\mathbf{S}\in\\{0,1\\}^{K\times N}$ has element $s_{kn}$ in row $k$ and
column $n$.

###  II-B Problem Properties

We identify two major properties of (2) that must be reflected in the proposed
solution as follows.

####  II-B1 Hierarchical Permutation-Equivariance

Problem (2) is _HPE_ : On the AP level, if we permute APs, the user assignment
of every AP should be permuted in the same way, as illustrated in Figure 2(a),
where $\mathbf{s}_{n}$ is the $n$-th column of $\mathbf{S}$. On the user level
of every AP, if we permute users, the user assignment of every AP should be
permuted in the same way, as illustrated in Figure 2(b). This property does
not exist in conventional NN architectures. In Section IV, we propose a GNN
with specially designed architectures for node and edge processing such that
the GNN is HPE.

AP 1InputAP 2AP 3$\mathbf{s}_{1}$Output$\mathbf{s}_{2}$$\mathbf{s}_{3}$AP
1InputAP 3AP 2$\mathbf{s}_{1}$Output$\mathbf{s}_{3}$$\mathbf{s}_{2}$

(a) Permutation-equivariance on AP level.

User 1InputUser 2User 3$s_{1n}$Output$s_{2n}$$s_{3n}$User 1InputUser 3User
2$s_{1n}$Output$s_{3n}$$s_{2n}$

(b) Permutation-equivariance on user level of AP $n$.

Figure 2: The hierarchical permutation-equivariance. If the inputs are
permuted, the outputs should be permuted in the same way automatically.

####  II-B2 Non-differentiable Combinatorial Optimization

Constraint (2b) makes (2) a combinatorial problem. With existing methods, this
problem can be tackled either through exhaustive search for small problem
sizes or via heuristic approaches with suboptimal performance. Neither of
these approaches is ideal. Problem (2) is also difficult to solve with NNs
because NN training is based on gradients. However, this problem is non-
differentiable. In Section V, we propose a training process that relaxes (2)
to a continuous one, and then uses a penalty and ALM to find solutions to the
original problem.

##  III General Approach of the Unsupervised Machine Learning for Optimization

We present the general approach of unsupervised ML for optimization as
follows. Given effective channel gains $\mathbf{G}\in\mathbb{R}^{K\times N}$,
with element $g_{kn}$ in row $k$ and column $n$ as defined in (1), we look for
a solution $\mathbf{S}$ that maximizes objective $f$ defined in (2a), which is
fully determined by $\mathbf{G}$ and $\mathbf{S}$, and can be written as
$f(\mathbf{G},\mathbf{S})$. We define an NN $N_{\theta}$, which is
parameterized by $\theta$, i.e., $\theta$ contains all trainable weights and
biases of $N_{\theta}$. $N_{\theta}$ maps from $\mathbf{G}$ to $\mathbf{S}$,
i.e., $\mathbf{S}=N_{\theta}(\mathbf{G})$. The objective can then be written
as $f(\mathbf{G},\mathbf{S})=f(\mathbf{G},N_{\theta}(\mathbf{G});\theta)$,
where it is emphasized that $f$ depends on $\theta$. We collect massive
amounts of realizations of $\mathbf{G}$ in a training set $\mathcal{D}$ and
formulate the problem as

$\max_{\theta}J=\sum_{\mathbf{G}\in\mathcal{D}}f(\mathbf{G},N_{\theta}(\mathbf{G});\theta).$
(3)

In this way, $N_{\theta}$ is optimized for the ensemble of
$\mathbf{G}\in\mathcal{D}$ (_training_) using gradient ascent:

$\theta\leftarrow\theta+\eta\nabla_{\theta}J,$ (4)

where $\eta$ is learning rate. If $N_{\theta}$ is well trained,
$\mathbf{S}^{\prime}=N_{\theta}(\mathbf{G^{\prime}})$ is also a good solution
for $\mathbf{G}^{\prime}\notin\mathcal{D}$ (_testing_), like a human uses
experience to solve new problems of the same type222A complete retraining is
only required when the input is fundamentally changed, e.g., when there is a
change of deployment environment. [20].

Compared to analytical methods, which optimize $\mathbf{S}$ itself given
$\mathbf{G}$, unsupervised ML optimizes the mapping from $\mathbf{G}$ to
$\mathbf{S}$. This has the following two advantages:

  1. 1.

Low complexity in operation: If analytical optimization is used, it is
performed after receiving $\mathbf{G}$. If its complexity is high, a real-time
application may be difficult. In contrast, with ML-based optimization,
$N_{\theta}$ is trained prior to operation. During operation, using the
trained $N_{\theta}$ to make inferences on a new $\mathbf{G}^{\prime}$ is much
less complex than training or analytical optimization.

  2. 2.

Scalability: As will be described in Section IV, the same mapping is applied
to all APs and users. This is a significant advantage in large network
scenarios with many APs and users. By optimizing the mapping $N_{\theta}$
instead of $\mathbf{S}$, the complexity becomes independent of the numbers of
APs and users. This achieves a high scalability, which is a major challenge in
CF mMIMO.

Formulation (3) presents a general approach. In the following sections, we
propose specialized GNN architecture and objective function according to the
problem properties described in Section II-B.

##  IV Hierarchically Permutation-Equivariant Graph Neural Network

The GNN structure is applied for CF mMIMO optimization due to its scalability
and decentralized nature. A GNN performs inference on a graph consisting of
nodes and edges, as illustrated in Figure 3. In our scenario, each node in the
graph represents an AP, and the edges between nodes indicate cooperation
between APs.

AP 1AP 2AP
3$\mathbf{e}_{12}$$\mathbf{e}_{21}$$\mathbf{e}_{31}$$\mathbf{e}_{13}$$\mathbf{e}_{23}$$\mathbf{e}_{32}$
Figure 3: An illustration of a graph. Each node in the graph represents an AP
and edges between nodes indicate cooperation between APs.

In layer $i$ of the GNN, inference at node $n$ is performed as

$\mathbf{F}_{n}^{(i+1)}=\gamma^{(i)}\left(\mathbf{F}_{n}^{(i)},\oplus_{m\in\mathcal{N}(n)}\phi^{(i)}\left(\mathbf{F}_{m}^{(i)},\mathbf{e}_{mn}\right)\right),$
(5)

where $\mathbf{F}_{n}^{(i+1)}$ is the node feature of node $n$ as the output
of layer $i$, $\gamma^{(i)}$ is an NN as node processing unit of layer $i$,
$\mathbf{F}_{n}^{(i)}$ is the node feature of node $n$ as the input of layer
$i$, $\oplus_{m\in\mathcal{N}(n)}$ is an aggregation function over all
neighboring nodes of node $n$ (in this work, averaging is used), $\phi^{(i)}$
is an NN as edge processing unit of layer $i$, $\mathbf{e}_{mn}$ is the
constant edge feature of the edge from node $m$ to node $n$. Inference (5) can
be interpreted as three steps: In Step 1, node $m$ uses $\phi^{(i)}$ to
generate a message for its neighboring node $n$ using its node feature
$\mathbf{F}_{m}^{(i)}$ and edge feature $\mathbf{e}_{mn}$. In Step 2, all
messages from neighbors of node $n$ are aggregated using operator $\oplus$. In
Step 3, node $n$ uses $\gamma^{(i)}$ to generate its node feature
$\mathbf{F}_{n}^{(i+1)}$ as output of layer $i$ according to
$\mathbf{F}_{n}^{(i)}$ and the aggregated message computed in Step 2.

Compared to the generic GNN model [21], a key difference in our model (5) is
that $\phi^{(i)}$ does not take $\mathbf{F}_{n}^{(i)}$ as an argument, such
that the message from node $m$ to node $n$ is computed with information
locally available at node $m$. Therefore, data traffic in fronthaul between
APs is significantly less than in the generic GNN model.

Note that $\gamma^{(i)}$ and $\phi^{(i)}$ are applied to all nodes and edges,
respectively. Therefore, the number of trainable GNN parameters is independent
of the number of nodes and edges, resulting in good scalability.

Although (5) ensures permutation-equivariance at the AP level, it does not
achieve permutation-equivariance at the user level when using a conventional
NN architecture, such as a fully connected NN, for $\phi$ and $\gamma$ [21].
In this work, we employ a dedicated architecture similar to the RISnet in
[22], which is permutation-equivariant to users. Its number of trainable
parameters is also independent of the number of users.

Both $\gamma$ and $\phi$ in all layers of the GNN (therefore without
superscripts) have two layers. The inputs of $\gamma$ and $\phi$ are two
dimensional matrices, where the first dimension is the feature and the second
dimension is the user, i.e., every user has a feature vector as a column of
the input matrix. Since the decision whether to serve a user depends on both
this user (the _current user_) and other users, we define an information
processing unit c for the current user and an information processing unit a
for all users. The output feature of Layer 1 for user $k$ is computed as

$\mathbf{f}_{k,2}=\begin{pmatrix}\text{ReLU}(\mathbf{W}^{\texttt{c}}_{1}\mathbf{f}_{k,1}+\mathbf{b}_{1}^{\texttt{c}})\\\
\big{(}\sum_{k^{\prime}}\text{ReLU}(\mathbf{W}^{\texttt{a}}_{1}\mathbf{f}_{k^{\prime},1}+\mathbf{b}_{1}^{\texttt{a}})\big{)}\big{/}K\\\
\end{pmatrix},$ (6)

where $\mathbf{W}_{1}^{\texttt{c}}$ and $\mathbf{b}_{1}^{\texttt{c}}$ are
weights and bias of information processing unit c on Layer 1, respectively,
and $\mathbf{W}_{1}^{\texttt{a}}$ and $\mathbf{b}_{1}^{\texttt{a}}$ are
weights and bias of information processing unit a on Layer 1, respectively. It
can be seen that the output for user $k$ using information processing unit c
depends only on the input feature of user $k$, whereas the output for user $k$
using information processing unit a is averaged over all users. In the second
layer, we do not use two information processing units to reduce data traffic
in fronthaul. The output feature of user $k$ is computed as

$\mathbf{f}_{k,3}=\text{ReLU}(\mathbf{W}_{2}\mathbf{f}_{k,2}+\mathbf{b}_{2}),$
(7)

where $\mathbf{W}_{2}$, $\mathbf{b}_{2}$ are weights and bias of Layer 2,
respectively. The above-described information processing of $\gamma$ and
$\phi$ is illustrated in Figure 4.

UsersInputInput matrixInfo. proc.Output matrixcaUsersFeature cFeature a

(a) First layer

UsersFeature cFeature aInput matrixInfo. proc.Output matrix Info. proc.
UsersOutput

(b) Second layer

Figure 4: Information processing (info. proc.) of $\gamma$ and $\phi$.

We let every AP make $U$ decisions (serving at most $U$ users according to
(2c)). Following the principle of recurrent neural networks (RNNs), we run the
GNN $U$ times. At the $u$-th run, a relaxed decision of assignment
$s_{kn}^{(u)}\in[0,1]$ between AP $n$ and user $k$ is made. The activation
function softmax is applied in the final layer such that
$s_{kn}^{(u)}\in[0,1]$ for all $n$, $k$ and $u$, and
$\sum_{k=1}^{K}s_{kn}^{(u)}=1$ for all $n$ and $u$. The input of user $k$ and
AP $n$ at the $u$-th run comprises the channel gain $g_{kn}$ and gap to
fulfill constraint (2d) after $u$ runs, i.e.,
$\text{ReLU}(L-\sum_{\mu=1}^{u}s_{kn}^{(\mu)})$.

##  V Training with Augmented Lagrangian Method

In this section, we propose a training process for Problem (2) using ALM.
Unlike conventional applications of ALM, which tune optimization variables
directly, we use ALM to formulate the objective, then use the Adam optimizer
to optimize the NN parameter $\theta$. During training, $s_{kn}$ is relaxed to
a continuous variable between 0 and 1. Constraint (2c) is guaranteed by the
recurrent process described in Section IV. Constraints (2b) and (2d) appear in
the objective as penalty terms. Using ALM, the objective function is defined
as

$\begin{split}g&=f-\lambda_{1}\sum_{k=1}^{K}\text{ReLU}\left(L-\sum_{n=1}^{N}s_{kn}\right)-\frac{1}{2}\lambda_{2}\sum_{n=1}^{N}p_{n}\\\
&\qquad-\nu_{1}\sum_{k=1}^{K}\text{ReLU}\left(L-\sum_{n=1}^{N}s_{kn}\right)^{2}-\frac{1}{2}\nu_{2}\sum_{n=1}^{N}p_{n}^{2},\end{split}$
(8)

where $\lambda_{j}$ and $\nu_{j}$ are the factors of the linear and quadratic
penalty terms, respectively ($j=1,2$), $p_{n}$ is the penalty term of AP $n$
for (2b) inspired by the information entropy:

$p_{n}=-\sum_{k=1}^{K}s_{kn}\log(s_{kn}).$ (9)

Note that $p_{n}$ is equal to 0 if $s_{kn}=0$ or $1$. Otherwise $p_{n}>0$.

The training is described in Algorithm 1, where $\Delta\nu$ is a predefined
constant changing rate of $\nu_{1}$ and $\nu_{2}$.

Algorithm 1 Training using ALM

1:Set $\lambda_{j}=\nu_{j}=0$ for $j=1,2$.

2:repeat

3: Update GNN with a gradient ascent step.

4:until Convergence

5:repeat

6:
$\lambda_{1}\leftarrow\lambda_{1}+\nu_{1}\sum_{k=1}^{K}\text{ReLU}\left(L-\sum_{n=1}^{N}s_{kn}\right)$

7: $\nu_{1}\leftarrow\nu_{1}+\Delta\nu$

8: repeat

9: Update GNN with a gradient ascent step.

10: until Convergence

11:until $\sum_{k=1}^{K}\text{ReLU}\left(L-\sum_{n=1}^{N}s_{kn}\right)=0$.

12:repeat

13:
$\lambda_{2}\leftarrow\lambda_{2}+\nu_{2}\sum_{k=1}^{K}\text{ReLU}\left(L-\sum_{n=1}^{N}s_{kn}\right)$

14: $\nu_{2}\leftarrow\nu_{2}+\Delta\nu$

15: repeat

16: Update GNN with a gradient ascent step.

17: until Convergence

18:until $-\sum_{k=1}^{K}s_{kn}\log(s_{kn})=0$.

Note that we added constraints (2d) and (9) one after another. This is
different from canonical ALM, which updates $\lambda_{j}$ and $\nu_{j}$
simultaneously for all constraints. However, training results show that the
proposed method guarantees the constraints better.

##  VI Training and Testing Results

For evaluating the training and testing performance, we define a small and a
large network scenario, where APs are uniformly placed (deterministically) and
users are randomly distributed, as shown in Figure 5. The user positions in
the training and testing sets are generated independently but with the same
uniform distributions. The path loss is assumed to be Rician distributed,
where the expectation is reciprocal to the distance from AP to user, and the
variance is constant. Important scenario parameters are listed in Table I.

$20$$40$$60$$80$$20$$40$$60$$80$APUser (a) 5 APs and 4 users

$0$$200$$400$$600$$800$$1{,}000$$0$$200$$400$$600$$800$$1{,}000$ (b) 20 APs
and 15 users

Figure 5: Visualization of AP-user assignment in two scenarios. Table I: Important scenario parameters. Parameter |  Value
---|---
Number of APs |  20 (large network scenario),
5 (small network scenario)
Number of users |  15 (large network scenario),
4 (small network scenario)
Minimum number of serving APs per user $L$ |  2
Maximum number of served users per AP $U$ |  2

We perform testing on a testing set (which is disjoint from the training set)
after every iteration of training. Figure 6 shows realized curves of the sum
rate for both training and testing phases. The sum rate in the testing phase
improves quickly to
$3.70\text{\,}\mathrm{bit}\text{/}\mathrm{Hz}\text{/}\mathrm{s}$ without
constraints, and then decreases slightly to
$3.67\text{\,}\mathrm{bit}\text{/}\mathrm{Hz}\text{/}\mathrm{s}$ due to the
constraints. The training curve fluctuates while the testing curve remains
smooth due to the fact that a random batch is sampled from the training set in
each iteration, whereas the same testing set is used for every testing
iteration. Despite the fluctuation in the training curve, the good match
between the training and testing curves on average indicates that there is no
overfitting. This is because the same information processing units are applied
to all AP-user connections (see Section IV). Therefore, the GNN is relatively
small compared to the problem size.

$0$$10{,}000$$20{,}000$$30{,}000$$40{,}000$$50{,}000$$60{,}000$$70{,}000$$80{,}000$$90{,}000$$100{,}000$$1$$2$$3$$4$IterationSum
rate (bit/Hz/s)TrainingTesting Figure 6: Sum-rate in training and testing over
iterations.

Figure 7 shows the realized connection penalty of users (2d) during both
training and testing. After introducing the penalty for constraint (2d) (i.e.,
Line 5 in Algorithm 1), the connection penalty is reduced to 0. After
introducing the penalty for discreteness (9) (i.e., Line 11 in Algorithm 1),
the connection penalty slightly increases before slowly decreasing to 0. We
confirm that constraint (2d) is satisfied in all 1024 testing samples at the
end of training.

$0$$10{,}000$$20{,}000$$30{,}000$$40{,}000$$50{,}000$$60{,}000$$70{,}000$$80{,}000$$90{,}000$$100{,}000$$0$$5$$10$IterationConnection
penaltyTrainingTesting Figure 7: Penalty term of constraint (2d) in training
and testing over iterations.

The training and testing curves of the discreteness penalty (9) are shown in
Figure 8. After introducing the penalty for discreteness (i.e., Line 11 in
Algorithm 1), the discreteness penalty reduces to 0. Therefore, the solution
to the relaxed problem is valid for the original problem (2).

$0$$10{,}000$$20{,}000$$30{,}000$$40{,}000$$50{,}000$$60{,}000$$70{,}000$$80{,}000$$90{,}000$$100{,}000$$0$$20$$40$IterationDiscreteness
penaltyTrainingTesting Figure 8: Discreteness penalty in training and testing
over iterations.

Two visualizations of small and large network scenarios are shown in Figure 5,
with dashed lines indicating the AP-user assignment. The assignments match our
intuition because users are always served by nearby APs, subject to
constraints (2c) and (2d).

Table II presents a comparison of sum spectral efficiency obtained with the
proposed method, an upper bound (found via exhaustive search), the average
performance of 100 random assignments, and the generalized serial dictatorship
(GSD) algorithm [23]. When performing an exhaustive search, we evaluate
$\binom{K}{U}^{N}$ possible assignments while considering constraints (2b) and
(2c). This results in $7776$ and $2.6\times 10^{40}$ assignments in the small
and large network scenarios, respectively. The proposed method achieves a
performance very close to the upper bound in the small network scenario,
whereas the upper bound cannot be calculated in the large network scenario due
to limited memory and time. In both scenarios, the proposed method outperforms
the random assignment and GSD algorithm. Additionally, the proposed method
provides a solution within a few milliseconds with reduced data traffic in
fronthaul (see Section IV), and without a centralized processing unit,
indicating potential for real-time application.

Table II: Achieved sum spectral efficiencies (bit/Hz/s) of proposed method and baselines. | Small scenario | Large scenario
---|---|---
Proposed | 1.14 | 3.67
Upper bound | 1.15 | Not available
Random | 0.60 | 0.51
GSD | 1.06 | 3.34

##  VII Conclusion

Due to challenging mmWave channel and hardware characteristics, the mmWave CF
mMIMO optimization problem is primarily an assignment problem between APs and
users. In this work, we have proposed an unsupervised ML approach. In
particular, we have utilized a customized GNN to achieve scalability,
distributed inference, hierarchical permutation-equivariance, and reduced data
traffic in fronthaul. To solve the combinatorial assignment problem, we first
relaxed it to a continuous problem. Next, we introduced a penalty inspired by
information entropy and use ALM to enforce discreteness. Testing results
indicate that the proposed method outperforms both random assignment and GSD
algorithm. In a small network scenario, its performance is very close to the
upper bound found via exhaustive search. In a large network scenario, where
exhaustive search is not feasible due to memory and time limitations, the
proposed method still achieves good performance. The proposed method could be
generalized to other large-scale combinatorial problems in the future. The
combination of ML and analytical ALM can also be applied to other constrained
optimization problems.

Code and data of this work will be made public on Github if the paper is
accepted.

## References

  * [1] Hussein A Ammar, Raviraj Adve, Shahram Shahbazpanahi, Gary Boudreau and Kothapalli Venkata Srinivas  “User-centric cell-free massive MIMO networks: A survey of opportunities, challenges and solutions”  In _IEEE Commun. Surv. Tutor._ 24.1 IEEE, 2021, pp. 611–652
  * [2] Usman Tahir Virk and Katsuyuki Haneda  “Modeling human blockage at 5G millimeter-wave frequencies”  In _IEEE Trans. Antennas Propag._ 68.3 IEEE, 2019, pp. 2256–2266
  * [3] Thomas L Marzetta, Erik G Larsson and Hong Yang  “Fundamentals of Massive MIMO”  Cambridge University Press, 2016
  * [4] Zhiyang Li, Fabian Göttsch, Siyao Li, Ming Chen and Giuseppe Caire  “Joint Fronthaul Load Balancing and Computation Resource Allocation in Cell-Free User-Centric Massive MIMO Networks”, 2023  arXiv:2310.14911 [cs.IT]
  * [5] Emil Björnson and Luca Sanguinetti  “Scalable cell-free massive MIMO systems”  In _IEEE Trans. Commun._ 68.7 IEEE, 2020, pp. 4247–4261
  * [6] Manobendu Sarker and Abraham O Fapojuwo  “Access Point-User Association and Auction Algorithm-Based Pilot Assignment Schemes for Cell-Free Massive MIMO Systems”  In _IEEE Syst. J._ 17.3 IEEE, 2023, pp. 4301–4312  DOI: 10.1109/JSYST.2023.3255865
  * [7] Thang X Vu, Symeon Chatzinotas, Shahram ShahbazPanahi and Bjorn Ottersten  “Joint power allocation and access point selection for cell-free massive MIMO”  In _2020 IEEE Int. Conf. Commun. (ICC)_ , 2020  IEEE
  * [8] Mahmoud Zaher, Emil Björnson and Marina Petrova  “Soft handover procedures in mmWave cell-free massive MIMO networks”  In _IEEE Trans. Wirel. Commun._ IEEE, 2023  DOI: 10.1109/TWC.2023.3330199
  * [9] Heng Liu, Jiayi Zhang, Shi Jin and Bo Ai  “Graph coloring based pilot assignment for cell-free massive MIMO systems”  In _IEEE Trans. Veh. Technol._ 69.8 IEEE, 2020, pp. 9180–9184
  * [10] Yoshua Bengio, Andrea Lodi and Antoine Prouvost  “Machine learning for combinatorial optimization: a methodological tour d’horizon”  In _Eur. J. Oper. Res._ 290.2 Elsevier, 2021, pp. 405–421
  * [11] Mamoun Guenach, Ali A Gorji and André Bourdoux  “A deep neural architecture for real-time access point scheduling in uplink cell-free massive MIMO”  In _IEEE Trans. Wirel. Commun._ 21.3 IEEE, 2021, pp. 1529–1541
  * [12] Elias Khalil, Hanjun Dai, Yuyu Zhang, Bistra Dilkina and Le Song  “Learning combinatorial optimization algorithms over graphs”  In _Adv. Neural Inf. Process. Syst._ 30, 2017, pp. 6348–6358
  * [13] Myeung Suk Oh, Anindya Bijoy Das, Seyyedali Hosseinalipour, Taejoon Kim, David J Love and Christopher G Brinton  “A decentralized pilot assignment algorithm for scalable O-RAN cell-free massive MIMO”  In _IEEE J. Sel. Areas Commun._ 42.2 IEEE, 2024, pp. 373–388  DOI: 10.1109/JSAC.2023.3336154
  * [14] Michel Verleysen and Damien François  “The Curse of Dimensionality in Data Mining and Time Series Prediction”  In _Computational Intelligence and Bioinspired Systems_ , Lecture Notes in Computer Science 3512  Springer, 2005, pp. 758–770  DOI: 10.1007/11494669˙93
  * [15] Vismika Ranasinghe, Nandana Rajatheva and Matti Latva-aho  “Graph neural network based access point selection for cell-free massive MIMO systems”  In _2021 IEEE Glob. Commun. Conf. (GLOBECOM)_ IEEE, 2021
  * [16] Bohan Li, Lie-Liang Yang, Robert G Maunder, Songlin Sun and Pei Xiao  “Heterogeneous graph neural network for power allocation in multicarrier-division duplex cell-free massive MIMO systems”  In _IEEE Trans. Wirel. Commun._ 23.2 IEEE, 2024, pp. 962–977  DOI: 10.1109/TWC.2023.3284263
  * [17] Ramprasad Raghunath, Bile Peng and Eduard A Jorswieck  “Energy-Efficient Power Allocation in Cell-Free Massive MIMO via Graph Neural Networks”  In _2024 IEEE Int. Conf. Mach. Learn. Commun. Netw. (ICMLCN)_ IEEE, 2024
  * [18] Lou Salaün, Hong Yang, Shashwat Mishra and Chung Shue Chen  “A GNN Approach for Cell-Free Massive MIMO”  In _2022 IEEE Glob. Commun. Conf. (GLOBECOM)_ IEEE, 2022, pp. 3053–3058
  * [19] Ryo Umagami, Yu Ono, Yusuke Mukuta and Tatsuya Harada  “HiPerformer: Hierarchically Permutation-Equivariant Transformer for Time Series Forecasting”, 2023  arXiv:2305.08073 [cs.LG]
  * [20] Wei Yu, Foad Sohrabi and Tao Jiang  “Role of Deep Learning in Wireless Communications”  In _IEEE BITS Inf. Theory Mag._ 2.2 IEEE, 2022, pp. 56–72  DOI: 10.1109/MBITS.2022.3212978
  * [21] Matthias Fey and Jan Eric Lenssen  “Fast graph representation learning with PyTorch Geometric”, 2019  arXiv:1903.02428 [cs.LG]
  * [22] Bile Peng, Karl-Ludwig Besser, Ramprasad Raghunath, Vahid Jamali and Eduard A. Jorswieck  “RISnet: a Scalable approach for reconfigurable intelligent surface optimization with partial CSI”  In _2023 IEEE Glob. Commun. Conf. (GLOBECOM)_ IEEE, 2023, pp. 4810–4816  DOI: 10.1109/GLOBECOM54140.2023.10437049
  * [23] Katarína Cechlárová, Pavlos Eirinakis, Tamás Fleiner, Dimitrios Magos, Ioannis Mourtos and Eva Potpinková  “Pareto optimality in many-to-many matching problems”  In _Discret. Optim._ 14 Elsevier, 2014, pp. 160–169