# A tunable quantum random number generator based on a fiber-optical Sagnac
interferometer

Joakim Argillander1, Alvaro Alarcón1 and Guilherme B. Xavier1 1 Institutionen
för Systemteknik, Linköping University, 581 83 Linköping, Sweden
<EMAIL_ADDRESS>
###### Abstract

Quantum random number generators (QRNG) are based on the naturally random
measurement results performed on individual quantum systems. Here, we
demonstrate a branching-path photonic QRNG implemented with a Sagnac
interferometer with a tunable splitting ratio. The fine-tuning of the
splitting ratio allows us to maximize the entropy of the generated sequence of
random numbers and effectively compensate for tolerances in the components. By
producing single-photons from attenuated telecom laser pulses, and employing
commercially-available components we are able to generate a sequence of more
than 2 gigabytes of random numbers with an average entropy of 7.99 bits/byte
directly from the raw measured data. Furthermore, our sequence passes
randomness tests from both the NIST and Dieharder statistical test suites,
thus certifying its randomness. Our scheme shows an alternative design of
QRNGs based on the dynamic adjustment of the uniformity of the produced random
sequence, which is relevant for the construction of modern generators that
rely on independent real-time testing of its performance.

  * January 2022

Keywords: quantum random number generation, tunable beamsplitter, fiber-optic
Sagnac interferometer

##  1 Introduction

Quantum information science is an interdisciplinary field where computer
science meets physics, and where long lasting principles of information
transfer, cryptography and secret sharing are challenged [1]. Quantum random
number generators (QRNG) are mature devices that directly benefit from the
properties of quantum systems in order to produce truly random numbers [2]. By
harnessing the fundamental randomness from the outputs of quantum mechanical
measurements it is possible to generate sequences of fundamentally
unpredictable random numbers. A popular implementation for QRNGs is based on
performing measurements on single-photons. Typically, in these photonic QRNGs,
one degree of freedom of a photon is used to prepare a superposition state of
the form $\ket{\psi}=\alpha\ket{0}+\beta\ket{1}$, where
$\left|\alpha\right|^{2}+\left|\beta\right|^{2}=1$, which is then projected
onto the basis formed by $\ket{0}$ and $\ket{1}$. By letting one of the basis
vectors represent a binary “0” and the other a binary “1”, we can utilize the
fact that measurements of an equally weighted superposition state yield
uniformly distributed random measurement outcomes to generate a random
sequence by performing repeated measurements.

One type of photonic QRNG is known as a branching-path generator [2, 3, 4, 5,
6], and is typically implemented either using polarized light and polarizing
beamsplitters, or with unpolarized light and conventional beamsplitters.
Regardless of the scheme, it is of utmost importance that the beamsplitters’
splitting ratio is $50:50$ for transmitted and reflected light respectively,
in order to guarantee an unbiased sequence. In practice, due to manufacturing
limitations, the splitting ratio will not be perfectly balanced, and this
leads to non-uniformity of the bit distribution of the generated sequence.
Furthermore, the splitting ratio may also change over time, due to temperature
fluctuations and other environmental factors, as well as component
degradation. Branching-path generators may also suffer from detector
imbalances, where two identical detectors used at each of the beamsplitter’s
outputs have slightly different detection efficiencies that also lead to a
bias in the sequence [2]. Such a bias reduces the unpredictability of the
sequence, and also so does its usability, not only for cryptography, but also
for uses like Monte-Carlo simulations and gambling. For cryptography
specifically, the security is directly proportional to the randomness of the
secret key, and therefore a low-entropy sequence increases the probability
that an attacker is able to successfully guess the key, which thus reduces the
effective strength of the cipher [7].

In this article, we demonstrate a tunable-ratio branching-path photonic QRNG
by implementing a dynamically tunable beamsplitter using a fiber-optical
Sagnac interferometer. The tunable beamsplitting function is implemented with
an active phase modulator placed in a Sagnac loop working as a single-photon
router [8]. By modulating the phase of a light packet propagating in the
clockwise direction vs the counter-clockwise one, we are able to tune the
output probabilities from the Sagnac interferometer through a change in the
amplitude of the driving signal fed to the phase modulator. This allows us to
change the entropy of the generated sequence as we wish. Both outputs are then
detected through a time-division multiplexing scheme to allow detection with
only one single-photon detector. Interferometers in a Sagnac-type
configuration have previously been used to modulate polarization [9, 10] and
intensity [11] in order to prepare states for quantum key distribution
experiments. While building on the same principle, our system is the first
demonstration of interferometric routing of single photons for quantum random
number generation. In this work we are able to generate a continuous stream of
random numbers over 36 hours of continuous operation at an average rate of 131
kbit/s. The generation rate and entropy are highly stable over time due to the
intrinsic phase stability of the Sagnac interferometer configuration. We
finally demonstrate the randomness of the sequence by passing it successfully
through the stringent statistical test suites NIST 800-22 [12] and Dieharder
[13] aimed at assessing random number generators. A critical issue in QRNGs,
especially the ones aimed at cryptographic applications, is whether the user
can trust the device. A solution to this problem is performing independent
real-time testing of the internal devices of the QRNG, usually through the
user changing either the prepared quantum states, or a measurement parameter
[14].  In order to test the device and achieve high confidence that the
generated sequence is private, a high interferometric visibility is necessary
[15]. The dynamic tunability provided by our setup is therefore useful to
tackle this issue, and is particularly appealing to recently demonstrated
measurement-device-independent QRNGs [16, 15]. Our results show that Sagnac
interferometers can be successfully employed as tunable beamsplitters for the
construction of quantum random number generators, thus opening up for
implementation of state-of-the-art measurement-device independent quantum
random number generators.

##  2 The experiment

The principle for a QRNG scheme using a TBS is based on a source of strongly
attenuated pulses, a beamsplitter with a variable splitting ratio and single-
photon detection of the two outputs of the beamsplitter (figure 1). A
detection at one output is assigned the bit “0”, while “1” is assigned to the
other one. The entire system works synchronously, driven by Field Programmable
Gate Array (FPGA) electronics, which also performs the data acquisition.  If
both outputs are simultaneously detected (either from two photons generated
together, or from noise sources), the FPGA discards that result in order to
remove the effects of after-pulsing and multiphoton events. Similarly, if no
photons are detected for a generated pulse (due to absorption) no result is
recorded. The generated bits are then streamed from the FPGA to a personal
computer for randomness extraction and storage. The bitstream then passes
through statistical tests to certify its randomness.

Figure 1:  Principle of a tunable-ratio beamsplitter for quantum random number
generation. Light from a photon source is split using a tunable beamsplitter
(TBS). Detectors (D1, D2) placed at the outputs of the TBS register photon
events and record a binary “0” or a binary “1” depending on whether the photon
impinges on D1 or D2. The choice of assignment of “0”/“1” to D1/D2 is
arbitrary. Figure 2: Experimental setup. The CW-laser is first fed to an
intensity modulator (IM) that chops the laser into 20 ns pulses, which are
then attenuated using a variable attenuator (ATT) to weak coherent pulses
(WCP) with an average photon probability of approximately $10$ photons per
pulse. The WCPs are then fed to a Sagnac interferometer that acts as a tunable
beamsplitter where the superposition
$\ket{\psi}=1/\sqrt{2}\left(\ket{\textrm{CW}}+ie^{i\phi}\ket{\textrm{CCW}}\right)$
is prepared, and CW and CCW denote the photon state of propagation inside the
Sagnac loop. At the output of the Sagnac loop, we use a  150 meter fiber-optic
delay line ($L_{TB}$) to add $\approx 750$ ns of delay to one term,
effectively time-bin multiplexing the two outputs from the Sagnac
interferometer into an _early_ and a _late_ time-bin for detection using one
single-photon detector.

The key aspect in our experiment is the use of a tunable beamsplitter (figure
1). Although in principle it is possible to create one by adding extra
attenuators in each output of a conventional beamsplitter, or by using a half-
wave plate in tandem with a polarizing beamsplitter, these options do not
provide possibilities of fast tuning. Furthermore, in the case of the first
option, extra losses are added to the system. Here we pursue an
interferometric approach to implement tunability. For instance, a Mach-Zehnder
interferometer can realize a tunable beamsplitter, with the tunability
provided by setting the relative phase $\phi$ between the interferometer arms
[17]. The output probability of the single-photon is then given by a
sinusoidal function of $\phi$. However, the drawback of using conventional
Mach-Zehnder interferometers is that they suffer from intrinsic phase
instability, and therefore require either mechanical isolation (i.e. dampened
tables, environmental protection, etc) [18] or active stabilization [19, 20].
In our work we follow the approach of [8] where we employ a Sagnac
interferometer to provide the same functionality, with the added benefit of
intrinsic phase stability. In our Sagnac interferometer a high-speed telecom
fiber-pigtailed phase modulator is placed to give a relative phase shift
$\phi$ between the two propagation directions (figure 2). This is done by
synchronizing the phase shift with the wave packet, such the phase shift is
only applied in one direction by the modulator.

