A non-adapted sparse approximation of PDEs with stochastic inputs

Introduction

Realistic analysis and design of complex engineering systems require not only a fine understanding and modeling of the underlying physics and their interactions, but also a significant recognition of intrinsic uncertainties and their influences on the quantities of interest. Uncertainty Quantification (UQ) is an emerging discipline that aims at addressing the latter issue; it aims at meaningful characterization of uncertainties in the physical models from the available measurements and efficient propagation of these uncertainties for a quantitative validation of model predictions.

Despite recent growing interests in UQ of complex systems, it remains a grand challenge to efficiently propagate uncertainties through systems characterized by a large number of uncertain sources where the so-called curse-of-dimensionality is yet an unsolved problem. Additionally, development of non-intrusive uncertainty propagation techniques is of essence as the analysis of complex multi-disciplinary systems often requires the use of sophisticated coupled deterministic solvers in which one cannot readily intrude to set up the necessary propagation infrastructure.

Sampling methods such as the Monte Carlo simulation and its several variants had been utilized for a long time as the primary scheme for uncertainty propagation. However, it is well understood that these methods are generally inefficient for large-scale systems due to their slow rate of convergence. There has been an increasing recent interest in developing alternative numerical methods that are more efficient than the Monte Carlo techniques. Most notably, the stochastic Galerkin schemes using polynomial chaos (PC) bases (Ghanem03, ; Deb01, ; Xiu02, ; Babuska04, ; Wan05, ) have been successfully applied to a variety of engineering problems and are extremely useful when the number of uncertain parameters is not large. In their original form, the stochastic Galerkin schemes are intrusive, as one has to modify the deterministic solvers for their implementation. Stochastic collocation schemes (Tatang95, ; Mathelin03, ; Xiu05a, ; Babuska07a, ; Nobile08a, ) belong to a different class of methods that rely upon (isotropic) sparse grid integration/interpolation in the stochastic space of the problem to reduce the curse-of-dimensionality associated with the conventional tensor-product integration/interpolation rules. As their construction is primarily based on the input parameter space, the computational cost of both stochastic Galerkin and collocation techniques increases rapidly for large number of independent input uncertainties.

More recently, efforts have been made to construct solution-adaptive uncertainty propagation techniques that exploit any structures in the solution to decrease the computational cost. Among them are the multi-scale model reduction of (Doostan07, ) and the sparse decomposition of (Todor07a, ; Bieri09a, ; Bieri09b, ; Bieri09c, ; Blatman10, ) for the stochastic Galerkin technique, anisotropic and adaptive sparse grids of (Nobile08b, ; Ma09a, ) for the stochastic collocation scheme, and low-rank solution approximations of (Nouy07, ; Nouy08, ; Doostan09, ).

In the present study, we are interested in cases where the quantity of interest is sparse at the stochastic level, i.e., it can be accurately represented with only few terms when linearly expanded into a stochastic, e.g., polynomial chaos, basis. Interestingly, sparsity is salient in the analysis of high-dimensional problems where the number of energetic basis functions (those with large coefficients) is small relative to the cardinality of the full basis. For instance, it has been shown in (Todor07a, ; Bieri09a, ) that, under some mild conditions, solutions to linear elliptic stochastic PDEs with high-dimensional random coefficients admit sparse representations with respect to the PC basis. Consequently, an approach based on a zero-dimensional algebraic stochastic problem has been proposed in (Bieri09a, ) to detect the sparsity pattern, which then guides the stochastic Galerkin analysis of the original problem. Moreover, a “quasi”-best $N$-term approximation for a class of elliptic stochastic PDEs has been proposed in (Bieri09c, ).

In this work, using concentration of measure inequalities and compressive sampling techniques, we derive a method for PC expansion of sparse solutions to stochastic PDEs. The proposed method is

Non-intrusive: it is based on the direct random sampling of the PDE solutions. This sampling can be done by using any legacy code for the deterministic problem as a black box.

Non-adapted: it does not tailor the sampling process to identify the important dimensions at the stochastic level

Provably convergent: we obtain probabilistic bounds on the approximation error proving the stability and convergence of the method.

Well-suited to problems with high-dimensional random inputs.

Compressive sampling is an emerging direction in signal processing that aims at recovering sparse signals accurately (or even exactly) from a small number of their random projections (Chen98, ; Chen01a, ; Candes06a, ; Donoho06b, ; Candes06b, ; Candes06c, ; Candes07a, ; Cohen09a, ; Bruckstein09, ). A sparse signal is simply a signal that has only few significant coefficients when linearly expanded into a basis, e.g., $\{\psi_{\bm{\alpha}}\}$.

For sufficiently sparse signals, the number of samples needed for a successful recovery is typically less than what is required by the Shannon-Nyquist sampling principle. Generally speaking, a successful signal reconstruction by compressive sampling is conditioned upon:

Incoherent random projections of the signal.

A square-measurable stochastic function $u(\omega)$, defined on a suitable probability space $(\Omega,\mathcal{F},\mathcal{P})$ can be expanded into a mean-squared convergent series of the chaos polynomial bases, i.e., $u(\omega)\approx\sum_{\bm{\alpha}}c_{\bm{\alpha}}\psi_{\bm{\alpha}}(\omega)$, with some cardinality $P$. The stochastic function $u(\omega)$ is then sparse in PC basis $\{\psi_{\bm{\alpha}}\}$, if only a small fraction of coefficients $c_{\bm{\alpha}}$ are significant. In this case, under certain conditions, the sparse PC coefficients $\bm{c}$ may be computed accurately and robustly using only $N\ll P$ random samples of $u(\omega)$ via compressive sampling. Given $N$ random samples of $u(\omega)$, compressive sampling aims at finding the sparsest (or nearly sparsest) coefficients $\bm{c}$ from an optimization problem of the form

