100100 2022

[1]Lißner Julian

[1]Data Analytis in Engineering, University of Stuttgart, Universitätsstraße
32, Stuttgart, 70569, Baden Württemberg, Germany

# Hybrid machine-learned homogenization: Bayesian data mining and
convolutional neural networks

<EMAIL_ADDRESS>Fritzen Felix<EMAIL_ADDRESS>*

###### Abstract

Beyond the generally deployed features for microstructure property prediction
this study aims to improve the machine learned prediction by developing novel
feature descriptors. Therefore, Bayesian infused data mining is conducted to
acquire samples containing characteristics inexplicable to the current feature
set, and suitable feature descriptors to describe these characteristics are
proposed. The iterative development of feature descriptors resulted in 37
novel features, being able to reduce the prediction error by roughly one
third. To further improve the predictive model, convolutional neural networks
(Conv Nets) are deployed to generate auxiliary features in a supervised
machine learning manner. The Conv Nets were able to outperform the feature
based approach. A key ingredient for that is a newly proposed data
augmentation scheme and the development of so-called deep inception modules. A
combination of the feature based approach and the convolutional neural network
leads to a hybrid neural network: A parallel deployment of the both neural
network archetypes in a single model achieved a relative rooted mean squared
error below 1%, more than halving the error compared to prior models operating
on the same data. The hybrid neural network was found powerful enough to be
extended to predict variable material parameters, from a low to high phase
contrast, while allowing for arbitrary microstructure geometry at the same
time.

###### keywords:

microstructure homogenization, convolutional neural networks, feature
engineering, bayesian neural networks, machine learning

##  1 Introduction

High performance materials are of great interest to the industry due to their
capabilities and scope of application, e.g., in aerospace applications or for
batteries, even though their development and manufacturing process is a highly
challenging task. The tailoring of specific materials can be conducted by,
e.g., tuning its microstructure to optimize its properties given specific
requirements. The development process can be significantly boosted by
replacing experimental tests through simulations, which are carried out by
taking the microscopic geometric information of the materials into account
(miehe2002, ; beyerlein2008dislocation, ). Some simulation methods are able to
directly operate on the 3D image representation of the microstructure, e.g.,
obtained from a computed tomography (CT) scan, recently achieving improvements
with respect to computational speed keshav2022fft . The high resolution of the
image representation renders even such efficient methods infeasible to
evaluate in a many query context by investigating numerous (e.g., in-silico
generated) microstructured materials, when optimizing for specific material
behaviour.

Figure 1: The general outline of this paper is presented in figure format. On
the left, the iterative approach of Bayesian assisted data mining is
displayed, and the newly developed features thereof are used in a hybrid
neural network (on the right). The backbone of the hybrid neural network are
deep inception modules in parralel to a bypass using the volume fraction and
the feature regressor.

Machine learning is a suitable tool to further reduce the computational cost
for material development: It can discover complex relationships by studying
the available data and replacing the costly simulations with computationally
affordable operations. Various machine learning algorithms are regularly
deployed in two actively researched fields of microstructure modeling, namely
microstructure reconstruction and microstructure property prediction. In
microstructure reconstruction/synthesis, the topology of an original
microstructured image is adjusted to optimize for selected material behaviour
under constraints (kumar2016markov, ; cang2018improving, ;
seibert2022descriptor, ), tailoring for specific microstructures to serve a
particular purpose. The field of microstructure property prediction is often
applied to a broader spectrum of microstructured materials lissner2019data ;
ford2021machine . It targets the efficient and accurate prediction of the
behaviour of inhomogeneous materials via, e.g., artificial neural networks
(ANN).

The present paper falls into a subclass of the latter category, specifically
microstructure homogenization, which directly predicts effective material
properties from given microstructural image data. Current state of the art
methods often deploy a blend of different unsupervised and supervised machine
learning methods, where the Principal Component Analysis (PCA) is used in
combination of the 2-Point Correlation Function (2PCF) (brough2017materials,
), using the principal scores as input for a supervised machine learned
regressor to conduct the microstructure property linkage, ranging from
polynomial regression (marshall2021autonomous, ) to artificial neural networks
(lissner2019data, ; farizhandi2022processing, ). Lately, the usage of
convolutional neural networks (Conv Nets) has gained popularity in both
outlined research fields (lubbers2017inferring, ; gayon2020pores, ). Conv Nets
are used to directly predict the effective material response
(liu2022correlation, ), or even in combination with the PCA where the reduced
representation of the target values is used to predict stress strain
relationships (yang2020prediction, ). Conv Nets have the advantage over
classical regressors in the sense that they are better suited to an extension
for the prediction of full field solution in voxel (3D pixels) representation,
since the data is often given in image representation, for the microstructured
material as well as the full field response.

The previous study of the authors (lissner2019data, ) has deployed a method
akin to the PCA and used the derived features in various state of the art
machine learned regressors, finding that the accuracy of the approach is
limited, which is confirmed in different studies, e.g., (ford2021machine, ;
fernandez2020generation, ). Improvements beyond the 2PCF have been attempted
by considering partial higher order correlation functions, e.g.,
(fast2011formulation, ). However, the feature identification and computation
becomes increasingly costly while yielding rather limited improvements with
respect to accuracy. In the search for a better prediction, the idea of
auxiliary features came to mind, which led us to data mining (mikut2011data,
). This constitutes the first block of this work (section 2.2): We
systematically categorize the underlying data by using Bayesian neural
networks (tipping2003bayesian, ). We evaluate the aleatoric uncertainty which
can hint at a lack of feature knowledge for certain samples. These samples can
then be examined systematically in order to engineer additional features while
being machine-guided during the feature engineering process. The thereby
identified recurrent characteristics across multiple samples were then
quantified by novel feature descriptors, while keeping computational
efficiency and physical interpretability in consideration.

The second block of this study considers Conv Nets (section 2.3), which are
able to derive machine learned features via supervised learning. Further
improvements with respect to Conv Nets are found by building upon the so
called inception modules (szegedy2015going, ) and, therefrom, developing a
deep inception module. The latter is explicitely designed to be able to
capture features at different length scales within a single microstructural
image. Ultimately, we combine the handcrafted, machine-guided features with
features derived from Conv Nets into a hybrid neural network. The proposed
model more than halved the prediction error compared to previous studies
(lissner2019data, ). This hybrid model is further extended to predict the
effective material properties for variable material parameters, allowing for
variable phase contrast ranging from $\frac{1}{2}$ up to $\frac{1}{100}$ while
considering variable microstructure material characteristics, e.g., with the
volume fraction ranging from 20-80%.

The data darus1151  used for development, training and validation purposes, as
well as the python code codebase  are made freely available.

##  2 Data and Methods

###  2.1 Data overview

Figure 2: The investigated data is shown, where a few RVE images of the
training set are shown on the left. The gray foreground represents the
inclusion phase, and the black background represents the matrix phase. On the
right, one component of the target values is shown for all samples out of the
training set.

In the machine learning focused manuscript, some emphasis has to be given to
the data. This manuscript deals with the prediction of material properties
from bi-phasic microstructured materials. The simplification of a
Representative Volume Element (RVE) is introduced, where it is assumed that a
single frame of the microstructure suffices to characterize the material
behaviour of the macroscopic material torquato2013random  while using periodic
boundary conditions. A compact overview of the data is given in Fig. 2, where
some exemplary images of the microstructure are plotted on the left. These
images, i.e., RVEs of the microstructured material serve as the input data to
our algorithm and will be denoted $\underline{\underline{\rm RVE}}$ in the
following. The target values of the machine learning algorithm are the
components of the effective heat conduction tensor
$\underline{\underline{\bar{\kappa}}}$. One Component of
$\underline{\underline{\bar{\kappa}}}$ is plotted in Fig. 2 b), considering
low conducting inclusions with a phase contrast of $R=5$. A more detailed
explanation of the data is given in the appendix A. Generally, the symmetric
effective heat conduction tensor will be given in de-dimensionalized Mandel
notation

$\bar{\underline{\kappa}}=\begin{bmatrix}\bar{\kappa}_{11}\\\
\bar{\kappa}_{22}\\\ \sqrt{2}\bar{\kappa}_{12}\end{bmatrix}\,.$ (1)

The objective of this paper is to find a machine learned model which is able
to accurately predict the effective material properties of the presented
microstructures, i.e.

$\bar{\underline{\kappa}}=f({\rm RVE})$ (2)

where the machine learned model $f(\cdot)$ directly operates on the image data
of the microstructured material (${\rm RVE}$) or uses features $\underline{x}$
directly extracted from the image. The constraint to our algorithm is to
solely rely on the given images of the RVEs, without any further information
such as, e.g., the number or shape of inclusions. In this manuscript we first
aim to find optimal features $\underline{x}$ and a suitable artificial neural
network $f(\cdot)$ accurately mapping the complex image based relationship
between features and the sought-after outputs.

###  2.2 Bayesian Modeling via artificial neural networks

####  2.2.1 Aleatoric and epistemic uncertainty

In supervised machine learning the training and testing data is always given
with deterministic target values such that error measures are computable.
However, during inference, reliable error measures are generally not
accessible. Consequently, having access to a measure of model confidence
indicating the likelihood for low or high prediction errors during inference
is advantageous.

Bayesian modeling approaches uncertainty quantification via aleatoric and
epistemic uncertainties (kendall2017uncertainties, ; depeweg2019modeling, ).
The epistemic uncertainty aims to represent the uncertainty arising based on
the lack of knowledge about the mapping from input to output values. In the
machine learning context it arises from the model being unable to recover the
original function, if it even exists.

The aleatoric uncertainty arises from the uncertainty in the system or data
and aims to represent noise in the input data, i.e., different output values
for similar or even identical input values. Considering, for instance Fig. 2
b) from a machine learning viewpoint, the variation of the material response
could be interpreted as noisy when considering only the volume fraction as
input. These variations, however, are not due to noise, but arise due to
different phenomena inexplicable when trying to predict the response
exclusively by the materials volume fraction. Consequently, the aleatoric
uncertainty can also be regarded as a measure of explainability in the data
given the underlying (possibly incomplete) feature set.

Since our main interest lies in the development of novel features based on the
examination of data, i.e., data mining, we employ the aleatoric uncertainty
and use it to detect samples which contain characteristics not explicable by
the current feature set. The aleatoric uncertainty is modeled using the
tensorflow probability distribution library (dillon2017tensorflow, ) via
artificial neural networks by assuming that the priors over the predicted data
follow a normal distribution. Practically speaking, the aleatoric uncertainty
is modeled by the prediction of the mean and standard deviation
$\hat{\mu}_{i},\hat{\sigma}_{i}$ which define a normal distribution of the
$i$th output component instead of a single deterministic value. The neural
networks predicted normal distribution is given as

$f(\underline{x})=\frac{1}{\sqrt{2\pi\hat{\sigma}_{i}^{2}}}\exp{\Big{(}\frac{\hat{\mu}_{i}-y_{i}}{2\hat{\sigma}_{i}^{2}}\Big{)}}\,=p(y_{i}|\underline{x})\,,$
(3)

with the model $f$ using the feature vector $\underline{x}$ to predict the
parameters $\hat{\mu}_{i}$ and $\hat{\sigma}_{i}$ for each component $y_{i}$
independently.