We employ a continuous wave (CW) distributed feedback telecom laser diode,
centered at $\lambda$ = 1546.92 nm, connected to variable optical attenuators
(ATT) to produce weak coherent states. The FPGA board generates electrical
pulses at a repetition rate of $250$ kHz in order to drive an intensity
modulator which chops the CW input light into optical pulses. The electronic
pulses are shaped to $20$ ns width, as limited by the employed electronics
driver before being fed to the intensity modulator. The attenuated optical
pulses pass through an optical circulator before entering one of the input
ports of a bidirectional $50:50$ fiber coupler. The output ports form the
Sagnac loop, together with a phase modulator and a 50 m long optical fiber
spool ($L_{s}$) that works as a delay line (figure 2). The delay line is
needed to ensure there is sufficient time separation for the phase modulation
signal to act on the wave packet propagating in only in one direction and not
in the opposing one, thus creating the necessary relative phase shift to
generate the different output probabilities in the interferometer outputs [8].
Not shown in figure 2 are two manual polarization controllers added in the
Sagnac loop, where one is used to align the pulses’ polarization with the
input polarization of the phase modulator, while the second one is used to
ensure that the two counter-propagating terms have the same polarization at
the recombining beamsplitter.

Inside the Sagnac loop the attenuated pulses can be represented as a coherent
superposition of the two counter-propagating paths as
$\ket{\psi}=1/\sqrt{2}(\ket{\textrm{CW}}+ie^{i\phi}\ket{\textrm{CCW}})$ where
$\ket{\textrm{CW}}$ and $\ket{\textrm{CCW}}$ represent the clockwise and
counter-clockwise propagation paths respectively. The detection probabilities
at the outputs A and B of the Sagnac loop are proportional to
$\cos^{2}(\frac{\phi}{2})$ and $\sin^{2}(\frac{\phi}{2})$. For voltages
corresponding to $\phi=\pi/2,3\pi/2$ the Sagnac works as a beamsplitter with a
perfectly balanced splitting ratio. The phase shift $\phi$ is proportional to
the modulator’s driving voltage, thus allowing us to precisely control the
splitting ratio between 0 and 1 in the ideal case. We fine-tune the desired
ratio by applying a small offset $\delta V_{\phi}$ to the theoretical voltage
value $V_{\phi}$. Therefore, the voltage we apply to achieve a balanced
splitting ratio is $\overline{V}_{\pi/2}=V_{\pi/2}+\delta V_{\pi/2}$.  Due to
insertion losses and propagation losses in the time-bin multiplexer, a
balanced splitting ratio between the early and late time bin does not imply a
50:50 probability for a photon to emerge at the outputs A:B of the Sagnac
interferometer. The two outputs of the Sagnac (denoted A and B in figure 2)
are then routed to the detection system (with output A passing through the
circulator to decouple it from the attenuated pulses coming from the source,
propagating in to the Sagnac loop). We opted to use one SPD and mapped each
output path to an early or a late time-bin, with both time slots detected
sequentially by the SPD. The time-bin separation is created by a 750 ns time
delay (a $\approx$150 m optical fiber spool) placed at the A output path from
the Sagnac loop. The two paths are then recombined with a 50:50 bidirectional
fiber coupler, and one of its outputs is connected to the SPD (idQuantique
id210), operating in Geiger mode with 10% detection efficiency, 2.5 ns gate
width and approximately $10^{-5}$ dark count probability per gate. The SPD is
externally gated by the FPGA by sending pairs of pulses, each separated by the
time-bin delay, that are synchronized with the signal sent to the intensity
modulator that generates the photon pulses. Each of these pulses correspond to
a gate in the _early_ and _late_ bin respectively. The output signals from the
SPD are also read by the FPGA for data acquisition. We set the deadtime of the
SPD to 500 ns, so that it is shorter than the time separation between the
_early_ and _late_ time-bins. We then let a detection of a photon in the
_early_ time-bin represent a “0” and a detection in the _late_ time-bin
represent a “1”.

##  3 Results

First we demonstrate the variable splitting ratio of the Sagnac interferometer
acting as a tunable beamsplitter. We set the attenuators such that we maximize
the detection probability at the detector for both time-bins, while still
operating in a linear regime (i.e. the detector is not saturating). This is
observed from figure 3a, where the number of detection events per integration
time fits well with the theoretical prediction that output A [B] from the
Sagnac loop should follow a $\cos^{2}(\frac{\phi}{2})$
[$\sin^{2}(\frac{\phi}{2})$] distribution. We also observe from these results
that the Sagnac interferometer is operating correctly, and complementary
output behaviour of the two outputs is observed, as expected.

(a)

(b)

Figure 3:  a) Characterization of detected photon events measured in the
_early_ and _late_ time-bin, as a function of applied voltage to the phase
modulator. The error bars are calculated assuming Poissonian photon statistics
($\sigma=\sqrt{N}$ for a photon count of $N$). Also plotted is a fit of the
theoretical detection probabilities that are proportional to $\cos^{2}$
[$\sin^{2}$] for the _early_ [_late_] time-bin, respectively. b) Raw entropy
(before extraction) as a function of applied voltage to the phase modulator
Also shown in the inset is a plot of the entropy versus applied voltage to the
phase modulator with a smaller step size.

We then wish to optimize the operating point of the Sagnac loop for optimal
random number generation. We scan the driving voltage of the phase modulator
and measure the Shannon entropy of the generated bits when considered as
$8$-bit symbols. The Shannon entropy $\mathcal{H}$ is defined as
$\mathcal{H}=-\sum_{i=0}^{255}p_{i}\log_{2}{(p_{i})}$, where $p_{i}$ is the
probability of finding the $i_{th}$ 8-bit sequence. We increase the voltage
$V_{\phi}$ applied to the phase modulator in the range $0\leq V_{\phi}\leq
4.2$ in steps of $200$ mV and measure the entropy, as can be seen in figure
3b. An inset is shown in figure 3b where we also sweep the voltage between
$2-2.3$ volts in steps of $20$ mV, and likewise measure $\mathcal{H}$ to
identify the optimal operating point where the entropy is maximized. It can be
clearly seen that the maximum obtained entropy (7.98 bits/byte) occurs for an
applied driving voltage of around $2.15$ V.

(a)

(b)

Figure 4:  Characterization of the stablity of the raw bitrate a) and entropy
stability b) over a period of 36 hours. Each of the data points represents the
average over 30 minutes.

(a)

(b)

Figure 5:  Results from statistical testing of randomness from NIST STS
800-22[12], a) and Dieharder[13] b). In figure 5a we plot the proportion of
tests that have passed, along with the confidence interval that the proportion
must fall within. For figure 5b, each colored vertical line corresponds to a
p-value from the tests. Lines in green are p-values that fall outside of the
confidence bound $10^{-2}<p<10^{0}$. The red lines are the resulting p-value
of a Kolmogorov-Smirnov test done on test results that yield multiple
p-values.

With the entropy maximized, we then proceed to acquire a long continuous
sequence of bits to feed to the statistical randomness tests, and also to
demonstrate the robustness of our implementation. The raw bit sequence is
buffered in the FPGA in blocks of 32 kilobytes and then transmitted over an
Ethernet connection to a personal computer for storage. Before assessing the
sequence, it is passed through a Toeplitz randomness extraction procedure in
order to remove any remaining non-uniformities from the raw data. To extract
randomness, we split the generated sequence into $N$ sequences of length $n$,
and multiply each subsequence with an $n\times m$ binary Toeplitz matrix which
yields an output of $m$ hashed bits [21]. We choose to use $n=400$, and as
defined in [21], we obtain $m=\mathcal{H}_{\mathrm{min}}-2\log_{2}\varepsilon$
where $\mathcal{H}_{\mathrm{min}}$ is the min-entropy of 8-bit strings which
is defined as
$\mathcal{H}_{\mathrm{min}}=-\log_{2}\left[\max_{x\in\left\\{0,1\right\\}^{8}}\mathrm{Pr}\left\\{X=x\right\\}\right]$.
We measure a min-entropy of $7.7451$ from the acquired sequence before
extraction. We choose a security parameter $\varepsilon=2^{-100}$ which is
derived from the leftover hash lemma [22]. We seed the Toeplitz matrix $T$ at
the beginning of the extraction procedure with $n+m-1$ bits taken from the raw
sequence, as it is sufficient to only specify the first row and first column.
The matrix $T$ can be reused for each of the subsequences, and we obtain the
extracted sequence of bits by concatenating the results of the multiplications
$Tn_{1},Tn_{2},Tn_{3},...,T_{N}$, where the subscript represents the
subsequence index [21, 22]. The extraction procedure is lossy, and for our
choice of parameters, we are able to achieve an efficiency of $\approx 50\%$.
After extraction we calculate the min-entropy of the sequence to be $7.9982$
bits.

