Probabilistic Reduced-Order Modeling for Stochastic Partial Differential Equations

TL;DR

This paper introduces a Bayesian, data-driven coarse-grained model for stochastic PDEs that learns microstructural features to efficiently predict fine-scale responses with quantified uncertainty.

Contribution

It develops a physically inspired probabilistic model that captures microstructural features for efficient and uncertain-aware predictions of stochastic PDEs.

Findings

Learned low-dimensional microstructural features predict fine-scale responses.

Model provides probabilistic predictions with quantified epistemic uncertainty.

Sparsity priors enhance feature selection and model efficiency.

Abstract

We discuss a Bayesian formulation to coarse-graining (CG) of PDEs where the coefficients (e.g. material parameters) exhibit random, fine scale variability. The direct solution to such problems requires grids that are small enough to resolve this fine scale variability which unavoidably requires the repeated solution of very large systems of algebraic equations. We establish a physically inspired, data-driven coarse-grained model which learns a low- dimensional set of microstructural features that are predictive of the fine-grained model (FG) response. Once learned, those features provide a sharp distribution over the coarse scale effec- tive coefficients of the PDE that are most suitable for prediction of the fine scale model output. This ultimately allows to replace the computationally expensive FG by a generative proba- bilistic model based on evaluating the much cheaper CG several times. Sparsity enforcing pri- ors further increase predictive efficiency and reveal microstructural features that are important in predicting the FG response. Moreover, the model yields probabilistic rather than single-point predictions, which enables the quantification of the unavoidable epistemic uncertainty that is present due to the information loss that occurs during the coarse-graining process.

Figures

Equations

$$\mathcal{L}\left(\mathbf{x},\lambda (\mathbf{x},\xi ))u(\mathbf{x},\lambda (\mathbf{x},\xi )\right)=f(\mathbf{x}),\text{+B.C.}$$

$${\mathbf{r}}_f({\mathbf{U}}_f;{\mathbf{\lambda }}_f)=0,$$

$$p({\mathbf{\lambda }}_f)=\int \delta ({\mathbf{\lambda }}_f-{\mathbf{\lambda }}_f(\xi ))p(\xi )d\xi$$

$$p({\mathbf{U}}_f)=\int p({\mathbf{U}}_f|{\mathbf{\lambda }}_f)p({\mathbf{\lambda }}_f)d{\mathbf{\lambda }}_f,$$

$$\begin{array}{rl}\bar{p}({\mathbf{U}}_f|{\mathbf{\lambda }}_f)& =\int \bar{p}({\mathbf{U}}_f,{\mathbf{\lambda }}_c|{\mathbf{\lambda }}_f)d{\mathbf{\lambda }}_c\\ & =\int \underset{\text{decoder}}{\underbrace{\bar{p}({\mathbf{U}}_f|{\mathbf{\lambda }}_c)}}\underset{\text{encoder}}{\underbrace{\bar{p}({\mathbf{\lambda }}_c|{\mathbf{\lambda }}_f)}}d{\mathbf{\lambda }}_c.\end{array}$$

$${\mathbf{r}}_c({\mathbf{U}}_c;{\mathbf{\lambda }}_c)=0,$$

$$\bar{p}({\mathbf{U}}_f|{\mathbf{\lambda }}_f,\mathbf{\theta })=\int p_{cf}({\mathbf{U}}_f|{\mathbf{U}}_c({\mathbf{\lambda }}_c),{\mathbf{\theta }}_{cf})p_c({\mathbf{\lambda }}_c|{\mathbf{\lambda }}_f,{\mathbf{\theta }}_c)d{\mathbf{\lambda }}_c,$$

$$\begin{array}{rl}\text{KL}\left(p({\mathbf{U}}_f,{\mathbf{\lambda }}_f)||\bar{p}({\mathbf{U}}_f,{\mathbf{\lambda }}_f)\right)& =\text{KL}\left(p({\mathbf{U}}_f,{\mathbf{\lambda }}_f)||\bar{p}({\mathbf{U}}_f|{\mathbf{\lambda }}_f,\mathbf{\theta })p({\mathbf{\lambda }}_f)\right)\\ & =\int p({\mathbf{U}}_f,{\mathbf{\lambda }}_f)\log \left(\frac{p({\mathbf{U}}_f,{\mathbf{\lambda }}_f)}{\bar{p}({\mathbf{U}}_f|{\mathbf{\lambda }}_f,\mathbf{\theta })p({\mathbf{\lambda }}_f)}\right)d{\mathbf{U}}_{f}d{\mathbf{\lambda }}_f\\ & =-\int p({\mathbf{U}}_f,{\mathbf{\lambda }}_f)\log \bar{p}({\mathbf{U}}_f|{\mathbf{\lambda }}_f,\mathbf{\theta })d{\mathbf{U}}_{f}d{\mathbf{\lambda }}_f+H(p({\mathbf{U}}_f,{\mathbf{\lambda }}_f))\\ & \approx -\frac{1}{N}\sum _{i=1}^N\log \bar{p}({\mathbf{U}}_f^{(i)}|{\mathbf{\lambda }}_f^{(i)},\mathbf{\theta })+H(p({\mathbf{U}}_f,{\mathbf{\lambda }}_f)),\end{array}$$

$${\mathbf{\lambda }}_f^{(i)}\sim p({\mathbf{\lambda }}_f),{\mathbf{U}}_f^{(i)}={\mathbf{U}}_f({\mathbf{\lambda }}_f^{(i)}),$$

$$p({\mathbf{\theta }}_c,{\mathbf{\theta }}_{cf}|{\mathcal{D}}_N)\propto e^{L({\mathcal{D}}_N|{\mathbf{\theta }}_c,{\mathbf{\theta }}_{cf})}p({\mathbf{\theta }}_c,{\mathbf{\theta }}_{cf}).$$

$$\begin{array}{rl}[{\mathbf{\theta }}_c^{\ast },{\mathbf{\theta }}_{cf}^{\ast }]& =\arg \underset{{\mathbf{\theta }}_c,{\mathbf{\theta }}_{cf}}{\max }{e}^{L({\mathcal{D}}_N|{\mathbf{\theta }}_c,{\mathbf{\theta }}_{cf})}p({\mathbf{\theta }}_c,{\mathbf{\theta }}_{cf})\\ & =\arg \underset{{\mathbf{\theta }}_c,{\mathbf{\theta }}_{cf}}{\max }\left(L({\mathcal{D}}_N|{\mathbf{\theta }}_c,{\mathbf{\theta }}_{cf})+\log p({\mathbf{\theta }}_c,{\mathbf{\theta }}_{cf})\right)\\ & =\arg \underset{{\mathbf{\theta }}_c,{\mathbf{\theta }}_{cf}}{\max }\left(\sum _{i=1}^N\log \bar{p}({\mathbf{U}}_f^{(i)}|{\mathbf{\lambda }}_f^{(i)},{\mathbf{\theta }}_c,{\mathbf{\theta }}_{cf})+\log p({\mathbf{\theta }}_c,{\mathbf{\theta }}_{cf})\right).\end{array}$$

