Bayesian Numerical Homogenization

Bayesian Numerical Analysis

This paper is inspired by a curious (and, perhaps, overlooked) link between Bayesian Inference and Numerical Analysis , known as Bayesian Numerical Analysis , and that can be traced back to Poincaré’s course on Probability Theory . We will recall Diaconis’ compelling example as an illustration of this link.

As a prototypical example, consider the numerical homogenization of the PDE

Recall that numerical homogenization concerns the approximation of the solution space of (1.1) with a finite-dimensional space. Although classical homogenization concepts might be present in some instances of this problem , one of the main objectives of numerical homogenization is to achieve a numerical approximation of the solution space of (1.1) with arbitrary rough coefficients, i.e., in particular, without the assumptions found in classical homogenization, such as scale separation, ergodicity at fine scales and $\epsilon$-sequences of operators. In this situation, piecewise linear finite-elements can perform arbitrarily badly and the numerical approximation of the solution space involves the identification of accurate basis elements adapted to the microstructure $a(x)$ .

As for the identification of quadrature rules in numerical analysis, the identification of accurate basis elements in numerical homogenization has been based on a difficult process of scientific investigation. Let us now turn our attention to the Bayesian approach to this problem. An immediate question is where do we place the prior? (1) If the prior is placed on $u$ then posterior values do not see (depend on) the microstructure. (2) If the prior is placed on $a$ then the microstructure becomes random whereas our purpose is the numerical homogenization of a given deterministic microstructure. Let us also note that the randomization of the microstructure, as investigated by Polynomial Chaos Approximation/Stochastic Expansion methods , does not lead to the simplification seen after homogenization but to increased complexity with the dimension of input stochastic variables (although Stochastic Expansion methods have been used successfully to beat Monte-Carlo sampling they do not lead to averaging results seen in homogenization). (3) If the prior is placed on $g$ then the noise propagates through the microstructure and the posterior value of $u$ contains that information.

This observation motivates us to place the prior on the source term $g$ in (1.1), e.g., replace it by white noise (i.e. a centered Gaussian field $\xi(x)$ on $\Omega$ with covariance function $\delta(x-y)$) and consider the stochastic PDE

where the functions $\phi_{i}$ are Rough Polyharmonic Splines (RPS) which have been identified as accurate basis elements for the numerical homogenization of (1.1) having noteworthy variational, optimal recovery and localization properties. The discovery of these Rough Polyharmonic Splines has required a significant amount of work and trial and errors but here, they are identified after a single step of Bayesian conditioning.

This observation motivates us to investigate what the same process of Bayesian conditioning would give under different priors and under other observations than the values of $u$ at individual points (we will consider data formed by the values of a finite number of linear functions of $u$). In particular, we will use this link between Bayesian Inference and Numerical Homogenization to identify bases for the numerical homogenization of arbitrary linear integro-differential equations. Our purpose is to show that this link is generic and could in principle be used, beyond numerical homogenization, as a guiding principle for the coarse-graining of multi-scale systems. The Bayesian approach to this problem is to (1) Put a prior on the degrees of freedom of the system (2) Select a finite number of coarse variables (3) Compute the posterior value of the state of the system conditioned on the coarse variables.

General setup

Let $\mathcal{L}$ and $\mathcal{B}$ be linear integro-differential operators on $\Omega$ and $\partial\Omega$ such that (1) $(\mathcal{L},\mathcal{B}):\mathcal{H}(\Omega)\rightarrow\mathcal{H}_{\mathcal{L}}(\Omega)\times\mathcal{H}_{\mathcal{B}}(\partial\Omega)$, where $\mathcal{H}(\Omega)$, $\mathcal{H}_{\mathcal{L}}(\Omega)$ and $\mathcal{H}_{\mathcal{B}}(\partial\Omega)$ are Hilbert spaces of Generalized functions on $\Omega$ and $\partial\Omega$ (2) $\mathcal{H}_{\mathcal{L}}(\Omega)$ contains $L^{2}(\Omega)$ and $\mathcal{H}(\Omega)$ is contained in $L^{2}(\Omega)$.

Consider the integro-differential equation

As with (1.1) the numerical homogenization of (2.1) will require the assumption that $g$ belongs to a strict subspace of $\mathcal{H}_{\mathcal{L}}(\Omega)$.

