# Active learning and molecular dynamics simulations to find high melting
temperature alloys

David E. Farachea,∗, Juan C. Verduzcoa,∗, Zachary D. McClurea,
Saaketh Desaib, and Alejandro Strachana

( a School of Materials Engineering and Birck Nanotechnology Center
Purdue University, West Lafayette, Indiana, 47907 USA
b Center for Integrated Nanotechnologies
Sandia National Laboratories, Albuquerque, New Mexico, 87123 USA
* Authors contributed equally to this work. )

###### Abstract

Active learning (AL) can drastically accelerate materials discovery; its power
has been shown in various classes of materials and target properties. Prior
efforts have used machine learning models for the optimal selection of
physical experiments or physics-based simulations. However, the latter efforts
have been mostly limited to the use of electronic structure calculations and
properties that can be obtained at the unit cell level and with negligible
noise. We couple AL with molecular dynamics simulations to identify multiple
principal component alloys (MPCAs) with high melting temperatures. Building on
cloud computing services through nanoHUB, we present a fully autonomous
workflow for the efficient exploration of the high dimensional compositional
space of MPCAs. We characterize how uncertainties arising from the stochastic
nature of the simulations and the acquisition functions used to select
simulations affect the convergence of the approach. Interestingly, we find
that relatively short simulations with significant uncertainties can be used
to efficiently find the desired alloys as the random forest models used for AL
average out fluctuations.

Keywords \- multiple principal component alloys (MPCAs); active learning;
uncertainty quantification;

##  1 Introduction

The combination of experiments, physics-based simulations, and data science
tools has been shown to have the potential to accelerate discovery of novel
materials with optimized properties [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12].
Active learning (AL), a subset of machine learning in which models learn
dynamically, has gathered significant interest, both from the point of basic
science [13] and for commercial applications [14]. AL models analyze existing
data and formulate queries to acquire additional information towards a design
goal. AL workflows start with a model trained with an initial set of data and
evaluate the expected gain towards an objective function of all possible new
experiments within a design space [15]. Top candidates are characterized (e.g.
by performing an experiment) and the outcome is added to the existing dataset.
Given this additional information, a new cycle is started. With this iterative
process, illustrated in Figure 1, the model becomes more accurate in regions
of interest within the design space. In order to identify the next best query,
AL uses selection strategies known as information acquisition functions. These
strategies differ from each other in the relative balance between exploitation
and exploration included in their mathematical formulation. Exploitation
favors cases expected to maximize the objective function while exploration
focuses on areas of high uncertainty.

Figure 1: Schematic representation of the AL iterative optimization. Starting
from existing data, a predictive model is trained with machine learning. The
model is then applied to all candidate materials in the design space to assess
their performance, including uncertainties. An acquisition function is used to
select the next material(s) to be tested (either via experiments or
simulations). If the new material does not satisfy the design conditions, it
is added back to the existing data set and the cycle re-initiated.

AL has been used for a wide range of applications from natural language
processing [16], reaction screening for pharmaceutical applications [1], and
multiscale modeling [2, 17]. In materials science, AL has been used to
accelerate discovery of materials with desired properties by coupling it with
experiments [18, 4, 19, 20, 21, 22] and physics-based simulations [5, 6, 7].
In addition, AL workflows paired with existing closed data sets has shown the
ability of these models to reduce the number of queries needed to identify the
best candidate [9, 10, 11, 12].

Kusne et al. [3] developed a closed-loop autonomous algorithm for materials
discovery and used it to explore the Ge-Sb-Te ternary system in search for an
optimal phase-change material with the highest optical bandgap difference
between phases. With an iterative approach using experimental ellipsometry
data, they identified a new composition, Ge4Sb6Te7, with nearly three times
the optical bandgap of Ge2Sb2Te5, a material used in random-access memory
devices. Xue et al. [4] demonstrated an adaptive design strategy for the
search of NiTi-based shape memory alloys with low thermal hysteresis. Starting
with 26 materials synthesized their optimization algorithm was used to screen
a space of $\sim$800,000 candidate materials. This screening resulted in 36
alloys queried and tested using differential scanning calorimetry (DSC) with
several showing lower thermal hysteresis than the starting training set. For
additional examples, see Refs. [19, 20, 21, 22].

Coupling AL workflows with physics-based simulations can further reduce the
costs and time required to discover new materials by minimizing the number of
experimental trials required. Tran and Ulissi [5] explored a high dimensional
space of density functional theory (DFT) results on intermetallics in search
for optimized catalyst materials for CO2 reduction. Their approach produced a
set of 54 promising materials from $\sim$1499 candidates in the Materials
Project [6]. Out of these materials, some have been since confirmed
experimentally. Seko et al. [7] reported a virtual screening of 54,779
materials looking for compounds of low lattice thermal conductivity. They used
AL to select appropriate descriptors for a Gaussian process regression (GPR)
model used to get information on the large dataset. Their screening yielded
221 materials with very low conductivity. After some constraints were imposed,
19 of them were characterized using DFT calculations for possible applications
as thermoelectric materials.

Work involving the combination of AL or similar optimization methods and
materials simulations has been mostly limited to DFT calculations and to
properties that can be obtained from relatively small simulation cells without
noise. While DFT often offers a good balance between accuracy and
computational efficiency, it remains limited to a subset of materials
properties. Fortunately, materials models across scales are available [23] and
significant progress in their coupling across scales has occurred over the
last few decades. Examples range from crystal plasticity models [24, 10] to
interatomic potentials from ab-initio simulations [25, 26, 27, 28, 29, 17].
Physics-based materials models across scales open the possibility of
significantly expanding the reach of AL approaches.

Molecular dynamics (MD) is an important rung in multiscale modeling that
connects ab initio physics to the meso- and macro-scales. MD has been paired
with optimization methods like reinforcement learning, including Monte Carlo
tree search (MCTS) [30, 31]. For example, Patra et al. [32] paired MCTS with
MD to find a copolymer compatibilizer, and Loeffler et al. [33] applied the
methods to study defect structures in metal dichalcogenides.

In this work, we couple AL with MD simulations to find multiple principal
component alloys (MPCAs) with high predicted melting temperature for a model
CrCoCuFeNi system. The main objective of this paper is to demonstrate the
potential of coupling of AL workflows with MD simulations and characterizing
the effect of the stochastic and noisy nature of MD results on AL workflows.
As discussed in Section 2.2, given the uncertainties associated with the
predictions of melting temperature with the current potential, the composition
with the highest melting temperature predicted by the model may not match
experiments. This paper is organized as follows: Section 2 introduces the
overall approach, design space, and simulation details. Section 3 describes
the AL runs exploring the effects of (a) different information acquisition
functions and (b) different MD simulation times. It also discusses how
uncertainty of the model evolves as the active learning workflow progresses.
In Section 4 we summarize and draw final conclusions of the paper. Finally, in
Section 5 we discuss availability of our models, datasets and workflows
through nanoHUB.

##  2 Methods

###  2.1 Problem statement

MPCAs are a new class of materials, which includes high entropy and complex
concentrated alloys, where four of more elements are combined in nearly equal
atomic percentages. Interest in these materials has grown due to several
desirable properties. Properties such as high strength at elevated
temperatures, radiation resistance [34], and high melting temperatures make
some MPCAs attractive for high temperature and extreme applications [35, 36,
37]. Experimental determination of this melting temperature is challenging
because of their complex processing [38] protocols, chemical reactions, and
phase separation [39, 40]. For these reasons, using computational methods to
determine the melting temperature of these alloys is highly desirable.
Fortunately, melting temperature can be accurately obtained from first
principles, and the first calculations for MPCAs are emerging [41].

The high-dimensional space of possible compositions prevents brute-force
approaches even for relatively cheap simulations. Therefore, an efficient way
to explore this space is essential. In this work, we pair AL with MD
simulations to find MPCAs with the highest possible temperature in a model
system. AL in the context of materials design seeks to reduce the number of
experiments needed to find an optimal candidate out of a pool of untested
contenders in an unexplored space. Within the AL scheme, predictions of the
surrogate machine learning model, along with sample-wise uncertainties, are
used to identify the next candidate to be tested, via acquisition functions.
In this work, we make use of the FUELS (Forests with Uncertainty Estimates for
Learning Sequentially) framework proposed by Ling et al. [11] based on random
forests and a supervised machine learning algorithm.

####  2.1.1 Initial and candidate space

We explore 5-component FCC alloys incorporating Cr, Co, Cu, Fe, and Ni.
Iterative AL models require an initial subset of entries to be evaluated to
serve as the prior for the ML model. Our initial set of alloys was selected to
be around the equiatomic composition Cr20Co20Cu20Fe20Ni20. Alloys included in
this initial set deviate by 10 at.% from this composition, such that each
element composes 10 at.% to 30 at.% of the alloy. Using steps of 10 at. %, we
get 39 data points. Melting temperatures were calculated via MD simulations,
as described below, using the two-phase coexistence method. This initial set
was the starting point for all of our AL runs.

For our exploration space, we relaxed the constraints from the initial set to
allow more deviations from the equiatomic such that any of the elements is
between 0 at.% and 50 at.%. This allows for exploration of binaries, ternaries
and quaternaries and ultimately, using steps of 10 at. %, results in 554 data
points that are not included in our initial subset. We recognize that this
grid spacing is coarse, but it is appropriate for an initial exploration. Our
iterative AL model will start with training on our initial space, and the
candidate space will serve as the pool from which it will be drawing
candidates.

###  2.2 MD simulations of melting

To predict the melting temperature of FCC MPCAs, we made use of MD simulations
and the two-phase coexistence method [42, 43]. To describe atomic
interactions, we used a many-body, embedded-atom-method interatomic potential
developed by Farkas and Caro [44]. The potential was designed for FCC Cr-Co-
Cu-Fe-Ni alloys with near equiatomic compositions. Its parameterization
focused on reproduction of elastic constants, vacancy formation energy,
stacking fault energy, surface energies and relative phase stability to ensure
that FCC is the stable structure for all five components. However, this
potential has not been trained or tested for melting temperature calculations.

Several techniques exist to compute the melting temperature of materials using
MD simulations. The most direct approach involves heating a crystalline sample
until it melts and then cooling the melt until it recrystallizes. This results
in overheating and undercooling, requiring a-posteriori corrections [45]. To
avoid these issues, we predict melting temperatures using a phase coexistence
method [42, 43]. In this method, a temperature at which the crystal and liquid
phases of the material of interest coexist is established.

