Hybrid machine-learned homogenization: Bayesian data mining and
  convolutional neural networks
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.
Julian Lißner, Felix Fritzen
# Hybrid machine-learned homogenization: Bayesian data mining and convolutional neural networks

###### 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.

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 [1; 2]. 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 [3]. 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.

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 [4, 5, 6], 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 [7, 8]. 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) [9], using the principal scores as input for a supervised machine learned regressor to conduct the microstructure property linkage, ranging from polynomial regression [10] to artificial neural networks [7, 11]. Lately, the usage of convolutional neural networks (Conv Nets) has gained popularity in both outlined research fields [12, 13]. Conv Nets are used to directly predict the effective material response [14], or even in combination with the PCA where the reduced representation of the target values is used to predict stress strain relationships [15]. 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 [7] 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., [8, 16]. Improvements beyond the 2PCF have been attempted by considering partial higher order correlation functions, e.g., [17]. 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 [18]. This constitutes the first block of this work (section 2.2): We systematically categorize the underlying data by using Bayesian neural networks [19]. 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 [20] 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 [7]. 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 [21] used for development, training and validation purposes, as well as the python code [22] are made freely available.
## 2 Data and Methods

### Data overview

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 [23] 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 RVE in the following. The target values of the machine learning algorithm are the components of the effective heat conduction tensor \(\frac{\bar{\kappa}}{\equiv}\). One Component of \(\frac{\bar{\kappa}}{\equiv}\) 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{\mathbf{\bar{\mathbf{\varepsilon}}}}}=\begin{bmatrix}\bar{ \kappa}_{11}\\ \bar{\kappa}_{22}\\ \sqrt{2}\bar{\kappa}_{12}\end{bmatrix}\,. \tag{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{\mathbf{\bar{\mathbf{\varepsilon}}}}}=f(\text{RVE}) \tag{2}\]