where $\|\bm{Wc}\|_{s}$, with $s=\{0,1\}$ and some positive diagonal weight matrix $\bm{W}$, is a measure of the sparsity of $\bm{c}$ and $\|\bm{\Psi}\bm{c}-\bm{u}\|_{2}$ is a measure of the accuracy of the truncated PC expansion in estimating the $u(\omega)$ samples. The $N$-vector $\bm{u}$ contains the independent random samples of $u(\omega)$ and the rows of the $N\times P$ matrix $\bm{\Psi}$ consist of the corresponding samples of the PC basis $\{\psi_{\bm{\alpha}}\}$.

Throughout the rest of this manuscript, we will elaborate on the formulation of the compressive sampling problem (1) and the required conditions under which it leads to an accurate and stable approximation of an arbitrary sparse stochastic function as well as sparse solutions to linear elliptic stochastic PDEs. Although we choose to study this particular class of stochastic PDEs, we stress that the proposed algorithms and theoretical developments are far more general and may be readily applied to recover sparse solution of other stochastic systems.

In Section 2, we describe the setup of the problem of interest. We then, in Section 3.1, briefly overview the spectral stochastic discretization of the random functions using PC basis. The main contribution of this work on sparse approximation of stochastic PDEs using compressive sampling is then introduced in Sections 3.2, 3.4, and 3.5. Sections 3.6 and 3.7 discuss some of the implementation details of the present technique. To demonstrate the accuracy and efficiency of the proposed procedures, in Section 4, we perform two numerical experiments on a $1$-$D$ (in space) linear elliptic stochastic PDE with high-dimensional random diffusion coefficients.

Problem setup

$\mathcal{P}-a.s.\;\omega\in\Omega$. The diffusion coefficient $a(\bm{x},\omega)$ is a stochastic function defined on $(\Omega,\mathcal{F},\mathcal{P})$ and is the source of uncertainty in (2). We assume that $a(\bm{x},\omega)$ is specified by a truncated Karhunen-Loève-“like” expansion

where $(\lambda_{i},\phi_{i})$, $i=1,\cdots,d$, are the eigenpairs of the covariance function $C_{aa}(\bm{x}_{1},\bm{x}_{2})\in L_{2}(\mathcal{D}\times\mathcal{D})$ of $a(\bm{x},\omega)$ and $\bar{a}(\bm{x})$ is the mean of $a(\bm{x},\omega)$. We further assume that $a(\bm{x},\omega)$ satisfies the following conditions:

A-I. For all $\bm{x}\in\mathcal{D}$, there exists constants $a_{\min}$ and $a_{\max}$ such that

A-II. The covariance function $C_{aa}(\bm{x}_{1},\bm{x}_{2})$ is piecewise analytic on $\mathcal{D}\times\mathcal{D}$ (Schwab06a, ; Bieri09a, ), implying that there exist real constants $c_{1}$ and $c_{2}$ such that for $i=1,\cdots,d$,

A-III. The random variables $\{y_{k}(\omega)\}_{k=1}^{d}$ are independent and uniformly distributed on $\Gamma_{k}:=$, $k=1,\cdots,d$, with probability distribution function $\rho_{k}(y_{k})=1/2$ defined over $\Gamma_{k}$. The joint probability distribution function of the random vector $\bm{y}:=(y_{1},\cdots,y_{d})$ is then given by $\rho(\bm{y}):=\prod_{k=1}^{d}\rho_{k}(y_{k})$.

Remark: The algorithm proposed in the paper requires the existence of a sparse solution. The only role of assumption A-II is to guarantee the existence of such a sparse solution. It is not necessary for the application and the validity of the proposed algorithm. In particular, if the coefficient $a$ is only essentially bounded, the proposed algorithm will be accurate as long as a sparse approximation exists.

Given the finite-dimensional uncertainty representation in (3), the solution $u(\bm{x},\omega)$ of (2) also admits a finite-dimensional representation, i.e.,

where $\Gamma:=\prod_{k=1}^{d}\Gamma_{k}$.

In what follows, we first briefly outline the Legendre spectral stochastic discretization of $u(\bm{x},\bm{y})$ and consequently introduce our approach based on compressive sampling to obtain such a discretization.

Numerical approach

In the context of the spectral stochastic methods (Ghanem03, ; Deb01, ; Xiu02, ; Babuska04, ), the solution $u(\bm{x},\bm{y})$ of (2) is represented by an infinite series of the form

We here assume that the univariate Legendre polynomials $\psi_{\alpha_{i}}(y_{i})$ are also normalized such that

The exact generalized Fourier coefficients $c_{\bm{\alpha}}(\bm{x})$ in (8), referred to as the PC coefficients, are computed by the projection of $u(\bm{x},\bm{y})$ onto each basis function $\psi_{\bm{\alpha}}(\bm{y})$,

where the set of multi-indices $\Lambda_{p,d}$ is

