A physics-aware, probabilistic machine learning framework for coarse-graining high-dimensional systems in the Small Data regime

TL;DR

This paper introduces a physics-aware, probabilistic machine learning framework designed for effective coarse-graining of high-dimensional physical systems in scenarios with limited data, leveraging physical insights and variational inference.

Contribution

It presents a novel probabilistic ML approach that combines physical structure and contextual knowledge to perform coarse-graining with small datasets, including model refinement capabilities.

Findings

Performs well under extrapolative conditions

Unifies dimensionality reduction and model calibration

Provides probabilistic estimates of information loss

Abstract

The automated construction of coarse-grained models represents a pivotal component in computer simulation of physical systems and is a key enabler in various analysis and design tasks related to uncertainty quantification. Pertinent methods are severely inhibited by the high-dimension of the parametric input and the limited number of training input/output pairs that can be generated when computationally demanding forward models are considered. Such cases are frequently encountered in the modeling of random heterogeneous media where the scale of the microstructure necessitates the use of high-dimensional random vectors and very fine discretizations of the governing equations. The present paper proposes a probabilistic Machine Learning framework that is capable of operating in the presence of Small Data by exploiting aspects of the physical structure of the problem as well as contextual knowledge. As a result, it can perform comparably well under extrapolative conditions. It unifies the tasks of dimensionality and model-order reduction through an encoder-decoder scheme that simultaneously identifies a sparse set of salient lower-dimensional microstructural features and calibrates an inexpensive, coarse-grained model which is predictive of the output. Information loss is accounted for and quantified in the form of probabilistic predictive estimates. The learning engine is based on Stochastic Variational Inference. We demonstrate how the variational objectives can be used not only to train the coarse-grained model, but also to suggest refinements that lead to improved predictions.

Figures

Tables

Table 1. Table 1: Computational complexity of training and prediction stages. The scaling with dim ​ ( 𝝀 f ) dim subscript 𝝀 𝑓 \text{dim}(\boldsymbol{\lambda}_{f}) is dependent on the applied set of feature functions 𝝋 m ​ ( 𝝀 f ) subscript 𝝋 𝑚 subscript 𝝀 𝑓 \boldsymbol{\varphi}_{m}(\boldsymbol{\lambda}_{f}) and reaches from 𝒪 ​ ( 1 ) 𝒪 1 \mathcal{O}(1) (e.g. for φ m ​ ( 𝝀 f ) = 1 subscript 𝜑 𝑚 subscript 𝝀 𝑓 1 \varphi_{m}(\boldsymbol{\lambda}_{f})=1 ) to 𝒪 ​ ( ( dim ​ ( 𝝀 f ) ) 2 ) 𝒪 superscript dim subscript 𝝀 𝑓 2 \mathcal{O}\left((\text{dim}(\boldsymbol{\lambda}_{f}))^{2}\right) (e.g. φ m ​ ( 𝝀 f ) = subscript 𝜑 𝑚 subscript 𝝀 𝑓 absent \varphi_{m}(\boldsymbol{\lambda}_{f})= mean mutual distance of exclusion centers).

Quantity Phase	$N$	$\text{dim}\left({𝒖}_f\right)$	$\text{dim}\left({𝝀}_f\right)$	$\text{dim}\left({𝒖}_c\right)$	$\text{dim}\left({𝝀}_c\right)$	$\text{dim}\left({\widetilde{𝜽}}_{c,m}\right)$
training	$N$	$\text{dim}\left({𝒖}_f\right)$	$1\dots$ ${\left(\text{dim}\left({𝝀}_f\right)\right)}^2$	$\text{dim}\left({𝒖}_c\right)$	$\text{dim}\left({𝝀}_c\right)$	$\text{dim}\left({\widetilde{𝜽}}_{c,m}\right)$
prediction	$1$	$\text{dim}\left({𝒖}_f\right)$	$1\dots$ ${\left(\text{dim}\left({𝝀}_f\right)\right)}^2$	$\text{dim}\left({𝒖}_c\right)$	$\text{dim}\left({𝝀}_c\right)$	$\text{dim}\left({\widetilde{𝜽}}_{c,m}\right)$

Table 2. Table 2: Error metrics R 2 superscript 𝑅 2 R^{2} and M ​ L ​ L 𝑀 𝐿 𝐿 MLL corresponding to the cross-prediction example in section 3.4 using distinct boundary conditions 𝒂 𝒂 \boldsymbol{a} and 𝒂 ~ ~ 𝒂 \tilde{\boldsymbol{a}} . Metrics are evaluated and averaged over 3 separate training sets of N = 128 𝑁 128 N=128 size (indicated errors are due to variation of training data) and tested on N test = 1024 subscript 𝑁 test 1024 N_{\text{test}}=1024 microstructural samples. Only slight deterioration is observed if boundary conditions of the training and test sets differ.

trained on prediction on
	$0.930\pm 1.3\cdot {10}^{-3}$	$0.863\pm 1.1\cdot {10}^{-2}$
	$0.887\pm 7.5\cdot {10}^{-3}$	$0.903\pm 7.0\cdot {10}^{-3}$

Table 3. Table 3: Set of 150 feature functions φ 𝜑 \varphi applied in the numerical examples of section 3 .

Index $j$	Function ${\phi }_{jm}$	Comment
1	constant	${\phi }_j=1$
2	pore fraction in $\Omega$	pore fraction evaluated on full domain $\Omega$
3	$\log$ pore fraction in $\Omega$
4–7	(pore fraction)0.5…0.5…2.5 in $\Omega$
8	$\exp$(pore fraction)
9	$\log$ SCA in $\Omega$	$\log$ self-cons. approx. Torquato (2001), sec. 18.1.2, inf. contrast limit
10	Maxwell-Approximation in $\Omega$	inf. contrast limit, see Torquato (2001), sec. 18.2.1
11	$\log$ Maxwell-Approximation
12–17	$\log$ chord length density in $\Omega$, $d=0.05\dots 0$	see Torquato (2001) sec. 6.2.4
18	interface area in $\Omega$
19	$\log$ interface area in $\Omega$
20–27	$|\log$ interface area${|}^{3/2,1/2,1/3,1/4,1/5,2,3,4}$ in $\Omega$
28–31	(interface area)1/3,1/4,1/5,2 in $\Omega$
32	mean distance edge in $\Omega$	measured from excl. edge to edge
33	$\log$ mean distance edge in $\Omega$	measured from excl. edge to edge
34	${\log }^2$ mean distance edge in $\Omega$	measured from excl. edge to edge
35	${\log }^3$ mean distance edge in $\Omega$	measured from excl. edge to edge
36	mean distance center in $\Omega$	measured from excl. center to center
37	min. distance center in $\Omega$	measured from excl. center to center
38	$\log$ min. distance center in $\Omega$	measured from excl. center to center
39	${\log }^2$ min. distance center in $\Omega$	measured from excl. center to center
40–43	lineal path in $\Omega$ for $d=0.025,0.01,0.005,0.002$	see Torquato (2001), sec. 2.4
44–47	$\log$ lineal path in $\Omega$ for $d=0.025,0.01,0.005,0.002$
48	void nearest-neighbor pdf, $d=0$, in $\Omega$	see Torquato (2001), sec. 2.8
49	$\log$ void nearest-neighbor pdf, $d=0$, in $\Omega$	see Torquato (2001). sec. 2.8
50	pore size density, $d=0$, in $\Omega$	see Torquato (2001), sec. 2.6
51	$\log$ pore size density, $d=0$, in $\Omega$	see Torquato (2001), sec. 2.6
52	mean chord length in $\Omega$	see Torquato (2001), sec. 2.5
53	$\log$ mean chord length in $\Omega$	see Torquato (2001), sec. 2.5
54	$\exp$ mean chord length in $\Omega$	see Torquato (2001), sec. 2.5
55	(mean chord length)0.5 in $\Omega$	see Torquato (2001), sec. 2.5
56–58	in ${\Omega }_m$	expected exclusion radii moments
59–61	in ${\Omega }_m$	$\log$ expected exclusion radii moments
62	pore fraction in ${\Omega }_m$
63	$\log$ pore fraction in ${\Omega }_m$
64	$\exp$ pore fraction in ${\Omega }_m$
65–68	(pore fraction)0.5,1.5,2,2.5 in ${\Omega }_m$
69	$\log$ self-consistent approximation (inf. contrast) in ${\Omega }_m$
70	Maxwell Approximation in ${\Omega }_m$
71	$\log$ Maxwell Approximation in ${\Omega }_m$
72–78	$\log$ chord length dens. in ${\Omega }_m$, $d=(50,25,12.5,6.25,3,1.5,0)\cdot {10}^{-4}$
79	interface area in ${\Omega }_m$
80–88	$|\log$ interface area${|}^{1,1.5,2,3,4,5,1/2,1/3,1/4}$ in ${\Omega }_m$
89–93	$|$ interface area${|}^{1/2,1/3,1/4,1/5,2}$ in ${\Omega }_m$
94	mean distance edge in ${\Omega }_m$	measured from excl. edge to edge
95	$\log$ mean distance edge in ${\Omega }_m$	measured from excl. edge to edge
96	max distance edge in ${\Omega }_m$	measured from excl. edge to edge
97	$\log$ max distance edge in ${\Omega }_m$	measured from excl. edge to edge
98	std distance edge in ${\Omega }_m$	measured from excl. edge to edge
99	$\log$ std distance edge in ${\Omega }_m$	measured from excl. edge to edge
100–104	square well potential, width = $(1,2,3,4,5)\cdot {10}^{-2}$ in ${\Omega }_m$
105	${\log }^2$ mean distance in ${\Omega }_m$	measured from excl. edge to edge
106	${\log }^3$ mean distance in ${\Omega }_m$	measured from excl. edge to edge
107	mean distance center in ${\Omega }_m$	measured from excl. center to center
108	min. distance center in ${\Omega }_m$	measured from excl. center to center
109	$\log$ min. distance center in ${\Omega }_m$	measured from excl. center to center
110	${\log }^2$ min. distance center in ${\Omega }_m$	measured from excl. center to center
111–114	lineal path in ${\Omega }_m$ for $d=0.025,0.01,0.005,0.002$
115–118	$\log$ lineal path in ${\Omega }_m$ for $d=0.025,0.01,0.005,0.002$
119	void nearest-neighbor pdf, $d=0$, in ${\Omega }_m$	see Torquato (2001), sec. 2.8
120	$\log$ void nearest-neighbor pdf, $d=0$, in ${\Omega }_m$	see Torquato (2001). sec. 2.8
121	pore size density, $d=0$, in ${\Omega }_m$	see Torquato (2001), sec. 2.6
122	$\log$ pore size density, $d=0$, in ${\Omega }_m$	see Torquato (2001), sec. 2.6
123	mean chord length in ${\Omega }_m$	see Torquato (2001), sec. 2.5
124	$\log$ mean chord length in ${\Omega }_m$	see Torquato (2001), sec. 2.5
125	$\exp$ mean chord length in ${\Omega }_m$	see Torquato (2001), sec. 2.5
126	(mean chord length)0.5 in ${\Omega }_m$	see Torquato (2001), sec. 2.5
127–129	in ${\Omega }_m$	expected exclusion radii moments
130–132	in ${\Omega }_m$	$\log$ expected exclusion radii moments
133	length scale of exp. approx. to lin. path in ${\Omega }_m$
134	mean of euclidean dist. transform in ${\Omega }_m$	see Soille (1999)
135	variance of euclidean dist. transform in ${\Omega }_m$	see Soille (1999)
136	max. of euclidean dist. transform in ${\Omega }_m$	see Soille (1999)
137	mean of chessboard dist. transform in ${\Omega }_m$	see Soille (1999)
138	variance of chessboard dist. transform in ${\Omega }_m$	see Soille (1999)
139	max. of chessboard dist. transform in ${\Omega }_m$	see Soille (1999)
140	mean of cityblock dist. transform in ${\Omega }_m$	see Soille (1999)
141	variance of cityblock dist. transform in ${\Omega }_m$	see Soille (1999)
142	max. of cityblock dist. transform in ${\Omega }_m$	see Soille (1999)
143–145	Gauss lin. filt. $d=2,5,10$ pixels in ${\Omega }_m$	see Grigo and Koutsourelakis (2019)
146	Ising energy in ${\Omega }_m$
147	shortest connected path, $x$-dir., Euclidean, in ${\Omega }_m$	see Soille (1999)
148	shortest connected path, $y$-dir., Euclidean, in ${\Omega }_m$	see Soille (1999)
149	shortest connected path, $x$-dir., cityblock, in ${\Omega }_m$	see Soille (1999)
150	shortest connected path, $y$-dir., cityblock, in ${\Omega }_m$	see Soille (1999)

Equations

$$\begin{split}p(\boldsymbol{u}_{f}|\boldsymbol{\lambda}_{f},\boldsymbol{\theta}_{cf},\boldsymbol{\theta}_{c})&=\int\underbrace{p_{cf}(\boldsymbol{u}_{f}|{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{u}_{c}},\boldsymbol{\theta}_{cf})}_{\text{decoder}}\underbrace{p_{\text{CGM}}(\boldsymbol{u}_{c}|\boldsymbol{\lambda}_{c},\boldsymbol{\theta}_{c})}_{\text{CGM}}\underbrace{p_{c}(\boldsymbol{\lambda}_{c}|\boldsymbol{\lambda}_{f},\boldsymbol{\theta}_{c})}_{\text{encoder}}d\boldsymbol{\lambda}_{c}d\boldsymbol{u}_{c}.\end{split}$$

$$\begin{array}{rl}p({\mathbf{u}}_f|{\mathbf{\lambda }}_f,{\mathbf{\theta }}_{cf},{\mathbf{\theta }}_c)& =\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.\end{array}$$

$${\nabla }_{\mathbf{x}}P-\mu {\Delta }_{\mathbf{x}}\mathbf{V}$$

$${\nabla }_{\mathbf{x}}\cdot \mathbf{V}$$

$$\mathbf{V}$$

$$\mathbf{t}$$

$${\mathbf{u}}_f({\mathbf{\lambda }}_f):{\mathbf{\lambda }}_f\mapsto {\mathbf{u}}_f$$

$${\ell }_f\ll r_0\ll L$$

$${\left(\frac{r_0}{L}\right)}^2\ll 1,r_0\gtrsim 5{\ell }_f$$

$$\bar{P}=\frac{1}{|{\Omega }_{\text{f, RVE}}|}{\int }_{{\Omega }_{f,\text{RVE}}}P(\mathbf{x})dV,\bar{\mathbf{V}}=\frac{1}{|{\Omega }_{\text{f, RVE}}|}{\int }_{{\Omega }_{f,\text{RVE}}}\mathbf{V}(\mathbf{x})dV,$$