We will assume that $\mathcal{L}$ and $\mathcal{B}$ are such that (2.1) (1) admits a unique solution in $\mathcal{H}(\Omega)$ (2) and a Green’s function $G$. Recall that $G$ is defined as the solution of

where $\delta(\cdot-y)$ is the Delta mass of dirac at the point $y$.

Note that for the prototypical example (1.1) we have

Our purpose is to identify a good basis for the numerical homogenization or coarse-graining of (2.1).

Bayesian Numerical Homogenization

Our Bayesian approach to the numerical homogenization of (2.1) is to replace the source term $g$ by a Gaussian field $\xi$. More precisely we introduce $\xi$, a centered Gaussian field on $\Omega$ with covariance function

and consider the stochastic integro-differential equation

Write $(\mathcal{L}^{*},\mathcal{B}^{*})$ the adjoint of $(\mathcal{L},\mathcal{B})$ with respect to the (scalar) product defined on $\mathcal{H}(\Omega)$ by \big{\langle}u,v\big{\rangle}_{L^{2}}:=\int_{\Omega}u(x)v(x)\,dx. Observe that $G(y,x)$ (the transpose of $G(x,y)$ with respect to the scalar product \big{\langle}\cdot,\cdot\big{\rangle}_{L^{2}}) is the Green’s function of $(\mathcal{L}^{*},\mathcal{B}^{*})$ (the complex conjugation of the Green’s function is not required to define its adjoint because the scalar product is bilinear and not sesquilinear). Observe that if $\xi$ is white noise (i.e. $\Lambda(x-y)=\delta(x-y)$) then

which is the Kernel of $\mathcal{L}^{*}\mathcal{L}$, i.e., $\mathcal{L}^{*}\mathcal{L}\Gamma(x,y)=\delta(x-y)$.

Since $\mathcal{L}$ and $\mathcal{B}$ are linear operators, $u$ is a linear function of $\xi$ and is therefore a Gaussian field. Moreover its covariance function is given by

Beyond Bayesian Homogenization, equations with random right hand side can also be of interest in practical applications, for instance in the modeling of the electrostatics in nanoscale field-effect sensors, where fluctuations arise from random charge concentrations .

We will show that the choice of the noise $\Lambda$ can be determined by the regularity of the source term $g$ in the right hand side of (2.1). More precisely if $\xi$ is white noise ($\Lambda(x,y)=\delta(x-y)$) then the resulting accuracy estimates will be obtained under the assumption that $g\in L^{2}(\Omega)$ and as a function of $\|g\|_{L^{2}(\Omega)}$.

If $\xi$ is not white noise (i.e. if its covariance function is not $\delta(x-y)$) then we assume that there exists two linear integro-differential operators $\mathcal{L}_{\Lambda}$ and $\mathcal{B}_{\Lambda}$ such that $\xi$ is the stochastic solution of the following equation with white noise $\xi^{\prime}$ as the source term:

In what follows, if $\xi$ is not white noise then we assume it to be obtained as in (3.6) and the resulting accuracy estimates will be obtained under the assumption that $\mathcal{L}_{\Lambda}g\in L^{2}(\Omega)$ and as a function of $\|\mathcal{L}_{\Lambda}g\|_{L^{2}(\Omega)}$. A prototypical example corresponds to the situation where $\xi$ is obtained as the regularization of white noise via a power of the Laplace Dirichlet operator on $\Omega$ and this allows us to identify optimal recovery bases under the assumption that $g\in H^{s}(\Omega)$ with $s\geq 0$ or $s<0$.

2 Identification of basis elements via conditioning

Let $N$ be a strictly positive integer. Our Bayesian approach is based on the conditioning of the solution of (3.2) posterior to the observation of $N$ linear functions of $u(x)$, expressed as

where $\psi_{1},\ldots,\psi_{N}$ are $N$ linearly independent generalized functions (distributions) on $\Omega$ such that for all $i$

Examples of $\psi_{i}$ include masses of Dirac ($\psi_{i}(x)=\delta(x-x_{i})$), indicator functions of subsets of $\Omega$ and elements of $L^{1}(\Omega)$. Let $\Theta$ be the $N\times N$ symmetric matrix defined by