The first step is to generate an initial sample with a liquid and solid phase
in contact, we achieve this with the use of two thermostats during the sample
preparation step. We start from the four-atom FCC unit cell with a lattice
parameter of 3.56 Å and replicate it 8 x 8 x 18 to create a periodic
simulation cell. The third direction (z) is normal to the solid-liquid
interface. Atom types are assigned based on the desired composition randomly
to each lattice site. The entire structure is relaxed via energy minimization
with respect to cell parameters and atomic positions. To create initial
liquid/solid structure the cell is then divided in halves along z and two
thermostats are used to create the two phases. The two regions are
equilibrated at their initial temperatures for 10 ps using isothermal,
isobaric MD simulations (NPT ensemble) with a Berendsen thermostat [46] and
Nose-Hoover barostat [47]. This first step seeks to create an initial cell
with both liquid and solid and the two temperatures need to chosen
appropriately (the temperature of the liquid half should be higher than the
melting temperature). A snapshot of the system after this initial step is
shown in Figure 2. Clearly, this initial sample is not in equilibrium, and we
follow the initial step with an isoenthalpic, isobaric (NPH) simulation where
the liquid and solid regions are allowed to exchange energy freely and come to
equilibrium. After steady state is achieved and the simulation cell has a
uniform temperature, if the sample contains both liquid and solid the
simulation temperature corresponds to the melting temperature of the system.
The use of the NPH ensemble allows the temperature of the system to evolve
towards the melting temperature. If the sample is initially too hot the liquid
phase will grow at the expense of the solid, the heat of fusion will
automatically cool the sample down towards the melting temperature; the
reverse happens if the sample is initially too cold.

Figure 2: MD snapshot of the simulation cell divided into a liquid and solid
region for an equiatomic alloy of Cr,Co,Cu,Fe and Ni. Each color represents a
different element.

A successful run requires the coexistence of liquid and solid after steady
state has been achieved. These conditions are checked by analyzing the atomic
structure and the time evolution of the instantaneous temperature of the
system. Phase fractions of solid and liquid are calculated using the
polyhedral template matching (PTM) algorithm from OVITO [48, 49]. Coexistence
is defined as the fraction of liquid being between 35 and 65 at.% and the FCC
crystal also between 35 and 65 at.% in the final snapshot of the simulation.
We also analyze thermodynamic data from the last 60% of the simulation to
assess steady state, we require the absolute value of the slope of the
instantaneous temperature vs. time to be less than 1 K/ps. Finally, we also
compute the 95% confidence interval on the temperature calculation. Melting
temperatures for our initial set ran with 100ps of simulation time.

###  2.3 Uncertainties and noise in the MD predictions

The description of interatomic forces is at the heart of MD simulations and
determines the predicted melting temperatures. The potential used in this work
was not developed or widely tested for melting temperature calculations and we
are unaware of any potential that has. Thus, our focus is on understanding the
use of MD in ML workflows and not in accurate predictions of melting
temperatures. Nevertheless, for completeness, Figure 3 compares the predicted
melting temperatures with experimental values [50] for various alloys and
single-elements [51]. Our simulations assume random FCC alloys since that is
what the AL workflow uses, this does not match the experiments in all cases.
We included experimental melting temperatures for single elements for
comparison. All MD simulations all arrange structures as FCC. For iron, we
show the predicted melting temperature of FCC Fe since the value is relevant
for the high-melting temperature alloys found by the AL workflow.
Interestingly, for Cobalt (Co), the MD simulation with this potential achieved
coexistence by creating a solid BCC structure rather than FCC. We find that
there is a consistent overestimation of melting temperatures, except for Cu
and Cr.

Figure 3: Comparison for MD simulated and experimental melting temperatures
for alloys composed of elements included in the potential. Dashed line
indicates a match between the values. Experimental values for MPCAs taken from
Wu et al [50]. Experimental values for single-elements taken from NIST
[51].Filled symbols represent alloys reported as FCC crystal structure, open
symbols represent non-FCC crystal structures.

It is important to note that the generation of the initial atomic structure is
done stochastically, and temperature is obtained as the time-average of a
fluctuating quantity. Thus, the resulting melting temperature from independent
runs will vary due to sample-to-sample variability (which can be reduced
increasing the simulation cell size) and time averaging of instantaneous
temperature (which can be reduced by running longer simulations and also by
increasing system size). One of the goals of this paper is to assess what
level of noise can be tolerated in an AL workflow.

To quantify the variability in predicted melting temperatures, we chose one of
our optimal alloy candidates. We performed 36 simulations with the same
overall composition but independent atomic configurations and initial
velocities. Figure 4 shows the resulting distributions of predicted melting
temperatures for the three simulation times used in this study, the inset
shows the mean and uncertainty estimates. Figure 4 shows that, as expected,
repeated experiments result in mean values relatively independent of the
simulation time (2470-2480 K) and reducing the simulation time results in
broader distributions (remember that temperature is associated with the time-
averaged kinetic energy per degree of freedom).

Figure 4: Distributions of MD simulated melting temperatures for a MPCA of 50%
Fe and 50% Co at different simulation times. Distributions indicate a mean
temperature around 2473 K and narrow down as the simulation time increases.
Table 1: Mean and uncertainty estimates for the distributions in Figure 4 of
MD simulated melting temperatures for a MPCA of 50% Fe and 50% Co at different
simulation times. These distributions were constructed by running 36
simulations with combinations of 6 different atomic arrangements and 6
different atomic velocities initializations.

Time Per Simulation (ps) |  Mean MD Simulated Tm (K) |  Standard deviation MD Simulated Tm (K)
---|---|---
50 | 2479 | 17
100 | 2476 | 13
200 | 2470 | 8

###  2.4 Random forests and uncertainties

Several regression models, including neural networks, gaussian processes, and
random forests have been used in AL workflows. Neural networks have stark
advantages in image recognition and other high-dimensionality problems and
have seen widespread adoption in a myriad of problems and disciplines.
However, uncertainty quantification on their applications for regression and
classification is remains an active area of research [52]. Random forests, on
the other hand, are often preferred when dealing with tabulated data with
well-defined and limited descriptors [53]. Random forests also have the
ability to yield descriptor importance for interpretability [54] and there has
been significant work on quantification of uncertainties [11, 55, 56].

In this work, we selected random forests as uncertainties are critical to make
use of information acquisition functions within AL and to deal with noisy data
present within our simulations. We use the lolopy library [57] to develop RF
models for melting temperature of MPCAs using composition and derived
descriptors as inputs. The models provide an expectation value and sample-wise
uncertainties that are fed to various acquisition functions to select the next
experiment. Uncertainty estimates are obtained based on the jackknife after
bootstrap and infinitesimal jackknife variance with Monte Carlo sampling
correction, shown in Equation 1, as done in Ling et al. [11] based on the work
of Efron [55] and Wager [56].

$\sigma_{i}^{2}(x)=Cov_{j}[n_{i,j},t_{j}(x)]^{2}+[\bar{t}_{-_{i}}(x)-\bar{t}(x)]^{2}-ev/T$
(1)

where $\sigma_{i}^{2}(x)$ is the variance at point $x$ due to training point
$i$, $n_{i,j}$ is the number of times point $i$ was used to train tree $j$,
$\bar{t}_{-_{i}}(x)$ is the average of the prediction over the trees that were
fit without using point $i$, $\bar{t}(x)$ average of the prediction over all
of the trees,$t_{j}(x)$ prediction of tree index $j$ for point $x$, $e$ is
Euler’s number and $v$ is the variance over all trees. Uncertainty estimates
are computed as sum of contributions from each estimator for each training
point, as seen in Equation 2:

$\sigma[M(x)]=\sqrt{\sum\limits_{i=1}^{S}max[\sigma^{2}_{i}(x),\omega]+\widetilde{\sigma}^{2}(x)}$
(2)

where $\omega$, the noise threshold, which is equivalent to the
$min\sigma^{2}(x_{i})$, the magnitude of the minimum variance, $[M(x_{i})]$ is
the model prediction of point $x_{i}$, $\tilde{\sigma}^{2}(x)$ is an explicit
bias model, and S is the number of training points. This explicit bias model
and uncertainty estimates coincide with those described in Ling et al. [11]

###  2.5 Materials Descriptors and Model Hyperparameters

For the RF models described in subsection 2.4 to predict melting temperature,
it is vital to represent our alloys with descriptors that can serve as inputs.
A brute-force approach would use composition as the sole inputs but, given the
relative scarcity of the data in materials applications, enhancing the inputs
with descriptors or features expected to correlate with the outputs is highly
beneficial. Previous work by Zhang et al. [58] in which they used a
classification algorithm to predict the phases of alloys outlined a materials
descriptor space and a genetic algorithm strategy to down select a subset of
the pool of properties. The complete set included 70 properties based on the
molar average value and mismatch value for each elemental property of the
alloy. Their work showed that RFs outperform other models in the prediction
task, but their algorithm selects subsets limited to four descriptors.
Similarly, work by Zhou et al. [59] shows that machine learning models can
help the exploration of phase design of high entropy alloys. Their model uses
13 descriptors based on the average and standard deviation of properties like
atomic radius, melting temperature, mixing enthalpy, mixing entropy,
electronegativity, bulk modulus and valence electron concentration.

In this work, atomic properties for elements contained in the alloys are
queried from properties available online on Pymatgen[60] and used to generate
a unique fingerprint of each alloy using a weight-average rule. The final set
of descriptors was based on the strongest correlation with melting temperature
of our initial set, measured with Pearson correlation coefficients. This can
be seen in Supplemental Information (SI) as Figure S1. However, it nicely
reflects some of the candidates proposed by Zhang et al. [58] as it includes
properties such as: melting point, electronegativity, boiling temperature,
atomic mass, atomic radii, and density. Another subset of the highly
correlated descriptors relates to mechanical and thermal properties and
includes Young’s modulus, Poisson’s ratio, hardness, and coefficient of
thermal expansion. Finally, electrical resistivity was added to the
descriptors as it showed a high correlation to melting temperature.

To test our initial model, we split the initial set in an 80/20 train/test
split. In Figure 5 we present a parity plot of the melting temperatures
predicted by our RF model and the melting temperatures obtained via MD
simulations.

Figure 5: Comparison for MD simulated and RF predicted melting temperatures
for the 39 MPCAs compositions in our initial set.

Given the small initial dataset, large deviations in MAE can occur due to the
random selection of the testing set. To properly determine the model’s
accuracy, the dataset was shuffled, split, and rerun 30 times through the
random forest. The average of the MAE for the RF predicted melting temperature
with respect to the MD simulated temperature was 36.44 K with an uncertainty
estimate of 12.2 K. A 10-fold cross validation analysis was used to create a
normalized residual distribution to determine calibration of the uncertainty.
Models with well calibrated uncertainties would produce a Gaussian
distribution with a mean of zero and unit standard deviation as seen in Figure
6.

Figure 6: Probability densities of normalized residuals of RF model computed
through tenfold cross-validation. Solid line is perfectly calibrated
uncertainties.

We optimized the number of estimators (trees) to 350 trees at the maximum
depth of each estimator, using the mean absolute error (MAE) as a metric to
measure the accuracy of our model predictions. Increasing the number of trees
reduces variance, while increasing the maximum depth can help reduce bias. Our
optimization stopped based on the idea that the larger quantity of either
would improve the RF until we observed diminishing returns. Plots for the
optimization of these parameters are included in the supplemental information
as Figure S2.

###  2.6 Acquisition functions

Acquisition functions are mathematical expressions designed to select which of
the possible tests (MD simulations in our case) to carry out next. This is
done by considering the model predicted value and uncertainty estimate of the
quantity of interest for all possible tests. For our case, we use the
predicted value and uncertainty estimates provided by the RF model. Various
acquisition functions have been proposed with different balances between
exploitation and exploration. Exploration functions select areas of high
uncertainty and, on the opposite end, functions based purely on exploitation
select candidates with the highest expected values.