$$\begin{array}{rl}{L}_i({\mathbf{\theta }}_c,{\mathbf{\theta }}_{cf})& =\log \bar{p}({\mathbf{U}}_f^{(i)}|{\mathbf{\lambda }}_f^{(i)},{\mathbf{\theta }}_c,{\mathbf{\theta }}_{cf})\\ & =\log \int p_{cf}({\mathbf{U}}_f^{(i)}|{\mathbf{U}}_c({\mathbf{\lambda }}_c^{(i)}),{\mathbf{\theta }}_{cf})p_c({\mathbf{\lambda }}_c^{(i)}|{\mathbf{\lambda }}_f^{(i)},{\mathbf{\theta }}_c)d{\mathbf{\lambda }}_c^{(i)},\end{array}$$

$$\begin{array}{rl}{L}_i({\mathbf{\theta }}_c,{\mathbf{\theta }}_{cf})& =\log \bar{p}({\mathbf{U}}_f^{(i)}|{\mathbf{\lambda }}_f^{(i)},{\mathbf{\theta }}_c,{\mathbf{\theta }}_{cf})\\ & =\log \int p_{cf}({\mathbf{U}}_f^{(i)}|{\mathbf{U}}_c({\mathbf{\lambda }}_c^{(i)}),{\mathbf{\theta }}_{cf})p_c({\mathbf{\lambda }}_c^{(i)}|{\mathbf{\lambda }}_f^{(i)},{\mathbf{\theta }}_c)d{\mathbf{\lambda }}_c^{(i)}\\ & =\log \int q_i({\mathbf{\lambda }}_c^{(i)})\frac{p_{cf}({\mathbf{U}}_f^{(i)}|{\mathbf{U}}_c({\mathbf{\lambda }}_c^{(i)}),{\mathbf{\theta }}_{cf})p_c({\mathbf{\lambda }}_c^{(i)}|{\mathbf{\lambda }}_f^{(i)},{\mathbf{\theta }}_c)}{q_i({\mathbf{\lambda }}_c^{(i)})}d{\mathbf{\lambda }}_c^{(i)}\\ & \ge \int q_i({\mathbf{\lambda }}_c^{(i)})\log \left(\frac{p_{cf}({\mathbf{U}}_f^{(i)}|{\mathbf{U}}_c({\mathbf{\lambda }}_c^{(i)}),{\mathbf{\theta }}_{cf})p_c({\mathbf{\lambda }}_c^{(i)}|{\mathbf{\lambda }}_f^{(i)},{\mathbf{\theta }}_c)}{q_i({\mathbf{\lambda }}_c^{(i)})}\right)d{\mathbf{\lambda }}_c^{(i)}\text{(Jensen)}\\ & ={\mathcal{F}}_i(q_i;{\mathbf{\theta }}_c,{\mathbf{\theta }}_{cf}),\end{array}$$

$$\begin{array}{rl}\log p({\mathbf{\theta }}_c,{\mathbf{\theta }}_{cf}|{\mathcal{D}}_N)& =\log L({\mathcal{D}}_N|{\mathbf{\theta }}_c,{\mathbf{\theta }}_{cf})+\log p({\mathbf{\theta }}_c,{\mathbf{\theta }}_{cf})\\ & =\sum _{i=1}^N\log L_i({\mathbf{\theta }}_c,{\mathbf{\theta }}_{cf})+\log p({\mathbf{\theta }}_c,{\mathbf{\theta }}_{cf})\\ & \ge \sum _{i=1}^N{\mathcal{F}}_i(q_i;{\mathbf{\theta }}_c,{\mathbf{\theta }}_{cf})+\log p({\mathbf{\theta }}_c,{\mathbf{\theta }}_{cf})\\ & =\mathcal{F}\left({\left\{q_i\right\}}_{i=1}^N,{\mathbf{\theta }}_c,{\mathbf{\theta }}_{cf}\right)+\log p({\mathbf{\theta }}_c,{\mathbf{\theta }}_{cf}).\end{array}$$

$$q_i({\mathbf{\lambda }}_c^{(i)})\propto p_{cf}({\mathbf{U}}_f^{(i)}|{\mathbf{U}}_c({\mathbf{\lambda }}_c^{(i)}),{\mathbf{\theta }}_{cf})p_c({\mathbf{\lambda }}_c^{(i)}|{\mathbf{\lambda }}_f^{(i)},{\mathbf{\theta }}_c)$$

$${\nabla }_{\mathbf{x}}(-\lambda (\mathbf{x},\xi (\mathbf{x})){\nabla }_{\mathbf{x}}U(\mathbf{x},\xi (\mathbf{x})))=0$$