While recording the raw data, we also measure the entropy of each streamed 32
kilobyte block, and as can be seen in figure 4a, we achieve an average entropy
of $7.9925\pm 0.0014$ bits/byte over $36$ hours, which demonstrates
exceptional phase stability of the QRNG. Likewise, we also measure the raw
bitrate when streaming the data, and we achieve an average bitrate of
$131.8\pm 2.3$ kbps (figure 4b). The bitrate is mainly limited by the delays
in the electronics that amplify the signals from the FPGA in order to drive
the phase modulator, and in the triggering electronics which is necessary to
synchronize the generation of pulses with the detection. By improving the
driver electronics, the bitrate is only limited by the separation of the time
bins.

After generating $2.2$ gigabytes of data, and extracting it, we pass the
sequence through the NIST 800-22 Statistical Test Suite[12] (NIST test) and
the Dieharder[13] test suite. Both test suites subject the random sequence to
a number of statistical tests aimed at assessing the randomness of any given
sequence. We run the NIST tests with $1000$ bitstreams of length $1$ Mb. To
deem any given test as passed, we form a confidence bound as defined in [12]
for the proportion of bitstreams that must yield a p-value $p\geq 10^{-2}$
using the formula $\hat{p}\pm 3\sqrt{\frac{\alpha(1-\alpha)}{n}}$ where
$\alpha=10^{-2}$, $n=1000$ is the number of bitstreams and $\hat{p}=1-\alpha$.
This gives us a confidence interval of $0.980\leq r<0.999$, and if the ratio
of tests that have passed falls outside this range, there is evidence that the
sequence is not random [12]. Figure 5a shows the proportions of tests that
pass, for any given test category, along with the significance level.

The Dieharder test suite yields a verdict (”Fail”, ”Weak”, ”Pass”) for each
test, and a test is deemed to pass if the p-value associated with it is in the
interval $10^{-2}<p<1-10^{-2}$. As can be seen in figure 5b, there are
parameterized tests where each test yields multiple p-values (one for each
value of the parameter). We plot all outputted p-values, and compute a
resulting p-value from a Kolmogorov-Smirnov test of all the p-values. As is
well-known, the p-values should be uniformly distributed under the null
hypothesis that the sequence is indeed random, and the Kolmogorov-Smirnov test
is used to eliminate false rejection of the null hypothesis (Type-I error)
[23]. For three of the tests, we observe that the resulting p-values fall
outside the confidence bound, and for those values of the parameter, the test
verdict is ”weak”. Since the Dieharder suite is a statistical test, it is
probable that there will be occasional ”weak” results. If the generator truly
is weak, then it will certainly fail other tests [13].

Figure 6:  Achieved visibilities for the states $\ket{\Psi}$ and $\ket{\Phi}$
that are expected to yield deterministic measurement results, together with
the visibilities for the balanced state $\ket{\Omega}$ used for randomness
generation.

To further demonstrate that the tunability can be used as a building block in
measurement-device independent protocols we run a prepare-and-measure type
protocol where we use the TBS to prepare either of the orthogonal states
$\ket{\Psi}$ and $\ket{\Phi}$ which are expected to yield deterministic
measurement results, and a state $\ket{\Omega}$ which yields random (but
balanced) measurement outcomes. We prepare the states by applying different
voltages to the PM inside the Sagnac for each state. From figure 3a we
determine the voltages to be $V_{\Psi}=0\textrm{ V}$, $V_{\Phi}=4.2\textrm{
V}$ and $V_{\Omega}=2.15\textrm{ V}$. For each acquired 32 kilobyte block of
data we then randomly choose between preparing either of the states with a
probability
$\mathrm{Pr}\left\\{\ket{\Psi}\right\\}=\mathrm{Pr}\left\\{\ket{\Phi}\right\\}=0.005$
and $\mathrm{Pr}\left\\{\ket{\Omega}\right\\}=0.99$, where
$\mathrm{Pr}\left\\{\ket{\bullet}\right\\}$ is the probability to prepare a
certain quantum state, in order to demonstrate the switching of states
necessary to implement measurement-device independent QRNG protocols. In
figure 6 we show the recorded visibilities for each of the states as a
function of time, where the visibility $\nu_{x}\in\left[0,1\right]$ for the
state $x\in\left\\{\ket{\Psi},\ket{\Phi},\ket{\Omega}\right\\}$ is defined
$\nu_{x}=\left|\frac{n_{\mathrm{early}}-n_{\mathrm{late}}}{n_{\mathrm{early}}+n_{\mathrm{late}}}\right|$
where $n_{\mathrm{early}}$ and $n_{\mathrm{late}}$ are the detector counts in
the early and late bin respectively. We are able to achieve a mean visibility
of $\nu_{\Psi}=0.997\pm 3.3\cdot 10^{-4}$, $\nu_{\Phi}=0.990\pm 3.1\cdot
10^{-3}$ and $\nu_{\Omega}=0.004\cdot 10^{-2}\pm 7.7\cdot 10^{-3}$ for each of
the three states without subtracting dark counts.

##  4 Conclusions

We have proposed and experimentally demonstrated a photonic quantum random
number generator based on a dynamically tunable beamsplitter. The random
numbers are generated depending on which output from the beamsplitter the
single-photons are detected. By tuning the beamsplitting ratio we are able to
maximize the generated entropy of the random sequence. This also helps
maximize the entropy in spite of differences in losses following the two
outputs. We implemented the tunable beamsplitter through a Sagnac
interferometer, in which the tunability is given by a relative phase between
the two counter-propagating paths, applied with an electro-optical telecom
phase modulator. This technique allows ultra-fast tunability, due to the short
response time of the phase modulator.

We are able to generate a continuous stream of random numbers over 36 hours of
continuous operation, showing an average raw generation rate of 131 kbit/s and
an average generated entropy of 7.9925 bits per byte in the raw sequence, thus
demonstrating a stable source of entropy.  As the bitrate is mainly limited by
internal delays in the employed driving electronics, by improving the driver
circuitry, a considerable bitrate gain can be achieved. By reducing these
delays, we would also reduce the lower limits on the time-bin spacing and thus
the repetition rate can be increased. Our setup is robust and built using only
off-the-shelf telecom fiber components, thus also showing the practicability
of our setup. Our results can aid in future designs of dynamically tunable
quantum random number generators, and serve as a foundation for measurement-
device quantum random number generation due to ease of changing the
beamsplitting ratio, through the appropriate selection of the applied phase in
the Sagnac loop.

Our system also demonstrates that the tunable beamsplitter enables realization
of measurement-device independent quantum random number generators thanks to
the high visibility when measuring orthogonal states, as well the fact that
the visibilities show minimal degradation over time.

## Acknowledgments

We acknowledge financial support from CENIIT Linköping University, the Swedish
Research Council (VR 2017-04470), QuantERA grant SECRET (VR 2019-00392), the
Knut and Alice Wallenberg Foundation through the Wallenberg Center for Quantum
Technology (WACQT).

## Author Contributions

J.A. and A.A. conceived the work together with G.B.X. J.A. developed the
electronic control system, and A.A. developed the optical system. Data
analysis was done by J.A. with assistance from G.B.X. The work was coordinated
by G.B.X. All authors contributed to the writing of the manuscript.

## Competing Interests

The authors declare no competing interests.

## References