A purely exploitative approach would be maximizing the mean predicted value;
this maximum mean (MM) function can be written as:

$MM:x^{*}=argmax\;\;E[M(x_{i})]$ (3)

Such functions can easily get trapped in local maxima, and some degree of
exploration is often desirable. Upper confidence bound (UCB) queries the
sample with the maximum predicted value plus its uncertainty estimate. For
this study, the adjustable parameter K is set to unity (K = 1).

$UCB:x^{*}=argmax\;\;(E[M(x_{i})]+K*\sigma[M(x_{i})])$ (4)

Maximum likelihood of improvement (MLI) chooses the sample with the highest
probability of surpassing the current best previously evaluated material. This
probability can be computed as the Z-score for predictions that are normally
distributed. Therefore, it is represented by the difference between the
expected value of the prediction and the value of the current best case,
$x_{best}$, over the uncertainty estimate.

$MLI:x^{*}=argmax\;\;\frac{E[M(x_{i})]-E[M(x_{best})]}{\sigma[M(x_{i})]}$ (5)

Maximum expected improvement (MEI) [61] works by modeling our knowledge of the
prediction as a normal distribution. It uses the model’s mean prediction and
uncertainty estimate to draw a normal probability density function at each
point. The improvement can then be understood as the probability that the
function at this value surpasses our current maximum. Intuitively, this
relates to the tail of the distribution crossing our current maximum. In the
absence of uncertainties, the function maximum expected improvement (MEI) [61,
62], reduces to maximizing the expected values. We note that the MM function
was incorrectly labeled as MEI in Refs. [12, 11].

$\begin{split}MEI:x^{*}=argmax\;\;\rho\;(E[M(x_{i})]-E[M(x_{best})],\sigma[M(x_{i})])\\\
with\;\;\rho\;(z,s)=\begin{cases}s\phi^{\prime}(\frac{z}{s})+z\phi(\frac{z}{s})&s>0\\\
max(z,0)&s=0\end{cases}\end{split}$ (6)

Finally, maximum uncertainty (MU) is at the extreme end of exploration and
only focuses on the candidates with the highest uncertainty in their
predictions. MU has little incentive to find the top performer and finds its
best use when trying to fine-tune ML models for surrogate-based optimization.

$MU:x^{*}=argmax\;\;\sigma[M(x_{i})]$ (7)

In these Equations, $x^{*}$ marks the composition selected by the function to
be tested next. $x_{i}$ are the possible experiments to run (our search
space). $x_{best}$ is the current best performer in the training set.
$E[M(x_{i})]$ is the expectation value of the prediction at point $x$. This
expectation value is equal to $\frac{1}{N}\sum\limits_{j}^{n_{T}}t_{j}(x)$
where $t_{j}(x)$ is the prediction of tree $j$ for point $x$ and $n_{T}$ is
the total number of trees. Within MEI, $\phi$ represents the Gaussian
cumulative distribution function, and $\phi^{\prime}$ is the Gaussian
probability density function. Finally, “arg max” (arguments of the maxima)
operation returns the sample material where the function is maximized for the
quantity of interest.

##  3 Active learning: exploring high melting temperature alloys

As described above, we started with knowledge of the melting temperature of 39
alloys and use AL to search for the alloy with the highest melting temperature
in a 5-dimensional compositional space with 554 possible alloys. The
performance of these acquisition functions was assessed from AL workflows with
a budget of 40 experiments, mimicking time and resource constraints in real
world applications. Here an experiment is equal to the convergence of a
composition of an alloy which may require several simulations. As previously
described, convergence is defined when the simulation outputs an alloy that
reached temperature equilibrium and the fraction of atoms in each of the
phases needs to account for 35% to 65% of the system. At each step of the
iterative process, we use the RF models to evaluate each acquisition function
over all unexplored alloys and select the alloy that maximizes the selected
function. A new MD simulation is launched with that alloy and the resulting
data is added to the data set and the processes is restarted.

The MD simulation is run for the desired time and convergence, or lack
thereof, is determined with the criteria in Section 2.2. Simulation inputs are
box length n=18, MD simulation time, and temperatures of the liquid and solid
regions. Additional inputs include seeds for pseudo-random number generation
for initialization of the initial random structure and atomic velocities. The
initial temperatures for the liquid and solid regions from the melting
temperature predicted by the RF (Tm) as follows: Tliquid,t=0 = 1.25 Tm ;
Tsolid,t=0 = 0.5 Tm. If a converged melting temperature calculation is not
reached, the initial temperatures are adjusted, and a new simulation launched.
The initial temperatures are increased or decreased by 5% if the resulting
structural analysis indicated little liquid or solid left, respectively. This
process continues until convergence is achieved. Once convergence is reached
the value will be added into the known set, the model will make new
predictions, and the process will repeat until our design goal is reached or
our experiment budget is exhausted.

###  3.1 Active Learning results

To explore the combination of AL with MD simulations we tested six selection
strategies (five information acquisition functions described in Section 2.6
and a random sampling baseline). To assess the effect the noise in the data we
considered three simulation times 50, 100 and 200 ps. Each of the 18 runs used
the same initial set as described in section 2.1.1.

Figures 7, 8 and 9 show the melting temperature, both predicted by the RF
(filled symbols) and the one obtained from the subsequent MD simulation (open
symbols), as a function of experiment number for the various acquisition
functions and simulation times. The insets in Figures 7, 8 and 9, show the
composition of the alloys with the highest predicted melting temperatures.
Tables 2, 3 and 4 summarize the various AL outcomes: for the alloy with the
highest melting temperature within the budget of 40 possible compositions for
each case, we list the temperature predicted by the RF and the actual value
obtained from the MD simulations, the uncertainty in the RF model, number of
experiments required to locate the alloy with the highest melting temperature
and its composition.

Figure 7: Performance of different acquisition functions in a 40-experiment
budget AL run with 50 picoseconds MD simulation time. All functions start from
identical initial sets. Open symbols represent MD simulated melting
temperatures. The filled symbols with error bars represent RF predicted
melting temperatures. Xx represents other elements, Cu and Cr. The insert
includes a close-up on the best performing MPCAs compositions. These
compositions contain high quantities of Fe. Figure 8: Performance of different
acquisition functions in a 40-experiment budget AL run with 100 picoseconds MD
simulation time. All functions start from identical initial sets. Open symbols
represent MD simulated melting temperatures. The filled symbols with error
bars represent RF predicted melting temperatures. Xx represents other
elements, Cu and Cr. The insert includes a close-up on the best performing
MPCAs compositions. These compositions contain high quantities of Fe. Figure
9: Performance of different acquisition functions in a 40-experiment budget AL
run with 200 picoseconds MD simulation time. All functions start from
identical initial sets. Open symbols represent MD simulated melting
temperatures. The filled symbols with error bars represent RF predicted
melting temperatures. Xx represents other elements, Cu and Cr. The insert
includes a close-up on the best performing MPCAs compositions. These
compositions contain high quantities of Fe. Table 2: Best results from
information acquisition functions running simulations for 100 picoseconds. RF
prediction melting temperature and uncertainty taken when best composition
selected for each acquisition function.

Acquisition function |  MD simulated Tm (K) |  RF predicted Tm (K) |  RF uncertainty (K) | Experiment | Composition
---|---|---|---|---|---
MM | 2472 | 2345 | 71 | 31 | Fe 50%, Ni 50%
UCB | 2498 | 2416 | 63 | 10 | Co 10%, Fe 50%, Ni 40%
MLI | 2502 | 2348 | 98 | 12 | Fe 50%, Ni 50%
MEI | 2492 | 2403 | 82 | 12 | Co 20%, Fe 50%, Ni 30%
MU | 2282 | 2137 | 118 | - | Co 40%, Fe 20%, Ni 40%
RAND | 2382 | 2202 | 52 | - | Cu 10%, Fe 50%, Ni 40%

Table 3: Best results from information acquisition functions running
simulations for 100 picoseconds. RF prediction melting temperature and
uncertainty taken when composition selected selected for each acquisition
function.

Acquisition function |  MD simulated Tm (K) |  RF predicted Tm (K) |  RF uncertainty (K) | Experiment | Composition
---|---|---|---|---|---
MM | 2489 | 2460 | 36 | 24 | Fe 50%, Ni 50%
UCB | 2521 | 2422 | 78 | 15 | Fe 50%, Ni 50%
MLI | 2539 | 2455 | 60 | 15 | Fe 50%, Ni 50%
MEI | 2486 | 2415 | 87 | 12 | Co 30%, Fe 50%, Ni 20%
MU | 2299 | 2160 | 117 | - | Co 20%, Cu 20%, Fe 50%, Ni 10%
RAND | 2382 | 2202 | 52 | - | Cu 10%, Fe 50%, Ni 40%

Table 4: Best results from information acquisition functions running
simulations for 200 picoseconds. RF prediction melting temperature and
uncertainty taken when composition selected selected for each acquisition
function.

Acquisition function |  MD simulated Tm (K) |  RF predicted Tm (K) |  RF uncertainty (K) | Experiment | Composition
---|---|---|---|---|---
MM | 2485 | 2461 | 29 | 27 | Fe 50%, Ni 50%
UCB | 2487 | 2417 | 81 | 17 | Co 30%, Fe 50%, Ni 20%
MLI | 2480 | 2413 | 85 | 14 | Co 10%, Fe 50%, Ni 40%
MEI | 2476 | 2418 | 91 | 9 | Co 10%, Fe 50%, Ni 50%
MU | 2285 | 2160 | 118 | - | Co 40%, Cu 10%, Fe 30%, Ni 20%
RAND | 2432 | 2360 | 80 | - | Co 40%, Fe 40%, Ni 20%

Quite interestingly, all acquisition functions that take exploitation into
consideration (MLI, UCB, MEI, and MM) find the alloys with the highest melting
temperatures in few experiments, regardless of the simulation time used. This
is quite remarkable given the significant variability in output for the 50 ps
long simulations. We find that highest melting temperature alloys consist of
50 at.% Fe with the remaining 50 at.% distributed between Ni and Co; this set
of alloys are predicted to have nearly identical melting temperatures, see
Figure 10. As expected, the MU function and random exploration do not explore
high-melting temperature compositions within the 40 experiment budget. The
graph seen for MM, UCB, MLI, and MEI has been made for MU and random and are
in the SI as S3, S4.

Figure 10: Predicted melting temperature best performing alloys with
composition Fe50CoxNi50-x. Mean values and uncertainty estimates over 36
independent runs are shown.

In all our tests, we find that the acquisition functions that combine both
exploration and exploitation (UCB, MLI, MEI) consistently outperform the
purely exploitative function (MM), which can spend resources in local maxima.
As seen in Tables 2-4 UCB, MLI, and MEI find the optimal composition within
less than 17 experiments while MM found it past 27. While in our example MM
finds the family of best performer candidates, it takes approximately twice as
many experiments as the other functions.

###  3.2 Noisy data and RF uncertainty

We find that the AL workflow can find best performer alloys in approximately
the same number of iterations regardless of the simulation time use. Between 9
and 17 simulations were required for the UCB, MEI, and MLI strategies. As
mentioned above, this is despite the significant level of noise, especially in
the 50 ps runs. We also find that while longer MD simulations do not result in
a reduction in the number of experiments required to achieve the design goal,
it results in improved accuracy of the RF prediction. Tables 2, 3 and 4 show
that the workflow with 50 ps runs significantly underestimates the temperature
of the optimal alloy and the uncertainty estimates are overly optimistic.