$$\lambda (\mathbf{x},\xi (\mathbf{x}))=\{\begin{array}{l}{\lambda }_{\text{hi}},\text{if}\xi (\mathbf{x})>c,\\ {\lambda }_{\text{lo}},\text{otherwise}\end{array}$$

$$\text{cov}({\mathbf{x}}_i,{\mathbf{x}}_j)=\exp \left\{-\frac{|{\mathbf{x}}_i-{\mathbf{x}}_j{|}^2}{l^2}\right\}.$$

$${\mathbf{z}}_c=\log {\mathbf{\lambda }}_c.$$

$$z_{c,k}=\sum _{j=1}^{N_{\text{features}}}{\theta }_{c,j}{\chi }_j({\mathbf{\lambda }}_{f,k})+{\sigma }_k{Z}_k,Z_k\sim \mathcal{N}(0,1),$$

$$p_c({\mathbf{z}}_c|{\mathbf{\lambda }}_f,{\mathbf{\theta }}_c,\mathbf{\sigma })=\mathcal{N}({\mathbf{z}}_c|\mathbf{\Phi }({\mathbf{\lambda }}_f){\mathbf{\theta }}_c,\text{diag}({\mathbf{\sigma }}^2)),$$

$$p_{cf}({\mathbf{U}}_f|{\mathbf{U}}_c({\mathbf{z}}_c),{\mathbf{\theta }}_{cf})=\mathcal{N}({\mathbf{U}}_f|\mathbf{W}{\mathbf{U}}_c({\mathbf{z}}_c),\mathbf{S}),$$

$$\begin{array}{rl}\log q_i^{(t+1)}({\mathbf{z}}_c^{(i)})& \propto \frac{1}{2}\sum _{k=1}^{N_{\text{el}}^c}\log \left(({\sigma }_k^{(t)}{)}^{-2}\right)-\frac{1}{2}\sum _{k=1}^{N_{\text{el}}^c}({\sigma }_k^{(t)}{)}^{-2}{\left({\mathbf{z}}_c^{(i)}-\mathbf{\Phi }({\mathbf{\lambda }}_f^{(i)}){\mathbf{\theta }}_c^{(t)}\right)}_k^2\\ & -\frac{1}{2}\sum _{j=1}^{N_{\text{el}}^f}\log s_j^{(t)}-\frac{1}{2}\sum _{j=1}^{N_{\text{el}}^f}{s}_j^{(t)}{\left({\mathbf{T}}_f^{(i)}-{\mathbf{W}}^{(t)}{\mathbf{T}}_c({\mathbf{z}}_c^{(i)})\right)}_j^2,\end{array}$$

$$\frac{\partial \mathcal{F}}{\partial {\sigma }_k^{-2}}$$

$${\nabla }_{{\mathbf{\theta }}_c}\mathcal{F}$$

$$s_j^{(t+1)}=\frac{1}{N}\sum _{i=1}^N{⟨{\left({\mathbf{T}}_f^{(i)}-\mathbf{W}{\mathbf{T}}_c({\mathbf{z}}_c^{(i)})\right)}_j^2⟩}_{q_i^{(t)}}.$$

$${\lambda }_{\text{eff}}={\lambda }_{\text{mat}}\frac{{\lambda }_{\text{mat}}+{\lambda }_{\text{inc}}+{\varphi }_{\text{inc}}({\lambda }_{\text{inc}}-{\lambda }_{\text{mat}})}{{\lambda }_{\text{mat}}+{\lambda }_{\text{inc}}-{\varphi }_{\text{inc}}({\lambda }_{\text{inc}}-{\lambda }_{\text{mat}})},$$

$${\lambda }_{\text{eff}}=\frac{\alpha +\sqrt{{\alpha }^2+4{\lambda }_{\text{mat}}{\lambda }_{\text{inc}}}}{2},\alpha ={\lambda }_{\text{mat}}(2{\varphi }_{\text{mat}}-1)+{\lambda }_{\text{inc}}(2{\varphi }_{\text{inc}}-1).$$

$$(1-{\varphi }_{\text{inc}})\frac{d}{d{\varphi }_{\text{inc}}}{\lambda }_{\text{eff}}({\varphi }_{\text{inc}})=2{\lambda }_{\text{eff}}({\varphi }_{\text{inc}})\frac{{\lambda }_{\text{inc}}-{\lambda }_{\text{eff}}({\varphi }_{\text{inc}})}{{\lambda }_{\text{inc}}+{\lambda }_{\text{eff}}({\varphi }_{\text{inc}})},$$

$$\left(\frac{{\lambda }_{\text{inc}}-{\lambda }_{\text{eff}}}{{\lambda }_{\text{inc}}-{\lambda }_{\text{mat}}}\right)\sqrt{\frac{{\lambda }_{\text{mat}}}{{\lambda }_{\text{eff}}}}=1-{\varphi }_{\text{inc}}.$$

$$\log p({\mathbf{\theta }}_c)=\log \frac{\sqrt{\gamma }}{2}-\sqrt{\gamma }\sum _{i=1}^{N_{{\mathbf{\theta }}_c}}\left|{\theta }_{c,i}\right|,$$

Full text

We establish a physically inspired, data-driven coarse-grained model which learns a low-dimensional set of microstructural features that are predictive of the fine-grained model (FG) response. Once learned, those features provide a sharp distribution over the coarse scale effective coefficients of the PDE that are most suitable for prediction of the fine scale model output.

This ultimately allows to replace the computationally expensive FG by a generative probabilistic model based on evaluating the much cheaper CG several times. Sparsity enforcing priors further increase predictive efficiency and reveal microstructural features that are important in predicting the FG response. Moreover, the model yields probabilistic rather than single-point predictions, which enables the quantification of the unavoidable epistemic uncertainty that is present due to the information loss that occurs during the coarse-graining process.

\heading

Constantin Grigo and Phaedon-Stelios Koutsourelakis

1 INTRODUCTION

Many engineering design problems such as flow in porous media and mechanical properties of composite materials require simulations that are capable of resolving the microstructure of the underlying medium. If the material components under consideration exhibit fine-scale heterogeneity, popular discretization schemes (e.g. finite elements) yield very large systems of algebraic equations. Pertinent solution strategies at best (e.g. multigrid methods) scale linearly with the dimension of the unknown state vector. Despite the ongoing improvements in computer hardware, repeated solutions of such problems, as is required in the context of uncertainty quantification (UQ), poses insurmountable difficulties.

It is obvious that viable strategies for such problems, as well as a host of other deterministic problems where repeated evaluations are needed such as inverse, control/design problems etc, should focus on methods that exhibit sublinear complexity with respect to the dimension of the original problem. In the context of UQ a popular and general such strategy involves the use of surrogate models or emulators which attempt to learn the input-output map implied by the fine-grained model. Such models, e.g. Gaussian Processes [26], polynomial chaos expansions [8], neural nets [2] and many more, are trained on a finite set of fine-grained model runs. Nevertheless, their performance is seriously impeded by the curse of dimensionality, i.e. they usually become inaccurate for input dimensions larger than a few tens or hundreds, or equivalently, the number of FG runs required to achieve an acceptable level of accuracy grows exponentially fast with the input dimension.

Alternative strategies for high-dimensional problems make use of multi-fidelity models [10, 16] as inexpensive predictors of the FG output. As shown in [12], lower-fidelity models whose output deviates significantly from that of the FG can still yield accurate estimates with significant computational savings, as long as the outputs of the models exhibit statistical dependence. In the case of PDEs where finite elements are employed as the FG, multi-fidelity solvers can be simply obtained by using coarser discretizations in space or time. While linear and nonlinear dimensionality reduction techniques are suitable for dealing with high-dimensional inputs [14], it is known which of the microstructural features are actually predictive of FG outputs [22].

The model proposed in the present paper attempts to address this question. By using a two-component Bayesian network, we are able to predict fine-grained model outputs based on only a finite number of training data runs and a repeated solution of a much coarser model. Uncertainties can be easily quantified as our model leads to probabilistic rather than point-like predictions.

2 THE FINE-GRAINED MODEL

Let $\left(\Omega,\mathcal{F},P\right)$ be a probability space. Let $\mathcal{H}$ be the Hilbert space of functions defined over the domain $\mathcal{D}$ over which the physical problem is defined. We consider problems in the context of heterogeneous media which exhibit properties given by a random process $\lambda(\boldsymbol{x},\xi)$ defined over the product space $\mathcal{D}\times\Omega$. The corresponding stochastic PDE may be written as

[TABLE]

where $\mathcal{L}$ is a stochastic differential operator and $\boldsymbol{x}\in\mathcal{D}$, $\xi\in\Omega$ are elements of the physical domain and the sample space, respectively. Discretization of the random process

$\lambda(\boldsymbol{x},\xi)\underbrace{\longrightarrow}_{\textrm{discretize}}\boldsymbol{\lambda}_{f}\in\mathbb{R}^{n_{\boldsymbol{\lambda}_{f}}}$ as well as the governing equation leads to a system of $n_{f}$ (potentially nonlinear) algebraic equations, which can be written in residual form as

[TABLE]

where $\boldsymbol{U}_{f}(\boldsymbol{\lambda}_{f})\in\mathbb{R}^{n_{f}}$ is the $n_{f}$-dimensional discretized solution vector for a given $\boldsymbol{\lambda}_{f}$ and

$\boldsymbol{r}_{f}:\mathbb{R}^{n_{f}}\times\mathbb{R}^{n_{\boldsymbol{\lambda}_{f}}}\to\mathbb{R}^{n_{f}}$ the discretized residual vector. It is the model described by equation (2) which is denoted as the fine-grained model (FG) henceforth.

3 A GENERATIVE BAYESIAN SURROGATE MODEL

Let

[TABLE]

be the density of $\boldsymbol{\lambda}_{f}$. The density of the fine-scale response $\boldsymbol{U}_{f}$ is then given by

[TABLE]

where the conditional density $p(\boldsymbol{U}_{f}|\boldsymbol{\lambda}_{f})$ degenerates to a $\delta(\boldsymbol{U}_{f}-\boldsymbol{U}_{f}(\boldsymbol{\lambda}_{f}))$ when the only uncertainties in the problem are due to $\boldsymbol{\lambda}_{f}$.

The objective of this paper is to approximate this input-output map implied by $\boldsymbol{U}_{f}(\boldsymbol{\lambda}_{f})$, or equivalently in terms of probability densities, the conditional distribution $p(\boldsymbol{U}_{f}|\boldsymbol{\lambda}_{f})$. The latter case can also account for problems where additional sources of uncertainty are present and the input-output map is stochastic. To that end, we introduce a coarse-grained model (CG) leading to an approximate distribution $\bar{p}(\boldsymbol{U}_{f}|\boldsymbol{\lambda}_{f})$ which will be trained on a limited number of FG solutions $\mathcal{D}_{N}=\left\{\boldsymbol{\lambda}_{f}^{(i)},\boldsymbol{U}_{f}^{(i)}(\boldsymbol{\lambda}_{f}^{(i)})\right\}_{i=1}^{N}$.

Our approximate model $\bar{p}(\boldsymbol{U}_{f}|\boldsymbol{\lambda}_{f})$ employs a set of latent [3], reduced, collective variables which we denote by $\boldsymbol{\lambda}_{c}\in\mathbb{R}^{n_{\boldsymbol{\lambda}_{c}}}$ for reasons that will be apparent in the sequel, such that

[TABLE]

As it can be understood from the equation above, the latent variables $\boldsymbol{\lambda}_{c}$ act as a probabilistic filter (encoder) on the FG input $\boldsymbol{\lambda}_{f}$, retaining the features necessary for predicting the FG output $\boldsymbol{U}_{f}$. In order for $\bar{p}(\boldsymbol{U}_{f}|\boldsymbol{\lambda}_{f})$ to approximate well ${p}(\boldsymbol{U}_{f}|\boldsymbol{\lambda}_{f})$, the latent variables $\boldsymbol{\lambda}_{c}$ should not simply be the outcome of a dimensionality reduction on $\boldsymbol{\lambda}_{f}$. Even if $\boldsymbol{\lambda}_{f}$ is amenable to such a dimensionality reduction, it is not necessary that the $\boldsymbol{\lambda}_{c}$ found would be predictive of $\boldsymbol{U}_{f}$. Posed differently, it is not important that $\boldsymbol{\lambda}_{c}$ provides a high-fidelity encoding of $\boldsymbol{\lambda}_{f}$ but it suffices that it is capable of providing a good prediction of the corresponding $\boldsymbol{U}_{f}(\boldsymbol{\lambda}_{f})$.

The aforementioned desiderata do not unambiguously define the form of the encoding/ decoding densities in (5) nor the type/dimension of the latent variables $\boldsymbol{\lambda}_{c}$. In order to retain some of the physical and mathematical structure of the FG, we propose employing a coarsened, discretized version of the original continuous equation (1). In residual form this can again be written as

[TABLE]

where $\boldsymbol{U}_{c}\in\mathbb{R}^{n_{c}}$ is the $n_{c}$-dimensional ($n_{c}\ll n_{f}$) discretized solution and $\boldsymbol{r}_{c}:\mathbb{R}^{n_{f}}\times\mathbb{R}^{n_{\boldsymbol{\lambda}_{c}}}\to\mathbb{R}^{n_{c}}$ is the discretized residual vector. Due to the significant discrepancy in dimensions $n_{c}\ll n_{f}$, the cost of solving the CG in (6) is negligible compared to the FG in (2).

It is clear that $\boldsymbol{\lambda}_{c}$ plays the role of effective/equivalent properties but it is not obvious (except for some limiting cases where homogenization results can be invoked [23]) how these should depend on the fine-scale input $\boldsymbol{\lambda}_{f}$ nor how the solution $\boldsymbol{U}_{c}(\boldsymbol{\lambda}_{c})$ of the CG should relate to $\boldsymbol{U}_{f}(\boldsymbol{\lambda}_{f})$ in (2). Furthermore, it is important to recognize a priori that the use of the reduced variables $\boldsymbol{\lambda}_{c}$ in combination with the coarse model in (6) would in general imply some information loss. The latter should introduce an additional source of uncertainty in the predictions of the fine-scale output $\boldsymbol{U}_{f}$ [25], even in the limit of infinite training data. For that purpose and in agreement with (5), we propose a generative probabilistic model composed of the following two densities:

• •

A probabilistic mapping from $\boldsymbol{\lambda}_{f}$ to $\boldsymbol{\lambda}_{c}$, which determines the effective properties $\boldsymbol{\lambda}_{c}$ given $\boldsymbol{\lambda}_{f}$. We write this as $p_{c}(\boldsymbol{\lambda}_{c}|\boldsymbol{\lambda}_{f},\boldsymbol{\theta}_{c})$ where $\boldsymbol{\theta}_{c}$ denotes a set of model parameters,

• •

A coarse-to-fine map $p_{cf}(\boldsymbol{U}_{f}|\boldsymbol{U}_{c},\boldsymbol{\theta}_{cf})$, which is the PDF of the FG output $\boldsymbol{U}_{f}$ given the output $\boldsymbol{U}_{c}$ of the CG. It is parametrized by $\boldsymbol{\theta}_{cf}$.

This model is illustrated graphically in Figure 1. The density $p_{c}$ encodes $\boldsymbol{\lambda}_{f}$ into $\boldsymbol{\lambda}_{c}$ and the coarse-to-fine map $p_{cf}$ plays the role of a decoder, i.e. given the CG output $\boldsymbol{U}_{c}$, it predicts $\boldsymbol{U}_{f}$.

Using the abbreviated notation $\boldsymbol{\theta}=[\boldsymbol{\theta}_{c},\boldsymbol{\theta}_{cf}]$, from (5) we obtain

[TABLE]

where $\boldsymbol{U}_{c}(\boldsymbol{\lambda}_{c})$ is the solution vector to equation (6). The previous discussion suggests the following generative process for drawing samples from $\bar{p}(\boldsymbol{U}_{f}|\boldsymbol{\lambda}_{f},\boldsymbol{\theta})$ i.e. predicting the FG output $\boldsymbol{U}_{f}$ given a FG input $\boldsymbol{\lambda}_{f}$,

• •

draw a sample $\boldsymbol{\lambda}_{c}\sim p_{c}(\boldsymbol{\lambda}_{c}|\boldsymbol{\lambda}_{f},\boldsymbol{\theta}_{c})$,

• •

solve the CG to obtain $\boldsymbol{U}_{c}(\boldsymbol{\lambda}_{c})$,

• •

draw a sample $\boldsymbol{U}_{f}\sim p_{cf}(\boldsymbol{U}_{f}|\boldsymbol{U}_{c}(\boldsymbol{\lambda}_{c}),\boldsymbol{\theta}_{cf})$.

3.1 Model training

In order to train the model described above, it is a reasonable strategy to minimize the Kullback-Leibler divergence [5] between the target density $p(\boldsymbol{U}_{f},\boldsymbol{\lambda}_{f})=\delta(\boldsymbol{U}_{f}-\boldsymbol{U}_{f}(\boldsymbol{\lambda}_{f}))$ and $\bar{p}(\boldsymbol{U}_{f}|\boldsymbol{\lambda}_{f},\boldsymbol{\theta})$. As these are conditional distributions, the KL divergence would depend on $\boldsymbol{\lambda}_{f}$. In order to calibrate the model for the $\boldsymbol{\lambda}_{f}$ values that are of significance, we operate on the augmented densities $p(\boldsymbol{U}_{f},\lambda_{f})=p(\boldsymbol{U}_{f}|\boldsymbol{\lambda}_{f})p(\boldsymbol{\lambda}_{f})$ and $\bar{p}(\boldsymbol{U}_{f},\boldsymbol{\lambda}_{f}|\boldsymbol{\theta})=\bar{p}(\boldsymbol{U}_{f}|\boldsymbol{\lambda}_{f},\boldsymbol{\theta})p(\boldsymbol{\lambda}_{f})$, where $p(\boldsymbol{\lambda}_{f})$ is defined by (3). In particular, we propose minimizing with respect to $\boldsymbol{\theta}$

[TABLE]

where $N$ is the number of training samples drawn from $p(\boldsymbol{U}_{f},\boldsymbol{\lambda}_{f})$, i.e.

[TABLE]

and $H(p(\boldsymbol{U}_{f},\boldsymbol{\lambda}_{f}))$ is the entropy of $p(\boldsymbol{U}_{f},\boldsymbol{\lambda}_{f}))$ that is nevertheless independent of the model parameters $\boldsymbol{\theta}$. It is obvious from the final expression in (8) that $\sum_{i=1}^{N}\log\bar{p}(\boldsymbol{U}_{f}^{(i)}|\boldsymbol{\lambda}_{f}^{(i)},\boldsymbol{\theta}_{c},\boldsymbol{\theta}_{cf})$ is a $\log$-likelihood function of the data $\mathcal{D}_{N}$ which we denote by $L(\mathcal{D}_{N}|\boldsymbol{\theta}_{c},\boldsymbol{\theta}_{cf})$. In a fully Bayesian setting, this can be complemented with a prior $p(\boldsymbol{\theta}_{c},\boldsymbol{\theta}_{cf})$ leading to the posterior