Thus, the neural network has to predict twice the number of output variables
when modeling the aleatoric uncertainty as described. The quantification of
the aleatoric uncertainty during training enters through the loss $\Phi_{i}$
for target value $y_{i}$, which is derived by minimizing the Kullback-Leibler
divergence depeweg2019modeling :

$\Phi_{i}\big{(}y_{i},{\cal
N}(\hat{\mu}_{i},\hat{\sigma}^{2}_{i})\big{)}=-\log\big{(}\frac{1}{2}p(y_{i}\lvert{\cal
N}(\hat{\mu}_{i},\hat{\sigma}^{2}_{i})+s)\big{)},$ (4)

where $p(y\lvert{\cal N}(\hat{\mu}_{i},\hat{\sigma}^{2}_{i}))$ is the
probability of observing $y_{i}$ given the current distribution ${\cal
N}(\hat{\mu}_{i},\hat{\sigma}^{2}_{i})$. Since the logarithm quickly tends to
$-\infty$ for very small values in $p\left(y\lvert{\cal
N}(\hat{\mu}_{i},\hat{\sigma}^{2}_{i})\right)$, the shift parameter $s>0$ has
been introduced to stabilize the training. This greatly improved convergence
behaviour of the model for $s=0.25$.

For a more detailed explanation of the theoretical background and derivation
of Bayesian neural networks the authors refer to the literature, e.g.,
(goan2020Bayesian, ; tipping2003bayesian, ), presenting a more generous
outline to bayesian modeling and kabir2018neural  who reviews different
modeling methods in uncertainty quantification. An example application using
bayesian modeling is presented in kendall2017uncertainties .

####  2.2.2 Use of Bayesian neural networks in feature engineering

As previously outlined, the aleatoric uncertainty can be used to indicate lack
of feature knowledge that hinders more accurate predictions. For instance, the
volume fraction alone is unable to consider particle shapes and their
orientation distribution (Fig. 2 b)), which is reflected in the aleatoric
uncertainty. Investigating these samples of high aleatoric uncertainty, it can
be seen that each of these samples contains one or several characteristics
inexplicable by the current feature set. Exploiting this information, our data
mining approach is guided by the Bayesian neural network (BNN). First, we
filter out a subset of our data by considering only samples of high aleatoric
uncertainty. With this subset at hand, we aim at finding suitable feature
descriptors to quantify the apparent characteristics.

Since a major motivation for machine learning is about gaining computational
efficiency, the selection of input features for the machine learned model
should bear computational efficiency in mind. This motivates us to develop
novel features which are either obtainable through computationally cheap
operations, or derived via the convolution operation.

The effect of the discrete convolution operation will be briefly motivated by
considering the 1D Sobel operator $\underline{k}^{\rm S}$

$\underline{k}^{\rm S}=\frac{1}{2}\begin{bmatrix}-1,&0,&1\end{bmatrix}$ (5)

which is applied to an arbitrary 1D signal
$\underline{s}\in{\mathbb{R}}^{n_{x}}$ via a discrete convolution: for any
admissible $i$, i.e., $1<i<n_{x}$

$\begin{array}[]{rl}\underline{f}=\underline{s}\ast\underline{k}^{\rm S}\to
f_{i}=&\sum\limits_{j=-1}^{1}s_{i+j}\,\,k^{\rm S}_{j+2}\\\\[11.38109pt]
=&\frac{1}{2}(s_{i+1}-s_{i-1})\,\,,\end{array}$ (6)

where the gradient information is obtained in the feature map $\underline{f}$
(c.f. central difference).

Since the kernel can be arbitrarily chosen, if one was to replace the Sobel
operator $\underline{k}^{S}$ with a normalized constant vector, one would
recover the result of the moving average in the feature map $\underline{f}$.

Figure 3: The convolution of one RVE with two different kernels is shown. The
superscripts ${\rm S,R}$ of $\underline{\underline{k}}$ represent the Sobel
filter and the searched rectangle, respectively. The feature map
$\underline{\underline{F}}$ to the right differs significantly when applying
different kernels to the same image.

Analogously, the discrete convolution of a 2d signal, i.e., of an image, is
conducted by adding a second dimension to the data and therefore a second
summation to the operation (equation (6)). A graphical overview of the
convolution operation is given a little later in section 2.3 in the
introduction of convolutional neural networks. Since the convolution is
evaluated at each admissible index position with a double summation, the
operation in itself becomes costly. The cost can be drastically reduced by
conducting the convolution in Fourier space when following the convolutional
theorem hunt1971matrix , i.e.

$\underline{\underline{F}}=\underline{\underline{A}}\ast\underline{\underline{k}}={\cal
F}^{-1}\big{(}\,{\cal F}(\underline{\underline{A}})\cdot{\cal
F}(\underline{\underline{k}})\big{)}\,.$ (7)

The discrete Fourier Transform induces periodicity, which is rather favorable
in our case where the microstructure, i.e., the RVE, and the fields are
periodic.

Figure 4: Iteration 1: A few samples which led to the highest relative
prediction error and highest aleatoric uncertainties are shown. These samples
were found by a model trained using only the reduced coefficients and were
used for the first iteration of feature development.

Similar to the motivation in 1d is the resulting feature map
$\underline{\underline{F}}$ depends on the kernel $\underline{\underline{k}}$,
even when applying it to the same image $\underline{\underline{A}}$. This
makes the convolution operation exceedingly flexible for the quantification of
characteristics. For instance, it could even be used to detect specific shapes
by interpreting the feature map as a match indicator of the kernel in the
image (Fig. 3). Thus, we develop feature descriptors by developing specific
kernels and let the BNN guide us in the process.

We suggest the following general strategy: the entire set of currently
available features is used to train a BNN. Thereafter, an additional data set,
i.e., a test set consisting of unseen data, is predicted and the subset of
samples of high prediction error and high aleatoric uncertainty (Fig. 5 inside
the red box) are investigated to identify common phenomena by the data
scientist/engineer. An efficient feature descriptor to quantify the patterns
is proposed. The new features are subsequently added to the existing feature
set, which are used to train a novel BNN. Therewith, the procedure was
repeated and samples of high aleatoric uncertainty were investigated,
postulating an iterative feature engineering approach. This is graphically
illustrated in Fig. 1 (left), and can be transferred to any application using
human interpretable data.

Figure 5: The relative prediction errors related to the aleatoric uncertainty
are shown for the models prediction of the test set in the first iteration.
Samples close to the top right corner (located inside the red box) are
displayed in Fig. 4 and used for feature engineering Figure 6: Iteration 2: A
few samples which led to the highest relative prediction error and highest
aleatoric uncertainties are shown. These samples were found by a model trained
with the reduced coefficients and the band feature coeffcients. These images
are studied in the second iteration of feature development.

The initial BNN was trained using the reduced coefficients of the two point
correlation function (2PCF), as described in (lissner2019data, ). In order to
isolate samples of high aleatoric uncertainty and high prediction error, we
consider a relative error measure

$\|e_{\rm
rel}\|_{2}=\frac{\|\underline{\bar{\kappa}}-\underline{\hat{\kappa}}\|_{2}}{\|\underline{\bar{\kappa}}\|_{2}}$
(8)

with the target value $\underline{\bar{\kappa}}$ and the prediction
$\underline{\hat{\kappa}}\leftarrow f(\underline{\underline{\rm RVE}})$, which
was taken as the mean $\hat{\mu}_{i}$ from the predicted normal distribution.
The error-uncertainty relationship in the first iteration is displayed in Fig.
5, where the entire test set was predicted. A few samples in the approximate
region of the top right corner of Fig. 5 are shown in Fig. 4.

After close investigation of these samples, it was found that many of these
samples of high aleatoric uncertainty share a diagonal structure, more
strikingly do they contain connected regions or even percolation in certain
directions. This led to the invention of the band features (Fig. 7), which aim
to detect linear connectivity. The band features have some resemblance with
the lineal path function (lu1992lineal, ). The lineal path function has been
previously adopted in multiple studies, e.g., being used as a feature
kalidindi2011microstructure , or as a descriptor to generate statistically
similar RVEs (scheunemann2015design, ; balzani2014construction, ). The lineal
path function computes the probability that a line segment of specific length
under a certain angle lies fully within the inclusion phase. It is quite
costly to compute for multiple angles and different lengths of the line
segment.

Figure 7: Each line on the left in a) denotes one scanned direction of the
band feature detectors. The determination of the band feature in vertical
direction is shown for two images in b) and c), where the full lines measures
the inclusion phase and the dashed line measures the matrix phase. The feature
detector $\underline{\underline{b}}_{d}$ in b) and c) are drawn in their
approximate position for feature determination.

The newly proposed band features are computed with a single convolution per
direction $d$ as

$p^{\rm
I}_{d}=\max(\underline{\underline{b}}_{d}\ast\underline{\underline{\rm{RVE}}})$
(9)

where $\underline{\underline{b}}_{d}$ denotes the detection band in the
indexed direction $d$. Note that each detector $\underline{\underline{b}}_{d}$
is normalized such that ${p^{\rm I}_{d}\in[0,1]}$. The maximum operation in
(9) virtually fixes the evaluation of the band feature on a single spot in the
RVE, i.e., the location where the band feature detector is within the
inclusion phase for the highest amount.

In addition to the linear connectivity found in the inclusion phase (Fig. 4),
the absence of inclusions in the respective direction is captured by

$\begin{array}[]{rl}\widehat{p}^{\rm
M}_{d}&=1-\min(\underline{\underline{b}}_{d}\ast\underline{\underline{\rm{RVE}}})\\\
&=\max\big{(}\underline{\underline{b}}_{d}\ast(1-\underline{\underline{\rm
RVE}})\big{)}\,,\end{array}$ (10)

where $1-\underline{\underline{\rm RVE}}$ is a phase inversion of the image in
the bi-phasic setting. Generally, ${(p^{\rm I}_{d}-1)\neq p^{\rm M}_{d}}$, due
to the minimum/maximum operation.

The band feature detectors introduce hyperparameters, one being the directions
to be considered with the band feature detectors (see Fig. 7 a) ), and one
being the width of the band feature detector, which controls the minimum
connection width. In this study, we chose the width of the band feature
detectors to be 4px, and the angle increments of the band feature detectors
are set to be $\frac{\pi}{8}$ for both phases.

Remark: An extension of the approach to multiphase materials can be gained by
using the band features on the individual phase indicator functions.

The novel band feature coefficients have been used to enrich the existing
feature set, and a new model was trained. After convergence, the model is once
again used to predict the test set and samples of high aleatoric uncertainty
are displayed in Fig. 6. It can be seen that the samples leading to high
aleatoric uncertainty now often contain dispersed inclusions which are spread
over the whole RVE. It can also be seen that the found inclusions are
comparatively small with some overlap. If multiple particles are connected
this often leads to non-convex inclusion clusters. Consequently, we seek
features that are able to characterize inclusion dispersal and inclusion
shape.

Figure 8: Two artificial microstructures are shown with their respective
effective heat conductivity to the right. On the top/left the respective
reduced line projections are graphically given. The RVEs only vary in the
horizontal position of half the inclusions but differ noticeably in their
effective response. Figure 9: Iteration 3: A selection of samples which led to
the high relative prediction error and high aleatoric uncertainties are shown.
These samples were found by a model trained using all the previously mentioned
features.

