Multilevel Sampling

Collaborator: Max Gunzburger (Florida State)

Once the probability distributions of the uncertain parameters in a PDE are known, they can be used to perform stochastic simulations. A statistical description of the model output \(u(x,\omega)\), or some related physical quantity, such as a local spatial average or a boundary flux, usually takes the form of a set of statistical quantities of interest \(Q\), such as measures of central tendency, spread, and correlation, as well as exceedence probabilities and confidence intervals. The vast majority of such statistical quantities of interest \(Q\) can be written in the form of a stochastic integral, \[ Q = E[v] := \int_\Gamma v(x,y)\rho(y)\;dy, \] where \(v = G(u)\) for some mapping \(G\), and \(Y\) is a random vector with density function \(\rho:\Gamma\rightarrow [0,\infty)\). Stochastic sampling methods, such as Monte Carlo-, quasi Monte Carlo, and sparse grid stochastic collocation, are arguably the most directly parallelizable and least intrusive means for estimating \(Q\), and are computed as a weighted sum of the form \(I_\eta[v_h] := \sum_{i=1}^\eta w_i v_h(y_i)\). Here, \(v_h(y)\) is the spatial approximation of the sample path \(v(y)\) whose accuracy depends on an underlying discretization parameter \(h\). Since the computation of each sample path requires the numerical solution of a PDE, the overall computation can carry a considerable cost, especially if reliable statistics are sought, requiring a potentially large number of samples computed at high levels of spatial accuracy.

Multilevel sampling methods improve upon the efficiency of traditional sampling schemes without compromising on accuracy and parallelizability, by dynamically incorporating the model’s physical discretization into the sampling procedure through the use of a hierarchy of discretization models. Let \(\{h_\ell\}_{\ell=0}^L\) be a sequence of spatial discretization parameters giving an increasing level of accuracy, where \(h_L\) is chosen to achieve a desired level of spatial accuracy. We can now write sample paths of this fine-scale approximation \(v_{h_L}\) of \(v\) as the sum of an initial coarse scale approximation plus a series of correction terms, i.e. \[ v_{h_{L}} = v_{h_0} + \sum_{\ell=1}^L (v_{h_{\ell}}-v_{h_{\ell-1}}) =: \sum_{\ell=0}^L \triangle v_\ell, \ \text{ where }\ \triangle v_\ell := \left\{\begin{array}{ll} v_{h_0} & \text{if \(\ell = 0\)}\\ v_{h_{\ell}}-v_{h_{\ell-1}} & \text{if \(\ell > 0\)}\end{array}\right. . \] The multilevel estimate of \(Q\) is then given by the sum \begin{equation} \widehat{Q}_{\{h_\ell\},\{\eta_\ell\}}^{\mathrm{ML}} := \sum_{\ell=0}^L I_{\eta_\ell}[\triangle v_{h_\ell}], \label{eqn:qml} \end{equation} where the sample sizes \(\eta_\ell\) may be chosen separately for each spatial refinement level \(\ell\). In particular, we can choose \(\eta_\ell\) so as to minimize the overall computational cost, subject to the error satisfying \(\displaystyle\|Q-\widehat{Q}_{\{h_\ell\},\{\eta_\ell\}}^{\mathrm{ML}}\|<\varepsilon\) for a given tolerance \(\varepsilon>0\). In principle, we would expect the solution of this optimal allocation sub-problem to prescribe drawing the bulk of samples cheaply at coarser refinement levels, while sampling more sparingly at finer ones. For Monte Carlo sampling, the sampling error is fairly simple and the sub-problem can be solved explicitly. This is not necessarily the case for other, more complicated sampling schemes. In VanWyk2014, we develop and analyze a multilevel sampling scheme for stochastic collocation methods.

These methods are often more efficient, depending on factors such as the stochastic dimension, the complexity of the underlying parameter space, and the regularity of \(v\) as a function of \(y\). Several difficulties arise in solving the allocation sub-problem for sparse grids: The convergence rate of the associated sampling error depends on the smoothness of \(v_h\) as a function of \(y\) and is therefore seldom known a priori. Moreover, since the sample sizes are related to fixed grids, there are additional constraints on what sizes are admissible. We nevertheless derive explicit formulae for \(\eta_0,...,\eta_L\) that give considerable cost savings, regardless of the rate at which the cost of generating each sample path increases as \(h_\ell \rightarrow 0\). Such savings are possible due to the fact that \(\triangle v_\ell \rightarrow 0\) as \(h_\ell\rightarrow 0\), rendering the absolute error estimates of the fine-scale corrections in \eqref{eqn:qml} less significant to the overall error than coarse scale ones, thus allowing for smaller sample sizes.