Here, $\|\bm{\alpha}\|_{1}=\sum_{i=1}^{d}\bm{\alpha}_{i}$ and $\|\bm{\alpha}\|_{0}=\#\{i:\alpha_{i}>0\}$ are the total order (degree) and dimensionality of the basis function $\psi_{\bm{\alpha}}(\bm{y})$, respectively. The approximation is then refined by increasing $p$ to achieve a given target accuracy. Under assumptions A-I, A-II, and A-III stated in Section 2, the solution $u(\bm{x},\bm{y})$ is analytic with respect to the random variables $\{y_{i}\}_{i=1}^{d}$ (see (Babuska07a, )), and as $p$ increases, the approximation (12) converges exponentially fast in the mean-squares sense (Babuska04, ; Babuska07a, ; Bieri09a, ).

Definition (Sparsity) The solution $u(\bm{x},\bm{y})$ is said to be (nearly) sparse if only a small fraction of coefficients $c_{\bm{\alpha}}(\bm{x})$ in (12) are dominant and contribute to the solution statistics.

As will be described in Section 3.2, a sparse solution $u(\bm{x},\bm{y})$ may be accurately recovered using $N\ll P$ random samples $\{u(\bm{x},\bm{y}_{i})\}_{i=1}^{N}$ using compressive sampling techniques. This has to be compared, for instance, with the least-squares regression-type techniques, (Hosder06, ), that normally require $N\gg P$ samples for an accurate recovery.

2 Sparse recovery using compressive sampling

Compressive sampling is an emerging theory in the field of signal and image processing (Chen98, ; Chen01a, ; Candes06a, ; Donoho06b, ; Candes06b, ; Candes06c, ; Candes07a, ; Cohen09a, ; Bruckstein09, ). It hinges around the idea that a set of incomplete random observations of a sparse signal can be used to accurately, or even exactly, recover the signal (provided that the basis in which the signal is sparse is known). In particular, the number of such observations may be much smaller than the cardinality of the signal. In the context of problem (2), compressive sampling may be interpreted as follows. The solution $u(\bm{x},\bm{y})$ that is sparse, in the sense of Lemma 3.2 defined in Section 3.3, can be accurately recovered using $N\ll P$ random samples $\{u(\bm{x},\bm{y}_{i})\}_{i=1}^{N}$, where $P$ is the cardinality of the Legendre PC basis $\{\bm{\psi}_{\bm{\alpha}}(\bm{y})\}$. We next elaborate on the above statement and address how such a sparse reconstruction is achieved and under what conditions it is successful.

Let $\{u(\bm{y}_{i})\}_{i=1}^{N}$ be i.i.d. random samples of $u(\bm{x},\bm{y})$ for a fixed point $\bm{x}$ in $\mathcal{D}$. For the time being, let us assume that the $p$th-order PC basis $\{\bm{\psi}_{\bm{\alpha}}(\bm{y})\}$ is a complete basis to expand $u(\bm{y})$; we will relax this assumption as we proceed. Given pairs of $\{\bm{y}_{i}\}_{i=1}^{N}$ and $\{u(\bm{y}_{i})\}_{i=1}^{N}$, we write

We are interested in the case that the number $N$ of solution samples is much smaller than the unknown PC coefficients $P$, i.e., $N\ll P$. Without any additional constraints on $\bm{c}$, the underdetermined linear system (16) is ill-posed and, in general, has infinitely many solutions. When $\bm{c}$ is sparse; that is, only a small fraction of the coefficients $c_{\bm{\alpha}}$ are significant, the problem (16) may be regularized to ensure a well-posed solution. Such a regularization may be achieved by seeking a solution $\bm{c}$ with the minimum number of non-zeros. This can be formulated in the optimization problem

In general, the $p$th-order PC basis is not complete for the exact representation of $u(\bm{y})$; therefore, we have to account for the truncation error. This can be accommodated in $(P_{0})$ and $(P_{1})$ by allowing a non-zero residual in the constraint $\bm{\Psi c}=\bm{u}$. Therefore, as in Sections 3.2.1 and 3.2.3 of (Bruckstein09, ), the proposed algorithms in this paper are error-tolerant versions of $(P_{0})$ and $(P_{1})$, with error tolerance $\delta$, i.e.,

respectively. The latter problem is named Basis Pursuit Denoising (BPDN) in (Chen98, ) and may be solved using techniques from quadratic programming. We leave the discussion on the available algorithms for solving problems $(P_{1,\delta})$ and $(P_{0,\delta})$ to Section 3.7. Instead, we henceforth delineate on sufficient conditions under which the BPDN problem $(P_{1,\delta})$ leads to a successful Legendre PC expansion of a general essentially bounded sparse stochastic function $u(\bm{y})$ and, subsequently, the sparse solution $u(\bm{x},\bm{y})$ to the problem (2). Our results are extensions of those in (Donoho06a, ; Bruckstein09, ), adapted to the case where the measurement matrix $\bm{\Psi}$ consists of random evaluations of the Legendre PC basis $\{\psi_{\bm{\alpha}}\}$. With slight differences that will be remarked accordingly, similar results hold for the case of the $(P_{0,\delta})$ problem.

Let $u(\bm{y})$ be an essentially bounded function of i.i.d. random variables $\bm{y}:=(y_{1},\cdots,y_{d})$ uniformly distributed on $\Gamma:=^{d}$. Define

with $S:=|\Lambda_{p,d}^{\epsilon}|$, then with probability

(on the $N$ samples $\{u(\bm{y}_{i})\}_{i=1}^{N}$) and for some constants $c_{1}$ and $c_{2}$, the solution $u_{p}^{1,\delta}$ must obey