Physically motivated the next feature descriptor aims to quantify global flux
hindrance, i.e., if there are multiple inclusions which form a disconnected
barrier in a certain direction. This measure is computed by reducing the
2d-image to a 1d-line and noting if there is at least one voxel containing the
inclusion phase. An example is illustrated in Fig. 8, where the full line (on
the top left) indicates the disconnected barrier. Taking the average of the
line, we obtain a scalar valued feature for each direction. This is
implemented in x- and y-direction via (given in python pseudocode)

$w_{i}=\text{mean\Big{(} sum( $\underline{\underline{\rm RVE}}$,
axis$=$}i)\geq 1\Big{)}$ (11)

where $\geq 1$ checks if there is at least 1 pixel found and the sum operator
is conducted in axis/dimension $i$, following the syntax of numpy.sum.

To quantify inclusion dispersal and approximate size of inclusions the RVE was
subdivided into a grid of local regions/cells. The local volume fraction
$c_{IJ}$ ($I,J$ in the coarse grid) was computed, coinciding with average
pooling (gholamalinezhad2020pooling, ) (cf. section 2.3). The local volume
fraction on the coarse grid yields a spatial distribution of the relative
amount of inclusion phase, see Fig. 10 left.

Figure 10: The process of extracting the distribution of volume fraction/edges
is shown above, reducing the $400\times 400$ image data to a $8\times 8$
distribution with the chosen parameters. On the left the volume fraction
distribution is shown, and on the right the edge distribution for $k^{2}$ is
shown.

Taking the mean of $\underline{\underline{c}}$ the global volume fraction is
obtained, i.e., average pooling is volume preserving. New insights can be
gained by taking the standard deviation
$\sigma^{f}(\underline{\underline{c}})$ and skewness
$\widetilde{\mu}^{f}_{3}(\underline{\underline{c}})$ of the distribution,
which approximately represents the size of inclusion clusters. Additionally,
the number of cells containing a certain volume fraction were counted.
Firstly, cells containing only one phase were counted, i.e., cells with the
local volume fraction of 0 or 1. The remaining cells were further subdivided
into thirds, such that the cells were counted as

$\underline{c}=\frac{1}{n_{\rm
cells}}\left\\{\begin{array}[]{lcrccr}\\#\big{(}&&f_{\rm
cell}&<&\epsilon&\big{)}\\\\[2.0pt] \\#\big{(}&\epsilon&\leq f_{\rm
cell}&<&\frac{1}{3}-\epsilon&\big{)}\\\\[2.0pt]
\\#\big{(}&\frac{1}{3}-\epsilon&\leq f_{\rm
cell}&<&\frac{2}{3}-\epsilon&\big{)}\\\\[2.0pt]
\\#\big{(}&\frac{2}{3}-\epsilon&\leq f_{\rm
cell}&<&1-\epsilon&\big{)}\\\\[2.0pt] \\#\big{(}&1-\epsilon&\leq f_{\rm
cell}&&&\big{)}\end{array}\right.,$ (12)

where $\\#(\bullet)$ denotes the operation counting the number of cells
fulfilling the condition, $f_{\rm cell}$ being the local cell volume fraction
($c_{IJ}$) and $\epsilon\leq\frac{1}{n_{x}\cdot n_{y}}$ a numerical tolerance.
Here one more hyperparameter is introduced, namely the window size of the
average pooling (we chose $\frac{1}{8}n_{x}\times\frac{1}{8}n_{y}$ in this
study), which determines the size of the coarse grain voxels. In addition to
that one could also consider a different partitioning for the counting
operation where we fixed it to be divided into equal thirds.

Table 1: The proposed features, are summarized in the table. The described phenomena of each feature as well as number of features and hyperparameters is shown. Note that when changing some hyperparameters the number of features may also change. feature name |  number of features |  number of hyperparameters |  captured phenomena |  related equation
---|---|---|---|---
band features | 16 | 2 | linear phase connectivity |  (9), (10)
global directional mean | 2 | 0 | global flux hindrance |  (11)
volume fraction distribution | 7 | 2 | inclusion size and dispersal |  (12)
directional edge distribution | 12 | 1 | inclusion shape and size |  (13),(14)

Another feature descriptor is proposed to further quantify the size and
sphericity of the inclusions through the surface information, since the ratio
of $\frac{\rm inclusion\,area}{\rm inclusion\,perimeter}$ assist in the
estimation. The ratio grows with inclusion size and is generally larger for
circles than rectangles. The surface information can be found via a
convolution, where edge detectors like the Sobel operator (vernon1991machine,
) are able to detect the (directional) surfaces of inclusions. The feature map
of edges $\underline{\underline{E}}^{e}$ can be computed as

$\underline{\underline{E}}^{e}=\lvert{\underline{\underline{k}}^{e}}\ast\underline{\underline{\rm
RVE}}\lvert,\qquad e\in\\{1,2,3,4\\}$ (13)

for different edge detectors $\underline{\underline{k}}^{e}$. The chosen edge
detectors are the horizontal, vertical and a diagonal Sobel filters, i.e.,

$\begin{array}[]{rl@{\qquad}rl}\underline{\underline{k}}^{1}&=\begin{bmatrix}-1&0&1\end{bmatrix}\,,&\underline{\underline{k}}^{2}&=\begin{bmatrix}0&-0.5&-1\\\
0.5&0&-0.5\\\ 1&0.5&0\end{bmatrix}\,,\\\\[17.07164pt]
\underline{\underline{k}}^{3}&=\begin{bmatrix}-1&-0.5&0\\\ -0.5&0&0.5\\\
0&0.5&1\end{bmatrix}\,,&\underline{\underline{k}}^{4}&=\begin{bmatrix}-1\\\
0\\\ 1\end{bmatrix}\,.\end{array}$ (14)

The absolute value is taken in (13) since the orientation of the surface
normal is irrelevant for our investigations. Similar to the local volume
fraction the edge feature maps are processed with average pooling, leading to
a coarse grid representation ${\underline{\underline{S}}^{e}}$. The mean
$\mu^{E}(\underline{\underline{S}}^{e})$, standard deviation
$\sigma^{E}(\underline{\underline{S}}^{e})$ and skewness
$\widetilde{\mu}_{3}^{E}(\underline{\underline{S}}^{e})$ is taken, leading to
three additional features per edge detector. The same window size as for the
local volume fraction $\underline{\underline{c}}$ was used, which constitutes
one hyperparameter for this feature. Additionally, other edge detectors, e.g.,
the Laplacian, could be considered.

With the auxiliary features at hand a new BNN model was trained for the next
iteration, with all of the previously used and newly introduced features. Once
again the RVE samples of high aleatoric uncertainty and high prediction errors
ares shown in Fig. 9. As can be seen in most of the samples, clusters of
inclusions are formed with narrow connections, leading often to non-convex
edges. Similarly, inclusions do almost form a cluster with a few pixels
distance in between them. No efficient method to quantify these phenomena was
found, and other features which have been tested did not notably improve the
models prediction. Thus, the feature engineering was halted at the third
iteration of feature engineering.

In total $37$ new features have been proposed by utilizing four different
feature descriptors, which are compactly summarized in tab. 1. Each of the
features is either obtainable via convolution, or obtainable by
computationally efficient operations.

###  2.3 Convolutional neural networks

Convolutional neural networks (Conv Nets) (o2015introduction, ) are derived
from feed forward neural networks (basheer2000artificial, ; lissner2019data, )
to make the classification/regression of high-resolution images tractable via
weight sharing (lecun1989generalization, ). The evaluation of Conv Nets is
similar to dense feed forward neural networks and conducted in one forward
pass. The matrix multiplication used in dense feed forward neural networks is
replaced by a convolution

$\underline{\underline{A}}_{l+1}=f_{l}(\underline{\underline{W}}_{l}\ast\underline{\underline{A}}_{l}+b_{l})\,,$
(15)

with $\underline{\underline{A}}_{l+1}$ being the image output of layer $l$,
$\underline{\underline{W}}_{l}$ denoting the kernel, $b_{l}$ the element wise
added bias/offset, and $f_{l}$ the activation function. Here the kernel
$\underline{\underline{W}}_{l}$ and the bias $b_{l}$ are the trainable
parameters of the neural network.

Figure 11: A schematic overview of a vanilla convolutional neural networks is
shown. The supposed feature extraction is conducted in the convolutional
layers (displayed as the light blue rectangles), where each rectangle
symbolizes one channel. The last convolutional layer is followed by a
flattening operation and a subsequent dense feed forward neural network
(displayed with the green neurons).

One major advantage of Conv Nets over dense neural networks is that the small
kernel (e.g. $3\times 3$ is a popular choice) is used to scan over the entire
image, processing the entire information with only few parameters. To further
improve information processing, multiple channels per layer can be specified
(illustrated through multiple slices in Fig. 11), where each channel
represents a different feature map (of same resolution). This leads to kernel
matrices being of size ${m_{y}\times m_{x}\times n_{\rm in}\times n_{\rm
out}}$, having different kernel weights for each input-output channel
combination. Thus, each output channel of the current layer is obtained by
processing all input channels with independent kernels, which results are
generally averaged.

Figure 12: A convolution using a $3\times 3$ kernel (displayed in the middle)
with stride 2 applied to a $7\times 7$ image (left) is shown. The resulting
image of size $3\times 3$ after convolution is given to the right. Each
colored box denotes the discrete positions of the kernel for the convolution
operation, matched on input and output image.

To keep the computational overhead feasible, downscaling of the image
resolution is conducted as more channels are added, which is generally
achieved via stride or pooling. The stride is implemented through the
convolution, denoting the stepping width of kernel evaluation interval, i.e.,
if stride=2 then the kernel gets evaluated on every other position, and a
downscaling of approximately factor two is achieved along each dimension (Fig.
12).

Figure 13: Average and max pooling is exemplarily displayed for an arbitrary
$9\times 9$ image on the left, where the discrete values after the operation
are shown in the respective colors. The red rectangles highlight the spatial
field considered for each output value. On the right ${10\times 10}$ average
and max pooling was applied for an exemplary microstructure image. The pooling
’grid’ is shown on the image in red, where each grid corresponds to one output
pixel.

A pooling layer replaces the convolution operation with e.g. the maximum or
average operation (Fig. 13), achieving downscaling analogously through the
stride. In a general case the stride in pooling is set equal to the size of
the pooling kernel. Interpreting the different pooling operations one can see
that max pooling singles out the peak activation of the kernel in the local
neighbourhood, whereas average pooling reflects the average occurrence of the
feature in a local region. After a pooling layer, the number of channels is
kept constant meaning that each channel is pooled individually.

Another special operation in Conv Nets are ${1\times 1}$ convolutions, which
are used to control/reduce the number of channels while generally retaining
the spatial resolution (with a stride$=1$). Additionally, a nonlinear
activation function can be deployed after the ${1\times 1}$ convolution.