To better understand how the models gain knowledge in the presence of noise we
explore the evolution of the model predictions for Fe50CoxNi50-x alloys, the
top performers, at various stages during the AL process. Figure 11 shows mean
and uncertainty estimate predictions at four stages of an MLI workflow for the
50 ps case; results for 100 and 200 ps simulations show similar trends and are
included in the SI as Figures S5, S6. Open circles denote MD simulations not
already explored at the corresponding cycle and filled circles represent the
same values for compositions that have been explored. We note that noise in
the data affects both mean predictions and uncertainty estimates and can thus
have non-trivial effects on AL decisions. Figure 11 shows that even with only
the initial 39 datapoints (Iteration 0), the model predicts relatively high
melting temperature for these alloys, comparable to the highest melting
temperature materials in the initial set, see Figure 7. As the AL workflow
with MLI explore compositions, the mean prediction for the selected family
improves from 11 (a) to (b) as the mean moves towards the true MD temperatures
of the alloys. At step 12, Fig. 11 (d), the model predicts high melting
temperatures for the alloys with x=0, 10, 20 % and x=0 is selected for
simulation. In subsequent steps, the MLI algorithm hones in the rest of the
family and improves the model accuracy. When comparing the results for 50 ps
runs with those of 100 and 200 ps (SI Figures S5 and S6) we observe similar
paths to the optimal alloys with a decrease in uncertainties for longer
simulation times. This is because the variability in the MD simulations is
relatively small compared to the differences in melting temperatures across
alloys.

Figure 11: Red lines represent the predicted mean melting temperature and
shaded region represents the uncertainty estimates for predictions on
compositions as a function of Co content in Fe50CoxNi50-x at various stages of
the AL workflow. Whereby x starts at 0% a goes to 50% with 1% step size.
Results correspond to the MLI function with a 50 ps MD simulation time. Open
circles represent MD-simulated values unknown to the model at the time and
filled symbols represent the values included in the model.

##  4 Conclusions

The combination of machine learning and physics-based simulations holds great
promise to guide the experimental search and accelerate the development of
novel materials. The use of surrogate-based models that can make use of
computational simulations instead of experiments can save resources and time
in design cycles. We find the molecular dynamics simulations, even under
significant stochastic noise, can be a powerful tool to be used in active
learning efforts.

We developed an AL workflow with MD simulations for the discovery of MPCAs
materials for high-temperature applications. Starting from a limited initial
set, we created a set of descriptors and a RF model to use as a surrogate for
the optimization of our AL model. We determined that RF models are capable of
handling multi-dimensional spaces with a limited initial dataset. We paired
our RF model with a cloud-based simulation environment for automatic querying
of melting temperatures from MD simulations with an interatomic potential for
FCC alloys. We analyzed the effect of MD simulation time in the overall
uncertainty quantification when running such calculations and showed that it
has an important effect in capturing the sample-to-sample variability caused
by the stochastic nature of MD. We compared the performance of four different
acquisition functions in a closed, but high-dimensional, search space. We
found that the best performing acquisition functions, MLI, MEI, and UCB, are
combinations of exploitation and exploration.

The workflow we developed can be adapted toward the optimization of other
material properties and to the use of other sources of information other than
computational simulations. We provide the analysis and all other relevant
information in a public SimTool on nanoHUB.org.

##  5 Models and Data Availability

The entire simulation workflow to calculate melting temperature, including
pre- and post-analysis was implemented as a SimTool [8] in nanoHUB [63] and is
available for online simulation [64]. SimTools [8] are complete simulation
workflows with verified inputs and outputs that enable workflows and the data
they generate to be findable, accessible, and reproducible. An additional
advantage of using SimTools is that every simulation performed in nanoHUB is
cached, so that given identical input parameters, cached results are retrieved
rather than re-running a simulation. The stand-alone SimTool [64] and the AL
workflow are available for online simulations in nanoHUB [65].

## Acknowledgements

This effort was supported by the US National Science Foundation
(DMREF-1922316). We acknowledge computational resources from nanoHUB and
Purdue University through the Network for Computational Nanotechnology. J. C.
V. thanks the Science and Technology Council of México (Consejo Nacional de
Ciencia y Tecnología, CONACYT) for partial financial support of this research.

## Ethics declarations

### Conflict of interest

On behalf of all authors, the corresponding author states that there is no
conflict of interest.

### Author Information

Corresponding Author

Alejandro Strachan \- School of Materials Engineering and Birck Nanotechnology
Center, Purdue University, West Lafayette, Indiana 47907; Email:
<EMAIL_ADDRESS>