Remark: Based on the conditions (22) and (24), the number $N$ of random samples has to grow like $P^{4c_{p,d}}\ln P$ and also proportional to the number of dominant coefficients $S=|\Lambda_{p,d}^{\epsilon}|$. Given any order $p$ of the PC expansion, for sufficiently high-dimensional problems, the constant $c_{p,d}<1/4$ (see Lemma 3.5 and Fig. 1), thus justifying $N\ll P$. In fact, the conditions (22) and (24) are too pessimistic; in practice, the number of random samples required for an accurate recovery is much smaller than the theoretical value in (22). We will elaborate on this statement in Section 3.5.

3 Sparsity of the solution u​(𝒙,𝒚)𝑢𝒙𝒚u(\bm{x},\bm{y})

Notice that the accurate recovery of $u(\bm{y})$ is conditioned upon the existence of a sparse PC expansion $u_{p}^{0}$ (see Theorem 3.1). In fact, this assumption may not hold for an arbitrary stochastic function $u(\bm{y})$, as all the elements of the basis set $\{\psi_{\bm{\alpha}}(\bm{y})\}$ may be important. In this case, our sparse approximation still converges to the actual solution but, perhaps, not using as few as $N\ll P$ random solution samples.

We will now summarize the results of (Todor07a, ; Bieri09a, ) on the sparsity of the Legendre PC expansion of the solution $u(\bm{x},\bm{y})$ to the problem (2). Alternative to the $p$th-order truncated PC expansion of (12), one may ideally seek a proper index set $\Lambda_{p,d}^{\epsilon}\subseteq\Lambda_{p,d}$, with sufficiently large $p$, such that for a given accuracy $\epsilon$

where $\Lambda_{p,d}$ is defined in (13). Such a reduction in the number of basis functions in (28) is possible as, given the accuracy $\epsilon$, the effective dimensionality $\nu$ of $u(\bm{x},\bm{y})$ in $\Gamma$ is potentially smaller than the apparent dimensionality $d$. More precisely, under assumptions A-I, A-II, and A-III stated in Section 2, the analyses of (Todor07a, ; Bieri09a, ) imply that the discretization of $u(\bm{x},\bm{y})$ using a sparse index set $\Lambda_{p,\nu}$,

preserves the exponential decay of the approximation error in the $H_{0}^{1}\left(\mathcal{D},L^{\infty}(\Gamma)\right)$ sense. For the sake of completeness, we cite this from (Bieri09a, ) in the following lemma.

Given assumptions A-I, A-II, and A-III in Section 2, there exist constants $c_{1},c_{2},c_{3},c_{4}>0$, depending only on $a(\bm{x},\omega)$ and $f(\bm{x})$ but independent of $d,p,\nu$, such that

In particular, for $d\geq c_{d}|\ln\epsilon|^{1/\kappa}$, choosing

where $u_{p_{\epsilon},\nu_{\epsilon}}$ is now defined on a sparse index set

for some arbitrary large $\rho>0$ and constants $c_{d},c_{p}$, and $c_{\nu}$ independent of $d,p_{\epsilon}$, and $\nu_{\epsilon}$ (Bieri09a, ).

In practice, the sparse set $\Lambda_{p,d}^{\epsilon}$ in (27) (or equivalently $\Lambda_{p_{\epsilon},\nu_{\epsilon}}$ in (33)) is not known a priori. In (Bieri09a, ), an approach based on an algebraic purely-stochastic problem is proposed to adaptively identify $\Lambda_{p,d}^{\epsilon}$. Having done this, the coefficients of the the spectral modes are computed via the (intrusive) stochastic Galerkin scheme (Ghanem03, ; Xiu02, ). Alternatively, in this work, we apply our sparse approximation using $(P_{1,\delta})$ and $(P_{0,\delta})$ to compute $u(\bm{x},\bm{y})$. The implementation of $(P_{1,\delta})$ and $(P_{0,\delta})$ is non-intrusive; only random samples of the solution are needed. Moreover, we do not adapt the sampling process to identify the important dimensions at the stochastic level; therefore, our constructions are non-adapted.

Throughout the rest of the present paper, we focus our attention on the case of the stochastic PDE (2) whose solution is provably sparse. The statement of Theorem 3.1 can be specialized to the approximation of the sparse solution to the stochastic PDE (2) using $(P_{1,\delta})$ (or $(P_{0,\delta})$) as follows.

Combining Lemma 3.2 with Theorem 3.1 leads to the following theorem.

There exists constants $c_{1},c_{2},c_{3},c_{4},c_{5}$ independent from $p,d,\kappa,N$ such that if $\lceil c_{2}d^{\kappa}\rceil\leq p$ and $\lceil c_{3}d^{\kappa/(\kappa+1)}\rceil\leq d$,

the solution $u_{p}^{1,\delta}$ must obey

The ability of problems $(P_{1,\delta})$ and $(P_{0,\delta})$ in accurately approximating the sparse PC coefficients $\bm{c}$ in (12), hence the solution $u(\bm{x},\bm{y})$, depends on two main factors: $i)$ the sparsity of the PC coefficients $\bm{c}$ and $ii)$ the mutual coherence of the measurement matrix $\bm{\Psi}$. In fact, the number $N$ of random solution samples required for a successful sparse approximation is dictated by these two factors. While sparsity is a characteristic of the solution of interest $u(\bm{x},\bm{y})$, the mutual coherence of the measurement matrix $\bm{\Psi}$ is universal as it only depends on the choice of PC basis $\{\psi_{\bm{\alpha}}(\bm{y})\}$ and the sampling process from which $\bm{\Psi}$ is assembled.

