Wavelet-Fourier CORSING techniques for multi-dimensional advection-diffusion-reaction equations

TL;DR

This paper introduces a wavelet-Fourier CORSING method for efficiently solving multi-dimensional advection-diffusion-reaction equations, combining compressed sensing with Petrov-Galerkin techniques to achieve accurate coefficient recovery.

Contribution

It provides the first rigorous error bounds and implementation strategies for multi-dimensional wavelet-Fourier CORSING methods in PDEs.

Findings

Achieves accurate approximation of PDE solutions using sparse recovery.

Provides stability and robustness demonstrated through numerical experiments.

Develops new estimates for local a-coherence in multi-dimensional settings.

Abstract

We present and analyze a novel wavelet-Fourier technique for the numerical treatment of multidimensional advection-diffusion-reaction equations based on the CORSING (COmpRessed SolvING) paradigm. Combining the Petrov-Galerkin technique with the compressed sensing approach, the proposed method is able to approximate the largest coefficients of the solution with respect to a biorthogonal wavelet basis. Namely, we assemble a compressed discretization based on randomized subsampling of the Fourier test space and we employ sparse recovery techniques to approximate the solution to the PDE. In this paper, we provide the first rigorous recovery error bounds and effective recipes for the implementation of the CORSING technique in the multi-dimensional setting. Our theoretical analysis relies on new estimates for the local a-coherence, which measures interferences between wavelet and Fourier basis functions with respect to the metric induced by the PDE operator. The stability and robustness of the proposed scheme is shown by numerical illustrations in the one-, two-, and three-dimensional case.

Equations

$$\mathcal{D}=(0,1{)}^n\subseteq {\mathbb{R}}^n.$$

$$(u,v{)}_k:=\sum _{\begin{array}{c}\mathbf{\alpha }\in [k{]}_0^n\\ {\alpha }_1+\cdots +{\alpha }_n\le k\end{array}}{\int }_{\mathcal{D}}{D}^{\mathbf{\alpha }}u(\mathbf{x})\overline{D^{\mathbf{\alpha }}v(\mathbf{x})}\text{d}\mathbf{x},$$

$$|u{|}_{H^k(\mathcal{D})}^2:=\sum _{\begin{array}{c}\mathbf{\alpha }\in [k{]}_0^n\\ {\alpha }_1+\cdots +{\alpha }_n=k\end{array}}\Vert D^{\mathbf{\alpha }}u{\Vert }_{L^2(\mathcal{D})}^2.$$

$$H_{\text{per}}^k(\mathcal{D}):={\mathrm{clos}}_{\Vert \cdot {\Vert }_{H^k(\mathcal{D})}}(C_{\text{per}}^{\infty }(\mathcal{D})),$$

$$C_{\text{per}}^{\infty }(\mathcal{D}):=\{v{|}_{\mathcal{D}}:v\in C^{\infty }({\mathbb{R}}^n),v(\mathbf{x}+{\mathbf{e}}_i)=v(\mathbf{x}),\forall \mathbf{x}\in {\mathbb{R}}^n,\forall i\in [n]\}.$$

$$u(\mathbf{x})=\sum _{\mathbf{r}\in {\mathbb{Z}}^n}{c}_{\mathbf{r}}{e}^{2\pi \text{i}\mathbf{r}\cdot \mathbf{x}},\text{with }{c}_{\mathbf{r}}:={\int }_{\mathcal{D}}u(\mathbf{x})e^{-2\pi \text{i}\mathbf{r}\cdot \mathbf{x}}\text{d}\mathbf{x}.$$

$$\Vert u{\Vert }_{H^k(\mathcal{D})}^2\sim \sum _{\mathbf{r}\in {\mathbb{Z}}^n}(1+|\mathbf{r}{|}^{2k})|c_{\mathbf{r}}{|}^2,\forall u\in H_{\text{per}}^k(\mathcal{D}),\forall k\in {\mathbb{N}}_0.$$

$$a(u,v):={\int }_{\mathcal{D}}\left[\eta (\mathbf{x})\nabla u(\mathbf{x})\cdot \overline{\nabla v(\mathbf{x})}+\mathbf{\beta }(\mathbf{x})\cdot \nabla u(\mathbf{x})\overline{v(\mathbf{x})}+\rho (\mathbf{x})u(\mathbf{x})\overline{v(\mathbf{x})}\right]\text{d}\mathbf{x},$$

$$\text{find }u\in H_{\text{per}}^1(\mathcal{D}):a(u,v)=(f,v),\forall v\in H_{\text{per}}^1(\mathcal{D}),$$

$$|a(u,u)|\ge \alpha \Vert u{\Vert }_{H^1(\mathcal{D})}^2,\forall u\in H_{\text{per}}^1(\mathcal{D}),$$

$$|a(u,v)|\le \mathcal{A}\Vert u{\Vert }_{H^1(\mathcal{D})}\Vert v{\Vert }_{H^1(\mathcal{D})},\forall u,v\in H_{\text{per}}^1(\mathcal{D}).$$

$$\eta \ge {\eta }_{\min },-\frac{1}{2}\nabla \cdot \mathbf{\beta }+\rho \ge \zeta ,\text{a.e. in }\mathcal{D}.$$

$$\Vert u{\Vert }_{H^1(\mathcal{D})}\le \frac{1}{\alpha }\Vert f{\Vert }_{H^{-1}(\mathcal{D})},$$

$$\mathcal{A}\le \max \left\{\Vert \eta {\Vert }_{L^{\infty }(\mathcal{D})},\underset{\mathbf{x}\in \mathcal{D}}{\sup }\Vert \mathbf{\beta }(\mathbf{x}){\Vert }_2,\Vert \rho {\Vert }_{L^{\infty }(\mathcal{D})}\right\}.$$

$$\Psi :={\Phi }_{{\ell }_0}\cup \bigcup _{\ell ={\ell }_0}^{L-1}{\Psi }_{\ell },$$

$${\Psi }^{\text{per}}:={\Phi }_{{\ell }_0}^{\text{per}}\cup \bigcup _{\ell ={\ell }_0}^{L-1}{\Psi }_{\ell }^{\text{per}},$$

$$\mathbb{Z}/(2^{\ell }\mathbb{Z})\equiv \{0,1,\dots ,2^{\ell }-1\},\forall \ell \in {\mathbb{N}}_0.$$

$${\psi }_{{\ell }_0-1,k}\equiv {\phi }_{{\ell }_0,k},\forall k\in \mathbb{Z}.$$

$${\mathbf{a}}_{[-1:1]}=[\frac{1}{2},1,\frac{1}{2}],{\tilde{\mathbf{a}}}_{[-2:2]}=[-\frac{1}{4},\frac{1}{2},\frac{3}{2},\frac{1}{2},-\frac{1}{4}],$$

$${\mathbf{b}}_{[-1:3]}=[\frac{1}{4},\frac{1}{2},-\frac{3}{2},\frac{1}{2},\frac{1}{4}],{\tilde{\mathbf{b}}}_{[0:2]}=[\frac{1}{2},-1,\frac{1}{2}].$$

$${\psi }_{\ell ,k}^{\text{per}}{|}_{\mathrm{supp}({\psi }_{\ell ,k})}\equiv {\psi }_{\ell ,k}{|}_{\mathrm{supp}({\psi }_{\ell ,k})},\forall \ell \ge {\ell }_0-1,\forall k\in \mathbb{Z}/(2^{\ell }\mathbb{Z}).$$

$${\psi }_{\mathbf{\ell },\mathbf{k}}^{\text{ani}}:={\psi }_{{\ell }_1,k_1}\otimes \cdots \otimes {\psi }_{{\ell }_n,k_n},$$

$$\mathbb{Z}/(2^{\mathbf{\ell }}\mathbb{Z}):=\mathbb{Z}/(2^{{\ell }_1}\mathbb{Z})\times \cdots \times \mathbb{Z}/(2^{{\ell }_n}\mathbb{Z})\equiv \prod _{j=1}^n\{0,1,\dots ,2^{{\ell }_j-1}\}.$$

$${\mathcal{J}}^{\text{ani}}:=\{(\mathbf{\ell },\mathbf{k}):\mathbf{\ell }\in {\mathbb{N}}^n,{\ell }_0-1\le \mathbf{\ell }<L,\mathbf{k}\in \mathbb{Z}/(2^{\mathbf{\ell }}\mathbb{Z})\},$$

$${\Psi }^{\text{ani}}=\{{\psi }_{\mathbf{j}}^{\text{ani}}{\}}_{\mathbf{j}\in {\mathcal{J}}^{\text{ani}}},$$