## References

  * [1] Michael Nielsen and Isaac L Chuang  “Quantum Computation and Quantum Information: 10th Anniversary Edition”  In _Cambridge University Press, Cambridge_ Cambridge University Press, 2010  DOI: 10.1017/CBO9780511976667
  * [2] Miguel Herrero-Collantes and Juan Carlos Garcia-Escartin  “Quantum random number generators”  In _Rev. Mod. Phys._ 89.1, 2017  DOI: 10.1103/RevModPhys.89.015004
  * [3] Thomas Jennewein et al.  “A fast and compact quantum random number generator”  In _Rev. Sci. Instrum._ 71.4, 2000, pp. 1675–1680  DOI: 10.1063/1.1150518
  * [4] André Stefanov et al.  “Optical quantum random number generator”  In _J. Mod. Opt._ 47.4 Taylor & Francis, 2000, pp. 595–598  DOI: 10.1080/09500340008233380
  * [5] P.. Wang, G.. Long and Y.. Li  “Scheme for a quantum random number generator”  In _J. Appl. Phys._ 100.5, 2006, pp. 056107  DOI: 10.1063/1.2338830
  * [6] Makoto Naruse et al.  “Single-photon decision maker”  In _Sci. Rep._ 5.1, 2015, pp. 13253  DOI: 10.1038/srep13253
  * [7] Shafi Goldwasser and Silvio Micali  “Probabilistic encryption”  In _J. Comput. Syst. Sci._ 28.2 Academic Press, 1984, pp. 270–299  DOI: 10.1016/0022-0000(84)90070-9
  * [8] Alvaro Alarcón et al.  “Polarization-independent single-photon switch based on a fiber-optical Sagnac interferometer for quantum communication networks”  In _Opt. Express_ 28.22, 2020, pp. 33731  DOI: 10.1364/oe.408637
  * [9] Costantino Agnesi et al.  “All-fiber self-compensating polarization encoder for quantum key distribution”  In _Opt. Lett._ 44.10, 2019, pp. 2398  DOI: 10.1364/ol.44.002398
  * [10] Marco Avesani et al.  “Deployment-ready quantum key distribution over a classical network infrastructure in Padua”  In _J. Light. Technol._ Institute of ElectricalElectronics Engineers Inc., 2021  DOI: 10.1109/JLT.2021.3130447
  * [11] G L Roberts et al.  “Patterning-effect mitigating intensity modulator for secure decoy-state quantum key distribution”  In _Opt. Lett._ 43.20, 2018, pp. 5110  DOI: 10.1364/ol.43.005110
  * [12] Andrew Rukhin, Juan Soto and James Nechvatal  “A Statistical Test Suite for Random and Pseudorandom Number Generators for Cryptographic Applications”  In _Nist Special Publication_ 22, 2010  DOI: 10.6028/NIST.SP.800-22r1a
  * [13] Robert G. Brown  “Dieharder: A random number test suite”  In _Dieharder: A random number test suite_ , 2011
  * [14] Xiongfeng Ma et al.  “Quantum random number generation”  In _npj Quantum Inf._ 2.1, 2016, pp. 16021  DOI: 10.1038/npjqi.2016.21
  * [15] Jaime Cariñe et al.  “Multi-core fiber integrated multi-port beam splitters for quantum information processing”  In _Optica_ 7.5 OSA, 2020, pp. 542–550  DOI: 10.1364/OPTICA.388912
  * [16] You-Qi Nie et al.  “Experimental measurement-device-independent quantum random-number generation”  In _Phys. Rev. A_ 94 American Physical Society, 2016, pp. 060301  DOI: 10.1103/PhysRevA.94.060301
  * [17] Xiaosong Ma et al.  “A high-speed tunable beam splitter for feed-forward photonic quantum information processing”  In _Opt. Express_ 19.23 OSA, 2011, pp. 22723–22730  DOI: 10.1364/OE.19.022723
  * [18] Michal Mičuda et al.  “Highly stable polarization independent Mach-Zehnder interferometer”  In _Rev. Sci. Instrum._ 85.8 American Institute of PhysicsAIP, 2014, pp. 083103  DOI: 10.1063/1.4891702
  * [19] G.. Xavier and J.. Weid  “Stable single-photon interference in a 1 km fiber-optic Mach–Zehnder interferometer with continuous phase adjustment”  In _Opt. Lett._ 36.10 Optical Society of America, 2011, pp. 1764  DOI: 10.1364/ol.36.001764
  * [20] F Reynaud, J J Alleman and P Connes  “Interferometric control of fiber lengths for a coherent telescope array”  In _Appl. Opt._ 31.19, 1992, pp. 3736  DOI: 10.1364/ao.31.003736
  * [21] Bing Qi  “True randomness from an incoherent source”  In _Rev. Sci. Instrum._ 88.11 American Institute of Physics Inc., 2017, pp. 113101  DOI: 10.1063/1.4986048
  * [22] Xiongfeng Ma et al.  “Postprocessing for quantum random-number generators: Entropy evaluation and randomness extraction”  In _Phys. Rev. A: At. Mol. Opt. Phys._ 87.6, 2013, pp. 62327  DOI: 10.1103/PhysRevA.87.062327
  * [23] H.. Hung, Robert T. O’Neill, Peter Bauer and Karl Kohne  “The Behavior of the P-Value When the Alternative Hypothesis is True”  In _Biometrics_ 53.1 JSTOR, 1997, pp. 11  DOI: 10.2307/2533093


††thanks: Corresponding author<EMAIL_ADDRESS>
# Enhancing material property prediction with ensemble deep graph
convolutional networks

Chowdhury Mohammad Abid Rahman  Lane Department of Computer Science and
Electrical Engineering, West Virginia University  Ghadendra Bhandari
Department of Physics and Astronomy, West Virginia University  Nasser M
Nasrabadi  Lane Department of Computer Science and Electrical Engineering,
West Virginia University  Aldo H. Romero  Department of Physics and Astronomy,
West Virginia University  Prashnna K. Gyawali  Lane Department of Computer
Science and Electrical Engineering, West Virginia University

###### Abstract

Machine learning (ML) models have emerged as powerful tools for accelerating
materials discovery and design by enabling accurate predictions of properties
from compositional and structural data. These capabilities are vital for
developing advanced technologies across fields such as energy, electronics,
and biomedicine, potentially reducing the time and resources needed for new
material exploration and promoting rapid innovation cycles. Recent efforts
have focused on employing advanced ML algorithms, including deep learning -
based graph neural network, for property prediction. Additionally, ensemble
models have proven to enhance the generalizability and robustness of ML and
DL. However, the use of such ensemble strategies in deep graph networks for
material property prediction remains underexplored. Our research provides an
in-depth evaluation of ensemble strategies in deep learning - based graph
neural network, specifically targeting material property prediction tasks. By
testing the Crystal Graph Convolutional Neural Network (CGCNN) and its
multitask version, MT-CGCNN, we demonstrated that ensemble techniques,
especially prediction averaging, substantially improve precision beyond
traditional metrics for key properties like formation energy per atom ($\Delta
E^{f}$), band gap ($E_{g}$) and density ($\rho$) in 33,990 stable inorganic
materials. These findings support the broader application of ensemble methods
to enhance predictive accuracy in the field.

††preprint: APS/123-QED

##  I Introduction

Predicting material crystal properties involves forecasting the chemical and
physical traits of crystalline materials based on their molecular and atomic
structures. This task is vital for fields like electronics, medicine,
aeronautics, and energy storage management [1, 2, 3, 4]. Accurately predicting
these properties from compositional and structural data is instrumental in the
discovery of new materials for advanced technologies. Using computational
methods and data-driven strategies, researchers can efficiently explore and
optimize material designs, avoiding the slow and costly process of
experimental trial-and-error [5]. Understanding how a material’s crystalline
structure impacts its properties requires a blend of computational and
experimental investigations. Although density functional theory (DFT) [6, 7]
is a well-established method, it is often perceived as computationally complex
and time-intensive. Moreover, the quest for materials with specific properties
within an extensive material search space poses challenges, and delays
progress in swiftly evolving domains like medical science, aeronautics, and
energy engineering, where accuracy and speed are paramount.

Machine learning (ML)-based models have emerged as a promising solution to
this challenge. These models can match the accuracy of DFT calculations and
support rapid material discovery, thanks to the growing material databases [8,
9, 10]. By harnessing the capabilities of ML algorithms and the growing
abundance of data in material repositories, these models can effectively
navigate the vast material landscape and pinpoint promising candidates with
desired properties. This data-driven approach for predicting material crystal
properties has shown promise and gained significant attention for its accuracy
and unparalleled speed [11].

However, the intricate atomic arrangements and the intrinsic correlations
between structure and properties present formidable challenges for ML models
to encode pertinent structural information accurately, mainly because of the
number of atoms involved and the internal degrees of freedom within the
crystal structures. Thus, effective representation capturing spatial
relationships and periodic boundary conditions within a crystal lattice
becomes challenging. Furthermore, a model adept with one crystal structure
might falter with another, given crystals’ inherent periodicity and symmetry.
Also, traditional ML models often fail to incorporate the nuanced knowledge of
unit cells and their repetitive nature, an essential aspect of crystallography
[12].

Moreover, for ML models the challenge of representing crystal systems, which
vary widely in size, arises because these models typically require input data
in the form of fixed-length vectors. To address this, researchers have
developed two main approaches. The first involves manually creating fixed-
length feature vectors based on basic material properties [13, 14], which,
though effective, necessitates tailored designs for each property being
predicted. The second approach uses symmetry-invariant transformations [15] of
the atomic coordinates to standardize input data, which, while solving the
issue of variable crystal sizes, complicates model interpretability due to the
complexity of these transformations. Historically, AI-driven material science
has focused on creating custom descriptors for predicting material properties,
utilizing expert knowledge for specific applications. However, these custom
descriptors often lack versatility beyond their initial scope.