## References

  * [1] Natalie S Eyke, William H Green, and Klavs F Jensen.  Iterative experimental design based on active machine learning reduces the experimental burden associated with reaction screening.  Reaction Chemistry & Engineering, 5(10):1963–1972, 2020.
  * [2] Lifei Zhao, Zhen Li, Bruce Caswell, Jie Ouyang, and George Em Karniadakis.  Active learning of constitutive relation from mesoscopic dynamics for macroscopic modeling of non-newtonian flows.  Journal of Computational Physics, 363:116–127, 2018.
  * [3] A Gilad Kusne, Heshan Yu, Changming Wu, Huairuo Zhang, Jason Hattrick-Simpers, Brian DeCost, Suchismita Sarker, Corey Oses, Cormac Toher, Stefano Curtarolo, et al.  On-the-fly closed-loop materials discovery via bayesian active learning.  Nature communications, 11(1):1–11, 2020.
  * [4] Dezhen Xue, Prasanna V Balachandran, John Hogden, James Theiler, Deqing Xue, and Turab Lookman.  Accelerated search for materials with targeted properties by adaptive design.  Nature communications, 7(1):1–9, 2016.
  * [5] Kevin Tran and Zachary W Ulissi.  Active learning across intermetallics to guide discovery of electrocatalysts for co2 reduction and h2 evolution.  Nature Catalysis, 1(9):696–703, 2018.
  * [6] Anubhav Jain, Shyue Ping Ong, Geoffroy Hautier, Wei Chen, William Davidson Richards, Stephen Dacek, Shreyas Cholia, Dan Gunter, David Skinner, Gerbrand Ceder, et al.  Commentary: The materials project: A materials genome approach to accelerating materials innovation.  APL materials, 1(1):011002, 2013.
  * [7] Atsuto Seko, Atsushi Togo, Hiroyuki Hayashi, Koji Tsuda, Laurent Chaput, and Isao Tanaka.  Prediction of low-thermal-conductivity compounds with first-principles anharmonic lattice-dynamics calculations and bayesian optimization.  Physical review letters, 115(20):205901, 2015.
  * [8] Saaketh Desai, Steven Clark, and Alejandro Strachan.  Introduction to simtools, Apr 2020.  URL: https://nanohub.org/resources/introtosimtools, DOI: 10.21981/4Z7T-X415.
  * [9] Prasanna V Balachandran, Dezhen Xue, James Theiler, John Hogden, and Turab Lookman.  Adaptive strategies for materials design using uncertainties.  Scientific reports, 6(1):1–9, 2016.
  * [10] Chiho Kim, Anand Chandrasekaran, Anurag Jha, and Rampi Ramprasad.  Active-learning and materials design: The example of high glass transition temperature polymers.  MRS Communications, 9(3):860–866, 2019.
  * [11] Julia Ling, Maxwell Hutchinson, Erin Antono, Sean Paradiso, and Bryce Meredig.  High-dimensional materials and process optimization using data-driven experimental design with well-calibrated uncertainty estimates.  Integrating Materials and Manufacturing Innovation, 6(3):207–217, 2017.
  * [12] Juan C Verduzco, Ernesto E Marinero, and Alejandro Strachan.  An active learning approach for the design of doped LLZO ceramic garnets for battery applications.  Integrating Materials and Manufacturing Innovation, pages 1–12, 2021.
  * [13] National Science and Technology Council (US).  Materials genome initiative for global competitiveness.  Executive Office of the President, National Science and Technology Council, 2011.
  * [14] John Wise, Alexandra Grebe de Barron, Andrea Splendiani, Beeta Balali-Mood, Drashtti Vasant, Eric Little, Gaspare Mellino, Ian Harrow, Ian Smith, Jan Taubert, et al.  Implementation and relevance of fair data principles in biopharmaceutical r&d.  Drug discovery today, 24(4):933–938, 2019.
  * [15] Burr Settles.  Synthesis Lectures on Artificial Intelligence and Machine Learning: Active Learning.  Morgan & Claypool Publishers, San Rafael, United States, 2011.
  * [16] Yanyao Shen, Hyokun Yun, Zachary C Lipton, Yakov Kronrod, and Animashree Anandkumar.  Deep active learning for named entity recognition.  arXiv preprint arXiv:1707.05928, 2017.
  * [17] Pilsun Yoo, Michael Sakano, Saaketh Desai, Md Mahbubul Islam, Peilin Liao, and Alejandro Strachan.  Neural network reactive force field for C, H, N, and O systems.  npj Computational Materials, 7(1):1–10, 2021.
  * [18] Aaron Gilad Kusne, Tieren Gao, Apurva Mehta, Liqin Ke, Manh Cuong Nguyen, Kai-Ming Ho, Vladimir Antropov, Cai-Zhuang Wang, Matthew J Kramer, Christian Long, et al.  On-the-fly machine-learning for high-throughput experiments: search for rare-earth-free permanent magnets.  Scientific reports, 4(1):1–7, 2014.
  * [19] Pavel Nikolaev, Daylond Hooper, Frederick Webber, Rahul Rao, Kevin Decker, Michael Krein, Jason Poleski, Rick Barto, and Benji Maruyama.  Autonomy in materials research: a case study in carbon nanotube growth.  npj Computational Materials, 2(1):1–6, 2016.
  * [20] Fang Ren, Logan Ward, Travis Williams, Kevin J Laws, Christopher Wolverton, Jason Hattrick-Simpers, and Apurva Mehta.  Accelerated discovery of metallic glasses through iteration of machine learning and high-throughput experiments.  Science advances, 4(4):eaaq1566, 2018.
  * [21] Marcus M Noack, Kevin G Yager, Masafumi Fukuto, Gregory S Doerk, Ruipeng Li, and James A Sethian.  A kriging-based approach to autonomous experimentation with applications to x-ray scattering.  Scientific reports, 9(1):1–19, 2019.
  * [22] Shin Kiyohara, Hiromi Oda, Koji Tsuda, and Teruyasu Mizoguchi.  Acceleration of stable interface structure searching using a kriging approach.  Japanese Journal of Applied Physics, 55(4):045502, 2016.
  * [23] Erik Van Der Giessen, Peter A Schultz, Nicolas Bertin, Vasily V Bulatov, Wei Cai, Gábor Csányi, Stephen M Foiles, Marc GD Geers, Carlos González, Markus Hütter, et al.  Roadmap on multiscale materials modeling.  Modelling and Simulation in Materials Science and Engineering, 28(4):043001, 2020.
  * [24] Alberto M Cuitiño, Laurent Stainier, Guofeng Wang, Alejandro Strachan, Tahir Çağin, William A Goddard, and Michael Ortiz.  A multiscale approach for modeling crystalline solids.  Journal of computer-aided materials design, 8(2):127–149, 2001\.
  * [25] Alejandro Strachan, Tahir Çağin, and William A Goddard III.  Phase diagram of mgo from density-functional theory and molecular-dynamics simulations.  Physical Review B, 60(22):15084, 1999.
  * [26] Adri CT Van Duin, Siddharth Dasgupta, Francois Lorant, and William A Goddard.  Reaxff: a reactive force field for hydrocarbons.  The Journal of Physical Chemistry A, 105(41):9396–9409, 2001.
  * [27] Albert P Bartók, Mike C Payne, Risi Kondor, and Gábor Csányi.  Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons.  Physical review letters, 104(13):136403, 2010.
  * [28] Aidan P Thompson, Laura P Swiler, Christian R Trott, Stephen M Foiles, and Garritt J Tucker.  Spectral neighbor analysis method for automated generation of quantum-accurate interatomic potentials.  Journal of Computational Physics, 285:316–330, 2015.
  * [29] Jörg Behler and Michele Parrinello.  Generalized neural-network representation of high-dimensional potential-energy surfaces.  Physical review letters, 98(14):146401, 2007.
  * [30] Thaer M. Dieb, Shenghong Ju, Kazuki Yoshizoe, Zhufeng Hou, Junichiro Shiomi, and Koji Tsuda.  Mdts: automatic complex materials design using monte carlo tree search.  Science and technology of advanced materials, 18(1):498–503, 2017\.
  * [31] Thaer M Dieb, Shenghong Ju, Junichiro Shiomi, and Koji Tsuda.  Monte carlo tree search for materials design and discovery.  MRS Communications, 9(2):532–536, 2019.
  * [32] Tarak K Patra, Troy D Loeffler, and Subramanian KRS Sankaranarayanan.  Accelerating copolymer inverse design using monte carlo tree search.  Nanoscale, 12(46):23653–23662, 2020.
  * [33] Troy D Loeffler, Suvo Banik, Tarak K Patra, Michael Sternberg, and Subramanian KRS Sankaranarayanan.  Reinforcement learning in discrete action space applied to inverse defect design.  Journal of Physics Communications, 5(3):031001, 2021.
  * [34] Stéphane Gorsse, Jean-Philippe Couzinié, and Daniel B Miracle.  From high-entropy alloys to complex concentrated alloys.  Comptes Rendus Physique, 19(8):721–736, 2018.
  * [35] Yiping Lu, Yong Dong, Sheng Guo, Li Jiang, Huijun Kang, Tongmin Wang, Bin Wen, Zhijun Wang, Jinchuan Jie, Zhiqiang Cao, et al.  A promising new class of high-temperature alloys: eutectic high-entropy alloys.  Scientific reports, 4(1):1–5, 2014.
  * [36] ZQ Xu, ZL Ma, M Wang, YW Chen, YD Tan, and XW Cheng.  Design of novel low-density refractory high entropy alloys for high-temperature applications.  Materials Science and Engineering: A, 755:318–322, 2019.
  * [37] Sathiyamoorthi Praveen and Hyoung Seop Kim.  High-entropy alloys: Potential candidates for high-temperature applications–an overview.  Advanced Engineering Materials, 20(1):1700645, 2018.
  * [38] O.N. Senkov, G.B. Wilks, D.B. Miracle, C.P. Chuang, and P.K. Liaw.  Refractory high-entropy alloys.  Intermetallics, 18(9):1758–1765, 2010.
  * [39] WL Wang, L Hu, SB Luo, LJ Meng, DL Geng, and B Wei.  Liquid phase separation and rapid dendritic growth of high-entropy CoCrCuFeNi alloy.  Intermetallics, 77:41–45, 2016.
  * [40] Feng He, Zhijun Wang, Qingfeng Wu, Junjie Li, Jincheng Wang, and CT Liu.  Phase separation of metastable CoCrFeNi high entropy alloy at intermediate temperatures.  Scripta Materialia, 126:15–19, 2017.
  * [41] Qi-Jun Hong, Jan Schroers, Douglas Hofmann, Stefano Curtarolo, Mark Asta, and Axel van de Walle.  Theoretical prediction of high melting temperature for a Mo–Ru–Ta–W HCP multiprincipal element alloy.  npj Computational Materials, 7(1):1–4, 2021.
  * [42] James R Morris, CZ Wang, KM Ho, and CT Chan.  Melting line of aluminum from simulations of coexisting phases.  Physical Review B, 49(5):3109, 1994.
  * [43] James R Morris and Xueyu Song.  The melting lines of model systems calculated from coexistence simulations.  The Journal of chemical physics, 116(21):9352–9358, 2002.
  * [44] Diana Farkas and Alfredo Caro.  Model interatomic potentials and lattice strain in a high-entropy alloy.  Journal of Materials Research, 33(19):3218–3225, 2018.
  * [45] Sheng-Nian Luo, Thomas J Ahrens, Tahir Çağın, Alejandro Strachan, William A Goddard III, and Damian C Swift.  Maximum superheating and undercooling: Systematics, molecular dynamics simulations, and dynamic experiments.  Physical Review B, 68(13):134206, 2003.
  * [46] A. S. Lemak and N. K. Balabaev.  On The Berendsen Thermostat.  Molecular Simulation, 13(3):177–187, 1994.
  * [47] Peter J. Daivis and Denis J. Evans.  Comparison of constant pressure and constant volume nonequilibrium simulations of sheared model decane.  J. Chem. Phys, 100, 1994.
  * [48] Peter Mahler Larsen, Søren Schmidt, and Jakob Schiøtz.  Robust structural identification via polyhedral template matching.  Modelling and Simulation in Materials Science and Engineering, 24(5):055007, 2016.
  * [49] Alexander Stukowski.  Visualization and analysis of atomistic simulation data with OVITO–the open visualization tool.  Modelling and Simulation in Materials Science and Engineering, 18(1):015012, 2009.
  * [50] Z. Wu, H. Bei, F. Otto, G.M. Pharr, and E.P. George.  Recovery, recrystallization, grain growth and phase stability of a family of fcc-structured multi-component equiatomic solid solution alloys.  Intermetallics, 46:131–140, 2014.
  * [51] Juan J Valencia and Peter N Quested.  Thermophysical properties.  NIST, 2013.  Table 5: Thermal expansion of selected pure metals at temperatures close to melting.
  * [52] Rohit K Tripathy and Ilias Bilionis.  Deep UQ: Learning deep neural network surrogate models for high dimensional uncertainty quantification.  Journal of computational physics, 375:565–588, 2018.
  * [53] Tianqi Chen and Carlos Guestrin.  Xgboost: A scalable tree boosting system.  In Proceedings of the 22nd ACM sigkdd international conference on knowledge discovery and data mining, pages 785–794, 2016.
  * [54] Leo Breiman.  Random forests.  Machine learning, 45(1):5–32, 2001.
  * [55] Bradley Efron.  Model Selection Estimation and Bootstrap Smoothing.  Division of Biostatistics, Stanford University, 2012.
  * [56] Stefan Wager, Trevor Hastie, and Bradley Efron.  Confidence intervals for random forests: The jackknife and the infinitesimal jackknife.  The Journal of Machine Learning Research, 15(1):1625–1651, 2014\.
  * [57] M Hutchinson.  Citrine informatics: Lolo, 2016.  URL: https://github.com/CitrineInformatics/ lolo.
  * [58] Yan Zhang, Cheng Wen, Changxin Wang, Stoichko Antonov, Dezhen Xue, Yang Bai, and Yanjing Su.  Phase prediction in high entropy alloys with a rational selection of materials descriptors and machine learning models.  Acta Materialia, 185:528–539, 2020.
  * [59] Ziqing Zhou, Yeju Zhou, Quanfeng He, Zhaoyi Ding, Fucheng Li, and Yong Yang.  Machine learning guided appraisal and exploration of phase design for high entropy alloys.  npj Computational Materials, 5(1):1–9, 2019.
  * [60] Shyue Ping Ong, William Davidson Richards, Anubhav Jain, Geoffroy Hautier, Michael Kocher, Shreyas Cholia, Dan Gunter, Vincent L Chevrier, Kristin A Persson, and Gerbrand Ceder.  Python materials genomics (pymatgen): A robust, open-source python library for materials analysis.  Computational Materials Science, 68:314–319, 2013.
  * [61] Donald R Jones, Matthias Schonlau, and William J Welch.  Efficient global optimization of expensive black-box functions.  Journal of Global optimization, 13(4):455–492, 1998.
  * [62] Emmanuel Vazquez and Julien Bect.  Convergence properties of the expected improvement algorithm with fixed mean and covariance functions.  Journal of Statistical Planning and inference, 140(11):3088–3095, 2010.
  * [63] Alejandro Strachan, Gerhard Klimeck, and Mark Lundstrom.  Cyber-enabled simulations in nanoscale science and engineering.  Computing in Science & Engineering, 12(2):12–17, 2010.
  * [64] Zachary D McClure, Saaketh Desai, and Alejandro Strachan.  High entropy alloy melting point calculation, Mar 2020.  URL: https://nanohub.org/resources/meltheas, DOI: 10.21981/W5VD-T039.
  * [65] Juan Carlos Verduzco Gastelum, David Enrique Farache, Zachary D McClure, Saaketh Desai, and Alejandro Strachan.  Active learning workflow for mpcas, Oct 2021.  URL: https://nanohub.org/resources/activemeltheas, DOI: 10.21981/NK7E-HA16.


# A Linear Bayesian Learning Receiver Scheme for Massive MIMO Systems
††thanks: An extended version of this paper including all proofs is in
preparation [1]. This research was supported by an Australian Research Council
Discovery Early Career Research Award (DE150101704) and the Research Training
Program Stipend from The University of Sydney.

Alva Kosasih1, Wibowo Hardjawana1, Branka Vucetic1, Chao-Kai Wen2
1Centre of Excellence in Telecommunications, The University of Sydney, Sydney,
Australia.
2Institute of Communications Engineering, National Sun Yat-sen University,
Kaohsiung, Taiwan.
Email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Much stringent reliability and processing latency requirements in ultra-
reliable-low-latency-communication (URLLC) traffic make the design of linear
massive multiple-input-multiple-output (M-MIMO) receivers becomes very
challenging. Recently, Bayesian concept has been used to increase the
detection reliability in minimum-mean-square-error (MMSE) linear receivers.
However, the latency processing time is a major concern due to the exponential
complexity of matrix inversion operations in MMSE schemes. This paper proposes
an iterative M-MIMO receiver that is developed by using a Bayesian concept and
a parallel interference cancellation (PIC) scheme, referred to as a linear
Bayesian learning (LBL) receiver. PIC has a linear complexity as it uses a
combination of maximum ratio combining (MRC) and decision statistic combining
(DSC) schemes to avoid matrix inversion operations. Simulation results show
that the bit-error-rate (BER) and latency processing performances of the
proposed receiver outperform the ones of MMSE and best Bayesian-based
receivers by minimum $2$ dB and $19$ times for various M-MIMO system
configurations.

