Reconstructing high-dimensional Hilbert-valued functions via compressed sensing

TL;DR

This paper introduces a new sparse polynomial method for approximating high-dimensional Hilbert-valued functions, with theoretical guarantees and an algorithm, applied to parameterized PDEs with stochastic inputs.

Contribution

It develops a mixed-norm $$ regularization technique for joint sparse recovery of Hilbert-valued functions, extending compressed sensing theory to infinite-dimensional settings.

Findings

Proves recovery guarantees for Hilbert-valued functions.

Provides an algorithm with strong convergence in infinite dimensions.

Demonstrates minimal sample complexity through numerical experiments.

Abstract

We present and analyze a novel sparse polynomial technique for approximating high-dimensional Hilbert-valued functions, with application to parameterized partial differential equations (PDEs) with deterministic and stochastic inputs. Our theoretical framework treats the function approximation problem as a joint sparse recovery problem, where the set of jointly sparse vectors is possibly infinite. To achieve the simultaneous reconstruction of Hilbert-valued functions in both parametric domain and Hilbert space, we propose a novel mixed-norm based $\ell_1$ regularization method that exploits both energy and sparsity. Our approach requires extensions of concepts such as the restricted isometry and null space properties, allowing us to prove recovery guarantees for sparse Hilbert-valued function reconstruction. We complement the enclosed theory with an algorithm for Hilbert-valued recovery, based on standard forward-backward algorithm, meanwhile establishing its strong convergence in the considered infinite-dimensional setting. Finally, we demonstrate the minimal sample complexity requirements of our approach, relative to other popular methods, with numerical experiments approximating the solutions of high-dimensional parameterized elliptic PDEs.

Equations

$$u=\sum _{\mathbf{\nu }\in \mathcal{F}}{\mathbf{c}}_{\mathbf{\nu }}{\Psi }_{\mathbf{\nu }}\approx \sum _{\mathbf{\nu }\in \mathcal{J}}{\mathbf{c}}_{\mathbf{\nu }}{\Psi }_{\mathbf{\nu }}=:u_{\mathcal{J}}.$$

$$\Vert u-u_{\mathcal{J}}{\Vert }_{L^2}^2=\Vert u_{{\mathcal{J}}^c}{\Vert }_{L^2}^2=\sum _{\mathbf{\nu }\in \mathcal{F}\setminus \mathcal{J}}\Vert {\mathbf{c}}_{\mathbf{\nu }}{\Vert }_{\mathcal{V}}^2$$

$$\mathcal{L}(u(\cdot ,\mathbf{y}),\mathbf{y})=0,\text{ in }D,$$

$$\mathbf{u}=\mathbf{A}\mathbf{z},\mathbf{z}\in {\mathcal{V}}^N,$$

$${\text{min}}_{\mathbf{z}\in {\mathcal{V}}^N}\Vert \mathbf{z}{\Vert }_{\mathcal{V},1}\text{subject to}\Vert \mathbf{A}\mathbf{z}-\mathbf{u}{\Vert }_{\mathcal{V},2}\le \frac{\eta }{\sqrt{m}},$$

$${\text{min}}_{\mathbf{z}\in {\mathcal{V}}^N}\Vert \mathbf{z}{\Vert }_{\mathcal{V},1}+\frac{\mu }{2}\Vert \mathbf{A}\mathbf{z}-\mathbf{u}{\Vert }_{\mathcal{V},2}^2,$$

$${\mathbf{c}}_{\mathbf{\nu }}=\sum _{r\in \mathbb{N}}{c}_{\mathbf{\nu },r}{\varphi }_r,\text{with }{c}_{\mathbf{\nu },r}\in \mathbb{R}.$$

$$\Vert \mathbf{c}{\Vert }_{\mathcal{V},q}={\left(\sum _{\mathbf{\nu }\in \mathcal{J}}\Vert {\mathbf{c}}_{\mathbf{\nu }}{\Vert }_{\mathcal{V}}^q\right)}^{1/q}$$