To address these challenges in material property prediction, the recent
adoption of deep learning (DL) has shown significant promise [16]. These
networks are adept at learning data distributions as embeddings, serving as
effective feature descriptors for predicting input data characteristics.
Unlike traditional grid-like image representations, crystal structures—with
their node-like atoms and bond-like edges—are ideally represented as graph-
based structures. This has highlighted the suitability of graph neural
networks (GNN) for modeling crystal structures, leveraging their natural
composition and bonding structures. GNNs thus facilitate the automatic
extraction of optimal representations from atoms and bonds by representing
these materials in deep networks.

While current DL models, including advanced GNNs, have successfully integrated
complex structural, geometrical, and topological features for predicting
material properties, a critical aspect often overlooked is the comprehensive
exploration of their training dynamics. Typically, the quest for the lowest
validation loss serves as a proxy for identifying the optimal model. However,
due to the highly non-convex nature of training deep neural networks, this
lowest validation loss does not necessarily correspond to the true optimal
point [17, 18]. Other regions within the loss landscape might capture the
relationship between material structure and properties differently. There has
been limited focus on understanding model behavior beyond the conventional
point of lowest validation loss. This oversight suggests we might not fully
grasp the potential and versatility of these models in capturing structure-to-
property relationships.

In this research, we propose a critical yet often overlooked hypothesis:
optimal model performance may not reside solely at the singular point of
lowest validation loss. Instead, it could be spread across multiple valleys
within the loss terrain. Our objective is to highlight these hidden,
underexplored models that could offer more accurate predictions and a deeper
understanding of material attributes. Furthermore, this approach may reveal
models that provide a better balance between variance and bias, underscoring
the necessity of examining the overall loss landscape to fully understand the
adaptability and efficiency of deep learning models.

Thus, we investigate various regions within the loss landscape where models
still perform satisfactorily. This exploration has yielded insights into
models that are robust and generalize well to new data. We propose combining
these different models to create a unified ensemble model. First, we analyze
the property prediction performance of a prominent GNN-based material
prediction model, CGCNN [19], and its multi-task variant, MT-CGCNN [20],
focusing on three widely studied properties: formation energy per atom,
bandgap, and density. We then examine the impact of ensemble techniques on
model performance, specifically the model average ensemble and prediction
average ensemble methods, to support our hypothesis of identifying influential
models within adjacent areas of the loss landscape. Finally, we conduct an
extensive evaluation to understand the ensemble’s effect on prediction
performance across the spectrum of properties. Overall, our contributions
include:

  1. 1.

We introduce ensemble techniques to improve the material property prediction
capabilities of prominent GNN-based methods, including CGCNN [19], and MT-
CGCNN [20].

  2. 2.

We conduct comprehensive experiments on three widely studied properties:
formation energy per atom, bandgap, and density.

  3. 3.

We assess the impact of ensemble models across a wide spectrum of material
properties, highlighting the effectiveness of ensemble-based approaches in
extreme test conditions.

##  II Related Work

Not having regular grid-like structures such as images or 1-D signals, CNN
could not be the automatic choice for studying material structure from DL
point of view. Rather its irregular structural shape made it a suitable
candidate for graph representation and graph neural network (GNN), where atoms
are considered nodes and atomic bonds edges. Therefore, convolving on the
graph structure which is converted from the actual atomic structure of the
material was a prominent choice for the researchers. Crystal Graph Convolution
Neural Network (CGCNN) [19] and SchNet[21] first proposed this graph
representation and utilized raw features like atom type and atomic bond
distance to predict material properties comparable with DFT computed values.
Although these models performed very well in terms of predicting material’s
properties they also exhibit inevitable and notable challenges due to their
reliance on distance-based message-passing mechanisms [22]. Firstly,
neglecting many-body and directional information can overlook critical aspects
important for understanding material properties. Secondly, the reliance on
nearest neighbors to define graph edges could misrepresent key interactions
due to the ambiguity of chemical bonding. Lastly, the models are limited by
their receptive field, compounded by issues like over-smoothing and over-
squashing, which restrict their ability to account for the long-range or
global influence of structures on properties.

The more recent DL models for material property prediction such as Atomistic
Line Graph Neural Network (ALIGNN) [11], iCGCNN [23],MatErials Graph Network
(MEGNet) [12], Orbital Graph Convolution Neural Network (OGCNN) [24], DimeNet
[25], GemNet [26] and Geometric-information-enhanced Crystal Graph Neural
Network (Geo-CGNN) [27] thus tried to incorporate more geometrical information
like bond angle, orbital interaction, body order information, directional
information, distance vector to outperform previously proposed distance-based
models. Some models proposed attention mechanisms and self supervised learning
(SSL) as well to choose the relative importance of features in predicting
material properties like Matformer [28], Equiformer [29], GATGNN [30],
Crystal-Twins [31] and DSSL [32].

It is a well-known fact that Deep learning models for complicated tasks
navigate a high-dimensional space to minimize a function that quantifies the
‘loss’ or error between the actual data and the expected results. We refer to
this optimization landscape as the ‘loss landscape’ [33, 34, 18]. It is
frequently represented visually as a surface or landscape with hills and
valleys. The loss landscapes of deep neural networks are extremely intricate
and non-convex. This indicates that while the model may converge to numerous
local minima (valleys), not all of them will result in the best solution
(global minimum). The training epoch in which the model achieves the lowest
validation loss is referred to as the ‘best-validated epoch’. Nevertheless,
the model at this epoch may not truly represent the best generalizable model
because of the complexity of the loss landscape and elements like overfitting
and identical loss with differences in functional space [35, 36, 37]. Though
they may have a little larger or the same loss, models from other epochs may
perform better on unknown data or have superior generalization ability for
having differences in loss dynamics. Therefore, our focus in this work is not
on the GNN models or associating features with them to strengthen material
property prediction but on creating a generalized framework of ensemble models
based on the cross-validation loss trajectory that might yield better
generalization of the prediction task with improved accuracy of prediction.

The comprehensive approach of any ensemble technique is to accumulate a set of
models using a function. It has widely been used to strengthen the performance
of ML and DL models. There have been a number of approaches that can be
followed to ensemble the power of various models in DL. For example, diversity
in base classifiers for ensembles is achieved through two distinct approaches,
depending on whether the ensemble is composed of models of the same type
(homogeneous) or different types (heterogeneous) [38]. Different data fusion
methods like max-voting, average voting, weighted average voting, and meta-
learning have been proposed [39]. Also, four types of ensemble techniques are
normally adopted in literature such as bagging, boosting, sorting, and
decision fusion [38]. Different ensemble techniques are already being used to
improve the overall performance of the single-base models in various arenas of
deep learning. Although the ensemble is a highly studied topic, it is mostly
used in fields like speech recognition [40], health-care [41], natural
language processing [42], and computer vision [43]. Compared to these works,
our work is the first to consider the ensemble strategies for material
property prediction tasks, and our ensemble framework is simple and considers
the ensembling by aggregating models across different training stages.

##  III Methodology

In this section, we first provide background details about the GCNN framework
for material property prediction, and then introduce our ensemble framework to
achieve enhanced predictions.

###  III.1 Preliminary: GCNN

Graph Convolutional Neural Networks (GCNNs) have emerged as a powerful tool in
materials science, enabling researchers to analyze and predict the properties
and behaviors of materials in a novel and efficient manner. Unlike traditional
convolutional neural networks that process grid-like data (e.g., images),
GCNNs are designed to handle graph-structured data, which is intrinsic to the
representation of atomic and molecular structures in materials science. These
models exploit the graph structure of materials, where nodes can represent
atoms and edges can denote chemical bonds or spatial relationships. By doing
so, GCNNs can capture both the local and global structural information of
materials.

Figure 1: Overview of GCNN in material property prediction tasks. Initially,
the crystal structure is created from information files. Then, the crystal
graph is constructed from the structure. GCN layers and pooling layers are
utilized to obtain crystal embeddings, after which fully connected (FC) layers
are employed to predict properties. For a single-task head, a single FC layer
is used, while for the prediction of multiple properties, multiple FC layers
are utilized.

Crystal Graph Convolutional Neural Networks (CGCNN) [19] and SchNet [21]
represent the two most prominent graph neural network architectures tailored
for material science applications. These models refine atom representations
within a structure by considering the types and bond lengths of neighboring
atoms. Subsequently, they aggregate these updated atom-level representations
to form a comprehensive representation of the entire structure. In CGCNN, the
crystal architecture is represented as a graph that accommodates the details
of atoms and their bonds with neighbors, and a graph convolution network is
constructed on such graph to attain the representations useful to material
property prediction. The architecture can comprise of single-task head [19] or
multi-task [20] head depending on the application. We present the overall
diagram combining the single- and multi-task setup in Fig. 1. In the presented
network, the atom feature encoding vector can be noted as $v_{i}$ where i is
an atom or node in the graph. Similarly, an edge or atomic bond among
neighbors is denoted as $(i,j)_{k}$ representing $k^{th}$ bond between atoms
$i$ and $j$, and its feature vector is $u_{(i,j)_{k}}$. The goal here is to
update the atom feature vector $v_{i}$ by iterative convolution with neighbors
and bonds (1) and generate a comprehensive feature vector for the crystal
structure by pooling (2).