With these operations at hand the convolutional layers can be built, where the
schematic overview of a Conv Net in Fig. 11 has four convolutional layers,
where each layer except the first one has a stride$>$1 and an increasing
number of output channels per layer. After the last convolutional layer, the
derived features in the spatial image representation are flattened into a
1d-vector and a dense feed forward neural network is deployed to conduct
regression/classification. Hence, the convolutional channels are often
interpreted as feature extractors, whereas the dense neural network is
utilized for the prediction based on the features extracted by the
convolutional layers. Such neural networks will be referred to as generic Conv
Nets below.

#### Padding in Conv Nets

The convolutional kernels derive features based on neighbourhood information
and consequently, on the boundary of the image where the neighbourhood is
undefined, a loss of information is introduced. To enable the Conv Net to
correctly consider the full neighbourhood relationship, periodic padding
before each convolution operation was deployed. This has been implemented in
tensorflow similar to (schubert2019circular, ) and is made publicly available
in (codebase, ), where also the implementation for the data augmentation
scheme is found.

Figure 14: The original RVE (inside red box) is padded by half the kernel
size. The padding is highlighted in orange. The blue rectangles show the
considered region while evaluating the convolution with a kernel at the dot
position for two exemplary convolution locations.

For illustration a periodically padded image is shown in Fig. 14, where it can
be seen that the padding ensures the spatial resolution to stay constant after
the convolution operation, since the kernel evaluation is only defined within
the red box of original spatial resolution. Note that in a deep Conv Net
periodic padding is deployed in each layer.

#### Data augmentation

Data augmentation aims at increasing the size of the training set, which
inherently has a regularizing effect kauderer2017quantifying . One study has
successfully applied cutout, which goes as far as to simply gray out larger
regions of the input image during training (devries2017improved, ). We aim to
virtually increase the size of the training set by translating the frame of
the RVE, which does not alter the macroscopic material due to assertion of
periodicity (Fig. 15). This additionally assists the Conv Net in learning the
translational invariance of the RVE frames, and can generate up to $400^{2}$
(i.e. image resolution) snapshots/samples of the input image per RVE. For
memory considerations, this is implemented ’online’ during training, where 50%
of the RVE frames are randomly translated every 10-th epoch.

Figure 15: The blue box highlights the originally generated RVE which
characterizes the macroscopic material by a repeated continuation in each
direction. The green box outlines a completely different image of the same RVE
obtained from translation in the periodic medium.

###  2.4 Deep inception module

In the context of microstructure homogenization, the material behaviour is
impacted by various factors, e.g., locations with narrow gaps between
inclusions, as well as large inclusion clusters. In the special case of
percolation, it is relevant that one inclusion cluster stretches over the
entire RVE, often in a curved manner, which also has to be detected in order
to make accurate predictions of the materials behaviour. Thus, different
characteristics with different relative size in the image have to be detected
by the convolutional layers.

In order to detect these various features, we deploy parallel branches of
differently sized kernels. Thereby, we design the Conv Net to capture
differently sized characteristics. Something similar, namely the inception
module has been previously implemented by szegedy2015going . Their original
intention was to reduce computational overhead and memory usage by increasing
width instead of depth in Conv Nets. Their resulting model, which only used
inception modules (Fig. 16), achieved state of the art results.

Figure 16: The original inception module is graphically illustrated where
multiple convolution operations are conducted between the layers. The graphic
layout is copied from (szegedy2015going, ) fig.2 and slightly modified.

The inception modules were implemented by replacing the convolution operation
between two layers by multiple parallel convolutions of differently sized
kernels, where each of the convolution operations has padding and a stride of
1, since the feature maps are concatenated channel wise at the end of the
module. To reduce the spatial resolution in the deep Conv Net,
szegedy2015going  used pooling between multiple inception modules.

Building upon the idea of multiple parallel convolutions, we increase the
depth in each parallel convolution branch, and propose the deep inception
modules (Fig. 17). The constraint of stride$=1$ in the convolutional layers is
dropped, and the number of operations in each branch can also be flexibly
adjusted. Then, we design each branch to capture different sized phenomena,
e.g., by a preceding average pooling/coarse graining***Note that large average
pooling introduces only a minor loss in information, since edge information is
mostly retained when using float values (c.f. Fig. 13), we can easily increase
the receptive field of the first $5\times 5$ convolution to be $25\times 25$
pixels (Fig. 17 rightmost branch). Additionally, we ensure that there exist
also branches which deploy convolution operations directly using the raw image
information, to capture small sized effects within the RVE (Fig. 17 leftmost
branch).

Figure 17: One deep inception module is shown where multiple deeper
convolution branches are conducted in parallel. Each kernel is quadratic and
the dimensions are read as size/stride $\cdot n_{\rm channels}$.

The only remaining constraint in the deep inception modules is that the
downsampling factor through stride/pooling has to match in each branch if the
channels are concatenated after the convolution operation. In Fig. 17, each
branch individually achieves a downsampling of a factor of $40$, and the
channels of each branch are concatenated at the end.

###  2.5 The hybrid neural network

Figure 18: The developed hybrid neural network utilizes parts of a generic
dense feed forward neural networks as well as parts of a Conv net with deep
inception modules. The prediction of the effective property is obtained by
summing up the prediction of each sub-part of the model.

So far two different model archetypes, i.e., the dense feed forward neural
network (FFNN) using the handcrafted features and the convolutional neural
network (Conv Net) are introduced to predict the homogenized material
property. Both model archetypes are able to outperform the reference model
lissner2019data  and we aim to obtain further improvements by combining the
FFNN and Conv Net in parallel. To additionally support the prediction, a
bypass using only the volume fraction, which is the single most impactful
variable in homogenization, has been implemented to serve as a baseline
prediction. The resulting hybrid neural network is shown in Fig. 18, which
consists of three different contributions to the models prediction. The models
prediction is given as

$\tilde{\underline{\kappa}}=\tilde{\underline{\kappa}}_{\rm
vol}+\Delta(\tilde{\underline{\kappa}}_{\rm
features}+\tilde{\underline{\kappa}}_{\rm ConvNet})\,,$ (16)

where each predicted subpart $\hat{\underline{\kappa}}_{\bullet}$ denotes one
subbranch of the hybrid neural network. The idea is that each branch considers
increasingly high level features in order to predict the effective heat
conductivity, which is governed by complex geometrical effects. To assist the
models convergence, a multistage training is implemented. The entire scheme is
graphically illustrated in Fig. 19.

Figure 19: The multiple stages of the hybrid models training are shown. If a
model subpart is transparent in a stage, then gradients were not computed
thereof. The line connections represent which part of the model contributed to
the prediction in the current stage. Each stage except for stage 2 is trained
to convergence.

Firstly, branch (I) using the volume fraction is trained independently until
convergence (stage 1). After convergence the parameters of (I) are frozen, but
the branch will contribute to all subsequent predictions. As a next step, the
entire model is trained for only a few epochs (stage 2), to move the model in
the proximity of a local minimum. Since the handcrafted features (section
2.2.2) remain static during training, branch (II) using these features is
trained until convergence (Stage 3), which trainable parameters are frozen
thereafter. In the next step the last branch (III) containing the Conv Net is
trained until convergence (stage 4), while the feature branch (II) contributes
to the prediction. The motivation is that the Conv Net can capture high level
features which explain the remaining variations uncaptured by the handcrafted
features. In the final step, the interdependent parameters are finetuned by
once again training the entire model until convergence (except for branch
(I)).

Figure 20: A few samples out of the benchmark set are displayed. Note that the
combination of different inclusions within one RVE is unseen to the machine
learning algorithm and represents truly unseen data. Some samples of the
training data were shown in Fig. 2.

##  3 Results

During the process of finding the best model, i.e., the hybrid neural network,
a multitude of neural networks has been trained in the step by step procedure.
The results of each intermediate step are compactly summarized in tab. 2 and
visualized via a R2 plot in the appendix 31. Further details on the different
prediction contributions, i.e., the handcrafted features and the convolutional
neural networks (Conv Nets) are given separately. In general, for the
different network archetypes, different losses have been deployed. Regarding
model parameter optimization, the same state of the art optimizers have been
used across the different models. This consists of ADAM gradient back
propagation with weight decay loshchilov2017decoupled , early stopping and a
constant learning rate with default hyperparameters. The model training has
been implemented using the tensorflow API tensorflow2015whitepaper , and the
code for training and post processing can be found in codebase .

All models have been trained using the same data outlined in section 2.1
(elaborated in appendix A). Validation was done using a benchmark dataset of
completely unseen inclusion shapes whilst a mix of inclusions within a single
RVE is allowed (Fig. 20). As an intuitively interpretable metric we use the
relative error measure:

$e=100\cdot\sqrt{\frac{1}{n_{\rm samples}\cdot
n_{\kappa}}\sum\limits_{i=1}^{n_{\rm
samples}}\sum\limits_{j=1}^{n_{\kappa}}\frac{(y_{ij}-\hat{y}_{ij})^{2}}{y_{ij}^{2}}}\,,$
(17)

which coincides with the error measure introduced in (8) and will be denoted
with rel. $\sqrt{\text{MSE}}$ in the future. Here $y_{ij}$ denotes the target
values, $\hat{y}_{ij}$ denotes the models prediction for all samples $i$ and
all components $j$ (of the heat conduction tensor).

###  3.1 Feature engineering and selection

During the feature engineering approach, in each iteration step a new Bayesian
neural network (BNN) was trained using the same architecture and the Bayesian
loss described in (4). The architecture is given in the appendix 5. The
features which have been added in each iteration step were motivated and
described in section 2.2.2. The actual improvements per step are given in tab.
2 and are graphically shown in Fig. 22. In Fig. 22, the density blur
highlights the locations where most of the model predictions coincide, and it
can be nicely seen that the center of it, i.e., the brightest spot, shifts
closer to the origin in each iteration. It can also be seen that the mean
prediction error improves whilst the average predicted uncertainty also
decreases. Note that in the first iteration the coefficients of the 2-point
correlation function (2PCF) has been used lissner2019data , in the second
iteration the band features were added, and in the third iteration the
remainder of the features was added, consisting of the volume fraction
distribution, the directional edge distribution and the global directional
mean (c.f tab. 1).

Figure 21: The absolute value of the Pearson correlation scores for all
features as well as the effective thermal properties of 3000 samples are given
above for the phase contrast of $R=5$ and $R=100$. Absolute values greater
than 0.85 are blacked out to indicate high correlation. The marked labels on
the ordinate note which feature is displayed up to the label (read from top to
bottom, ordered as introduced in section 2.2.2), and the abscissa up to which
index the current feature is given in the feature vector. The six features
beyond the blue line denote the target values, i.e., the effective heat
conductivity for the different phase contrasts. Figure 22: A density plot is
shown which relates the aleatoric uncertainty to the rel. $\sqrt{\text{MSE}}$
in each iteration step, overlayed with a scatterplot of the actual
predictions. Samples further away from the brightest spot in the density plot
are displayed larger for better visibility. The ’x’ marks the mean values.
Note that the abscissa of the left most plot is differently scaled.

As has been previously stated, any tests on additional features features
yielded little to no improvements. Since the improvements from the second to
the third iteration were smaller than in the first iteration, a feature
selection process was initiated to spot if the introduced features are
redundant or do not positively contribute to the prediction. The correlation
matrix in Fig. 21 shows the absolute values of the Pearson correlation scores,
indicates that the features are suitably uncorrelated, though some
correlations remain within each feature class.