$$\bar{\mathbf{V}}$$

$${\nabla }_{\mathbf{x}}\cdot \bar{\mathbf{V}}$$

$$\bar{\mathbf{V}}\cdot \mathbf{n}$$

$$\bar{P}$$

$$\mathbf{K}(\mathbf{x},{\mathbf{\lambda }}_c)=\sum _{m=1}^{N_{\text{cells, c}}}{1}_{{\Omega }_m}(\mathbf{x}){\mathbf{K}}_m({\mathbf{\lambda }}_c),$$

$${\mathbf{K}}_m({\mathbf{\lambda }}_c)={\mathbf{K}}_m({\lambda }_{c,m})=e^{{\lambda }_{c,m}}\mathbf{I},$$

$${\mathbf{u}}_c({\mathbf{\lambda }}_c):{\mathbf{\lambda }}_c\mapsto {\mathbf{u}}_c.$$

$${\lambda }_{c,m}=\sum _{j=1}^{N_{\text{features},m}}{\tilde{\theta }}_{c,jm}{\phi }_{jm}({\mathbf{\lambda }}_f)+{\sigma }_{c,m}{Z}_m,Z_m\sim \mathcal{N}(0,1),$$

$$p_c({\mathbf{\lambda }}_c|{\mathbf{\lambda }}_f,{\mathbf{\theta }}_c)=\prod _{m=1}^{\dim ({\mathbf{\lambda }}_c)}\mathcal{N}({\lambda }_{c,m}|{\tilde{\mathbf{\theta }}}_{c,m}^T{\mathbf{\phi }}_m({\mathbf{\lambda }}_f),{\sigma }_{c,m}^2),$$

$$u_{f,i}=\sum _{j\in G^{(c)}}{W}_{ij}{u}_{c,j}({\mathbf{\lambda }}_c)+{\tau }_{cf,i}^{-1/2}{Z}_i,Z_i\sim \mathcal{N}(0,1),$$

$$p_{cf}({\mathbf{u}}_f|{\mathbf{u}}_c({\mathbf{\lambda }}_c),{\mathbf{\theta }}_{cf})=\mathcal{N}({\mathbf{u}}_f|\mathbf{W}{\mathbf{u}}_c({\mathbf{\lambda }}_c),\text{diag}[{\mathbf{\tau }}_{cf}^{-1}]),$$

$$p({\tilde{\mathbf{\theta }}}_c|\mathbf{\gamma })=\prod _{m=1}^{\dim ({\mathbf{\lambda }}_c)}\prod _{j=1}^{N_{\text{features},m}}\mathcal{N}({\tilde{\theta }}_{c,jm}|0,{\gamma }_{jm}^{-1}),$$

$$p({\gamma }_{jm})=Gamma({\gamma }_{jm}|a,b)=b^a{\gamma }_{jm}^{a-1}{e}^{-b{\gamma }_{jm}}/\Gamma (a).$$

$$p({\tau }_{c,m})$$

$$p({\tau }_{cf,i})$$

$$\begin{array}{rl}\mathcal{L}({\mathbf{\theta }}_{cf},{\mathbf{\theta }}_c)=\prod _{n=1}^N& p({\mathbf{u}}_f^{(n)}|{\mathbf{\lambda }}_f^{(n)},{\mathbf{\theta }}_{cf},{\mathbf{\theta }}_c)\\ =\prod _{n=1}^N& \int \mathcal{N}({\mathbf{u}}_f^{(n)}|\mathbf{W}{\mathbf{u}}_c({\mathbf{\lambda }}_c^{(n)}),\text{diag}[{\mathbf{\tau }}_{cf}{]}^{-1})\cdot \\ \prod _{m=1}^{\dim ({\mathbf{\lambda }}_c)}& \left[\mathcal{N}({\lambda }_{c,m}^{(n)}|{\tilde{\mathbf{\theta }}}_{c,m}^T{\mathbf{\phi }}_m({\mathbf{\lambda }}_f^{(n)}),{\tau }_{c,m}^{-1})\right]d{\mathbf{\lambda }}_c^{(n)},\end{array}$$

$$\begin{array}{rl}p({\tilde{\mathbf{\theta }}}_c,{\mathbf{\tau }}_c,{\mathbf{\tau }}_{cf}|\mathcal{D})\propto \prod _{n=1}^N& [\int \mathcal{N}({\mathbf{u}}_f^{(n)}|\mathbf{W}{\mathbf{u}}_c({\mathbf{\lambda }}_c^{(n)}),\text{diag}[{\mathbf{\tau }}_{cf}{]}^{-1})\\ \prod _{m=1}^{\dim ({\mathbf{\lambda }}_c)}& \left[\mathcal{N}({\lambda }_{c,m}^{(n)}|{\tilde{\mathbf{\theta }}}_{c,m}^T{\mathbf{\phi }}_m({\mathbf{\lambda }}_f^{(n)}),{\tau }_{c,m}^{-1})\right]d{\mathbf{\lambda }}_c^{(n)}]\\ & \cdot \int p({\tilde{\mathbf{\theta }}}_c|\mathbf{\gamma })p(\mathbf{\gamma })d\mathbf{\gamma }\cdot p({\mathbf{\tau }}_c)p({\mathbf{\tau }}_{cf})\end{array}$$

$$\begin{gathered}p(\boldsymbol{\theta}|\mathcal{D})=\cfrac{p(\boldsymbol{\theta},\mathcal{D})}{p(\mathcal{D})}=\frac{1}{p(\mathcal{D})}\prod_{n=1}^{N}\Bigg{[}\mathcal{N}(\boldsymbol{u}_{f}^{(n)}|\boldsymbol{W}\boldsymbol{u}_{c}(\boldsymbol{\lambda}_{c}^{(n)}),\text{diag}[\boldsymbol{\tau}_{cf}]^{-1})\cdot\\ \cdot\prod_{m=1}^{\dim(\boldsymbol{\lambda}_{c})}\left[\mathcal{N}(\lambda_{c,m}^{(n)}|\tilde{\boldsymbol{\theta}}_{c,m}^{T}\boldsymbol{\varphi}_{m}(\boldsymbol{\lambda}_{f}^{(n)}),\tau_{c,m}^{-1})\right]\Bigg{]}~{}p(\tilde{\boldsymbol{\theta}}_{c}|\boldsymbol{\gamma})p(\boldsymbol{\gamma})p(\boldsymbol{\tau}_{c})p(\boldsymbol{\tau}_{cf})\end{gathered}$$

$$\text{KL}(Q(\mathbf{\theta })||p(\mathbf{\theta }|\mathcal{D}))=-\int Q(\mathbf{\theta })\log \frac{p(\mathbf{\theta }|\mathcal{D})}{Q(\mathbf{\theta })}d\mathbf{\theta }.$$

$$\log p(\mathcal{D})=\log \int p(\mathbf{\theta },\mathcal{D})d\mathbf{\theta }=\mathcal{F}(Q)+\text{KL}(Q(\mathbf{\theta })||p(\mathbf{\theta }|\mathcal{D})),$$

$$\mathcal{F}(Q)=\int Q(\mathbf{\theta })\log \frac{p(\mathbf{\theta },\mathcal{D})}{Q(\mathbf{\theta })}d\mathbf{\theta }.$$

$$\begin{array}{rl}& Q({\left\{{\mathbf{\lambda }}_c^{(n)}\right\}}_{n=1}^N,{\tilde{\mathbf{\theta }}}_c,{\mathbf{\tau }}_c,{\mathbf{\tau }}_{cf},\mathbf{\gamma })=\\ & =\prod _{m=1}^{\text{dim}({\mathbf{\lambda }}_c)}\left[Q_{{\tau }_{c,m}}({\tau }_{c,m})\prod _{j=1}^{N_{\text{features},m}}\left[Q_{{\tilde{\theta }}_{c,jm}}({\tilde{\theta }}_{c,jm})Q_{{\gamma }_{jm}}({\gamma }_{jm})\right]\right]\\ & \prod _{i=1}^{\text{dim}({\mathbf{u}}_f)}\left[Q_{{\tau }_{cf,i}}({\tau }_{cf,i})\right]\prod _{n=1}^N{Q}_{{\mathbf{\lambda }}_c^{(n)}}({\mathbf{\lambda }}_c^{(n)})\\ & \approx p({\left\{{\mathbf{\lambda }}_c^{(n)}\right\}}_{n=1}^N,{\tilde{\mathbf{\theta }}}_c,{\mathbf{\tau }}_c,{\mathbf{\tau }}_{cf},\mathbf{\gamma }|\mathcal{D}).\end{array}$$

Full text

††

Source code as well as examples can be found at: github.com/congriUQ/physics_aware_surrogate

Keywords: Uncertainty Quantification, Bayesian Inference, Coarse-Graining, Variational Inference, Random Media, Multi-Physics Models

1 Introduction

The present paper is concerned with the data-driven construction of coarse-grained descriptions and predictive surrogate models in the context of random heterogeneous materials. This is of paramount importance for carrying out various uncertainty quantification (UQ) tasks, both in terms of forward and backward uncertainty propagation. More generally, the development of efficient surrogates is instrumental in enabling many-query applications in the presence of complex models of physical and engineering systems. The significant progress that has been materialized in the fields of statistical or machine learning Ghahramani (2015) in the last few years, particularly in the context of Deep Learning LeCun et al. (2015); Goodfellow et al. (2016) has found its way in physical systems where training data for the surrogates are generated by forward simulations of the reference, high-fidelity model. We note however that such problems exhibit (at least) two critical differences as compared to typical machine learning applications. Firstly, the dimension of the (random) inputs and outputs encountered in complex physical models is generally much larger than speech- or image-related data and associated tasks. Furthermore, the dependencies between these two sets of variables cannot be discovered with brute-force dimensionality reduction techniques. Secondly and most importantly, these are not Big Data problems. The abundance of data that can be rightfully assumed in many data-driven contexts, is contrary to the nature of the problem we are trying to address, that is learning surrogate models, or more generally enabling UQ tasks, with the least possible high-fidelity output data. Such regimes could be more appropriately be referred to as Small Data or better yet Tall Data in consideration of the dimensionality above Koutsourelakis et al. (2016).

The development of surrogates for UQ purposes in the context of continuum thermodynamics where pertinent models are based on PDEs and ODEs has a long history and some of the most well-studied methods have been based on (generalized) Polynomial Chaos expansions (gPC) Wiener (1938), Ghanem and Spanos (1991), Xiu and Karniadakis (2002) which have become popular in the last decade also due to the emergence of non-intrusive, sparse-grid stochastic collocation approaches Xiu and Hesthaven (2005); Ma and Zabaras (2009); Lin and Tartakovsky (2009). Despite recent progress in connecting gPC with domain decomposition methods Lin et al. (2010); Tipireddy et al. (2017, 2018), these approaches typically struggle in problems with high-dimensional stochastic input, as is the case with random heterogeneous materials Torquato and Lu (1993) encountered in problems such as flow through porous media Bilionis et al. (2013); Atkinson and Zabaras (2019); Mo et al. (2018) or electromagnetic wave propagation through random permittivity fields Chauvière et al. (2006).

Another well-established strategy makes use of statistical couplings between the reference model and less expensive, lower-fidelity models in order to enable *multi-fidelity * methods O’Hagan and Kennedy (2000). In the context of PDEs, such lower-fidelity or reduced-order models may readily be constructed by utilizing potentially much coarser spatio-temporal discretizations which, in conjunction with appropriate statistical models, are capable of providing fully probabilistic predictions Koutsourelakis (2009); Perdikaris et al. (2017). The most striking benefit of such methods is that they automatically retain some of the underlying physics of the high-fidelity model and are therefore able to produce accurate predictions even when the data is scarce and the stochastic input dimensionality is high. Another category of algorithms falling under the class of reduced-order models is given by reduced-basis (RB) methods Hesthaven et al. (2016); Quarteroni et al. (2016) where, based on a small set of ‘snapshots’ i.e. forward model evaluations, the solution space dimensionality is reduced by projection onto the principal ‘snapshot’ directions. The reduced-order solution coefficients are then found either by standard Galerkin projection Rowley et al. (2004); Cui et al. (2015); Afkham and Hesthaven (2017) or regression-based strategies Yu and Hesthaven (2018); Guo and Hesthaven (2018); Hesthaven and Ubbiali (2018). Apart from issues of efficiency and stability, both gPC and RB approaches in their standard form are generally treated in a non-Bayesian way and therefore only yield point estimates instead of full predictive posterior distributions.

A more recent trend is to view surrogate modeling as a supervised learning problem that can be handled by pertinent statistical learning tools, e.g. Gaussian Process (GP) regression Rasmussen and Williams (2006); Bilionis and Zabaras (2017) which provide closed-form predictive distributions by default. Although recent advances have been made towards multi-fidelity data fusion O’Hagan and Kennedy (2000); Raissi et al. (2017); Perdikaris et al. (2015) and incorporation of background physics Raissi and Karniadakis (2018); Yang et al. (2018); Lee et al. (2018); Tipireddy and Tartakovsky (2018) via Gaussian Processes, the poor scaling behavior with stochastic input dimension remains one of the main challenges for GP models.

More recently, deep neural networks (DNNs) LeCun et al. (2015); Goodfellow et al. (2016) have found their way into surrogate modeling of complex computer codes Tripathy and Bilionis (2018); Zhu and Zabaras (2018); Mo et al. (2018); Yang and Perdikaris (2019). One of the most promising developments in the integration of such tools in physical modeling are physics-informed neural networks Raissi et al. (2017, 2019); Tartakovsky et al. (2018) which are trained by minimization of a loss function augmented by the residuals of the governing equations Raissi et al. (2017) and Raissi et al. (2017). In Zhu et al. (2019), this methodology is generalized and extended to cases without any labeled data, i.e. without computing the output of the reference, high-fidelity model. However, such models transfer much of the numerical burden from the generation of training data to the repeated evaluation of the residual during model training. While promising results have been obtained, we note that several challenges remain with regards to the incorporation of prior physical knowledge in the form of symmetries/invariances etc. into DNNs as well as their capabilities under extrapolative conditions.

In this paper, we present a novel, fully Bayesian strategy for surrogate modeling designed for problems characterized by

• •

high input dimension, i.e. the “curse of dimensionality” – in the problems examined, we consider input data (random material microstructures) of effective dimensionality $\gtrsim 10^{4}$; and

• •