###### Index Terms:

Massive MIMO, PIC, DSC, Bayesian learning, low complexity, URLLC.

##  I Introduction

Massive multiple-input-multiple-output (M-MIMO) technology has been proposed
to support ultra reliability and low latency (URLLC) data transmissions. The
reduction of minimum transmission time intervals (TTI) and bit-error-rate
(before coding) requirements from $15$ ms and $10^{-3}$ to $1$ ms and
$10^{-5}$ in 4G and 5G NR URLLC [2, 3, 4, 5, 6, 7], respectively, pose a
challenging problem in developing an M-MIMO receiver. Two types of M-MIMO
receivers have been investigated in the literature; classical and Bayesian
learning receivers.

Classical receivers can be divided into two categories, non-linear and linear
receivers. While non-linear receivers e.g. maximum likelihood (ML) receiver
[8] can achieve an optimal symbol detection reliability, it suffers a high
computational complexity due to an exhaustive search operation to find the ML
combinations of user symbols. This leads to a long detection processing time
and thus a high latency. The second type is the linear M-MIMO receivers; 1)
minimum mean square error (MMSE) and zero forcing (ZF) receivers [9]; and 2)
iterative parallel interference cancellation (PIC) receivers that use a
combination of maximum ratio combining (MRC) and decision statistic combining
(DSC) schemes [10]. MMSE and ZF receivers both rely on a matrix inversion
operation to cancel multiple user interference which leads to an exponential
increase in computational complexity with the number of antennas. In contrast,
PIC receivers avoid the use of a matrix inversion operation. Specifically, PIC
receivers iteratively performs MRC to estimate multiple user symbols used to
reconstruct interfering symbols. The interfering symbols are then subtracted
from the received signal to recover the desired symbols, in a parallel manner.
Consequently, PIC receivers have a linear computational complexity and thus a
much lower latency. Note that all above linear receivers exhibit poor
detection reliability when compared to an ML receiver.

Recently, Bayesian learning concept has been introduced to reduce the
reliability performance gap between ML and linear receivers. These receivers
are referred to as the Bayesian learning receivers [11, 12, 13, 14, 15, 16].
The Bayesian learning concept is used to incorporate detection probability
measures when estimating the detected symbols from the received signals[17].
The best Bayesian learning receivers in term of bit-error-rate (BER)
performance combines the Bayesian learning concept with MMSE scheme, referred
to as expectation propagation (EP) receivers [11, 12, 13, 14]. Despite a
significant performance improvement compared to the MMSE receivers, matrix
inversion operations are still required. This results in a exponentially
latency processing at the receivers. To avoid the matrix inversion operation,
another Bayesian based receiver referred to as approximate message passing
(AMP) receiver [15, 16] has been proposed, albeit at the cost of higher BER as
compared to EP receivers. Note also that all Bayesian M-MIMO receivers above
rely on the learning parameters that needs to be searched by using a trial and
error process for different environments.

We propose a novel iterative M-MIMO receiver referred to as linear Bayesian
learning (LBL) receiver to cater for higher reliability and lower latency
requirements in URLLC traffic. We first build the system model wherein an
M-MIMO receiver at the base station is used to detect the symbols sent by
multiple users. We assume that the channel estimates for different users are
available and have been calculated by using different processes at the base
stations. The developed LBL receiver consists of three modules; Bayesian
symbol observation (BSO), Bayesian symbol estimate (BSE) and DSC modules. In
the BSO module, we apply the MRC scheme to the received signals in order to
get the symbol estimates, referred to as the observed symbols. For each
observed symbol, PIC scheme is then used to remove its interference. The
symbols variance is also calculated. In the BSE module, the observed symbols
and their variances are then utilised to construct maximum likelihood Gaussian
distribution functions. The soft symbol estimates are then calculated based on
the likelihood functions. These estimates are used to compute symbol error
between estimations and observations. In the DSC module, the value of symbol
errors in current and previous iterations are used to calculate the symbol
estimates. The process is then repeated iteratively and finally the DSC
outputs are taken as the symbol estimates. The simulation results show that
the BER and latency performances of the proposed LBL receiver outperform
existing classical and Bayesian M-MIMO receivers in [9, 10, 15] by minimum $2$
dB and $19$ times for various M-MIMO system configurations. The main
contributions of this paper are as follows

  * •

First Bayesian M-MIMO iterative receiver that uses PIC scheme [10]. This leads
to an elimination of the matrix inversion operations or approximations used in
several advanced iterative receivers [11, 18, 13, 14]. As a result, linearly
computational latency processing and a near-optimal performance are achieved.

  * •

First Bayesian M-MIMO iterative receiver that derives the learning parameters
directly from the symbol errors between estimations and observations in
different iterations via the DSC scheme. This is in contrast to the trial and
error process used by many Bayesian receivers [11, 13, 16, 14] to find the
learning parameters.

  * •

We perform an analysis on the maximum number of users that can be supported by
the proposed receiver in exchange of a low complexity signal processing. This
analysis has not been discussed in the literature.

Notations—$\mathbf{I}$ denotes a proper size identity matrix. For any matrix
$\mathbf{A}$, $\mathbf{A}^{T}$ is the transpose of $\mathbf{A}$,
$\mathbf{A}^{H}$ is the conjugate transpose of $\mathbf{A}$, and ${\sf
tr}(\mathbf{A})$ denotes the traces of $\mathbf{A}$. ${\rm diag}({\bf Q})$
refers to the operation to force the non-diagonal elements of matrix ${\bf Q}$
to zero. Conversely, ${\rm off}({\bf Q})$ refers to the operation to force the
diagonal elements of matrix ${\bf Q}$ to zero. ${{\sf E}}[{\bf x}]$ is the
mean of random vector ${\bf x}$ and ${\sf Var}[{\bf x}]={{\sf
E}}\big{[}\left({\bf x}-{{\sf E}}[{\bf x}]\right)^{2}\big{]}$ is its variance.
${\mathcal{N}}(x_{k},c_{k};v_{k})$ represents a complex single variate
Gaussian distribution of random variable $x_{k}$ with mean $c_{k}$ and
variance $v_{k}$. By letting ${\bf x}=[x_{1},\cdots,x_{K}]^{T},{\bf
c}=[c_{1},\cdots,c_{K}]^{T}$, the multivariate Gaussian distribution of random
vector ${\bf x}$ is denoted as

${\mathcal{N}}({\bf x},{\bf
c};{\boldsymbol{\Sigma}})=\frac{1}{\sqrt{(2\pi)^{K}|{\boldsymbol{\Sigma}}|}}{\rm
exp}\bigg{(}-\frac{1}{2}{\left({{\bf x}-{\bf
c}}\right)^{T}{\boldsymbol{\Sigma}}^{-1}\left({{\bf x}-{\bf
c}}\right)}\bigg{)}$ (1)

${\boldsymbol{\Sigma}}$ is a covariance matrix, and $|{\boldsymbol{\Sigma}}|$
is the determinant of the covariance matrix ${\boldsymbol{\Sigma}}$.

Figure 1: The uplink M-MIMO system. Figure 2:  System model of the LBL
receiver.

##  II System Model

We consider an M-MIMO receiver that receives uplink signals from $K$ users,
each with a single antenna, as depicted in Fig. 1. The receiver is equipped
with a large number of antennas $N>>K$. Each user first maps its information
bit stream to a symbol $x_{k}$ that belongs to a constellation point of
$M$-QAM, $\Omega_{k}=[s_{1},\dots,s_{M}]$ where ${\bf
x}=[x_{1},\cdots,x_{K}]^{T}$, $s_{m}$ is the $m$-th symbol in the $M$-QAM
symbol constellations, and $k=1,\dots,K$. The average symbol energy is
$E_{x}={\sf E}\\{|x_{k}|^{2}\\}$. In the receiver side, the received signal at
the M-MIMO receiver, $\mathbb{y}=[y_{1}\ldots y_{N}]^{T}$ where $y_{n}$ is the
received signal at antenna $n$ can then be written as

${\bf y}={\bf H}{\bf x}+{\bf n}$ (2)

where, ${\bf H}\in\mathbb{C}^{N\times K}$ is the coefficients of complex
Rayleigh wireless channels between $N$ antennas at M-MIMO receiver and $K$
users, ${\bf n}\in\mathbb{C}^{N}$ denotes the additive white Gaussian noise
(AWGN) with zero mean and covariance matrix $\sigma^{2}{\bf I}$.

The function used to cancel multiple user interference and recover the symbols
sent by $K$ users, $\mathbb{x}$ from the received signals, $\mathbb{y}$ when
MMSE receiver is used, is given as

${{\bf x}}\propto\left({\bf H}^{H}{\bf H}+\sigma^{2}\mathbf{I}\right)^{-1}{\bf
H}^{H}{\bf y}.$ (3)

Note that the matrix inversion operation used in (3) is costly as its
complexity increases exponentially with the number of antennas. In contrast,
when the iterative PIC receiver is used there is no matrix inversion operation
performed. At the M-MIMO receiver, the symbols recovery from $\mathbb{y}$ in
iteration $t$, ${\bf x}^{(t)}$ is given as

${\bf x}^{(t)}\propto{\rm diag}^{(-1)}({\bf H}^{H}{\bf H})({\bf H}^{H}{\bf
y}-{\rm off}({\bf H}^{H}{\bf H})\tilde{{\bf x}}^{(t-1)})$ (4)

where $\tilde{{\bf
x}}^{(t-1)}=[\tilde{x}_{1}^{(t-1)},\cdots,\tilde{x}_{K}^{(t-1)}]^{T}$ are the
symbol estimates from previous iteration $t-1$. Note that to estimate ${\bf
x}^{(t+1)}$, we set $\tilde{{\bf x}}^{(t)}={\bf x}^{(t)}$.

##  III Expectation Propagation Receiver

In this section, we briefly review the best Bayesian based receiver in the
literature, Expectation Propagation (EP) receivers [11, 12, 13, 14, 19, 20],
which is a combination of Bayesian and MMSE concepts. The main idea of the EP
is to iteratively approximate the distribution of a random transmitted symbol
vector, ${\bf x}$ by using a Gaussian probability distribution function (PDF)
approximation based on the received signals ${\bf y}$, $p^{(t)}({\bf x}|{\bf
y})$, and a pair of tuning parameters,
$({\boldsymbol{\lambda}}^{(t)},{\boldsymbol{\gamma}}^{(t)})$ which are
obtained from the exponential family distributions [17]. $p^{(t)}({\bf x}|{\bf
y})$ is given as

$\displaystyle p^{(t)}({\bf x}|{\bf y})=\mathcal{N}\left({\bf
x},{{\boldsymbol{\mu}}}^{(t)}_{\rm obs};{\boldsymbol{\Sigma}}^{(t)}_{\rm
obs}\right)$ (5)

with its mean and variance are obtained from the received signal observation,
${\bf y}$,

$\displaystyle{\boldsymbol{\mu}}^{(t)}_{\rm
obs}={{\boldsymbol{\Sigma}}}^{(t)}_{\rm obs}{\left(\sigma^{-2}{\bf
H}^{H}\mathbf{y}+{\boldsymbol{\gamma}}^{(t)}\right)}$ (6a)
$\displaystyle{{\boldsymbol{\Sigma}}}^{(t)}_{\rm obs}={\left(\sigma^{-2}{\bf
H}^{H}{\bf H}+{\boldsymbol{\lambda}}^{-1(t)}\right)}^{-1}.$ (6b)