In Section 3.3, based on the analysis of (Bieri09a, ), we rationalized the sparsity of $u(\bm{x},\bm{y})$ with respect to the Legendre PC basis. We now give the definition of the mutual coherence of $\bm{\Psi}$ and discuss its role in our sparse approximation using $(P_{1,\delta})$ and $(P_{0,\delta})$.

In plain words, the mutual coherence is a measure of how close to orthogonal a matrix is. Clearly, for any general matrix $\bm{\Psi}$,

where the lower bound is achieved, for instance, by unitary matrices. However, for the case of $N<P$, the mutual coherence $\mu(\bm{\Psi})$ is strictly positive. It is well understood that measurement matrices with smaller mutual coherence have a better ability to recover a sparse solution using compressive sampling techniques, e.g., see Lemma 3.6. Therefore, we shall proceed to examine the mutual coherence of the random measurement matrix $\bm{\Psi}$ in (16). We first observe that, by the orthogonality of the Legendre PC basis, the mutual coherence $\mu(\bm{\Psi})$ converges to zero almost surely for asymptotically large random sample sizes $N$. However, it is essential for our purpose to $i)$ investigate if a desirably small $\mu(\bm{\Psi})$ can be achieved by a sample size $N\ll P$ and $ii)$ quantify how large $\mu(\bm{\Psi})$ can get for a finite $N$. These are addressed in the following theorem.

Figure 1 illustrates the decay of $c_{p,d}$, for several values of $p$, as a function of $d$. Based on Theorem 3.4, for cases where the number $d$ of random variables $\bm{y}$ is large enough such that $c_{p,d}<1/4$, it is sufficient to have $N\sim\mathcal{O}(16P^{4c_{p,d}}\ln P)\ll P$ to keep $\mu(\bm{\Psi})$ bounded from above with a large probability. Notice that such a requirement on $c_{p,d}$ is particularly suited to high-dimensional problems.

Remark: We observe that, given the choice of $r$ in (40), the upper bound on $\mu(\bm{\Psi})$ in (41) decays like $1/\sqrt{N}$ for asymptotically large $N$, which is consistent with the Central Limit Theorem.

In order to prove Theorem 3.4, we first need to compute the maximum of the Legendre PC basis functions $\psi_{\bm{\alpha}}(\bm{y})$. This is given in the following lemma.

Let $\{\psi_{\bm{\alpha}}(\bm{y})\}$ be the Legendre polynomial chaos basis of total order $p$ in $d$ i.i.d uniform random variables $\bm{y}$ (as defined in (13)) and with cardinality $P$. Then,

Proof of Theorem 3.4. The mutual coherence $\mu(\bm{\Psi})$ is

Given the independence of samples $\{\bm{y}_{i}\}_{i=1}^{N}$ and using the McDiarmid’s inequality, we obtain

Using Lemma 3.5, we may probabilistically bound the numerator in (44) as

for some $\zeta>1$, we arrive at the statement of the Theorem 3.4. $\square$

To summarize, we observe that with large probability, the mutual coherence $\mu(\bm{\Psi})$ of the measurement matrix $\bm{\Psi}$ in (16) can be arbitrarily bounded from above by increasing the number $N$ of independent random solution samples. Moreover, given the discussions of Section 3.3, we know that the solution to problem (2) is sparse in the Legendre PC basis. These are the two key factors affecting the stability and accuracy of our sparse approximation.

Let $u_{p}^{0}(\bm{y}):=\sum_{\bm{\alpha}\in\Lambda_{p,d}}c_{\bm{\alpha}}^{0}\psi_{\bm{\alpha}}(\bm{y})$ with $c_{\bm{\alpha}}^{0}=0$ for $\bm{\alpha}\notin\Lambda_{p_{\epsilon},\nu_{\epsilon}}$ be the sparse Legendre PC approximation of $u(\bm{x},\bm{y})$ at a spatial point $\bm{x}$ where the sparse index set $\Lambda_{p_{\epsilon},\nu_{\epsilon}}$ is defined in (33). Assume that the vector of PC coefficients $\bm{c}^{0}$ satisfies the sparsity condition

Proof. Using Theorem 3.1 of (Donoho06a, ), we obtain that if $\bm{c}^{0}$ satisfies the sparsity condition $\|\bm{c}^{0}\|_{0}<(1+1/\mu(\bm{\Psi}))/4$, then

Remark: The error bound in (50) is not tight; in fact, the actual error is significantly smaller than the upper bound given in (54). More importantly, according to (Donoho06a, ), the sparsity condition (49) is unnecessarily too restrictive. In practice, both far milder sparsity conditions are needed and much better actual errors are achieved.

Remark: We will later use the sparsity condition (49) to derive the sufficient condition (22) (together with (24)) on the number $N$ of random samples needed for a successful recovery. As the condition (49) is too restrictive, the theoretical lower bound on $N$ given in (22) and (24) is too pessimistic.

Remark: According to Lemma 3.6, we do not need to know a priori the sparse index set $\Lambda_{p_{\epsilon},\nu_{\epsilon}}$; only the sparsity condition (49) is required.

Remark: We stated Lemma 3.6 for the case where the sparsity of the PC expansion is due to the fact that the effective dimensionality $\nu_{\epsilon}$ is potentially smaller than $d$. However, as far as the stability condition (49) is satisfied, similar stability results are valid for situations where dominant basis are defined over all the dimensions.