[TABLE]

It is up to the analyst if predictions using equation (7) are carried out using point estimates of $\boldsymbol{\theta}$ (e.g. maximum likelihood (MLE) or maximum a posteriori (MAP)) or if a fully Bayesian approach is followed by averaging over the posterior $p(\boldsymbol{\theta}_{c},\boldsymbol{\theta}_{cf}|\mathcal{D}_{N})$. The latter has the added advantage of quantifying the uncertainty introduced due the finite training data $N$. We pursue the former in the following as it is computationally more efficient.

3.1.1 Maximizing the posterior

Our objective is to find $\boldsymbol{\theta}^{*}=[\boldsymbol{\theta}_{c}^{*},\boldsymbol{\theta}_{cf}^{*}]$ which maximizes the posterior given in equation (10), i.e.

[TABLE]

The main difficulty in this optimization problem arises from the log-likelihood term which involves the log of an analytically intractable integral with respect to $\boldsymbol{\lambda}_{c}$ since

[TABLE]

where we note that an independent copy of $\boldsymbol{\lambda}_{c}^{(i)}$ pertains to each data point $i$. Due to this integration, typical deterministic optimization algorithms are not applicable.

However, as $\boldsymbol{\lambda}_{c}$ appears as a latent variable, we may resort to the well-known Expectation-Maximization algorithm [6]. Using Jensen’s inequality, we establish a lower bound on every single term $L_{i}$ of the sum in the log-likelihood $L(\mathcal{D}_{N}|\boldsymbol{\theta}_{c},\boldsymbol{\theta}_{cf})$ by employing an arbitrary density $q_{i}(\boldsymbol{\lambda}_{c}^{(i)})$ (different for each sample $i$) as

