{'fraction_non_alphanumeric': 0.028704255937361412, 'fraction_numerical': 0.008013335168448257, 'mean_word_length': 4.686260869565217, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 23, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A heteroencoder architecture for prediction of failure locations in porous metals using variational inference

In this work we employ an encoder-decoder convolutional neural network to predict the failure locations of porous metal tension specimens based only on their initial porosities. The process we model is complex, with a progression from initial void nucleation, to saturation, and ultimately failure. The objective of predicting failure locations presents an extreme case of class imbalance since most of the material in the specimens do not fail. In response to this challenge, we develop and demonstrate the effectiveness of data- and loss-based regularization methods. Since there is considerable sensitivity of the failure location to the particular configuration of voids, we also use variational inference to provide uncertainties for the neural network predictions. We connect the deterministic and Bayesian convolutional neural networks at a theoretical level to explain how variational inference regularizes the training and predictions. We demonstrate that the resulting predicted variances are effective in ranking the locations that are most likely to fail in any given specimen.


Introduction
Additive manufacturing (AM) [1] and topology optimization [2,3] have opened a new era for design of complex and optimized components; however, fundamental issues with the porosity of metal AM parts persist. Experimental assessment tools exist, such as computed tomography (CT) and other non-destructive screening methods, but they need to be coupled to models for quantitative predictions of performance and reliability. Typically there is large uncertainty in predictions since failure is highly sensitive to small flaws, such as voids, and is, by its nature, an unstable process. Since failure is usually localized, indications of critical flaws can be somewhat imprecise and yet useful. Direct simulation of ductile failure is computationally expensive due the complexity of the models and the issues with material stability. Here we design a convolutional neural network to efficiently predict failure locations given initial porosity at a suitable resolution. Since the output image is of different character than the input image, as opposed to the widely-employed autoencoder, the proposed architecture is more accurately described as a heteroencoder [4,5].
Convolutional neural networks (CNNs) [6] and similar neural networks have been used for a variety of related tasks in physics and engineering. For instance, Garland et al. [7] used a CNN as a screening tool for predicting the mechanical properties of additively manufactured lattices. Li et al. [8] employed an encoder-decoder CNN to predict the state of a reaction-diffusion system. Nie et al. [9] used two different CNN architectures to predict the stress field of cantilever beams. CNNs have been widely employed in homogenization and related multi-fidelity applications. Frankel et al. [10,11,12] have employed pixel and graph based convolutional neural networks for homogenization and field predictions. Teichert and Garikipati [13] have used multi-fidelity learning and surrogate optimization as alternates to phase field modeling to predict precipitate morphologies in alloys. Zhang and Garikipati [14] also have explored dense neural networks (DNNs) and CNNs with multiresolution learning of the free energy and homogenized the nonlinear elastic response of evolving microstructures. Teichert et al. [15,16,17] have developed a class of integrable deep neural networks and active learning in the context of scale bridging in materials physics. Johnson et al. [18] used a CNN to predict the peak load sustained by simulated tensile specimens with different porosity realizations. They also demonstrated that the network trained to J2 plasticity response was predictive outside the training conditions. Zhu and Zabaras [19] employed Bayesian CNNs to predict the response of stochastic permeability data. They show that Stein variational gradient descent can outperform standard variational inference (VI) techniques and Gaussian processes. CNNs have also been applied in the more general context of solving partial differential equations (PDEs). Zhang and Garikipati [20] have employed both deterministic and Bayesian neural networks for the direct solution of elliptic PDEs. Physics informed neural networks (PINNs) [21] were developed as general framework for learning from and predicting field data from PDEs. Muhammad et al. [22] used a PINNslike framework to predict strain fields of experimental tensile specimens using processing and local microstructural information as inputs. Digital image correlation data was used to provide strains for training. Particularly related to the present work, Buda et al. [23] discussed and investigated methods to alleviate class imbalance with CNNs.
The contributions of this work to predicting failure with CNNs are manyfold. We develop a means to ameliorate class imbalance inherent in predicting single/few points of failure in a larger domain through data transformations and loss re-weightings. We also recast the categorical problem of failure/no-failure at a voxel level to a regression problem and then use precision and recall as scoring metrics that encode the usefulness of predictions. We relate perturbative sensitivity analysis of the system response to the CNN and the Bayesian CNN (BCNN) predictions. We connect the formulation of the calibration of deterministic and Bayesian convolutional neural networks at a theoretical level to explain how variational inference regularizes training and prediction. Furthermore we provide concrete results of the value of using variational inference (VI) with the CNN to provide confidence in the predictions and indications of the epistemic uncertainty of the model. We then use the predicted mean and variance of the resulting BCNN to effectively rank the likely failure locations.
In Sec.2 we describe the phenomenology of the failure process of ductile porous metals as represented by a well-calibrated physics-based model. Then, in Sec.3, we develop the neural network architecture as well as the data transformations and loss manipulations needed for useful predictions. In Sec.4 we demonstrate the performance of the proposed CNN. Then, in Sec. 5, we augment the CNN architecture with VI and connect the Bayesian augmentation to the deterministic model. In Sec. 6 we illustrate the advantages of the Bayesian network. We conclude, in Sec. 7, with a summary of our findings and some additional observations in the context of future work.

Ductile failure of porous metals
To emulate a real failure process we employed a traditional physics-based model calibrated to experimental data consisting of porosity imaging and stress-strain response curves. The physical model has two levels of porosity: (a) larger pores that are represented explicitly in the mesh, and (b) smaller pores represented implicitly in the viscoplastic constitutive model with damage evolution. The distinction between the two is the explicit pores were visible with computed tomography (CT) [24,25], while the implicit porosity was inferred in the calibration process [26]. The initial explicit pores ϕ(X), initial background void density φ(X, 0), and loading ε(t) interact to create an evolving damage field φ(X, t) and, ultimately, failure locations Φ(X). As in the high-throughput experiments used to generate the calibration data, we simulate the rectangular gauge section of a tensile specimen loaded to failure.

Constitutive Model
We employ a previously calibrated model of AM 17-4 PH stainless steel [26] to represent the material behavior of the porous tension specimens in this study. This plasticity model is affected by accumulation of damage that results in local ductile failure and ultimately structural instability with sufficient loading in the 1 mm × 1 mm × 4mm porous tension specimens.
In the constitutive model, stress σ follows a linear elastic rule: where C is the isotropic elastic modulus tensor with components which depend on Young's modulus E and Poisson's ratio v, is the total strain, and p is the plastic strain. The evolution of the void fraction, φ, effects a reduction of stiffness of the material. The mechanical behavior evolves according to a coupled set of ordinary differential equations. The plastic strain p increases according to: with p (t = 0) = 0 and where Y is the yield strength, s = dev σ is the deviatoric stress, and σ vm ≡ 3 2 s · s is the von Mises stress. The evolution of the isotropic hardening variable κ is governed by: with κ(t = 0) = κ 0 and where H is the hardening modulus, R is the recovery coefficient, and˙ p = 2 3˙ p ·˙ p is the equivalent plastic strain rate. The damage due to implicit void volume fraction φ evolves according to: which includes nucleation and growth of voids, where the damage exponent m, and the volume of newly nucleated voids ν 0 , are additional model parameters. To model nucleation, the void concentration η, evolves according to: where p = 1/3 σ · I is the pressure, J 2 = 1/2 tr s 2 = 1/3 σ 2 vm and J 3 = 1/3 tr s 3 . Lastly, once the void fraction φ exceeds a threshold φ max the material is considered completely failed and acts as an explicitly represented pore.
The calibrated parameters are given in Table 1 and additional details can be found in Ref. [26].

Generation of porosity realizations
As in our previous work [26] we utilized a 20×20×80 structured mesh of 8 node hexahedral elements. On this grid we generated the initial explicit porosity field via a Karhunen-Loève (KL) process with ≈ 12,000 modes calibrated to the average porosityφ = 0.0008 and spatial power-exponential correlation function fitted to the empirical correlation information as extracted from the available CT scans. Since the observed porosity is confined to clusters that consist of only a few voxels, the correlation length for the porosity process is short. Hence many modes are needed to resolve the high frequencies and capture 99.9% of the energy. Full details are given in Ref. [26].

Boundary value problem
Tensile loading was imposed with minimal displacement boundary conditions on the top (high z) and bottom (low z) surfaces of the specimen such that these surfaces remained planar but could change dimension laterally. A nominally quasi-static strain rate,ε = 0.002/s, was applied through these boundary conditions until structural failure.

Phenomenology of the failure process
Damage evolution leading to failure progresses through three stages in these specimens. Fig. 1 shows the response of a representative sample. In stage I (ε ≈ [0:2]%) there is a smooth rise of maximum damage in the sample from the background ϕ 0 =0.008 with no material failure. In stage II (ε ≈ [2:4]%) there is a plateau of maximum damage at φ/φ max ≈ 0.6 as sites nucleate damage, harden and then alternate sites start to damage. This stage ends when all the likely/vulnerable sites saturate. In stage III (ε ≈ [4:5]%) a single most vulnerable site (or in some rare cases a few clustered sites) reaches failure and the specimen quickly fails with a horizontal ductile fracture surface. This progression has some qualitative similarity with the well-known phenomenon of the evolution of Lüder bands [27]. The location along the tensile axis, i.e. , the z-location, of the first element(s) to reach failure is highly correlated with the ultimate z-location of failure as measured at the final load step. In fact, the difference in z-location of the first and final failure is on average only about 0.04 mm, or 1% of the specimen length or about one element, which indicates the first elements to fail are the most likely cause of ultimate failure. The histograms of damage for a sequence of strains are shown in Fig. 2. The panels of Fig. 2 for different realizations illustrate that the three stage progression is representative of the ensemble. Lastly, the bundle of response curves shown in Fig. 3 demonstrates that there is considerable variability within the ensemble. The progression of maximum damage in the samples over time is similar, but the rise, plateaus, and final failures vary considerably with the different pore configurations. Also once an element fails (at φ/φ max = 1 locally) the rise to structural failure is rapid for all samples but the strain at which it occurs, and the void volume fraction, φ, at which the structure ultimately fails vary substantially.

Sensitivity of failure to pore locations
Although the damage criterion used to define the failure location typically provides a single region of failure, there are multiple regions of high damage that compete in early stages of the loading (stage II). This suggests that small variations in the explicit porosity could lead to failure occurring in alternate regions of the specimen. To explore this sensitivity, we carry out one voxel perturbations of the porosity field. We randomly select a void voxel on the edge of a pore within a given distance from the original failure location and swap it with a randomly selected solid voxel anywhere in the specimen. This process slightly decreases the size of the pore in the neighborhood of original failure while maintaining the same overall porosity, but potentially alters the porous surface area. Loading the modified porosity structure results in a failure location that may or may not be close to the failure location of original structure which reflects the sensitivity of the failure location to the porosity field. This process was repeated 100 times for nine realizations and 500 times for three realizations (realizations (g), (h), and (l) in Fig. 4). The locus of failure locations are shown along with the explicit porosity of the unperturbed realization in Fig. 4. The porosity is projected onto the x-z-plane by summing the number of voided voxels along the y-axis and plotting the result as a grayscale value, where black indicates no pores. The failure locations are shown using colored squares, where multiple failures with the same x and z values are indicated using the color scale next to each plot. Some realizations, such as (a), (c), and (f) exhibit very little sensitivity to the random perturbation of a single void voxel. Others, e.g. (b), (i), (k), demonstrate two or three regions of failure, while realization (g) displays the most sensitivity, with failure occurring at several different locations within the specimen.

A convolutional neural network model of localized failure
Our goal in using neural networks is to predict likely failure locations, which is a binary categorization Φ(X) (failed/not failed), given the binary initial porosity field ϕ(X). Initial failure is defined as the first time at which the damage φ(X) in any element/voxel at X reaches φ max , at which point the element will start to act as a void and initiate structural instability. The 20×20×80 grid that discretizes the domain makes this a high, yet finite, dimensional problem. The emphasis on predicting the failure locations leads to stark class imbalance, where the failed to not failed ratio is on the order of 1:10 4 , since only a few of the elements fail in any particular mesh. The fact that a prediction of no failure anywhere scores nearly as well in mean error as a perfect prediction if there are only a few failures illustrates the class imbalance problem. So even relative accuracies of 0.1-0.01 % would be insufficient to reliably predict one element failing in a mesh of 32,000 elements.
As a first step to alleviating this issue, we recast this categorical prediction problem into a regression problem of predicting the underlying damage field φ(X, t fail ) at time of failure t fail . Although the damage fields are still non-smooth, the spatial correlations in the damage field make the regression task easier; in effect this recasting is a physically motivated regularization. Recall, the physics-based direct numerical simulations generate training data consisting of an input binary porosity field ϕ(X) mapped to an output damage field φ(X, t fail ) at the time of first failure t fail . Since the transition of φ → φ max in a particular location rapidly leads to complete structural failure, we have a natural, direct failure definition and a binary failure location field Φ(X) (output). So a high dimensional, two category, Φ ∈ {0, 1}, problem with data D = {ϕ, Φ} has a direct connection to the regression problem with D = {ϕ, φ}. To obtain a dataset D with 110,712 full-field instances, we generated 18,452 realizations and augmented this ensemble with 6-fold symmetry transforms (π/2 rotations about the tensile axis and reflections in each of the cardinal directions). In addition to the regularization of recasting the categorical problem, predicting a continuous, albeit irregular, field in the regression problem allows for data manipulation (Sec. 3.1) and loss re-weighting (Sec. 3.2) to further regularize the problem and promote accurate prediction of the failure locations. In this mode, failure locations are recovered post-prediction by applying a threshold to the predicted damage field.

Data transformations
As can be seen in Figs. 1 and 4 the damage fields can have isolated high damage voxels as well as more significant clusters that correlate with the ultimate failure locations. Based on this observation, we applied smoothing, contrast, and normalization transformations to condition the data and improve the CNN training without sacrificing the ability to predict failure locations.
We first applied a Gaussian filter to the damage field φ to obtain: with k 0 being the width of Gaussian kernel g, X being the distance between the current filter location and a neighboring location, and * being convolution. The use of a spatial smoothing filter was motivated by the fact that damage fields in metals are typically more smooth and regularly distributed than the direct numerical simulation results, since the viscoplastic model in Sec. 2.1 lacks an intrinsic length-scale. Furthermore, convolution such as Eq. (6) has the effect of removing high frequencies, which is compatible with the encoderdecoder network (Sec. 3.3) we employed that distills the full resolution images to a reduced dimensionality latent space. It is also consistent with the fact that the KL process (Sec. 2.2) that generated the porosity realizations is dominated by the low frequency components. For this dataset k 0 < 1.0 mm had a negligible smoothing effect, whereas k 0 > 2.0 mm removed important features (clusters of high damage); hence, kernel width k 0 = 1.5 mm provided a good balance between the two effects. Next, we applied the softmax function voxel-wise to each voxel I in the damage field φ of a realization:φ with c 0 being the contrast parameter that increases the difference between low and high values. The softmax transformation is commonly used in categorization and image processing. In our application, a larger c 0 results in higher contrast between high and low damage values, i.e. bringing the damage field closer to a binary field, and emphases the clusters of high damage. We used c 0 = 5.0 based on preliminary studies. Note that that the contrast and smoothing parameters could be treated as hyperparameters in the CNN training but the fact that extremes in their values lead to trivial predictions presents challenges which are left for future work. Finally, we normalize the resulting damage fieldφ per realization to obtain: As discussed in Sec. 2.4 and shown in Fig. 1, initially damage can rise at multiple sites due to nucleation and hardening, but ultimately only a most vulnerable site in a given porosity field will reach failure. Failure typically occurs at one or a few voxels with the maximum damage value φ max which are surrounded by a cluster of high damage. Since Gaussian smoothing results in a weighted damage value that accounts for the values in the surrounding neighborhood, ideally smoothing will: (a) remove isolated damage voxels that developed in stage II but did not lead to ultimate failure in stage III and (b) cluster the voxels surrounding the voxel that ultimately fails. The contrast transform operates locally but also suppresses low damage values and enhances high values. As shown in Fig. 5, the smoothing treatment emphasizes the region where failure ultimately occurs in stage III. The effectiveness of the different data transformations in promoting training and accurate predictions is discussed Sec. 4.2.
pores damage field transformed damage field

Weighted loss functions
As stated, predicting failure voxel locations with a convolutional neural network from porosity realizations can be naturally recast as a damage field regression problem. Hence, we use the mean squared error (MSE) of the predicted damage field as the optimization objective function for the regression problem: where φ = {φ I } is a vector of damage values at voxels for a specific sample, N s is the number of samples in the training set, N v the number of voxels in a damage field realization, and φ s ,φ s are the true and predicted fields on the grid for sample s, respectively. Since the damage field at the pore locations is not defined, we omit voxel locations corresponding to the initial voids in Eq. (9).
A primary issue with this standard loss function in the present context is that the damage fields are unbalanced in terms of the number of voxels with small damage values versus the number of voxels where there is significant damage. This discrepancy causes a CNN trained to this loss to learn the mapping for the low damage regions accurately at the expense of the high damage regions which are of greater importance. The distribution of voxels and regression errors depicted as a function of damage value in Fig. 6a  The frequency normalized error of the trained CNN predictions in Fig. 6b illustrate the effect damage value frequency has on the solution found by the optimizer. Gradient descent finds a local minimum where the reconstruction error is small around the highest frequency damage value and increases quickly away from this value. This can be understood by partitioning the terms of the MSE loss into two bins according to the true damage voxel value: where |D low | |D high | and D low contains all the relatively undamaged voxels (those below a threshold ≈ φ max ) and D high the few voxels that approach failure at φ max . Hence, given relatively uniform errors across each sample, the contribution to the gradient of the predominant D low data is likely to be significantly larger than the contribution of D high . This leads the optimizer to "overlearn" the high frequency low damage values.
This issue in our regression problem is related to class-imbalance in classification problems [28,29,30,23,31] where training data is unevenly distributed among the classes leading to solutions which are inaccurate for classes with fewer training examples. Common techniques to address this issue are loss-based techniques, where the loss function is modified with a weighting scheme to balance contributions from different classes [23], and data-based techniques, where sampling strategies are used to balance the training data [23]. In this work, we use the loss weighting approach as the imbalance exists in each field realization so that common sampling methods would not be effective. Weightings of the MSE loss take the form: where w(φ) is the per-voxel weighting function which depends only on the true damage values φ. Several possible weighting functions are depicted in Fig. 7: uniform weighting, logistic weighting, inverse histogram (IH) weighting based on the ensemble, and inverse histogram per realization (IH pr ) weighting. The optimization procedure responds to the weighting by finding a solution with a different distribution of the accuracy per damage value shown in Fig. 7. By weighting the high-frequency true values with smaller weights, the regression accuracy can be made more uniform across the range of damage values. While using a weighting inversely proportional to the damage value frequency, either per realization or for the whole ensemble, seems to create the most even distribution of prediction accuracy, the logistic weighting provides a similar effect on the loss. A comparable effective treatment to inverse histogram weighting could likely be achieved by optimal choice of the logistic functions parameters. To avoid this parameter tuning, we focus on the inverse histogram weightings which are generated directly from the training data. Note that a drawback of reweighting is that by downweighting the contribution of smaller damage values, the loss function penalizes poor small damage value predictions less which leads to more false positives. Although predicting the low background damage is not the primary goal, this behavior is a consideration in choosing the optimal weight function w. Related discussion will be given in Sec. 4.3.

A heteroencoder convolutional neural network
In Sec.2.2 we described a reduced-dimensional representation of the explicit porosity process based on a Karhunen-Loève expansion (KLE). It is optimal in the sense that given a fixed number of modes, no other linear decomposition can capture as much energy of the random porosity process as the KL expansion. In this setting, we decided to capture 99.9% of the energy which requires ≈12,000 KLE modes. In other terms, the random porosity process can be optimally captured using ≈12,000 parameters, namely coefficients of orthogonal modes of a linear decomposition. In the context of our modeling challenge, KLE, being a linear approximation, can be used to project the input signal onto ≈12,000 vectors. This is a form of encoding of the high-dimensional porosity realizations which serve as inputs to regression problem. Our problem could be modeled with a linear encoder-decoder where W 1 is an encoding linear transform/matrix, W 2 is a decoding one and b is a bias correction. The KL process already provides the optimal version of W 1 , being the collection of all "linear" KLE modes. However, this is not an auto-encoder since the output (damage field, φ) is not an approximation of the input (porosity field, ϕ). This linear model, which we call a heteroencoder [5], attempts to reduce the input porosity to a smaller latent space then expands back to the high-dimensional output, namely the damage field. The reduction is formally achievable given the dimensionality of the KLE approximation of the porosity process used to generate the input ϕ.
Starting with the KLE as a baseline for dimensionality reduction from the naïve dimensionality of the porosity on the grid, the dimensionality of the latent space can be further reduced using a non-linear latent variable model. The idea is similar to the use of non-linear principal component analysis (PCA) methods, and uses neural network training procedures to generate nonlinear features to achieve further dimensionality reduction of the latent variable space [32]. This can be enacted using networks with nonlinear activation functions.
In general, an output field often pertains to a small number of functionals of the input process (or latent variables thereof) [33], which justifies further reduction of the number of requisite latent variables. In our context, the output field to be modeled relates to regions of high damage near failed elements which should depend on a subset of, rather than all, the latent variables. The convolutional architecture and the spatial smoothing of the input porosity process (achieved in a data pre-processing step) helps further reduce the requisite number of latent variables by diminishing high-frequency content in the input. Lastly, we note that the use of softmax to define the voxel-wise vector of damage field values,φ I , confers an attention network form to the overarching framework [34].
Therefore, starting with a "linear" KLE model and associated linear encoder-decoder with a latent variable space of dimensionality ≈12,000, we propose a non-linear heteroencoder with a total of ≈7,000 latent variables (refer to Fig. 8). This heteroencoder CNN is parametrically parsimonious compared to fully connected feed-forward networks and CNNs that have all layers operate at full resolution, hence the ≈110,000 samples we have for training and validation are sufficient.
Training was carried out with a 7:1:2 train/validation/test split of the available data. An Adam optimizer [35] was employed with a learning rate of 10 −3 over 300 epochs. A batch size of 128 was used for the deterministic networks, while a larger batch size of 256 was used for the probabilistic networks (Sec. 5) to reduce the variance of the gradient estimations. The weights were randomly initialized for deterministic networks, while a warm-start strategy is adopted for mean parameters of the Bayesian networks by initializing them from a corresponding trained, deterministic model.

Evaluation of the convolutional neural network prediction performance
In the previous sections, transformations of the damage field and modifications of the MSE loss were motivated as strategies for improving the ability of the CNN to predict the relatively sparse regions displaying high damage. We investigate these choices in a systematic manner to compare their effect on CNN performance and to reveal a more global picture of how the CNN model depends on these hyperparameters. (20,20,

Post-processing metrics
While the MSE provides a standard measure of predictive accuracy for regression problems that is easily understood and formulated as a differentiable loss function, it does not completely capture some of the engineering aspects of performance. Given the consequences of failure, it is more important to capture the locations of clusters of large damage values than it is to match the exact damage value predictions across the entire field. Also predictions of failure in close proximity to the true failure locations are valuable in that they indicate clusters of porosity that are most susceptible to failure. This motivates the consideration of alternative post-prediction performance metrics to evaluate the quality of the model after training.
Overlap metrics There are a number of metrics which provide a standard measure of similarity between two sets which appear frequently in binary classification tasks such as image segmentation and information retrieval. Three such metrics, relevant to our work, are precision, recall, and overlap which are defined as: Overlap = T P T P + F P + F N where TP, FP, and FN denote true positives, false positives, and false negatives, respectively. True positives are defined as the voxels which are predicted to fail by the CNN and also were actual failed voxels in the physical simulations that provided the data. Set-wise: TP = Fail true ∩ Fail pred . These measure different aspects of similarity between true and predicted sets and are illustrated in Fig. 9. Note that in the present context these sets (and the associated spatial clusters) depend on a selected damage threshold value.

True negatives
True positives False negatives False positives Intuitively, precision reflects the percentage of positive predictions which are accurate, while recall measures the percentage of true values which are captured by the prediction. The F-score combines the precision and recall metrics to create an overall measure of accuracy from the harmonic mean of these two quantities. On the other hand, overlap measures the extent to which the true and predicted regions coincide. The overlap metric is also known as the Jaccard index, originally defined to measure the similarity of sets or, equivalently, binary vectors. The Sorenson-Dice coefficient is a similar metric related to the Jaccard index by a monotonic transformation such that the two metrics induce the same ordering on sets and can be considered equivalent [36]. For our goal of predicting all the likely failure locations recall is the most important metric, since moderate over-prediction of failure provides a useful and conservative estimate of reliability.
Cluster overlap metric The overlap metrics based on voxel-wise comparison serve as strict measurements to evaluate the performance of the proposed CNN. Considering the continuous nature of the damage field, the damage locations from CNN predictions often appear as spatially contiguous clusters for sufficiently high threshold values. To have an evaluation of CNN performance that accommodates a sense of cluster proximity, we also calculate a relaxed overlap metric based on damage clusters instead of true and predicted damage values in the same voxel. In this metric, we consider a predicted damage cluster as a true positive if it overlaps with a true damage cluster by at least one voxel. Conversely, it is considered to be a false positive if it does not overlap with any true damage clusters. Likewise, a true damage cluster is considered to be a false negative if it does not overlap with any predicted damage cluster.

Cross-validation of model variants
A 5-fold cross validation study of various choices of data transformation and loss weighting functions was carried out and evaluated according the the overlap metrics discussed in the previous section. Variants of the CNN architecture shown in Fig. 8 were trained for each combination of the following choices:

Data transformation
Notation Loss re-weightings Notation None id None id Softmax σ Inverse histogram IH Gaussian filter g Inverse histogram per realization IH pr Softmax • Gaussian filter t The precision and recall metrics for CNNs based on all 12 combinations are plotted as a function of damage value threshold according to the standard overlap metrics in Fig. 10a and according to the cluster overlap metric in Fig. 10b. Note we only show results for the cluster overlap metric for the range of interest, i.e. where the threshold is high enough to discern distinct damage clusters and regions of likely failure, due to the computational expense of computing the cluster overlap metric.
Observe that recall and precision both improve as the data becomes increasingly smooth from the data transformation. In particular, the least smooth data is produced by softmax (σ) which tends to amplify gaps in the underlying field. For the unweighted MSE, this displays the worst performance across the entire range of damage values. The softmax and Gaussian filter transformation (t) and Gaussian filter (g) provide relatively smooth data the highest precision and recall values. Both transformations are generally comparable and seem to marginally outperform one another in certain damage value regions. This suggests  that the smoothness induced by convolution is the most significant factor in improving CNN predictions on the damage field. As expected, the loss weighting schemes combat class imbalance by emphasizing lowfrequency, high-damage regions and result in an improvement to the recall and precision metrics. Recall improves across the entire damage range while precision improves only at the upper end. This is likely due to precision decreasing with false positives which increase with loss weightings as discussed in Sec. 3.2.
Note that the metrics display a general pattern of evolution. Initially, there is a constant region at the maximal metric value of one. This is followed by a region of rapid decrease which then slows and, possibly, reverses direction, sometimes forming a cusp in the graph. After this, there is a local maximum and the function then monotonically decreases.
To explain these phenomena, we consider Fig. 11. For each data transformation, Fig. 11 displays the damage value frequency distribution for true and predicted damage fields as well as the proportion of damage values above a given threshold for the unweighted MSE loss. The initial plateaus are due to the prediction problem being degenerate for low thresholds, i.e. essentially all voxels fail. The cusps occur at the frequency peaks in the true damage fields suggesting that the trend reversal is because of the data frequency dependent performance of the CNN. This is further supported by the observation that the loss reweightings move the local maxima toward the higher end of the damage value region in Fig. 10.
Note that the predicted frequency distribution does not exactly capture the shape of the true frequency distribution. This discrepancy is amplified at the true damage field frequency peaks which are overemphasized by the predictions as an effect of class-imbalance. This is reflected in the second row of Fig. 11 which shows the volume of superlevel sets and reflects how quickly the damage fields transition to sparse clusters as the damage threshold increases. Here, the superlevel sets reflect the collection of high damage clusters present after a thresholding operation transforms the continuous damage field to a binary field. The predicted damage fields make a delayed, but faster transition to sparse clusters which accounts for the initial rapid decrease in the overlap metrics since the two fields will temporarily display significantly different volumes. This shows that the peaks of the true damage fields are smoother than those of the predicted fields resulting in different structure to the level set evolution. Similar behavior is seen with each data transformation except the identity transformation where the transition of both the true and predicted level sets nearly coincide.

Discussion of porosity sensitivity and predicted false positives
The CNN predicted results of realizations corresponding to those in Fig. 4(a,g,h,l), are presented in Fig. 12, where the pores, the true damage clusters and CNN predicted damage clusters for a threshold value of 0.8, and the failure locations from the sensitivity analysis are shown. This CNN is trained with the softmax and Gaussian filter data transformation using an inverse histogram per realization weighted loss. Fig. 12 shows that all the true damage clusters in these four realizations overlap with CNN predicted damage clusters. This is consistent with the cluster recall curve in Fig. 10, which has a value above 0.9 at a threshold of 0.8 for the selected data transformation and loss weighting. Even though the cluster precision value in Fig. 10 for a threshold of 0.8 is around 0.3, we observe from Fig. 12 that some false positives indicate potential failure locations that were captured by the sensitivity analysis. The CNN model appears to not only learn the complex nonlinear mapping between the pore distribution and damage initiation location (high recall score), but also learn subtle relationships between the two fields that were only revealed by the sensitivity analysis with expensive physics-based direct numerical simulations.

A Bayesian convolutional model of localized failure
We investigated the addition of uncertainty quantification (UQ) to the hetereoencoder CNN model described in Sec. 3 using variational inference (VI) to approximate the Bayesian posterior distribution of the weights given training data. This addition provides: (a) an assessment of the sufficiency of the training data volume as reflected in the distribution of uncertainty, and (b) a mitigation of possible overfitting by the regularization provided by Bayesian prior distributions.

Uncertainty quantification with Bayesian inference
In general, model predictions can be endowed with uncertainty estimates by allowing the model parameters to be random variables whose distributions are calibrated to the training data. In the present context, we let w denote the parameters, both the weights and biases, of the neural network transformation, written as NN w (ϕ), and let D = {ϕ i , φ i } denote the training data where i ∈ {1, N s } and N s is the number of data points. The data is assumed to be generated by the neural network transformation and additive Gaussian noise such Figure 12: CNN prediction of four selected realizations corresponding to Fig. 4(a,g,h,l) for a threshold value of 0.8. Red: true damage, blue: CNN prediction, yellow: failure locations from the sensitivity analysis, gray wireframe: pore location.
that φ = NN w (ϕ) + = φ pred + where ∼ N (0, σ 2 ) captures model discrepancy. Note that we make a distinction between the exact and noisy outputs φ pred and φ, respectively. The noise accounts for the approximate nature of the model NN w and defines the likelihood distribution p(D | w) over the data given the weights. The posterior probability of the parameters given the training data is provided by Bayes' rule: where p(w) is the prior and p(D) is the evidence. Uncertainty in the parameters as captured by p(w | D) can then be pushed forward through the model to induce uncertainty in the outputs resulting in the pushed forward posterior distribution (PFP) [37]. We denote this as p(φ pred | ϕ, D), which defines the distribution over outputs φ pred given a particular input ϕ. The PFP distribution is not analytically tractable so we employ a Monte Carlo scheme to estimate its first two moments based on sampling the weights w ∼ p(w | D) according to the posterior.

Approximating posteriors with variational inference
The model evidence integral in the denominator of Eq. (16) is typically intractable to evaluate, which necessitates methods to approximate the posterior (16) and to sample from the PFP. Briefly, there are a variety of techniques to enable Bayesian calibration/inference with a NN, e.g. variational inference (VI) [38], stochastic gradient descent (SGD)-based [39], and Dropout [40] schemes. Variational inference minimizes a measure of discrepancy, such as Kullback-Liebler divergence between a tractable parametric family F of probability distributions and the actual posterior p(w|D) in order to select the optimal member q θ of F. Here θ is a vector parameterizing the family in the sense that F = {q θ | θ ∈ Θ ⊆ R n }. This gives us a closed form approximation of the true posterior whose fidelity is determined by the chosen family of distributions F. By sampling from this approximate posterior, we can form Monte Carlo estimates of the PFP distribution.
In this work, we use stochastic variational inference (SVI) which formulates a differentiable cost function equivalent to (17) and then carries out a form of gradient descent optimization. The evidence lower-bound (ELBO) is a standard cost function in VI. The ELBO L(θ) differs from the KL-divergence by a constant, the log of the evidence log p(D): A detailed derivation of Eq. (18) can be found in Ref. [38]. From Eq. (18), it is clear that we can minimize the KL-divergence by maximizing the ELBO or, equivalently, minimizing the negative of the ELBO: Note that the ELBO no longer references the intractable posterior p(w | D) and hence can be approximated in an optimization procedure.
To provide an understanding of the ELBO loss we relate it to the MSE loss used in deterministic CNN training, refer to Sec. 3. We first focus on the (second) data-dependent term of Eq. (19) and assume the likelihood p(D | w) takes the form of a mean field Gaussian with means given by the model predictions NN w (ϕ s ) and homoscedastic noise σ. Note that we treat σ as an additional optimization parameter that is calibrated along with θ in the training of our Bayesian CNN models. We can express the likelihood component as an expectation of the deterministic MSE with respect to the approximate posterior q θ (w). Note N v is the dimension of the damage field and c is a constant with respect to parameters θ and noise σ 2 . The (first) KL-divergence term in Eq. (19) involves the prior and exhibits no data-dependence. It provides a regularization that drives the posterior toward the prior. This suggests a connection between the ELBO cost and regularized least squares, similar to maximum likelihood estimation. We can illustrate the connection by considering the case where the model is defined by NN W (x) = Wx, i.e. , the model is simply a linear transformation by the weight matrix W ∈ R n×n , as in PCA or a simple network model. In this case, we can derive a closed form for the ELBO loss function. Let q θ (W) = N (W | µ q , Σ q ) and p(W) = N (W | µ p , Σ p ) be matrix valued distributions denoting the approximate posterior and prior, respectively, where µ q ∈ R n×n is the matrix of means and Σ q ∈ R n 2 ×n 2 is the corresponding covariance. Together, these constitute the variational parameters θ = (µ q , Σ q ) ∈ R n 2 +n 4 . Parameters µ p ∈ R n×n and Σ p ∈ R n 2 ×n 2 are defined similarly for the prior p(W) but are constants in the optimization process. Given these assumptions, the ELBO has the form: and Y = [φ 1 · · · φ Ns ] are data matrices consisting of the N s samples of porosity inputs and damage outputs, respectively. The matrix V is defined by V ij = Σ kikj given the definition E[W ij W kl ] = Σ ijkl + µ ij µ kl with Σ ijkl being the components of Σ q . Again, c is a constant with respect to the variational parameters θ. Note the details of this derivation, as well as a 1-dimensional example, are given in App. A.
Observe that with respect to the means, the problem has the form of least squares 1 2σ 2 tr{(Y−µ q X) T (Y−µ q X)} with a quadratic regularization term (µ p −µ q ) T Σ −1 p (µ p −µ q ) arising from the KL-divergence D KL (q θ (w | D) p(w)) portion of the ELBO. The noisedependent term log det(Σ p Σ −1 q ) + tr(Σ −1 p Σ q ) also arises from the KL-divergence and is a convex function attaining a unique minimum at Σ q = Σ p . This term, along with the quadratic regularization, tend to drive the posterior to the prior. The loss component 1 2σ 2 tr{VXX T } arises from the likelihood portion of the ELBO and is characterized by the input data dependent matrix XX T . This term is uniquely minimized by Σ q = 0 (or when V is orthogonal to XX T , when XX T is not full rank). It has the effect of limiting the magnitude of the variance, which minimizes the propagation of noise through the model W and improves the accuracy of predictions.
In the general, nonlinear case when the ELBO lacks a closed form expression, we can gain some intuition from (20) when the surrogate posterior is a Gaussian distribution q θ = N (w | µ q , Σ q ) by letting MSE(w) = 1 Ns Ns s=1 (φ s − NN w (ϕ s )) 2 so that (20) becomes where Hence the likelihood portion of the loss is a function of θ = (µ q , Σ q ) and takes the form of a convolution in µ q with a Gaussian kernel of covariance Σ q . This has the effect of smoothing the MSE to differing extents over different regions of parameter space. We expect smoothing to approximately preserve local minima so that they are carried over from the MSE to the ELBO.

Training of the Bayesian neural network
To train the model, we use SVI which utilizes gradient descent optimization of the ELBO cost function, Eq. (19). At each iteration, a Monte Carlo scheme is used to estimate the gradient of the ELBO. One way to obtain a gradient estimator is to reparametrize random variables to allow the gradient of the ELBO to be expressed as an expectation with respect to a probability density which can be estimated with standard Monte Carlo [38]. Due to the scale of neural networks, it is only feasible to use a small number of Monte Carlo samples. This estimator is then subject to large variance, so we use Flipout [41] to reduce the variance without having to use additional random samples. The Flipout method decorrelates gradient samples within a batch by generating pseudo-independent perturbations of the samples via random matrices.

Challenges with variational inference
A common challenge in optimization is non-convexity of the loss function given a model that is non-linear in its parameters. This leads to multiple local minima such that local optimizers may get stuck at a non-optimal solution, depending on initial conditions. It then becomes necessary to employ strategies to lead the optimizer to a better local minimum such as initialization procedures or transformations of the loss function. The form of the ELBO (19) indicates that, while the KL-divergence term is convex for Gaussian distributions, the likelihood term in (20) will inherit non-convexity from the MSE loss. Indeed, Refs. [38,42,43] report issues with SVI converging to poor local minima and employ strategies such as annealing to mitigate this. In preliminary studies of the BCNN, we observed SVI converging to a model which predicted a uniform field at approximately the highest-frequency damage value independent of the input. The uniform field solution was found to also be present in the deterministic CNN, which confirmed that the likelihood term does appear to inherit local optima from the MSE as suggested by (22). By re-weighting the BCNN loss or using an informed initialization scheme, the uniform field local minimum could be avoided suggesting that the existence of this solution is a result of class imbalance. Namely, the uniform field attains a low MSE value because the damage fields are approximately uniform in that the high-frequency damage values far outnumber other damage values. Viewing the likelihood term in the form (22) suggests that local optima in the likelihood will have mean values approximately equal to the weights representing local optima in the MSE. Hence, to avoid poor local minima, we initialize the weight means µ q of the BCNN with the weights of the trained CNN in a "warm-start" scheme.

Evaluation of the Bayesian convolutional neural network model predictions
The BCNN provides a mean prediction, which can be thought of as a point estimate of the random weights and defines a deterministic mean transformation of ϕ to φ. The variance of the weights induces a random distribution of outputs φ around the mean providing uncertainty quantification on top of deterministic predictions. Both the mean and the variance of the pushed forward posterior distribution over outputs can be estimated through Monte Carlo sampling of the pushed forward posterior distribution. Each sample is generated by first sampling the weights according to their approximate Bayesian posterior (16) obtained through variational inference. A forward pass with this random sample of the weights is used to generate a sample of the posterior distribution. The sample mean and variance are then used to estimate the actual actual mean and variance.

Mean predictions
To analyze the mean predictions, the BCNN was compared to the cross-validation study from Sec. 4 and evaluated according to the recall and precision metrics. Two different likelihood loss weighting schemes were used analogous to the unweighted and the inverse histogram weighted versions of the MSE. Given the results of Sec. 4.2, only the softmax and Gaussian filter transformation were considered. The results are displayed in Fig. 13 and show that for the unweighted loss, the mean predictions of the BCNN slightly outperform all of the deterministic CNNs. This indicates that the inherent regularization provided by Bayesian neural networks is improving the quality of the solution found by the optimizer. The BCNN trained using an inverse histogram weighting scheme for the likelihood term of the ELBO performs significantly worse than the best of the inverse histogram weighted CNNs. For the unweighted loss, the calibrated noise standard deviation was σ = 0.065, while for the weighted loss it was significantly higher. This suggests that weighting schemes may have a different effect in the ELBO loss despite the connection to MSE.

Uncertainty quantification
In the previous section, we established that the mean predictions of the unweighted BCNN perform similarly to the CNNs and even display improvement in the precision and recall metrics compared to the unweighted CNN. In addition to the mean prediction, the BCNN also provides useful information about uncertainty in the network predictions given the amount of training data. In this section, the uncertainty in the posterior predictive distribution is analyzed and used to define a method to classify damage clusters according to their likelihood of failure.
Uncertainty in the posterior predictive distribution The outputs of the BCNN are stochastic samples from a posterior predictive distribution which reflect the variance or uncertainty around the network's mean prediction. A well-known phenomenon in Bayesian problems is the convergence of the posterior to a Dirac delta, a distribution with no uncertainty, around a particular parameter value with increasing amounts of data. Similarly, in simple regression problems with Bayesian networks, it can be seen that uncertainty in the outputs correlates with data sparsity in the domain, where sparse regions lead to higher uncertainty [44]. From these observations, we expect the uncertainty in the damage field BCNN to reflect data sparsity. This suggests viewing the uncertainty as a function of damage value, similar to the loss distributions in the previous section. This quantity is plotted in Fig. 14 where the variance across field voxels is binned according to the damage value. Also displayed is the frequency distribution of the predicted damage field. Examination of Fig. 14 indicates that the class imbalance aspect of the problem is affecting the uncertainty. Comparison of the uncertainty distribution with the frequency distribution shows the damage value regions of high frequency exhibit a smaller level of uncertainty and vice versa. Also note that the inverse histogram weighting of the log likelihood loss causes an overall increase in the uncertainty across the range of damage values. In the limit of w(φ) weighting any subset of data with zero weight, the size of our training set decreases. Hence, down-weighting the high frequency regions of the damage spectrum is effectively reducing the amount of data we are using to train which results in a general increase in the amount of uncertainty. This is consistent with what is expected in Bayesian inference. Also observe that the uncertainty for the inverse histogram BCNN decreases toward larger damage values reflecting the emphasis that the inverse histogram weighting places on the high damage region.
Ranking damage clusters with UQ Now that we have a characterization of the uncertainty the BCNN predictions is capturing, we use both the mean and variance of the posterior predictive outputs to analyze the damage clusters. A natural approach is to threshold the damage field to obtain multiple high damage value clusters (superlevel sets) and then classify them based on an appropriate metric representing the likelihood of failure. To this end, we would like a metric which depends on the mean and variance of each cluster's Monte Carlo predictions. A metric motivated by reliability analysis which satisfies this condition is the probability mass above a given threshold. This allows ranking of the clusters according to their average damage value in a manner sensitive to the distribution around this average. A cluster whose distribution shows a significant likelihood of producing a stochastic realization below its average value will be ranked below a cluster with the same average with lower variance. To illustrate this concept the empirical CDFs of three clusters from the same damage field samples are shown in Fig. 15. The empirical distributions coincide at the lower and upper limits of the damage threshold but display a consistent ranking across an intermediate region of the damage threshold where the CDFs for all are less than one and greater than zero. This suggests that the probability mass metric is fairly invariant with respect to choice of threshold and provides a robust method for ranking clusters. This ranking scheme was carried out on realizations (a), (c), (g), and (j) from the sensitivity analysis in Fig. 4 with results displayed in Fig. 16. Realizations (a), (c) display low sensitive and have only a single, small cluster of perturbed failures. Observe that the probability mass ranking successfully determines the damage clusters containing these perturbed failures as the most likely failure regions. Similarly, realization (j) has two perturbed failure regions and these clusters are ranked more highly than the clusters which do not coincide with perturbed failure locations. Realization (g) is the most sensitive and each damage cluster coincides with a cluster of perturbed failures. The probability mass ranking reflects this by assigning similar rankings to all four clusters indicating that the failure likelihoods are more difficult to distinguish in this case.
In this ranking scheme, a connected components algorithm [45] designed for 3-dimensional fields was used. It employs several strategies which provide optimal time complexity for the algorithm. For more general clustering, we can represent voxel connectivity as a graph and use graph-based connected component algorithms such as depth-first-search or union-find. Figure 16: Classification of damage clusters in realizations (a), (c), (g), and (j) using the probability mass above a threshold. Red indicates regions of higher mass while yellow represents the perturbed failure locations.

Discussion and conclusion
We demonstrated that a hetereoencoder comprised of an encoder of the input and a decoder that provides predictions of an output with different character is an effective representation of the porosity-to-damage-at-failure map for a porous ductile AM specimen. The architecture proved useful in conjunction with data and loss transformations designed to alleviate the extreme class imbalance and spatial irregularity of the data, and yet preserve the ac-curacy of the predictions of failure locations. We were also able to illustrate the utility of adding variational inference to this architecture. Not only does it regularize the training/optimization of the network, VI provides indications of whether the data is sufficient to train the network beyond the usual comparison of training and test errors and more elaborate cross-validation. The UQ provided by the BCNN architecture also proved useful in the pragmatic task of ranking the location most likely to fail. The methods we developed have direct application to reliability assessment based on non-destructive CT scans and structural optimization given the inherent porosity of the AM process.
As a final observation to complement the examination of the uncertainty in the pushed forward posterior distribution of the BCNN in Sec.6.2, we studied the posterior distribution of the network parameters. In particular, we looked at the overall uncertainty in the distribution of the weights as a function of each layer of the BCNN. For each layer, we form a mean coefficient of variation (CV) from i σ i / i µ i , the mean of standard deviations of the weights σ i over the mean of the means of the weights µ i . Fig. 17 shows that on average, the uncertainty in both the convolutional kernel and bias parameters increases with depth in the network. We offer two explanations for this phenomenon. The first is that uncertainty in the initial layers of the network is propagated through a composition of nonlinear functions which may display high sensitivity to inputs. Large uncertainty in the initial layers would lead to highly unpredictable outputs and thus, poor predictions. Hence the uncertainty distribution may reflect effects of noise propagation. The second explanation is based on the similarity of the BCNN architecture to an encoder-decoder which possesses a low-dimensional intermediate dimensional layer, which, in our case, is the result of the third convolutional layer of the stack of six. Decoding transformations (layers 4-6) can be seen as allowing for larger variance as they are essentially adding information in reconstructing the output at the full resolution.
This observation has ramifications to network designs in that it may be sufficient to limit stochastic behavior to the final layer or layers of a neural network. This would allow for UQ at significantly reduced computational cost. It also motivates the study of how the uncertainty distribution across a network may correlate to the occurrence of different types of transformations such as the encoding and decoding portion of an autoencoder or heteroencoder.
The issue of choosing a threshold to predict regions likely to fail remains. This suggests post-training tuning based on validation and other means of subverting or optimizing the threshold for accuracy. A number of methods for segmenting images exist [46,47,48,49] and could be adapted to the failure prediction problem to improve the proposed methodology.

A Derivation of the ELBO for a linear model
Herein, we provide supplementary details for derivation of the ELBO cost function Eq. (21) associated with the mean field approximation of the parameter posterior for the linear transformation case described in Sec. 5.2. We subsequently provide the ELBO for the single parameter linear model case as a form of validation.

A.1 Derivation for the N-dimensional case
Recall that the cost function for variational inference takes the form of the negative ELBO as in Eq. (19). The KL-divergence term has a well-known analytical formula so we focus on deriving an expression for the likelihood term. Starting with the mean-field normal form of the likelihood, we expand this term to where Y = [φ 1 · · · φ Ns ], X = [ϕ 1 · · · ϕ Ns ] are output and input data matrices. In this case, the input-output transformation defined by the neural network model is linear such that NN W (ϕ) = Wϕ with W ∈ R n×n being parameters which we assume are normally distributed W ∼ N (µ q , Σ q ). Note that we have assumed, without loss of generality, that the input and output vectors are of the same dimensionality n = N v . The last term c = NsNv 2 log(2π) is a constant with respect to the variational parameters θ and σ 2 , which are both optimized during training. Focusing on the matrix expression, we obtain To simplify notation, we will denote the covariance matrix by Σ q ≡ Σ, where we have dropped the q in anticipation of subscript indices. Typically, the covariance is thought of as a matrix Σ ∈ R n 2 ×n 2 over the random vector vec(µ q ) ∈ R n 2 where vec(·) is the operation which maps an n × m matrix into a column vector in R nm . Here, we adopt the alternative perspective of defining the covariance as a rank-4 tensor so that we have E q θ (W) [W ⊗ W] = Σ + µ q ⊗ µ q or, equivalently, in component notation, E q θ (W) [W ij W jk ] = Σ ijkl + µ ij µ kl . To obtain E q θ (W) W T W , we can contract this expression, with implied Einstein summation convention, to get E q θ (W) W T W ij = Σ kikj + µ ki µ kj = V ij + µ ki µ kj where we have defined V ij = Σ kikj . Returning to matrix notation, we have 1 2σ 2 tr{E q θ (W) W T W XX T } = 1 2σ 2 tr{(V + µ T q µ q )XX T } = 1 2σ 2 tr{X T µ T q µ q X} + 1 2σ 2 tr{VXX T } so that the likelihood term can finally be written in terms of known matrix quantities Hence, our full negative ELBO cost function is as in Eq. (21).

A.2 Reduction to the 1-dimensional case
To gain intuition into how the various terms contribute to the loss function, consider the case where W = w ∈ R 1×1 , i.e. , our matrix is just scaling a univariate input by a constant to produce a scalar output. In this single parameter case, the negative ELBO cost function becomes where q θ (w) = N (w | µ q , σ 2 q ), p(w) = N (w | µ p , σ 2 p ) are the surrogate posterior and the prior on w, respectively.
To check the correctness of this formula, we show that we can express the negative ELBO in the 1-dimensional case as L(µ l , σ l , µ p , σ p , µ q , σ q ), where µ l , σ l are the mean and standard deviation of the likelihood probability density function of parameter w. Bayesian inference does not distinguish between the prior and likelihood in the sense that p(w | D) = p(D | w)p(w) p(D) = p(w)p(D | w) p(D) . Hence, as VI will find the exact Gaussian posterior in this case, we should expect that the optimal solution µ * q , σ * q is invariant under exchanging the prior and likelihood parameters (mean and standard deviation), i.e. , µ * q , σ * q = arg min µq,σq L(µ l , σ l , µ p , σ p , µ q , σ q ) = arg min µq,σq L(µ p , σ p , µ l , σ l , µ q , σ q ) Letting y = [φ 1 · · · φ Ns ], x = [ϕ 1 · · · ϕ Ns ], recall that Bayesian linear regression with a uniform prior is defined by the log likelihood χ 2 = − 1 2σ 2 (y − wx) T (y − wx). We can extract the likelihood mean and variance from the first and second order derivatives of the log likelihood [50]: so that µ l = y T x x T x and σ 2 l = σ 2 x T x . Based on these expressions, the first term in the negative ELBO can be written as in which the first term is a complete square of the difference between the likelihood and variational posterior means. We also note that Hence the univariate ELBO can be re-expressed as L(µ l , σ l , µ p , σ p , µ q , σ q ) = 1 2σ 2 where c = Ns 2 log σ 2 + 1 2σ 2 l y T y x T x − µ 2 l + 1 2 log σ 2 p is a constant with respect to the variational parameters µ q , σ q . Therefore, the optimal solution is which is preserved under the interchange of the prior and likelihood parameters, as expected.