Linear reconstructions and the analysis of the stable sampling rate

TL;DR

This paper analyzes the stable sampling rate (SSR) in linear reconstruction methods, showing its linearity with Walsh functions and wavelets, and providing exact SSR calculations for certain wavelets to improve reconstruction accuracy.

Contribution

It demonstrates the linearity of SSR in binary Walsh sampling and wavelet reconstruction, and derives exact SSR values for specific wavelets, enhancing understanding of sampling efficiency.

Findings

SSR is linear with Walsh functions and wavelet reconstruction.

Exact SSR can be determined for certain wavelets.

Results enable sharper performance estimates for reconstruction methods.

Abstract

The theory of sampling and the reconstruction of data has a wide range of applications and a rich collection of techniques. For many methods a core problem is the estimation of the number of samples needed in order to secure a stable and accurate reconstruction. This can often be controlled by the Stable Sampling Rate (SSR). In this paper we discuss the SSR and how it is crucial for two key linear methods in sampling theory: generalized sampling and the recently developed Parametrized Background Data Weak (PBDW) method. Both of these approaches rely on estimates of the SSR in order to be accurate. In many areas of signal and image processing binary samples are crucial and such samples, which can be modelled by Walsh functions, are the core of our analysis. As we show, the SSR is linear when considering binary sampling with Walsh functions and wavelet reconstruction. Moreover, for certain wavelets it is possible to determine the SSR exactly, allowing sharp estimates for the performance of the methods.

Figures

Equations

$$\cos (\omega ({\mathcal{R}}_N,{\mathcal{S}}_M)):=\underset{r\in {\mathcal{R}}_N,\Vert r\Vert =1}{\inf }\Vert P_{{\mathcal{S}}_M}r\Vert .$$

$$\mu ({\mathcal{R}}_N,{\mathcal{S}}_M)=1/\cos (\omega ({\mathcal{R}}_N,{\mathcal{S}}_M))\in [1,\infty ],$$

$$\Theta (N,\theta )=\min \left\{M\in \mathbb{N}:\mu ({\mathcal{R}}_N,{\mathcal{S}}_M)<\theta \right\}.$$

$$\mathrm{Wal}(n,x)=(-1{)}^{{\sum }_{i\in \mathbb{Z}}(n_i+n_{i+1})x_{-i-1}}.$$

$$\mathrm{Wal}(n,x)=\underset{k=1}{\overset{d}{⨂}}\mathrm{Wal}(n_k,x_k).$$

$${\mathcal{S}}_M=\mathrm{span}\left\{\mathrm{Wal}(n,\cdot ),n=(n_k{)}_{k=1,\dots ,d},n_k=1,\dots ,m,k=1,\dots ,d\right\}.$$

$$\text{stackon}[2pt]f\text{stretchto}\text{scalerel}\ast [\text{widthof}f]\bigwedge {0.5ex}_{\text{belowbaseline}[-2pt]W}(n)=⟨f(\cdot ),\mathrm{Wal}(n,\cdot )⟩={\int }_{[0,1{]}^d}f(x)\mathrm{Wal}(n,x)dx,n\in {\mathbb{R}}^d.$$

$${\psi }_{R,j}(x):=2^{R/2}\psi (2^{R}x-j)\text{ and }{\varphi }_{R,j}(x):=2^{R/2}\varphi (2^{R}x-j),$$