$v^{(t+1)}_{i}=\text{Conv}\left(v^{(t)}_{i},v^{(t)}_{j},u_{(i,j)k}\right),\quad(i,j)\in\mathcal{G}.$
(1)
$v_{c}=\text{Pool}\left(v^{(0)}_{0},v^{(0)}_{1},\ldots,v^{(0)}_{N},\ldots,v^{(R)}_{N}\right)$
(2) $\min_{\mathbf{W}}J(y,f(\mathbf{C};\mathbf{W}))$ (3)

The training procedure involves minimizing a cost function $J(y;\hat{y})$
where $y$ is the DFT computed value and $\hat{y}$ is the prediction of the
model. Therefore, CGCNN can be considered a function that tries to approximate
the actual property value $y$ by mapping a crystal $C$ employing weights $W$
as shown in Eq. (3).

###  III.2 Ensemble models

In this section, we present an ensemble framework to enhance the results
obtained from GCNN networks. Our central hypothesis is that in the training of
deep neural networks, selecting a single model based solely on the lowest
validation error may not always guarantee the most comprehensive learning of
all features within the dataset. This limitation can be attributed to the
highly non-convex nature of the loss landscape that characterizes neural
network optimization [44]. In such a complex terrain, the path to minimizing
loss involves numerous local minima and saddle points. This suggests that
multiple models—each residing in different areas of the loss landscape—could
perform similarly well on the validation set, albeit potentially capturing
different aspects of the data [17].

Toward this, we propose ensemble framework for achieving ensemble models from
the training dynamics of GCNNs, that capitalizes on the temporal evolution of
model parameters across training epochs. Consider the training process to span
a fixed number of epochs, $T$, during which each epoch yields a candidate
model characterized by unique properties and attributes reflective of its
learning state at that point in time. Let, $f(\mathbf{x};\Theta_{t})$
represent the GCNN model at some epoch $t$. Here $x$ is the input data and
$\Theta_{t}$ represents the model’s parameters (weight and biases) at epoch
$t$. By running the training for $T$ epochs, we generate the sequence of
models $\\{f(x;\Theta_{1}),f(x;\Theta_{2}),\ldots,f(x;\Theta_{T})\\}$, and for
each model $f(\mathbf{x};\Theta_{t})$, we compute the MSE on a given
validation set as:

$\text{MSE}_{t}=\frac{1}{N_{\text{val}}}\sum_{i=1}^{N_{\text{val}}}\left(y_{i}-f(\mathbf{x}_{i};\Theta_{t})\right)^{2}$
(4)

where $N_{val}$ is the number of sample in the validation set, $y_{i}$ is the
true value and $f(x_{i};\Theta_{t})$ is the predicted value for the
$i^{\text{th}}$ sample in the validation set by the model at epoch $t$. Using
the metric from Eq. 4, we sort the models based on their MSE, selecting the
top $n$ models with the lowest MSE values. Let
$\Theta^{(1)},\Theta^{(2)},\ldots,\Theta^{(n)}$ be the parameters of these
top-$n$ models after sorting, and we present two different strategies for
aggregating these $n$ models.

####  III.2.1 Prediction-based ensemble modeling

For prediction-based ensemble modeling, we first calculate the prediction for
each $\Theta^{(t)}$ model, within the top-$n$ as:

$\hat{y}_{t}=f(x;\Theta_{t})$ (5)

and create prediction ensemble as:

$\bar{y}_{\text{prediction-ensemble}}=\frac{1}{n}\sum_{i=1}^{n}\hat{y}_{n}$
(6)

####  III.2.2 Model-based ensemble modeling

Here, we first aggregate top-$n$ models together by creating an ensemble model
$\Theta_{avg}$:

$\Theta_{\text{avg}}=\frac{1}{n}\sum_{j=1}^{n}\Theta^{(j)}$ (7)

The final prediction of the ensemble model for a new input $x$ is then given
as:

$\hat{y}_{\text{model-ensemble}}=f(x;\Theta_{\text{avg}})$ (8)

We present the overall schematic of our proposed ensemble-based framework in
Fig. 2. On the left, we illustrate the prediction ensemble, and on the right
panel, the model ensemble framework is depicted.

Figure 2: Overview of ensemble strategies: (a) prediction averaging ensemble
technique and (b) the model averaging ensemble technique.

##  IV Experiment

In this section, we discuss the data we have used for experiments, the
implications of choosing the properties we worked with, the outline of the
experimental setup, and the results with their interpretation, significance,
and implications.

###  IV.1 Data

All the models in this study were trained on the dataset from the Materials
Project [8]. The Materials Project is a multi-institutional, multinational
effort to compute the properties of all inorganic materials and provide the
data and associated analysis algorithms to every materials researcher free of
charge. One of the largest and most popular three-dimensional (3D) materials
datasets in the materials science field, the Materials Project collection
covers approximately 155,000 materials, including a wide variety of inorganic
compounds. This broad coverage guarantees a representative and varied sampling
of material kinds, extending the generalizability of our results. Moreover,
the MP dataset has been successfully used in numerous studies to create and
evaluate predictive models for a range of material properties. The
trustworthiness of the dataset and the accomplishments of earlier studies
employing MP data highlight its suitability for verifying our methods. For
this study, we worked with 33,990 stable materials, which refers to the set of
materials that, under standard conditions, have a low energy state and are
hence thermodynamically favorable to exist in their current form. An indicator
of stability is the ‘energy above hull’ metric, which shows how much energy
would be released by each atom in the event that the material changes into the
most stable phase combination. For several reasons, it is critical for this
work to concentrate on stable materials. Stable materials are more important
for real-world applications because they are thermodynamically favored,
meaning they occur naturally or can be created with less energy. Focusing on
these materials allows the research to directly target materials that are
useful in energy storage, electronics, and catalysis, among other real-world
applications. To demonstrate the efficacy of the proposed ensemble framework,
we focused on three distinct properties: Formation energy per atom ($\Delta
E^{f}$), Density ($\rho$), and Bandgap ($E_{g}$). These three properties are
important because formation energy determines the thermodynamic stability of a
material, bandgap signifies whether a material is an insulator, conductor, or
semiconductor, and density defines the stiffness.

###  IV.2 Experimental setup

For exploring the efficacy of our proposed ensemble model, we primarily
consider CGCNN [19] and its multi-task extension, MT-CGCNN [20], as base
models, and apply our proposed ensembling to construct the ensemble framework.
Furthermore, we vary the number of convolutional layers within CGCNN and MT-
CGCNN to create two separate versions of the network for each category.
Throughout our experiments and results, we refer to them as
CGCNN${}_{\text{1}}$ (number of convolutional layers, $n{\text{c}}$ = 3),
CGCNN${}_{\text{2}}$ ($n{\text{c}}$ = 5), MT-CGCNN${}_{\text{1}}$
($n{\text{c}}$ = 3), and MT-CGCNN${}_{\text{2}}$ ($n{\text{c}}$ = 5). For the
MT-CGCNN models, our multi-task objective involved predicting all three
properties together using three separate heads, and optimization weights of
1.5 were used for $\Delta E^{f}$, 3 for $E_{g}$, and 1.5 for $\rho$. Although
we present both prediction and model ensemble frameworks, we found the
prediction ensemble to yield the best results. Therefore, we utilize the
prediction ensemble for all our analyses and comparisons against the Best-val
model. However, it is important to note that we also compare the performance
between prediction and model-based ensembles to establish the efficacy of the
prediction average ensemble over the latter.

Across all models, the length of the atom feature vector was set to 64, and
the hidden feature vector to 128. For all experiments, MSE was used as the
loss function, SGD as the optimizer, and a fixed learning rate of 0.01 was
applied. We utilized an NVIDIA GeForce GTX TITAN X graphics processing unit
(GPU) with 12 GB of memory for all model training and evaluation tasks. The
training, validation, and test data were randomly selected from 33990 data
points and were kept consistent across all experiments, employing a 70-10-20
distribution strategy, and the batch size was uniformly set at 256. This
random selection ensures that our model is not biased towards any specific
subset of data, thereby enhancing its ability to generalize to unseen data.
Our sampling strategy did not follow any particular distribution, ensuring
that the training, validation, and test sets were representative of the
overall dataset. By not constraining the sampling process to a specific
distribution, we mitigate the risk of overfitting to particular patterns in
the training data, thus improving the model’s robustness. All models were run
for 100 epochs, and up to 50 models were considered for the ensemble. To
determine the best model among the baseline models, and to select several
models for creating the ensemble, MSE loss on validation data was used.
However, for reporting results, we used MAE on the test dataset to follow
standard practice in the literature.