[TABLE]

Hence, the $\log$-posterior in (10) can be lower bounded by

[TABLE]

The introduction of the auxiliary densities $q_{i}$ suggests the following recursive procedure [4] for maximizing the log-posterior:

• E-step:

Given some $\boldsymbol{\theta}_{c}^{(t)},\boldsymbol{\theta}_{cf}^{(t)}$ in iteration $t$, find the auxiliary densities $q_{i}^{(t+1)}$ that maximize $\mathcal{F}$,

• M-step:

Given $q_{i}^{(t+1)}$, find the parameters $\boldsymbol{\theta}_{c}^{(t+1)},\boldsymbol{\theta}_{cf}^{(t+1)}$ that maximize $\mathcal{F}$.

It can be readily shown that the optimal $q_{i}$ is given by

[TABLE]

with which the inequality in (13) becomes an equality. In fact, both E- and M-steps can be relaxed to find suboptimal $q^{(t)}_{i},\boldsymbol{\theta}_{c}^{(t)},\boldsymbol{\theta}_{cf}^{(t)}$, which enables the application of approximate schemes such as e.g. Variational Inference (VI) [24, 9]. For the M-step, we may resort to any (stochastic) optimization algorithm [17] or, on occasion, closed-form updates might also be feasible.

4 SAMPLE PROBLEM: 2D STATIONARY HEAT EQUATION

As a sample problem, we consider the 2D stationary heat equation on the unit square $[0,1]^{2}$

[TABLE]

where $U(\boldsymbol{x},\xi(\boldsymbol{x}))$ represents the temperature field. For the boundary conditions, we fix the temperature $U$ in the upper left corner (see Figure 2) to $-50$ and prescribe the heat flux $\boldsymbol{Q}(\boldsymbol{x})=\left(\begin{smallmatrix}150-30y\\ 100-30x\end{smallmatrix}\right)$ on the remaining domain boundary $\partial\mathcal{D}$.