small number of training data $N$, i.e. forward or fine grained model (FGM) input/output pairs that can practicably be generated in order to train any sort of reduced-order model or surrogate – in the examples presented here, we use $N\lesssim 100$ training samples.

In the context of random media and in order to achieve the aforementioned objectives, it is of primary importance (a) to extract only the microstructural features pertinent for reconstruction of the fine-grained model (FGM) output and (b) to construct a surrogate that retains as much as possible the underlying physics of the problem. Whilst point (a) should effectively reduce the input dimension, (b) may exclude a myriad of potential model outputs violating essential physical principles of the FGM. Consequently, the core unit of the surrogate model established here is based on a coarse-grained model (CGM) operating on simplified physics valid on a larger scale than the characteristic length of the microstructural features in the FGM. The CGM serves as a stencil that automatically retains the primary physical characteristics of the FGM. It is parametrized by a latent, low-dimensional set of variables that locally encode the effective material properties of the underlying random heterogeneous microstructure and are found by appropriate sparsity-inducing statistical learning Schöberl et al. (2017). Finally, a decoder distribution maps the CGM outputs to the FGM solution space, yielding fully probabilistic model predictions. Instead of viewing the encoding step (dimension reduction) separately from the surrogate modeling task Ma and Zabaras (2011); Xing et al. (2016), both are integrated and trained simultaneously, ensuring that the encoded, latent representation of the microstructure is maximally predictive of the sought FGM response instead of the original fine-scale microstructure itself. Model training is carried out using modern Stochastic Variational Inference (SVI) techniques Paisley et al. (2012); Hoffman et al. (2013) in combination with efficient stochastic optimization algorithms Kingma and Ba (2014). A hierarchical Bayesian treatment of the model parameters obviates any need for fine-tuning or additional input by the analyst Bishop and Tipping (2000).

The rest of the paper is organized as follows. Section 2 lays the methodological foundations for the proposed model, including a presentation of the pertinent physical structure, associated state variables and governing equations. Moreover, a detailed definition of both encoding and decoding steps is given, as well as a clear specification of the prior model and the training/prediction algorithms. In section 3, we provide numerical evidence that the proposed model can lead to accurate predictions even when trained on a few tens of FGM data. Furthermore, we present a methodological framework for adaptively adjusting the CGM complexity in order to enhance its predictive ability. Section 4 summarizes the main findings and presents further extensions and possible applications of the proposed model.

2 Methodology

We present a general, probabilistic framework to coarse-graining which we specialize in the subsequent sections. The starting point is an expensive, deterministic, fine-grained model (FGM) with high-dimensional, random inputs $\boldsymbol{\lambda}_{f}$ and high-dimensional outputs $\boldsymbol{u}_{f}=\boldsymbol{u}_{f}(\boldsymbol{\lambda}_{f})$. Due to the high dimensionality $\dim(\boldsymbol{\lambda}_{f})$ of the input uncertainties $\boldsymbol{\lambda}_{f}$ and the small number of forward evaluations that can practically be performed, any naive regression approach, i.e. learning a direct map from $\boldsymbol{\lambda}_{f}$ to $\boldsymbol{u}_{f}$ is doomed to fail.

We seek instead a lower-dimensional, effective representation $\boldsymbol{\lambda}_{c}$ of the high-dimensional stochastic inputs $\boldsymbol{\lambda}_{f}$ which encodes as much as possible information on the respective FGM output $\boldsymbol{u}_{f}$, see Figure 1. We note that this differs from classical dimension reduction approaches in the sense that the encoded representation $\boldsymbol{\lambda}_{c}$ should be maximally predictive of the FGM output $\boldsymbol{u}_{f}(\boldsymbol{\lambda}_{f})$ – and not for the reconstruction of $\boldsymbol{\lambda}_{f}$ itself Tishby et al. (1999). Rather than assigning an abstract, statistical meaning to $\boldsymbol{\lambda}_{c}$, we associate them with a less-expensive, physical, coarse-grained model (CGM). In particular, the latent variables $\boldsymbol{\lambda}_{c}$ represent the inputs of the CGM. We denote by $\boldsymbol{u}_{c}=\boldsymbol{u}_{c}(\boldsymbol{\lambda}_{c})$ the corresponding CGM outputs which we subsequently attempt to link with the FGM response $\boldsymbol{u}_{f}(\boldsymbol{\lambda}_{f})$, see Figure 1. The CGM effectively represents the bottleneck through which information that the FGM input $\boldsymbol{\lambda}_{f}$ provides about the FGM output $\boldsymbol{u}_{f}$ is squeezed. While several possibilities exist, independently of the final choice for the CGM, information loss will generally take place during the coarse-graining process, which in turn results in predictive uncertainty that we aim at quantifying. Moreover, tools for comparing different CGMs as well as strategies of refining the CGM adaptively in order to improve predictive accuracy are developed in a later section.

The general framework described above can be effected by a model that integrates the following three key components (see Figure 1 and Grigo and Koutsourelakis (2017, 2019)):

• •

Encoding - Dimension reduction: A mapping from the high - dimensional microstructural description $\boldsymbol{\lambda}_{f}$ to a much lower-dimensional quantity $\boldsymbol{\lambda}_{c}$ retaining as much as possible information on the corresponding response $\boldsymbol{u}_{f}(\boldsymbol{\lambda}_{f})$ – this mapping is carried out probabilistically by the conditional distribution $p_{c}(\boldsymbol{\lambda}_{c}|\boldsymbol{\lambda}_{f},\boldsymbol{\theta}_{c})$ parametrized by $\boldsymbol{\theta}_{c}$.

• •

Model reduction: A skeleton of a coarse-grained model (CGM) that employs the dimension-reduced quantities $\boldsymbol{\lambda}_{c}$ as input and yields an output $\boldsymbol{u}_{c}$. Such a model can be based on simplified physics and/or coarser spatial discretization. In our case, this is constructed by employing a coarse discretization of the Darcy-flow governing equations as it will be explained in detail later. We denote by $p_{\text{CGM}}(\boldsymbol{u}_{c}|\boldsymbol{\lambda}_{c})$ the generally stochastic input-output map implied.

• •

Decoding - Reconstruction: A probabilistic mapping from the CGM output $\boldsymbol{u}_{c}$ back onto the original FGM output $\boldsymbol{u}_{f}$ – this mapping is mediated by the distribution $p_{cf}(\boldsymbol{u}_{f}|\boldsymbol{u}_{c},\boldsymbol{\theta}_{cf})$ parametrized by $\boldsymbol{\theta}_{cf}$.

The combination of the aforementioned three densities yields

[TABLE]

Assuming a deterministic CGM, i.e. $p_{\text{CGM}}(\boldsymbol{u}_{c}|\boldsymbol{\lambda}_{c})=\delta(\boldsymbol{u}_{c}-\boldsymbol{u}_{c}(\boldsymbol{\lambda}_{c}))$, this simplifies to

[TABLE]

The resulting parametrized density can be viewed as the likelihood of a FGM input/output pair $\{\boldsymbol{\lambda}_{f},\boldsymbol{u}_{f}\}$ and is thus the key quantity in the proposed probabilistic surrogate. We note that the three constitutive densities are highly modular and can be adapted to the particulars of the FGM and the CGM adopted. In the following, the FGM consists of a Stokes-flow simulator through a random heterogeneous medium, see section 2.1. The CGM is described by Darcy’s equations of flow through a permeable medium and is explained in section 2.2. We present the essential components of the proposed model as well as the specialization of the variables and densities above in section 2.3. In section 2.4, we present the learning engine that is capable of training the model with small data and in section 2.5, we discuss how the trained model can be used to produce probabilistic predictions. We conclude this methodological section with a discussion on the numerical complexity (section 2.6) as well as a novel algorithm for the adaptive refinement of the coarse-grained model (section 2.7).

2.1 The fine-grained model: Stokes flow through random porous media

Flow of incompressible Newtonian fluids through random porous media is characterized by low fluid velocities and is therefore dominated by viscous forces. In such regimes, the steady-state Navier-Stokes equations reduce to Stokes equations of momentum and mass conservation,

[TABLE]

where $P$ and $\boldsymbol{V}$ are pressure and velocity fields, $\boldsymbol{f}$ is an external force field, $\Omega_{f}=\Omega\backslash\Omega_{s}$ the pore part of the domain $\Omega=\Omega_{f}\cup\Omega_{s}$, where $\Omega_{s}$ denotes the part of impermeable solid material, and $\mu$ the fluid viscosity which we may set to $\mu=1$ for convenience.

On the solid-fluid interface $\Gamma_{\text{int}}$, the no-slip boundary condition

[TABLE]

is imposed, whereas on the boundary $\Gamma_{\text{ext}}=\partial\Omega=\partial\Omega_{f}\backslash\Gamma_{\text{int}}=\Gamma_{P}\cup\Gamma_{\boldsymbol{V}}$ of the macroscopic system, we apply

[TABLE]

where $\boldsymbol{t}$ is the Cauchy traction and $\boldsymbol{n}$ the unit outward normal.

For the remainder of the paper, we consider a unit square domain $\Omega=[0,1]^{2}\subset\mathbb{R}^{2}$ with randomly distributed non-overlapping polydisperse spherical exclusions, see Figure 2 for a schematic representation. Such microstructures exhibit topologically fully connected pore spaces $\Omega_{f}$ and lead to unique solutions of (3) for the case where $\Gamma_{\text{ext}}=\partial\Omega$.

For the numerical solution of Equation (3), we employ triangular meshes with standard Taylor-Hood elements (i.e. quadratic and linear shape functions for velocity and pressure, respectively), implemented using the FEniCS finite element software package Logg et al. (2012). For generic two-phase media, $\boldsymbol{\lambda}_{f}$ consists of a list of binary variables representing the phase of each pixel/voxel. For the particular case considered of randomly distributed non-overlapping polydisperse spherical exclusions, one could also represent the microstructure with a list of exclusion center coordinates and the corresponding radii. We denote by $\boldsymbol{u}_{f}$ the discretized solution vector which can consist of pressures or velocities and write

[TABLE]

to indicate the map from the microstructure $\boldsymbol{\lambda}_{f}$ to the PDE response of interest.

Computational considerations

Typically, the structure of the porous material is very complex with a characteristic length-scale $\ell$ much smaller than the macroscopic length scale $L=1$ of the problem domain $\Omega$. This large difference in scales means that (a) from a modeling point of view, a probabilistic description would be more realistic for the porous microstructure, and (b) the numerical solution of (3) requires discretizations fine enough to resolve all microstructural details. Moreover, the uncertainty of microstructures $\boldsymbol{\lambda}_{f}\sim p(\boldsymbol{\lambda}_{f})$ requires repeated evaluation of the numerical solver (4) to obtain accurate estimates (e.g. with Monte Carlo) of any statistic of interest relating to the output which quickly becomes impracticable due to the significant cost associated with each FGM run.

2.2 The coarse-grained model: Darcy equations for diffusion through a permeable medium

The macroscopic Darcy-like behavior of Stokes flow through random porous media (3) is a long-standing result from homogenization theory first shown for periodic microstructures Sanchez-Palencia (1980); Tartar (1980), later generalized for connected solid phase matrix Allaire (1989), the non-stationary process Allaire (1992) and non-periodic media Whitaker (1986). In classical homogenization, the crucial step is to find a local boundary value problem on a Representative Volume Element (RVE) Sandström et al. (2013); Marchenko and Khruslov (2006); Tartar (2010). The RVE size $r_{0}$, i.e. the radius of the averaging volume, must be large compared to the microstructural length scales $\ell_{f},\ell_{s}$. On the other hand, the volume averaged quantities, i.e. the intrinsic phase averages $\bar{P}$ and $\bar{\boldsymbol{V}}$, should show sufficient variation over the problem domain $\Omega$, i.e. $L$ should be large compared to $r_{0}$. To summarize, Stokes to Darcy convergence is observed in the limit of

[TABLE]

which, according to Whitaker (1986), may slightly be relaxed to the conditions

[TABLE]

to safely assume Darcy-type flow. Using the definition for the intrinsic phase averages of pressure and velocity fields $P$ and $\boldsymbol{V}$

[TABLE]

where $|\Omega_{f,\text{RVE}}|$ is the volume of the fluid part of the RVE, the effective Darcy constitutive behavior can be written as

[TABLE]

with the unit outward normal vector $\boldsymbol{n}$, the viscosity $\mu=1$ and some unknown permeability tensor $\boldsymbol{K}$. Determining/estimating the effective permeability field $\boldsymbol{K}$ remains a substantial computational problem which in homogenization theory is usually approached by solving an RVE subscale problem Sandström et al. (2013). We emphasize that full Stokes/Darcy equivalence can only be assumed in the limit of infinite scale separation as defined by Equation (5). However, this does not rule out the possibility to use the Darcy equations as a stencil of a machine learning model applied to fine-grained Stokes flow data even far away from the scale separation limit, as we will see in the experiments section 3.

2.2.1 The Darcy permeability tensor field $\boldsymbol{K}=\boldsymbol{K}(\boldsymbol{x},\boldsymbol{\lambda}_{c})$

The material parameters of the CGM which provide the necessary closure pertain to the permeability tensor field $\boldsymbol{K}=\boldsymbol{K}(\boldsymbol{x},\boldsymbol{\lambda}_{c})$. We subdivide the problem domain $\Omega$ into $N_{\text{cells},c}=\text{dim}(\boldsymbol{\lambda}_{c})$ of non-overlapping subregions/ cells $\left\{\Omega_{m}\right\}_{m=1}^{N_{\text{cells, c}}}$ such that $\cup_{m=1}^{N_{\text{cells, c}}}\Omega_{m}=\Omega$ and $\Omega_{m_{1}}\cap\Omega_{m_{2}}=\emptyset$ for $m_{1}\neq m_{2}$. In most cases examined in the experimental section 3, the subregions $\Omega_{m}$ are squares of equal size, as can be seen in the bottom row of Figure 11. Within each subregion $\Omega_{m}$, a constant permeability tensor $\boldsymbol{K}$ is assumed, i.e.

[TABLE]

where $\mathbb{1}_{\Omega_{m}}(\boldsymbol{x})$ is the indicator function of subregion $\Omega_{m}$. The local permeability tensors $\boldsymbol{K}_{m}(\boldsymbol{\lambda}_{c})$ are positive definite matrices which, in the simplest case based on the isotropy assumption, take the form

[TABLE]