Table 1: Comparison of Best-val and Ensemble for CGCNN1,2 (20 models) | CGCNN${}_{\text{1}}$ | CGCNN${}_{\text{2}}$ |
---|---|---|---
| Best-val | Ensemble | Best-val | Ensemble |
Formation Energy ($\Delta E^{f}$) | 0.058 |  0.054 (6.90% $\uparrow$) | 0.060 |  0.055 (8.33% $\uparrow$) |
BG ($E_{g}$) | 0.322 |  0.301 (6.90% $\uparrow$) | 0.293 |  0.278 (5.12% $\uparrow$) |
Density ($\rho$) | 0.134 |  0.128 (4.47% $\uparrow$)) |  0.145 | 0.146 (0.69% $\downarrow$) |
Table 2: Comparison of Best-val and Ensemble for MT-CGCNN${}_{\text{1},\text{2}}$ (40 models) | MT-CGCNN${}_{\text{1}}$ | MT-CGCNN${}_{\text{2}}$ |
---|---|---|---
| Best-val | Ensemble | Best-val | Ensemble |
Formation Energy ($\Delta E^{f}$) | 0.082 |  0.073 (11% $\uparrow$) | 0.081 |  0.076 (6.17% $\uparrow$) |
BG ($E_{g}$) | 0.316 |  0.307 (2.85% $\uparrow$) | 0.294 |  0.280 (4.76% $\uparrow$) |
Density ($\rho$) | 0.216 |  0.192 (11.11% $\uparrow$) | 0.205 |  0.190 (7.32% $\uparrow$) |
Figure 3: Results for Prediction Ensemble on Single-task CGCNN and Mulit-task
CGCNN.

##  V Result

We present our main results in Table 1 and Table 2. For results in both
tables, we present outcomes from the prediction ensemble for a fixed number of
epochs (20 for single-task and 40 for multi-task) as determined by validation
performance. In Table 1, we compare the ensemble framework (prediction
ensemble) against CGCNN${}_{\text{1}}$ and CGCNN${}_{\text{2}}$ across all
three properties. We found that, out of six different settings, our proposed
ensemble framework achieved better results in five settings. Moreover, in one
of the settings where the standard approach achieved a better result, the gap
between our results and the baseline was the smallest (0.69%) compared to
other margins (4% to 8%). In Table 2, we observed enhanced results by our
proposed framework in all six different settings, with the improvement margin
up to 11%.

Figure 4: Prediction ensemble CGCNN${}_{\text{1}}$ across different groups of
data distribution for three properties.

In Fig. 3, we analyze the effect of the prediction ensemble across a different
number of models. The leftmost point, representing the number of models used
to create an ensemble as one, is the best validation model and represents the
standard practice of validation using a single model. Compared to that single
point, every result to the right represents our proposed approach of using an
ensemble-based framework. As seen in both single-task (left panel) and multi-
task (right panel) and across all three properties, the effect of ensembling
for enhancing property prediction is quite evident. It should be noted that in
some cases, such as the formation energy per atom for the single-task model,
the behavior of the ensemble framework appears not to be as effective after a
certain number, but we can still see that the ensemble framework is always
better or similar to the best validation model. This also underscores the
importance of the ensemble framework in achieving enhanced prediction results.

Figure 5: Prediction ensemble MT-CGCNN${}_{\text{1}}$ across different groups
of data distribution for three properties.

We also explore the impact of the ensemble-based framework across different
regions of the material property spectrum. This analysis is crucial for a
comprehensive understanding of the various properties involved in our study.
For example, specific regions of the bandgap determine a material’s
suitability as a conductor, semiconductor, or insulator. Grasping the range of
formation energy per atom is vital, as it reflects the stability and
synthesizability of materials. Materials that are thermodynamically stable,
and therefore more likely to occur in nature or be successfully synthesized in
the laboratory, are distinguished by lower formation energies. Additionally,
pinpointing the extremes within density ranges aids in assessing the
durability of high-density materials or the practicality of lightweight
insulators with low density.

We present the results in Fig.4 and Fig.5, where we partitioned the test data
from the 10th to the 90th percentile in both top-bottom (left panel) and
bottom-top (right panel) distributions. This approach helps to identify
performance differences across different regions of the property spectrum for
CGCNN${}_{\text{1}}$ and MT-CGCNN${{}_{\text{1}}}$ models. It is observed that
the ensemble model, in all instances, aligns with the trend of the original
single best model, exhibiting a reduced MAE value or improved accuracy across
all percentiles of data distribution in most scenarios. Notably, significant
improvement is observed in certain regions, for example, for the bottom 10% of
bandgap materials in multi-task models.

Finally, although our experiment clearly demonstrated the benefit of a
prediction ensemble over a model ensemble, in Fig. 6, we include an analysis
in this paper that compares the performance between a prediction ensemble and
a model ensemble for both single-task and multi-task frameworks for band gap.
As shown in the figure, although the performances of both ensemble approaches
appear similar in the multi-task settings (right panel), for the single-task
model, we observe that the model ensemble resulted in the worst performance,
even when compared to the single best validation model.

Figure 6: Result comparison for Prediction Ensemble and Model average Ensemble
on Single-task CGCNN1 (a) and MT-CGCNN1 (b) for Band gap.

##  VI Conclusion

In this paper, we explore the impact of an ensemble framework on the task of
material prediction. Our proposed framework is both simple and effective,
leveraging the loss landscape of deep neural network training without
requiring computationally expensive ensemble strategies. We present two types
of ensembles: prediction-based and model-based. The former involves
aggregating predictions across different models to form the ensemble model,
while the latter aggregates the models first, then generates a single
prediction as the ensemble result. Our analysis shows that the prediction
ensemble consistently outperforms the model ensemble. As a result, we
conducted a comprehensive analysis across various models (single-task vs.
multi-task), architectures (variations in GCNN depth), and properties, finding
that the prediction ensemble almost always improves the predictive performance
over the single best validation model. We also examined the efficacy of our
proposed framework across different property spectra for test data
distribution. Overall, our extensive analysis demonstrates the robustness of
the proposed framework in generating enhanced predictive results. Future work
will focus on investigating systematic approaches for calculating the number
of models to select candidate models for the ensemble framework, extending the
analysis to include other GNN models, and on expanding our results to 2D
datasets.

###### Acknowledgements.

This research is partially supported by West Virginia Higher Education Policy
Commission’s Research Challenge Grant Program 2023 and CITeR NSF/IUCRC center
(23F-02W).