$$\equiv {\left(\sum _{\mathbf{\nu }\in \mathcal{J}}\Vert {\hat{\mathbf{c}}}_{\mathbf{\nu }}{\Vert }_2^q\right)}^{1/q}=\Vert \hat{\mathbf{c}}{\Vert }_{2,q},$$

$$\Vert {\mathbf{z}}_S{\Vert }_{\mathcal{V},2}\le \frac{\rho }{\sqrt{s}}\Vert {\mathbf{z}}_{S^c}{\Vert }_{\mathcal{V},1}+\tau \Vert \mathbf{A}\mathbf{z}{\Vert }_{\mathcal{V},2}$$

$$\forall \mathbf{z}\in {\mathcal{V}}^N,\forall S\subset [N]\text{with}\#(S)\le s.$$

$$(1-{\delta }_{2s})\Vert \mathbf{z}{\Vert }_2^2\le \Vert \mathbf{A}\mathbf{z}{\Vert }_2^2\le (1+{\delta }_{2s})\Vert \mathbf{z}{\Vert }_2^2,$$

$$\forall \mathbf{z}\in {\mathbb{R}}^N,\mathbf{z}2s\text{-sparse},$$

$$m\ge C{\Theta }^{2}s\max \{{\log }^2({\Theta }^{2}s)\log (N),$$

$$\log ({\Theta }^{2}s)\log (\log ({\Theta }^{2}s)N^{\log (s)})\},$$

$${\text{minimize}}_{\mathbf{z}\in {\mathcal{V}}^N}\Vert \mathbf{z}{\Vert }_{\mathcal{V},1}\text{s.t.}\Vert \mathbf{A}\mathbf{z}-\mathbf{u}{\Vert }_{\mathcal{V},2}\le \frac{\eta }{\sqrt{m}}$$