where $\boldsymbol{I}$ is the identity matrix. We note that the number of cells $N_{\text{cells},c}$ does not need to coincide with the finite elements employed for the solution of the CGM. The former control the dimension of the reduced representation whereas the latter determine the computational cost in the solution of the CGM. Naturally, different discretizations or representations of the coarse model’s spatially-varying properties can be adopted, e.g. $\boldsymbol{\lambda}_{c}$ could relate to basis functions’ coefficients in an appropriate expansion of $\boldsymbol{K}(\boldsymbol{x},\boldsymbol{\lambda}_{c})$.

Analogously to the FGM in Equation (4), the Darcy-based, coarse-grained model (CGM) is defined as the deterministic mapping

[TABLE]

Remarks

• •

The low-dimensional latent variables $\boldsymbol{\lambda}_{c}$ encode the high dimensional descriptors $\boldsymbol{\lambda}_{f}$ of the fine scale porous domain of the Stokes flow FGM as an effective, Darcy-type permeability random field $\boldsymbol{K}=\boldsymbol{K}(\boldsymbol{x},\boldsymbol{\lambda}_{c})$. We re-emphasize that it is extraneous if the latent representation $\boldsymbol{\lambda}_{c}$ is a faithful decoding of its fine scale analogue $\boldsymbol{\lambda}_{f}$ for the purpose of accurately reconstructing $\boldsymbol{\lambda}_{f}$ directly. Instead, $\boldsymbol{\lambda}_{c}$ must be maximally predictive (when solving the inexpensive CGM $\boldsymbol{u}_{c}(\boldsymbol{\lambda}_{c})$) of the reconstruction of the FGM output $\boldsymbol{u}_{f}$. Hence, $\boldsymbol{\lambda}_{c}$ may be completely different from the reduced coordinates identified by typical dimensionality reduction techniques applied directly on $\boldsymbol{\lambda}_{f}$. The encoder density $p_{c}(\boldsymbol{\lambda}_{c}|\boldsymbol{\lambda}_{f},\boldsymbol{\theta}_{c})$ should therefore be trained in such a way that it extracts all microstructural features from $\boldsymbol{\lambda}_{f}$ that are relevant in predicting the FGM output $\boldsymbol{u}_{f}$.

• •

The decoder density p_{cf}(\boldsymbol{u}_{f}|{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{u}_{c}(\boldsymbol{\lambda}_{c})},\boldsymbol{\theta}_{cf}) attempts to reconstruct the FGM Stokes flow response $\boldsymbol{u}_{f}$ given the solution $\boldsymbol{u}_{c}(\boldsymbol{\lambda}_{c})$ of the Darcy-based CGM Equation (8). While a general parametrization can be employed (e.g. making use of (deep) neural networks), this would necessitate either huge amounts of data or very strong regularization assumptions. We exploit instead the spatial character of the problem, i.e. if the FGM and CGM outputs $u_{f,i},u_{c,j}$ are directly related to spatial locations $\boldsymbol{x}_{i},\boldsymbol{x}_{j}$ (e.g. if they represent pressure or velocity responses at these points), it is natural to assume that components $u_{c,j}$ corresponding to points $\boldsymbol{x}_{j}$ in the vicinity of $\boldsymbol{x}_{i}$ are most relevant for reconstruction of $u_{f,i}$ corresponding to point $\boldsymbol{x}_{i}$.

• •

Independently of the values assigned to the parameters $\boldsymbol{\theta}_{c},\boldsymbol{\theta}_{cf}$, predictions produced by the surrogate will be probabilistic as reflected by the density in Equation (2). This can be attributed to the epistemic uncertainty due to the unavoidable information loss during the coarse-graining process as there is an upper bound on the mutual information $I(\boldsymbol{\lambda}_{c},\boldsymbol{u}_{f})<I(\boldsymbol{\lambda}_{f},\boldsymbol{u}_{f})$ when $\dim(\boldsymbol{\lambda}_{c})\ll\dim(\boldsymbol{\lambda}_{f})$ and no redundancies in $\boldsymbol{\lambda}_{f}$ Tishby et al. (1999).

The following subsections are devoted to giving a clear description of the main building blocks $p_{c}$ and $p_{cf}$ as well as of the prior model we apply on the model parameters $\boldsymbol{\theta}_{c}$ and $\boldsymbol{\theta}_{cf}$.

2.3 Components of the proposed model

2.3.1 The encoder distribution $p_{c}$

The encoder distribution $p_{c}$ that effects the mapping from the high-dimensional $\boldsymbol{\lambda}_{f}$ to the lower-dimensional $\boldsymbol{\lambda}_{c}$ can assume various forms. Any successful strategy should be capable of extracting appropriate features from $\boldsymbol{\lambda}_{f}$ or, more generally, learning the right representation for predicting $\boldsymbol{\lambda}_{c}$. An obvious choice would be to exploit the expressive power of (deep) Neural Nets (NNs) and, given the spatial character of $\boldsymbol{\lambda}_{f}$, Convolutional Neural Networks (CNNs) in particular Fukushima (1980); Lecun et al. (1998); Krizhevsky et al. (2012). A common misconception is that Deep Learning techniques are devoid of any feature engineering. While deep NNs are highly expressive, contain a lot of learnable parameters and reduce the need for feature engineering, it is not true that data pre-processing or feature extraction are totally irrelevant333https://github.com/tensorflow/transform. Several components, such as the number, size and type of layers or e.g. the type and number of pooling layers in CNNs are in fact hard-coded, feature-extraction mechanisms. More importantly, the basic premise in successful applications of these models is the availability of large amounts of data, from which the necessary information can be extracted. As mentioned earlier however, we operate in the Small Data regime.

Based on this, we make use of a large set of physically-motivated, feature functions $\varphi_{jm}(\boldsymbol{\lambda}_{f})$ (discussed in detail in the sequel) in the context of the following model:

[TABLE]

In principle, the type and number of features can differ with every latent space component $m$. These feature functions are linearly combined with coefficients $\tilde{\boldsymbol{\theta}}_{c,m}=\left\{\tilde{\theta}_{c,jm}\right\}_{j=1}^{N_{\text{features},m}}$ (potentially different for every component $m$ as well) and augmented with residual Gaussian white noise of variance $\sigma_{c,m}^{2}$, so that the set of free parameters for the distribution $p_{c}$ is $\boldsymbol{\theta}_{c}=\left\{\tilde{\boldsymbol{\theta}}_{c,m},\sigma_{c,m}^{2}\right\}_{m=1}^{\dim(\boldsymbol{\lambda}_{c})}$ 444We use the notation $\boldsymbol{\Sigma}_{c}=\text{diag}[\boldsymbol{\sigma}_{c}^{2}]$, where $\boldsymbol{\sigma}_{c}^{2}=\left\{\sigma_{c,m}^{2}\right\}_{m=1}^{\dim{\boldsymbol{\lambda}_{c}}}$ is the vector of variances, wherever it is more convenient..

The distribution for $p_{c}$ is thus

[TABLE]

with $\boldsymbol{\varphi}_{m}(\boldsymbol{\lambda}_{f})=\left\{\varphi_{jm}(\boldsymbol{\lambda}_{f})\right\}_{j=1}^{N_{\text{features},m}}$ being the vector of feature functions used to predict the latent space component $\boldsymbol{\lambda}_{c,m}$.

Feature functions:

The most critical aspect in the expressivity of the encoder pertains to the feature functions $\varphi_{m}(\boldsymbol{\lambda}_{f})$ employed. In this work, we make use of morphological measures Torquato (2001); Lu and Torquato (1992); Torquato and Lu (1993); Lowell et al. (2006), quantities from classical fluid dynamics Sutera and Skalak (1993); Kozeny (1927); Carman (1937) or other physics Archie (1947), image analysis Soille (1999) and low-dimensional autoencoder representations Doersch (2016). Some important classes of possible feature functions are depicted in Figure 3 and briefly explained in the following list:

• •

Pore fraction: The fraction of pore space $|\Omega_{f}|/|\Omega|$ over a certain region and several nonlinear functions thereof (e.g. $\exp,\log,$ certain powers);

• •

Interface area: The interface area between matrix- and exclusion phase over a certain region and nonlinear functions thereof;

• •

Effective medium approximations: There exist several approximation formula for effective material properties in diffusion problems (Bruggeman formula Bruggeman (1935), Maxwell-Garnett approximation Landauer (1978)) which can be generalized to the infinite contrast limit of matrix/exclusion phase permeabilities;

• •

Chord length density: The density of lengths of segments of an infinitely long line being in either matrix or exclusion phase, see blue line with red dots in Figure 3. A ‘segment’ is a part of the line between two neighboring dots – the feature function is given by evaluating the chord length density for several lengths $d$, see 3 for the length scales used in the experiments;

• •

Statistics of exclusion radii: Mean, variance, max, min, etc. of exclusion radii;

• •

Statistics of mutual exclusion distances: Mean, variance, max, min, etc. of mutual exclusion distances (center-to-center or edge-to-edge). See the orange lines in Figure 3 for illustration;

• •

Pore-size density: Density of the distance from a random point in $\Omega$ to the nearest matrix/exclusion interface. For illustration see the red/gray circles in Figure 3. Nearest-neighbor functions Torquato (2001) are closely related. Again, the pore-size density is to be evaluated for several lengths $d$, see 3;

• •

Lineal-path function: The probability of a random line of length $d$ being purely in matrix (or exclusion) phase. See the green lines in Figure 3 for illustration and 3 for lengths $d$ used;

• •

2-point correlations: The probability, that two random points of distance $d$ lie both in matrix (exclusion) phase – see 3 for lengths $d$;

• •

Distance transforms: For any point in $\Omega$, compute the distance to the closest exclusion (matrix phase), see right part of Figure 3. Features can be constructed from statistics of the distance transform.

Ultimately, several hundreds of such functions can be deployed – in Section 3, we use the 150 listed in Table 3. In order to avoid overfitting in the Small Data regime we are operating (i.e. $N\lesssim 100$), we make use of a fully Bayesian, sparsity-inducing prior model Bishop and Tipping (2000) capable of automatically controlling model complexity by pruning all irrelevant features and identifying the most salient ones, see section 2.3.3.

2.3.2 The decoder distribution $p_{cf}$

The probabilistic coarse-to-fine mapping $p_{cf}(\boldsymbol{u}_{f}|\boldsymbol{u}_{c},\boldsymbol{\theta}_{cf})$ should take into account the spatial characteristics of the problem, e.g. that FGM and CGM outputs $u_{f,i},u_{c,j}$ may be given by physical quantities evaluated at spatial locations $\boldsymbol{x}_{i},\boldsymbol{x}_{j}$. For simplicity, let us assume that $u_{f,i}=P(\boldsymbol{x}_{i})$ is the pressure field solution of the FGM Equation (3) interpolated at a regular fine scale grid $G^{(f)}$, $\boldsymbol{x}_{i}\in G^{(f)}$, and $u_{c,j}=\bar{P}(\boldsymbol{x}_{j})$ its Darcy counterpart from Equation (8) evaluated at a much coarser grid $G^{(c)}$, $\boldsymbol{x}_{j}\in G^{(c)}$. Then the model for $p_{cf}$ could look like

[TABLE]

where $\boldsymbol{W}$ should have the structure of an interpolation matrix from $G^{(c)}$ to $G^{(f)}$ and $\text{diag}[\boldsymbol{\tau}_{cf}]$ is a diagonal precision matrix to be learned from the data. The decoder distribution $p_{cf}$ can thus be written as

[TABLE]

with model parameters $\boldsymbol{\theta}_{cf}=\left\{\boldsymbol{W},\boldsymbol{\tau}_{cf}\right\}$. We emphasize that many different possibilities exist in terms of the parametrization and architecture of $p_{cf}$. One should note here that, in contrast to $p_{c}$ which involves a map from the high-dimensional $\boldsymbol{\lambda}_{f}$ to a lower-dimensional $\boldsymbol{\lambda}_{c}$, $p_{cf}$ encompasses a map from the lower-dimensional $\boldsymbol{u}_{c}$ to the high-dimensional $\boldsymbol{u}_{f}$ which is obviously much more manageable. As with $p_{c}$, the option of NNs comes to mind and given the spatial character of both $\boldsymbol{u}_{c}$ and $\boldsymbol{u}_{f}$ CNNs could be introduced. We have again decided against this option, due to the Small Data setting of the problem. Given the lower dimension of the inputs $\boldsymbol{u}_{c}$, nonlinear models employing polynomials, kernels etc could be employed. An interesting alternative involves Bayesian non-parametric extensions. For instance, a Gaussian Process regression could be deployed for the probabilistic mapping $u_{c,j}=\bar{P}(\boldsymbol{x}_{j})\mapsto u_{f,i}=P(\boldsymbol{x}_{i})$, exploiting the spatial structure of the problem by making use of a stationary covariance kernel of the form $\text{cov}(u_{f,i},u_{c,j})=\text{cov}(P(\boldsymbol{x}_{i}),\bar{P}(\boldsymbol{x}_{j}))=\text{cov}(\boldsymbol{x}_{i},\boldsymbol{x}_{j})=\text{cov}(|\boldsymbol{x}_{i}-\boldsymbol{x}_{j}|)$. Alternatively, each $u_{f,i}$ can be expressed using an independent Gaussian Process in a fashion analogous to GP-LVMs Lawrence (2005); Titsias and Lawrence (2010). Such an extension is not pursued here due to the complexities associated with the inference and prediction steps Damianou et al. (2015).

2.3.3 The prior model

We advocate a fully Bayesian formulation where the (approximate) posterior densities on the model parameters $\boldsymbol{\theta}_{c},\boldsymbol{\theta}_{cf}$ may be computed. In view of the Small Data setting, we employ priors that can safeguard against overfitting, that do not require any fine-tuning or have ad hoc hyperparameters, but also promote interpretability of the results. To that end, in order to induce sparsity in the feature functions of $p_{c}$ given in Equation (13), we make use of the Automatic Relevance Determination (ARD) model Faul (2003) corresponding to

[TABLE]

where the precision hyperparameters $\gamma_{jm}$ are equipped with a conjugate Gamma hyperprior

[TABLE]

This hierarchical model has been widely studied and the capability of pruning irrelevant parameters $\theta_{c,jm}$ has been shown when small values for the hyperparameters are employed555In our investigations, we set $a=b=c=d=e=f=10^{-10}$ to avoid division by zero without affecting the converged values. Bishop and Tipping (2000); Tipping (2001).

Similarly, we apply conjugate, uninformative Gamma hyperpriors666See footnote 5 on the precision parameters $\boldsymbol{\tau}_{c},\boldsymbol{\tau}_{cf}$ in the encoder/decoder distributions $p_{c}$ and $p_{cf}$, i.e.

[TABLE]