## References

  * Shen  _et al._ [2022] L. Shen, J. Zhou, T. Yang, M. Yang, and Y. P. Feng, High-throughput computational discovery and intelligent design of two-dimensional functional materials for various applications, Acc. Mater. Res. 3, 572 (2022).
  * Shahzad  _et al._ [2024] K. Shahzad, A. I. Mardare, and A. W. Hassel, Accelerating materials discovery: combinatorial synthesis, high-throughput characterization, and computational advances, Science and Technology of Advanced Materials: Methods 4, 2292486 (2024).
  * Wei  _et al._ [2019] J. Wei, X. Chu, X.-Y. Sun, K. Xu, H.-X. Deng, J. Chen, Z. Wei, and M. Lei, Machine learning in materials science, InfoMat 1, 338 (2019).
  * Liu  _et al._ [2017] Y. Liu, T. Zhao, W. Ju, and S. Shi, Materials discovery and design using machine learning, Journal of Materiomics 3, 159 (2017).
  * Kauwe  _et al._ [2020] S. K. Kauwe, T. Welker, and T. D. Sparks, Extracting knowledge from dft: experimental band gap predictions through ensemble learning, Integr Mater Manuf Innov 9, 213 (2020).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Self-consistent equations including exchange and correlation effects, Physical review 140, A1133 (1965).
  * Hohenberg and Kohn [1964] P. Hohenberg and W. Kohn, Inhomogeneous electron gas, Physical review 136, B864 (1964).
  * Jain  _et al._ [2013] A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, _et al._ , Commentary: The materials project: A materials genome approach to accelerating materials innovation, APL materials 1, 10.1063/1.4812323 (2013).
  * Choudhary  _et al._ [2020] K. Choudhary, K. F. Garrity, A. C. Reid, B. DeCost, A. J. Biacchi, A. R. Hight Walker, Z. Trautt, J. Hattrick-Simpers, A. G. Kusne, A. Centrone, _et al._ , The joint automated repository for various integrated simulations (jarvis) for data-driven materials design, npj Comput. Mater. 6, 173 (2020).
  * Kirklin  _et al._ [2015] S. Kirklin, J. E. Saal, B. Meredig, A. Thompson, J. W. Doak, M. Aykol, S. Rühl, and C. Wolverton, The open quantum materials database (oqmd): assessing the accuracy of dft formation energies, npj Comput. Mater. 1, 1 (2015).
  * Choudhary and DeCost [2021] K. Choudhary and B. DeCost, Atomistic line graph neural network for improved materials property predictions, npj Comput. Mater. 7, 185 (2021).
  * Chen  _et al._ [2019] C. Chen, W. Ye, Y. Zuo, C. Zheng, and S. P. Ong, Graph networks as a universal machine learning framework for molecules and crystals, Chem. Mater. 31, 3564 (2019).
  * Isayev  _et al._ [2017] O. Isayev, C. Oses, C. Toher, E. Gossett, S. Curtarolo, and A. Tropsha, Universal fragment descriptors for predicting properties of inorganic crystals, Nat. Commun. 8, 15679 (2017).
  * Xue  _et al._ [2016] D. Xue, P. V. Balachandran, J. Hogden, J. Theiler, D. Xue, and T. Lookman, Accelerated search for materials with targeted properties by adaptive design, Nat. Commun. 7, 1 (2016).
  * Seko  _et al._ [2017] A. Seko, H. Hayashi, K. Nakayama, A. Takahashi, and I. Tanaka, Representation of compounds for machine-learning prediction of physical properties, Phys. Rev. B 95, 144110 (2017).
  * Cao  _et al._ [2019] Z. Cao, Y. Dan, Z. Xiong, C. Niu, X. Li, S. Qian, and J. Hu, Convolutional neural networks for crystal material property prediction using hybrid orbital-field matrix and magpie descriptors, Crystals 9, 191 (2019).
  * Gyawali  _et al._ [2022] P. K. Gyawali, X. Liu, J. Zou, and Z. He, Ensembling improves stability and power of feature selection for deep learning models, in  _Machine Learning in Computational Biology_ (PMLR, 2022) pp. 33–45.
  * Garipov  _et al._ [2018] T. Garipov, P. Izmailov, D. Podoprikhin, D. P. Vetrov, and A. G. Wilson, Loss surfaces, mode connectivity, and fast ensembling of dnns, Advances in neural information processing systems 31 (2018).
  * Xie and Grossman [2018] T. Xie and J. C. Grossman, Crystal graph convolutional neural networks for an accurate and interpretable prediction of material properties, Phys. Rev. Lett. 120, 145301 (2018).
  * Sanyal  _et al._ [2018] S. Sanyal, J. Balachandran, N. Yadati, A. Kumar, P. Rajagopalan, S. Sanyal, and P. Talukdar, Mt-cgcnn: Integrating crystal graph convolutional neural network with multitask learning for material property prediction, arXiv preprint arXiv:1811.05660 10.48550/arXiv.1811.05660 (2018).
  * Schütt  _et al._ [2018] K. T. Schütt, H. E. Sauceda, P.-J. Kindermans, A. Tkatchenko, and K.-R. Müller, Schnet–a deep learning architecture for molecules and materials, J. Chem. Phys. 148, 10.1063/1.5019779 (2018).
  * Gong  _et al._ [2023] S. Gong, K. Yan, T. Xie, Y. Shao-Horn, R. Gomez-Bombarelli, S. Ji, and J. C. Grossman, Examining graph neural networks for crystal structures: limitations and opportunities for capturing periodicity, Sci. Adv. 9, eadi3245 (2023).
  * Park and Wolverton [2020] C. W. Park and C. Wolverton, Developing an improved crystal graph convolutional neural network framework for accelerated materials discovery, Phys. Rev. Mater. 4, 063801 (2020).
  * Karamad  _et al._ [2020] M. Karamad, R. Magar, Y. Shi, S. Siahrostami, I. D. Gates, and A. B. Farimani, Orbital graph convolutional neural network for material property prediction, Phys. Rev. Mater. 4, 093801 (2020).
  * Gasteiger  _et al._ [2020] J. Gasteiger, J. Groß, and S. Günnemann, Directional message passing for molecular graphs, arXiv preprint arXiv:2003.03123 10.48550/arXiv.2003.03123 (2020).
  * Gasteiger  _et al._ [2021] J. Gasteiger, F. Becker, and S. Günnemann, Gemnet: Universal directional graph neural networks for molecules, Advances in Neural Information Processing Systems 34, 6790 (2021).
  * Cheng  _et al._ [2021] J. Cheng, C. Zhang, and L. Dong, A geometric-information-enhanced crystal graph network for predicting properties of materials, Commun. Mater. 2, 92 (2021).
  * Yan  _et al._ [2022] K. Yan, Y. Liu, Y. Lin, and S. Ji, Periodic graph transformers for crystal material property prediction, Advances in Neural Information Processing Systems 35, 15066 (2022).
  * Liao and Smidt [2022] Y.-L. Liao and T. Smidt, Equiformer: Equivariant graph attention transformer for 3d atomistic graphs, arXiv preprint arXiv:2206.11990 10.48550/arXiv.2206.11990 (2022).
  * Louis  _et al._ [2020] S.-Y. Louis, Y. Zhao, A. Nasiri, X. Wang, Y. Song, F. Liu, and J. Hu, Graph convolutional neural networks with global attention for improved materials property prediction, Phys. Chem. Chem. Phys. 22, 18141 (2020).
  * Magar  _et al._ [2022] R. Magar, Y. Wang, and A. Barati Farimani, Crystal twins: self-supervised learning for crystalline material property prediction, NPJ. Comput. Mater. 8, 231 (2022).
  * Fu  _et al._ [2024] N. Fu, L. Wei, and J. Hu, Physics-guided dual self-supervised learning for structure-based material property prediction, J. Phys. Chem. Lett. 15, 2841 (2024).
  * Li  _et al._ [2018] H. Li, Z. Xu, G. Taylor, C. Studer, and T. Goldstein, Visualizing the loss landscape of neural nets, Advances in neural information processing systems 31, (2018).
  * Goodfellow  _et al._ [2016] I. Goodfellow, Y. Bengio, and A. Courville,  _Deep learning_ (MIT press, 2016).
  * Draxler  _et al._ [2018] F. Draxler, K. Veschgini, M. Salmhofer, and F. Hamprecht, Essentially no barriers in neural network energy landscape, in  _International conference on machine learning_ (PMLR, 2018) pp. 1309–1318.
  * Fort and Jastrzebski [2019] S. Fort and S. Jastrzebski, Large scale structure of neural network loss landscapes, Advances in Neural Information Processing Systems 32 (2019).
  * Fort  _et al._ [2019] S. Fort, H. Hu, and B. Lakshminarayanan, Deep ensembles: A loss landscape perspective, arXiv preprint arXiv:1912.02757 https://doi.org/10.48550/arXiv.1912.02757 (2019).
  * Ganaie  _et al._ [2022] M. A. Ganaie, M. Hu, A. K. Malik, M. Tanveer, and P. N. Suganthan, Ensemble deep learning: A review, Eng. Appl. Artif. Intell. 115, 105151 (2022).
  * Mohammed and Kora [2023] A. Mohammed and R. Kora, A comprehensive review on ensemble deep learning: Opportunities and challenges, J. King Saud Univ. Comput. Inf. Sci. 35, 757 (2023).
  * Li  _et al._ [2017] S. Li, X. Lu, S. Sakai, M. Mimura, and T. Kawahara, Semi-supervised ensemble dnn acoustic model training, in  _2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ (IEEE, 2017) pp. 5270–5274.
  * Tanveer  _et al._ [2021] M. Tanveer, A. H. Rashid, M. Ganaie, M. Reza, I. Razzak, and K.-L. Hua, Classification of alzheimer’s disease using ensemble of deep neural networks trained through transfer learning, IEEE J. Biomed. Health. Inf. 26, 1453 (2021).
  * Elnagar  _et al._ [2020] A. Elnagar, R. Al-Debsi, and O. Einea, Arabic text classification using deep learning models, Inf. Process. Manage. 57, 102121 (2020).
  * Roshan  _et al._ [2024] S. Roshan, J. Tanha, M. Zarrin, A. F. Babaei, H. Nikkhah, and Z. Jafari, A deep ensemble medical image segmentation with novel sampling method and loss function, Comput. Biol. Med. , 108305 (2024).
  * Cooper [2018] Y. Cooper, The loss landscape of overparameterized neural networks, arXiv preprint arXiv:1804.10200 https://doi.org/10.48550/arXiv.1804.10200 (2018).