We consider a binary random medium whose conductivity $\lambda(\boldsymbol{x})$ can take the values $\lambda_{\textrm{hi}}$ and $\lambda_{\textrm{lo}}$. To define such a field we consider transformations of a zero-mean Gaussian process $\xi(\boldsymbol{x})$ of the form [18, 13]

[TABLE]

where the thresholding constant $c$ is selected so as to achieve the target volume fraction $\phi_{\textrm{hi}}$ (or equivalently $\phi_{\textrm{lo}}=1-\phi_{\textrm{hi}}$) of the material with conductivity $\lambda_{\textrm{hi}}$ (or equivalently of $\lambda_{\textrm{lo}}$). In the following we use an isotropic squared-exponential covariance function for $\xi(\boldsymbol{x})$ of the form

[TABLE]

In the following studies, we used values $l\approx 0.01$.

Due to the small correlation length $l$, we discretize the SPDE in Equation (16) using $256\times 256$ standard quadrilateral finite elements, leading to a linear system of $N_{\textrm{eq}}=\dim(\boldsymbol{U}_{f})-1=257\times 257-1=66048$ algebraic equations. We choose the discretization mesh of the random process $\lambda(\boldsymbol{x},\xi(\boldsymbol{x}))$ to coincide with the finite element discretization mesh of the SPDE, i.e. $\dim(\boldsymbol{\lambda}_{f})=256\times 256=65536$.

Samples $\mathcal{D}_{N}=\{\lambda_{f}^{(i)},\boldsymbol{U}_{f}^{(i)}\}_{i=1}^{N}$ are readily obtained by simulating realizations of the discretized Gaussian field, transforming them according to (17) and solving the discretized SPDE. Three such samples can be seen in Figure 2.

4.1 Model specifications

We define a coarse model employing $n_{\boldsymbol{\lambda}_{c}}$ quadrilateral elements, the conductivities of which are given by the vector $\boldsymbol{\lambda}_{c}$. Since these need to be strictly positive, we operate instead on $\boldsymbol{z}_{c}$ defined as

[TABLE]

For each element $k=1,\ldots,n_{\boldsymbol{\lambda}_{c}}$ of the coarse model/mesh, we assume that

[TABLE]

where $\boldsymbol{\lambda}_{f,k}$ is the subset of the vector $\boldsymbol{\lambda}_{f}$ that belongs to coarse element $k$ and $\chi_{j}$ some feature functions of the underlying $\boldsymbol{\lambda}_{f,k}$ that we specify below. In a more compact form, we can write

[TABLE]

where $\boldsymbol{\Phi}(\boldsymbol{\lambda}_{f})$ is an $n_{\boldsymbol{\lambda}_{c}}\times N_{\textrm{features}}$ design matrix with $\Phi_{kj}(\boldsymbol{\lambda}_{f})=\chi_{j}(\boldsymbol{\lambda}_{f,k})$ and $\textrm{diag}(\sigma^{2})$ is a diagonal covariance matrix with components $\sigma_{k}^{2}$.

For the coarse-to-fine map $p_{cf}$, we postulate the relation

[TABLE]

where $\boldsymbol{\theta}_{cf}=(\boldsymbol{W},\boldsymbol{S})$ are parameters to be learned from the data. The matrix $\boldsymbol{W}$ is of dimension $n_{f}\times n_{c}$ and $\boldsymbol{S}$ is a positive definite covariance matrix of size $n_{f}\times n_{f}$. To reduce the large amount of free parameters in the model, we fix $\boldsymbol{W}$ to express coarse model’s shape functions. Furthermore, we assume that $\boldsymbol{S}=\textrm{diag}(\boldsymbol{s})$ where $\boldsymbol{s}$ is the $n_{f}$-dimensional vector of diagonal entries of the diagonal matrix $\boldsymbol{S}$. The aforementioned expression implies, on average, a linear relation between the fine and coarse model outputs, which is supplemented by the residual Gaussian noise implied by $\boldsymbol{S}$. The latter expresses the uncertainty in predicting the FG output when only the CG solution is available. We note that the relation in Equation (20) (i.e. $\boldsymbol{\theta}_{c}$ and $\sigma_{k}^{2}$) will be adjusted during training so that the model implied in Equation (22) represents the data as good as possible.

From equation (15), we have that

[TABLE]

where with $(.)_{i}$, we mean the $i$-th component of the vector in brackets. For use in the M-step, we compute the gradients

[TABLE]

where $\boldsymbol{\Sigma}=\textrm{diag}(\sigma_{1}^{2},\ldots,\sigma^{2}_{n_{\boldsymbol{\lambda}_{c}}})$. We can readily compute the roots to find the optimal $\boldsymbol{\Sigma},\boldsymbol{\theta}_{c}$. Furthermore, from$\frac{\partial\mathcal{F}}{\partial s_{j}}=0$ we get

[TABLE]

4.2 Feature functions $\chi$

The framework advocated allows for any form and any number of feature functions $\chi_{j}$ in (20). Naturally such a selection can be guided by physical insight [11]. In practice, the feature functions we are using can be roughly classified into two different groups:

4.2.1 Effective-medium approximations

We consider existing effective-medium approximation formula that can be found in the literature [23] as ingredients for construction of feature functions $\chi_{j}$ in equation (20). The majority of commonly used such features only retain low-order topological information. only The following approximations to the effective property $\lambda_{\textrm{eff}}$ can be used as building blocks for feature functions $\chi$ in the sense that we can transform them nonlinearly such that $\chi(\boldsymbol{\lambda}_{f})=f(\lambda_{\textrm{eff}}(\boldsymbol{\lambda}_{f}))$. In particular, as we are modeling $\boldsymbol{z}_{c}=\log\boldsymbol{\lambda}_{c}$, we include feature functions of type $\chi(\boldsymbol{\lambda}_{f})=\log(\lambda_{\textrm{eff}}(\boldsymbol{\lambda}_{f}))$.

Maxwell-Garnett approximation (MGA)

The Maxwell-Garnett approximation is assumed to be valid if the microstructure consists of a matrix phase $\lambda_{\textrm{mat}}$ and spherical inclusions $\lambda_{\textrm{inc}}$ that are small and dilute enough such that interactions between them can be neglected. another Moreover, the heat flux far away from any inclusion is assumed to be constant. Under such conditions, the effective conductivity can be approximated in 2D as

[TABLE]

where the volume fraction $\phi_{i}(\boldsymbol{\lambda}_{f})$ is given by the fraction of phase $i$ elements in the binary vector $\boldsymbol{\lambda}_{f}$.

Self-consistent approximation (SCA)

The self-consistent approximation (or Bruggeman

formula) was originally developed for effective electrical properties of random microstructures. It also considers non-interacting spherical inclusions and follows from the assumption that perturbations of the electric field due to the inclusions average to [math]. In 2D, the formula reads

[TABLE]

Note that the SCA exhibits phase inversion symmetry.