A vague, uninformative prior could also be used for $\boldsymbol{W}$ but as this was found to have no measurable effect on the results and in order to lighten the notation, we completely omit it.

2.4 Model training

The probabilistic graphical model described thus far is depicted in Figure 4. Given a dataset consisting of FGM input/output pairs (i.e. Stokes-flow model simulations) $\mathcal{D}=\left\{\boldsymbol{\lambda}_{f}^{(n)},\boldsymbol{u}_{f}^{(n)}\right\}_{n=1}^{N}$, the model likelihood can be written as (see Equation (2))

[TABLE]

where $\text{diag}[\boldsymbol{\tau}_{c}]=\text{diag}[\boldsymbol{\sigma}_{c}^{-2}]$ is a diagonal precision matrix and $\boldsymbol{\lambda}_{c}^{(n)}$ are the latent variables that encode the coarse-grained material properties for each sample $n$. Given the aforementioned priors, the posterior of the model parameters is

[TABLE]

where $\tilde{\boldsymbol{\theta}}_{c}=\left\{\tilde{\boldsymbol{\theta}}_{c,m}\right\}_{m=1}^{\text{dim}(\boldsymbol{\lambda}_{c})}$. The posterior over all variables $\boldsymbol{\theta}=\left\{\left\{\boldsymbol{\lambda}_{c}^{(n)}\right\}_{n=1}^{N},\right.$ $\left.\tilde{\boldsymbol{\theta}}_{c},\boldsymbol{\tau}_{c},\boldsymbol{\tau}_{cf},\boldsymbol{\gamma}\vphantom{\left\{\boldsymbol{\lambda}_{c}^{(n)}\right\}_{n=1}^{N}}\right\}$ is given by

[TABLE]

where $p(\mathcal{D})$ is the model evidence of the postulated coarse-graining model Gelman et al. (2003). This is a pivotal quantity in model validation that as we show later will serve as the main driver for the refinement/enhancement of the CGM.

Evaluating the integrals in Equation (21) is analytically intractable, making necessary the use of suitable (approximate) inference methods that are capable of giving fast and accurate estimates of the above posterior. In this work, we employ Stochastic Variational Inference (SVI, Paisley et al. (2012); Hoffman et al. (2013)) methods which produce closed-form approximations of the true posterior $p(\boldsymbol{\theta}|\mathcal{D})$ and simultaneously of the model evidence $p(\mathcal{D})$. In contrast to sampling-based procedures (e.g. MCMC, SMC), SVI yields biased estimates at the benefit of computational efficiency and computable convergence objectives in the form of the Kullback-Leibler (KL) divergence between the approximation $Q(\boldsymbol{\theta})$ and the true posterior $p(\boldsymbol{\theta}|\mathcal{D})$ Blei et al. (2017):

[TABLE]

We note that the $\log$ evidence $\log p(\mathcal{D})$ can be decomposed as

[TABLE]

where

[TABLE]

Since $\text{KL}(Q(\boldsymbol{\theta})||p(\boldsymbol{\theta}|\mathcal{D}))\geq 0$, $\mathcal{F}(Q)$ provides a rigorous lower bound to $\log p(\mathcal{D})$, called the evidence lower bound (ELBO). Hence minimizing $\text{KL}(Q(\boldsymbol{\theta})||p(\boldsymbol{\theta}|\mathcal{D}))$ $\geq 0$ is equivalent to maximizing the ELBO $\mathcal{F}(Q)$.

We employ a mean-field approximation to the full posterior based on the following decomposition:

[TABLE]

If $\boldsymbol{\theta}_{k}$ denotes an arbitrary subset of the parameters $\boldsymbol{\theta}$ above (i.e. if $\boldsymbol{\theta}_{k}$ is effectively one of the component variables $\tilde{\theta}_{c,jm},\tau_{c,m},\gamma_{jm},\tau_{cf,i},\boldsymbol{\lambda}_{c}^{(n)}$), then the optimal $Q_{\boldsymbol{\theta}_{k}}(\boldsymbol{\theta}_{k})$ can be found by setting the first order variation of $\mathcal{F}(Q)$ to zero yielding Beal and Ghahramani (2003)

[TABLE]

where $\left<~{}\cdot~{}\right>_{\ell\neq k}$ denotes expectation w.r.t. all $Q_{\boldsymbol{\theta}_{\ell}}$’s except for $Q_{\boldsymbol{\theta}_{k}}$. We emphasize that every $Q_{\boldsymbol{\theta}_{k}}$ implicitly depends on all other $Q_{\boldsymbol{\theta}_{\ell}}$’s, which means we need to self-consistently cycle over all $k$ and update $Q_{\boldsymbol{\theta}_{k}}$ given expected values w.r.t. all $Q_{\boldsymbol{\theta}_{\ell}},~{}\ell\neq k$, until convergence is attained. We provide complete details of the derivations in A and summarize here the closed-form update equations

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\left<~{}\cdot~{}\right>$ denotes expected values w.r.t. $Q(\left\{\boldsymbol{\lambda}_{c}^{(n)}\right\}_{n=1}^{N},\tilde{\boldsymbol{\theta}}_{c},\boldsymbol{\tau}_{c},\boldsymbol{\tau}_{cf},\boldsymbol{\gamma})$ as specified in (LABEL:eq:meanfield). Due to the choice of conjugate priors, most expected values are given in closed form,

[TABLE]

The optimal form of the remaining approximate posteriors, i.e. $Q_{\boldsymbol{\lambda}_{c}^{(n)}}$, is given by (based on Equation (27))

[TABLE]

which is analytically intractable due to the expectations with respect to the output of the CGM $\boldsymbol{u}_{c}(\boldsymbol{\lambda}_{c}^{(n)})$ (i.e. of the Darcy-flow simulator). Instead of approximating the expectations above by Monte Carlo, we retain the variational character of the algorithm by employing approximations to the optimal $Q_{\boldsymbol{\lambda}_{c}^{(n)}}$ which take the form of multivariate Gaussians with a diagonal covariance, i.e.

[TABLE]

and find the optimal values of the parameters $\boldsymbol{\mu}_{\boldsymbol{\lambda}_{c}}^{(n)},\boldsymbol{\sigma}_{\boldsymbol{\lambda}_{c}}^{(n)}$ via black-box variational inference Ranganath et al. (2014); Hoffman et al. (2013); Paisley et al. (2012), i.e. by minimizing the KL-divergence between $\tilde{Q}_{\boldsymbol{\lambda}_{c}^{(n)}}$ and the optimal $Q_{\boldsymbol{\lambda}_{c}^{(n)}}$ in Equation (33),

[TABLE]

The computation of the objective above as well as gradients with respect to $\boldsymbol{\mu}_{\boldsymbol{\lambda}_{c}}^{(n)},(\boldsymbol{\sigma}_{\boldsymbol{\lambda}_{c}}^{(n)})^{2}$ involves expectations of the CGM outputs as well as their respective gradients $\partial\boldsymbol{u}_{c}(\boldsymbol{\lambda}_{c}^{(n)})/\partial\boldsymbol{\lambda}_{c}^{(n)}$. The latter are efficiently computed using adjoint formulations (see e.g. Heinkenschloss (2008); Constantine et al. (2014)). In order to minimize the noise in the expectations involved, we apply the the reparametrization trick Kingma and Welling (2013) and supply the stochastic gradients to the adaptive moment estimation optimizer (ADAM Kingma and Ba (2014)). We note that this stochastic optimization problem needs to be run for every data point $n$ in every training iteration, but is entirely parallelizable in $n$. After optimization, the expected values $\left<\lambda_{c,m}^{(n)}\right>,\left<(\lambda_{c,m}^{(n)})^{2}\right>$ are readily given by the first and second moments of the resulting Gaussian specified in Equation (34), whereas the expected values $\left<u_{c,k}(\boldsymbol{\lambda}_{c}^{(n)})\right>,\left<u_{c,k}^{2}(\boldsymbol{\lambda}_{c}^{(n)})\right>$ can efficiently be obtained by direct Monte Carlo. Full details about the performed black-box VI are given in A.1. All steps of the training process are summarized in Algorithm 1.

2.5 Model predictions

A main feature of the coarse-grained model presented in this work is the capability to produce probabilistic predictions that quantify uncertainty both due to limited training data and due to limited model complexity, i.e. the information loss happening during the coarse-graining process going $\boldsymbol{\lambda}_{f}\to\boldsymbol{\lambda}_{c}\to\boldsymbol{u}_{c}\to\boldsymbol{u}_{f}$. In particular, given training data $\mathcal{D}$ and a new FGM input $\boldsymbol{\lambda}_{f}$ (which is not included in the training data $\mathcal{D}$), the predictive posterior density $p_{\text{pred}}(\boldsymbol{u}_{f}|\boldsymbol{\lambda}_{f},\mathcal{D})$ for the corresponding FGM output $\boldsymbol{u}_{f}$ is

[TABLE]

with the decoding/encoding distributions $p_{c},p_{cf}$ as described in sections 2.3.1, 2.3.2. We make use of the variational approximation (section 2.4) to the posterior

[TABLE]

in order to rewrite the predictive posterior as

[TABLE]

By drawing samples from the approximate posterior, the integration above can be performed efficiently with Monte Carlo since the cost of each sample is very small as it entails solely CGM output evaluations. We note though that some of the integrations involved in the computation of lower-order statistics can also be done (semi)-analytically. In particular, we note that the predictive posterior mean $\boldsymbol{\mu}_{\text{pred}}(\boldsymbol{\lambda}_{f})$ is given by:

[TABLE]

where $\left<\boldsymbol{u}_{c}(\boldsymbol{\lambda}_{c})\right>$ denotes the predictive posterior average of the CGM output computed with Monte Carlo by ancestral sampling. Similarly, one can compute the posterior predictive covariance of the FGM output vector $\boldsymbol{u}_{f}$, which we denote by $\text{Cov}_{\text{pred}}(\boldsymbol{\lambda}_{f})$, as

[TABLE]

where $\left<\left(\boldsymbol{u}_{c}(\boldsymbol{\lambda}_{f})-\left<\boldsymbol{u}_{c}(\boldsymbol{\lambda}_{f})\right>\right)\left(\boldsymbol{u}_{c}(\boldsymbol{\lambda}_{f})-\left<\boldsymbol{u}_{c}(\boldsymbol{\lambda}_{f})\right>\right)^{T}\right>$ denotes the predictive posterior covariance of the CGM output vector $\boldsymbol{u}_{c}$ which can be estimated by Monte Carlo and $\left<\tau_{cf,i}^{-1}\right>$ correspond to the posterior averages of reciprocal precisions which can be computed in closed form based on Equation (30).

The computational steps to generate predictive samples from the generally non-Gaussian $p_{\text{pred}}(\boldsymbol{u}_{f}|\boldsymbol{\lambda}_{f},\mathcal{D})$ for an unseen microstructural input $\boldsymbol{\lambda}_{f}$ is summarized in Algorithm 2.

The training (offline) and prediction (online) stages are represented graphically in the workflow diagram in Figure 5.

2.5.1 Model performance metrics

Apart from computational efficiency, the key objective of any surrogate model is its predictive accuracy, i.e. its closeness to the true FGM solution $\boldsymbol{u}_{f}(\boldsymbol{\lambda}_{f})$. In order to assess this, we consider separate test datasets $\mathcal{D}_{\text{test}}=$ $\left\{\boldsymbol{\lambda}_{f}^{(n)},\boldsymbol{u}_{f}^{(n)}\right\}_{n=1}^{N_{\text{test}}}$, where the $\boldsymbol{\lambda}_{f}^{(n)}\sim p(\boldsymbol{\lambda}_{f})$ are drawn from the same distribution $p(\boldsymbol{\lambda}_{f})$ that was used for the training set $\mathcal{D}$. We report metrics that reflect both the accuracy of point estimates of the predictive posterior (e.g. $\boldsymbol{\mu}_{\text{pred}}(\boldsymbol{\lambda}_{f})$ in Equation (39)) as well as the whole density $p_{\text{pred}}(\boldsymbol{u}_{f}|\boldsymbol{\lambda}_{f},\mathcal{D})$ (Equation (38)).

In particular, we compute the coefficient of determination $R^{2}$ Zhang (2017)

[TABLE]

where $\bar{\boldsymbol{u}}_{f}=\frac{1}{N_{\text{test}}}\sum_{n=1}^{N_{\text{test}}}\boldsymbol{u}_{f}^{(n)}$ is the sample mean of the FGM output $\boldsymbol{u}_{f}$ over the test dataset. The denominator normalizes the point deviation in the numerator by the actual variability of the true outputs. In the case of perfect prediction, $\boldsymbol{\mu}_{\text{pred}}(\boldsymbol{\lambda}_{f}^{(n)})=\boldsymbol{u}_{f}^{(n)}~{}\forall n$, the second term in (41) vanishes and $R^{2}=1$ which is the largest value it can attain. Smaller values of $R^{2}$ indicate either larger deviation of the mean prediction from the truth or smaller variability in the $\boldsymbol{u}_{f}$ components.

The second metric we employ is the mean log likelihood $MLL$ Zhu and Zabaras (2018) given by

[TABLE]

We note that when $\log p_{\text{pred}}(\boldsymbol{u}_{f}|\boldsymbol{\lambda}_{f},\mathcal{D})$ degenerates to a discrete density centered at the true $\boldsymbol{u}_{f}(\boldsymbol{\lambda}_{f}),~{}\forall\boldsymbol{\lambda}_{f}$ i.e. not only it predicts perfectly the FGM output but there is no predictive uncertainty, then the $MLL$ attains the highest possible value of [math]. In the other extreme if for one or more of the test samples $\boldsymbol{\lambda}_{f}^{(n)}$, $p_{\text{pred}}\to 0$ at the true $\boldsymbol{u}_{f}(\boldsymbol{\lambda}_{f})$, then $MLL$ attains its lowest possible value of $-\infty$. Hence, $MLL$ attempts to measure the quality of the whole predictive posterior $p_{\text{pred}}$. Given that the integration in Equation (38) is analytically intractable and no closed-form expression for $p_{\text{pred}}$ can be established, we evaluate $MLL$ under the assumption that $p_{\text{pred}}$ can be adequately approximated by a multivariate Gaussian with mean $\boldsymbol{\mu}_{\text{pred}}$ (Equation (39)) and a diagonal covariance which consists of the diagonal entries of $\text{Cov}_{\text{pred}}(\boldsymbol{\lambda}_{f})$ (Equation (40)) which we represent with the vector $\boldsymbol{\sigma}^{2}_{\text{pred}}(\boldsymbol{\lambda}_{f})$. In this case:

[TABLE]

In the numerical experiments discussed in section 3, $N_{\text{test}}=1024$ is used, leading to negligible Monte Carlo error due to variation of test samples.

2.6 Numerical complexity analysis

For analysis of the numerical complexity of the proposed model, it is vital to distinguish between training (offline) and prediction (online) phases, see also Figure 5. There are basically six quantities that are relevant for the computational cost of both stages. These are the number of training data $N$, the dimension of the FGM input $\text{dim}(\boldsymbol{\lambda}_{f})$, the dimension of the FGM output $\text{dim}(\boldsymbol{u}_{f})$, the dimension of CGM input $\text{dim}(\boldsymbol{\lambda}_{c})$, the dimension of the CGM output $\text{dim}(\boldsymbol{u}_{c})$ and the number of feature functions per latent space variable $\text{dim}(\tilde{\boldsymbol{\theta}}_{c})$.

As it can be seen in Algorithm 1, the complexity of the training stage depends linearly on the number of training data $N$ (for-loop), but is fully parallelizable over every data point. The prediction stage only sees the final (approximate) posterior $p(\tilde{\boldsymbol{\theta}}_{c},\boldsymbol{\tau}_{c},\boldsymbol{\tau}_{cf}|\mathcal{D})\approx Q_{\tilde{\boldsymbol{\theta}}_{c}}(\tilde{\boldsymbol{\theta}}_{c})Q_{\boldsymbol{\tau}_{c}}(\boldsymbol{\tau}_{c})Q_{\boldsymbol{\tau}_{cf}}(\boldsymbol{\tau}_{cf})$, so that the scaling is independent of $N$. Both training and prediction phases scale linearly with the fine scale output dimension $\text{dim}(\boldsymbol{u}_{f})$, since they only involve vector operations. For the same reason, the scaling with respect to $\text{dim}(\boldsymbol{\lambda}_{c})$ is linear. Scaling with $\text{dim}(\boldsymbol{u}_{c})$ depends on the CGM solver i.e. in the worst case, for the linear Darcy model with a direct solver, the scaling is $\mathcal{O}\left(\text{dim}((\boldsymbol{u}_{c}))^{3}\right)$. We apply a sparse banded solver, which scales as $\mathcal{O}\left(\text{dim}((\boldsymbol{u}_{c}))\right)$. Scaling with the dimension of microstructural inputs $\text{dim}(\boldsymbol{\lambda}_{f})$ is dependent on the set of applied feature functions $\boldsymbol{\varphi}_{m}$ and ranges from $\mathcal{O}(1)$ (e.g. constant feature, i.e. $\varphi(\boldsymbol{\lambda}_{f})=1$) to $\mathcal{O}\left((\text{dim}(\boldsymbol{\lambda}_{f}))^{2}\right)$ (e.g. mean mutual distance of circular exclusions). Finally, scaling with respect to $\text{dim}(\tilde{\boldsymbol{\theta}}_{c,m})$ is $\mathcal{O}\left((\text{dim}(\tilde{\boldsymbol{\theta}}_{c,m}))\right)$ see Equation (31). Furthermore, in the prediction stage, it is possible to omit components of $\tilde{\boldsymbol{\theta}}_{c,m}$ that have been pruned out by the applied sparsity prior, thereby avoiding the cost connected to the evaluation of corresponding feature functions. The scaling characteristics of training and prediction phases is summarized in table 1.

2.7 Automatic adaptive refinement

As already mentioned in section 2.4, the evidence lower bound (ELBO) is an approximation to the model’s $\log-$evidence. The latter balances the model’s fit to the data with the model’s complexity and serves as a parsimonious metric of how well the model can explain the training data Murray and Ghahramani (2005); Rasmussen and Ghahramani (2001). This quantity can therefore be used not only for scoring different competing models which have been trained on the same data, but also as the objective to any refinement of an existing model. According to Equation (25), the ELBO is given by the expected value

[TABLE]

where $p(\boldsymbol{\theta},\mathcal{D})$ is given in Equation (22) and $Q(\boldsymbol{\theta})$ is the variational approximation to the posterior $p(\boldsymbol{\theta}|\mathcal{D})$ as defined in (LABEL:eq:meanfield). Since all of the expected values needed in the right side of the above equation are estimated during model training, the ELBO $\mathcal{F}(Q)$ can be evaluated explicitly and be used to also monitor convergence behavior during training. In order to keep this section uncluttered, the final expression for the ELBO is given in D.

While one can envision many ways of refining the coarse-grained model, we focus here on the latent variables $\boldsymbol{\lambda}_{c}$. These constitute the bottleneck through which information that the FGM input $\boldsymbol{\lambda}_{f}$ provides about the FGM output $\boldsymbol{u}_{f}$ is squeezed. Physically, they relate to the effective permeability of the Darcy-based CGM and assuming a model according to equations (9)–(10), they are inherently tied with the piecewise-constant discretization adopted, i.e. each component $\lambda_{c,m}$ relates to the Darcy permeability in a square subregion $\Omega_{m}$ of the problem domain.

The question we would like to address is the following: If one was to refine the piecewise-constant representation of the Darcy permeability field in equations (9)–(10) by subdividing one square cell $\Omega_{m}$ of the existing partition into smaller subregions, then which $\Omega_{m}$ would yield the best possible model? We note that such a refinement would lead to an increase in the dimension $\text{dim}(\boldsymbol{\lambda}_{c})$ of the bottleneck variables which, in principle, should lead to more information form $\boldsymbol{\lambda}_{f}$ being retained.

Given a fixed refinement mechanism of an $\Omega_{m}$, say by subdividing into 4 equal squares (in two-dimensions), a brute force solution strategy would be to refine one of the existing $\Omega_{m}$ at a time, train the new model, compute its ELBO and compare with the others. Obviously such a strategy is computationally cumbersome particularly when the number of existing subregions $N_{\text{cells,c}}$ is high. The cost would be even greater if one was to consider all possible combinations of two or more subregions that could be refined at the same time.

We thus propose a procedure based on the decomposition of the ELBO $\mathcal{F}(Q)$ into contributions from different subregions $\Omega_{m}$ (and therefore $\lambda_{c,m}$) according to $\mathcal{F}(Q)=\sum_{m=1}^{N_{\text{cells,c}}}\mathcal{F}_{m}(Q)+\mathcal{H}(Q)$, where with $\mathcal{H}(Q)$ we denote terms that cannot be assigned to any of the subregions $\Omega_{m}$. As the goal is ultimately to increase the ELBO of the model, we propose selecting for refinement the cell $m$ with the lowest contribution $\mathcal{F}_{m}(Q)$ and therefore contributes the least in explaining the data. This strategy is analogous to the one advocated in Ghahramani and Beal (2000) for refining components in a mixture of factor analyzers. Hence we use $\mathcal{F}_{m}(Q)$ as the scoring function of each subregion/cell. The derivation of $\mathcal{F}_{m}(Q)$ is contained in D. and we include here only the final expression which, for a model with tied $\gamma_{jm}=\gamma_{j}$, is

[TABLE]

where $(\sigma_{\lambda_{c,m}}^{(n)})^{2}$ is the posterior variance of $\lambda_{c,m}^{(n)}$ (see equations (34) and (35)), $\tilde{c},\tilde{d}_{m}$ are given in Equation (29) and $\sigma_{\tilde{\theta}_{c,jm}}^{2}$ is the posterior variance of the parameters $\theta_{c,jm}$ given in Equation (31).

Keeping in mind that the subregion/cell $\Omega_{m}$ with the smallest $\mathcal{F}_{m}$ is proposed for refinement, the scoring function defined by Equation (45) favours cells where the posterior variances $(\sigma_{\lambda_{c,m}}^{(n)})^{2}$ of the $\lambda_{c,m}^{(n)}$ tend to be small, i.e. tight (approximate) posterior distribution $Q_{\lambda_{c,m}^{(n)}}(\lambda_{c,m}^{(n)})$, whereas at the same time, the expected decoder misfit $\sum_{n=1}^{N}\left<\left(\lambda_{c,m}^{(n)}-\tilde{\boldsymbol{\theta}}_{c,m}^{T}\boldsymbol{\varphi}_{m}(\boldsymbol{\lambda}_{f}^{(n)})\right)^{2}\right>$ (see Equation (29)) should be large but there should be high certainty about the feature coefficients $\tilde{\boldsymbol{\theta}}_{c,m}$, i.e. low $\sigma_{\tilde{\theta}_{c,jm}}^{2}$. We finally note that once the cell $m$ with the lowest $\mathcal{F}_{m}(Q)$ has been identified and split, training is restarted with all model parameters initialized to their previously converged values. We illustrate these steps by example in section 3.6.

3 Numerical examples

The numerical experiments carried out in this paper are aimed at replacing the expensive Stokes flow FE simulation defined by Equations (3)–(4) by an inexpensive, Bayesian, coarse-grained model built around the Darcy-flow skeleton as presented in sections 2.2–2.6. We begin this section with a clear specification of the FGM data generation and the model distributions $p_{c},p_{cf}$. Subsequently, we show that close to the homogenization limit (i.e. infinite scale separation) defined by Equation (5), the proposed coarse-grained model is capable to accurately reconstruct the FGM solution by learning (from a small number of FGM training samples) the requisite effective properties without ever solving any of the equations homogenization theory describes. Thereafter, we investigate the coarse-grained model’s predictive quality as a function of the number of training data $N$ and the latent space dimensionality $\text{dim}(\boldsymbol{\lambda}_{c})$ for microstructural data far from the homogenization limit (without scale separation). We show that the model indeed exhibits feature sparsity in learning the effective Darcy permeability field $\boldsymbol{K}(\boldsymbol{x},\boldsymbol{\lambda}_{c})$. We then apply the surrogate to an uncertainty propagation (UP) problem and finally demonstrate the performance of the proposed adaptive refinement objective and associated algorithm. Source code as well as examples can be found at: https://github.com/congriUQ/physics_aware_surrogate

3.1 Experimental setup

3.1.1 Fine scale microstructural data

Throughout the following experiments, we use randomly distributed, non - overlapping polydisperse spherical exclusions from a unit square domain $\Omega=[0,1]^{2}$ as depicted in Figure 2. The total number of exclusions $N^{(n)}_{\text{ex}}$ of a material sample $n$ as well as the radii $r_{\text{ex},i}^{(n)},i=1,\ldots,N^{(n)}_{\text{ex}}$, are sampled from a log-normal distribution,

[TABLE]

The coordinates of the center of each exclusion $\boldsymbol{x}_{\text{ex},i}^{(n)}$ are drawn from a Gaussian process with squared exponential kernel of length scale $l_{\boldsymbol{x}}$ warped with the logistic sigmoid function $S(z)=1/(1+e^{-l_{s}z})$,

[TABLE]

where $N_{\rho^{(n)}}=\int\rho^{(n)}_{\boldsymbol{x}_{\text{ex}}}(\boldsymbol{x})\boldsymbol{x}$ and a squared exponential covariance $k_{\boldsymbol{x}}=$

$\exp\left\{-\left|\boldsymbol{x}-\boldsymbol{x}^{\prime}\right|^{2}/l_{\boldsymbol{x}}^{2}\right\}$777Note that the resulting distribution is distorted by the fact that a sampled exclusion is rejected if it overlaps with an exclusion which has been inserted before.. Also, the (location-dependent) $\log$-normal mean of exclusion radii $\mu_{r}(\boldsymbol{x})$ is drawn from a Gaussian process for every sample according to

[TABLE]

with $k_{r}=\exp\left\{-\left|\boldsymbol{x}-\boldsymbol{x}^{\prime}\right|^{2}/l_{r}^{2}\right\}$ the corresponding covariance kernel. Four microstructural samples generated according to the described distributions are depicted in Figure 6. After generation of a microstructure, the domain is pixelized into a grid of size $256\times 256$. The corresponding binary vector representing the phase of each pixel makes up the FGM input $\boldsymbol{\lambda}_{f}$, i.e.:

[TABLE]

The sampled microstructres are passed to the FEniCS mshr Logg et al. (2012) module to generate a triangular finite element mesh containing roughly $256\times 256=65,536$ vertices. The output of interest, i.e. $\boldsymbol{u}_{f}$, consists of the the interpolated pressure values on a regular grid of $129\times 129$ i.e.:

[TABLE]

Depending on the total number of exclusions $N_{\text{ex}}^{(n)}$, generation of a single mesh takes in between $\sim 2$ hours (for $N_{\text{ex}}^{(n)}\approx 1,000$) up to 10 days (for $N_{\text{ex}}^{(n)}\approx 20,000$) of computation time on a single Intel Xeon E5-2620 (2.00 GHz) CPU. The average solution for each FGM (i.e. Stokes-flow) PDE-solve takes 1512s $\pm$ 5.7s on the same hardware.

We use boundary conditions on the FGM pressure $P_{bc}(\boldsymbol{x})$ and velocity fields $\boldsymbol{V}_{bc}(\boldsymbol{x})$ of the form

[TABLE]

where $\Gamma_{P}=\boldsymbol{0}$ denotes the origin of the domain, i.e. the point $\boldsymbol{x}=\boldsymbol{0}$. The boundary condition coefficients $a_{x},a_{y}$ and $a_{xy}$ may assume different values in the following experiments and are therefore specified later.

3.1.2 The CGM $\boldsymbol{u}_{c}(\boldsymbol{\lambda}_{c})$

For the Darcy-type reduced order model $\boldsymbol{u}_{c}(\boldsymbol{\lambda}_{c})$, we use a finite-element solver to Equations (8) as described in section 2.2, operating on a regular $16\times 16$ square grid endowed with bilinear shape functions. The coarse-grained random permeability field $\boldsymbol{K}(\boldsymbol{x},\boldsymbol{\lambda}_{c})$ of the CGM Darcy solver is assumed to have the form defined by equations (9)–(10) throughout this section.

3.1.3 The encoder distribution $p_{c}$

For the encoder distribution $p_{c}(\boldsymbol{\lambda}_{c}|\boldsymbol{\lambda}_{f},\boldsymbol{\theta}_{c})$, we assume a model as defined by equations (12)–(13), with the restriction that

[TABLE]