$${\vartheta }_{\ell ,k,e}:=\{\begin{array}{ll}{\phi }_{\ell ,k},& \text{if }e=0,\\ {\psi }_{\ell ,k},& \text{if }e=1.\end{array}$$

$${\psi }_{\ell ,\mathbf{k},\mathbf{e}}^{\text{iso}}:={\vartheta }_{\ell ,k_1,e_1}\otimes \cdots \otimes {\vartheta }_{\ell ,k_n,e_n}.$$

$${\mathcal{J}}^{\text{iso}}:=\{(\ell ,\mathbf{k},\mathbf{e}):\ell \in \mathbb{N},{\ell }_0\le \ell <L,\mathbf{k}\in (\mathbb{Z}/(2^{\ell }\mathbb{Z}){)}^n,\mathbf{e}\in \{0,1{\}}^n\},$$

$${\Psi }^{\text{iso}}:=\{{\psi }_{\mathbf{j}}^{\text{iso}}{\}}_{\mathbf{j}\in {\mathcal{J}}^{\text{iso}}},$$

$$\bigg{\|}\sum_{{\boldsymbol{j}}\in\mathcal{J}^{\star}}c_{{\boldsymbol{j}}}\psi_{{\boldsymbol{j}}}^{\star}\bigg{\|}_{H^{s}(\mathcal{D})}\sim\|D_{\star}^{s}{\boldsymbol{c}}\|_{2},\quad\forall\star\in\{\textnormal{ani},\textnormal{iso}\},$$

$$(D_{\star }{)}_{\mathbf{j},{\mathbf{j}}^{\prime }}=\{\begin{array}{ll}{2}^{\Vert \mathbf{\ell }{\Vert }_{\infty }}{\delta }_{\mathbf{j},{\mathbf{j}}^{\prime }},& \text{if }\star =\text{ani},\\ 2^{\ell }{\delta }_{\mathbf{j},{\mathbf{j}}^{\prime }},& \text{if }\star =\text{iso},\end{array}\forall \mathbf{j},{\mathbf{j}}^{\prime }\in {\mathcal{J}}^{\star },$$

Full text

1 Introduction

This paper deals with the theoretical analysis and the numerical implementation of a recently-introduced paradigm in numerical approximation of Partial Differential Equations (PDEs), named COmpRessed SolvING (in short, CORSING). The CORSING method has been proposed and studied in [8, 11, 12] for the solution of linear PDEs set in Hilbert spaces, and combines the Petrov-Galerkin discretization techniques with compressed sensing [14, 20]. Assuming the sparsity of the solution with respect to a suitable trial function basis, the idea is to build a reduced Petrov-Galerkin discretization of the weak formulation of the problem by considering a randomly subsampled test subspace, and then to recover a sparse approximation to the solution via sparse recovery techniques, such as $\ell^{1}$-minimization or the greedy Orthogonal Matching Pursuit (OMP) algorithm. As discussed in [12], the main advantages of CORSING with respect to adaptive finite element techniques for PDEs are that no a posteriori error estimators are needed and that assembly and recovery via OMP are fully parallelizable.

In this paper, we focus on the scenario of multi-dimensional advection-diffusion-reaction (ADR) equations on a torus, and we employ biorthogonal wavelets as trial functions basis and a Fourier basis as test functions. In this way, we introduce a hybrid Wavelet-Fourier technique named CORSING $\mathcal{WF}$, which is able to approximate the largest wavelet coefficients of the solution to the PDE by sampling randomly the Fourier test space, using a suitable probability measure.

It is worth noticing that the results showed here can be extended to the nonperiodic case by considering biorthogonal wavelets on the interval (and on the hypercube) (see [18, 25, 33]) as the trial functions basis and the sine basis $\{\sin(\pi kx)\}_{k\in\mathbb{N}}$ (or a tensorized version of it) as the test functions basis. However, we decided to stick to the periodic case to make the theoretical exposition free of an excessive quantity of technical details regarding the construction of wavelets at the boundary. In this respect, the present work should be considered as a first step towards the setup of CORSING in practical applications. We also observe that the construction of wavelets over general domains is not a straightforward task. Hence, the CORSING $\mathcal{WF}$ shares this difficulty with classic wavelets methods for PDEs. We refer to [13] for an extension and implementation of the CORSING method to domains with more general geometries via the IsoGeometric Analysis principle [22].

The choice of biorthogonal wavelet functions is motivated by the need of working with trial and test functions bases that satisfy the Riesz basis property [12]. Roughly speaking, this assumption guarantees a control on the condition number of the Petrov-Galerkin discretization matrix, which is crucial for a successful application of the compressed sensing paradigm (see, e.g., [8, Theorem 1.21] or [2, Theorem 3.6]). In [12], the authors considered the hierarchical basis of hat functions and the sine functions basis as trial and test basis, respectively. While this choice ensures the Riesz basis property in the one-dimensional setting, tensorizing these two bases breaks the Riesz basis property in the multi-dimensional framework. Indeed, the tensorization of the sine functions basis is a Riesz basis (up to suitable diagonal rescaling) in any dimension, while the Riesz property is broken for the tensorized hierarchical basis of hat functions. However, thanks to the so-called norm equivalence property, tensorized biorthogonal wavelets form a Riesz basis in any dimension, hence providing a remedy for this issue.

Although compressed sensing is becoming a standard paradigm for signal processing applications, understanding its full potential and limitations in scientific computing is still object of active work. This paper moves a step forward in this direction. In a fast-growing literature, it is worth mentioning here the applications of compressed sensing to numerical approximation of high-dimensional functions (see, e.g., [3, 4, 15, 28]) and of parametric PDEs, with special emphasis on uncertainty quantification (see, e.g., [7, 21, 27, 34]). In these cases, a smooth function, which can be the quantity of interest of a parametric PDE, is approximated with respect to a global sparsity basis like orthogonal polynomials by means of random pointwise observations. The PDE solver is a black box used to evaluate the quantity of interest for different values of the parameters, and compressed sensing is performed outside the black box. Our focus is different, since the CORSING method takes advantage of the compressed sensing paradigm inside the black box, i.e., to solve the PDE itself given a particular choice of the parameters.

The compressed sensing principle has also been recently employed for the efficient numerical approximation of diffusion equations via Sturm-Liouville spectral collocation in [9]. Finally, we observe that wavelet-Fourier techniques are widely used in signal processing applications (see, e.g., [5, 6, 23] for theoretical contributions in this direction). Yet, to the best of our knowledge, this paper is the first comprehensive study of this type of techniques in the context of numerical approximation of multidimensional PDEs.

1.1 Main contributions

The main contribution of this paper is threefold. First, we present and study the first hybrid wavelet-Fourier discretization for ADR equations in arbitrary dimension on the torus based on Petrov-Galerkin discretization and on compressed sensing, named CORSING $\mathcal{WF}$ (see Section 3). Moreover, in Section 4 we show the applicability of the theoretical analysis in [12] to $n$-dimensional ADR equations with constant coefficients for $n\geq 1$ and with nonconstant coefficients for $n=1$. Finally, in Section 5 we provide a Matlab${}^{\text{\textregistered}}$ implementation of the CORSING $\mathcal{WF}$ method for $n$-dimensional ADR equations, with $n=1,2,3$.

In view of the aforementioned contributions, we will focus on the following three key technical issues necessary to implement and quantify the performance of the CORSING $\mathcal{WF}$ method (more details on these three issues are given in Section 3.4):

• (i)

Find a suitable truncation condition on the Fourier test space in order to guarantee stability of the resulting Petrov-Galerkin discretization.

• (ii)

Give lower bounds to the sampling complexity, i.e., the minimum number of randomly selected Fourier test functions needed to recover the $s$ dominant coefficients of the solution in the wavelet expansion.

• (iii)

Provide explicit expressions for the probability distribution on the Fourier test space needed for the random selection of the basis functions.

In order to address issues (i), (ii), and (iii), we will take advantage of the general framework given in [12] for the analysis of CORSING, based on the so-called local $a$-coherence, which can be interpreted as a measure of the angle between the wavelet trial functions and the Fourier test functions with respect to the metric induced by the sesquilinear form associated with the ADR problem (see Definition 3.3). In view of this, our main efforts will be aimed at producing upper bounds to the local $a$-coherence in Section 4 for the specific ADR problems addressed. In particular, we provide local $a$-coherence upper bounds for the 1D case with nonconstant coefficients in Theorem 4.2, and for the multi-dimensional case with constant coefficients, when employing anisotropic and isotropic tensor product wavelets in Theorems 4.7 and 4.12, respectively. As a consequence, we derive explicit and computable answers to (i), (ii), and (iii) and provide recovery error guarantees for the CORSING $\mathcal{WF}$ method in Theorem 4.4 (1D case), Theorem 4.9 (multi-dimensional anisotropic case), and Theorem 4.12 (multi-dimensional isotropic case). Numerical results in Section 5 confirm the theoretical findings.

2 Problem setting

In this section, we recall some basics on periodic Sobolev spaces and discuss ADR problems in weak form with periodic boundary conditions.

Notation

We denote by $\mathbb{N}:=\{1,2,3,\ldots\}$, $\mathbb{N}_{0}:=\{0\}\cup\mathbb{N}$, and $\mathbb{Z}:=\mathbb{N}_{0}\cup(-\mathbb{N})$. We define $[k]:=\{1,\ldots,k\}$ and $[k]_{0}:=\{0\}\cup[k]$. The notation $X\lesssim Y$ means $X\leq CY$, with $C>0$ a constant independent of $X$ and $Y$; $X\sim Y$ means that $X\lesssim Y$ and $X\gtrsim Y$ hold simultaneously. By $X\propto Y$, we understand that there exists a constant $C>0$ independent of $X$ and $Y$ such that $X=CY$. Given a multi-index ${\boldsymbol{r}}$, we denote its 2-norm by $|{\boldsymbol{r}}|$. For every $z\in\mathbb{C}$, $|z|$ is its modulus, $\overline{z}$ is its complex conjugate, and $\operatorname{Re}(z)$ and $\operatorname{Im}(z)$ denote its real and imaginary part, respectively. Given a vector ${\boldsymbol{x}}\in\mathbb{C}^{n}$ and $1\leq p<+\infty$, then $\|{\boldsymbol{x}}\|_{p}:=(\sum_{j=1}^{n}|x_{j}|^{p})^{1/p}$, $\operatorname{supp}({\boldsymbol{x}}):=\{j\in[n]:x_{j}\neq 0\}$, and $\|{\boldsymbol{x}}\|_{0}:=|\operatorname{supp}({\boldsymbol{x}})|$. Inequalities between vectors in $\mathbb{R}^{n}$ have to be read componentwise; for example, ${\boldsymbol{x}}\leq{\boldsymbol{y}}$ means $x_{i}\leq y_{i}$ for every $i$. The vectors of the canonical basis of $\mathbb{C}^{n}$ are denoted by ${\boldsymbol{e}}_{1},\ldots,{\boldsymbol{e}}_{n}$ and ${\boldsymbol{x}}\cdot{\boldsymbol{y}}:=x_{1}\overline{y}_{1}+\cdots+x_{n}\overline{y}_{n}$ is the standard inner product of $\mathbb{C}^{n}$. Given a matrix $A\in\mathbb{C}^{m\times n}$, $A^{*}$ denotes its conjugate transpose. The set of sequences indexed by integer multi-indices is denoted by $\ell(\mathbb{Z}^{n}):=\{(x_{{\boldsymbol{j}}})_{{\boldsymbol{j}}\in\mathbb{Z}^{n}}:x_{{\boldsymbol{j}}}\in\mathbb{C},\,\forall{\boldsymbol{j}}\in\mathbb{Z}^{n}\}$.

2.1 Sobolev spaces

We start by recalling some standard notions about periodic Sobolev spaces. Let $n\in\mathbb{N}$ and consider the domain

[TABLE]

Given $k\in\mathbb{N}_{0}$, let $H^{k}(\mathcal{D})=H^{k}(\mathcal{D};\mathbb{C})$ be the Sobolev space of order $k$ of complex-valued functions over $\mathcal{D}$, being understood that $H^{0}(\mathcal{D})=L^{2}(\mathcal{D})$. Moreover, we denote the $H^{k}(\mathcal{D})$-inner product by

[TABLE]

where $D^{{\boldsymbol{\alpha}}}:=\frac{\partial^{\alpha_{1}+\cdots+\alpha_{n}}}{\partial x_{1}^{\alpha_{1}}\cdots\partial x_{n}^{\alpha_{n}}}$ is the derivative in the sense of distributions and $[k]_{0}^{n}:=[k]_{0}\times\cdots\times[k]_{0}$ ($n$ times). For the sake of simplicity, we use the shorthand notation $(\cdot,\cdot):=(\cdot,\cdot)_{0}$ for the $L^{2}(\mathcal{D})$-inner product. The $H^{k}(\mathcal{D})$-norm is defined as $\|\cdot\|_{H^{k}(\mathcal{D})}^{2}=(\cdot,\cdot)_{k}$ and the $H^{k}(\mathcal{D})$-seminorm is given by

[TABLE]

The periodic Sobolev space of order $k$, $H^{k}_{\textup{per}}(\mathcal{D})\subseteq H^{k}(\mathcal{D})$, is then defined as

[TABLE]

where

[TABLE]

Notice that $H^{0}_{\textup{per}}(\mathcal{D})\equiv L^{2}(\mathcal{D})$ since $C^{\infty}_{\textup{per}}(\mathcal{D})$ is dense in $L^{2}(\mathcal{D})$. Now, let $u\in H^{k}_{\textup{per}}(\mathcal{D})$ and consider its Fourier series expansion

[TABLE]

Then, the following norm equivalence holds:

[TABLE]

Moreover, we define $H^{-1}_{\textup{per}}(\mathcal{D}):=[H^{1}_{\textup{per}}(\mathcal{D})]^{*}$, where the superscript ∗ denotes the dual space. For the proofs of these results and for more details about periodic Sobolev spaces, we refer the reader to [1, 31, 32].

2.2 Advection-diffusion-reaction problems

Consider the sesquilinear form $a:H^{1}_{\textup{per}}(\mathcal{D})\times H^{1}_{\textup{per}}(\mathcal{D})\to\mathbb{C}$ defined as

[TABLE]

where $\eta,\rho:\mathcal{D}\to\mathbb{R}$ are the diffusion and reaction coefficients, respectively, and ${\boldsymbol{\beta}}:\mathcal{D}\to\mathbb{R}^{n}$ is the advective field. Then, the weak formulation of the periodic ADR equation reads

[TABLE]

where $f:\mathcal{D}\to\mathbb{R}$ is a forcing term. Although we are interested in the case of real-valued coefficients $\eta$, $\rho$, ${\boldsymbol{\beta}}$, and $f$, the bilinear form $a(\cdot,\cdot)$ is defined over complex-valued Hilbert spaces to allow us the use of the Fourier basis in Section 3.

We recall that the coercivity constant of $a(\cdot,\cdot)$ is the largest $\alpha>0$ such that

[TABLE]

and the continuity constant of $a(\cdot,\cdot)$ is the smallest constant $\mathcal{A}>0$ such that

[TABLE]

In the following proposition, we provide conditions on the coefficients $\eta,{\boldsymbol{\beta}},\rho$ sufficient to ensure the well-posedness of problem (2.5). A proof of this result is provided in [10].

Proposition 2.1** (Well-posedness of the periodic ADR problem).**

Let $\eta,\rho\in L^{\infty}(\mathcal{D})$ and ${\boldsymbol{\beta}}\in[L^{\infty}(\mathcal{D})]^{n}$ be such that ${\boldsymbol{\beta}}$ is one-periodic with respect to each variable. Moreover, assume that there exist two constants $\eta_{\min},\zeta>0$ such that

[TABLE]

Then, for every $f\in H^{-1}_{\textup{per}}(\mathcal{D})$, the weak problem (2.5) admits a unique solution $u$ such that

[TABLE]

with $\alpha\geq\min\{\eta_{\min},\zeta\}$ the coercivity constant of $a(\cdot,\cdot)$. Moreover, the continuity constant of $a(\cdot,\cdot)$ satisfies

[TABLE]

Remark 2.2* (Diffusion equation).*

Notice that Proposition 2.1 does not encompass the purely diffusive case, i.e. ${\boldsymbol{\beta}}\equiv{\boldsymbol{0}}$ and $\rho\equiv 0$. Indeed, in this case we only have uniqueness of the solution up to constants. This issue can be fixed by assuming $u,v\in H^{1}_{\textup{per}}(\mathcal{D})/\mathbb{R}$ in the weak formulation (2.5). $\blacksquare$

3 CORSING Wavelet-Fourier

In order to solve problem (2.5), we describe the CORSING $\mathcal{WF}$ (Wavelet-Fourier) method.

3.1 Trial functions: wavelets

In this section, we present the biorthogonal wavelets on the periodic interval and on the periodic hypercube that will be employed as trial functions of the Petrov-Galerkin discretization. We refer to [10] and to [18, 25, 33] for technical details about wavelet construction. A crucial property of the biorthogonal wavelets is that, after suitable normalization, they are a Riesz basis for $H^{1}_{\textup{per}}(\mathcal{D})$, thanks to the so-called norm equivalence property (see Theorem 3.2).

Biorthogonal wavelets on the periodic interval.

Given $\ell_{0},L\in\mathbb{N}_{0}$ with $\ell_{0}<L$, we consider biorthogonal B-spline wavelet basis on the real line

[TABLE]

where $\Phi_{\ell}=\{\varphi_{\ell,k}\}_{k\in\mathbb{Z}}$ are scaling functions and $\Psi_{\ell}=\{\psi_{\ell,k}\}_{k\in\mathbb{Z}}$ are wavelet functions (for more details, see Figure 1, Setting 3.1 and [10]).

Notice that the dependence of $\Psi$ on the levels $\ell_{0}$ and $L$ is understood. For $L=\infty$, $\Psi$ is a basis for $L^{2}(\mathbb{R})$. In order to build a basis for $L^{2}(\mathcal{D})$, we resort to periodization (see [10]). Similarly, we denote by

[TABLE]

where $\Phi_{\ell}^{\textup{per}}=\{\varphi_{\ell,k}^{\textup{per}}\}_{k\in\mathbb{Z}/(2^{\ell}\mathbb{Z})}$ are periodized scaling functions and $\Psi_{\ell}^{\textup{per}}=\{\psi_{\ell,k}^{\textup{per}}\}_{k\in\mathbb{Z}/(2^{\ell}\mathbb{Z})}$ are periodized wavelet functions. In particular, $\mathbb{Z}/(2^{\ell}\mathbb{Z})$ denotes the ring of integers modulo $2^{\ell}$ that coincides with the set of canonical representatives, i.e.,

[TABLE]

Assuming the coarsest level $\ell_{0}\in\mathbb{N}_{0}$ to be fixed, we also introduce the following notation:

[TABLE]

Setting 3.1** (1D biorthogonal B-spline wavelets [17]).**

We consider 1D biorthogonal B-spline wavelets of order $(d,\widetilde{d})=(2,2)$, corresponding to primal and dual filters

[TABLE]

Moreover, we assume scaling functions $\Phi_{\ell}^{\text{per}}$ and wavelets $\Psi_{\ell}^{\text{per}}$ to be normalized with respect to the $L^{2}(\mathcal{D})$-norm (i.e., $\|\varphi_{\ell,k}\|_{L^{2}(\mathcal{D})}\sim\|\psi_{\ell,k}\|_{L^{2}(\mathcal{D})}\sim 1$, for every $\ell\in\mathbb{N}_{0}$ and $k\in\mathbb{Z}/(2^{\ell}\mathbb{Z})$) and $\ell_{0}\geq 2$.

When wavelets are built as in Setting 3.1, the corresponding $\Psi^{\textup{per}}$ is a basis for $H^{1}_{\textup{per}}(\mathcal{D})$ when $L=\infty$. In fact, for $L\in\mathbb{N}$, $\operatorname{span}(\Psi^{\textup{per}})$ coincides with the space of functions in $H^{1}_{\textup{per}}(\mathcal{D})$ that are continuous and piecewise linear on the uniform grid $2^{-L}\mathbb{Z}\cap[0,1]$. Moreover, choosing $\ell_{0}\geq 2$ guarantees that555In particular, the assumption $\ell_{0}\geq 2$ ensures that the periodization $\psi_{\ell,k}^{\textup{per}}$ of $\psi_{\ell,k}$ (defined as $\psi_{\ell,k}^{\textup{per}}(x)=2^{\ell/2}\sum_{j\in\mathbb{Z}}\psi(2^{\ell}(x+j)-k)$, for every $x\in\mathbb{R}$ – see also [10]) is the sum of terms with disjoint support. On the contrary, let us assume, for example, $\ell_{0}=1$. Then, since $\operatorname{supp}(\psi_{1,1})\subseteq(0,3/2)$ the periodization $\psi_{1,1}^{\textup{per}}$ is built by summing two overlapping terms on $[0,1/2]$. In particular, $\psi_{1,1}^{\textup{per}}$ is constant over $[0,1/2]$, whereas $\psi_{1,1}$ is not. This shows that condition (3.5) is not satisfied for $\ell_{0}=1$.

[TABLE]

For the sake of simplicity, the apex ${}^{\textup{per}}$ will be omitted in the subsequent developments. In fact, from now on we will always assume to be in the periodic setting.

In order to generalize this construction to arbitrary dimension $n>1$, we consider anisotropic and isotropic tensor product wavelets (see Figure 2). Here we just recall the basic definitions and we refer the reader to [18, 25, 33] for a more detailed discussion.

Anisotropic tensor product wavelets.

Given a $\ell_{0}\in\mathbb{N}$ and multi-indices ${\boldsymbol{\ell}}\in\mathbb{N}^{n}$, with ${\boldsymbol{\ell}}\geq\ell_{0}-1$ and ${\boldsymbol{k}}\in\mathbb{Z}^{n}$, a first straightforward way to obtain a multi-dimensional basis is by tensorizing the 1D wavelet basis $\Psi$ with itself $n$ times, namely,

[TABLE]

with $\psi_{\ell,k}=\psi_{\ell,k}^{\textup{per}}$. For every level multi-index ${\boldsymbol{\ell}}\in\mathbb{N}^{n}$, the spatial multi-index ${\boldsymbol{k}}$ takes values in $\mathbb{Z}/(2^{\boldsymbol{\ell}}\mathbb{Z})$, where

[TABLE]

Therefore, fixing $L\in\mathbb{N}$ with $L>\ell_{0}$ and defining the multi-index set

[TABLE]

we have

[TABLE]

where we set ${\boldsymbol{j}}=({\boldsymbol{\ell}},{\boldsymbol{k}})$.

Isotropic tensor product wavelets.

An alternative (and less straightforward) way to define a multi-dimensional wavelet basis is by isotropic tensorization. In this case, we need the following auxiliary notation:

[TABLE]

Then, given $\ell\in\mathbb{N}$ with $\ell\geq\ell_{0}$, ${\boldsymbol{k}}\in\mathbb{Z}^{n}$, and ${\boldsymbol{e}}\in\{0,1\}^{n}$, we define

[TABLE]

Then, fixing $L\in\mathbb{N}$ with $L>\ell_{0}$ and defining the multi-index set

[TABLE]

we have

[TABLE]

where we set ${\boldsymbol{j}}=(\ell,{\boldsymbol{k}},{\boldsymbol{e}})$.

The main difference between anisotropic and isotropic tensor product structure is visualized in Figures 2 and 3 (in the 2D case).

In particular, the anisotropic tensorization blends all the dyadic levels together, whereas the isotropic tensorization only combines basis functions of the same level.

Riesz basis property.

The biorthogonal wavelets employed here satisfy the norm equivalence property, presented in the following theorem. We refer to [10] for the proof.

Theorem 3.2** (Norm equivalence).**

Let $\Psi\in\{\Psi^{\textnormal{ani}},\Psi^{\textnormal{iso}}\}$ be the tensor product wavelet basis defined as in (3.8) or (3.12) from periodized 1D biorthogonal B-spline wavelets of order $(d,\widetilde{d})$. Then, the following norm equivalence holds:

[TABLE]

for every $s\in[0,d-1]$, where ${\boldsymbol{c}}=({\boldsymbol{c}})_{{\boldsymbol{\ell}},{\boldsymbol{k}}}$ and $D_{\star}\in\{D_{\textnormal{ani}},D_{\textnormal{iso}}\}$ is the diagonal matrix defined by

[TABLE]

where $\delta_{{\boldsymbol{j}},{\boldsymbol{j}}^{\prime}}$ is Kronecker’s delta function and where $\mathcal{J}^{\textnormal{ani}}$ and $\mathcal{J}^{\textnormal{iso}}$ are defined as in (3.7) and (3.11), respectively.

In Setting 3.1, Theorem 3.2 implies that a suitable weighted $\ell^{2}$-norm of the wavelet coefficients is equivalent to the $H^{1}(\mathcal{D})$-norm. In particular, using the fact that

[TABLE]

we can rescale the basis functions of $\Psi^{\textnormal{ani}}$ and $\Psi^{\textnormal{iso}}$ and normalize them with respect to the $H^{1}(\mathcal{D})$-norm and obtain two Riesz bases with respect to the $H^{1}(\mathcal{D})$-norm. More details on the norm equivalence property are provided in [10].

3.2 Test functions: Fourier basis

Consider the one-periodic Fourier basis functions in dimension one and their tensorized version in the $n$-dimensional case

[TABLE]

It is easy to verify that $\{\xi_{\boldsymbol{q}}:{\boldsymbol{q}}\in\mathbb{Z}^{n}\}$ is an orthonormal basis for $L^{2}(\mathcal{D})$ and that its elements are orthogonal with respect to the $H^{1}(\mathcal{D})$-inner product. Moreover, we have

[TABLE]

Given $R\in\mathbb{N}$, let us consider the following finite multi-index set:

[TABLE]

Then, we define the Fourier basis as

[TABLE]

In particular, $\Xi$ is a Riesz basis with respect to the $H^{1}(\mathcal{D})$-norm.

3.3 The CORSING Wavelet-Fourier method

We are now in a position to introduce the CORSING $\mathcal{WF}$ (Wavelet-Fourier) method.

Normalization with respect to the $H^{1}(\mathcal{D})$-norm.

Let $\Psi=\{\psi_{{\boldsymbol{j}}}\}_{{\boldsymbol{j}}\in\mathcal{J}}$ be a tensor product of periodized biorthogonal B-spline wavelets (or simply a family of periodized biorthogonal B-spline wavelets for $n=1$) defined as in Section 3.1. In particular, $\Psi\in\{\Psi^{\textnormal{ani}},\Psi^{\textnormal{iso}}\}$ and $\mathcal{J}\in\{\mathcal{J}^{\textnormal{ani}},\mathcal{J}^{\textnormal{iso}}\}$. Let $\Xi=\{\xi_{{\boldsymbol{q}}}\}_{{\boldsymbol{q}}\in\mathcal{Q}}$ be the Fourier basis in Section 3.2. Then, we normalize both trial and test functions with respect to the $H^{1}(\mathcal{D})$-norm, namely,

[TABLE]

such that

[TABLE]

In view of (3.14), (3.15), and (3.17), this normalization can be realized by defining

[TABLE]

Petrov-Galerkin discretization.

We consider a Petrov-Galerkin discretization of (2.5) associated with the trial basis $\widehat{\Psi}$ and test basis $\widehat{\Xi}$ (see [26] for an introduction to the Petrov-Galerkin method). The stiffness matrix $B\in\mathbb{C}^{M\times N}$, with $M,N\in\mathbb{N}$ and $M\geq N$, and the load vector ${\boldsymbol{g}}\in\mathbb{C}^{M}$ are defined as

[TABLE]

where

[TABLE]

and the resulting Petrov-Galerkin discretization is given by the linear system

[TABLE]

with ${\boldsymbol{v}}\in\mathbb{C}^{N}$ the vector of the unknowns.

The CORSING $\mathcal{WF}$ method.

The next step is to reduce the dimensionality of the Petrov-Galerkin discretization (3.22) by random subsampling. Given a probability distribution ${\boldsymbol{p}}\in\mathbb{R}^{M}$ over $\mathcal{Q}$, we draw $m\ll M$ multi-indices ${\boldsymbol{\tau}}_{1},\ldots,{\boldsymbol{\tau}}_{m}\in\mathcal{Q}$ i.i.d.  randomly according to

[TABLE]

Then, we define the CORSING stiffness matrix $A\in\mathbb{C}^{m\times N}$ and load vector ${\boldsymbol{f}}\in\mathbb{C}^{m}$ as

[TABLE]

The CORSING reduced discretization corresponds to the underdetermined linear system

[TABLE]

with ${\boldsymbol{z}}\in\mathbb{C}^{N}$ the vector of the unknowns. Then, we consider the diagonal preconditioner $D\in\mathbb{C}^{m\times m}$, defined as666The diagonal preconditioner $D$ is chosen in such a way that $\mathbb{E}[(DA)^{*}DA]=B^{*}B$. For further details, see [12].

[TABLE]

where $\delta_{i,k}$ is Kronecker’s delta function. Given a target sparsity level $s\in\mathbb{N}$, with $s\ll N$, we consider the following optimization problem:

[TABLE]

Although NP-hard, (3.26) can be approximated by the orthogonal matching pursuit (OMP) algorithm. Finally, the CORSING solution is given by

[TABLE]

where $\widehat{{\boldsymbol{u}}}=(\widehat{u}_{{\boldsymbol{j}}})_{{\boldsymbol{j}}\in\mathcal{J}}\in\mathbb{C}^{N}$ is an approximate solution to (3.26) computed via OMP.

Note that the sampling complexity $m$ can be reduced in order to avoid repeated indices among the $\tau_{i}$’s. In that case, the preconditioner $D$ has to be slightly modified (see [12, Remark 3.9]).

3.4 Towards a recovery error analysis

We are now able to restate issues (i), (ii), and (iii) in Section 1.1 in a rigorous way. Assuming to fix $s,N\in\mathbb{N}$ as parameters chosen by the user, with $s\ll N$, the CORSING procedure outlined above depends on the following three choices:

1. (i)

Choose $R=R(s,N)$, defining the size of $\mathcal{Q}$ and, consequently, the test space dimension $M=R^{n}$ of the Petrov-Galerkin discretization (3.22). 2. (ii)

Choose the number $m=m(s,N,M)$ of random samples, depending sublinearly on $N$ and $M$ (in order to have dimensionality reduction from (3.22) to (3.24)); 3. (iii)

Define the probability distribution ${\boldsymbol{p}}\in\mathbb{R}^{M}$ on the test multi-index space $\mathcal{Q}$.

We observe that the choice of $s$ and $N$ theoretically depends on the best $s$-term approximation error decay rate for the target function class (containing $u$). This, in turn, could be studied using results from nonlinear approximation theory (see, e.g., [19] or [24, Chapter 9]).

In order to solve issues (i), (ii), and (iii), we resort to the theoretical analysis carried out in [12], based on the following

Definition 3.3** (Local $a$-coherence).**

The local $a$-coherence of $\widehat{\Psi}$ with respect to $\widehat{\Xi}$ is a sequence ${\boldsymbol{\mu}}\in\ell(\mathbb{Z}^{n})$ defined by

[TABLE]

Since the CORSING solution $\widehat{u}$ is $s$-sparse with respect to $\Psi$ (or, equivalently, to $\widehat{\Psi}$), the corresponding best possible accuracy is the best $s$-term approximation error

[TABLE]

We specialize [12, Theorem 3.15] to the CORSING $\mathcal{WF}$ setting, providing a recovery error estimate in expectation. The fact that $\widehat{\Psi}$ and $\widehat{\Xi}$ are Riesz bases with respect to the $H^{1}(\mathcal{D})$-norm (guaranteed by Theorem 3.2) is required to apply this result, and this justifies the use of biorthogonal wavelets when $n>1$. Indeed, tensorizing the hierachical basis of hat functions as in [8, 11, 12] does not give the Riesz basis assumption in the multi-dimensional case. It is also possible to state an analogous result in probability instead of expectation (see [12, Theorems 3.16 and 3.18]).

Theorem 3.4** (CORSING $\mathcal{WF}$ recovery in expectation).**

Let $n,s,L\in\mathbb{N}$, with $s\ll N=2^{nL}$, $K>0$ be such that $\|u\|_{H^{1}(\mathcal{D})}\leq K$, where $u$ is the unique solution to (2.5), and assume to have an upper bound to the local $a$-coherence ${\boldsymbol{\mu}}\in\ell(\mathbb{Z}^{n})$, i.e., there exists a sequence ${\boldsymbol{\nu}}\in\ell(\mathbb{Z}^{n})$ such that

[TABLE]

Choose $R\in\mathbb{N}$ (or, equivalently, $\mathcal{Q}$) such that the truncation condition777Let us clarify a small difference between [12, Theorem 3.15] and Theorem 3.4. In [12, Theorem 3.15], ${\boldsymbol{\nu}}$ is an upper bound to ${\boldsymbol{\mu}}|_{\mathcal{Q}}$ and the truncation condition (3.31) involves ${\boldsymbol{\mu}}$ instead of ${\boldsymbol{\nu}}$. Therefore, the truncation condition of Theorem 3.4 implies the truncation condition of [12, Theorem 3.15]. However, this does not make any difference since in practice the truncation condition of [12, Theorem 3.15] is verified using an upper bound to ${\boldsymbol{\mu}}$, and not ${\boldsymbol{\mu}}$ itself.

[TABLE]

holds and where $\mathcal{Q}^{c}:=\mathbb{Z}^{n}\setminus\mathcal{Q}$. Then, for every $0<\varepsilon<1$, the CORSING solution $\widehat{u}\in H^{1}_{\textup{per}}(\mathcal{D})$ exactly solving (3.26) satisfies

[TABLE]

where $\alpha$ and $\mathcal{A}$ are the inf-sup and the continuity constants respectively associated with $a(\cdot,\cdot)$, provided that

[TABLE]

and that the drawings $\tau_{1},\ldots,\tau_{m}\in\mathcal{Q}$ are i.i.d. according to the probability distribution

[TABLE]

Some considerations are in order:

• •

Relation (3.33) corresponds to a quadratic scaling of $m$ with respect to $s$ (up to logarithmic factors and up to the quantity $\|{\boldsymbol{\nu}}|_{\mathcal{Q}}\|_{1}$). In practice, a linear dependence of $m$ on $s$ (up to logarithmic factors) seems to be sufficient (see the [12, Section 5.4]). A different theoretical analysis carried out in [8, Section 3.2.5] based on the concept of restricted isometry property seems to confirm this conjecture up to rescaling $DA$ by a suitable factor depending on $\mathcal{A}$ and on the true Riesz constant of $\Psi$, hidden in the norm equivalence (3.13) with $s=1$. We have preferred to employ this slightly suboptimal result to avoid this technical rescaling issue.

• •

A necessary condition for (3.31) is $\|{\boldsymbol{\nu}}\|_{1}<\infty$. This will always be the case in the applications discussed in this paper.

• •

In order to have an actual compression of the Petrov-Galerkin discretization, $\|{\boldsymbol{\nu}}|_{\mathcal{Q}}\|_{1}$ has to depend sublinearly on $N$ and $M$ in (3.33).

• •

Notice that knowing an upper bound $K$ to $\|u\|_{H^{1}(\mathcal{D})}$ is not a restrictive hypothesis in view of the a priori estimate (2.8).

• •

The hypothesis that $\widehat{{\boldsymbol{u}}}$ solves (3.26) exactly does not take into account the approximation error due to the OMP algorithm, which can be included by resorting to the restricted isometry property analysis.

• •

Condition (3.31) is sufficient to guarantee the so-called $s$-sparse restricted inf-sup property for the Petrov-Galerkin discretization (3.22). This is a variant of the classical inf-sup (or Ladyzhenskaya-Babuška-Brezzi) condition, adapted to the sparse case (see [12, Section 3.3] for more details).

Starting from Theorem 3.4, we can tackle issues (i), (ii), and (iii) in Section 1.1:

1. (i)

Choose $R=R(s,N)$ large enough to satisfy (3.31). 2. (ii)

Choose $m=m(s,N)$ according to (3.33). 3. (iii)

Choose ${\boldsymbol{p}}\in\mathbb{R}^{M}$ according to (3.34).

Our next goal is to find an upper bound ${\boldsymbol{\nu}}$ to the local $a$-coherence ${\boldsymbol{\mu}}$ as in (3.30) such that:

• (a)

It is possible to find an explicit formula for $R=R(s,N)$ such that (3.31) is satisfied.

• (b)

The quantity $\|{\boldsymbol{\nu}}|_{\mathcal{Q}}\|_{1}$ depends sublinearly on $M$ and $N$.

4 Local $a$-coherence estimates

This section is the technical core of the article. We extend the results in [12] by deriving upper bounds to the local $a$-coherence defined in (3.28) for ADR equations with nonconstant coefficients in dimension one (Section 4.1) and with constant coefficients in arbitrary dimension (Section 4.2).

4.1 The 1D ADR problem with nonconstant coefficients

We start by showing some auxiliary inequalities involving inner products between biorthogonal B-spline wavelets and the Fourier basis functions, and their respective derivatives. Note that in the following statement, the basis functions are normalized with respect to the $L^{2}(\mathcal{D})$-norm.

Lemma 4.1** (Auxiliary inequalities, $n=1$).**

In Setting 3.1, the following inequalities hold for every $\ell\geq\ell_{0}$, $k\in\mathbb{Z}/(2^{\ell}\mathbb{Z})$, $q\in\mathbb{Z}\setminus\{0\}$, $(\alpha_{1},\alpha_{2})\in\{0,1\}^{2}$, and $\gamma\in[0,2]$:

[TABLE]

Moreover, for $q=0$, we have

[TABLE]

Proof.

If $q=0$, equations (4.3)-(4.5) can be verified via direct computation, being $\xi_{0}\equiv 1$. Now, we analyze the case $q\neq 0$. Considering the case $(\alpha_{1},\alpha_{2})=(1,1)$, thanks to hypothesis (3.5), we can directly compute

[TABLE]

The inequality $\sin^{2}(x)\leq|x|^{\gamma}$, which holds for every $x\in\mathbb{R}$ and $\gamma\in[0,2]$ (see Figure 4), yields

[TABLE]

Moreover, thanks again to hypothesis (3.5), employing (4.10), the discrete Cauchy-Schwarz inequality, and the definition of the mother wavelet (see also [10])

[TABLE]

we obtain

[TABLE]

This concludes the case $(\alpha_{1},\alpha_{2})=(1,1)$. The case $(\alpha_{1},\alpha_{2})\neq(1,1)$ can be addressed using integration by parts since $\xi_{q}^{\prime}=(2\pi\text{i}q)\xi_{q}$.

∎

We are now in a position to estimate the local $a$-coherence of the wavelet basis with respect to the Fourier basis for 1D ADR equations with nonconstant coefficients.

Theorem 4.2** (Local $a$-coherence upper bound, $n=1$).**

In Setting 3.1, for $\eta,\beta\in H^{1}_{\textup{per}}(\mathcal{D})$ and $\rho\in L^{2}(\mathcal{D})$, the local $a$-coherence in (3.28) can be bounded from above as

[TABLE]

Proof.

First, let us consider the case $q=0$. We have $|(\eta\psi_{\ell,k}^{\prime},\xi_{0}^{\prime})|=|(\eta\psi_{\ell,k}^{\prime},0)|=0$,

[TABLE]

and

[TABLE]

Therefore,

[TABLE]

which, by maximization over $\ell$ and $k$, implies the estimate for $\mu_{0}$.

Now, let $q\neq 0$. The idea is to expand the diffusion, advection, and reaction terms with respect to the Fourier basis

[TABLE]

with $\eta_{r}:=(\eta,\xi_{r})$, $\beta_{r}:=(\beta,\xi_{r})$, and $\rho_{r}:=(\rho,\xi_{r})$, for every $r\in\mathbb{Z}$. The decay of the coefficients $(\eta_{r})_{r\in\mathbb{Z}}$, $(\beta_{r})_{r\in\mathbb{Z}}$, and $(\rho_{r})_{r\in\mathbb{Z}}$ is strictly linked with the Sobolev regularity of $\eta$, $\beta$, and $\rho$, respectively, thanks to the norm equivalence (2.3). The estimate of $\mu_{q}$ is divided into four parts. In the first three parts, we assess the impact of the terms $\mu$, $\beta$, and $\rho$ on the final estimate separately. Then, we combine these estimates in the fourth part.

Part I: diffusion term $\eta$ ($q\neq 0$).

We start by considering the diffusion term

[TABLE]

Recalling (2.3), the first factor is bounded from above by $\|\eta\|_{H^{1}(\mathcal{D})}^{2}$. Now, exploiting standard algebraic properties of the Fourier basis, performing integration by parts, and using (4.2) and (4.4), we obtain the following estimate for every $\gamma\in[0,2]$:

[TABLE]

where

[TABLE]

Now, we study the asymptotic behavior of $S_{y}(q)$ for $y=1$, and $y=2$ (corresponding to $\gamma=1/2$, and $\gamma=0$, respectively).888Choosing $\gamma=1/2$ and $\gamma=0$ leads to the estimates $S_{1}(q)\lesssim 1/|q|$ and $S_{2}(q)\lesssim 1/|q|^{2}$, respectively. These, in turn, imply two upper bounds to $\mu_{q}$: the first one independent of $N$ and decaying linearly with respect to $q$, the second one linear in $N$ but decaying quadratically with respect to $q$. These two properties will be crucial to answer issues (i) and (ii) in Theorem 4.4.

We start by considering $S_{1}(q)$, when $q>0$. We have the splitting

[TABLE]

We study the three sums separately, to show that each term can be bounded from above by $1/q$, up to a constant. Indeed,

[TABLE]

The first inequality employed to bound $S_{1,2}(q)$ relies on a property of the function $g(r):=1/[(q-r)(1+r^{2})]$ such that for every $q\geq 2$, $g$ admits only one stationary point $r^{*}=\frac{1}{3}(q-\sqrt{q^{2}-3})$ in the open interval $(0,q-1)$ such that $g(r^{*})\leq g(q-1)\leq g(0)$, and $g$ decreases monotonically in $[0,r^{*}]$ and increases monotonically in $[r^{*},q-1]$. In particular, this implies

[TABLE]

Also, notice that when $q=1$ the term $S_{1,2}(q)$ is equal to 1.

Observing that $S_{y}(q)$ is even with respect to $q$, we conclude that $S_{1}(q)\lesssim 1/|q|$, for every $q\neq 0$.

We carry out a similar analysis for $S_{2}(q)$. For $q>0$, we have

[TABLE]

Using arguments similar to those employed for $S_{1}(q)$, we obtain

[TABLE]

Using again that $S_{y}(q)$ is even with respect to $q$, we conclude that $S_{2}(q)\lesssim 1/q^{2}$, for every $q\neq 0$.999In view of Lemma 4.1, we conjecture that $S_{y}(q)\lesssim 1/|q|^{y}$ for every $y\in[-4,2]$ (corresponding to $\gamma\in[0,2]$ and $y=2(1-\gamma)$). Nevertheless, proving this rigorously is not straightforward, since the terms corresponding to $S_{1,2}(q)$ and $S_{2,2}(q)$ become very difficult to analyze.

Recalling relations (4.16) and (4.17), we have

[TABLE]

Part II: advection term $\beta$ ($q\neq 0$).

Analogously to the diffusion case, we have

[TABLE]

Estimating the sum in the right-hand side as before, we obtain

[TABLE]

Part III: reaction term $\rho$ ($q\neq 0$).

We deal with the nonconstant reaction term in an analogous way. Recalling Lemma 4.1 and the norm equivalence (2.3) with $k=0$, we have

[TABLE]

for every $\gamma\in[0,2]$, where we have employed (4.1)-(4.2) with $(\alpha_{1},\alpha_{2})=(0,0)$ and $\xi_{q-r}$ for $r\neq q$ and (4.5) with $\xi_{q-r}=\xi_{0}$ for $r=q$. Choosing $\gamma=1/2$, we see that

[TABLE]

since $\sum_{r\in\mathbb{Z}\setminus\{q\}}|q-r|^{-3}\lesssim 1$, for every $q\neq 0$.

Part IV: conclusion ($q\neq 0$).

Combining (4.30), (4.32), and (4.34) finally yields

[TABLE]

As a consequence, normalizing the trial and test functions with respect to the $H^{1}(\mathcal{D})$-norm and using that $\|\psi_{\ell,k}\|_{H^{1}(\mathcal{D})}\sim 2^{\ell}$ and that $\|\xi_{q}\|_{H^{1}(\mathcal{D})}\sim|q|$ (recall (3.14) and (3.17)), we obtain

[TABLE]

This completes the proof. ∎

Remark 4.3* (Sharper upper bound for $\rho\in H^{1}_{\textup{per}}(\mathcal{D})$).*

It is not difficult to show that an upper bound for $\mu_{q}$ as in Theorem 4.2 holds when $\rho\in H^{1}_{\textup{per}}(\mathcal{D})$, with the following estimate for the local $a$-coherence when $q\neq 0$:

[TABLE]

Notice that (4.37) generalizes [12, Proposition 4.4] (in particular, we refer to [12, equation (145)]) for 1D ADR equations with constant coefficients and nonperiodic boundary conditions. $\blacksquare$

Finally, we have the CORSING $\mathcal{WF}$ recovery theorem that provides an answer to the three items (i), (ii), and (iii) in Section 3.4 for the 1D case.

Theorem 4.4** (CORSING $\mathcal{WF}$ recovery).**

In Setting 3.1 and under the same hypotheses as in Theorem 3.4, let $n=1$, $\eta,\beta\in H^{1}_{\textup{per}}(\mathcal{D})$ and $\rho\in L^{2}(\mathcal{D})$. Then, provided

• (i)

$R\sim CsN$,

• (ii)

$m\gtrsim Cs(s\ln(eN/(2s))+\ln(2s/\varepsilon))(\ln N+\ln s+\ln C)$,

• (iii)

$p_{q}\propto\begin{cases}1,&q=0,\\ \min\left\{\frac{N}{q^{2}},\frac{1}{|q|}\right\},&q\neq 0,\end{cases}$**

where

[TABLE]

the CORSING $\mathcal{WF}$ method recovers the best $s$-term approximation to $u$ in expectation, in the sense of estimate (3.32).

Proof.

Employing Theorem 4.2, we choose

[TABLE]

As a consequence, we have $\|{\boldsymbol{\nu}}\|_{1}<+\infty$. Then, to ensure condition (3.31), using that $\nu_{q}\leq CN/q^{2}$ for $q\neq 0$ we estimate

[TABLE]

Therefore, to ensure (3.31) we let $R=R(s,N)\sim CsN$. Moreover, using that $\nu_{q}\leq C/|q|$ for $q\neq 0$, we see that

[TABLE]

which depends sublinearly on $M$ and $N$, as desired. ∎

4.2 The multi-dimensional case

We consider issues (i), (ii), and (iii) in Section 3.4 for the CORSING $\mathcal{WF}$ for multi-dimensional ADR equations with constant coefficients. In Section 4.2.1 we analyze the case of anisotropic tensor product wavelet, while in Section 4.2.2 we deal with the isotropic case.

4.2.1 Anisotropic tensor product wavelets

We provide local $a$-coherence upper bounds (Theorem 4.7) and a recovery result (Theorem 4.9) for the CORSING $\mathcal{WF}$ method with anisotropic tensor product wavelets.

We start by proving a technical result analogous to Lemma 4.1. To shorten notations, we introduce

[TABLE]

Moreover, we denote ${\boldsymbol{1}}=(1,1,\ldots,1)$, ${\boldsymbol{2}}=(2,2,\ldots,2)$, and so on.

Lemma 4.5** (Auxiliary inequalities, anisotropic wavelets).**

In Setting 3.1, let $n>1$, ${\boldsymbol{\ell}}\in\mathbb{N}^{n}$, with ${\boldsymbol{\ell}}\geq\ell_{0}-1$, ${\boldsymbol{k}}\in\mathbb{Z}/(2^{{\boldsymbol{\ell}}}\mathbb{Z})$, and ${\boldsymbol{q}}\in\mathbb{Z}^{n}$. Moreover, define101010Note that, according to (3.2), $\operatorname{scal}({\boldsymbol{\ell}})$ is the set of indices $j$ such that $\psi_{\ell_{j},k_{j}}$ is a scaling function.

[TABLE]

Then, it follows

• •

if ${\boldsymbol{q}}={\boldsymbol{0}}$, we have

[TABLE]

• •

if ${\boldsymbol{q}}\neq{\boldsymbol{0}}$ and $\operatorname{zero}({\boldsymbol{q}})\subseteq\operatorname{scal}({\boldsymbol{\ell}})$, then, for every $\boldsymbol{\gamma}\in[0,2]^{\|{\boldsymbol{q}}\|_{0}}$, it holds

[TABLE]

where $\widehat{{\boldsymbol{x}}}:={\boldsymbol{x}}|_{\operatorname{supp}({\boldsymbol{q}})}\in\mathbb{R}^{\|{\boldsymbol{q}}\|_{0}}$ for every ${\boldsymbol{x}}\in\mathbb{R}^{n}$ and the inequalities hide constants that depend exponentially on $n$;

• •

if ${\boldsymbol{q}}\neq{\boldsymbol{0}}$ and $\operatorname{zero}({\boldsymbol{q}})\not\subseteq\operatorname{scal}({\boldsymbol{\ell}})$, it holds

[TABLE]

Proof.

In the case ${\boldsymbol{q}}={\boldsymbol{0}}$, the equalities (4.42), (4.43), and (4.44) are a direct consequence of (4.3), (4.4), and (4.5). Let us consider ${\boldsymbol{q}}\neq{\boldsymbol{0}}$. Then, we organize the proof discussing two cases: $\|{\boldsymbol{q}}\|_{0}=n$ and $\|{\boldsymbol{q}}\|_{0}<n$.

Case $\|{\boldsymbol{q}}\|_{0}=n$.

Due to the tensorized form of the trial and the test basis functions, the following relations hold:

[TABLE]

Plugging relations (4.1) and (4.2) into (4.51) we see that, for every $\boldsymbol{\gamma}\in[0,2]^{n}$, we have

[TABLE]

Similarly, plugging (4.1) and (4.2) into (4.52), we obtain

[TABLE]

Finally, plugging (4.1) and (4.2) into (4.53), it follows

[TABLE]

The relations above prove (4.45), (4.46), and (4.47).

Case $\|{\boldsymbol{q}}\|_{0}<n$.

Let us consider the diffusion term $(\nabla\psi^{\textnormal{ani}}_{{\boldsymbol{\ell}},{\boldsymbol{k}}},\nabla\xi_{{\boldsymbol{q}}})$.

First, assume that $\operatorname{zero}({\boldsymbol{q}})\nsubseteq\operatorname{scal}({\boldsymbol{\ell}})$ (notice also that $\operatorname{zero}({\boldsymbol{q}})$ is nonempty since $\|{\boldsymbol{q}}\|_{0}<n$). In this case, we can pick an index $j_{0}\in\operatorname{zero}({\boldsymbol{q}})\setminus\operatorname{scal}({\boldsymbol{\ell}})$, i.e., such that $q_{j_{0}}=0$ and that $\ell_{j_{0}}\geq\ell_{0}$ (that is such that $\psi_{\ell_{j_{0}},k_{j_{0}}}$ is a wavelet function and not a scaling function). Combining (4.51) with relations (4.3) and (4.5) yields

[TABLE]

This proves (4.48) ((4.49) and (4.50) are shown analogously). As a consequence, the only possibility for $(\nabla\psi^{\textnormal{ani}}_{{\boldsymbol{\ell}},{\boldsymbol{k}}},\nabla\xi_{{\boldsymbol{q}}})$ to be nonzero is to have $\operatorname{zero}({\boldsymbol{q}})\subseteq\operatorname{scal}({\boldsymbol{\ell}})$. In this case, by splitting the sum above and employing (4.3), we obtain

[TABLE]

Now, splitting the product and using (4.5) yields

[TABLE]

In order to prove (4.45), it is sufficient to apply an argument analogous to the case $\|{\boldsymbol{q}}\|_{0}=n$, where the set $[n]$ is replaced with $\operatorname{supp}({\boldsymbol{q}})$, ${\boldsymbol{q}}$ with $\widehat{{\boldsymbol{q}}}$ and ${\boldsymbol{\ell}}$ with $\widehat{{\boldsymbol{\ell}}}$.

Similar arguments lead to (4.46), (4.47). ∎

Remark 4.6* (Curse of dimensionality).*

It is worth stressing that estimates (4.42)–(4.50) are affected by the curse of dimensionality, since they hide constants that blow up exponentially with $n$. This makes the CORSING $\mathcal{WF}$ approach applicable only for moderate values of $n$. This is not a major problem for fluid-dynamics applications, where the physical domain has always dimension $n\leq 3$. $\blacksquare$

Equipped with the auxiliary inequalities of Lemma 4.5, we are now in a position to provide upper bounds to the local $a$-coherence in the multi-dimensional case for anisotropic tensor product wavelets.

Theorem 4.7** (Local $a$-coherence upper bound, anisotropic wavelets).**

In Setting 3.1, let $n>1$. Then, for every $\eta,\rho\in\mathbb{R}$ and ${\boldsymbol{\beta}}\in\mathbb{R}^{n}$, the following upper bounds hold:

[TABLE]

where the inequalities hide constants depending exponentially on $n$.

Proof.

Let us assume ${\boldsymbol{q}}={\boldsymbol{0}}$. Recalling (4.42) and (4.43), we have

[TABLE]

Employing (4.44) and recalling (3.14) and (3.17), we obtain

[TABLE]

which, in turn, implies (4.61).

When ${\boldsymbol{q}}\neq{\boldsymbol{0}}$, we consider the cases $\|{\boldsymbol{q}}\|_{0}=n$ and $\|{\boldsymbol{q}}\|_{0}<n$.

Case $\|{\boldsymbol{q}}\|_{0}=n$.

Employing the auxiliary inequalities (4.45)–(4.47) (notice that in this case $\operatorname{zero}({\boldsymbol{q}})=\emptyset\subseteq\operatorname{scal}({\boldsymbol{\ell}})$) and recalling relations (3.14) and (3.17) on the $H^{1}(\mathcal{D})$-norm of $\psi^{\textnormal{ani}}_{{\boldsymbol{\ell}},{\boldsymbol{k}}}$ and of $\xi_{{\boldsymbol{q}}}$, for every $\boldsymbol{\gamma}\in[0,2]^{n}$ and ${\boldsymbol{q}}\neq{\boldsymbol{0}}$, we have

[TABLE]

For any fixed ${\boldsymbol{q}}\in\mathbb{Z}^{n}\setminus\{{\boldsymbol{0}}\}$, we make two different choices of $\boldsymbol{\gamma}=\boldsymbol{\gamma}({\boldsymbol{q}})$ that, in turn, generate two upper bounds to the local $a$-coherence.

We start by choosing $\boldsymbol{\gamma}={\boldsymbol{0}}$, in order to let the upper bound decay as fast as possible in ${\boldsymbol{q}}$. Squaring (4.63) and considering the maximum over ${\boldsymbol{j}}\in\mathcal{J}$, we obtain

[TABLE]

Finally, by noticing that

[TABLE]

we obtain the upper bound corresponding to the first argument of the minimum in (4.62), for $\|{\boldsymbol{q}}\|_{0}=n$.

Now, we find a second upper bound to $\mu_{\boldsymbol{q}}$ that is independent of $L$, but decays more slowly with respect to ${\boldsymbol{q}}$. Consider an index $j_{\infty}\in[n]$ such that $|q_{j_{\infty}}|=\|{\boldsymbol{q}}\|_{\infty}$. We define $\widetilde{\boldsymbol{\gamma}}=\widetilde{\boldsymbol{\gamma}}({\boldsymbol{q}})$ componentwise as

[TABLE]

and choose ${\boldsymbol{\gamma}}=\widetilde{{\boldsymbol{\gamma}}}$. Considering the set $\triangle:=\{{\boldsymbol{\ell}}\in\mathbb{N}^{n}:\ell_{0}\leq{\boldsymbol{\ell}}<L,\,\ell_{1}\geq\ell_{2}\geq\cdots\geq\ell_{n}\}$ and defining $S_{n}$ as the permutation group of the set $[n]$, we have

[TABLE]

where we have exploited the identity $\{{\boldsymbol{\ell}}\in\mathbb{N}^{n}:\ell_{0}\leq{\boldsymbol{\ell}}<L\}=\bigcup_{\sigma\in S_{n}}\sigma(\triangle)$ and the fact that $\|{\boldsymbol{\ell}}\|_{1}$ and $\|{\boldsymbol{\ell}}\|_{\infty}$ are invariant with respect to permutations of the components of ${\boldsymbol{\ell}}$. Combining (4.63) with (4.69) (with ${\boldsymbol{\gamma}}=\widetilde{{\boldsymbol{\gamma}}}$), and observing that

[TABLE]

we obtain

[TABLE]

where the right hand side does not depend on ${\boldsymbol{\ell}}$. This relation corresponds to the second argument of the minimum in (4.62) when $\|{\boldsymbol{q}}\|_{0}=n$. The case $\|{\boldsymbol{q}}\|_{0}=n$ is hence concluded.

Case $0<\|{\boldsymbol{q}}\|_{0}<n$.

The argument is analogous to the case $\|{\boldsymbol{q}}\|_{0}=n$. We just need to replace $[n]$ with $\operatorname{supp}({\boldsymbol{q}})$, $n$ with $\|{\boldsymbol{q}}\|_{0}$, ${\boldsymbol{q}}$ with $\widehat{{\boldsymbol{q}}}$, ${\boldsymbol{\ell}}$ with $\widehat{{\boldsymbol{\ell}}}$, and ${\boldsymbol{\gamma}}$ with $\widehat{{\boldsymbol{\gamma}}}$. Moreover, notice that $\|\widehat{{\boldsymbol{q}}}\|_{2}=\|{\boldsymbol{q}}\|_{2}$, $\|\widehat{{\boldsymbol{q}}}\|_{\infty}=\|{\boldsymbol{q}}\|_{\infty}$, and $\|\widehat{{\boldsymbol{\ell}}}\|_{\infty}\leq\|{\boldsymbol{\ell}}\|_{\infty}$ (note that the last relation is not an equality since $\widehat{{\boldsymbol{\ell}}}={\boldsymbol{\ell}}|_{\operatorname{supp}({\boldsymbol{q}})}$). This concludes the proof. ∎

Remark 4.8* (Consistency with the 1D case.).*

The upper bounds of Theorem 4.2 are compatible with those in Theorem 4.7. Indeed, when the coefficients are constant, they all belong to $H^{1}_{\textup{per}}(\mathcal{D})$ and it is immediate to verify that (4.37) coincides with (4.62) when $n=1$. Moreover, (4.61) with $n=1$ yields $\mu_{0}\lesssim|\rho|^{2}2^{-3\ell_{0}}$, which is sharper than the upper bound $\mu_{0}\lesssim|\rho|^{2}2^{-2\ell_{0}}$ implied by Theorem 4.2 when the coefficients are constant. This small discrepancy is due to the upper bound $|(\rho\psi_{\ell,k},\xi_{0})|\leq\|\rho\|_{L^{2}(\mathcal{D})}\|\psi_{\ell,k}\|_{L^{2}(\mathcal{D})}$ employed in the proof of Theorem 4.2, which is not sharp when $\rho$ is constant. $\blacksquare$

The local $a$-coherence estimates of Theorem 4.7 answer issues (i), (ii), and (iii) in Section 3.4 and lead to the following

Theorem 4.9** (CORSING $\mathcal{WF}$ recovery, anisotropic wavelets).**

In Setting 3.1 and under the same assumptions as in Theorem 3.4, let $n>1$ and $\eta,\rho\in\mathbb{R}$ and ${\boldsymbol{\beta}}\in\mathbb{R}^{n}$. Then, provided

• (i)

$R\sim CsN^{3-\frac{2}{n}}$,

• (ii)

$m\gtrsim C2^{n\ell_{0}}s(s\ln(eN/(2s))+\ln(2s/\varepsilon))(\ln N+\ln s+\ln C)^{n}$,

• (iii)

p_{{\boldsymbol{q}}}\propto\begin{cases}2^{-(2+n)\ell_{0}},&{\boldsymbol{q}}={\boldsymbol{0}},\\ 2^{-(n-\|{\boldsymbol{q}}\|_{0})\ell_{0}}\min\Big{\{}\frac{2^{(3\|{\boldsymbol{q}}\|_{0}-2)L}\|{\boldsymbol{q}}\|_{2}^{2}}{|\widehat{{\boldsymbol{q}}}|^{{\boldsymbol{4}}}},\frac{\|{\boldsymbol{q}}\|_{2}^{2}}{\|{\boldsymbol{q}}\|_{\infty}^{2}|\widehat{{\boldsymbol{q}}}|^{{\boldsymbol{1}}}}\Big{\}},&{\boldsymbol{q}}\neq{\boldsymbol{0}},\end{cases}**

where $N=2^{nL}$ and $C=|\eta|^{2}+\|{\boldsymbol{\beta}}\|_{2}^{2}+|\rho|^{2}$, the CORSING $\mathcal{WF}$ method with $\Psi=\Psi^{\textnormal{ani}}$ recovers the best $s$-term approximation to $u$ in expectation in the sense of estimate (3.32).

Proof.

Let us consider the upper bound ${\boldsymbol{\nu}}$ defined according to (4.61) and (4.62). The definition of $p_{{\boldsymbol{q}}}$ in (iii) directly follows from the definition of $\nu_{{\boldsymbol{q}}}$. Now, we derive condition (i) by estimating the tail $\|{\boldsymbol{\nu}}|_{\mathcal{Q}^{c}}\|_{1}$ and using the upper bound corresponding to the first argument of the minimum in (4.62). Letting $Q:=\lfloor R/2\rfloor$ and splitting the sum involved in the 1-norm with respect to the sparsity levels of ${\boldsymbol{q}}$, we obtain

[TABLE]

Now, we analyze $T(s)$. It is easy to verify that, for every $s\in[n]$, it holds

[TABLE]

where $X_{k}^{{\boldsymbol{t}}}\subseteq\mathbb{Z}^{s}$ is the set of multi-indices having the $k^{\text{th}}$ component larger than $Q$ and with sign pattern ${\boldsymbol{t}}$ (see Figure 5 for $s=2$ and $Q=2$).

Since the function ${\boldsymbol{r}}\mapsto\|{\boldsymbol{r}}\|_{2}^{2}/|{\boldsymbol{r}}|^{\boldsymbol{4}}$ is invariant with respect to $\operatorname{sign}({\boldsymbol{r}})$ and to permutations of the components of ${\boldsymbol{r}}$, we can just consider the set $X_{1}^{\mathbf{1}}$. Moreover, the number of possible sets $X_{k}^{\boldsymbol{t}}$ in $\mathbb{Z}^{s}$ is $2^{s}s$. Therefore, we estimate

[TABLE]

where $c:=\sum_{k\in\mathbb{N}}1/k^{4}<\infty$ and where the last inequality involves a constant depending exponentially on $s$ and hence on $n$ (since $s\leq n$). Plugging the estimate for $T(s)$ above into (4.73) yields

[TABLE]

Condition (i) is obtained from $s\|{\boldsymbol{\nu}}|_{\mathcal{Q}^{c}}\|_{1}\lesssim 1$. Of course, the above inequality also shows implicitly that $\|{\boldsymbol{\nu}}\|_{1}<+\infty$.

Finally, in order to prove (ii), we estimate $\|{\boldsymbol{\nu}}|_{\mathcal{Q}}\|_{1}$. Employing the upper bound corresponding to (4.61) and to the second argument of the minimum in (4.62), and recalling that $\|{\boldsymbol{q}}\|_{2}^{2}\leq\|{\boldsymbol{q}}\|_{0}\|{\boldsymbol{q}}\|_{\infty}^{2}$ for every ${\boldsymbol{q}}\in\mathbb{Z}^{n}$, we obtain

[TABLE]

which proves the theorem. ∎

4.2.2 Isotropic tensor product wavelets

As in the previous sections, we provide local $a$-coherence upper bounds (Theorem 4.11) and a recovery result (Theorem 4.12) for the CORSING $\mathcal{WF}$ method.

We start by proving auxiliary inequalities analogous to those in Lemma 4.5. We skip the proof, which follows the same arguments as in Lemma 4.5.

Lemma 4.10** (Auxiliary inequalities, isotropic wavelets).**

In Setting 3.1, let $n>1$, $\ell\in\mathbb{N}$, with $\ell\geq\ell_{0}$, ${\boldsymbol{k}}\in(\mathbb{Z}/(2^{\ell}\mathbb{Z}))^{n}$, ${\boldsymbol{e}}\in\{0,1\}^{n}$, and ${\boldsymbol{q}}\in\mathbb{Z}^{n}$. Moreover, define $\operatorname{zero}({\boldsymbol{q}})$ as in (4.41). Then, the following inequalities hold:

• •

If ${\boldsymbol{q}}={\boldsymbol{0}}$, we have

[TABLE]

• •

if ${\boldsymbol{q}}\neq{\boldsymbol{0}}$ and $\operatorname{zero}({\boldsymbol{q}})\subseteq\operatorname{zero}({\boldsymbol{e}})$, then, for every $\boldsymbol{\gamma}\in[0,2]^{\|{\boldsymbol{q}}\|_{0}}$, it holds

[TABLE]

where $\widehat{{\boldsymbol{x}}}:={\boldsymbol{x}}|_{\operatorname{supp}({\boldsymbol{q}})}\in\mathbb{R}^{\|{\boldsymbol{q}}\|_{0}}$ for every ${\boldsymbol{x}}\in\mathbb{R}^{n}$ and the above inequalities hide constants depending exponentially on $n$;

• •

if ${\boldsymbol{q}}\neq{\boldsymbol{0}}$ and $\operatorname{zero}({\boldsymbol{q}})\not\subseteq\operatorname{zero}({\boldsymbol{e}})$, it holds

[TABLE]

In the following theorem, we provide upper bounds to the local $a$-coherence for isotropic tensor product wavelets. We only outline a sketch of its proof, which is analogous to that of Theorem 4.7.

Theorem 4.11** (Local $a$-coherence upper bound, $n>1$, isotropic wavelets).**

In Setting 3.1, let $n>1$. Then, for every $\eta,\rho\in\mathbb{R}$ and ${\boldsymbol{\beta}}\in\mathbb{R}^{n}$, the following upper bounds hold:

[TABLE]

where the inequalities hide constants depending exponentially on $n$.

Proof.

The upper bound (4.92) to $\mu_{\boldsymbol{0}}$ is easy to verify. Let us consider ${\boldsymbol{q}}\neq{\boldsymbol{0}}$. Employing the auxiliary inequalities of Lemma 4.10 and using arguments analogous to those in Theorem 4.7, we obtain, for every ${\boldsymbol{\gamma}}\in[0,2]^{n}$,

[TABLE]

Analogously to Theorem 4.7, we obtain the upper bound corresponding to the first argument of the minimum in (4.93) by choosing $\widehat{{\boldsymbol{\gamma}}}={\boldsymbol{0}}$ in the inequality above. Notice in particular that

[TABLE]

where, in the second expression, we have used that $2^{2(-\frac{n}{2}+\|{\boldsymbol{q}}\|_{0}-1)\ell}\leq 1$ when ${-\frac{n}{2}+\|{\boldsymbol{q}}\|_{0}-1}$ is negative. The upper bound corresponding to the second argument of the minimum in (4.93) is proved by letting ${\boldsymbol{\gamma}}=\widetilde{{\boldsymbol{\gamma}}}$, as in (4.66). Note that, in this case, we have $\|\widehat{{\boldsymbol{\gamma}}}\|_{1}=\frac{3}{2}\|{\boldsymbol{q}}\|_{0}-1$ and, consequently,

[TABLE]

This concludes the proof. ∎

Finally, we obtain the CORSING $\mathcal{WF}$ recovery theorem for isotropic tensor product wavelets solving issues (i), (ii), and (iii) in Section 3.4. The proof is analogous to that of Theorem 4.9 and therefore will be omitted.

Theorem 4.12** (CORSING $\mathcal{WF}$ recovery, isotropic wavelets).**

In Setting 3.1 and under the same assumptions as in Theorem 3.4, let $n>1$ and let $\eta,\rho\in\mathbb{R}$ and ${\boldsymbol{\beta}}\in\mathbb{R}^{n}$. Then, provided

• (i)

$R\sim CsN^{3-\frac{2}{n}}$,

• (ii)

$m\gtrsim C2^{n\ell_{0}}s(s\ln(eN/(2s))+\ln(2s/\varepsilon))(\ln N+\ln s+\ln C)^{n}$,

• (iii)

p_{{\boldsymbol{q}}}\propto\begin{cases}2^{-(2+n)\ell_{0}},&{\boldsymbol{q}}={\boldsymbol{0}},\\ \min\Big{\{}\frac{(1+2^{2(-\frac{n}{2}+2\|{\boldsymbol{q}}\|_{0}-1)L})\|{\boldsymbol{q}}\|_{2}^{2}}{|\widehat{{\boldsymbol{q}}}|^{{\boldsymbol{4}}}},\frac{2^{-(n-\|{\boldsymbol{q}}\|_{0})\ell_{0}}\|{\boldsymbol{q}}\|_{2}^{2}}{\|{\boldsymbol{q}}\|_{\infty}^{2}|\widehat{{\boldsymbol{q}}}|^{{\boldsymbol{1}}}}\Big{\}},&{\boldsymbol{q}}\neq{\boldsymbol{0}},\end{cases}**

where $C=|\eta|^{2}+\|{\boldsymbol{\beta}}\|_{2}^{2}+|\rho|^{2}$, the CORSING $\mathcal{WF}$ method with $\Psi=\Psi^{\textnormal{iso}}$ recovers the best $s$-term approximation to $u$ in expectation in the sense of estimate (3.32).

5 Numerical assessment

In this section, numerically investigate the reliability and robustness of the CORSING $\mathcal{WF}$ approach in different dimensions. In Section 5.1, we consider a 1D ADR equation with constant and nonconstant coefficients. As predicted by the theory, we show that CORSING $\mathcal{WF}$ is robust and reliable also in the case of nonconstant coefficients. In Section 5.2, we consider the 2D case and compare the performance of isotropic and anisotropic wavelets in different case studies. Moreover, it turns out that the nonuniform subsampling strategy based on the local $a$-coherence significantly outperforms uniform random subsampling. Finally, in Section 5.3, we consider a 3D case.

All the numerical experiments have been performed in Matlab${}^{\text{\textregistered}}$ with the aid of OMP-Box for OMP [29, 30]. We have used Matlab${}^{\text{\textregistered}}$ R2017b version 9.3 64-bit on a MacBook Pro equipped with a 3 GHz Intel Core i7 processor and with 8 GB DDR3 RAM.

5.1 1D case with nonconstant diffusion

We consider a 1D equation with $\beta=0$, $\rho=1$ and let the diffusion coefficient vary. In particular, we consider $\eta(x_{1})\equiv 1$ and $\eta(x_{1})=1+0.5\sin(6\pi x_{1})$.

The Petrov-Galerkin stiffness matrix $B$.

We compute the stiffness matrix $B$ associated with the Petrov-Galerkin discretization and show the absolute value of the entries in Figure 6 (top row).

We set $L=9$, resulting in $N=512$, and choose $R=N$. We observe that the oscillations of the diffusion coefficient only impact the stiffness matrix “horizontally”. In particular, this does not impact the decay properties of the local $a$-coherence ${\boldsymbol{\mu}}$ (see Figure 6, bottom row). The Wavelet-Fourier Petrov-Galerkin discretization of the ADR equation gives rise to a matrix with a comparable structure with respect to the matrices in [6, Figure 4]. This qualitatively confirms that the proposed discretization is suitable for compressed sensing. We also point out that the condition number of $B\in\mathbb{C}^{512\times 512}$ is very small, being $20.4$ for the constant diffusion case and $28.2$ for the nonconstant diffusion case, to be compared with $10^{6}$ and $1.6\cdot 10^{6}$, respectively, when the trial and test functions are not normalized with respect to the $H^{1}(\mathcal{D})$-norm.

Best $s$-term approximation.

We consider the synthetic solution

[TABLE]

This solution is smooth, with a global support in $[0,1)$, and exhibits a bump close to the point $x=0.3$. Moreover, it is periodic up to machine precision (see Figure 7 (left)). In Figure 7 (right), we show the wavelet coefficients and highlight in red the largest 50 ones in absolute value.

The resulting relative best 50-term approximation error with respect to the $H^{1}(\mathcal{D})$-norm is

[TABLE]

where ${\boldsymbol{u}}_{1}$ is the vector of coefficients of $u_{1}$ with respect to the biorthogonal wavelet basis $\widehat{\Psi}$ (normalized with respect to the $H^{1}(\mathcal{D})$-norm), $\widetilde{{\boldsymbol{u}}}_{1}$ is the best 50-term approximation of ${\boldsymbol{u}}_{1}$, and $\widetilde{u}_{1}$ is the function corresponding to the wavelet coefficients in $\widetilde{{\boldsymbol{u}}}_{1}$.

Sensitivity of the recovery error to the number of test functions.

In Figure 8, we show the box plot of the relative error between $u_{1}$ and the CORSING approximation, $\widehat{u}_{1}$, with respect to the $H^{1}(\mathcal{D})$-norm as a function of the number $m$ of test functions.

We fix $s=50$ and let $m$ vary from $100$ to $500$ ($N=512$). The data are relative to 100 random runs of the CORSING procedure. We can appreciate that for both choices of $\eta$ CORSING is able to reach a good accuracy (less than twice the best $50$-term approximation error) for $m\geq 250$. The presence of a nonconstant diffusion term does not impact the performance of the method to a substantial extent. We only observe more outliers for the nonconstant diffusion.

Sensitivity of the recovery error to the sparsity.

In Figure 9, we show the relative CORSING error with respect to the $H^{1}(\mathcal{D})$-norm as a function of $s$ and compare it with the best $s$-term approximation error.

The box plot is relative to 100 runs of CORSING. For each value of $s$ varying between 5 and 50, we set $m=\lceil 2s\log(N)\rceil$. We remark that, for $s$ large enough, the recovery error exhibits the same decay rate as the best $s$-term approximation error. No striking difference can be detected varying the diffusion.

5.2 2D case

We consider a 2D ADR problem over $\mathcal{D}=(0,1)^{2}$ with constant coefficients $\mu=\rho=1$, and ${\boldsymbol{b}}=[1,1]^{T}$.

On the one hand, we compare the performance of anisotropic and isotropic wavelets on solutions that exhibit different spatial features. On the other hand, we show that nonuniform sampling strategy based on local $a$-coherence outperforms the uniform random subsampling.

Wavelet coefficients and best $s$-term approximation error.

We consider the following solutions:

[TABLE]

both periodic up to machine precision. The function $u_{2}$ exhibits two local Gaussian-shaped features, one isotropic around the point $(0.3,0.4)$ and the other anisotropic, around $(0.6,0.5)$. The function $u_{3}$ is purely anisotropic, having a Gaussian behavior along the $x_{1}$-direction and being constant along the $x_{2}$-direction. The functions $u_{2}$ and $u_{3}$ are shown in Figure 10 along with the corresponding anisotropic and isotropic $64\times 64$ wavelet coefficients (with respect to $H^{1}(\mathcal{D})$-normalized wavelets).

Letting $s=100$, the relative best $s$-term approximation error with respect to the $H^{1}(\mathcal{D})$-norm of is $10^{-1}$ (anisotropic wavelets) and $6.6\cdot 10^{-2}$ (isotropic wavelets) for $u_{2}$, and $8.2\cdot 10^{-3}$ (anisotropic wavelets) and $1.2\cdot 10^{-1}$ (isotropic wavelets) for $u_{3}$. As expected, anisotropic wavelets generate a very sparse representation of $u_{3}$. For $u_{2}$, the compression achieved by anisotropic and isotropic wavelets is comparable, in slight favor of isotropic wavelets.111111Taking advantage of the norm equivalence property, the $H^{1}(\mathcal{D})$-norm is approximated using the $\ell^{2}$-norm of the wavelet coefficients as in (5.2).

Sensitivity of the recovery error to the number of test functions.

We assess the performance of CORSING in the case of anisotropic and isotropic wavelets and compare uniform random subsampling (${\boldsymbol{p}}\propto{\boldsymbol{1}}$) with the nonuniform subsampling based on the local $a$-coherence upper bound

[TABLE]

which is obtained from the upper bounds in Theorems 4.7 and 4.11. In particular, in (4.62) and (4.93), we consider the second argument of the minimum and use that $2^{-(n-\|{\boldsymbol{q}}\|_{0})\ell_{0}}\leq 1$, while we use that $2^{-(2+n)\ell_{0}}\leq 1$ in (4.61) and (4.92). We set $\ell_{0}=2$, $L=6$ (corresponding to $N=2^{2L}=4096$). As for the test space, we fix $R=2^{L}$, corresponding to $M=N$. Although issue (i) in Theorems 4.9 and 4.12 suggests choosing $R\sim sN^{3-\frac{2}{n}}$, the choice $R=N$ turns out to be sufficient to have a well-conditioned Petrov-Galerkin discretization matrix $B$ in practice. We set $s=100$ and let $m=100,200,300,400,500$. For each value of $m$, we run $100$ tests of CORSING with uniform and nonuniform subsampling. We plot the relative recovery error measured with respect to the $H^{1}(\mathcal{D})$-norm as a function of the number of tests $m$ in Figures 11 and 12 for $u_{2}$ and $u_{3}$, respectively.

In the case of $u_{2}$, isotropic wavelets slightly outperform anisotropic wavelets. The benefit of nonuniform subsampling over uniform subsampling is evident both in terms of the probability of success (smaller boxes) and of accuracy. For $u_{3}$, anisotropic wavelets significantly outperform isotropic wavelets, thanks to the better compressibility of the solution. Moreover, uniform sampling fails to recover the solutions in both cases, whereas nonuniform sampling exhibits a convergent behavior. This experiment confirms the key role played by the local $a$-coherence for a successful implementation of the CORSING $\mathcal{WF}$ method.

5.3 The 3D case

We validate the CORSING $\mathcal{WF}$ method on a 3D ADR problem on $\mathcal{D}=(0,1)^{3}$ with constant coefficients $\mu=\rho=1$, ${\boldsymbol{b}}=[1,1,1]^{T}$. We consider the exact solution

[TABLE]

The function $u_{4}$ exhibits an anisotropic Gaussian-shaped feature centered at the point $(0.4,0.5,0.6)$. We compare anisotropic and isotropic wavelets.

Wavelet coefficients and best $s$-term approximation.

We fix $\ell_{0}=2$, $L=4$ (corresponding to a trial space of dimension $N=2^{3L}=4096$) and $s=200$. The wavelet coefficients and the best $s$-term approximation are shown in Figure 13.

The relative best $s$-term approximation error with respect to the $H^{1}(\mathcal{D})$-norm is $9.3\cdot 10^{-2}$ for the anisotropic wavelets and $2.5\cdot 10^{-2}$ for the isotropic wavelets. In this case, the isotropic tensorization is able to sparsify the function to a slightly better extent.

Sensitivity of the recovery error to the number of test functions.

We compare anisotropic and isotropic wavelets and uniform and nonuniform subsampling. We set $m=200,300,400,500,600$. The box plots corresponding to $100$ runs of CORSING are shown in Figure 14.

The performance of anisotropic and isotropic wavelets is similar. However, uniform subsampling is not able to recover the solution at all. Comparing this results with those of the 2D case (Figures 11 and 12), we note how using a bad probability measure (i.e., the uniform) deteriorates the performance of the method more heavily as the dimension of the domain increases. Finally, this experiment confirms that the theoretical analysis carried out in Section 4 turns out to be a useful tool for an effective implementation of the CORSING $\mathcal{WF}$ method.

6 Conclusions

We presented a wavelet-Fourier discretization technique for ADR equations based on the Petrov-Galerkin method and on the compressed sensing paradigm, called CORSING $\mathcal{WF}$. We carried out a theoretical analysis of the method, which hinges on the concept of local $a$-coherence and provides practical recipes for a successful implementation of the method. Numerical experiments confirm the robustness and reliability of the CORSING $\mathcal{WF}$ approach for $n$-dimensional ADR equations with $n=1,2,3$.

In particular, we showed that the method achieves a recovery error comparable to the best $s$-term approximation error, and that the sampling measure based on the local $a$-coherence proposed here is able to successfully exploit the sparsity of the exact solution in the discretization (in contrast to other randomization strategies such as uniform random subsampling).

Several open issues still remain to be investigated. First, one needs to understand whether the sampling measure proposed in this paper is or is not the “optimal” one (in some sense to be specified). On the practical and computational side there is still a lot of work to be done. Although we compared the accuracy of the CORSING $\mathcal{WF}$ solution with the best $s$-term approximation error (Figure 9), the computational cost of the CORSING $\mathcal{WF}$ procedure with OMP reconstruction scales linearly in $N$ (i.e., the dimension of the trial basis of wavelet functions). Yet, adaptive wavelet methods can recover the best $s$-term approximation error accuracy with an optimal computational cost of $O(s)$ flops. This is a crucial issue to address in order to understand what is the real impact of the CORSING $\mathcal{WF}$. In this direction, a line of research currently under investigation is the use of techniques for sublinear-time compressed sensing recently proposed in [16]. Finally, developing an effective and optimized implementation for CORSING $\mathcal{WF}$ that takes advantage of the wavelet transform, the Fourier transform, and of the tensor product structure of the basis functions is still an open issue that has to be tackled to implement CORSING $\mathcal{WF}$ in dimension $n>3$.

Acknowledgements

The first author acknowledges the support of the Postdoctoral Training Centre in Stochastics of the Pacifical Institute for the Mathematical Sciences (PIMS), the Centre for Advanced Modelling Science (CADMOS), and the Natural Sciences and Engineering Research Council of Canada through grant number 611675 for the financial support. The first author thanks Ben Adcock, Wolfgang Dahmen, and Holger Rauhut for very insightful discussions about approximation theory and compressed sensing. The fourth author acknowledges the research project GNCS-INdAM 2018 “Tecniche di Riduzione di Modello per le Applicazioni Mediche”, which partially supported this research. The authors would also like to thank the two anonymous reviewers for their helpful and constructive comments.