Notice that the normalized truncation error

Although our sparse approximations are point-wise in space, we are ultimately interested in deriving suitable global stability and error estimates for our sparse reconstructions. Such extensions are readily available from Lemma 3.6 and are stated in the following corollary.

We have now all the necessary tools to proceed with the proof of our main result stated in Theorem 3.3 which is primarily a direct consequence of Lemma 3.2, Theorem 3.4, and Corollary 3.7.

Proof of Theorem 3.3. We first note that, given the conditions of Lemma 3.2, the solution to problem (2) admits a sparse Legendre PC expansion $u_{p}$ with sparsity $S=|\Lambda_{p_{\epsilon},\nu_{\epsilon}}|\lesssim\epsilon^{-1/\rho}$ when $\epsilon\geq\exp\left(-\left(\frac{d}{c_{1}}\right)^{\kappa}\right)$ for some constants $c_{1}$ and (arbitrary) $\rho>0$. Notice that the sparse approximation $u_{p}$ has an accuracy better than $\epsilon$ in the $H_{0}^{1}(\mathcal{D},L^{\infty}(\Gamma))$ sense. Based on Theorem 3.4, it is sufficient to have random solution samples of size $N\geq 64P^{4c_{p,d}}(\ln P)S$, to meet the sparsity requirement $S=|\Lambda_{p_{\epsilon},\nu_{\epsilon}}|<(1+1/\mu(\bm{\Psi}))/4$, in Corollary 3.7, with probability at least $1-4P^{2-2S_{\max}}$ where $S_{\max}:=\frac{N}{64P^{4c_{p,d}}(\ln P)}$. On the other hand, given $N$ random samples of solution, we require $\epsilon\geq\frac{1}{S_{\max}^{\rho}}$ to satisfy the sparsity condition. Given Corollary 3.7 and using the triangular and Poincaré inequalities, with probability at least $1-P^{-8S_{\max}}$, we have

for all $\delta>0$. Moreover, by choosing $r=\frac{1}{4}\left\|u-u_{p}^{0}\right\|_{L^{2}\left(\mathcal{D},L^{\infty}(\Gamma)\right)}^{2}$ in (55),

with probability at least $1-P^{-8S_{\max}P^{4c_{p,d}}}$. Finally, by taking

we arrive at the statement of the Theorem 3.3. $\square$

6 Choosing the truncation error tolerance δ𝛿\delta

An important component of the sparse approximation using $(P_{1,\delta})$ and $(P_{0,\delta})$ is the selection of the truncation error tolerance $\delta$. Although the stability bounds given in Lemma 3.6 and Corollary 3.7 are valid for any $\delta\geq 0$, the actual error and the sparsity level of the solution to $(P_{1,\delta})$ and $(P_{0,\delta})$ depend on the choice of $\delta$. Ideally, we desire to choose $\delta\approx\|\bm{\Psi}\bm{c}^{0}-\bm{u}\|_{2}$; while larger values of $\delta$ deteriorate the accuracy of the approximation, as in Lemma 3.6, smaller choices of $\delta$ may result in over-fitting the solution samples and, thus, less sparse solutions. In practice, as the exact values of the PC coefficients $\bm{c}^{0}$ are not known, the exact values of the truncation error $\|\bm{\Psi}\bm{c}^{0}-\bm{u}\|$ and, consequently, $\delta$ are not known a priori. Therefore, $\delta$ has to be estimated, for instance, using statistical techniques such as the cross-validation (Boufounos07, ; Ward09, ).

In this work, we propose a heuristic cross-validation algorithm to estimate $\delta$. We first divide the $N$ available solution samples to $N_{r}$ reconstruction and $N_{v}$ validation samples such that $N=N_{r}+N_{v}$. The idea is to repeat the solution of $(P_{1,\delta})$ (or $(P_{0,\delta})$) on the reconstruction samples and with multiple values of truncation error tolerance $\delta_{r}$. We then set $\delta=\sqrt{\frac{N}{N_{r}}}\hat{\delta}_{r}$ in which $\hat{\delta}_{r}$ is such that the corresponding truncation error on the $N_{v}$ validation samples is minimum. This is simply motivated by the fact that the truncation error on the validation samples is large for values of $\delta_{r}$ considerably larger and smaller than $\|\bm{\Psi}\bm{c}^{0}-\bm{u}\|_{2}$ evaluated using the reconstruction samples. While the former is expected from the upper bound on the approximation error in Lemma 3.6, the latter is due to the over-fitting the reconstruction samples. The following exhibit outlines the estimation of $\delta$ using the above cross-validation approach:

Algorithm for cross-validation estimation of $\delta$: • Divide the $N$ solution samples to $N_{r}$ reconstruction and $N_{v}$ validation samples. • Choose multiple values for $\delta_{r}$ such that the exact truncation error $\|\bm{\Psi}\bm{c}^{0}-\bm{u}\|_{2}$ of the reconstruction samples is within the range of $\delta_{r}$ values. • For each value of $\delta_{r}$, solve $(P_{1,\delta})$ (or $(P_{0,\delta})$) on the $N_{r}$ reconstruction samples. • For each value of $\delta_{r}$, compute the truncation error $\delta_{v}:=\|\bm{\Psi}\bm{c}^{1,\delta_{r}}-\bm{u}\|_{2}$ (or $\delta_{v}:=\|\bm{\Psi}\bm{c}^{0,\delta_{r}}-\bm{u}\|_{2}$) of the $N_{v}$ validation samples. • Find the minimum value of $\delta_{v}$ and its corresponding $\hat{\delta}_{r}:=\delta_{r}$. • Set $\delta=\sqrt{\frac{N}{N_{r}}}\hat{\delta}_{r}$.