We then construct a Gaussian likelihood function to approximate the
distribution of received signal ${\bf y}$ for a given statistic of ${\bf x}$
in (6). This is denoted as $p^{(t)}({\bf y}|{\bf x})$,

$\displaystyle p^{(t)}({\bf y}|{\bf x})=\mathcal{N}\left({\bf x},{\bf x}_{\rm
ext}^{(t)};{\bf v}^{(t)}_{{\rm ext}}\right)$ (7)

where $\mathcal{N}\left({\bf x},{\bf x}_{\rm ext}^{(t)};{\bf v}^{(t)}_{{\rm
ext}}\right)=\prod_{k=1}^{K}\mathcal{N}\left(x_{k},x_{k,{\rm
ext}}^{(t)};v_{k,{\rm ext}}^{(t)}\right)$ and

$\displaystyle x_{k,{\rm ext}}^{(t)}=v_{k,{\rm
ext}}^{(t)}{\left(\frac{\mu_{k,{\rm obs}}^{(t)}}{\Sigma_{k,{\rm
obs}}^{(t)}}-\gamma^{(t)}_{k}\right)}$ (8a) $\displaystyle v_{{k,{\rm
ext}}}^{(t)}=\frac{\Sigma_{k,{\rm obs}}^{(t)}}{1-\Sigma_{k,{\rm
obs}}^{(t)}\lambda_{k}^{(t)}}.$ (8b)

Here, $x_{k,{\rm ext}}^{(t)}$, $\mu^{(t)}_{k,{\rm obs}}$ and
$\gamma^{(t)}_{k}$ are the $k-$th element of vectors ${\bf x}_{{\rm
ext}}^{(t)}$, ${\boldsymbol{\mu}}^{(t)}_{\rm obs}$ and
${\boldsymbol{\gamma}}^{(t)}$, respectively. $v_{{k,{\rm
ext}}}^{(t)},\Sigma^{(t)}_{k,{\rm obs}}$, and $\lambda_{k}^{(t)}$ are the
$k$-th diagonal element of matrices ${\bf v}_{{{\rm
ext}}}^{(t)},{\boldsymbol{\Sigma}}^{(t)}_{\rm obs}$, and
${\boldsymbol{\lambda}}^{(t)}$, respectively.

Now, we estimate the soft symbols of ${\bf x}$ based on the Gaussian
likelihood function in (7). The soft symbol estimates and its variance are
given as

$\hat{{\bf x}}^{(t)}={\sf E}{\left[{\bf x}|{\bf x}_{\rm ext}^{(t)},{\bf
v}_{\rm ext}^{(t)}\right]}=c\sum_{{\bf x}\in{\boldsymbol{\Omega}}}{\bf
x}\text{ }\mathcal{N}\left({\bf x},{\bf x}_{\rm ext}^{(t)};{\bf v}^{(t)}_{{\rm
ext}}\right)$ (9a) ${\bf v}^{(t)}={\sf Var}{\left[{\bf x}|{\bf x}_{\rm
ext}^{(t)},{\bf v}_{\rm ext}^{(t)}\right]}=\\\ {\sf E}\left[\left|{\bf x}-{\sf
E}\left[{\bf x}|{\bf x}_{\rm ext}^{(t)},{\bf v}_{\rm
ext}^{(t)}\right]\right|^{2}\right],$ (9b)

respectively, where $c$ is a normalisation constant to ensure that the
summation of $p^{(t)}({\bf y}|{\bf x})$ is $1$.

The iteration of the EP receiver is performed until the values of $\hat{{\bf
x}}^{(t)}$ and ${\bf v}^{(t)}$ in (9) are close to those of
${\boldsymbol{\mu}}^{(t)}_{\rm obs}$ and ${\boldsymbol{\Sigma}}^{(t)}_{\rm
obs}$ in (6), respectively. This is referred to as the moment matching
condition [11]. When the moment matching condition has not been satisfied, the
EP receiver recomputes the parameters
$\left({\boldsymbol{\gamma}}^{(t+1)},{\boldsymbol{\lambda}}^{(t+1)}\right)$ in
(6) by using (9) and a predefined learning parameter, $\beta$. Note that
$\beta$ is used to weight the parameters,
$\left({\boldsymbol{\gamma}}^{(t+1)},{\boldsymbol{\lambda}}^{(t+1)}\right)$ in
two consecutive iterations.

Their details are described in [11]. Although EP receivers can achieve near
optimal detection performance [11, 12, 13, 14, 19, 20], it suffers a highly
computational complexity as it performs the matrix inversion operation (6b),
in every iteration. Furthermore, the predefined learning parameter needs to be
searched beforehand by using the trial and error processes. These limitations
prohibit the deployment of EP receivers in real time systems.

##  IV Linear Bayesian Learning M-MIMO Receiver

In this section, we propose a novel linear Bayesian based receiver referred to
as a linear Bayesian learning (LBL) receiver which avoids the matrix inversion
operations in (11a) and the search of a learning parameter. The proposed
receiver is shown in Fig. 2. It consists of three modules, BSO module that
computes the probability distribution function (PDF) of observed symbols from
the received signals; BSE module that yields soft symbol estimates based on
the computed PDFs; DSC module refines the symbol estimates by using the BSE
outputs and returns the refined estimates to BSO module.

###  IV-A Bayesian Symbol Observation

The computation of PDF of symbols from the observed received signal is done by
treating ${\bf x}$ in (2) at each iteration $t$ as a random vector. Its mean
and variance are obtained from observation ${\bf y}$ based on (4). The PDF of
symbol vector ${\bf x}$ for a given ${\bf y}$ at iteration $t$, $p^{(t)}({\bf
x}|{\bf y})$ is given as

$p^{(t)}({{\bf x}}|{\bf y})=\mathcal{N}\big{(}{\bf x},{\bf x}_{\rm
obs}^{(t)};{\boldsymbol{\Sigma}})$ (10)

where

$\displaystyle{\boldsymbol{\Sigma}}={\sf Var}[{\bf
x}]\propto\left({\sigma^{-2}{\rm diag}({\bf H}^{H}{\bf H})}\right)^{-1},$
(11a) $\displaystyle{\bf x}_{\rm obs}^{(t)}={{\sf E}}[{\bf
x}]={\boldsymbol{\Sigma}}\sigma^{-2}{\left({\bf H}^{H}\mathbf{y}-{\rm
off}({\bf H}^{H}{\bf H})\tilde{{\bf x}}^{(t-1)}\right)}$ (11b)

where $\Sigma_{k}$ is the $k$-th diagonal element of matrix
${\boldsymbol{\Sigma}}$, $x_{{\rm obs},k}^{(t)}$ is the $k$-th element of
vector ${\bf x}_{\rm obs}^{(t)}$. Note that ${\boldsymbol{\Sigma}}$ in (11a)
is a diagonal matrix and thus no matrix inversion operation is needed. This is
in contrast to the matrix inversion operations in EP receiver (6b). The
results, $\left({\bf x}_{\rm obs}^{(t)},{\boldsymbol{\Sigma}}\right)$ are then
forwarded to BSE module, as shown in Fig. 2.

###  IV-B Bayesian Symbol Estimator

In the BSE module, in Fig. 2, we compute the soft symbol estimate,
$\tilde{{\bf x}}^{(t)}$. We use (10) and (11) obtained from the BSO module to
write the likelihood probability function of symbol $x_{k}$ of user $k$,

$p(x_{k}|x_{{\rm
obs},k}^{(t)},\Sigma_{k})=c\times{\mathcal{N}}{\left(x_{k};x_{{\rm
obs},k}^{(t)},\Sigma_{k}\right)}$

where $c$ is a normalisation constant to ensure that
$\sum_{x_{k}\in\Omega_{k}}p(x_{k}|x_{{\rm obs},k}^{(t)},\Sigma_{k})=1$. The
soft symbol estimate, $\tilde{x}_{k}^{(t)}$ is then given as

$\tilde{x}_{k}^{(t)}={\sf E}{\left[x_{k}|x_{{\rm
obs},k}^{(t)},\Sigma_{k}\right]}=\sum_{x_{k}\in\Omega_{k}}x_{k}{p}{\left(x_{k}|x_{{\rm
obs},k}^{(t)},\Sigma_{k}\right)},$ (12)

$\tilde{{\bf x}}^{(t)}=[\tilde{x}_{1}^{(t)},\dots,\tilde{x}_{K}^{(t)}]$ and
$\Omega_{k}$ is the $k$-th constellations. The errors between estimations and
observations for $K$ users at iteration $t$, ${\bf V}^{(t)}$ is given as

${\bf V}^{(t)}=\\\ {\rm diag}\left(\left({\bf H}^{H}{\bf y}-({\bf H}^{H}{\bf
H})\tilde{{\bf x}}^{(t)}\right)\left({\bf H}^{H}{\bf y}-({\bf H}^{H}{\bf
H})\tilde{{\bf x}}^{(t)}\right)^{T}\right).$ (13)

The outputs of BSE module $\left(\tilde{{\bf x}}^{(t)},{\bf V}^{(t)}\right)$
from (12) and (13) are sent to the DSC module.

###  IV-C Decision Statistic Combining

In this section, we use the DSC scheme to approximate the learning parameter
based on the values of instantaneous errors in subsequent iterations. This is
shown in Fig. 2 and done by weighting the outputs of the BSE module,
$\tilde{{\bf x}}$ in the current and previous iterations based on the error
${\bf V}$ in (13), given as

${\bf x}_{\rm DSC}^{(t)}=\left({\bf V}^{(t)}+{\bf
V}^{(t-1)}\right)^{-1}\left({\bf V}^{(t-1)}\tilde{{\bf x}}^{(t)}+{\bf
V}^{(t)}\tilde{{\bf x}}^{(t-1)}\right)$ (14)

where ${{\bf x}}_{\rm DSC}^{(t)}=[x_{1,{\rm DSC}}^{(t)},\dots,x_{K,{\rm
DSC}}^{(t)}]$. The iterative process will stop if the following is achieved,

$\|{\bf x}_{\rm DSC}^{(t)}-{\bf x}_{\rm DSC}^{(t-1)}\|\leq\epsilon\lor
t=T_{\rm max}$ (15)

where $\epsilon$ is the minimum acceptable difference of the ${\bf
x}^{(t)}_{\rm DSC}$ in two consecutive iterations and $T_{\rm max}$ is the
maximum number of the iterations. If (15) is not satisfied, we use ${\bf
x}_{\rm DSC}^{(t)}$ as the inputs of BSO in the next iteration

$\tilde{{\bf x}}^{(t)}={{\bf x}}_{\rm DSC}^{(t)},$ (16)

else the symbol hard decision is then made by using the outputs of DSC in
(14),

${\bf x}_{\rm hard}=\operatornamewithlimits{argmin}\limits\|{\bf x}_{{\rm
DSC}}^{(T)}-{\boldsymbol{\Omega}}\|^{2}$ (17)