where the machine learned model \(f(\cdot)\) directly operates on the image data of the microstructured material (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.

### Bayesian Modeling via artificial neural networks

#### 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

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.
confidence indicating the likelihood for low or high prediction errors during inference is advantageous.

Bayesian modeling approaches uncertainty quantification via _aleatoric_ and _epistemic_ uncertainties [24, 25]. 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 [26] 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})\,, \tag{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 [25]:

\[\Phi_{i}\big{(}y_{i},\mathcal{N}(\hat{\mu}_{i},\hat{\sigma}_{i}^{2})\big{)}=- \log\big{(}\frac{1}{2}p(y_{i}|\mathcal{N}(\hat{\mu}_{i},\hat{\sigma}_{i}^{2})+ s)\big{)}, \tag{4}\]

where \(p(y|\mathcal{N}(\hat{\mu}_{i},\hat{\sigma}_{i}^{2}))\) is the probability of observing \(y_{i}\) given the current distribution \(\mathcal{N}(\hat{\mu}_{i},\hat{\sigma}_{i}^{2})\). Since the logarithm quickly tends to \(-\infty\) for very small values in \(p\left(y|\mathcal{N}(\hat{\mu}_{i},\hat{\sigma}_{i}^{2})\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., [19, 27], presenting a more generous outline to bayesian modeling and [28] who reviews different

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.
modeling methods in uncertainty quantification. An example application using bayesian modeling is presented in [24].

#### 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{\underline{\mbox{\it k}}}^{\mathrm{S}}\)

\[\underline{\underline{\mbox{\it k}}}^{\mathrm{S}}=\frac{1}{2}\left[-1,\;0,\; 1\right] \tag{5}\]

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

\[\begin{array}{rcl}\underline{\underline{\mbox{\it f}}}=\underline{\underline{ \mbox{\it s}}}*\underline{\underline{\mbox{\it k}}}^{\mathrm{S}}\to f_{i}&=& \sum\limits_{j=-1}^{1}s_{i+j}\ k_{j+2}^{\mathrm{S}}\\ &=&\frac{1}{2}(s_{i+1}-s_{i-1})\,\end{array} \tag{6}\]

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

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

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 [29], i.e.

\[\underline{\underline{\mbox{\it F}}}=\underline{\underline{\mbox{\it A}}}* \underline{\underline{\mbox{\it k}}}=\mathcal{F}^{-1}\big{(}\,\mathcal{F}( \underline{\underline{\mbox{\it A}}})\cdot\mathcal{F}(\underline{\underline{ \mbox{\it k}}})\big{)}\,. \tag{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.

Similar to the motivation in 1d is the resulting feature map \(\underline{\underline{\mbox{\it F}}}\) depends on the kernel \(\underline{\underline{\mbox{\it k}}}\), even when applying it to the same image \(\underline{\underline{\mbox{\it 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.

Figure 3: The convolution of one RVE with two different kernels is shown. The superscripts S, R of \(\underline{\underline{\mbox{\it k}}}\) represent the _Sobel_ filter and the searched _rectangle_, respectively. The feature map \(\underline{\underline{\mbox{\it F}}}\) to the right differs significantly when applying different kernels to the same image.
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.

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

\[\|e_{\mathrm{rel}}\|_{2}=\frac{\|\bar{\kappa}-\hat{\kappa}\|_{2}}{\|\bar{ \kappa}\|_{2}} \tag{8}\]

with the target value \(\bar{\kappa}\) and the prediction \(\hat{\kappa}\leftarrow f(\underline{\mathrm{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 [30]. The lineal path function has been previously adopted in multiple studies, e.g., being used as a feature [31], or as a descriptor to generate statistically similar RVEs [32, 33]. 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.

The newly proposed band features are computed with a single convolution per direction \(d\) as

\[p_{d}^{\mathrm{I}}=\max(\underline{b}_{d}*\underline{\mathrm{RVE}}) \tag{9}\]

where \(\underline{b}_{d}\) denotes the detection band in the indexed direction \(d\). Note that each detector \(\underline{b}_{d}\) is normalized such that \(p_{d}^{\mathrm{I}}\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

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.

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
by

\[\begin{array}{l}\widehat{p}_{d}^{\mathrm{M}}\ =1-\min(\underline{b}_{d}*\underline{ \mathrm{RVE}})\\ \ =\max\left(\underline{b}_{d}*(1-\underline{\mathrm{RVE}})\right),\end{array} \tag{10}\]

where \(1-\underline{\mathrm{RVE}}\) is a phase inversion of the image in the bi-phasic setting. Generally, \((p_{d}^{\mathrm{I}}-1)\neq p_{d}^{\mathrm{M}}\), 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.

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}=\mathrm{mean}\Big{(}\ \mathrm{sum}(\ \underline{\mathrm{RVE}},\ \mathrm{axis =}i)\geq 1\Big{)} \tag{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_[34] (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 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 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.

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{b}_{d}\) in b) and c) are drawn in their approximate position for feature determination.
Taking the mean of \(\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{c})\) and skewness \(\widetilde{\mu}_{3}^{f}(\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_{\text{cells}}}\left\{\begin{array}{l@{\quad}l}\# \big{(}\quad&f_{\text{cell}}\,<\quad\epsilon\quad\big{)}\\ \#\big{(}\quad\epsilon\,\leq f_{\text{cell}}\,<\,\frac{1}{3}-\epsilon\,\, \big{)}\\ \#\big{(}\,\,\frac{1}{3}-\epsilon\,\leq f_{\text{cell}}\,<\,\frac{2}{3}- \epsilon\,\,\big{)}\\ \#\big{(}\,\,\frac{1}{3}-\epsilon\,\leq f_{\text{cell}}\,<\,1-\epsilon\,\, \big{)}\\ \#\big{(}\,\,1-\epsilon\,\leq f_{\text{cell}}\,\,\big{)}\\ \end{array}\right., \tag{12}\]

where \(\#(\bullet)\) denotes the operation counting the number of cells fulfilling the condition, \(f_{\text{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.

Another feature descriptor is proposed to further quantify the size and sphericity of the inclusions through the surface information, since the ratio of \(\frac{\text{inclusion\,area}}{\text{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 [35] 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}=|\underline{k}^{e}*\underline{\text{RVE}}|, \qquad e\in\{1,2,3,4\} \tag{13}\]

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

\[\underline{\underline{k}}^{1} =\left[-1\,\,0\,\,1\right]\,, \underline{\underline{k}}^{2} =\left[\begin{matrix}0&-0.5&-1\\ 0.5&0&-0.5\\ 1&0.5&0\end{matrix}\right]\,, \tag{14}\] \[\underline{\underline{k}}^{3} =\left[\begin{matrix}-1&-0.5&0\\ -0.5&0&0.5\\ 0&0.5&1\end{matrix}\right]\,, \underline{\underline{k}}^{4} =\left[\begin{matrix}-1\\ 0\\ 1\end{matrix}\right]\,.\]

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{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 are 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

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.

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.
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.

### Convolutional neural networks

Convolutional neural networks (Conv Nets) [36] are derived from feed forward neural networks [7, 37] to make the classification/regression of high-resolution images tractable via weight sharing [38]. 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}* \underline{\underline{A}}_{l}+b_{l})\,, \tag{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.

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_{\text{in}}\times n_{\text{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.

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

\begin{table}
\begin{tabular}{|l|c|c|c|c|} \hline feature name & number of features & number of hyperparameters & captured phenomena & related equation \\ \hline 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) \\ \hline \end{tabular}
\end{table}
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.

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).

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.
on every other position, and a downscaling of approximately factor two is achieved along each dimension (Fig. 12).

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 [39] and is made publicly available in [22], where also the implementation for the data augmentation scheme is found.

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 [40]. One study has successfully applied

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.

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.
cutout, which goes as far as to simply gray out larger regions of the input image during training [41]. 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.

### 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 [20]. 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.

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, [20] 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*, 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).

Footnote *: 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)

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

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

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.
individually achieves a downsampling of a factor of 40, and the channels of each branch are concatenated at the end.

### The hybrid neural network

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 [7] 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}}_{\text{vol}}+\Delta( \tilde{\underline{\kappa}}_{\text{features}}+\tilde{\underline{\kappa}}_{ \text{ConvNet}})\,, \tag{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.

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)).

## 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 R\({}^{2}\) plot in the appendix C3. 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 [42], early stopping and a constant learning rate with default hyperparameters. The model training has been implemented using the tensorflow API [43], and

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_{\text{channels}}\).
the code for training and post processing can be found in [22].

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_{\mathrm{samples}}\cdot n_{\kappa}}\sum_{i=1}^{n_{ \mathrm{samples}}}\sum_{j=1}^{n_{\kappa}}\frac{(y_{ij}-\hat{y}_{ij})^{2}}{y_{ ij}^{2}}}\,, \tag{17}\]

which coincides with the error measure introduced in (8) and will be denoted with rel. \(\sqrt{\mathrm{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).

### 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 B1. 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 [7], 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 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.

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.
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 [44]. The deployed wrapper methods were the recursive feature elimination (RFE)[45] and a method deploying random forests [46]. Each of these methods found a slightly different order, presented as indices, given in arrays:

\[\begin{array}{l}\text{Pearson ranking:}\\ \left[\begin{array}{
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.

In a final step, multiple BNNs (with the same architecture tab. 1) 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.

### Convolutional and hybrid neural networks

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 B2. 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

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.

Figure 22: A density plot is shown which relates the aleatoric uncertainty to the rel. \(\sqrt{\mathrm{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.
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 B2.

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 over-fitting 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

\begin{table}
\begin{tabular}{|l|c|c|c|c|c|c|} \hline \multirow{2}{*}{Model type} & \multirow{2}{*}{mean relative error [\%]} & \multicolumn{2}{c|}{real.} & \multicolumn{1}{c|}{MSE} & \multicolumn{1}{c|}{trainable} \\  & & \(\bar{\kappa}_{11}\) & \(\bar{\kappa}_{22}\) & \(\bar{\kappa}_{11}\) & \(\bar{\kappa}_{22}\) & \multicolumn{1}{c|}{[\%]} & \multicolumn{1}{c|}{[-]} \\ \hline \hline 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 [7] & 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 \\ \hline 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 (0.1) & 1.20 (1.02) & 1.00 (0.81) & 0.91 (0.79) & 1.71 (1.38) & 5.33 (3.11) & 4 997 \\ \hline 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 & 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 \\ + 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 & \(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 \\ +data augmentation & & & & & & \\ inception Conv Net & & & & & & \\ + vol bypass & 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 \\ +data augmentation & & & & & & \\ + vol bypass & 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 \\ +data augmentation & & & & & & \\ \hline 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 \\ \hline \hline \end{tabular}
\end{table}
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.

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.
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.

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

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.

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.

#### 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{\mathrm{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. 14.

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.
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

\begin{table}
\begin{tabular}{|l|c|c|c c|c c|} \hline \begin{tabular}{c} RVE image \\ configuration \\ \end{tabular} & \multicolumn{2}{c|}{mean relative error [\%]} & \multicolumn{2}{c|}{median relative error [\%]} & \multicolumn{2}{c|}{
\begin{tabular}{c} relative \\ MSE [\(\cdot\)] \\ \end{tabular} } & \multicolumn{2}{c|}{\(\begin{bmatrix}\text{MSE}\\ 10^{-5}\) [\(\cdot\)] \\ \end{tabular} } \\  & \(\bar{\kappa}_{11}\) & \(\bar{\kappa}_{22}\) & \(\bar{\kappa}_{11}\) & \(\bar{\kappa}_{22}\) & & \\ \hline original & 0.730 & 0.742 & 0.552 & 0.562 & 1.027 & 1.915 \\ \hline rotated & 0.759 & 0.751 & 0.564 & 0.574 & 1.045 & 1.982 \\ \hline \end{tabular}
\end{table}
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.

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 \(\underline{\hat{\kappa}}\) as well as the actual target values \(\underline{\hat{\kappa}}\) are shown below each RVE.
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 [7] 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

\begin{table}
\begin{tabular}{|c|c c|c c c|c c|} \hline Phase & mean relative error [\%] & \multicolumn{2}{c|}{mean absolute error [-]} & \multicolumn{2}{c|}{relative} & \multicolumn{1}{c|}{MSE} \\ contrast & \multicolumn{2}{c|}{\(\bar{\kappa}_{\bullet}>0.2\)} & \multicolumn{2}{c|}{\(\bar{\kappa}_{\bullet}\leq 0.2\)} & \multicolumn{2}{c|}{\(\cdot 10^{-5}[-]\)} \\  & \(\bar{\kappa}_{11}\) & \(\bar{\kappa}_{22}\) & \(\bar{\kappa}_{11}\) & \(\bar{\kappa}_{22}\) & \(\bar{\kappa}_{12}\) & \\ \hline \hline
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 \\ \hline \end{tabular}
\end{table}
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. 14). The error measures are given at the discrete sampled phase contrasts the model was trained on.

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.
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.

AcknowledgmentsFunded 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 FR2702/8 - 406068690 and DFG-FR2702/10 - 517847245. We acknowledge the support by the Stuttgart Center for Simulation Science (SimTech).

## Declarations

* 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) [21]
* Code availability The code will be made publicly available in [22] upon publication.

## Appendix A Data overview

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

\[A\in\mathbb{A}^{n\times n},\quad A_{ij}=A(\mathbf{x}_{ij}),\quad\mathbf{x}_{ij}\in \Omega, \tag{11}\]

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,\ldots,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 conductivity2. 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

Footnote 2: Equivalently, this can be interpreted as a prediction of the permeability of the material which has the same underlying mathematical structure.

\[\kappa_{\rm M}=1\,\frac{\rm J}{\rm sKm},\quad\quad\quad\kappa_{\rm I}=\frac{ \kappa_{\rm M}}{R}\,, \tag{12}\]

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 [47] yield the symmetric effective heat conduction tensor \(\tilde{\underline{\tilde{\kappa}}}\) as

\[\tilde{\underline{\tilde{\kappa}}}=\begin{bmatrix}\tilde{\kappa}_{11}&\tilde{ \kappa}_{12}\\ \tilde{\kappa}_{12}&\tilde{\kappa}_{22}\end{bmatrix}\,, \tag{11}\]

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

\[\tilde{\underline{\tilde{\kappa}}}=\begin{bmatrix}\tilde{\kappa}_{11}\\ \tilde{\kappa}_{22}\\ \sqrt{2}\tilde{\kappa}_{12}\end{bmatrix}\,. \tag{12}\]

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

\[\tilde{\underline{\tilde{\kappa}}}=\frac{\tilde{\underline{\tilde{\kappa}}}}{ \kappa_{\text{M}}}\,, \tag{13}\]

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

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

\[\overline{X}=\frac{1}{n}\sum_{j=1}^{n}X_{j},\;\;\;\widetilde{X}_{i}=\frac{ \sqrt{n-1}(X_{i}-\overline{X})}{\sqrt{\sum_{j=1}^{n}(X_{j}- \overline{X})^{2}}}\,. \tag{14}\]

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. 11). 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

## Appendix C Full error plots