Differential effective medium (DEM)

From a first-order expansion in volume fraction of the effective conductivity in the dilute limit of spherical inclusions, one can deduce the differential equation [23]

[TABLE]

which can be integrated to

[TABLE]

We solve for $\lambda_{\textrm{eff}}$ and use it as a feature function $\chi$.

4.2.2 Morphology-describing features

Apart from the effective-medium approximations which only take into account the phase conductivities and volume fractions, we wish to have feature functions $\chi_{j}$ that more thoroughly describe the morphology of the underlying microstructure. Popular members of this class of microstructural features are the two-point correlation, the lineal-path function, the pore-size density or the specific surface to mention only a few of them [23].

We are however free to use any function $\chi_{j}:(\mathbb{R}^{+})^{\dim(\boldsymbol{\lambda}_{f,k})}\mapsto\mathbb{R}$ as a feature, no matter from which field or consideration it may originate. We thus make use of existing image processing features [15] as well as novel topology-describing functions. Some important examples are

Convex area of connected phase blobs

This feature identifies distinct, connected “blobs” of only high or low conducting phase pixels and computes the area of the convex hull to each blob. We then can e.g. use the mean or maximum value thereof as a feature function.

Blob extent

From an identified phase blob, we can compute its maximum extension in $x$- and $y$-directions. One can for instance take the mean or maximum of maximum extension among identified blobs as a feature.

Distance transformation functions

The distance transform of a binary image assigns a number to each pixel $i$ that is the distance from that pixel $i$ to the nearest nonzero pixel in the binary image, see right part of Figure 3. One can use either phase to correspond to nonzero in the binary image as well as different distance metrics. As a feature, one can e.g. take the mean or maximum of the distance transformed image.

Pixel-cross function

This feature counts the number of high or low conducting pixels one has to cross going on a straight line from boundary to boundary in $x$- or $y$-direction. One can again take e.g. mean, maximum or minimum values as the feature function outputs.

Straight path mean function

A further refinement of the latter function is to take generalized means instead of numbers of crossed pixels along straight lines from boundary to boundary. In particluar, we use harmonic, geometric and arithmetic means as features.

4.3 Sparsity priors

It is clear from the previous discussion that the number of feature functions is practically limitless. The more such $\chi_{j}$ one introduces, the more unknown parameters $\theta_{c,j}$ must be learned. From the modeling point of view, while ML estimates can always be found, it is desirable to have as clear of a distinction as possible between relevant and irrelevant features that could provide further insight as well being able to do so with the fewest possible training data available. For that purpose we advocate the use of sparsity-enforcing priors in $\boldsymbol{\theta}_{c}$ [1, 7, 19]. From a statistical perspective, this is also motivated by the bias-variance-tradeoff. Model prediction accuracy is adversely affected by two factors: One is noise in the training data (variance), the other is due to overly simple model assumptions (bias). Maximum-likelihood estimates of model parameters tend to have low bias but high variance, i.e. they accurately predict the training data but generalize poorly. To address this issue, a common Bayesian approach to control model complexity is the use of priors, which is the equivalent to regularization in frequentist formulations. A particularly appealing family of prior distributions is the Laplacian (or LASSO regression, [21]), as it sets redundant or unimportant predictors to exactly 0, thereby simplifying interpretation.

In particular, we use a prior on the coefficients $\boldsymbol{\theta}_{c}$ of the form

[TABLE]

with a hyper-parameter $\gamma$ which can be set by either applying some cross-validation scheme or more efficiently by minimization of Stein’s unbiased risk estimate [20, 27]. A straightforward implementation of this prior is described in [7].

5 RESULTS

5.1 Required training data

We use a total number of 306 feature functions $\chi_{j}$ and a Laplacian prior (31) with a hyperparameter $\gamma$ we find by cross-validation. We set the length scale parameter to $l=0.0781$ and the expected volume fraction of the high conducting phase to $\phi_{\textrm{hi}}=0.2$. For the contrast, we take $c=\frac{\lambda_{\textrm{hi}}}{\lambda_{\textrm{lo}}}=10$, where we set $\lambda_{\textrm{lo}}=1$. The full- and reduced-order models are both computed on regular quadrilateral finite element meshes of size $256\times 256$ for the FG and $2\times 2$, $4\times 4$ and $8\times 8$ for the CG.

Our goal is to measure the predictive capabilities of the described model. To that end, we compute the mean squared distance

[TABLE]

of the predictive mean $\left<\boldsymbol{U}_{f}^{(i)}\right>_{\bar{p}(\boldsymbol{U}_{f}^{(i)})|\boldsymbol{\lambda}_{f}^{(i)},\boldsymbol{\theta})}$ to the true FG output $\boldsymbol{U}_{f}^{(i)}$ on a test set of $N_{\textrm{test}}=256$ samples. Predictions are carried out by drawing 10,000 samples of $\boldsymbol{\lambda}_{c}$ from the learned distribution $p_{c}(\log\boldsymbol{\lambda}_{c}|\boldsymbol{\lambda}_{f},\boldsymbol{\theta}_{c}^{*},\boldsymbol{\sigma})=\mathcal{N}(\log(\boldsymbol{\lambda}_{c})|\boldsymbol{\Phi}(\boldsymbol{\lambda}_{f})\boldsymbol{\theta}_{c}^{*},\textrm{diag}(\boldsymbol{\sigma}^{2}))$, solving the coarse model $\boldsymbol{U}_{c}^{(i)}=\boldsymbol{U}_{c}(\boldsymbol{\lambda}_{c}^{(i)})$ and drawing a sample from $p_{cf}(\boldsymbol{U}_{f}^{(i)}|\boldsymbol{U}_{c}^{(i)},\boldsymbol{\theta}_{cf}^{*})=\mathcal{N}(\boldsymbol{U}_{f}^{(i)}|\boldsymbol{W}\boldsymbol{U}_{c}^{(i)},\boldsymbol{S}^{*})$ for every test data point. Monte Carlo noise is small enough to be neglected.

As a reference value for the computed error $d^{2}$, we compute the mean variance of the FG output

[TABLE]

where expectation values $\left<.\right>$ are estimated on a set of 1024 FG samples such that errors due to Monte Carlo can be neglected. Figure 4 shows the relative squared prediction error $\frac{d^{2}}{var(U_{f})}$ in dependence of the number of training data samples for different coarse mesh sizes $\dim(\boldsymbol{\lambda}_{c})=2\times 2,4\times 4$ and $8\times 8$. We observe that the predictive error converges to a finite value already for relatively few training data. This is due to the inevitable information loss during the coarsening process $\boldsymbol{\lambda}_{f}\rightarrow\boldsymbol{\lambda}_{c}$ as well as finite element discretization errors. In accordance with that, we see that the error drops with the dimension of the coarse mesh, $\dim(\boldsymbol{\lambda}_{c})$.

5.2 Activated microstructural feature functions

For the same volume fraction and microstructural length scale parameter as above