Note that (3.8) implies that if $u$ is the solution of (3.2) then

is a well defined center Gaussian random vector with covariance matrix $\Theta$.

We will from now on assume that the covariance function (3.1) is not degenerate in the sense that for $f\in\mathcal{H}(\Omega)$,

is zero if and only if $f$ is the null function. Note that if $\xi$ is obtained via (3.6) then $\|f\|_{\Lambda}^{2}=\|\mathcal{L}_{\Lambda}^{-1}f\|_{L^{2}(\Omega)}^{2}$ (writing $\mathcal{L}_{\Lambda}^{-1}f$ the solution of $\mathcal{L}_{\Lambda}u=f$ in $\Omega$ with $\mathcal{B}_{\Lambda}u=0$ on $\partial\Omega$) and the non-degeneracy of $\Lambda$ is equivalent to that of the operator $\mathcal{L}_{\Lambda}$.

Our motivation for using Gaussian noise in (3.2) lies in the fact that for Gaussian fields, conditional expected values can be computed via linear projection. Henceforth our approach is also akin to Gaussian filtering for numerical homogenization and the following Theorem shows that this approach allows for the identification of a (projection) basis $\phi_{i}$.

Let $u$ be the solution of (3.2) and $\Psi$ defined by (3.10), then

Furthermore, $u(x)$ conditioned on the value of $\Psi$ is a Gaussian random variable with mean (3.16) and variance

Since $u$ and $\Psi$ belong to the same Gaussian space, it follows that $u_{\Psi}$ is a linear function of $\Psi$ obtained by minimizing the mean squared error

If $\mathcal{L}$ and $\mathcal{B}$ correspond to the prototypical example (1.1) (see also Example 2.1), if $\xi$ is white noise (i.e. if its covariance matrix is $\Lambda(x,y)=\delta(x-y)$), and if the observable functions are masses of Diracs at points $x_{i}\in\Omega$ (and $d\leq 3$ which is required for (3.8)), then Theorem 3.5 implies (1.3) and the basis elements $\phi_{i}$ are the RPS elements of which are a generalization of Polyharmonic Splines to PDEs with rough coefficients. Recall that Polyharmonic splines can be traced back to the seminal work of Harder and Desmarais and Duchon .

Note also that according to Theorem (3.5) the process of Bayesian conditioning gives us the whole posterior distribution of $u(x)$ and not only its (conditional) expected value. In particular, the distribution of $u(x)$ conditioned on $u(x_{1}),\ldots,u(x_{N})$ is a Gaussian random variable with mean (1.3) and variance

and this observation can be used to compute the probability of deviation of the RPS interpolation from $u(x)$ by a given margin and guide the addition of interpolation points (note that $\sigma^{2}(x)=0$ at the interpolation points $x_{1},\ldots,x_{N}$).

We will show in Theorem 5.1 that $\sigma(x)$ also controls the pointwise error between the solution of the original integro-differential equation (2.1) and the approximation $\sum_{i=1}^{N}\phi_{i}(x)\int_{\Omega}u(y)\psi_{i}(y)\,dy$.

Variational properties of basis elements

In this section we will show that as for RPS , the basis elements $\phi_{i}$ from Bayesian Inference have remarkable variational and optimal recovery properties that can be used (1) for their practical computation (2) for the derivation of accuracy estimates.

In this subsection we will assume that $\xi$ is white noise (i.e. $\Lambda(x,y)=\delta(x-y)$). Define

and let \big{\langle}\cdot,\cdot\big{\rangle} be the (scalar) product on $V$ defined by: for $u,v\in V$,

Note in particular that \big{\langle}v,v\big{\rangle}=0 if and only if $v=0$ and we write

the corresponding norm (note that $\|v\|_{V}$ is a norm on $V$ because $\|v\|_{V}=0$ and $v\in V$ imply $\mathcal{L}v=0$ in $\Omega$ and $\mathcal{B}v=0$ on $\partial\Omega$ which leads to $v=0$ by the non-degeneracy of the operator $\mathcal{L}$).

If $\Gamma(x,x)<\infty$ then for $v\in V$ and $x\in\Omega$

and the space $V$ with the reproducing Kernel $\Gamma(x,y)$ forms a Reproducing Kernel Hilbert Space. In particular, for all $v\in V$