In the numerical experiments of Section 4, we repeat the above cross-validation algorithm for multiple replications of the reconstruction and validation samples. The estimate of $\delta=\sqrt{\frac{N}{N_{r}}}\hat{\delta}_{r}$ is then based on the value of $\hat{\delta}_{r}$ for which the average of the corresponding truncation errors $\delta_{v}$, over all replications of the validation samples, is minimum. This resulted in more accurate solutions in our numerical experiments.

7 Algorithms

There are several numerical algorithms for solving problems $(P_{0,\delta})$ and $(P_{1,\delta})$ each with different optimization kernel, computational complexity, and degree of accuracy. An in-depth discussion on the performance of these algorithms is outside the scope of the present work; however, below we name some of the available options for each problem and briefly describe the algorithms that have been utilized in our numerical experiments. For comprehensive discussions on this subject, the interested reader is referred to (Bruckstein09, ; Berg08, ; Figueiredo07, ; Becker09, ; Yang09, ; Tropp10a, ).

Problem $(P_{0,\delta})$: A brute force search through all possible support sets in order to identify the correct sparsity for the solution $\bm{c}^{0,\delta}$ of $(P_{0,\delta})$ is NP-hard and not practical. Greedy pursuit algorithms form a major class of schemes to tackle the solution of $(P_{0,\delta})$ with a tractable computational cost. Instead of performing an exhaustive search for the support of the sparse solution, these solvers successively find one or more components of the solution that result in the largest improvement in the approximation. Some of the standard greedy pursuit algorithms are Orthogonal Marching Pursuit (OMP) (Pati93, ; Davis97, ), Regularized OMP (ROMP)(Needell07, ), Stagewise OMP (StOMP) (Donoho06c, ), Compressive Sampling MP (CoSaMP) (Needell08a, ), Subspace Pursuit (Dai09, ), and Iterative Hard Thresholding (IHT) (Blumensath09, ). Under well-defined conditions, all of the above schemes provide stable and accurate solutions to $(P_{0,\delta})$ in a reasonable time.

Orthogonal Matching Pursuit (OMP) Algorithm: • Set $k=0$. – Set the initial solution $\bm{c}^{0,\delta,(0)}=\bm{0}$ and residual $\bm{r}^{(0)}=\bm{u}-\bm{\Psi}\bm{c}^{0,\delta,(0)}=\bm{u}$. – Set the solution support index set $\mathcal{I}^{(0)}=\emptyset$. • While $\|\bm{u}-\bm{\Psi c}^{0,\delta,(k)}\|_{2}>\delta$ perform: – For all $j\notin\mathcal{I}^{(k)}$ evaluate $\epsilon(j)=\|\bm{\psi}_{j}\alpha_{j}-\bm{r}^{(k)}\|_{2}$ with $\alpha_{j}=\bm{\psi}_{j}^{T}\bm{r}^{(k)}/\|\bm{\psi}_{j}\|_{2}^{2}$. – Set $k=k+1.$ – Update the support index set $\mathcal{I}^{(k)}=\mathcal{I}^{(k-1)}\bigcup\left\{\arg\min_{j}\epsilon(j)\right\}$. – Solve for $\bm{c}^{0,\delta,(k)}=\arg\min_{{\bm{c}}^{0,\delta}}\|\bm{u}-\bm{\Psi c}^{0,\delta}\|_{2}$ subject to $Support\{\bm{c}^{0,\delta}\}=\mathcal{I}^{(k)}$. – Update the residual $\bm{r}^{(k)}=\bm{u}-\bm{\Psi}\bm{c}^{0,\delta,(k)}$ • Output the solution $\bm{c}^{0,\delta}=\bm{c}^{0,\delta,(k)}$.

Although we chose OMP in our analysis, we note that further studies are needed to identify the most appropriate greedy algorithm for the purpose of this study.

In the next section, we explore some aspects of the proposed scheme through its application to a $1$-$D$ (in space) elliptic stochastic PDE with high-dimensional random diffusion coefficients.

Numerical examples

We consider the solution of a one-dimensional, i.e., $D=1$, version of problem (2),

where the stochastic diffusion coefficient $a(x,\omega)$ is given by the expansion

Here, $\{\lambda_{i}\}_{i=1}^{d}$ and $\{\phi_{i}(x)\}_{i=1}^{d}$ are, respectively, $d$ largest eigenvalues and the corresponding eigenfunctions of the Gaussian covariance kernel

in which $l_{c}$ is the correlation length of $a(x,\omega)$ that prescribes the decay of the spectrum of $C_{aa}$ in (60). Random variables $\{y_{i}(\omega)\}_{i=1}^{d}$ are assumed to be independent and uniformly distributed on $$. The coefficient$\sigma_{a}$controls the variability of$a(x,\omega)$.