where $\boldsymbol{\lambda}_{f}^{(m)}$ denotes the subset of the microstructure $\boldsymbol{\lambda}_{f}$ which corresponds to the permeability field cell $\Omega_{m}$. Equation (53) essentially states that the same set of features $\boldsymbol{\varphi}_{m}(\boldsymbol{\lambda}_{f})=\boldsymbol{\varphi}(\boldsymbol{\lambda}_{f})=\left\{\boldsymbol{\varphi}_{\text{glob}}(\boldsymbol{\lambda}_{f}),\boldsymbol{\varphi}_{\text{loc}}(\boldsymbol{\lambda}_{f}^{(m)})\right\}$ are used for every permeability field cell $m$, where $\boldsymbol{\varphi}_{\text{glob}},\boldsymbol{\varphi}_{\text{loc}}$ denote features evaluated on the global domain and on a local cell permeability field cell $m$, respectively. It is noted though that the coefficients $\tilde{\theta}_{c,jm}$ are generally different for different cells $m$. A Table of all feature functions used in the experiments is given in B.

Prior on $\tilde{\theta}_{c,jm}$

The prior defined in Equation (16) is restricted such that the precision hyperparameters $\gamma_{jm}=\gamma_{j}$, i.e. all permeability field cells $m$ share the same set of hyperparameters $\boldsymbol{\gamma}$. This leads to the same set of activated features, i.e. features with nonzero, but generally different $\tilde{\theta}_{c,jm}$ for different cells $m$, such that permeability cells are enabled to share information across the whole domain $\Omega$.

3.1.4 The decoder distribution $p_{cf}$

The decoder distribution we use in the experiments is defined by equations (14)–(15). As the FGM model output $\boldsymbol{u}_{f}$, we choose to use the Stokes flow pressure response $u_{f,i}=P(\boldsymbol{x}_{i})$ linearly interpolated onto a square rectangular grid $\boldsymbol{x}_{i}\in G^{(f)}=129\times 129$ in the unit square domain $\Omega=[0,1]^{2}$. Although parameter uncertainty on the coarse-to-fine projection matrix $\boldsymbol{W}$ could readily be included to the model by specifying a suitable prior distribution $p(\boldsymbol{W})$, it has proven to be beneficial (in terms of predictive quality) in the low-data regime $N\lesssim 100$ to fix $\boldsymbol{W}$ to the CGM Darcy solver shape function interpolant, i.e.

[TABLE]

where $\psi_{c,j}(\boldsymbol{x})$ is the Darcy solver finite element shape function associated with the CGM degree of freedom $u_{c,j}$. The precisions $\tau_{cf,i}$ remain a free parameter and give the inverse reconstruction variance for $u_{f,i}=P(\boldsymbol{x}_{i})$.

3.2 Case 1: Scale-separation (homogenization limit) & high data variability

To reveal the full potential of the proposed approach, we use input data that contains microstructures $\boldsymbol{\lambda}_{f}$ that fulfill the homogenization limit Equation (5) as much as it is computationally affordable, whilst having large variability, i.e. high $\sigma_{\text{ex}}^{2},\sigma_{r}^{2}$, see equations (46)–(47). We therefore sample FGM inputs $\boldsymbol{\lambda}_{f}$ from the model presented earlier with the following parameter values:

[TABLE]

The boundary conditions employed are of the form given in Equation (52) with

[TABLE]

We generate three training data sets of size $N=128$ and report the averages of the performance metrics defined in 2.5.1 over the three models trained. We find

[TABLE]

evaluated on a test set of size $N_{\text{test}}=1024$. Predictions for four random test samples are illustrated in Figure 7. We use $\text{dim}(\boldsymbol{\lambda}_{c})=16$ corresponding to a uniform $4\times 4$ square discretization of the problem domain $\Omega$.

In C, we report the results on the PCA decomposition of the FGM outputs $\boldsymbol{u}_{f}$. One observes a rapid decay of the corresponding eigenvalues in Figure 18 which would suggest that the problem could be addressed by methods that project the governing equations on the low-dimensional subspace spanned by the first few eigenvectors (e.g. reduced-basis (RB) methods Hesthaven et al. (2016); Quarteroni et al. (2016)). We would like to emphasize though that it would be unclear how the Stokes/Darcy homogenization process should be accomplished for the high-dimensional random porous domain described by $\boldsymbol{\lambda}_{f}$, how to find the RB solution coefficients and if a standard Galerkin projection procedure would indeed be computationally more efficient. Moreover, we note that a RB-based Darcy CGM could readily be included as a substitute for the $\boldsymbol{u}_{c}(\boldsymbol{\lambda}_{c})$ CGM as was defined in section 2.2, which might be a possible extension to investigate in the future. Apart from that, we emphasize that fast decay of PCA eigenmodes of the output data $\boldsymbol{u}_{f}$ does not affect the high effective dimension of inputs $\boldsymbol{\lambda}_{f}$.

3.3 Case 2: No scale separation (far from homogenization limit) and high output variability

In this setting, we consider microstructures that lead to PCA decompositions of the FGM outputs $\boldsymbol{u}_{f}$ that are not concentrated on a very small subset of eigenvectors. In particular we generate microstructures based on the following parameters:

[TABLE]

see Figure 18 for a PCA analysis of the output data $\left\{\boldsymbol{u}_{f}^{(n)}\right\}_{n=1}^{2048}$, revealing the effective output dimension. Boundary conditions are again fixed to $a_{x}=a_{y}=1,a_{xy}=0$. We consider three cases in terms of $\text{dim}(\boldsymbol{\lambda}_{c})$ which correspond to regular, square discretization of size $2\times 2$ (i.e. $\text{dim}(\boldsymbol{\lambda}_{c})=4$), $4\times 4$ (i.e. $\text{dim}(\boldsymbol{\lambda}_{c})=16$) and $8\times 8$ (i.e. $\text{dim}(\boldsymbol{\lambda}_{c})=64$). While the PDE discretization is independent of the representation of $\boldsymbol{\lambda}_{c}$, it needs to be able to resolve the corresponding Darcy permeability field. To that end, we employ a regular square mesh of $16\times 16$ bilinear elements for the CGM model, leading to $\text{dim}(\boldsymbol{u}_{c})=288$. One CGM evaluation requires $(1.04\pm 2)\times 10^{-3}~{}s$ of computational time on a single Intel Xeon E5-2620 (2.00 GHz) CPU, as opposed to the $1512\pm 5.7s$ needed for each FGM solve. For predictive purposes (section 2.5), we use $N_{\text{samples}}=100$ for the estimation e.g. of $\boldsymbol{\mu}_{\text{pred}},\boldsymbol{\sigma}_{\text{pred}}^{2}$ as defined in Equations (39) and (40).

Figure 8 depicts both error metrics $R^{2}$ and $MLL$ as a function of the number of training samples $N$ for the different $\text{dim}(\boldsymbol{\lambda}_{c})$ as described above. One observes that with only $N\approx 32$ training samples (i.e. FGM simulations), the asymptotic values for $R^{2}$ and $MLL$ are attained. Furthermore, it is interesting to note that the highest $\text{dim}(\boldsymbol{\lambda}_{c})=64$ achieves the second best $R^{2}$ score after $\text{dim}(\boldsymbol{\lambda}_{c})=16$. This could be attributed to the higher information loss taking place using the local feature functions $\varphi_{\text{loc}}$ on a finer grid in Equation (53).

Figure 9 shows three predictive means $\boldsymbol{\mu}_{\text{pred}}$ (blue) $\pm$ one predictive standard deviation $\boldsymbol{\sigma}_{\text{pred}}$ (transparent gray) on the same microstructure $\boldsymbol{\lambda}_{f}$, but using different number of training data $N=4,16,32$, computed with $\text{dim}(\boldsymbol{\lambda}_{c})=16$. One can observe how the true solution (colored surface) is better captured the more data $N$ is used for training.

3.3.1 Predictive variance $\boldsymbol{\sigma}^{2}_{\text{pred}}$ & $L^{2}$-error

In Figure 10, the predictive uncertainty $\boldsymbol{\sigma}_{\text{pred}}$ (top row) and the corresponding ‘true’ $L^{2}$-error $|\boldsymbol{u}_{f}(\boldsymbol{\lambda}_{f})-\boldsymbol{\mu}_{\text{pred}}(\boldsymbol{\lambda}_{f})|$ are plotted for four different microstructures $\boldsymbol{\lambda}_{f}$ drawn according to the same distribution as in section 3.3, using a $4\times 4$ permeability field discretization and $N=128$ training samples. To provoke different predictive errors in different regions in the domain, we use non-homogeneous flux boundary conditions of type (52) with $a_{x}=a_{y}=0,a_{xy}=-1$. It can be observed that $\boldsymbol{\sigma}_{\text{pred}}$ is indeed higher in samples/regions where the true $L^{2}$ error $|\boldsymbol{u}_{f}(\boldsymbol{\lambda}_{f})-\boldsymbol{\mu}_{\text{pred}}(\boldsymbol{\lambda}_{f})|$ is higher, whereby it is incapable to resolve errors on a smaller scale than the one defined by the effective permeability field discretization.

3.3.2 Effective Darcy permeabilities

Further insight on the numerical homogenization process taking place in the $\boldsymbol{\lambda}_{f}\mapsto\boldsymbol{\lambda}_{c}$ encoding process is provided by Figure 11, which shows the mean effective permeability field $\left<\boldsymbol{K}(\boldsymbol{x},\boldsymbol{\lambda}_{c})\right>_{Q_{\boldsymbol{\lambda}_{c}}}$ (bottom row, see section 3.1.2) and the corresponding microstructures (top) for four randomly chosen training samples of a $4\times 4$ effective model trained on $N=128$ samples using boundary conditions and microstructure distributions as in 3.3. It is clear that effective permeabilities are lower in cells with higher solid fraction/interface area.

Figure 12 provides information on the activated feature functions. These are determined by the posterior mean of the precision parameters $\gamma_{jm}=\gamma_{j}$ that do not converge to large values (according to Equation (31) when $\left<\gamma_{jm}\right>\rightarrow\infty$ then $\sigma_{\tilde{\theta}_{c,jm}}^{2}\rightarrow 0$). We note that, apart from the bias term, the features activated pertain to the self-consistent approximation Torquato (2001) (infinite contrast limit), the squared area of total solid/fluid interfaces, and void-phase lineal-path function which pertains to the connectivity of the void phase. In Figure 12, $\text{dim}(\boldsymbol{\lambda}_{c})=16$ is the dimension of the CGM’s constitutive parameter vector and a set of $N_{\text{features}}=\text{dim}(\boldsymbol{\gamma})=150$ feature functions $\varphi(\boldsymbol{\lambda}_{f})$ specified in B is used. The left side of the Figure 12 shows the average number of activated features (i.e. with nonzero $\left<\gamma_{j}\right>^{-1}$) as a function of the number of training data $N$, where error bars are due to randomization of training sets. As expected, the number of activated feature functions increases with the number of training data. The right part is an example of the final, converged values of $\left<\gamma_{j}\right>^{-1}$ for $N=128$ training samples, where 5 features are activated (explained in the Figure caption).

3.4 Predictions under different boundary conditions

A particularly useful quality of the proposed coarse-grained model is that boundary conditions can be imprinted deterministically onto the Darcy-based CGM as specified in equations (8) (c-d). As a result, a model trained on boundary conditions, say $\boldsymbol{a}$ (see section 3.1.1), can be used for predictions on different boundary conditions, say $\tilde{\boldsymbol{a}}\neq\boldsymbol{a}$. Good predictive performance in such cases would suggest the ability of the proposed coarse-grained model to extrapolate which presupposes that it has correctly encoded salient physical information. This is a desirable trait in all surrogates which cannot always be achieved with black-box statistical models that attempt to interpolate the input-output map Tripathy and Bilionis (2018); Zhu and Zabaras (2018).

As an example, we use the boundary conditions $\boldsymbol{a}=\{a_{x}=a_{y}=1,a_{xy}=0\}$ (as in section 3.3) and boundary conditions $\tilde{\boldsymbol{a}}=\{\tilde{a}_{x}=\tilde{a}_{y}=0,\tilde{a}_{xy}=-1\}$. We generate $N=128$ FGM input-output training pairs for each of the aforementioned boundary conditions and train a coarse-grained model with $\text{dim}(\boldsymbol{\lambda}_{c})=16$ as previously described. In Table 2, we report the performance metrics $R^{2}$ and $MLL$ as evaluated on test data generated with each of the above boundary conditions. We observe that the off-diagonal entries which correspond to different BCs for training and testing are comparable with the diagonal entries which correspond to identical BCs for training and testing. Indicative test samples and the corresponding predictions are shown in Figure 13. We observe again comparable predictive accuracy between the different combinations of training and test data.

We also perform full randomization of boundary conditions in a further scenario, where different boundary conditions for every single test and training datum are drawn. In particular, we draw $a_{x}\sim\mathcal{N}(0,1)$, $a_{y}\sim\mathcal{N}(0,1)$, $a_{xy}\sim\mathcal{N}(0,1)$ and generate $N=128$ training data (i.e. FGM input-output pairs). Figure 14 depicts predictions over four test cases with randomly selected BCs as above, where it is clearly seen that the surrogate is capable of producing accurate estimates both in terms of the mean as well as the breadth of the predictive uncertainty. The error metrics $R^{2}$ and $MLL$ are evaluated on a $N_{\text{test}}=1024$ test set and yield the following values:

[TABLE]

which are comparable to the values provided earlier.

3.5 Uncertainty propagation problem

The numerical experiment carried out in this section pertains to a simple uncertainty propagation (UP) problem. In a UP setting, one is interested in computing statistics of some quantity of interest (QoI) $f(\boldsymbol{u}_{f})$ of the FGM output $\boldsymbol{u}_{f}$ given some density $p(\boldsymbol{\lambda}_{f})$ for the FGM input $\boldsymbol{\lambda}_{f}$. For the sake of illustration we assume that the QoI is scalar and therefore its density can be written as

[TABLE]

where $p(\boldsymbol{u}_{f}|\boldsymbol{\lambda}_{f})=\delta(\boldsymbol{u}_{f}-\boldsymbol{u}_{f}(\boldsymbol{\lambda}_{f}))$ if $\boldsymbol{u}_{f}$ is computed by running the expensive FGM $\boldsymbol{u}_{f}(\boldsymbol{\lambda}_{f})$. One can approximate $p(f)$ with a histogram where the values in each bin are estimated with Monte Carlo (MC) sampling. Due to slow MC convergence rates and the high cost of the FGM forward solver, an obvious strategy is to replace the FGM by the probabilistic coarse-grained model proposed and more specifically the predictive posterior $p_{\text{pred}}(\boldsymbol{u}_{f}|\boldsymbol{\lambda}_{f},\mathcal{D})$ (section 2.5) as follows:

[TABLE]