where $T$ is the last number of iteration and ${\boldsymbol{\Omega}}$ is the
set of the constellations. Note that the use of (13) and (14) eliminates the
trial and error issues for finding optimal learning parameters used in most of
Bayesian learning iterative receivers, Eq. (37) and (38) in [11], Eqs. (15)
and (16) in[13], and Eq. (8) in [16]. The complete pseudo-code is shown in
Alg. 1.

Initialization: $\tilde{{\bf x}}^{(0)}=0,T_{\rm max}=10$;

while  _( 15) is not satisfied_ do

1\. Compute the mean and variance,
${\boldsymbol{\mu}}^{(t)},{\boldsymbol{\Sigma}}$, in (11)

2\. Compute the mean and variance, $\hat{{\bf x}}^{(t)},\hat{{\bf V}}^{(t)}$,
in (12),(13)

3\. Compute ${\bf x}_{\rm DSC}^{(t)}$, in (14)

4\. Execute (16)

end while

5\. Execute (17)

Algorithm 1 LBL Algorithm

##  V Complexity and Convergence Analysis

Figure 3: The Convergence behaviour of the LBL receiver.  TABLE I: Computational complexity comparison. Receiver | Complexity
---|---
LBL | $\mathcal{O}(NKT)$
AMP [15] | $\mathcal{O}(NKT)$
MMSE | $\mathcal{O}(N^{2}K+NK)$
EP in [11] | $\mathcal{O}((N^{2}K+NK)T)$
EP in [14] | $\mathcal{O}((NK^{2}+NK)T)$
ML [8] | $\mathcal{O}(M^{K})$

We first evaluate the computational complexity of the LBL receiver.The
computational complexity of the proposed receiver; the representative of best
linear receivers, MMSE [9] scheme; the representative of Bayesian receivers EP
[11, 14] and AMP [15] schemes; and the exhaustive search based ML scheme [8]
are tabulated in Table I. The table indicates that the computational
complexity of the proposed receiver increases linearly with the number of
antennas, $N$ and users, $K$ by avoiding matrix inversion operations. This is
in contrast to all other receivers whose computational complexity grow
exponentially with $N$ and/or $K$. We have also observed that the number of
iterations, $T$ needed for the EP, AMP, and LBL receivers are similar.
However, the convergence rates of the AMP and EP receivers are not shown in
this paper due to a space limitation. Therefore, the proposed receiver has a
significantly lower processing latency and thus more suitable for URLLC data
traffic. For instance, for $K=64,N=192,\text{and }T=10$, LBL receiver
complexity is 19 times lower than MMSE receiver, and 192 times lower than EP
receiver. We conclude that by eliminating the matrix inversion operations, we
can obtain a significantly lower computational complexity and thus much lower
latency.

To analyse the convergence rate of the proposed LBL receiver, we plot the
value of the difference between symbol estimates, ${\bf x}_{{\rm DSC}}$ in
iteration $t$ and $t-1$ versus the iteration $t$ for $100,000$ channel
realisations in Fig. 3, and $\epsilon=10^{-6}$. The figure shows that the
maximum number of iterations needed is $6$ for various system configurations,
implying a low latency. This also indicates that the iterations needed is
relatively insensitive to the system configurations.

##  VI Reliability Analysis

To analyse receivers’ reliability, we plot the BER performances of our
proposed receiver versus the linear receivers based on PIC [10] and MMSE [9]
schemes; Bayesian learning receivers based on EP [11, 14] and AMP[15] schemes;
an optimal ML receiver [8] for different system configurations. The same
modulation scheme and number of channel realisation as in Section V are used.
Note also that we perform trial and error process to search for the optimal
learning parameters for EP and AMP off-line. We found that the optimal
learning parameters in the EP and AMP receivers are $\beta=0.9$ and $0.2$,
respectively.

In Fig. 4a, we employ the number of receive antennas, $N=24$ and users, $K=8$.
It shows that the LBL receiver outperforms (MMSE and PIC) receivers for about
2 dB to achieve BER $10^{-4}$. In practice, this also implies that the radio
coverage of the proposed M-MIMO receiver is $26\%$ broader than those of the
existing linear receivers. Note also the similarity in term of BER performance
with the EP and ML receivers that have significantly higher computational
complexity and latency as explained in the previous section. In Fig. 4b, the
number of received antennas and users are increased to $N=192$ and $K=64$.
Similar behaviour is observed in relation to the BER performance in comparison
with the linear and Bayesian learning receivers. Note that in Fig. 4b we are
unable to draw ML receivers due to its high computational complexity. From
these facts, we conclude that the LBL receiver is more reliable than several
existing classical and advance iterative receivers.

In Fig. 5, we analyse the relationship between the number of users and used
received antennas in the proposed M-MIMO receiver. We first define the ratio
between the number of users and used received antennas as $\alpha$ in Fig. 5.
In Fig. 5a, we set a fixed number of antennas, $N=256$. The EP and LBL
receivers can achieve BER less than $10^{-4}$ when $\alpha<0.71$ and
$\alpha<0.65$, respectively. In Fig. 5b, we increase the number of antennas,
$N=1024$ and obtain $\alpha<0.61$ and $0.7$ for LBL and EP receivers,
respectively. This implies that the exchange of the low complexity signal
processing of the proposed receiver is the number of supported users.

(a)  $N=24,K=8$, $M=4.$ (b) $N=192,K=64$, $M=4.$

Figure 4: The BER performance of PIC, MMSE, AMP, EP, LBL, ML receivers in
several system configurations.

(a)  $N=256$ and $M=4.$ (b) $N=1024$ and $M=4.$

Figure 5: The BER performance of LBL and EP receivers with respect to the
ratio of users to antennas, $\alpha$.

##  VII Conclusion

We propose a novel Bayesian-based parallel interference cancellation receiver
referred to as the LBL receiver. Simulation results show that the BER
performance of the LBL receiver is very close to that of maximum likelihood
receiver while maintaining a linear latency processing time in contrast to
other existing schemes in the literature.

## References

  * [1] A. Kosasih, W. Hardjawana, C. K. Wen, and B.Vucetic, “Bayesian-based parallel interference cancellation receiver for massive MIMO systems,” _in preparation for a journal_.
  * [2] L. Lu, G. Y. Li, A. L. Swindlehurst, A. Ashikhmin, and R. Zhang, “An overview of massive MIMO: Benefits and challenges,” _IEEE J. Sel. Topics Signal Process._ , vol. 8, no. 5, pp. 742–758, Oct. 2014.
  * [3] Y. Huo, X. Dong, and W. Xu, “5G cellular user equipment: From theory to practical hardware design,” _IEEE Access_ , vol. 5, pp. 13 992–14 010, Jul. 2017.
  * [4] C. Sexton, N. J. Kaminski, J. M. Marquez-Barja, N. Marchetti, and L. A. DaSilva, “5G: Adaptable networks enabled by versatile radio access technologies,” _IEEE Commun. Surveys Tuts._ , vol. 19, no. 2, pp. 688–720, 2nd Quart. 2017.
  * [5] T. K. Vu, C. Liu, M. Bennis, M. Debbah, M. Latva-aho, and C. S. Hong, “Ultra-reliable and low latency communication in mmWave-enabled massive MIMO networks,” _IEEE Commun. Lett._ , vol. 21, no. 9, pp. 2041–2044, Sep. 2017.
  * [6] P. Popovski, Č. Stefanović, J. J. Nielsen, E. de Carvalho, M. Angjelichinoski, K. F. Trillingsgaard, and A. Bana, “Wireless access in ultra-reliable low-latency communication (URLLC),” _IEEE Trans. Commun._ , vol. 67, no. 8, pp. 5783–5801, Aug. 2019.
  * [7] A. Bana, G. Xu, E. D. Carvalho, and P. Popovski, “Ultra reliable low latency communications in massive multi-antenna systems,” in _Proc. 52nd Asilomar Conf. Signals, Syst. and Comput._ , USA, Oct. 2018, pp. 188–192.
  * [8] U. Fincke and M. Pohst, “Improved methods for calculating vectors of short length in a lattice, including a complexity analysis,” _Math. Comput._ , vol. 44, no. 170, p. 463–471, Apr. 1985.
  * [9] G. Caire, R. Muller, and T. Tanaka, “Iterative multiuser joint decoding: Optimal power allocation and low-complexity implementation,” _IEEE Trans. Inf. Theory_ , vol. 50, no. 9, p. 1950–1973, Sep. 2004.
  * [10] N. Aboutorab, W. Hardjawana, and B. Vucetic, “A new iterative doppler-assisted channel estimation joint with parallel ICI cancellation for high-mobility MIMO-OFDM systems,” _IEEE Trans. Veh. Technol._ , vol. 61, no. 4, pp. 1577–1589, May 2012.
  * [11] J. Céspedes, P. M. Olmos, M. Sánchez-Fernández, and F. Pérez-Cruz, “Expectation propagation detection for high-order high-dimensional MIMO systems,” _IEEE Trans. Commun._ , vol. 62, no. 8, pp. 2840–2849, Aug. 2014.
  * [12] K. Takeuchi and C. K. Wen, “Rigorous dynamics of expectation-propagation signal detection via the conjugate gradient method,” in _Proc. IEEE 18th Int. Workshop Signal Process. Advances in Wireless Commun._ , Japan, July 2017, pp. 1–5.
  * [13] G. Yao, G. Yang, J. Hu, and C. Fei, “A low complexity expectation propagation detection for massive MIMO system,” in _Proc. IEEE Global Commun. Conf._ , UAE., Dec. 2018, pp. 1–6.
  * [14] H. Wang, A. Kosasih, C. K. Wen, J. Shi, and W. Hardjawana, “Expectation propagation detector for extra-large scale massive MIMO,” _submitted for publication_.
  * [15] C. Jeon, R. Ghods, A. Maleki, and C. Studer. (2018) Optimal data detection in large MIMO. [Online]. Available: https://arxiv.org/abs/1811.01917
  * [16] J. Zhang, C. K. Wen, S. Jin, and G. Y. Li, “Artificial intelligence-aided receiver for a CP-free OFDM system: Design, simulation, and experimental test,” _IEEE Access_ , vol. 7, pp. 58 901–58 914, May 2019.
  * [17] Minka, Thomas P, “Expectation propagation for approximate Bayesian inference,” in _Proc. 17th Conf. Uncertainty in Artificial Intell._ , USA, Aug. 2001, pp. 362–369.
  * [18] X. Yuan, J. Ma, and L. Ping, “Energy spreading transform based MIMO systems: Iterative equalization, evolution analysis, and precoder optimization,” _IEEE Trans. Wireless Commun._ , vol. 13, no. 9, pp. 5237–5250, Sept. 2014.
  * [19] Y. Zhang, Z. Wu, C. Li, Z. Zhang, X. You, and C. Zhang, “Expectation propagation detection with neumann-series approximation for massive MIMO,” in _Proc. IEEE Int. Workshop Signal Process. Syst._ , South Afr., Oct. 2018, pp. 59–64.
  * [20] X. Tan, Y. Ueng, Z. Zhang, X. You, and C. Zhang, “A low-complexity massive MIMO detection based on approximate expectation propagation,” _IEEE Trans. Veh. Technol._ , vol. 68, no. 8, pp. 7260–7272, Aug. 2019.