In a further investigation step, the features have been ranked by two filter
and wrapper methods. Filter methods estimate feature importance on plain data
observation, and wrapper methods rate the importance based on an intermediate
surrogate model which measures scores of prediction contribution. The deployed
filter methods used the Pearson correlation scores and the analysis of
variance (ANOVA) via F-scores (feir1974empirical, ). The deployed wrapper
methods were the recursive feature elimination (RFE)(guyon2002gene, ) and a
method deploying random forests (geurts2006extremely, ). Each of these methods
found a slightly different order, presented as indices, given in arrays:

$\begin{gathered}\text{\small{Pearson ranking:}}\\\
\begin{array}[]{lllllllllllllll}\large[&{\color[rgb]{0,0.25390625,0.5703125}{0}},&{\color[rgb]{0,0.25390625,0.5703125}{38}},&{\color[rgb]{0,0.74609375,1}{22}},&{\color[rgb]{0,0.74609375,1}{14}},&{\color[rgb]{0,0.74609375,1}{32}},&{\color[rgb]{0,0.25390625,0.5703125}{24}},&{\color[rgb]{0,0.25390625,0.5703125}{28}},&{\color[rgb]{0,0.74609375,1}{16}},&{\color[rgb]{0,0.25390625,0.5703125}{20}},&{\color[rgb]{0,0.25390625,0.5703125}{26}},&{\color[rgb]{0,0.74609375,1}{18}},&{\color[rgb]{0,0.74609375,1}{15}},&{\color[rgb]{0,0.74609375,1}{23}},&\\\
&17,&21,&25,&1,&29,&19,&27,&36,&35,&31,&30,&4,&33,&\\\
&49,&40,&46,&5,&43,&48,&39,&45,&42,&8,&50,&9,&\\\
&{\color[rgb]{1,0.36328125,0}{37}},&{\color[rgb]{1,0.36328125,0}{3}},&{\color[rgb]{0.578125,1,0}{11}},&{\color[rgb]{1,0.36328125,0}{2}},&{\color[rgb]{0.578125,1,0}{44}},&{\color[rgb]{0,0.25390625,0.5703125}{47}},&{\color[rgb]{0.578125,1,0}{10}},&{\color[rgb]{0,0.74609375,1}{12}},&{\color[rgb]{0,0.25390625,0.5703125}{34}},&{\color[rgb]{0,0.74609375,1}{41}},&{\color[rgb]{0,0.25390625,0.5703125}{7}},&{\color[rgb]{0.578125,1,0}{6}},&{\color[rgb]{0,0.25390625,0.5703125}{13}}\large]\end{array}\\\\[5.69054pt]
\text{\small{ANOVA ranking:}}\\\
\begin{array}[]{lllllllllllllll}\large[&{\color[rgb]{0,0.25390625,0.5703125}{0}},&{\color[rgb]{0,0.25390625,0.5703125}{38}},&{\color[rgb]{0,0.25390625,0.5703125}{22}},&{\color[rgb]{0,0.74609375,1}{15}},&{\color[rgb]{0,0.74609375,1}{14}},&{\color[rgb]{0,0.74609375,1}{32}},&{\color[rgb]{0,0.25390625,0.5703125}{24}},&{\color[rgb]{0,0.74609375,1}{23}},&{\color[rgb]{0,0.25390625,0.5703125}{28}},&{\color[rgb]{0,0.74609375,1}{16}},&{\color[rgb]{0,0.25390625,0.5703125}{26}},&{\color[rgb]{0,0.25390625,0.5703125}{20}},&{\color[rgb]{0,0.74609375,1}{18}},&\\\
&29,&21,&17,&25,&1,&19,&27,&36,&33,&35,&4,&30,&37,&\\\
&31,&2,&48,&5,&39,&42,&45,&46,&8,&40,&3,&11,&\\\
&{\color[rgb]{0,0.74609375,1}{49}},&{\color[rgb]{0.578125,1,0}{43}},&{\color[rgb]{0,0.25390625,0.5703125}{7}},&{\color[rgb]{0,0.74609375,1}{41}},&{\color[rgb]{1,0.36328125,0}{9}},&{\color[rgb]{0,0.25390625,0.5703125}{13}},&{\color[rgb]{0.578125,1,0}{10}},&{\color[rgb]{0,0.25390625,0.5703125}{34}},&{\color[rgb]{0.578125,1,0}{6}},&{\color[rgb]{0,0.74609375,1}{12}},&{\color[rgb]{1,0.36328125,0}{50}},&{\color[rgb]{0.578125,1,0}{44}},&{\color[rgb]{0,0.25390625,0.5703125}{47}}\large]\end{array}\\\\[5.69054pt]
\text{\small{forest ranking:}}\\\
\begin{array}[]{lllllllllllllll}\large[&{\color[rgb]{0,0.25390625,0.5703125}{0}},&{\color[rgb]{0,0.25390625,0.5703125}{38}},&{\color[rgb]{0,0.74609375,1}{23}},&{\color[rgb]{0,0.74609375,1}{14}},&{\color[rgb]{0,0.74609375,1}{22}},&{\color[rgb]{0,0.74609375,1}{32}},&{\color[rgb]{0,0.74609375,1}{15}},&{\color[rgb]{0,0.25390625,0.5703125}{20}},&{\color[rgb]{0,0.25390625,0.5703125}{24}},&{\color[rgb]{1,0.36328125,0}{25}},&{\color[rgb]{0.578125,1,0}{3}},&{\color[rgb]{0,0.25390625,0.5703125}{26}},&{\color[rgb]{0,0.25390625,0.5703125}{28}},&\\\
&18,&17,&16,&21,&29,&27,&2,&19,&33,&31,&30,&37,&4,&\\\
&1,&5,&50,&44,&8,&9,&35,&42,&48,&10,&6,&43,&\\\
&{\color[rgb]{0.578125,1,0}{11}},&{\color[rgb]{1,0.36328125,0}{39}},&{\color[rgb]{0,0.74609375,1}{41}},&{\color[rgb]{0,0.25390625,0.5703125}{7}},&{\color[rgb]{0,0.25390625,0.5703125}{13}},&{\color[rgb]{0.578125,1,0}{45}},&{\color[rgb]{0,0.25390625,0.5703125}{47}},&{\color[rgb]{0,0.74609375,1}{12}},&{\color[rgb]{0.578125,1,0}{36}},&{\color[rgb]{0.578125,1,0}{46}},&{\color[rgb]{1,0.36328125,0}{40}},&{\color[rgb]{0,0.74609375,1}{49}},&{\color[rgb]{0,0.25390625,0.5703125}{34}}\large]\end{array}\\\\[5.69054pt]
\text{\small{RFE ranking:}}\\\
\begin{array}[]{lllllllllllllll}\large[&{\color[rgb]{0,0.25390625,0.5703125}{0}},&{\color[rgb]{0,0.74609375,1}{18}},&{\color[rgb]{0.578125,1,0}{3}},&{\color[rgb]{0,0.25390625,0.5703125}{26}},&{\color[rgb]{0,0.25390625,0.5703125}{38}},&{\color[rgb]{0,0.25390625,0.5703125}{20}},&{\color[rgb]{0,0.25390625,0.5703125}{24}},&{\color[rgb]{0,0.25390625,0.5703125}{28}},&{\color[rgb]{1,0.36328125,0}{4}},&{\color[rgb]{0,0.74609375,1}{16}},&{\color[rgb]{1,0.36328125,0}{2}},&{\color[rgb]{1,0.36328125,0}{27}},&{\color[rgb]{1,0.36328125,0}{21}},&\\\
&14,&22,&19,&15,&23,&17,&39,&25,&1,&8,&31,&29,&50,&\\\
&5,&41,&30,&9,&40,&44,&6,&11,&10,&35,&12,&37,&\\\
&{\color[rgb]{0,0.25390625,0.5703125}{47}},&{\color[rgb]{0.578125,1,0}{43}},&{\color[rgb]{0.578125,1,0}{45}},&{\color[rgb]{1,0.36328125,0}{48}},&{\color[rgb]{0,0.25390625,0.5703125}{7}},&{\color[rgb]{0.578125,1,0}{46}},&{\color[rgb]{1,0.36328125,0}{32}},&{\color[rgb]{0,0.74609375,1}{49}},&{\color[rgb]{0,0.25390625,0.5703125}{13}},&{\color[rgb]{1,0.36328125,0}{33}},&{\color[rgb]{1,0.36328125,0}{42}},&{\color[rgb]{0,0.25390625,0.5703125}{34}},&{\color[rgb]{0.578125,1,0}{36}}\large]\end{array}\end{gathered}$
(18)

where Python notation is used and the indices start at 0. Note that the first
and the last 13 features are color coded to display agreement between the
feature selection methods. Features marked by dark blue have been ranked as
the top/bottom 13 features by all methods, light blue by three, green by two,
and orange by one. The utmost curious reader may compare the indices
represented above with the indices and corresponding feature displayed in the
correlation matrix (Fig. 21). Overall, the volume fraction and most of the
band features were rated as the most significant for all methods, and
generally some edge distribution features and reduced coefficients of the 2PCF
were commonly rated with the lowest scores.

Figure 23: The validation loss development with respect to the total number of features is shown. Each scatter shows the lowest validation loss with the best BNN out of five for the different feature selection methods. Table 2: Error measures are given for each step in the improvement process of the artificial neural networks. The full numbers present the error measures computed on the benchmark set and the values given in the brackets refer to the error measures computed on a test set affine to the training data. The different model types are separated via horizontal lines. Model type |  mean relative error [%] |  median relative error [%] |  rel. $\sqrt{\text{MSE}}$ [%] |  MSE $10^{-5}$ [-] |  trainable parameters [-]
---|---|---|---|---|---
| $\bar{\kappa}_{11}$ | $\bar{\kappa}_{22}$ | $\bar{\kappa}_{11}$ | $\bar{\kappa}_{22}$ |  |  |
volume fraction only |  6.99 (5.60) |  6.86 (5.81) |  5.89 (4.14) |  5.54 (4.24) |  8.91 (7.47) |  144.12 (91.2) | 241
reference model lissner2019data  |  1.92 (1.62) |  1.91 (1.63) |  1.37 (1.28) |  1.30 (1.24) |  2.39 (1.99) |  10.40 (6.47) | 241
feature iteration (1) |  1.40 (1.15) |  1.13 (1.19) |  1.04 (0.91) |  0.83 (0.96) |  1.90 (1.72) |  6.58 (4.84) | 4 997
feature iteration (2) |  1.26 (1.01) |  1.20 (1.02) |  1.00 (0.81) |  0.91 (0.79) |  1.71 (1.38) |  5.33 (3.11) | 4 997
generic Conv Net |  1.87 (1.91) |  1.81 (1.85) |  1.46 (1.55) |  1.41 (1.49) |  2.37 (2.29) |  10.21 (8.55) | 220 331
generic Conv Net \+ vol bypass |  1.70 (1.56) |  1.56 (1.46) |  1.21 (1.28) |  1.14 (1.21) |  4.91 (1.98) |  43.79 (6.43) | 220 359
inception Conv Net \+ vol bypass |  2.01 (1.90) |  2.11 (1.84) |  1.49 (1.51) |  1.65 (1.47) |  2.68 (2.43) |  13.07 (9.65) | 278 831
generic Conv Net \+ vol bypass +data augmentation |  $10^{8}$ (0.88) |  $10^{7}$ (0.90) |  0.65 (0.72) |  0.70 (0.71) |  $10^{8}$ (1.23) |  $10^{17}$ (2.45) | 220 359
inception Conv Net \+ vol bypass +data augmentation |  0.74 (0.77) |  0.79 (0.80) |  0.58 (0.62) |  0.62 (0.63) |  1.05 (1.04) |  2.00 (1.77) | 278 831
hybrid model |  0.73 (0.68) |  0.74 (0.71) |  0.55 (0.54) |  0.56 (0.56) |  1.03 (0.93) |  1.91 (1.42) | 283 703