Theorem 4.1 is a direct consequence of the fact that

and (by Cauchy-Schwartz inequality and \big{\langle}\Gamma(\cdot,x),\int_{\Omega}\Gamma(\cdot,x)\big{\rangle}=\Gamma(x,x))

and consider the following optimization problem over $V_{i}$:

$V_{i}$ is a non-empty closed affine subspace of $V$. Problem (4.9) is a strictly convex quadratic optimization problem over $V_{i}$. The unique minimizer of (4.9) is $\phi_{i}$ as defined by (3.18).

Let us first prove that $\phi_{i}\in V_{i}$. Let

First observe that for all $i\in\{1,\ldots,N\}$,

and $\mathcal{B}\theta_{i}(x)=0$ on $\partial\Omega$. Noting that \big{\|}\mathcal{L}\theta_{i}\big{\|}_{L^{2}(\Omega)}^{2}=\Theta_{i,i} we deduce from (3.8) that $\theta_{i}\in V$. We conclude from (3.18) and Lemma 3.4 that $\phi_{i}\in V$. Now observe that (3.9) implies that

where $\delta_{i,i}=1$ and $\delta_{i,j}=0$ for $j\not=i$. We conclude that $\phi_{i}\in V_{i}$ which implies that $V_{i}$ is non empty (it is easy to check that it is a closed affine sub-space of $V$).

Now let us prove that problem (4.9) is a strictly convex optimization problem over $V_{i}$. Let $v,w\in V_{i}$ such that $v\not=w$. Write for $\lambda\in$,

and we need to show that $f(\lambda)$ is a strictly convex function. Observing that

and noting that \big{\langle}v-w,v-w\big{\rangle}>0 (otherwise one would have $v=w$) we deduce that $f$ is strictly convex in $\lambda$. We conclude that (see, for example, [29, pp. 35, Proposition 1.2]) that Problem (4.9) is a strictly convex optimization problem over $V_{i}$ and that it admits a unique minimizer in $V_{i}$. We will postpone the proof of the fact that $\phi_{i}$ is the minimizer of (4.9) to the proof of Theorem 4.6. ∎

It is important to note that in practical (numerical) applications each element $\phi_{i}$ would be obtained by solving the quadratic optimization problem (4.9) rather than through the representation formula (3.18) because the identification of $\Gamma$ in (3.18) is more expensive than solving the linear systems associated with (4.9) (inverting a matrix is more expensive than solving a linear system). Note also that, if $u$ is the (stochastic) solution of (3.2), then $\phi_{i}$ is also equal to the expected value of $u(x)$ conditioned on $\int_{\Omega}u(x)\psi_{i}(x)=1$ and $\int_{\Omega}u(x)\psi_{j}(x)=0$ for $j\not=i$, i.e.

A simple calculation allows us to show that $\phi_{i}$ is also the solution of the following nested equations

Another simple calculation allows us to show that $\phi_{i}$ is also the solution of the following nested equations

Write $V_{0}$ the subset of $V$ defined by

The basis $\phi_{i}$ is orthorgonal to $V_{0}$ with respect to the product \big{\langle}\cdot,\cdot\big{\rangle}, i.e.

$\sum_{i=1}^{N}w_{i}\phi_{i}$ is the unique minimizer of \big{\langle}v,v\big{\rangle} over all $v\in V$ such that $\int_{\Omega}v(x)\psi_{i}(x)\,dx=w_{i}$.

For all $i\in\{1,\ldots,N\}$ and for all $v\in V$,

Theorem 4.6 and its proof is analogous to the optimal property of strictly conditionally positive definite kernels when used as interpolant solutions of the optimal recovery problem .

It follows that $\sum_{i=1}^{N}w_{i}\phi_{i}$ is the unique minimizer of \big{\langle}v,v\big{\rangle} over all $v\in V$ such that $\int_{\Omega}v(x)\psi_{i}(x)\,dx=w_{i}$. Note that this also implies that $\phi_{i}$ is the minimizer of (4.9). ∎

2 Non-white Gaussian noise

If $\xi$ is not white noise (i.e. $\Lambda(x,y)\not=\delta(x-y)$) then Theorem 4.1,Theorem 4.6 and Proposition 4.2 remain true provided that the definitions of the space $V$ and scalar product \big{\langle}\cdot,\cdot\big{\rangle} are changed to