We verify the accuracy and efficiency of the present sparse approximation schemes for both moderate and high-dimensional diffusion coefficient $a(x,\omega)$. These two cases are obtained, respectively, by assuming $(l_{c},d)=(1/5,14)$ and $(l_{c},d)=(1/14,40)$ in (60) and (59). We further assume that $\bar{a}=0.1$, $\sigma_{a}=0.03$ when $d=14$, and $\sigma_{a}=0.021$ when $d=40$. These choices ensure that all realizations of $a(x,\omega)$ are strictly positive on $\mathcal{D}=(0,1)$. Table LABEL:tab:parameters summarizes the assumed parameters for the two test cases.

For both cases, the spatial discretization is done by the Finite Element Method using quadratic elements. A mesh convergence analysis is performed to ensure that spatial discretization errors are inconsequential.

As elucidated in Section 3.5, the accuracy of our sparse reconstruction depends on the mutual coherence $\mu(\bm{\Psi})$, the sample size $N$, and the truncation error $\|\bm{\Psi c}-\bm{u}\|_{2}$ (hence $\delta$). In order to reduce the approximation error, we need to reduce $\|\bm{\Psi c}-\bm{u}\|_{2}$, which may be done by increasing $p$ and, therefore, $P$. However, with a fixed number $N$ of samples, an increase in $P$ may result in a larger mutual coherence and, thus, the degradation of the reconstruction accuracy. Therefore, in practice, we start by approximating the lower order PC expansions when $N$ is small and increase $p$ when larger number of samples become available. Notice that such an adaptivity with respect to the order $p$ is a natural way of refining the accuracy of PC expansions, for instance, when the intrusive stochastic Galerkin scheme is adopted (Ghanem03, ). In particular, in this example, for sample sizes $N=\{29,120\}$, we attempt to estimate the coefficients of the $3$rd-order Legendre PC expansion, i.e. $p=3$ and $P=680$. For larger sample sizes $N$, we also include the first $320$ basis function from the $4$th-order chaos, thus resulting in $P=1000$. Since all of the $4$th-order basis functions are not employed, we need to describe the ordering of our basis construction. We sort the elements of $\{\psi_{\bm{\alpha}}(\bm{y})\}$ such that, for any given order $p$, the random variables $y_{i}$ with smaller indices $i$ contribute first in the basis.

For each analysis, we estimate the truncation error tolerance $\delta$ based on the cross-validation algorithm described in Section 3.6. For each $N$, we use $N_{r}\approx 3N/4$ of the samples (reconstruction set) to compute the PC coefficients $\bm{c}^{1,\delta_{r}}$ and the rest of the samples (validation set) are used to evaluate the truncation error $\delta_{v}$. The cross-validation is performed for four replications of reconstruction and validation sample sets. We then find the value $\hat{\delta}_{r}$ that minimizes the average of $\delta_{v}$ over the four replications of the cross-validation samples. Given an estimate of the truncation error tolerance $\delta\approx\sqrt{4/3}\hat{\delta}_{r}$, we then use all $N$ samples to compute the coefficients $\bm{c}^{1,\delta}$.

The convergence of the mean, standard deviation, and root mean-squares of the approximation error for $u(0.5,\bm{y})$ is illustrated in Figs. 3 (a), (b), and (c), respectively. For the case of stochastic collocation, we apply sparse grid quadrature (cubature) integration rule to directly compute the mean and the standard deviation. The root mean-squares error of the Monte Carlo and the stochastic collocation solution are evaluated by estimating the corresponding PC coefficients using sampling and sparse grid quadrature integration, respectively, and then comparing them with the exact coefficients.

2 Case II: d=40

The objective of this example is to highlight that a sparse reconstruction may lead to significant computational savings for problems with high-dimensional random inputs. Similar to the analysis of Case I described in Section 4.1, we compute the solution statistics using multiple numbers of independent samples. More specifically, we evaluate the solution at $x=0.5$ for independent samples of size $N=\{81,200,400,600,800,1000\}$. The number of grid points in the level $l=1$ and $l=2$ of the Clenshaw-Curtis rule in dimension $d=40$ is $N=81$ and $N=3281$, respectively. To obtain a reference solution, the $3$rd order PC coefficients $\bm{c}$ of the solution at $x=0.5$ are computed using level $l=5$ stochastic collocation with the Clenshaw-Curtis rule.

For $N=\{81,200\}$ we only estimate the coefficients associated with the $2$nd-order PC expansion, i.e. $p=2$ and $P=861$. For larger sample sizes, we also include the first $639$ basis functions from the $3$th-order chaos, thus leading to $P=1500$. For each combination of $N$ and $p$, we estimate the truncation error $\delta$ using an identical cross-validation procedure described in Section 4.2. Figure 4 illustrates the estimation of PC coefficients of $u(0.5,\bm{y})$ with BPDN and OMP algorithms with $N=200$ and $N=1000$. We again note that the recovered solution from the BPDN algorithm is less sparse as compared to that of the OMP approach, although the over-estimated coefficients (mostly from the second order term) are indeed small. As the samples size $N$ is increased, we are naturally able to recover more dominant coefficients on the expansion. Figure 5 depicts the convergence of the statistics of the solution as functions of the sample size $N$ as well as one instance of the estimation of the truncation error tolerance $\delta$. The implementation details are similar to those described in Section 4.1 for the case of $d=14$.

Conclusion

Acknowledgments

The first author acknowledges the support of the United States Department of Energy under Stanford’s Predictive Science Academic Alliance Program (PSAAP) for the preliminary stages of his work. The second author acknowledges the support of the National Science Foundation via NSF grant CMMI-092600 and of the United States Department of Energy under Caltech’s Predictive Science Academic Alliance Program (PSAAP).

References