In a final step, multiple BNNs (with the same architecture tab. 5) have been
trained using only a subset of the available features. The BNNs used
[6,9,12,…,51] features in each step and were trained in a ’best out of 5’
setup, meaning that five BNN have been trained for each feature subset ordered
by each feature selection method. The achieved validation losses are given in
Fig. 23. Generally, there is very little deviation seen between the different
feature selection methods, and more strikingly, the loss decreases almost
monotonously with respect to the number of features. Thus, all 51 features
have been chosen for the subsequent hybrid neural network since all of the
proposed features positively contribute to the prediction.

###  3.2 Convolutional and hybrid neural networks

Figure 24: The plots display the target values and the respective hybrid
neural network prediction for some randomly selected samples and the samples
with the highest prediction error for each component of the heat conduction
tensor. Each plot shares the same legend displayed in the middle. The bar
plots at the bottom are computed using the absolute error and show the error
quantiles on the current volume fraction interval. Note that the error bars
consider the entire dataset, even the samples which were not plotted in favor
of a less cluttered and better readable plot.

For the convolutional neural networks (Conv Net) the aleatoric uncertainty was
dropped and only one deterministic value was predicted. The deployed loss for
model optimization and calibration was the mean squared error (MSE). At first,
a generic Conv Net was deployed to predict the effective heat conductivity
using the available image data. The full architecture of the Conv Net is given
in appendix 6. This model underperformed previous state of the art results at
first, however, it has been continuously improved to outperform even the
feature based approach. Since the volume fraction is the most relevant
parameter in homogenization, the plain Conv Net was improved by implementing
the bypass using the volume fraction (equivalent to the hybrid model in
section 2.5), such that the Conv Net was forced to dedicate its attention in
detecting high level features which quantify phenomena explaining variations
around the volume fraction. The next improvement was found by implementing the
deep inception modules. The deep inception modules were both applied to the
input image, and their output was flattened and concatenated into one dense
regressor. The full architecture of each intermediate neueral network is given
in appendix 30.

The attentive reader might have noted that the loss of the model deploying
deep inception modules on the benchmark set in tab. 2 is slightly larger than
for the generic Conv Net, however the model using deep inception modules
achieved a low training loss of $1.7\cdot 10^{-5}$, outperforming all previous
models by a factor of 4. Thus, to fully capitalize on the deep inception
modules, the dataset was enriched through the proposed data augmentation
scheme, by randomly translating 50% of the input images every 10-th epoch
during training. The data augmentation scheme has also been deployed for the
generic Conv Net, however leading to an accurate but strongly overfitting
model, giving prohibitive prediction errors for few samples. This can be best
seen in the discrepancy between the mean and median error on the test and
benchmark set in tab. 2.

Further improvements were only found by the implementation of the hybrid
neural network, adding only very few additional parameters to the inception
Conv Net. As discussed in section 2.5, the hybrid model did not always
converge to a good local minimum, however, after implementing the multistage
training, it reliably found an excellent minimum (Fig. 25). Note that for
almost every model the validation loss is lower than the training loss in Fig.
25, which is explicable by the data augmentation scheme.

Figure 25: The achieved validation (full lines) and training loss (dashed
lines) is shown after model convergence. The hybrid neural network is compared
with and without pretraining for five randomly initialized trainings, as well
as the Conv Net using deep inception modules.

A more detailed error analysis of the hybrid model is given in figure format
in Fig. 24. The plot shows only a few sample predictions but it is ensures
that the samples with the highest prediction errors are shown. As can be seen,
the model is incredibly accurate, and is even able to accurately predict
outliers of the offdiagonal component $\bar{\kappa}_{12}$.

#### Physical correctness

Figure 26: The prediction of the hybrid neural network is given for 50 random samples in 100 different configurations with respect to frame translation. Predictions and samples are color matched to make them distinguishable for similar volume fractions. Table 3: The error measures have been computed for the hybrid neural network once in the original configuration of the RVE and once in the rotated configuration for the entire benchmark dataset. In the rotated configuration the components $\bar{\kappa}_{11}$ and $\bar{\kappa}_{22}$ are flipped and the off-diagonal component $\bar{\kappa}_{12}$ changes its sign.  RVE image configuration |  mean relative error [%] |  median relative error [%] |  relative MSE [-] |  MSE $10^{-5}$ [-]
---|---|---|---|---
|  $\bar{\kappa}_{11}$ | $\bar{\kappa}_{22}$ |  $\bar{\kappa}_{11}$ | $\bar{\kappa}_{22}$ |  |
|  |  |  |  |  |
original | 0.730 | 0.742 | 0.552 | 0.562 | 1.027 | 1.915
rotated | 0.759 | 0.751 | 0.564 | 0.574 | 1.045 | 1.982

In order to further quantify the accuracy of the hybrid neural network, a
brief study regarding physical correctness is conducted. Starting out with the
first property of frame translation invariance, which is not enforced but the
data augmentation scheme is of assistance. Note that the previously mentioned
manually derived features are designed to be frame translational invariant,
with the minor exception of the volume fraction/directional edge distribution,
which fluctuates ever so slightly based on the location of the pooling grid.

Some errors of the hybrid neural network with respect to frame translation are
seen in Fig. 26, where it can be seen that there is no configuration of the
RVE in which the model is significantly worse than in any other. The model
does not yield the same prediction for every possible configuration, however,
there is only slight variations between the predictions of the different
configurations.

Figure 27: The prediction of two RVEs out of the benchmark dataset is given
for the original representation of the RVE as well as the rotated
configuration. The models prediction $\hat{\underline{\underline{\kappa}}}$ as
well as the actual target values $\bar{\underline{\underline{\kappa}}}$ are
shown below each RVE.

Additionally, a different geometric transform can be used to generate
auxiliary samples, i.e., by rotating the RVE frames (Fig. 27). After rotation,
the distinct values of the heat conduction tensor $\bar{\underline{\kappa}}$
are swapped and do not change, except for the sign in $\bar{\kappa}_{12}$.
Ideally, this property should also be learned by the hybrid model, which was,
similar to the translational invariance, not the case up to minor fluctuations
(tab. 3). There seems to be an ever so slight bias in favor of the prediction
of $\bar{\kappa}_{11}$, which implies a minor bias in the training data. The
model seems to perform slightly better for the original configuration of the
RVE, however, this effect is negligibly small.

Note that the rotation of the RVE frames could also be implemented for data
augmentation, however the features would have to be recomputed on every
augmentation step since they are, correctly so, not rotational invariant.

The hybrid model with data augmentation is able to correctly reflect physical
behaviour up to an acceptable prediction error.

Table 4: Averaged error measure over the entire benchmark dataset are shown for the model being able to predict variable phase contrast. Depending on the target value of $\bar{\kappa}_{\bullet}$, a relative or absolute error measure is computed, values denoted with - are not defined (c.f.Fig. 32). The error measures are given at the discrete sampled phase contrasts the model was trained on. Phase contrast |  mean relative error [%] $\bar{\kappa}_{\bullet}>0.2$ |  mean absolute error [-] $\bar{\kappa}_{\bullet}\leq 0.2$ |  relative $\sqrt{\text{MSE}}$ $[\%]$ |  MSE $\cdot 10^{-5}[-]$
---|---|---|---|---
|  $\bar{\kappa}_{11}$ |  $\bar{\kappa}_{22}$ | $\bar{\kappa}_{11}$ | $\bar{\kappa}_{22}$ | $\bar{\kappa}_{12}$ |  |
|  |  |  |  |  |  |
2 | 0.61 | 0.59 | - | - | 0.0021 | 0.70 | 2.37
5 | 1.01 | 0.97 | - | - | 0.0043 | 1.33 | 3.60
10 | 2.07 | 2.01 | 0.0048 | 0.0053 | 0.0065 | 2.74 | 9.56
20 | 2.78 | 2.69 | 0.0087 | 0.0111 | 0.0090 | 4.71 | 18.98
50 | 3.17 | 3.36 | 0.0139 | 0.0159 | 0.0121 | 8.40 | 34.37
100 | 3.46 | 3.74 | 0.0159 | 0.0190 | 0.0134 | 12.16 | 44.50
Figure 28: The prediction for variable phase contrast is given for the full
interpolation range of R=[2,100] using $\Delta R=1$. To the left the averaged
prediction over the entire set is give. The right three plots, which share the
same legend, show the average prediction in each component, as well as the
prediction of two distinct samples of the highest and lowest relative error.
The predictions are always given in full lines and the reference solution in
dashed lines. The phase contrasts, which the hybrid neural network has seen
during training, are highlighted by dots.

####  3.2.1 Variable phase contrast

To further test the capabilities of the potent hybrid model, it was trained to
predict the thermal behaviour for variable phase contrast ranging from $R=2$
to $R=100$, while considering insulating inclusions (previously, R=5 was
fixed). During training, the input data of the model did not change compared
to the previously discussed data, however, the output data of the model
changed, where the effective heat conductivity of all training samples got
computed at the discrete values of the phase contrast
$R\in\\{2,5,10,20,50,100\\}$. In order to inform the hybrid model of the
variable phase contrast, one extra input neuron was added in each branch of
the model (Fig. 18 the red neuron), leading to less than $500$ additional
parameters. The phase contrast input parameter was linearly scaled from 0 to
1. One further adjustment was made, i.e., instead of the MSE the rel.
$\sqrt{\text{MSE}}$ was used as a cost function. The resulting prediction
errors on the benchmark set are compactly summarized in tab. 4, and presented
more elaborately in the appendix in Fig. 32. There it can also be seen that
the prediction becomes increasingly challenging for higher phase contrasts,
which the model was still able to predict with a moderately low rel.
$\sqrt{\text{MSE}}$ of 12% for the highest phase contrast. The models accuracy
for $R=5$ did slightly deteriorate compared to the previously shown hybrid
model, however, in the tradeoff for generalization capabilities, where the
model is able to predict every phase contrast with low prediction errors.

One major advantage of training the different phase contrasts in a single
model is the interpolation capabilities, where the prediction errors are given
in Fig. 28 for $R\in{1,2,3,...,100}$. The plot shows that the smooth
interpolation is in general accurately given for almost all phase contrasts
within the training range with a minor exception for the short interval of
$R\approx[60,72]$, where the prediction errors are slightly larger.