where $\mathcal{L}_{\Lambda}$ and $\mathcal{B}_{\Lambda}$ are defined in (3.6).

Assume that $\Gamma(x,x)<\infty$. Let $v\in V$. It holds true that for $x\in\Omega$

where $\sigma^{2}(x)$ is the variance of $u(x)$ (solution of (3.2)) conditioned on $\int_{\Omega}u(y)\psi_{1}(y)\,dy,\ldots,\int_{\Omega}u(y)\psi_{N}(y)\,dy$ as defined by (3.19). In particular if $u$ is the solution of the original integro-differential equation (2.1), then

if $\phi_{i},\sigma$ are derived from white noise, and

if $\phi_{i},\sigma$ are derived from the noise with covariance function $\Lambda$ described in (3.6).

Let $v\in V$ and $x\in\Omega$. Using the reproducing kernel property of Theorem 4.1 we obtain that

Therefore, using Cauchy-Schwartz inequality

We conclude by expanding the right hand side of (5.5) and the definition $\phi_{i}(x)=\sum_{j=1}^{N}\Theta^{-1}_{i,j}\int_{\Omega}\Gamma(x,y)\psi_{i}(y)\,dy$. ∎

$\sigma^{2}(x)$ is also known as the Power function in radial basis function interpolation . The proof of Theorem 5.1 is similar to the one used to derive local error estimates for radial basis function interpolation of scattered data (see in which $\sigma^{2}(x)$ was referred to as the Kriging function, a terminology coming from geostatistics ).

2 ℋ​(Ω)ℋΩ\mathcal{H}(\Omega)-norm estimates

Let $V_{0}$ be the subset of $V$ defined by (4.20). Write

where $\|.\|_{\mathcal{H}(\Omega)}$ is the natural norm associated with the space on which the operator $\mathcal{L}$ is defined.

and $\rho(V_{0})$ is the smallest constant for which (5.7) holds for all $v\in V$.

Write v_{\Psi}(x):=\sum_{i=1}^{N}\phi_{i}(x)\big{(}\int_{\Omega}v(y)\psi_{i}(y)\,dy\big{)}

Observing that $v-v_{\Psi}$ belongs to $V_{0}$ implies that

Observe that Theorem 5.3 implies that if $u$ is the solution of the original integro-differential equation (2.1) and $\phi_{i},\sigma$ are derived from white noise, then

Similarly, if $\phi_{i},\sigma$ are derived from the noise with covariance function $\Lambda$ described in (3.6), then

If $\mathcal{L}$ and $\mathcal{B}$ correspond to the prototypical example (1.1) (Example 2.1), if $\xi$ is white noise, and if the observable functions are masses of Diracs at points $x_{i}\in\Omega$ (and $d\leq 3$), then ,

where $C$ depends only on $\lambda_{\min}(a),\lambda_{\max}(a)$ and where $\lambda_{\max}(a):=\sup_{x\in\Omega,l\not=0}l^{T}a(x)l/|l|^{2}$, $\lambda_{\min}(a):=\inf_{x\in\Omega,l\not=0}l^{T}a(x)l/|l|^{2}$ and $H$ is the mesh-norm

Let us also recall that the proof of (5.13) is based on the following Poincaré inequality (Lemma 3.1 of )

([60, Lemma 3.1 ]) Let $d\leq 3$ and $B_{1}$ be the open ball of center and radius $1$. There exists a finite strictly positive constant $C_{\lambda_{\min}(a),\lambda_{\max}(a)}$ such that for all $v\in\mathcal{H}^{1}(B_{1})$ such that $\operatorname{div}(a\nabla v)\in L^{2}(B_{1})$ it holds true that

We will recall the proof of this lemma (as presented in [60, Lemma 3.1 ]) for the sake of completeness. The proof is per absurdum. Note that since $d\leq 3$ the assumptions $v\in\mathcal{H}^{1}(B_{1})$ and $\operatorname{div}(a\nabla v)\in L^{2}(B_{1})$ imply the Hölder continuity of $v$ in $B_{1}$. Assume that (5.16) does not hold. Then there exists a sequence $v_{n}$ and a sequence $a_{n}^{\prime}$ whose maximum and minimum eigenvalues are uniformly bounded by $\lambda_{\min}(a)$ and $\lambda_{\max}(a)$ (we need to introduce that sequence because we want the constant in (5.16) to depend only $d,\lambda_{\min}(a),\lambda_{\max}(a)$) such that