$${\Vert u-\sum _{\mathbf{\nu }\in \mathcal{J}}{\mathbf{c}}_{\mathbf{\nu }}^{\#}{\Psi }_{\mathbf{\nu }}\Vert }_{L_{\varrho }^2(\mathcal{U};\mathcal{V})}\le C_{\epsilon }\sqrt{s}\exp \left[-{\left(\frac{\kappa s}{(1+\epsilon )}\right)}^{1/d}\right],$$

$$0\in T(\mathbf{x})$$

$${\mathbf{x}}^{k+1}=(I+\tau T_1{)}^{-1}(I-\tau T_2){\mathbf{x}}^k,$$

$$J_{\tau }:=(I+\tau T_1{)}^{-1},G_{\tau }:=(I-\tau T_2),\text{ and }{S}_{\tau }:=J_{\tau }\circ G_{\tau },$$

$$G_{\tau }(\mathbf{x})$$

$$J_{\tau ,j}(\mathbf{x})$$

$$\Vert G_{\tau }(\mathbf{v})-G_{\tau }(\mathbf{w}){\Vert }_{2,2}\le \Vert \mathbf{v}-\mathbf{w}{\Vert }_{2,2}\forall \mathbf{v},\mathbf{w}\in \mathcal{H},$$

$$\begin{array}{rl}& \Vert J_{\tau }({\mathbf{v}}_j)-J_{\tau }({\mathbf{w}}_j){\Vert }_2^2\\ & \le \Vert {\mathbf{v}}_j-{\mathbf{w}}_j{\Vert }_2^2-\Vert (I-J_{\tau }){\mathbf{v}}_j-(I-J_{\tau }){\mathbf{w}}_j{\Vert }_2^2.\end{array}$$

$$\displaystyle\left\{\begin{array}[]{rll}-\nabla\cdot\left(a(x,{\bm{y}})\nabla u(x,{\bm{y}})\right)&=f(x)&\textnormal{in }D,\\ u(x,{\bm{y}})&=0&\textnormal{on }\partial D.\end{array}\right.$$

$$a(x,\mathbf{y})$$

$${\zeta }_i$$

$${\vartheta }_i(x)$$

Full text

I Introduction

Underlying many successful applications of compressed sensing to problems in applied mathematics and the physical sciences is the fact that, for many practical problems, the target to be reconstructed possesses sufficient sparsity in order to enable unique solutions from systems that would otherwise be ill-posed. In the basic compressed sensing problem, the signal is an unknown vector ${\bm{c}}\in\mathbb{R}^{N}$, and the sensing process yields a measurement vector ${\bm{u}}\in\mathbb{R}^{m}$ that is formed by the product of ${\bm{c}}$ with a sensing matrix, i.e., ${\bm{u}}={\bm{A}}{\bm{c}}$, where ${\bm{A}}\in\mathbb{R}^{m\times N}$. The key observation is that when the signal ${\bm{c}}$ is sufficiently sparse, it can be uniquely determined from an underdetermined set of measurements ($m\ll N$), provided ${\bm{A}}$ satisfies certain additional properties. To overcome the NP-hardness of directly finding the sparsest ${\bm{c}}$ consistent with a given set of measurements, various greedy and convex relaxation strategies have been proposed and demonstrated, both empirically and theoretically, to have good reconstruction performance in a range of settings.

In the parameterized PDE literature, one typically seeks to approximate the parameter-to-solution map ${\bm{y}}\mapsto u({\bm{y}})$, taking values in a Hilbert space ${\cal V}$, by a truncation of its orthonormal expansion, i.e.,

[TABLE]

Here, ${\bm{y}}\in\mathcal{U}$ a high-dimensional tensor product domain, ${\cal J}$ is a finite multi-index set of cardinality $N=\#({\cal J})$ with ${\cal J}\subset{\cal F}:={\mathbb{N}}_{0}^{d}$, $(\Psi_{\bm{\nu}})_{{\bm{\nu}}\in\mathcal{J}}$ is an $L^{2}(\mathcal{U},d{\varrho})$-orthonormal basis, and ${\bm{c}}=({\bm{c}}_{\bm{\nu}})_{{\bm{\nu}}\in{\cal J}}\in\mathcal{V}^{N}$ are the Hilbert-valued coefficients to be computed. Generally, one selects ${\cal J}$ to be a set large enough to ensure that

[TABLE]

is minimal. However, this can lead to a less efficient approximation if ${\cal J}$ is not chosen carefully.

Over the course of the last decade, a series of works in the parameterized PDEs community (see [4] and the reference therein) have demonstrated that, under reasonable assumptions on the input data to the PDEs, the solutions are compressible, and hence well-represented by sparse expansions in given orthonormal systems. In other words, the solution vector ${\bm{c}}$ from (1) is sparse, and accurate reconstructions of the most important components of ${\bm{c}}$ are enough for satisfactory approximations of the PDE solutions. The ability of compressed sensing (CS) to exploit sparsity and allow far fewer samples than traditional approaches (e.g., Monte-Carlo, projection, interpolation) makes it a promising tool for such reconstruction problems. Therefore it is no surprise that CS-based polynomial approximation has attracted growing interest in the area of high-dimensional parameterized PDEs in recent years, [8, 13, 21, 17, 15, 12, 16].

However, there has been a critical mismatch between standard CS techniques and the problem of reconstructing Hilbert-valued solutions to parameterized PDEs: these methods do not enable direct recovery of the vector $({\bm{c}}_{\bm{\nu}})_{{\bm{\nu}}\in{\cal J}}\in{\cal V}^{N}$, i.e., a vector with Hilbert-valued components. Instead, CS-based polynomial approximation methods only allow the recovery of real or complex sparse vectors. Hence, in the context of parameterized PDEs, standard CS-based polynomial approximation methods do not yield approximations of the entire solution map ${\bm{y}}\mapsto u(\cdot,{\bm{y}})\in{\cal V}$, but only functionals of the solution, i.e., maps of the form ${\bm{y}}\mapsto G(u({\bm{y}}))\in\mathbb{R}$, where $G:{\cal V}\to\mathbb{R}$ is a functional of $u$. In particular, many of the existing works (cited above) perform point-wise reconstruction of the solution, i.e., reconstruction of the map ${\bm{y}}\mapsto G_{x^{\star}}(u({\bm{y}}))$ with $G_{x^{\star}}(u({\bm{y}}))=u(x^{\star},{\bm{y}})$ for a fixed $x^{\star}\in D$, with $D$ the physical domain of the PDE.

In this paper, we present an overview of latest progress in addressing this mismatch by developing a novel sparse recovery technique which enables direct reconstruction of Hilbert-valued vectors, [5, 6]. The theoretical framework, based on the problem of joint-sparse recovery, treats the recovery problem as a matrix recovery problem in which each row of the solution matrix may have infinitely many terms corresponding to the coefficients of a Hilbert-valued function in a given basis. Our regularization is performed with respect to a mixed norm $\|\cdot\|_{{\cal V},1}$, which is defined to be the $\ell_{1}$ norm of the vector $(\|{\bm{c}}_{\bm{\nu}}\|_{{\cal V}})_{{\bm{\nu}}\in{\cal J}}$. Often, the decay of the polynomial coefficients and tail bounds are estimated in global energy norms rather than pointwise over the physical domain $D$ [4]. As we will show, our choice of regularization enables us to prove convergence rates for global approximations in the same energy norms, with explicit computable coefficients. We present the theoretical support of this strategy via several extensions of compressed sensing concepts such as the restricted isometry property (RIP) and null space property (NSP) to the Hilbert-valued setting.

Our regularization problem is solved with a forward-backward splitting approach. We present a new strong convergence result for forward-backward splitting in a joint sparse recovery scenario where we assume neither the strict convexity of the fidelity functions nor the finite cardinality of the set of signals to construct. We note that most of similar strong convergence results along the line need at least one of these conditions. In relation to parameterized PDEs problem, this result means that forward-backward splitting strongly converges even before discretization in (infinite-dimensional) Hilbert space is introduced.

The rest of the paper proceeds as follows. In Section II, we introduce our approach of sparse regularization and provide a brief comparison to the problem of joint-sparse recovery. In Section III, we present theoretical guarantees for the direct reconstruction of Hilbert-valued functions through mixed norm regularization, as well as convergence estimates in case of parameterized PDEs. Section IV introduces a version of forward-backward algorithm adapted for Hilbert-valued recovery and provides its strong convergence result. Section V presents numerical results on applying the sparse regularization method to the solution of a parameterized elliptic PDE. Finally, Section VI concludes with some remarks on our approach and applications to a wider array of problems.

II Sparse regularization for parameterized PDEs

The sparse polynomial techniques proposed in this work are applicable to general parameterized PDE problems of the form: find $u(\cdot,\bm{y}):\overline{D}\times\mathcal{U}\to\mathbb{R}$ for all ${\bm{y}}\in\mathcal{U}$ such that

[TABLE]

where $\mathcal{L}$ is a differential operator defined on a spatio-temporal domain $D$. Our development to the general CS/joint-sparse polynomial approximation problem can be formulated as follows. One first generates $m$ samples ${\bm{y}}_{1},\ldots,{\bm{y}}_{m}$ in ${\cal U}$ independently from the orthogonalization measure $\varrho$ associated with $(\Psi_{\bm{\nu}})_{{\bm{\nu}}\in{\cal J}}$, for instance, uniform samples for the Legendre basis and Chebyshev samples for Chebyshev basis, see [9, Chapter 12], and solves the equation (3) at these samples to form the normalized output vector ${\bm{u}}:=(u({\bm{y}}_{1}),\ldots,u({\bm{y}}_{m}))/\sqrt{m}$ as well as the normalized sampling matrix $\bm{A}:=(\Psi_{\nu}({\bm{y}}_{i})/\sqrt{m})_{{{1\leq i\leq m},\,{\nu\in{\cal J}}}}$. Taking account that the true unknown coefficient ${\bm{c}}=({\bm{c}}_{{\bm{\nu}}})_{{\bm{\nu}}\in{\cal J}}$ approximately solves the linear system

[TABLE]

and further, ${\bm{c}}$ is compressible [4], it is reasonable to approximate ${\bm{c}}$ by ${\bm{c}}^{\#}$, the solution to the following problem

[TABLE]

where $\eta$ relates to an estimate of the tail (2); or the equivalent unconstrained convex minimization:

[TABLE]

for appropriately chosen $\mu>0$. Here, the norm $\|\cdot\|_{{\cal V},q}$ is defined for ${\bm{c}}\in{\cal V}^{N}$ as $\|{{\bm{c}}}\|_{{\cal V},q}:=(\sum_{{\bm{\nu}}\in{\cal J}}\|{\bm{c}}_{{\bm{\nu}}}\|^{q}_{{\cal V}})^{1/q}$. This is arguably the most natural extension of the $\ell_{1}$ minimization approach, traditionally for real and complex signal recovery, to the reconstruction of sparse generalized Hilbert-valued vectors. We denote our approach simultaneous compressed sensing (SCS). In depth description, analysis and application of SCS for solving parameterized PDEs are provided in [6].

Problem (5) can be related to the joint-sparse basis pursuit denoising problem as follows. Let $(\phi_{r})_{r\in\mathbb{N}}$ be an orthonormal basis of ${\cal V}$, then ${\bm{c}}_{{\bm{\nu}}}\in{\cal V}$ has unique representation

[TABLE]

Each coefficient ${\bm{c}}_{{\bm{\nu}}}$ corresponds to an $\mathbb{R}^{1\times\mathbb{N}}$ vector $(c_{{\bm{\nu}},1},c_{{\bm{\nu}},2},\ldots,c_{{\bm{\nu}},r},\ldots)^{\top}$, thus, ${\bm{c}}=({\bm{c}}_{{\bm{\nu}}})_{{\bm{\nu}}\in{\cal J}}$ is completely determined by the $\mathbb{R}^{N\times\mathbb{N}}$ matrix ${\widehat{\bm{c}}}=(c_{{\bm{\nu}},r})_{{\bm{\nu}}\in{\cal J},r\in\mathbb{N}}$. Furthermore,

[TABLE]

where the matrix norm $\|\cdot\|_{p,q}$ is defined as the $\ell_{q}$ norm of the vector $(\|\widehat{{\bm{c}}}_{{\bm{\nu}}}\|_{p})_{{\bm{\nu}}\in{\cal J}}$, implying the equivalence of (5) with the infinite-dimensional joint-sparse recovery problem. In practice, one needs to employ a discretization over $\mathcal{V}$ to be able to numerically solve (5) or (6). Any preferred method of spatial/temporal discretization may be used, e.g., the finite element, difference, and volume methods. However, our strong convergence result of forward-backward splitting in Section IV is applicable to (6) without any discretization.

III Error estimates for Hilbert-valued recovery

Straightforward extensions of concepts and results from compressed sensing and joint-sparse recovery can be made to ensure uniform recovery of Hilbert-valued signals via $\ell_{{\cal V},1}$-relaxation. Well-known concepts such as the null space property (NSP) and restricted isometry property (RIP) have Hilbert-valued counterparts. In this section, we review Hilbert-valued versions of the NSP and RIP to guarantee uniform recovery for the mixed norm regularization in the presence of noise or sparsity defects. We note that the extension of these concepts to the Hilbert-valued recovery setting does not require finite-dimensionality of the Hilbert space ${\cal V}$. For more detailed discussion and complete proofs of the results in this section, we refer the interested readers to [6].

First, we define a Hilbert-valued version of the $\ell_{2}$-robust NSP, which guarantees the reconstruction of vectors ${\bm{c}}\in\mathcal{V}^{N}$ (up to $s$ largest components and up to a noise level).

Definition 1** ($\ell_{{\cal V},2}$-robust null space property).**

The matrix ${\bm{A}}\in{\mathbb{R}}^{m\times N}$ is said to satisfy the $\ell_{{\cal V},2}$-robust null space property of order $s$ with constants $0<\rho<1$ and $\tau>0$ if

[TABLE]

An RIP-type condition is required to quantify the sample complexity of solving (5) to a prescribed accuracy. The following result establishes the implication of the $\ell_{{\cal V},2}$-robust NSP from the standard RIP.

Proposition 1**.**

Suppose that ${\cal V}$ is a separable Hilbert space, and that the matrix ${\bm{A}}\in{\mathbb{R}}^{m\times N}$ satisfies the RIP, that is

[TABLE]

with $\delta_{2s}<\frac{4}{\sqrt{41}}$. Then, ${\bm{A}}$ satisfies the $\ell_{{\cal V},2}$-robust NSP of order $s$ with constants $0<\rho<1$ and $\tau>0$ depending only on $\delta_{2s}$.

Proposition 1 implies that sample complexity results for solving the standard $\ell_{1}$ problem hold for the mixed norm problem (5) as well.

In the parameterized PDE context, with error estimates as in [19, Theorem 2] for quasi-optimal approximations, we are able to provide convergence rates for approximations to parameterized PDEs obtained through the mixed norm regularization. First, we assume that the solution $u$ has parametric expansion with coefficients $({\bm{c}}_{{\bm{\nu}}})_{{\bm{\nu}}\in{\cal F}}$ as in (1) satisfying $\|{\bm{c}}_{{\bm{\nu}}}\|_{\cal V}\lesssim e^{-b({\bm{\nu}})}$ for every ${\bm{\nu}}\in{\cal F}$, with $b({\bm{\nu}})$ obeying [19, Assumption 3]. For brevity, we do not detail that assumption herein, but remark that it is satisfied by most parameterized PDE models we are aware of.

Theorem 1**.**

For any $\varepsilon>0$, assume that the solution $u$ to (3) with parametric expansion (1) has coefficients $({\bm{c}}_{{\bm{\nu}}})_{{\bm{\nu}}\in{\cal F}}$ satisfying $\|{\bm{c}}_{\bm{\nu}}\|_{\cal V}\leq e^{-b({\bm{\nu}})}$ for every ${\bm{\nu}}\in{\cal F}$ with $b:[0,\infty)^{d}\to{\mathbb{R}}$ also satisfying [19, Assumption 3]. Denote by ${\cal J}_{s}$ the set of indices corresponding to the $s$ largest bounds of the sequence $(e^{-b({\bm{\nu}})})_{{\bm{\nu}}\in{\cal F}}$. Let ${\cal J}$ be such that ${\cal J}_{s}\subseteq{\cal J}$, and assume that the number of samples $m\in{\mathbb{N}}$ satisfies

[TABLE]

with $\Theta=\sup_{{\bm{\nu}}\in{\cal J}}\|\Psi_{\bm{\nu}}\|_{L^{\infty}({\cal U})}$ and $N=\#({\cal J})$. Then with probability $1-N^{-\log(s)}$, the solution ${\bm{c}}^{\#}$ of

[TABLE]

approximates $u$ with asymptotic error

[TABLE]

where $\kappa,{C}_{\varepsilon}>0$ are independent of $s$.

IV Forward-backward algorithm for Hilbert-valued recovery

In this section, we present a forward-backward splitting algorithm for solving problem (6) over the real Hilbert space ${\cal V}^{N}$. Assuming $\mu=1$ for simplicity, define $T_{1}$ and $T_{2}$ to be the subdifferential and gradient of the $\|\cdot\|_{{\cal V},1}$ and $\frac{1}{2}\|{\bm{A}}(\cdot)-{\bm{u}}\|_{{\cal V},2}^{2}$ parts of (6), respectively, and $T=T_{1}+T_{2}$. The algorithm can be derived as follows. Given $\tau>0$, we have

[TABLE]

The last identity in (12) leads to the forward-backward splitting algorithm: given initial guess ${\bm{x}}^{0}\in\mathcal{H}:=\{{\bm{x}}\in{\cal V}^{N}:\|{\bm{x}}\|_{{\cal V},2}<\infty\}$, compute

[TABLE]

where ${\bm{x}}^{k}$ denotes the approximation at $k$-th iterate. Letting

[TABLE]

then (13) can be written as ${\bm{x}}^{k+1}=S_{\tau}({\bm{x}}^{k})$. A straightforward derivation yields:

[TABLE]

One can observe that the forward operator $G_{\tau}$ resembles a step of gradient descent algorithm with stepsize $\tau$. The backward operator $J_{\tau}$, on the other hand, is a soft thresholding step associated with proximal point method.

Under a standard assumption on the stepsize $\tau$, namely $0<\tau<2/\|{\bm{A}}^{*}{\bm{A}}\|_{2}$, we can obtain some nonexpansive properties for $G_{\tau}$ and $J_{\tau}$. In particular, $G_{\tau}$ is nonexpansive, i.e.,

[TABLE]

and $J_{\tau}$ is row-wise firmly nonexpansive, i.e., $\forall{\bm{v}},{\bm{w}}\in{\cal H}$, $\forall j\in[N]:=\{1,2,\ldots,N\}$

[TABLE]

see [2, Chapter 4]. To prove the strong convergence of the forward-backward splitting algorithm for Hilbert-valued recovery in infinite-dimensional setting, we consider certain partition of the index set $[N]$ into two subsets $L$ and $E$, $L\subset(\text{supp}(\bm{x}^{*}))^{c}$ and $E\supset\text{supp}(\bm{x}^{*})$ where ${\bm{x}}^{*}$ is a solution to (6), as they require different treatments, see [11, 5]. With properties (16)-(17) in hand, our result is obtained in three steps:

Using a partitioning technique from [11], establish finite convergence ${\bm{x}}_{j}^{k}\to{\bm{x}}^{*}_{j}$ for every $j\in L$. 2. 2.

Establish an angular convergence result on the extended support set $E$. 3. 3.

Combine the known general weak convergence of the forward-backward algorithm, see, e.g., [2], with our angular convergence result to obtain convergence in norm on $E$, and hence strong convergence since $L\cup E=[N]$.

For more details on the above approach, see [5]. The established result can then be summarized as follows.

Theorem 2**.**

Let $0<\tau<2/\|{\bm{A}}^{*}{\bm{A}}\|_{2}$ and $({\bm{x}}^{k})$ be the sequence generated by the forward-backward iterations ${\bm{x}}^{k+1}=S_{\tau}({\bm{x}}^{k})$ starting from any ${\bm{x}}^{0}\in{\cal H}$. Then $({\bm{x}}^{k})$ converges strongly to an element ${\bm{x}}^{*}\in{\cal V}^{N}$ solving (6).

V Numerical experiments on parameterized elliptic PDE models

In this section we present numerical experiments demonstrating the efficiency of the proposed approach in approximating the solution of the following parameterized elliptic PDE: for all ${\bm{y}}\in{\cal U}$, find $u(\cdot,{\bm{y}}):D\to{\mathbb{R}}$ such that

[TABLE]

Here $D=[0,1]^{2}$, $f\equiv 1$, the diffusion coefficient is given by

[TABLE]

or its log transformed version $\log(a(x,{\bm{y}})-0.5)$, and $y_{i}$ are uniform random variables on $(-\sqrt{3},\sqrt{3})$. Convergence is compared against “highly-enriched” reference sparse grid stochastic collocation approximations which we denote $u_{h,\textnormal{ex}}(x,{\bm{y}})$, see [18]. All approximations (including the enriched reference solution) are computed on fixed finite element meshes $\mathcal{T}_{h}$, and our enriched SC approximation is computed using Clenshaw-Curtis (CC) abscissas with high level $L_{\textnormal{ex}}$. We compare performance using the relative errors of the expectation and standard deviation in the $H_{0}^{1}(D)$-norm.

In our plots and discussion, we use the following abbreviations. For the SCS method, we use “SCS-TD” to denote approximations obtained using the SCS method with the total degree (TD) subspace. For the stochastic Galerkin (SG) method, we use “SG-TD” to denote the SG approximation with TD subspace, see, e.g., [20, 10]. The stochastic collocation (SC) method with CC points and doubling growth rule is denoted “SC-CC,” see [14]. The Monte Carlo approximation is denoted “MC.”

We follow convention from [1], identifying stochastic degrees of freedom (SDOF) as the number of random sample points $m$ for MC and SCS, and sparse grid points $m_{L}$ for SC with level $L$. For SG, we define SDOF to be $N$, the cardinality of the index set used in construction. We include the SG method only to compare the $L^{2}({\cal U},d\varrho)$-optimal (w.r.t. SG SDOF) error of the Galerkin projection against the error of the sampling-based approximations.

We employ the orthonormal Legendre series for the parametric discretization of the SG and SCS methods. For the MC and SCS methods, the random samples $({\bm{y}}_{i})_{i=1}^{m}$ are drawn from the uniform distribution $\varrho$. We average the random sampling results over 24 trials, fixing the initial seed for the pseudorandom number generator on each trial, and then solving each trial’s problem with the same set of $m$ samples when plotting convergence.

Figures 1 & 2 display the achieved results. For highly anisotropic problems, the compressed-sensing based approach is able to naturally detect the underlying anisotropy in refinement. We note that, for the problem under consideration in Figure 2, we are unable to obtain an SG-TD approximation due to the difficulty of solving the nonlinearly coupled SG systems, see [7] for more details. We also expect that incorporation of anisotropic weighting schemes such as those considered in [1] and weighted $\ell_{1}$-minimization techniques such as those from [3] to improve results for the problems considered. We leave a further study of such improvements to a future work.

VI Concluding remarks

In this work we presented an overview of a novel sparse polynomial approximation technique, enabling global recovery of solutions to parameterized PDEs. Our approach, based on extensions of compressed sensing and joint-sparse recovery, treats the solution vector as an element of a tensor product of real Hilbert spaces ${\cal V}$. The key difference in our technique is the use of a mixed norm involving the energy norm of the associated PDE problem and the standard vector $\ell_{1}$ norm. Within this framework, we are able to prove uniform recovery results through straightforward extensions of concepts such as the restricted isometry and null space properties. Moreover, by combining extensions of error estimates for the standard basis pursuit denoising problem with quasi-optimal error estimates for sparse approximation of parameterized PDE systems, we are able to derive sub-exponential convergence of our method. For more details see [6].

We have also presented a modification of the standard forward-backward splitting algorithm for Hilbert-valued recovery. As the considered convex optimization problem is posed over an infinite dimensional Hilbert space, new techniques were needed to establish its convergence properties. By deriving a novel angular convergence result from the firm nonexpansiveness of the soft-thresholding operator, we are able to prove the strong convergence of the algorithm in the considered setting.

Finally, we presented numerical results on the application of our approach to the approximation of solutions to both affinely and non-affinely parameterized elliptic PDEs. We compare our results with those obtained using the stochastic Galerkin, stochastic collocation, and Monte Carlo methods. The achieved results are positive, highlighting a key benefit of compressed sensing-based approaches, namely the ability to detect underlying problem anisotropy in refinement.