To illustrate the performance of such a scheme, we use as the QoI the FGM pressure $P(\boldsymbol{x})$ at $\boldsymbol{x}=\begin{pmatrix}1,&1\end{pmatrix}^{T}$ (i.e. top right corner of the problem domain $\Omega$). The blue line in Figure 15 shows a kernel density estimate of the MC-based histogram of 10,000 FGM runs. The black line is based on a trained coarse-grained model ($\text{dim}(\boldsymbol{\lambda}_{c})=16$) with $N=32$ training samples. The gray shaded area of Figure 15 aims to capture predictive uncertainties about $p(f)$ due to limited data. To that end, instead of $p_{\text{pred}}(\boldsymbol{u}_{f}|\boldsymbol{\lambda}_{f},\mathcal{D})$ as defined in Equation (36), we generate samples of the model parameters $\tilde{\boldsymbol{\theta}}_{c}^{(s)},\boldsymbol{\tau}_{c}^{(s)}$ and $\boldsymbol{\tau}_{cf}^{(s)}$ from their approximate posteriors $Q_{\tilde{\boldsymbol{\theta}}_{c}},Q_{\boldsymbol{\tau}_{c}},Q_{\boldsymbol{\tau}_{cf}}$ and for each of those values we generate samples from

[TABLE]

The mean and $90\%$ credible intervals for the histogram value in each bin are subsequently computed and plotted in Figure 15. One can observe that the credible intervals envelop correctly even the tails of the true, MC-based histogram.

3.6 Illustration of automatic adaptive refinement

In this section, we illustrate the automatic adaptive refinement scheme of the coarse-grained model as explained in section 2.7. We consider as the base model $\text{dim}(\boldsymbol{\lambda}_{c})=4$ which corresponds to a regular $2\times 2$ square grid discretization of the Darcy permeability field (see top left of Figure 17). In order to better illustrate the capabilities of the proposed method, we draw the FGM microstructural samples, i.e. $\boldsymbol{\lambda}_{f}$, from a particular distribution and refer to the corresponding samples as ‘tiled’ microstructures. These exhibit a homogeneous distribution of circular exclusions (solid phase) over the whole domain except for the lower left quadrant, i.e. $\Omega_{ll}=[0,0.5]\times[0,0.5]$. In this region, exclusions are distributed homogeneously on sub-cells of size $0.125\times 0.125$, but with varying random volume fraction as can be seen in Figure 16 where four representative microstructures are depicted. Due to this localized inhomogeneity, it is natural to expect that refinements in $[0,0.5]\times[0,0.5]$ would be most promising. The baseline $\text{dim}(\boldsymbol{\lambda}_{c})=2\times 2$ model is trained on $N=32$ data generated by the Stokes-based FGM.

After training, the cell scoring function defined in equation (45) indicates which of the $2\times 2=4$ cells of constant permeability should be split into four new square sub-cells. After splitting that cell, training continues with all model parameters initialized to their previously found values. To facilitate the problem as much as possible, boundary conditions of type $a_{x}=1$, $a_{y}=a_{xy}=0$ are used, i.e. uniform expected velocity $\left<\boldsymbol{V}(\boldsymbol{x})\right>=\left(1,~{}0\right)^{T}$ all over the domain. Given these boundary conditions and the ‘tiled’ microstructures as described above, it should be beneficial to have a finer CGM permeability field discretization in $\Omega_{ll}$ than in the rest of the domain, where exclusions are distributed homogeneously.

Figure 17 depicts the process of 4 successive splits each of which corresponds to dividing the selected cell into four equal-sized squares. Hence the number of cells/subregions and therefore the dimension of the associated $\text{dim}(\boldsymbol{\lambda}_{c})$ will become $7,10,13$ and $16$ after each of the aforementioned splits. The left part of Figure 17 shows the cell scores $\mathcal{F}_{m}(Q)$ computed as in Equation (45) right before each split. The cell that is subsequently subdivided is the one with the lowest $\mathcal{F}_{m}(Q)$. As expected, the splits are concentrated in the lower-left quadrant $\Omega_{ll}$ and the right part of the Figure shows the corresponding discretization of the permeability field after the 4 splits. The tendency to refine close to the origin $\boldsymbol{x}=\boldsymbol{0}$ (lower left corner) and along the lower boundary $y=0$ may be explained by the fact that the microstructures are constructed such that solid phase exclusions have a minimum distance of $d=0.003$ to the domain boundary, i.e. there always exists a narrow path along the boundaries where fluid can flow. As boundary conditions imply flux from left to right, it is beneficial to refine close to the lower domain boundary. The lower left corner is further distinguished by the fact that the pressure boundary condition $P(\boldsymbol{x}=0)=0$ is implied here.

More importantly, the middle part of 17 shows the evolution of the ELBO over training iterations (blue line). The jumps that are observed after each split are a result of the improvement in the model which is to be expected due to the increased dimension $\text{dim}(\boldsymbol{\lambda}_{c})$. We also compare the ELBO score after 4 splits (where $\dim(\boldsymbol{\lambda}_{c})=16$) with a regular discretization of the Darcy permeability into $4\times 4$ squares in which $\text{dim}(\boldsymbol{\lambda}_{c})=16$ as well, i.e. the dimension of the information bottleneck variables is the same. We note that the adaptive refinement proposed leads to an improved final model.

4 Conclusions

We have presented a physics-informed, fully Bayesian surrogate model to predict response fields of stochastic partial differential equations with high-dimensional stochastic input as those arising in random heterogeneous media. Employing an encoder-decoder architecture with a coarse-grained model (CGM) based on simplified physics as a stencil at its core, it is possible to provide tight and accurate probabilistic predictions even when the training data i.e. fine-grained model simulations, is scarce (i.e. of the order of N\approx 10{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sim}100) and the stochastic input dimension is large (i.e. $\sim 10^{4}$ or larger). The encoding or coarse-graining step is capable of extracting the salient features of the input random field which are most predictive of the fine-scale PDE-response, thereby drastically reducing the stochastic input dimension in an optimal way. This low-dimensional, latent representation then serves as the input to the CGM, the output of which is basal for the probabilistic reconstruction of the fine scale output. Predictive uncertainty due to limited training data as well as due to the information loss unavoidably taking place during the coarse-graining step is accurately quantified by the fully probabilistic nature of the model. All parts of the model were trained simultaneously in a fully Bayesian way using Stochastic Variational Inference, without the need of any parameters to be specified by the analyst.

On the basis of the Stokes/Darcy similarity for flow through random porous media, we provided numerical evidence of the predictive capabilities under the aforementioned setting. Another important advantage of the proposed formulation is its ability to perform comparably well under extrapolative conditions. In the context of the models explored, we demonstrated this by carrying out predictions on data under a considerably different (or even randomized) boundary conditions than those used for training. Furthermore, we demonstrated one of the many possible application of the proposed model for Uncertainty Propagation problems where not only accurate estimates of output statistics but the confidence in these predictions is quantified. Finally, we outlined a method for automatically and adaptively refining the model and demonstrated its benefits in the context of the problems explored.

Applications of the proposed model can be found in any multi-query setting with expensive forward model evaluations, as is often the case in design-, optimization- or inverse problems. Several model extensions can be contemplated. For coarse-graining of the high-dimensional stochastic inputs, a deep learning framework could be adopted, potentially benefiting from recent advances in computer vision with e.g. convolutional neural networks (CNNs) or deep Gaussian process models. We emphasize however that significant regularization would be needed to make those work in the Small Data setting considered. Furthermore, one might think of a model that not only passes information through the CGM, but also entails ‘bypassing’ components, thereby effectively expanding the information bottleneck and allowing for more detailed reconstruction of fine scale output variability.

Appendix A Optimizing the variational distributions $Q$

To maximize the ELBO functional $\mathcal{F}(Q)$ in equation (25) w.r.t. the variational distributions $Q$ under the normalization constraint $\int Q_{\boldsymbol{\theta}_{k}}(\boldsymbol{\theta}_{k})d\boldsymbol{\theta}_{k}=1$, we maximize the functional $\mathcal{J}[Q;\boldsymbol{\xi}]$

[TABLE]

where $\boldsymbol{\theta}_{k}\subset\boldsymbol{\theta}=\left\{\left\{\boldsymbol{\lambda}_{c}^{(n)}\right\}_{n=1}^{N},\tilde{\boldsymbol{\theta}}_{c},\boldsymbol{\tau}_{c},\boldsymbol{\tau}_{cf},\boldsymbol{\gamma}\right\}$ is a subset of the model parameters $\boldsymbol{\theta}$ and $\boldsymbol{\xi}$ are Lagrangean multipliers. The first order variation $\delta\mathcal{J}$ w.r.t. $Q_{\boldsymbol{\theta}_{k}}$ is

[TABLE]

where $d\boldsymbol{\theta}_{-k}=\prod_{i\neq k}d\boldsymbol{\theta}_{i}$. As $\delta\mathcal{J}\stackrel{{\scriptstyle!}}{{=}}0$ for arbitrary $\delta Q_{\boldsymbol{\theta}_{k}}$,

[TABLE]

or

[TABLE]

where $\left<~{}.~{}\right>_{i\neq k}$ denotes the expected value w.r.t. $\prod_{i\neq k}\left(Q_{\boldsymbol{\theta}_{i}}(\boldsymbol{\theta}_{i})\right)$. The normalization constraints are recaptured by $\nabla_{\xi}\mathcal{J}\stackrel{{\scriptstyle!}}{{=}}0$ and determine the constant $\xi_{k}-1$ in Equation (66). Thus, it is found that

[TABLE]

A.1 Black-box variational inference for $Q_{\boldsymbol{\lambda}_{c}^{(n)}}$

We note that minimization of the KL divergence as given in equation (35) is equivalent to

[TABLE]

where $H\left(\tilde{Q}_{\boldsymbol{\lambda}_{c}^{(n)}}(\boldsymbol{\lambda}_{c}^{(n)})\right)$ is the Shannon entropy of $\tilde{Q}_{\boldsymbol{\lambda}_{c}^{(n)}}$. Expected values of terms in (68) depending on the CGM $\boldsymbol{u}_{c}(\boldsymbol{\lambda}_{c}^{(n)})$ are analytically intractable and gradients w.r.t. $\boldsymbol{\mu}_{\boldsymbol{\lambda}_{c}}^{(n)},\boldsymbol{\sigma}_{\boldsymbol{\lambda}_{c}}^{(n)}$ are therefore only accesible by noisy Monte Carlo estimates. To reduce the variance in those estimates, we apply the reparametrization trick Kingma and Welling (2013)

[TABLE]

and express Equation (68) as

[TABLE]

where ‘$\circ$’ denotes element-wise multiplication and we used the fact that

$H\left(\tilde{Q}_{\boldsymbol{\lambda}_{c}^{(n)}}(\boldsymbol{\lambda}_{c}^{(n)})\right)\vphantom{\sum_{m=1}^{\text{dim}(\boldsymbol{\lambda}_{c})}}\propto\sum_{m=1}^{\text{dim}(\boldsymbol{\lambda}_{c})}\log\sigma_{\boldsymbol{\lambda}_{c},m}^{(n)}$. Derivatives w.r.t. $\boldsymbol{\mu}_{\boldsymbol{\lambda}_{c}}^{(n)}$, $\boldsymbol{\sigma}_{\boldsymbol{\lambda}_{c}}^{(n)}$ can be computed using the chain rule as

[TABLE]

where $\partial\boldsymbol{\lambda}_{c}^{(n)}/\partial\boldsymbol{\mu}_{\boldsymbol{\lambda}_{c}}^{(n)}=\boldsymbol{I}$ and $\partial\boldsymbol{\lambda}_{c}^{(n)}/\partial\boldsymbol{\sigma}_{\boldsymbol{\lambda}_{c}}^{(n)}=\text{diag}[\boldsymbol{\epsilon}^{(n)}]$. To solve the optimization problem in (70), the above gradients are passed to the adaptive moment estimation optimizer (ADAM) using the default optimization parameters suggested in Kingma and Ba (2014).

Appendix B Feature functions

Table 3 shows a list of the 150 feature functions used in the numerical examples of section 3. Features 1–55 take the whole microstructure $\boldsymbol{\lambda}_{f}$ as input, features 56–150 use the subset $\boldsymbol{\lambda}_{f}^{(m)}$ pertaining to subdomain/cell $\Omega_{m}$ for which $\boldsymbol{K}_{m}(\boldsymbol{x},\boldsymbol{\lambda}_{c})=e^{\lambda_{c,m}}\boldsymbol{I}$.

Appendix C Proper orthogonal decomposition of output data $\boldsymbol{u}_{f}$

In section 3.2, it is mentioned that we observe very fast decay of POD eigenmodes for data with large variance in the pore space fraction $\text{vol}(\Omega_{f})/\text{vol}(\Omega)$, which translates to high $\sigma_{\text{ex}}^{2}$ as explained in section 3.1.1. We therefore plot the explained PCA variance, i.e. the cumulative sum of the $n$ largest variances of a PCA analysis performed on 2048 sample vectors $\mathcal{D}_{PCA}=\left\{\boldsymbol{u}_{f}(\boldsymbol{\lambda}_{f}^{(n)})\right\}_{n=1}^{2048}$ for the data used in sections 3.2 and 3.3 in figure 18. We note that this fast decay in POD eigenmodes observed in section 3.2 does not affect the effective input dimensionality of $\boldsymbol{\lambda}_{f}$. To show that, we discretize the 2048 microstructures from the distributions defined in sections 3.2, 3.3 to bitmap image vectors of $256\times 256=65,536$ dimensions and again plot the explained variance from the first $n$ components in the second row of figure 18. We observe that our claim of high-dimensional inputs is verified, as the explained variance is only slowly increasing with the number of POD components.

Appendix D Evidence Lower BOund (ELBO) $\mathcal{F}(Q)$

Using the expected values found in equations (28)–(32), omitting constants and canceling terms, we find

[TABLE]

as the most compact form of the evidence lower bound (ELBO) $\mathcal{F}(Q)$. We note that with the exception of the first term, the remaining ones are automatically separated into additive contributions from each cell $m$ or equivalently each $\lambda_{c,m}$. Hence, in order to arrive at the cell-scoring function $\mathcal{F}_{m}(Q)$ given in Equation (45), we pick out the terms in equation (74) which explicitly carry a cell index $m$. These are

[TABLE]

In the numerical experiments, we tie the precision parameters $\gamma_{jm}$ across different cells $m$ together i.e. $\gamma_{jm}=\gamma_{j}$ which allows sharing of information in terms of the microstructural features. As a result $b_{jm}=b_{j}$ and the corresponding contribution from each cell $m$ becomes:

[TABLE]