($\phi_{\textrm{hi}}=0.2,l=0.0781$), a varying contrast of $c=\left\{10,20,50\right\}$, a coarse model dimension of $\dim(\boldsymbol{\lambda}_{c})=4\times 4$ and a training set of $N=128$ samples, we find the optimal $\boldsymbol{\theta}_{c}^{*}$ shown in Figure 5.

Due to the application of a sparsity enforcing prior as described in section 4.3, we observe that most components of $\boldsymbol{\theta}_{c}^{*}$ are exactly 0. Comparability between different feature functions can be ensured by standardization or normalization of feature function outputs on the training data. For all three contrast values, we see that the three most important feature functions are given by the maximum convex area of blobs of conductivity $\lambda_{\textrm{hi}}$, the maximum $\log$ of the geometric mean of conductivities along a straight line from boundary to boundary in $y$-direction and the $\log$ of the self-consistent effective medium approximation as described above. The maximum convex area feature returns the largest convex area of all blobs found within a coarse element. convex area. It is an interesting to note that although, in the set of 306 feature functions $\chi_{j}$, the $(\max/\min/$mean/variance of) the blob area (for both high and low conducting phases) are also included, it is only the convex area which is activated.

In Figure 5, it is observed that with increasing contrast, the coefficient $\theta_{c,j}^{*}$ belonging to the $\log$ self-consistent approximation is decreasing in contrast to increasing values of $\theta_{c,j}^{*}$’s corresponding the $\max.$ high conducting convex area and the $\max.$ $\log$ geometric mean along $y$-direction. We believe that this is because, the higher the contrast $c$ is, the more the exact geometry and connectedness of the microstructure plays a role for predicting the effective properties. As the SCA only considers the volume fraction and the conductivities of both phases, it disregards such information.

5.3 Learned effective property $\boldsymbol{\lambda}_{c}$

Figure 6 shows three test samples of $\phi_{\textrm{hi}}=0.2,l=0.0781$ and $c=10$ (top row) along with the mode $\exp\left(\boldsymbol{\Phi}(\boldsymbol{\lambda}_{f})\boldsymbol{\theta}_{c}^{*}\right)$ of the learned distribution of the effective conductivity $p_{c}(\boldsymbol{\lambda}_{c}|\boldsymbol{\lambda}_{f},\boldsymbol{\theta}_{c}^{*})$, which is connected to the mean by $\left<\boldsymbol{\lambda}_{c}\right>_{p_{c}}=\exp\left(\boldsymbol{\Phi}(\boldsymbol{\lambda}_{f})\boldsymbol{\theta}_{c}^{*}+\frac{1}{2}\boldsymbol{\sigma}^{*2}\right)$. We emphasize again that the latent variable $\boldsymbol{\lambda}_{c}$ is not a lower-dimensional compression of $\boldsymbol{\lambda}_{f}$ with the objective of most accurately reconstructing $\boldsymbol{\lambda}_{f}$, but of providing good predictions of $\boldsymbol{U}_{f}(\boldsymbol{\lambda}_{f})$. Even though Figure 6 gives the impression of a simple local averaging relation between $\boldsymbol{\lambda}_{f}$ and $\boldsymbol{\lambda}_{c}$, this is not always be the case. In particular, $p_{c}(\boldsymbol{\lambda}_{c}|\boldsymbol{\lambda}_{f},\boldsymbol{\theta}_{c}^{*})$ was found to have non-vanishing probability mass for $\lambda_{c,i}<\lambda_{\textrm{lo}}$ or $\lambda_{c,i}>\lambda_{\textrm{hi}}$, especially for more general models where the coarse-to-fine mapping (22) is not fixed to be the shape function interpolation $\boldsymbol{W}$.

5.4 Predictive uncertainty

The predictive uncertainty is composed of the uncertainty in having an accurate encoding of $\boldsymbol{\lambda}_{f}$ in $\boldsymbol{\lambda}_{c}$ which is described by $\sigma_{k}^{2}$ in $p_{c}$ (21), as well as the uncertainty in the reconstruction process from $\boldsymbol{U}_{c}$ to $\boldsymbol{U}_{f}$, which is given by the diagonal covariance $\boldsymbol{S}=\textrm{diag}(\boldsymbol{s})$ in $p_{cf}$ (22). Both are shown after training ($\phi_{2}=0.2,l=0.0781,c=10$) in Figure 7. For $\boldsymbol{\sigma}^{*2}$, we observe that values in the corner elements always converge to very tight values, whereas some non-corner elements can converge to comparably large values. The exact location of these elements is data dependent. The coarse-to-fine reconstruction variances $\boldsymbol{s}^{*}$ is depicted on the right column of Figure 7. As expected, we see that the estimated coarse-to-fine reconstruction error is largest in the center of coarse elements i.e. at large distances from the coarse model finite element nodes.

5.5 Predictions

As in section 5.1, for a model with $N_{\textrm{train}}=32$ and $\dim(\boldsymbol{\lambda}_{c})=8\times 8$ and random microstructures with parameter values $\phi_{\textrm{hi}}=0.2,l=0.0781,c=10$, we consider predictions by sampling from $\bar{p}(\boldsymbol{U}_{f}|\boldsymbol{\lambda}_{f},\boldsymbol{\theta}_{cf},\boldsymbol{\theta}_{c})$ using 10,000 samples. The predictive histogram for the temperature $U_{f,lr}$ in the lower right corner of the domain can be seen in Figure 8. Figure 9 shows a surface plot of the true solution (colored), the predictive mean (blue) $\pm\sigma$ (gray). As can be seen, the true solution $\boldsymbol{U}_{f}$ is nicely included in $\bar{p}$ everywhere.

6 CONCLUSION

In this paper, we described a generative Bayesian model which is capable of giving probabilistic predictions to an expensive fine-grained model (FG) based only on a small finite number of training data and multiple solutions of a fast, but less accurate reduced order model (CG).

In particular, we consider the discretized solution of stochastic PDEs with random coefficients where the FG corresponds to a fine-scale discretization. Naturally, this comes along with a very high-dimensional vector of input uncertainties. The proposed model is capable of extracting the most relevant features of those input uncertainties and gives a mapping to a much lower dimensional space of effective properties (encoding). These lower dimensional effective properties serve as the input to the CG, which solves the PDE on a much coarser scale. The last step consists of a probabilistic reconstruction mapping from the coarse- to the fine-scale solution (decoding).

We demonstrated features and capabilities of the model proposed for a $2D$ steady-state heat problem, where the fine scale of the conductivity implies (upon discretization) a random input vector of dimension $256\times 256$ and the solution of a discretized system of equations of comparable size. In combination with a sparsity-enforcing prior, the proposed model identified the most salient features of the fine-scale conductivity field and allowed accurate predictions of the FG response using a CGs of size only $2\times 2,4\times 4$ and $8\times 8$. The predictive distribution always included the true FG solution as well as provided uncertainty bounds arising from the information loss taking place during the coarse-graining process.