Letting $w_{n}=\frac{v_{n}-v_{n}(0)}{\|v_{n}-v_{n}(0)\|_{L^{2}(B_{1})}}$ we obtain that $w_{n}(0)=0$, $\|w_{n}\|_{L^{2}(B_{1})}=1$ and

it follows that there exists a subsequence $w_{n_{j}}$ and a $w\in\mathcal{H}^{1}(B_{1})$ such that $w_{n_{j}}\rightharpoonup w$ weakly in $\mathcal{H}^{1}(B_{1})$ and $\nabla w_{n_{j}}\rightharpoonup\nabla w$ weakly in $L^{2}(B_{1})$. Using $\|\nabla w_{n}\|_{L^{2}(B_{1})}\leq 1/n$ we deduce that $\nabla w=0$ which implies that $w$ is a constant in $B_{1}$. Since by the Rellich–-Kondrachov theorem the embedding $\mathcal{H}^{1}(B_{1})\subset L^{2}(B_{1})$ is compact it follows from (5.19) that $w_{n_{j}}\rightarrow w$ strongly in $L^{2}(B_{1})$ which (using $\|w_{n}\|_{L^{2}(B_{1})}=1$) implies that $\|w\|_{L^{2}(B_{1})}=1$. Now (5.19) together with the fact that \big{\|}\operatorname{div}(a_{n}^{\prime}\nabla w_{n})\big{\|}_{L^{2}(B_{1})}^{2} is uniformly bounded and that $d\leq 3$ implies that $w_{n}$ is uniformly Hölder continuous on $B(0,\frac{1}{2})$ (see for instance ). This implies that $w$ is continuous in $B(0,\frac{1}{2})$ and that $w(0)=0$. This contradicts the fact that $w$ is a constant in $B_{1}$ with $\|w\|_{L^{2}(B_{1})}=1$. ∎

If $\mathcal{L}$ and $\mathcal{B}$ correspond to the prototypical example (1.1) (Example 2.1), if $\xi$ is white noise, and if the observable functions are indicator functions of Voronoï cells around points in $x_{i}\in\Omega$ or of tetrahedra of a regular tessellation of the points $x_{i}\in\Omega$ then (5.13) remains valid as a simple consequence of localized Poincaré inequalities. Indeed for $v\in V_{0}$, writing $C_{i}$ the Voronoï cells at the points $x_{i}\in\Omega$, we have (assuming $\Omega$ is the union of those Voronoï cells)

and we conclude by applying Poincaré’s inequality to the $L^{2}$-norm of $v$ within each cell $C_{i}$, i.e.

We will give the last example as a theorem.

Let $\mathcal{L}$ and $\mathcal{B}$ be as in the prototypical example (1.1) (Example 2.1) and let $\xi$ be white noise. Let $\psi_{1},\ldots,\psi_{N}$ be linearly independent generalized probability densities on $\Omega$ with (possibly overlapping) support $\operatorname{support}(\psi_{i})$. Define

where $C$ depends only on $\lambda_{\min}(a)$ and $\lambda_{\max}(a)$. Henceforth, for $u\in V$

Observe that if for all $i$ the support of $\psi_{i}$ is contained in a ball of center $x_{i}$ and radius $H^{\prime}$, then

in particular if the points $x_{i}$ have mesh norm $H^{\prime\prime}$ (see (5.14)) then $H\leq H^{\prime}+H^{\prime\prime}$.

The proof of (5.23) is simply based on the observation that if $v\in V_{0}$ then (since $\int_{\Omega}v(x)\psi_{i}(x)\,dx=0$) there exists $N$ points $y_{1},\ldots,y_{N}$ such that $v(y_{i})=0$ and the mesh norm of those points is bounded by $H$. Therefore we can apply the result of Example 5.1. ∎

Pseudo-algorithm

A simple pseudo-algorithmic description of the proposed framework for the numerical homogenization of (2.1) is as follows:

Select $N$ linearly independent (measurement) functions $\psi_{1},\ldots,\psi_{N}$ in $L^{2}(\Omega)$.

Let $\xi$ in (3.2) be a Gaussian field of mean and covariance function $\Lambda(x,y)$ (assumed to be non-degenerate, i.e. such that there exists an inverse covariance function $\Lambda^{-1}(x,y)$ with $\int_{\Omega^{2}}\Lambda(x,y)\Lambda^{-1}(y,z)\,dy=\delta(x-z)$).

The basis functions $\phi_{1},\ldots,\phi_{N}$ for the numerical homogenization of (2.1) are identified as (writing $u$ the solution of (3.2) and $\delta_{i,j}=1$ if $i=j$ and $\delta_{i,j}=0$ if $i\not=j$) the deterministic functions

Each $\phi_{i}$ can also be identified as the unique minimizer of

Under appropriate choice of the measurement functions $\psi_{i}$ and the covariance function $\Lambda(x,y)$, the basis functions $\phi_{i}$ can be computed by localizing the optimization problems (6.2) to subdomains of $\Omega$.

Statistical Decision Theory and Practical Applications

Another motivation for exploring Bayesian approximations of the solution space, lies in the decision theory/game theory approach to numerical homogenization. In this approach one looks at the numerical homogenization problem (1.1) as a repeated game where player B chooses a function $\theta$ of the linear measurements (data) $\int_{\Omega}u(x)\psi_{1}(x)\,dx,$ $\ldots,$ $\int_{\Omega}u(x)\psi_{N}(x)\,dx$ and player A chooses a source term $g$ in the unit ball of $L^{2}(\Omega)$. These two choices combine and form an error term

Player’s B objective is to minimize the error (7.1) while player’s A objective is to maximize it. A surprising result stemming from a generalization of Wald’s Decision Theory and Von Neumann’s Game Theory is that, although such games are deterministic, under weak regularity conditions, the optimal strategy for player $A$ is to play at random by placing an optimal probability distribution $\pi_{A}$ on the set of candidates for $g$ and, similarly, the best strategy for player $B$ is to assume that player A is playing at random and to use a function $\theta$ living in the Bayesian class (obtained by placing a prior $\pi_{B}$ on the set of candidates for $g$ and conditioning with respect to the measurements $\int_{\Omega}u(x)\psi_{i}(x)\,dx$).

Although the estimator employed by player B may be called Bayesian, the game described here is not (i.e. the choice of player A might be distinct from that of player B) and player B must solve a min max optimization problem over $\pi_{A}$ and $\pi_{B}$ to identify an optimal prior distribution for the Bayesian estimator (a careful choice of the prior also appears to be important due to the possible high sensitivity of posterior distributions ).

We refer to for (1) the complete description of the generalization of the Bayesian framework described here to the decision theory/information game formulation (described above) (2) practical (including numerical) applications of that generalized framework to the problems of finding numerical homogenization bases and fast solvers for (1.1). In that generalization, optimal numerical homogenization bases functions are obtained by selecting the prior distribution of $\xi$ (in (1.2)) to be that of a Gaussian field with mean zero and covariance function the operator (1.1) (i.e. such that for $f\in H^{1}_{0}(\Omega)$, $\int_{\Omega}f(x)\xi(x)\,dx$ is a Gaussian random variable of mean zero and variance $\int_{\Omega}(\nabla f(x))^{T}a(x)\nabla f(x)\,dx$). In particular shows how the identification of an optimal distribution for $\xi$ (in the Gaussian class) leads to the (automated) discovery of multigrid and multiresolution solvers for PDEs with rough coefficients.

The author gratefully acknowledges this work supported by the Air Force Office of Scientific Research under Award Number FA9550-12-1-0389 (Scientific Computation of Optimal Statistical Estimators) and the U.S. Department of Energy Office of Science, Office of Advanced Scientific Computing Research, through the Exascale Co-Design Center for Materials in Extreme Environments (ExMatEx, LANL Contract No DE-AC52-06NA25396, Caltech Subcontract Number 273448). The author also thanks Dongbin Xiu, Lei Zhang and Guillaume Bal for stimulating discussions and Leonid Berlyand for comments on the manuscript. The author also thanks two anonymous referees for valuable comments and suggestions.

References