##  4 Summary

The prediction of the homogenized response is improved by novel features which
are developed with Bayesian assisted data mining. The aleatoric uncertainty is
used to collect samples which contain characteristics inexplicable to the
current feature set, commonalities within these samples were found and
multiple feature descriptors were developed, efficiently quantifying these
characteristics. In addition to the manually engineered features convolutional
neural networks are deployed, where we propose deep inception modules which
are a priori designed to capture phenomena at different length scales within
the microstructural image data. The major advantage of the deep inception
modules appear to be the improved generalization capabilities and its
regularizing effect. To utilize the deep inception modules to their full
potential, a data augmentation scheme is presented, which can generate more
than 100 000 input samples per data point, without increasing memory
consumption.

The two different neural network archetypes utilizing the newly engineered
features and deep inception modules are combined into a hybrid neural network,
deploying the archetypes in parallel to yield its prediction. To improve
convergence behaviour, a multistage training is implemented. The resulting
model was able to more than half the prediction error compared to previous
state of the models.

After extending the model to predict variable phase contrast, the proposed
hybrid neural network was able to accurately predict the material response for
the challenging data. It performs just slightly worse on single phase contrast
compared to the model trained only on the respective data, even without
touching the network layout. Remarkably, the variable contrast model is still
majorly outperforming the reference model lissner2019data  by almost halving
the prediction error, while enabling for completely arbitrary inputs. The
model does also perform reasonably well on interpolation, yielding accurate
predictions even for phase contrasts the model has not been trained on.

##  5 Discussion

One of the key motivations of machine learning is its efficient evaluation
during inference. The overhead of computing all the features as well as the
evaluation of the convolutional neural network is noticeable. Since all of the
features, including the reduced coefficients, are computable in Fourier space
or trivially obtainable, the FFT has to be conducted only once for each
sample. Additionally, the convolutional kernels can be pre-computed and stored
in Fourier representation when the resolution is fixed. One computational
downside for the band features is that in each direction the IFFT has to be
taken, since the maximum and the minimum in real space are required. Similarly
to the edge distribution where the absolute value is required with its pixel
location. Consequently, these features are only obtainable with an additional
computational overhead. When timing the feature computation as well as the
prediction, the model outperformed the intrinsically fast Fourier accelerated
solvers by a factor of $\approx 10$, where the feature derivation made up
about 80% of the computational effort by taking $\approx 77$ seconds when
considering 1500 samples, whilst the hybrid models prediction took $\approx 9$
seconds. Note that the prediction can be evaluated for multiple phase
contrasts once the input feature vector is obtained.

Further computational speedups can be obtained with additional hyperparameter
tuning of the deep inception modules. However, the obtained improvement is
already significant, such that additional hyperparameter tuning has been
omitted in favor of more future research, especially when considering the
required time investment to explore the actual infinite range of possible
variations (inception module depth, width, tuning of each branch, etc.), and
the possibility of nested deep inception modules.

Regarding the prediction accuracy, the hybrid neural network which uses the
manually engineered features as well as the deep inception modules was able to
significantly outperform different models. Ultimately, it has failed to learn
exact physical behaviour, which is not too surprising since the prediction is
solely based of a neural network without any constraints.

Acknowledgments Funded by Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under Germany’s Excellence Strategy - EXC 2075 –
390740016. Contributions by Felix Fritzen are funded by Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) within the Heisenberg
program DFG-FR2702/8 - 406068690 and DFG-FR2702/10 - 517847245. We acknowledge
the support by the Stuttgart Center for Simulation Science (SimTech).

## Declarations

  * •

Funding
Funded by Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)
under Germany’s Excellence Strategy - EXC 2075 – 390740016. Contributions by
Felix Fritzen are funded by Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) within the Heisenberg program DFG-FR2702/8 - 406068690
and DFG-FR2702/10 - 517847245. We acknowledge the support by the Stuttgart
Center for Simulation Science (SimTech).

  * •

Conflict of interest/Competing interests
The authors declare no conflict of interest

  * •

Availability of data and materials

The data is publicly available in the data repository of the University of
Stuttgart (DaRUS) darus1151

  * •

Code availability

The code will be made publicly available in codebase  upon publication.

##  Appendix A Data overview

The binary image data $A$ is represented via voxels as a discrete field, such
that we have

$\displaystyle A$ $\displaystyle\in{\mathbb{A}}^{n\times n},$ $\displaystyle
A_{ij}$ $\displaystyle=A({\boldsymbol{x}}_{ij}),$
$\displaystyle{\boldsymbol{x}}_{ij}$ $\displaystyle\in\varOmega,$ (19)

where the value $A_{ij}$ denote pixel-wise constant values and the data range
${\mathbb{A}}\subset{\mathbb{N}}_{0}$ depends on the number of discrete
phases. In general, $A$ represents image data, more specifically the image
data representing the RVE. The analyzed microstructure images are given by
binary voxel data, i.e., ${\mathbb{A}}=\\{0,1\\}$, and the resolution is fixed
to $400\times 400$, i.e., $i,j\in[0,1,\dots,399]$ (using python notation). The
binary color space $\mathbb{A}=\\{0,1\\}$ is displayed as black in the matrix
phase of the material (value $0$), and the light gray foreground represents
the inclusions (value 1), as can be seen in Fig. 2. The material is exemplary
for metal-ceramic or polymer-glass composite materials.

This study focuses on the thermal behaviour of the material, i.e., on the
effective heat conductivity†††Equivalently, this can be interpreted as a
prediction of the permeability of the material which has the same underlying
mathematical structure.. The properties of the macroscopic material are
induced by the parameters of the phases in the RVE, where a phase contrast of
$R=5$ is considered such that we have

$\displaystyle\kappa_{\rm M}$ $\displaystyle=1\,\frac{\rm J}{{\rm sKm}},$
$\displaystyle\kappa_{\rm I}$ $\displaystyle=\frac{\kappa_{\rm M}}{R}\,,$ (20)

i.e., low-conducting inclusions are considered. Here the subscript M denotes
the matrix phase, and the subscript I denotes the inclusion phase. Simulations
of the heat flow conducted via FANS (leuschner2018fourier, ) yield the
symmetric effective heat conduction tensor
$\tilde{\underline{\underline{\kappa}}}$ as

$\tilde{\underline{\underline{\kappa}}}=\begin{bmatrix}\tilde{\kappa}_{11}&\tilde{\kappa}_{12}\\\
\tilde{\kappa}_{12}&\tilde{\kappa}_{22}\\\ \end{bmatrix}\,,$ (21)

in the 2D setting. Due to the symmetry the target values are represented in
Mandel notation

$\tilde{\underline{\kappa}}=\begin{bmatrix}\tilde{\kappa}_{11}\\\
\tilde{\kappa}_{22}\\\ \sqrt{2}\tilde{\kappa}_{12}\end{bmatrix}\,.$ (22)

To obtain de-dimensionalized output values the heat conduction tensor is
normalized by the material property of the matrix phase, i.e.,

$\bar{\underline{\kappa}}=\frac{\tilde{\underline{\kappa}}}{\kappa_{\rm
M}}\,,$ (23)

again, given in Mandel notation. The investigated data, generated by our in-
house algorithms, is made publicly available (darus1151, ). It hosts all of
the data used for the methodological developments.

Figure 29: The effective heat conductivity for every component of
$\bar{\underline{\kappa}}$ is related to the volume fraction for the phase
contrast $R=5$. The Voigt and Reuss bounds (hashin1962variational, ) are shown
with dashed red lines for reference. The inspected data covers the admissible
variation within the physical bounds pretty well.

Considering the training data for the neural networks, the input data of is
shifted and scaled such that each input feature of the neural network has zero
mean and unit standard deviation in order to ensure statistical equivalence
for all features. The transformation of samples $X_{i},i=1,\dots,n$ into
$\widetilde{X}_{i}$ of a feature $X$ is computed via

$\displaystyle\overline{X}$ $\displaystyle=\frac{1}{n}\sum_{j=1}^{n}X_{j},$
$\displaystyle\widetilde{X}_{i}$
$\displaystyle=\frac{\sqrt{n-1}(X_{i}-\overline{X})}{\sqrt{\displaystyle\sum_{j=1}^{n}(X_{j}-\overline{X})^{2}}}\,.$
(24)

This process is performed for individually for each feature, i.e., zero cross-
correlation of the inputs is asserted for simplicity. Regarding the target
values $\bar{\underline{\kappa}}$: the diagonal components in the heat
conduction tensor are well defined in $\tilde{\kappa}_{ii}\in(0,1)$ and the
off-diagonal component $\tilde{\kappa}_{ij}$ fluctuates around 0 with
relatively small values due to the positive definiteness of the tensor (c.f.
Fig. 29). The physically possible range of output values is consequently
machine learning friendly and, therefore, no additional scaling is performed.

The complete dataset contains 30000 samples, however, only 3000 samples were
used for training and validation purposes. The development of the methods were
supported by using an auxiliary test set containing 1500 unseen samples.
Additionally, a new set of microstructure images containing a mix of circular
and rectangular inclusions within a single RVE and, further, structures
composed of ellipsoidal inclusions are used as benchmark, constituting truly
unseen data.

##  Appendix B Neural network architectures

Table 5: The architecture of the Bayesian neural network used for the regression during feature selection is shown. The abbreviation BN stands for batch normalization. Every Dense connection used the selu activation function, except for the output layer which uses the identity function. BNN architecture | Dense 45 | BN | Dense 32 | BN
---|---|---|---|---
Dense 25 | BN |  Dense $3\cdot 2$ |
Table 6: The chosen architecture of the generic Conv Net is shown. The convolutional layers are read as ’filter size/stride $\times n_{\rm channels}$’. Dense layers simply show the number of output neurons of the layer. Every convolutional and dense layer have a selu activation functions, except for the last output layer which uses the identity function.  Convolutional layers | AvgPool 2/2 |  conv 11/4$\cdot 32$ | MaxPool 2/2 |  conv 7/3$\cdot$32
---|---|---|---|---
MaxPool 2/2 |  conv 5/3$\cdot 64$ | MaxPool 2/2 |  conv 3/2$\cdot 64$
MaxPool 2/2 |  conv 3/2$\cdot 96$ | MaxPool 2/2 | Flatten
Regressor | Dense 100 | BatchNorm | Dense 70 | BatchNorm
Dense 50 | BatchNorm | Dense 30 | Dense 3

Figure 30: The full layout of the deployed deep inception modules is shown.
Both modules are connected in parallel to the input image, as well as the
following dense regressor. These deep inception modules have been used for the
final hybrid model, which contained four hidden layers of (32-32-16-16) hidden
neurons using the selu activation function with batch normalization in
between.

##  Appendix C Full error plots

Figure 31: The R2 plots for each of the relevant models is shown for the
entire benchmark set and each component of the heat conduction tensor.
Relative errors are color coded, such that darker scatters indicate better
predictions. Figure 32: The predictions for all trained phase contrasts for
the variable model are given alike to Fig. 24. Not all available data points
have been plotted as scatters but are considered for the error bar computation
at the bottom.