$$V_r^b=\mathrm{span}\left\{{\varphi }_{r,j}:j=0,\dots ,2^r-p-1,{\varphi }_{r,j}^{\#}:j=2^r-p,\dots ,2^r-1\right\},$$

$$\bigcup _{r<R}{W}_r^b=V_R^b.$$

$${\mathcal{R}}_N:=V_R^b.$$

$${\mathcal{R}}_N=V_R^{b,d}:=V_R^b\otimes \dots \otimes V_R^b\text{ (d-times)}$$

$$\mathcal{F}(f,\zeta )=0,$$

$${\mathcal{R}}_0\subset {\mathcal{R}}_1\subset \dots \subset {\mathcal{R}}_N,\dim ({\mathcal{R}}_j)=j,$$

$$\mathcal{K}=\left\{f\in \mathcal{H}:\mathrm{dist}(f,{\mathcal{R}}_N)\le {ϵ}_N\right\}\text{ and }{\mathcal{H}}_s=\left\{f\in \mathcal{H}:P_{{\mathcal{S}}_M}f=s\right\}.$$

$$g^{\ast }={\mathrm{argmin}}_{g\in {\mathcal{R}}_N}||s-P_{{\mathcal{S}}_M}g|{|}^2,$$

$$A^{\ast }(P_{{\mathcal{S}}_M}f):=f^{\ast }=s+P_{{\mathcal{S}}_M^{\perp }}{g}^{\ast }.$$

$$||f-A(s)||\le C_A(s)\mathrm{dist}(f,{\mathcal{R}}_N),f\in {\mathcal{K}}_s.$$

$$E_A(W):=\underset{f\in W}{\sup }||f-A(P_{{\mathcal{S}}_M}f)||.$$

$$E(W):=\underset{A}{\inf }{E}_A(W).$$

$$||f-A^{\ast }(P_{{\mathcal{S}}_M}f)||\le \mu ({\mathcal{R}}_N,{\mathcal{S}}_M)\mathrm{dist}(f,{\mathcal{R}}_N).$$

$$||f-A^{\ast }(P_{{\mathcal{S}}_M}f)||\le \mu ({\mathcal{R}}_N,{\mathcal{S}}_M)\mathrm{dist}(f,{\mathcal{R}}_N\oplus ({\mathcal{S}}_M\cap {\mathcal{R}}_N^{\perp })).$$

$$⟨P_{{\mathcal{S}}_M}{G}_{N,M}(f),r_i⟩=⟨P_{{\mathcal{S}}_M}f,r_i⟩,i=1,\dots ,N.$$

$$U^{[N,M]}{\alpha }^{[N]}$$

$$U^{[N,M]}$$

$$||G_{N,M}(f)||\le \mu ({\mathcal{R}}_N,{\mathcal{S}}_M)||f||,$$

$$||f-P_{{\mathcal{R}}_N}f||\le ||f-G_{N,M}(f)||\le \mu ({\mathcal{R}}_N,{\mathcal{S}}_M)||f-P_{{\mathcal{R}}_N}f||.$$

$$\kappa ({\mathcal{R}}_N,{\mathcal{S}}_M)=\mu ({\mathcal{R}}_N,{\mathcal{S}}_M).$$

$$\underset{\alpha \in {\ell }_1}{\inf }||\alpha |{|}_{{\ell }_1}\text{ with }{P}_{[M]}U\alpha =P_{[M]}U\beta .$$

$$E{G}_M(f)=\sum _{i=1}^{\infty }{\alpha }_i{r}_i.$$

Full text

Key words and phrases : Sampling, Wavelets, binary measurements, Generalized sampling, linear reconstruction,

2010 AMS Mathematics Subject Classification : 94A20, 42C10, 42C40 (primary); 65R32, 94A08, 94A12 (secondary)

1. Introduction

Sampling theory is a core principle in image and signal processing as well as the mathematics of information, data science and inverse problems. Since the early results of Shannon [25, 38, 40] many techniques have been developed, and there are now a myriad of methods available. Moreover, the many applications, such as Magnetic Resonance Imaging (MRI) [30, 21], electron tomography [27, 28], lensless cameras, fluorescence microscopy [39, 36], X-ray computed tomography [14, 35], surface scattering [26] as well as parametrised PDEs [7, 10, 17], make the field well connected to different areas of the sciences.

A standard sampling model is as follows. We have an element $f\in\mathcal{H}$, where $\mathcal{H}$ is a separable Hilbert space, and the goal is to reconstruct an approximation to $f$ from a finite number of linear samples of the form $l_{i}(f)$, $i\in\mathbb{N}$. In particular, given that the $l_{i}$s are linear functionals, we measure the scalar product between $f$ and some sampling element $s_{i}\in\mathcal{H}$, $i\in\mathbb{N}$, i.e. $l_{i}(f)=\langle f,s_{i}\rangle$. It is important to note that the $l_{i}$s cannot be chosen freely, but are dictated by the modality of the sampling device, say an MRI scanner providing Fourier samples or a fluorescence microscope giving binary measurements modelled by Walsh coefficients. A natural question that arises from this setting is, what is the number of samples that is needed for a accurate and stable reconstruction? This can be made explicit by the stable sampling rate (SSR).

To define the SSR we first introduce the sampling space and the reconstruction space. We define the sampling space $\mathcal{S}=\overline{\operatorname{span}}\{s_{i}:i\in\mathbb{N}\}\subset\mathcal{H}$, meaning the closure of the span. In practice, one can only acquire a finite number of samples, therefore, we denote by $\mathcal{S}_{M}=\operatorname{span}\{s_{i}:i=1,\ldots,M\}$ the sampling space of the first $M$ elements. Similarly, the reconstruction space denoted by $\mathcal{R}$ is spanned by reconstruction functions $(r_{i})_{i\in\mathbb{N}}$, i.e. $\mathcal{R}=\overline{\operatorname{span}}\{r_{i}:i\in\mathbb{N}\}$. As in the sampling case, one has to restrict to a finite reconstruction space, which is denoted by $\mathcal{R}_{N}=\operatorname{span}\{r_{i}:i=1,\ldots,N\}$.

The key ingredient in the definition of the SSR is the subspace angle $\omega$ between the subspaces $\mathcal{R}_{N}$ and $\mathcal{S}_{M}$. In particular,

[TABLE]

The orthogonal projection onto the sampling space is denoted by $P_{\mathcal{S}_{M}}$. Mainly, one is interested in the reciprocal value

[TABLE]

which, as we will see below, plays a key role in all the error estimates of the two linear algorithms discussed here. Due to the definition of cosine, $\mu$ takes values in $[1,\infty]$. We can now define the stable sampling rate

[TABLE]

In particular, the SSR determines how many samples $M$ are needed when given $N$ reconstruction vectors in order to bound $\mu$ by $\theta$.

Although the SSR is developed mainly for linear methods there is a strong connection to non-linear techniques. For example in infinite-dimensional compressed sensing a key condition is the so-called balancing property [5], which is very similar to the SSR. When the balancing property is satisfied, the truncated matrix $P_{[N]}UP_{[M]}$ is close to an isometry. Hence, it also ensures the stability of the reconstruction.

The stable sampling rate has already been analysed for different settings. A prominent one is the reconstruction from Fourier measurements. In this case the sampling functions $s_{i}$ are the complex exponentials. For the reconstruction of wavelets it was shown in [3] that the SSR is linear. The other main measurements are binary. Combined with wavelets they also have a linear SSR [23].

The purpose of this paper is twofold:

• •

First, we want to compare two linear reconstruction methods: generalized sampling [4] and the PBDW approach based on data assimilation [10]. We also show how both linear methods completely rely on the SSR in order to be accurate. Additionally, the non-linear extension of generalized sampling, which is presented in [34], is included in the comparison.

• •

Second, we provide sharp results on the SSR when considering Walsh functions and Haar wavelets. This can be done by realising the common structure of the Walsh functions and the Haar wavelets. Although the SSR is linear when considering Walsh samples and wavelet reconstructions, sharp results on the SSR for arbitrary Daubechies wavelets are still open. The difficulty is that the higher order Daubechies wavelets share very little structural similarities with the Walsh functions.

2. Reconstruction Methods

In terms of reconstruction methods, there are three different properties that are often desired. The most important are obviously accuracy and stability. However, consistency, meaning that the reconstruction will yield the same samples as the true solution, is also often considered an advantage. Below, we will see how the SSR is crucial for the two former properties.

2.1. Reconstruction and Sampling Space

Throughout the paper $\mathcal{H}=L^{2}([0,1]^{d})$. Due to the fact that we are dealing with the $d$-dimensional case, we introduce multi-indices to make the notation more readable. Let $j=\left\{j_{1},\ldots,j_{d}\right\}\in\mathbb{N}^{d}$, $d\in\mathbb{N}$, be a multi-index. A natural number $n$ is in the context with a multi-index interpreted as a multi-index with the same entry, i.e. $n=\left\{n,\ldots,n\right\}$. Then, we define the addition of two multi-indices for $j,r\in\mathbb{N}^{d}$ by the point-wise addition, i.e. $j+r=\left\{j_{1}+r_{1},\ldots,j_{d}+r_{d}\right\}$ and the sum $\sum_{j=k}^{r}x_{j}:=\sum_{j_{1}=k_{1}}^{r_{1}}\ldots\sum_{j_{d}=k_{d}}^{r_{d}}x_{j_{1},\ldots,j_{d}}$, where $k,r\in\mathbb{N}^{d}$. The multiplication of a multi-index with a real number is understood point-wise as well as the division by a multi-index.

We now discuss the sampling space $\mathcal{S}_{M}$. As pointed out before, we are interested in binary measurements. Binary measurements come from applications such as fluorescence microscopy or single pixel cameras. Mathematically they are modelled by scalarproducts of our function of interest $f$ and sampling functions $s_{i}$, $i\in\mathbb{N}$, which take only the values $\left\{0,1\right\}$ or $\left\{-1,1\right\}$. Without loss of generality, we assume in the further scope of this paper that the functions take values $\left\{-1,1\right\}$. The former can be attained by taking an additional measurement with the constant function. For practical reasons it is desirable to have a fast transform for these measurements. Therefore, measurement matrices which consist of random matrices are not desirable. Additionally, it has been shown that the reconstruction quality is better with more structured sampling spaces. Therefore, we use Walsh functions to represent the measurements. Walsh functions obey the advantage that they correspond to a fast transform. Moreover, they are the kernels of the Hadamard transform and the analogue or dyadic version of exponential functions. Hence, they are a natural choice to model binary measurements. Now, we define the Walsh functions. The Walsh functions in higher dimensions can be represented by the tensor product of one-dimensional Walsh functions.

Definition 1** (Walsh function [19]).**

Let $n=\sum_{i\in\mathbb{Z}}n_{i}2^{i-1}$ with $n_{i}\in\left\{0,1\right\}$ be the dyadic expansion of $s\in\mathbb{R}_{+}$. Analogously, let $x=\sum_{i\in\mathbb{Z}}x_{i}2^{i-1}$ with $x_{i}\in\left\{0,1\right\}$. The generalized Walsh functions in $L^{2}([0,1])$ are given by

[TABLE]

We extend it to functions in $L^{2}([0,1]^{d})$ by the tensor product for $n=(n_{k})_{k=1,\ldots,d},x=(x_{k})_{k=1,\ldots,d}$

[TABLE]

Notice that the parameter $n$ represents the number of zero crossings of the function. For this reason, it is often referred to as sequency, which is similar to frequency in the case of exponential functions. The Walsh functions span the sampling space, i.e. for $M=m^{d},m\in\mathbb{N}$ we have

[TABLE]

Moreover, Walsh functions can be extended to negative inputs by $\operatorname{Wal}(-n,x)=\operatorname{Wal}(n,-x)=-\operatorname{Wal}(n,x)$.

With the help of the Walsh functions we can define the continuous Walsh transform of a function $f\in L^{2}([0,1]^{d})$ as in [19] almost everywhere by

[TABLE]

The Walsh functions have the following useful properties. First, they obey the scaling property, i.e. $\operatorname{Wal}(2^{j}n,x)=\operatorname{Wal}(n,2^{j}x)$ for all $j\in\mathbb{N}$ and $n,x\in\mathbb{R}$. Second, the multiplicative identity holds, this means $\operatorname{Wal}(n,x)\operatorname{Wal}(n,y)=\operatorname{Wal}(n,x\oplus y)$, where $\oplus$ is the dyadic addition, i.e. the element-wise addition modulo two. These properties are also easily transferred to the Walsh transform. For further information on Walsh functions and transforms see [8, 37, 13].

Direct inversion from a finite number of samples, both in the Fourier and Walsh case, may lead to substantial artefacts such as the Gibbs phenomenon in the Fourier case or block artefacts known from Walsh functions. This can be seen in the numerical experiments in Figures 3c, 4c and 5b. Therefore, it is important to consider reconstruction spaces $\mathcal{R}$ which represent the data in a better way so that a finite and low-dimensional subspace leads to a good reconstruction. In many different applications, such as image and signal processing and also representations of solution manifolds for PDEs [32], wavelets have become highly popular alternatives. The reconstruction space is spanned by reconstruction functions $r_{i}$, $i\in\mathbb{N}$. As it is not possible to deal numerically with an infinite number of samples, it also is not possible to reconstruct infinite number of coefficients. Therefore, we examine the reconstruction space $\mathcal{R}_{N}=\operatorname{span}\left\{r_{i}:i=1,\ldots,N\right\}$. When $\mathcal{R}_{N}$ is used as an approximation for the solution manifold of a PDE, we are also given the approximation error $\epsilon_{N}$.

In the following, we use wavelets as the reconstruction space due to their good time and frequency localisation. First, we study the one-dimensional case. Then, we get to higher dimensions. We use the common notation and denote the mother wavelet with $\psi$ and the corresponding scaling function with $\phi$. These functions are then scaled and translated. This results in the functions

[TABLE]

where $R,j\in\mathbb{Z}$. The wavelet space at a certain level $r$ is given by $W_{r}:=\operatorname{span}\left\{\psi_{r,j}:j\in\mathbb{Z}\right\}$ and the scaling space is given by $V_{r}:=\operatorname{span}\left\{\phi_{r,j}:j\in\mathbb{Z}\right\}$. Often one is interested in the representation of functions in $L^{2}([0,1])$ instead of $L^{2}(\mathbb{R})$. For this sake boundary corrected wavelets were introduced in [15]. We are focussing on the version presented in chapter $4$ of [15], because they still obey the same smoothness and vanishing moments properties as the wavelets they are derived from on the real line. Additionally, they also remain the multi-resolution analysis property. We now shortly repeat the construction. The scaling space for boundary corrected wavelets is spanned by the original scaling function $\phi$ and reflections around $1$ of the scaling function $\phi^{\#}$, i.e.

[TABLE]

for Daubechies wavelets of order $p$. The boundary corrected Daubechies wavelets still obey the multi-resolution analysis. Therefore, it is possible to represent the union of the wavelet spaces up to a certain level $R-1$ by the scaling space at this level, i.e.

[TABLE]

Hence, it is not necessary for the analysis to have a deeper look in the ordering of the wavelets and their construction for the boundary corrected version as long as we have that the number of coefficients equals the number of elements in that level, i.e. if $N=2^{R}$ for some $R\in\mathbb{N}$. The reconstruction space is then given by

[TABLE]

The higher-dimensional scaling spaces are constructed by the tensor product of the one-dimensional one, such that we get in $d$ dimensions

[TABLE]

for $N=2^{dR}$. Note that the corresponding wavelet space is not simply the tensor product of the one-dimensional wavelets. Fortunately, this is not often a problem in the analysis, since Daubechies wavelets obey the multi-resolution analysis. We will deal with the internal ordering for the one- and two-dimensional case for the Haar wavelets.

2.2. Reconstruction Techniques

In the following, we present two different interesting reconstruction methods, which are both optimal in their setting. We highlight some advantages and disadvantages. Particularly, we see that the performance of both methods depends highly on the subspace angle between the sampling and the reconstruction space. This gives rise to the discussion in §3.1 about the question for which sampling and reconstruction spaces the stable sampling rate is linear.

2.2.1. PBDW-method

In [7, 10, 17] the PBDW-method from [31] is analysed. In contrast to the prominent application for image and signal analysis this method arises from the application with PDEs. One tries to estimate a state $f$ of a physical system by solving a parametric family of PDEs

[TABLE]

where $\mathcal{F}$ is a differential operator. The parameter $\zeta$ may not be known exactly. Therefore, the value of $f$ cannot be attained by simply solving the PDE. Hence, other information is necessary. In most applications, one has access to linear measurements $l_{i}(f)$, $i=1,\ldots,M$ of the state $f$. This alone is not sufficient to estimate $f$ or more ambitiously even the parameter $\zeta$. Fortunately, one also has information about the PDE, which can be used to analyse the solution manifold $\mathcal{M}$. The solution manifold is usually quite complicated and not given directly. Hence, approximations are used. The method used in this context is the approximation by a sequence of nested finite subspaces

[TABLE]

where the approximation error is known to be $\epsilon_{j}$ for each subspace $\mathcal{R}_{j}$. There are a lot of different methods that allow us to construct the spaces $\mathcal{R}_{j}$, such as the reduced basis method [9, 12, 16, 41] and the use of wavelets [32].

The given setting leads to the goal of trying to merge the data driven and model based information, which leads to the concept of PBDW-method. Let the measurement $s:=P_{\mathcal{S}_{M}}f$ be given. The idea is to combine the information given by the measurements and the PDE. This means we are searching for an approximation $f^{*}\in\mathcal{K}_{s}$ where

[TABLE]

The intersection is then the space of possible solutions, i.e. $\mathcal{K}\cap\mathcal{H}_{s}=\mathcal{K}_{s}$. The aim is to reduce the distance between the approximation $f^{*}$ and the true solution $f$.

It was shown in [10] that the following approach introduced in [31] is optimal for this task. First, the minimization problem

[TABLE]

is solved. The solution to the reconstruction problem is then given by the mapping $A^{*}:\mathcal{S}_{M}\rightarrow\mathcal{H}$ defined by

[TABLE]

For the performance analysis, we want to compare the reconstruction algorithm $A^{*}$ with all other mappings $A:\mathcal{S}_{M}\rightarrow\mathcal{H}$. In the following, let $A$ represent any alternative reconstruction algorithm which also only depends on the information of $f$ in the sampling space $\mathcal{S}_{M}$.

We first examine the instance optimality. This means we analyse the error for a given measurement $s$, i.e.

[TABLE]

The algorithm $A$ which leads to the smallest constant $C_{A}(s)$ is called instance optimal. It is clear that the error scales with the distance of the element $f$ from the reconstruction space $\mathcal{R}_{N}$. Due to the fact that we normally do not know $s$ a priori, this estimate is not very helpful. Hence, one is interested in the performance for any kind of input $s\in\mathcal{S}_{M}$. Therefore, the performance of a recovery algorithm $A$ on any subset $W\subset\mathcal{H}$ is given by

[TABLE]

Taking the infimum over the whole class gives the class optimal performance on a given set $W$ defined by

[TABLE]

It is shown in [10] that the presented algorithm is both instance and class optimal and gives the estimate

[TABLE]

Hence, with this approach we do not get reasonable estimates, if $\mathcal{R}_{N}\cap\mathcal{S}_{M}^{\perp}\neq\left\{0\right\}$. Moreover, a detailed knowledge about the stable sampling rate is necessary to get a method which can be used in practice. A reconstruction method is of little help if we do not have good error bounds or if the condition number is too high. We will see in the next chapter that the condition number of this reconstruction method is also bounded if the stable sampling rate has been taken into account, see Equation (27). Next, it was shown in [31] that this estimate can even be improved to

[TABLE]

Note that (21) demonstrates how the PBDW-method is dependent on the SSR. Moreover, it was shown in [10] that the constant $\mu(\mathcal{R}_{N},\mathcal{S}_{M})$ cannot be improved. Thus, the estimate is sharp, which further demonstrates the importance of the SSR.

2.2.2. Generalized Sampling

After the examination of the concept of the PBDW-method, we want to study the different reconstruction technique generalized sampling. Here, the approach is a stable improvement of concepts as finite section methods [11, 20, 22, 29]. The main difference is that generalized sampling allows for different dimensions on $\mathcal{S}_{M}$ and $\mathcal{R}_{N}$, whereas in the finite section method they are always the same. However, the finite section method becomes a special case of generalized sampling when the dimensions of the sampling space and reconstruction space are equal. The method is defined now. Afterwards, we discuss how the equation can be attained by solving the equivalent least square problem.

Definition 2** ([1]).**

For $f\in\mathcal{H}$ and $N,M\in\mathbb{N}$ we define the reconstruction method of generalized sampling $G_{N,M}:\mathcal{H}\rightarrow\mathcal{R}_{N}$ by

[TABLE]

We also refer to $G_{N,M}(f)$ as the generalized sampling reconstruction of $f$.

We stress at this point that generalized sampling is also a linear reconstruction scheme. In particular, Equation (22) is equivalent to solving the following linear equation for $\alpha^{[N]}\in\mathbb{R}^{N}$.

[TABLE]

where

[TABLE]

and $u_{ij}=\langle r_{j},s_{i}\rangle,~{}l(f)^{[M]}=(l_{1}(f),\ldots,l_{M}(f))\in\mathbb{R}^{M}$. The matrix can be seen in Figure 1 for different sampling and reconstruction spaces. The reconstruction is given by $G_{N,M}(f)=\sum_{i=1}^{N}\alpha_{i}r_{i}$. An interesting fact about this is that the matrix $U^{[N,N]}$ is very ill-conditioned for most cases. As pointed out in [3], in the case of Fourier samples and wavelet reconstructions, the condition number grows exponentially in $N$. This means that this approach is only feasible due to the above mentioned allowance of a different number of samples and reconstructed coefficients. It can be shown that there exists a certain number of samples such that (22) obeys a solution.

Theorem 1** ([2]).**

Let $N\in\mathbb{N}$. Then, there exists $M_{0}\in\mathbb{N}$ such that for every $f\in\mathcal{H}$ Equation (22) has a unique solution $G_{N,M}(f)$ for all $M\geq M_{0}$. Moreover, the smallest $M_{0}$ is the least number such that $\cos(\omega(\mathcal{R}_{N},\mathcal{S}_{M_{0}}))>0$.

Just as for the PBDW-method, the performance bounds of generalized sampling were studied. Here we observe again that the reconstruction quality highly depends on the subspace angle and on the relation between the data and the reconstruction space.

Theorem 2** ([2]).**

Retaining the definitions and notations from this chapter, for all $f\in\mathcal{H}$ we have

[TABLE]

and

[TABLE]

In particular, these bounds are sharp.

Additionally, it was shown in [4] that the condition number $\kappa(\mathcal{R}_{N},\mathcal{S}_{M})$ of the reconstruction method can also be controlled by the SSR, i.e.

[TABLE]

Hence, not only the accuracy but also the stability is controlled by the SSR.

We conclude that along with mappings that map into the reconstruction space $\mathcal{R}_{N}$, generalized sampling is also optimal in terms of achieving a low condition number and high accuracy.

2.2.3. Non-linear extension of generalized sampling

In [34] a consistent approach to generalized sampling was introduced. This means that, in contrast to generalized sampling, the output of the reconstruction method has the same measurements in the sampling space as the input. Let $f\in\mathcal{R}$, then it can be represented as $f=\sum_{i=1}^{\infty}\beta_{i}r_{i}$. The measurements can be written as $l(f)^{[M]}=P_{[M]}U\beta$, where $P_{[M]}$ is the orthogonal projection onto space spanned by the first $M$ elements of the canonical bases of $\ell^{2}(\mathbb{N})$, and $U$ is defined as in (24). The introduced method solves the non-linear minimization problem

[TABLE]

The solution is given by

[TABLE]

The measurements of $EG_{M}(f)$ and $f$ are naturally equal, because of $P_{\mathcal{S}_{M}}EG_{M}(f)=P_{[M]}U\alpha=P_{[M]}U\beta=P_{\mathcal{S}_{M}}f$. Hence, this approach is consistent and maps into the reconstruction space. In contrast to generalized sampling, it is not necessary to decide the number of reconstructed coefficients a priori.

In the setting of arbitrary sampling and reconstruction spaces, it was shown that for every number of reconstructed coefficients $N$ there exists some $M_{0}$ such that the method reconstructs the data correct up to its first $N$ coefficients. Hence, this approach is convergent as the error is given by $\mathcal{O}(||P_{\mathcal{R}_{N}}^{\perp}f||)$, which goes to zero for $M\rightarrow\infty$. The speed of convergence is then dependent upon the sampling and reconstruction space. This speed was analysed in [34] for the case of Fourier measurements and wavelet reconstruction. Let the reconstruction space $\mathcal{R}$ be given by wavelets and the sampling space $\mathcal{S}$ be the space representing Fourier measurements. Then, for $\beta\in\ell_{1}(\mathbb{N})$ and $N\in\mathbb{N}$ the following holds:

1. (1)

If for some $D>0$ and $\gamma\geq 1$, the Fourier transform of the scaling function $\phi$ decays with

[TABLE]

then there exists some constant $C$ independent of $N$ (but dependent on $\gamma$ and $\epsilon$) such that for $M=CN^{1+1/(2\gamma-1)}$, any solution $\alpha$ to (28) satisfies

[TABLE] 2. (2)

If for $k=0,1,2$, for some $D>0$ and $\gamma\geq 1.5$, the Fourier transform of the scaling function $\phi$, the wavelet $\psi$ and their first two derivatives decays with

[TABLE]

then there exists some constant $C$ independent of $N$ (but dependent on $\phi$, $\psi$ and $\epsilon$) such that for $M=CN$, any solution $\alpha$ to (28) satisfies

[TABLE]

For Dauchechies wavelets among other wavelets, the assumptions of this result are satisfied by their construction. The first assumption is fulfilled by all Daubechies wavelets and the second one is fulfilled for Daubechies wavelets of $7$ or more vanishing moments. Numerical experiments suggest that even Daubechies wavelets with less vanishing moments might have a linear relationship [34].

2.2.4. Comparison

In this chapter we want to point out the similarity of generalized sampling and the PBDW-method in terms of optimality and the dependence on the SSR. Additionally, we want to discuss the differences concerning the output space and the consistency. We also include the non-linear extension of generalized sampling into this comparison and examine the overall robustness.

In Equation (21) and Theorem 2 we see that the error bounds of generalized sampling and the PBDW-method depend both on $\mu(\mathcal{R}_{N},\mathcal{S}_{M})$. Moreover, since both systems solve the same least square problem

[TABLE]

as part of the reconstruction, they also have the same condition number $\kappa(\mathcal{R}_{N},\mathcal{S}_{M})$. In more detail, we have that the PBDW-method first solves the generalized sampling problem with a solution $v^{*}$, see Equation (15). Then, the solution is tweaked to be consistent, i.e.

[TABLE]

Moreover, both methods have been proven to be optimal in their setting. This means that linear methods cannot outperform them. Additionally, the bounds are sharp. This underlines the importance of the analysis of the stable sampling rate. Therefore, we dedicate §3 to the investigation of the relation between the dimension of the sampling and the reconstruction space to bound the subspace angle.

Besides the similarities both methods take different parts into consideration for the construction of the algorithms. Generalized sampling ensures that the output is contained in the reconstruction space. This is very desirable when the signal is sparse in the reconstruction bases, because we do not get the artefacts from the measurements in the output signal. Nevertheless, the approach is not consistent. This means that the projection onto the sampling space of the input and output may not be equal.

This is overcome with the non-linear approach presented in [34]. For this gain in consistency one, unfortunately, has to compromise with the robustness. Stability is considered in the $\ell_{1}$ setting, and the used definition is equal to the existence of the condition number of the method. This means that the problem is well-conditioned in terms of solving an $\ell_{1}$ minimization problem rather than the reconstruction of $f$ from $M$ samples. Particularly, this means that the problem is not robust against noise.

Finally, the PBDW method is linear and consistent. This is possible due to the fact that the solution is not forced to stay in the reconstruction space $\mathcal{R}_{N}$. Remark that (21) shows that the error gets smaller with larger $M$ even though we keep the reconstruction space $\mathcal{R}_{N}$ the same. Hence, with increasing $\mathcal{S}_{M}$ we are leaving the reconstruction space $\mathcal{R}_{N}$ and get further away as $M$ increases. This can be desirable, if the artefacts from the sampling space are not too severe, as it allows us a mix of properties from the sampling and the reconstruction space. Nevertheless, if the function is very sparse in the reconstruction space and has a lot of artefacts in the sampling space, this approach leads to less impressive reconstructions. The impact is further illustrated in §4 Numerical experiments.

3. Stable Sampling Rate

3.1. Linearity of the Stable Sampling Rate

In §2.2 we saw that the subspace angle between the sampling and the reconstruction space controls the reconstruction accuracy. Therefore, one is interested in the relation between the number of samples and the number of reconstructed coefficients, such that the subspace angle is bounded. In detail, we are interested in the stable sampling rate:

[TABLE]

The stable sampling rate has been analysed for important cases which appear frequently in practice, i.e. for the Fourier-Wavelet [3], Fourier-Polynomial bases [24] and Walsh-Wavelet cases [23]. For the reconstruction with wavelets we get for both Fourier and Walsh sampling that the stable sampling rate is linear. This is the best relation one can wish for and means that the methods discussed above are, up to a constant, as good as if one could access the wavelet coefficients directly. In particular, for the Walsh case we have the following theorem.

Theorem 3** ([23]).**

Let $\mathcal{S}$ and $\mathcal{R}$ be the sampling and reconstruction space spanned by the d-dimensional Walsh functions and separable boundary wavelets respectively. Moreover, let $N=2^{dR}$ with $R\in\mathbb{N}$. Then for all $\theta\in(1,\infty)$ there exists $S_{\theta}$ such that for all $M\geq 2^{dR}S_{\theta}$ we have $\mu(\mathcal{R}_{N},\mathcal{S}_{M})\leq\theta$. In particular one gets $\Theta\leq S_{\theta}N$. Hence, the relation $\Theta(N;\theta)=\mathcal{O}(N)$ holds for all $\theta\in(1,\infty)$.

A natural question which arises from this theorem is whether it is possible to give sharp bounds on the constant $S_{\theta}$. In Figure 1 we can see the stable sampling rate for different Wavelets and the bound $\theta=2$. The slope $S_{\theta}$ is unknown in most cases and very difficult to find. This comes from the fact that for the majority of wavelets the reconstruction matrix is not perfectly block diagonal, as in 1d and 1f. Hence, one has to take the off diagonals into consideration. The numerics suggest that the slope is higher the further away the reconstruction matrix gets from block diagonal. Only for the case of Haar wavelets and Walsh functions we get that the reconstruction matrix is perfectly block diagonal. This can be seen in 1b. Note, that from the numerical example one may deduce that $S_{\theta}=1$. This is indeed the case, and the analysis detailed below establishes that $S_{\theta}=1$ for all $\theta\in(1,\infty).$

3.2. Sharpness for the Haar wavelet - Walsh case

The sharp bound on $S_{\theta}$ can be summarised in the following theorem.

Theorem 4**.**

Let the sampling space $\mathcal{S}$ be spanned by the Walsh functions and the reconstruction space $\mathcal{R}$ by the Haar wavelets in $L^{2}([0,1]^{d})$. If $N=2^{dR}$ for some $R\in\mathbb{N}$, then for every $\theta\in(1,\infty)$ we have that the stable sampling rate is the identity, i.e. $\Theta(N,\theta)=N$.

For the proof we first study the behaviour of Haar wavelets in one dimension under the Walsh transform. This way we also get a theoretical argument for the block structure that can be seen in the numerical implementation.

Lemma 1**.**

Let $\psi=\mathcal{X}_{[0,1/2]}-\mathcal{X}_{(1/2,1]}$ be the Haar wavelet. Then, we have that

[TABLE]

Proof.

For the scalarproduct we have

[TABLE]

where $\Delta_{k}^{p}=[2^{-p}k,2^{-p}(k+1))$. We know from [6] that the function $\operatorname{Wal}(n,x)$ for $2^{p}\leq n<2^{p+1}$ takes the value $+1$ on the interval $\Delta_{2k}^{p+1}$ or $\Delta_{2k+1}^{p+1}$ and $-1$ on the other one for $k=0,\ldots,2^{p}-1$. Now, we consider three different cases.

Case 1: $n<2^{R}$. There exist $r<R$ such that $2^{r}\leq n<2^{r+1}$. Then, the function $\operatorname{Wal}(n,x)$ is constant on the interval $\Delta^{r}_{k}$ for any $k=0,\ldots,2^{r}-1$. Note that for $j=0,\ldots,2^{R}-1$ we have that the rounding error is bounded as follows

[TABLE]

Then, we have the interval inclusion

[TABLE]

Moreover,

[TABLE]

Hence, $\operatorname{Wal}(n,x)$ takes the same value on $\Delta^{R+1}_{2j}$ and $\Delta^{R+1}_{2j+1}$. Therefore, the two integrals are equal and the scalarproduct vanishes.

Case 2: $2^{R}\leq n<2^{R+1}$. We have for $j=0,\ldots,2^{R}-1$ that $\operatorname{Wal}(n,x)$ is equal to $+1$ on $\Delta_{2j}^{R+1}$ or $\Delta_{2j+1}^{R+1}$ and $-1$ on the other. Therefore, we have either

[TABLE]

Or in the other case analogously

[TABLE]

Now, we are left with the last case.

Case 3: $n\geq 2^{R+1}$. There exists an integer $r\geq R+1$ such that $2^{r}\leq n<2^{r+1}$. This is similar to the first case. Moreover, we have for $j=0,\ldots,2^{R}-1$ that

[TABLE]

With the fact that $\operatorname{Wal}(n,x)$ takes the value $+1$ on one of the invervals $\Delta^{r}_{2l}$ and $\Delta^{r}_{2l+1}$ and $-1$ on the other, we have that the function $\operatorname{Wal}(n,x)$ takes the values $+1$ and $-1$ on half of the interval of $\Delta_{2j}^{R+1}$. Therefore, the integral vanishes. The same holds true for $\Delta_{2j+1}^{R+1}$ such that we get the desired result.

∎

Before we prove Theorem 4, we analyse the reconstruction matrix for the Haar wavelet - Walsh case in two dimensions. The number of functions which span the wavelet space in $d$ dimensions grows exponentially with $d$. Therefore, we restrict ourselves to two dimensions to underline the main idea. In Figure 2 we can see that the reconstruction matrix in two dimensions has an additional structure in each level. Similar to the one-dimensional case we have perfect block structure for the Haar case and nearly block structure for the higher order wavelets. For the analysis of this phenomena in the Haar wavelet - Walsh case we examine the definition and order of the two dimensional Haar wavelet. Note that this is necessary for the analysis of the reconstruction matrix, instead of the SSR. In two dimensions the wavelets are constructed by the tensor product. Due to the multi-resolution analysis we have that $V_{r}^{2}=(V_{r-1}\oplus W_{r-1})^{2}=V_{r-1}^{2}\oplus V_{r-1}\otimes W_{r-1}\oplus W_{r-1}\otimes V_{r-1}\oplus W_{r-1}^{2}$. Hence, the wavelets can be in three different spaces, $V_{r-1}\otimes W_{r-1}$, $W_{r-1}\otimes V_{r-1}$ and $W_{r-1}^{2}$. In the first space wavelets are constructed by the tensor product of the one-dimensional scaling function and the one-dimensional wavelet. For the second one the one-dimensional wavelet and the one-dimensional scaling function are combined by the tensor product. And finally, for the third space we take the tensor products of two one-dimensional wavelets. This results for $R\in\mathbb{N},0\leq j_{1},j_{2}\leq 2^{R}-1$ in

[TABLE]

The scaling function is simply the tensor product of two one-dimensional scaling functions:

[TABLE]

For the order of the reconstruction matrix, we first take the first level scaling function $\phi$. Then, we increase by levels, in each level $R$ we let first $j_{1}$ go from $0,\ldots,2^{R}-1$ and then $j_{2}=0,\ldots,2^{R}-1$. Finally, we let $l=1,\ldots,3$ such that we get for the order of the wavelets: $\phi$, $\psi_{R,0,0,1},\ldots,\psi_{R,2^{R}-1,0,1}$,$\psi_{R,0,1,1},\ldots$ $\psi_{R,2^{R}-1,1,1},$ $\ldots\psi_{R,2^{R}-1,2^{R}-1,1}$, $\psi_{R,0,0,2},\ldots,\psi_{R,2^{R}-1,2^{R}-1,3}$.

Due to the fact that the higher dimensional wavelets are constructed also by means of the scaling functions, it is necessary to analyse the decay rate for the scaling function as well. Moreover, this is also a main ingredient for the proof of Theorem 4 as we represent the union of the wavelet spaces by the scaling space.

Lemma 2**.**

Let $\phi=\mathcal{X}_{[0,1]}$ be the Haar scaling function. Then, we have that the Walsh transform obeys the following block and decay structure

[TABLE]

Proof.

The scalarproduct can be expressed as an integral over the interval $\Delta_{j}^{R}$.

[TABLE]

We look at the two different cases

Case 1: $n<2^{R}\quad$ Remember from Lemma 1 that $\operatorname{Wal}(n,x)$ is constant to $+1$ or $-1$ on the interval $\Delta_{j}^{R}$ for $j=0,\ldots,2^{R}-1$. Hence, we get that

[TABLE]

Case 2: $n\geq 2^{R}\quad$ This follows as in Case 3 of Lemma 1. With the difference that we are looking at the integral over the interval $\Delta_{j}^{R}$ instead of the two integrals $\Delta_{2j}^{R+1}$ and $\Delta_{2j+1}^{R+1}$. Nevertheless, they vanish for the same reason. ∎

With this in hand we can now state the structure of the reconstruction matrix in two dimensions.

Corollary 1**.**

Let $\psi_{R,j_{1},j_{2},l}$ be the Haar wavelet defined as in (49). Then, the Walsh transform has the following property for $0\leq j_{1},j_{2}\leq 2^{R}-1$

[TABLE]

[TABLE]

and for the third version

[TABLE]

Proof.

The proof follows directly from the tensor product structure and Theorem 1 and Lemma 2. ∎

After this detailed analysis of the behaviour of the wavelet and the scaling function under the Walsh transform, we are able to proof Theorem 4.

Proof of Theorem 4.

We want to analyse the subspace angle $\mu(\mathcal{R}_{N},\mathcal{S}_{M})$ for $N=M$. In detail, we are interested in bounding $\mu(\mathcal{R}_{M},\mathcal{S}_{M})$ by $\theta$ for all $\theta\in(1,\infty)$. Hence, we try to show that $\mu(\mathcal{R}_{M},\mathcal{S}_{M})=1$ or equally $1/\mu(\mathcal{R}_{M},\mathcal{S}_{M})=1$ for $M=2^{dR}$. Due to the fact that the circle $\left\{r\in\mathcal{R}_{M},||r||=1\right\}$ is compact, and the orthogonal projection is continuous there exists $r_{0}\in\mathcal{R}_{M}$, $||r_{0}||=1$ such that we have

[TABLE]

The minimal element $r_{0}$ can be represented by

[TABLE]

where the multi-index notation is used. Then, we have that

[TABLE]

With (51) we get that this sum vanishes. Hence,

[TABLE]

as desired. ∎

3.3. Approximation Rate

In this chapter we discuss the approximation rate for Walsh functions and wavelets. With approximation theory we get a good insight in the representation properties of bases. For a given orthonormal basis $\left\{r_{i}\right\}_{i\in\mathbb{N}}$ for $L^{2}([0,1^{d}])$ we have that for all $f\in L^{2}([0,1]^{d})$ it holds

[TABLE]

Unfortunately, in most applications we cannot store or access $\langle f,r_{i}\rangle$ for all $i\in\mathbb{N}$ but only for $i=1,\ldots,N$ for some $N\in\mathbb{N}$. Hence, instead of the true function $f$ we can only have an approximation $f_{N}=\sum_{i=1}^{N}\langle f,r_{i}\rangle r_{i}$. In the field of approximation theory one studies for different representation systems how good this estimate is. Particularly, one is interested in the error

[TABLE]

In general, representation systems are desirable with the property that this error decays very fast in $N$. The reason why is that we get a good representation from only a small amount of information. This results in more efficient compression or less measurement time.

In [6] it is shown that for all continuous functions in $L^{2}([0,1]^{d})$ the approximation error decays with $\mathcal{O}(N^{-2})$. The block artefacts are reflected in the poor approximation rate, and it underlines the need of reconstruction methods, which change the basis to achieve better approximation rates.

Daubechies wavelets obey this property. This is analysed in detail in [33]. Consider the representation with Daubechies wavelets of order $p$, and let the data $f$, which we try to represent, be in the Sobolev space $W^{\gamma}([0,1]^{d})$ for some $\gamma<p$. Then, we get an improved approximation rate of

[TABLE]

Hence, for all Daubechies wavelets of order $p\geq 3$ and all functions $f\in W^{\gamma}([0,1]^{d})$ for some $\gamma\geq 2$ we get an improved decay rate with Daubechies wavelets in contrast to the representation with Walsh functions.

From Theorem 3 we know that the stable sampling rate is linear for Walsh functions and wavelets and in the following chapter we will even see that the constant is reasonably low, i.e. $S_{2}=2$ for Daubechies $8$ wavelets. This means that with only a constant increase of measurements we get from a decay rate of $\mathcal{O}(N^{-2})$ to an approximation rate of $\mathcal{O}(N^{-2\gamma})$.

4. Numerical Experiments

In this section we want to underline the findings from §2.2.4, especially, the differences between generalized sampling and the PBDW technique with examples using both Fourier and Walsh samples. We emphasise that some of the examples are chosen particularly to highlight the differences between the two methods. Hence, the examples may not reflect typical practical scenarios. In Figure 3 we look at an extreme case and consider the task of recovering the Haar wavelet from Fourier coefficients. Generalized sampling will obviously recover the function perfectly, given that we choose the Haar basis for the reconstruction space. In this very special case, the recovered solution is consistent with the samples. Nevertheless, this is not true in general and only possible because one gets perfect recovery. In particular, generalized sampling will usually provide non-consistent solutions. The PBDW-method on the other hand is always consistent. The effect of this is that, even when the reconstruction space is fixed, the solution changes with the number of samples. Moreover, in this example it gets closer, as the number of samples increase, to the solution provided by simply truncating the Fourier series.

The next example in Figure 4 considers Walsh samples. The original function displayed in Figure 4a is continuous with two jump discontinuities. Moreover, the continuous part is very well represented with Daubechies $8$ wavelets, as demonstrated by the generalized sampling reconstruction in Figure 4b. Nevertheless, there are some artefacts at the jumps. Walsh functions represent the jumps better, but lead to heavy artefacts along the continuous part of the function 3c. Hence, both spaces have advantages and disadvantages and can represent the signal well in different areas. The PBDW-method allows us to take advantages from both spaces, as it does not force the solution to stay in the reconstruction space. This leads to a different reconstruction quality as in 4d and 4f. In this case we also get better results with more samples even though we reconstruct the same number of coefficients.

In the last example we used the code from [18] and consider reconstruction of images from Fourier samples. In Figure 5 we see that the PBDW-method provides good results in the $2$D setting and with Fourier measurements. The function is very smooth outside the discontinuous part. Therefore, it can be effectively represented with wavelets in this area. Nevertheless, there are some obvious artefacts around the discontinuity in Figure 5c. The artefacts that arise from the truncated Fourier transform are also easy to spot in Figure 5b and less clear in Figure 5e due to the increased number of samples. With the PBDW-method it is possible to merge the advantages of both systems and decrease the error, as seen in Figures 5d and 5f. We can also see that an increased number of samples leads to better performance of the PBDW-method even if the number of reconstructed coefficients stays the same.

ACKNOWLEDGEMENT

LT acknowledges support by the UK Engineering and Physical Sciences Research Council (EPSRC) grant EP/L016516/1 for the University of Cambridge Centre for Doctoral Training, the Cambridge Centre for Analysis. AH acknowledges support from a Royal Society University Research Fellowship.