## References

  *   * (1) Miehe, C.: Strain-driven homogenization of inelastic microstructures and composites based on an incremental variational formulation. International Journal for Numerical Methods in Engineering 55, 1285–1322 (2002). https://doi.org/10.1002/nme.515
  * (2) Beyerlein, I., Tomé, C.: A dislocation-based constitutive law for pure Zr including temperature effects. International Journal of Plasticity 24(5), 867–895 (2008). https://doi.org/10.1016/j.ijplas.2007.07.017
  * (3) Keshav, S., Fritzen, F., Kabel, M.: FFT-based Homogenization at Finite Strains using Composite Boxels (ComBo). arXiv preprint arXiv:2204.13624 (2022)
  * (4) Kumar, A., Nguyen, L., DeGraef, M., Sundararaghavan, V.: A Markov random field approach for microstructure synthesis. Modelling and Simulation in Materials Science and Engineering 24(3), 035015 (2016). https://doi.org/10.1088/0965-0393/24/3/035015
  * (5) Cang, R., Li, H., Yao, H., Jiao, Y., Ren, Y.: Improving direct physical properties prediction of heterogeneous materials from imaging data via convolutional neural network and a morphology-aware generative model. Computational Materials Science 150, 212–221 (2018). https://doi.org/10.1016/j.commatsci.2018.03.074
  * (6) Seibert, P., Raßloff, A., Ambati, M., Kästner, M.: Descriptor-based reconstruction of three-dimensional microstructures through gradient-based optimization. Acta Materialia, 117667 (2022). https://doi.org/10.1016/j.actamat.2022.117667
  * (7) Lißner, J., Fritzen, F.: Data-Driven Microstructure Property Relations. Mathematical and Computational Applications 24(2), 57 (2019). https://doi.org/10.3390/mca24020057
  * (8) Ford, E., Maneparambil, K., Rajan, S., Neithalath, N.: Machine learning-based accelerated property prediction of two-phase materials using microstructural descriptors and finite element analysis. Computational Materials Science 191, 110328 (2021). https://doi.org/10.1016/j.commatsci.2021.110328
  * (9) Brough, D.B., Wheeler, D., Kalidindi, S.R.: Materials knowledge systems in python—a data science framework for accelerated development of hierarchical materials. Integrating materials and manufacturing innovation 6(1), 36–53 (2017). https://doi.org/10.1007/s40192-017-0089-0
  * (10) Marshall, A., Kalidindi, S.R.: Autonomous development of a machine-learning model for the plastic response of two-phase composites from micromechanical finite element models. JOM 73(7), 2085–2095 (2021). https://doi.org/10.1007/s11837-021-04696-w
  * (11) Farizhandi, A.A.K., Mamivand, M.: Processing Time, Temperature, and Initial Chemical Composition Prediction from Materials Microstructure by Deep Network for Multiple Inputs and Fused Data. Materials & Design, 110799 (2022)
  * (12) Lubbers, N., Lookman, T., Barros, K.: Inferring low-dimensional microstructure representations using convolutional neural networks. Physical Review E 96(5), 052111 (2017). https://doi.org/10.1103/PhysRevE.96.052111
  * (13) Gayon-Lombardo, A., Mosser, L., Brandon, N.P., Cooper, S.J.: Pores for thought: generative adversarial networks for stochastic reconstruction of 3D multi-phase electrode microstructures with periodic boundaries. npj Computational Materials 6(1), 1–11 (2020). https://doi.org/10.1038/s41524-020-0340-7
  * (14) Liu, X., Zhou, S., Yan, Z., Zhong, Z., Shikazono, N., Hara, S.: Correlation between microstructures and macroscopic properties of nickel/yttria-stabilized zirconia (Ni-YSZ) anodes: Meso-scale modeling and deep learning with convolutional neural networks. Energy and AI 7, 100122 (2022). https://doi.org/10.1016/j.egyai.2021.100122
  * (15) Yang, C., Kim, Y., Ryu, S., Gu, G.X.: Prediction of composite microstructure stress-strain curves using convolutional neural networks. Materials & Design 189, 108509 (2020). https://doi.org/10.1016/j.matdes.2020.108509
  * (16) Fernández, M., Fritzen, F.: On the generation of periodic discrete structures with identical two-point correlation. Proceedings of the Royal Society A 476(2242), 20200568 (2020)
  * (17) Fast, T., Kalidindi, S.R.: Formulation and calibration of higher-order elastic localization relationships using the MKS approach. Acta Materialia 59(11), 4595–4605 (2011)
  * (18) Mikut, R., Reischl, M.: Data mining tools. Wiley interdisciplinary reviews: data mining and knowledge discovery 1(5), 431–443 (2011)
  * (19) Tipping, M.E.: Bayesian inference: An introduction to principles and practice in machine learning. In: Summer School on Machine Learning, pp. 41–62 (2003). https://doi.org/10.1007/978-3-540-28650-9_3. Springer
  * (20) Szegedy, C., Liu, W., Jia, Y., Sermanet, P., Reed, S., Anguelov, D., Erhan, D., Vanhoucke, V., Rabinovich, A.: Going deeper with convolutions. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 1–9 (2015)
  * (21) Lißner, J.: 2d Microstructure Data. https://doi.org/10.18419/darus-1151. https://doi.org/10.18419/darus-1151
  * (22) Lißner, J.: Hybrid neural network codebase. GitHub (2023). https://github.com/J-lissner/hybrid_neural_network
  * (23) Torquato, S.: Random Heterogeneous Materials: Microstructure and Macroscopic Properties vol. 16. Springer, ??? (2013)
  * (24) Kendall, A., Gal, Y.: What uncertainties do we need in bayesian deep learning for computer vision? Advances in neural information processing systems 30 (2017)
  * (25) Depeweg, S.: Modeling epistemic and aleatoric uncertainty with bayesian neural networks and latent variables. PhD thesis, Technische Universität München (2019). https://mediatum.ub.tum.de/1482483
  * (26) Dillon, J.V., Langmore, I., Tran, D., Brevdo, E., Vasudevan, S., Moore, D., Patton, B., Alemi, A., Hoffman, M., Saurous, R.A.: Tensorflow distributions. arXiv preprint arXiv:1711.10604 (2017)
  * (27) Goan, E., Fookes, C.: Bayesian neural networks: An introduction and survey. In: Case Studies in Applied Bayesian Data Science, pp. 45–87. Springer, ??? (2020). https://doi.org/10.1007/978-3-030-42553-1_3
  * (28) Kabir, H.D., Khosravi, A., Hosen, M.A., Nahavandi, S.: Neural network-based uncertainty quantification: A survey of methodologies and applications. IEEE access 6, 36218–36234 (2018). https://doi.org/10.1109/ACCESS.2018.2836917
  * (29) Hunt, B.: A matrix theory proof of the discrete convolution theorem. IEEE Transactions on Audio and Electroacoustics 19(4), 285–288 (1971). https://doi.org/10.1109/TAU.1971.1162202
  * (30) Lu, B., Torquato, S.: Lineal-path function for random heterogeneous materials. Physical Review A 45(2), 922 (1992). https://doi.org/10.1103/PhysRevA.45.922
  * (31) Kalidindi, S.R., Niezgoda, S.R., Salem, A.A.: Microstructure informatics using higher-order statistics and efficient data-mining protocols. Jom 63(4), 34–41 (2011). https://doi.org/10.1007/s11837-011-0057-7
  * (32) Scheunemann, L., Balzani, D., Brands, D., Schröder, J.: Design of 3D statistically similar representative volume elements based on Minkowski functionals. Mechanics of Materials 90, 185–201 (2015). https://doi.org/10.1016/j.mechmat.2015.03.005
  * (33) Balzani, D., Scheunemann, L., Brands, D., Schröder, J.: Construction of two-and three-dimensional statistically similar RVEs for coupled micro-macro simulations. Computational Mechanics 54(5), 1269–1284 (2014)
  * (34) Gholamalinezhad, H., Khosravi, H.: Pooling methods in deep neural networks, a review. arXiv preprint arXiv:2009.07485 (2020)
  * (35) Vernon, D.: Machine Vision: Automated Visual Inspection and Robot Vision. Prentice-Hall, Inc., ??? (1991)
  * (36) O’Shea, K., Nash, R.: An introduction to convolutional neural networks. arXiv preprint arXiv:1511.08458 (2015)
  * (37) Basheer, I.A., Hajmeer, M.: Artificial neural networks: fundamentals, computing, design, and application. Journal of microbiological methods 43(1), 3–31 (2000). https://doi.org/10.1016/S0167-7012(00)00201-3
  * (38) LeCun, Y., et al.: Generalization and network design strategies. Connectionism in perspective 19(143-155), 18 (1989)
  * (39) Schubert, S., Neubert, P., Pöschmann, J., Protzel, P.: Circular convolutional neural networks for panoramic images and laser data. In: 2019 IEEE Intelligent Vehicles Symposium (IV), pp. 653–660 (2019). IEEE
  * (40) Kauderer-Abrams, E.: Quantifying translation-invariance in convolutional neural networks. arXiv preprint arXiv:1801.01450 (2017)
  * (41) DeVries, T., Taylor, G.W.: Improved regularization of convolutional neural networks with cutout. arXiv preprint arXiv:1708.04552 (2017)
  * (42) Loshchilov, I., Hutter, F.: Decoupled weight decay regularization. arXiv preprint arXiv:1711.05101 (2017)
  * (43) Abadi, M., Agarwal, A., Barham, P., Brevdo, E., Chen, Z., Citro, C., Corrado, G.S., Davis, A., Dean, J., Devin, M., Ghemawat, S., Goodfellow, I., Harp, A., Irving, G., Isard, M., Yangqing, J., Jozefowicz, R., Kaiser, L., Kudlur, M., Levenberg, J., Mané, D., Monga, R., Moore, S., Murray, D., Olah, C., Schuster, M., Shlens, J., Steiner, B., Sutskever, I., Talwar, K., Tucker, P., Vanhoucke, V., Vasudevan, V., Viégas, F., Vinyals, O., Warden, P., Wattenberg, M., Wicke, M., Yu, Y., Zheng, X.: TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems. Software available from tensorflow.org (2015). https://www.tensorflow.org/
  * (44) Feir-Walsh, B.J., Toothaker, L.E.: An empirical comparison of the ANOVA F-test, normal scores test and Kruskal-Wallis test under violation of assumptions. Educational and Psychological Measurement 34(4), 789–799 (1974). https://doi.org/10.1177/001316447403400406
  * (45) Guyon, I., Weston, J., Barnhill, S., Vapnik, V.: Gene selection for cancer classification using support vector machines. Machine learning 46(1), 389–422 (2002). https://doi.org/10.1023/A:1012487302797
  * (46) Geurts, P., Ernst, D., Wehenkel, L.: Extremely randomized trees. Machine learning 63(1), 3–42 (2006). https://doi.org/10.1007/s10994-006-6226-1
  * (47) Leuschner, M., Fritzen, F.: Fourier-accelerated nodal solvers (FANS) for homogenization problems. Computational Mechanics 62(3), 359–392 (2018)
  * (48) Hashin, Z., Shtrikman, S.: A variational approach to the theory of the elastic behaviour of polycrystals. Journal of the Mechanics and Physics of Solids 10(4), 343–352 (1962). https://doi.org/10.1016/0022-5096(62)90005-4