Sampling and Reconstruction in Distinct Subspaces

# Sampling and Reconstruction in Distinct Subspaces Using Oblique Projections

Peter Berger Institute of Telecommunications
Vienna University of Technology,
Gusshausstrasse 25/389
A-1040 Vienna, Austria
Karlheinz Gröchenig Faculty of Mathematics
University of Vienna
Oskar-Morgenstern-Platz 1
A-1090 Vienna, Austria
and  Gerald Matz Institute of Telecommunications
Vienna University of Technology,
Gusshausstrasse 25/389
A-1040 Vienna, Austria
###### Abstract.

We study reconstruction operators on a Hilbert space that are exact on a given reconstruction subspace. Among those the reconstruction operator obtained by the least squares fit has the smallest operator norm, and therefore is most stable with respect to noisy measurements. We then construct the operator with the smallest possible quasi-optimality constant, which is the most stable with respect to a systematic error appearing before the sampling process (model uncertainty). We describe how to vary continuously between the two reconstruction methods, so that we can trade stability for quasi-optimality. As an application we study the reconstruction of a compactly supported function from nonuniform samples of its Fourier transform.

###### 2000 Mathematics Subject Classification:
Funding by the Austrian Science Fund (FWF) through grant NFN SISE (S10602) and P26273 - N25 and by the Vienna Science and Technology Fund (WWTF) through project VRG12-009 and project ICT15-119.

## 1. Introduction

### 1.1. The reconstruction problem

In this paper we treat the following sampling problem. Let be a separable Hilbert space over with inner product .We assume that we are given linear measurements , , of an unknown function . We call the sampling frame and the sampling space. Our goal is to approximate the function by an element in the reconstruction space with , by a series expansion from the given measurements. The main point is that in general the reconstruction space is distinct from the sampling space, whereas in classical frame theory these two spaces coincide.

### 1.2. Areas of application and related work

This type of sampling problem arises in many concrete applications and in the numerical modelling of infinite dimensional problems.

(i) Sampling of bandlimited functions. In [25] a bandlimited function is approximated from finitely many, nonuniform samples by means of a trigonometric polynomial. In this case the sampling space consists of the reproducing kernels and the reconstruction vectors are .

(ii) Inverse Polynomial Reconstruction Method. In this method one tries to approximate an algebraic polynomial or an analytic function from its Fourier samples. Thus the sampling space consists of vectors , and the reconstruction space consists of a suitable polynomial basis, usually the monomials , or the Legendre polynomials. This method claims to efficiently mitigate the Gibbs phenomenon  [38, 28, 29, 30], and, indeed, the modified inverse polynomial reconstruction method [27] leads to a numerically stable reconstruction when .

(iii) Fourier sampling. More generally, the goal is to approximate a compactly supported function in some smoothness class from its nonuniform Fourier samples . Thus the sampling space consists again of the functions . The reconstruction space depends on the signal model and on a priori information. If is smooth and belongs to a Besov space, then the reconstruction space may be taken to be a wavelet subspace. The problems of Fourier sampling have motivated Adcock and Hansen to revisit nonuniform sampling theory and to create the impressive and useful framework of generalized sampling [8, 7, 10, 6, 33].

(iv) Model reduction in parametric partial differential equations and the generalized empirical interpolation method. In general the solution manifold to a parametric partial differential equation is quite complicated, therefore it is approximated by finite-dimensional spaces . The Generalized Empirical Interpolation Method (GEIM) [34, 35] builds an interpolant in an -dimensional space based on the knowledge of physical measurements . In [36, 13] an extension based on a least squares method has been proposed, where the dimension of is smaller than the number of the measurements . A further generalization to Banach spaces is contained in [18]. The focus in [36, 13] lies in minimizing the error caused by the model mismatch. This is done by using a correction term outside of the reconstruction space, which means that (in contrast to our work) the reconstruction is allowed to be located outside of the reconstruction space. This approach is optimal in the absence of measurement noise [13].

In all these problems the canonical approximation or reconstruction is by means of a least squares fit, namely

 ~f=arg~{}ming∈T∑j∈Nwj∣∣⟨g,uj⟩H−dj∣∣2. (1)

The weights are usually chosen to be , but in many contexts is has turned out to be useful to use weights as a kind of cheap preconditioners. The use of adaptive weights in sampling theory goes back at least to [23, 24], and has become standard in the recent work on (Fourier) sampling, see for example [5, 25, 26, 40, 1, 3, 2, 4, 11].

### 1.3. The reconstruction operators

In this paper we restrict ourselves to the case where the approximation of the unknown function is obtained by a linear and bounded reconstruction operator . Thus the approximation from the data is given by . We will use two quantities to measure the quality of such a reconstruction operator. As a measure of stability with respect to measurement noise we use the operator norm . As a measure of stability with respect to model mismatch we follow [9] and use the so-called quasi-optimality constant (see Definition 2.1).

Let denote the orthogonal projection onto , the target function, and be the noise vector. Then the input data are given by the sequence , the reconstruction is , and the error is bounded by

 ∥f−Q((⟨f,uj⟩H+lj)j∈N)∥H⩽μ(Q)∥f−PTf∥H+∥Q∥op∥l∥2. (2)

### 1.4. Contributions

The error bound (2) raises several questions:

• Which operators admit an error bound of the form (2)?

• Under what circumstances does such an operator exist?

• Which operator has the smallest possible operator norm ?

• Which operator has the smallest possible quasi-optimality constant ?

• Is there a way to trade-off between quasi-optimality and operator norm?

Our objective is to answer these questions both in finite-dimensional and infinite dimensional Hilbert spaces. The results can be formulated conveniently in the language of frame theory.

(i) Characterization of all reconstruction operators. We characterize all reconstruction operators that admit an error estimate of the form (2). In fact, every dual frame of the set yields a reconstruction satisfying (2). Conversely, every reconstruction operator subject to (2) is the synthesis operator of a dual frame of . Note that (2) implies that such a reconstruction operator is exact on the reconstruction space, i.e., for all . Reconstruction operators fulfilling this property are called perfect. For a precise formulation see Theorem 2.2.

The important insight of  [9] is the connection between stability and the angle between the sampling space and the reconstruction space. We will see that a perfect reconstruction operator exists if and only if . It should also be mentioned that the reconstruction operators considered in this paper are a special case of pseudoframes [32].

(ii) Least squares approximation. As already mentioned, the canonical approximation of the data by a vector in is by a least squares fit. Let denote the analysis operator of the frame and let denote the vector containing the noisy measurements . Let the reconstruction operator be defined by the least squares fit

 Q1d=arg~{}ming∈T∑j∈N∣∣⟨g,uj⟩H−dj∣∣2=arg~{}ming∈T∥U∗g−d∥2. (3)

It is folklore that the least squares solution (3) is optimal in the absence of additional information on . Precise formulations of this optimality were proven in [9, Theorem 6.2.] (including even non-linear reconstructions) and in [12] (in abstract Hilbert space). We will show in addition (Theorem 3.1) that is the synthesis operator of the canonical dual frame of . Using this property we derive a simple proof for the statement that the operator has the smallest possible operator norm among all perfect reconstruction operators.

(iii) Minimizing the quasi-optimality constant. Let be the square root of the Moore-Penrose pseudoinverse of the Gramian of the sampling frame and consider the operator defined by

 Q0d=arg~{}ming∈T∥WU∗g−Wd∥2. (4)

We will show (Theorem 3.3) that has the smallest possible quasi-optimality constant. The reduction of the quasi-optimality constant is one of the motivations of weighted least squares, see [2, 1, 5, 3]. In (4) we go a step further and use the non-diagonal matrix as a weight for the least squares problem. From the point of view of linear algebra, may be seen as a preconditioner.

In [2, 1, 5, 3] and also [11, 24, 25, 26, 4, 23] the stability with respect to a bias in the measured object is considered, i.e., the reconstruction from (stated in terms of a frame inequality in the latter). In this context, is the most stable operator with respect to biased objects, see the discussion in Section 3.6.

(iv) Trading stability and quasi-stability. It is natural to ask whether one can mix between the two least squares problems (3) and (4). Let and and define by

 Qλd=arg~{}ming∈T∥Σ−12λU∗g−Σ−12λd∥2.

These reconstruction operators “interpolate” between (most stable with respect to noise) and (most stable with respect to model uncertainty). The parameter can be seen as a regularization parameter, or alternatively the matrix as version of the adaptive weights in sampling. In Theorem 3.7 and Lemma 3.8 we will study this class of reconstruction operators and derive several representations for .

(v) Fourier resampling — numerical experiments. In the last part we carry out a numerical comparison of the various reconstruction operators on the basis of the so-called resampling problem. We approximate a function with compact support from finitely many, nonuniform samples of its Fourier transform and then resample the Fourier transform on a regular grid. For this problem we test the performance of the reconstruction operators .

The paper is organized as follows: In Section 2 we introduce the frame theoretic background, discuss the angle between subspaces, and characterize all reconstruction operators satisfying the required stability estimate (2). In Section 3 we study the various least squares problems (3) and (4) and analyze several representations of the corresponding reconstruction operators. The section is complemented by general numerical considerations. Section 4 covers the numerical experiments on Fourier sampling. The brief appendix collects some standard facts about frames.

## 2. Classification of all reconstruction operators

We will use the language of frame theory throughout the whole paper. The Appendix contains a short list of basic definitions and well known facts from frame theory. For more details on this topic, see for instance [15].

Let us introduce some notation. To every set of measurement vectors in a Hilbert space (of finite or infinite dimension) we associate the synthesis operator defined formally by and the sampling space . The adjoint operator consists of the measurements and is called the analysis operator. The frame operator is and the Gramian is . With this notation, is a frame for , if there exist constants , such that for every

 A∥f∥2H⩽∥U∗f∥22⩽B∥f∥2H.

We always assume that is a frame for , thus is bounded from to and has closed range in . We use for the range of an operator and for its kernel (null space).

Likewise we assume that is a frame for the reconstruction space with synthesis operator and analysis operator . Thus

 C∥g∥2H⩽∥T∗g∥22⩽D∥g∥2H for g∈T.

Given a sequence of linear measurements , we try to find an approximation of in the subspace . Assuming that all occurring operartors are bounded, we investigate the class of reconstruction operators , such that is the desired reconstruction or approximation of . We use two metrics to quantify the stability of a reconstruction operator . As a measure for stability with respect to measurement noise we use the operator norm . In order to measure how well deals with the part of the function lying outside of the reconstruction space, we use the quasi-optimality constant from [9].

###### Definition 2.1.

Let and be the orthogonal projection onto . The quasi-optimality constant is the smallest number , such that

 ∥f−QU∗f∥H⩽μ∥f−PTf∥H,for all f∈H.

If we call a quasi-optimal operator. Since is the element of closest to , the quasi-optimality constant is a measure of how well performs in comparison to orthogonal projection . Note that for we have , thus a quasi-optimal reconstruction operator is perfect.

The following theorem characterizes all bounded quasi-optimal operators.

###### Theorem 2.2.

Let and be closed subspaces of , and be a Bessel sequence spanning the closed subspace . For an operator the following are equivalent.

1. There exist constants , such that for and

 ∥f−Q(U∗f+l)∥H⩽μ∥f−PTf∥H+β∥l∥2. (5)
2. for and is a bounded operator.

3. The sequence is a frame for . Let be a dual frame of , then is of the form

 Qc=∑j∈Ncjhj,

i.e., is the synthesis operator of some dual frame of .

4. The operator is bounded and is a bounded oblique projection onto .

Theorem 2.2 sets up a bijection between the class of reconstruction operators and the class of all dual frames of .

To prove Theorem 2.2, we need the concept of subspace angles. Among the many different definitions of the angle between subspaces (see [41, 39]) the following definition is most suitable for our analysis.

###### Definition 2.3.

Let and be closed subspaces of a Hilbert space . The subspace angle between and is defined as

 cos(φT,U)=infg∈T∥g∥H=1∥PUg∥H=infg∈T∥g∥H=1 sup∥u∥H=1u∈U|⟨g,u⟩H|. (6)

We observe that in general . For , and therefore . If , then and .

The following lemma collects the main properties of oblique projections and angles between subspaces.

###### Lemma 2.4.

Assume that and are closed subspaces of a Hilbert space . Then

1. if and only if and the direct sum (not necessarily orthogonal) is closed in .

2. If and is a closed subspace of , then the oblique projection with range and kernel is well defined and bounded on .

3. Let , , and let be the oblique projection with range and null space . Then

 ∥PT,W∥op=1cos(φT,W⊥)

and

 ∥f−PTf∥H⩽∥f−PT,Wf∥H⩽1cos(φT,W⊥)∥f−PTf∥H, (7)

for all . The upper bound in (7) is sharp.

Item (i) of Lemma 2.4 is stated in [42, Theorem 2.1], the proof of (ii) can be found in [14, Theorem 1], and for (iii) see [41], [14], and [9, Corollary 3.5].

Proof of Theorem 2.2 (i) (ii). Set and choose . Then (5) implies , since otherwise . Setting in (5) implies that is bounded.

(ii) (iii) Let be a bounded operator with for . Let be the standard basis of and let . Then . In particular for ,

 QU∗g=∑j∈N⟨g,PTuj⟩Hhj=g.

Since is bounded, is a Bessel sequence in . By assumption is a Bessel sequence in with Bessel bound and consequently

 ∑j∈N|⟨f,PTuj⟩H|2=∑j∈N|⟨PTf,uj⟩H|2⩽B∥PTf∥2H⩽B∥f∥2H.

Therefore is a dual frame of .

(iii) (iv) Let be a dual frame of and define by and . Since the range of is contained in and for , it follows that is onto and that . Since both and are bounded, is bounded.

(iv) (i). Let be a bounded operator, and let be a bounded oblique projection onto . Lemma 2.43 implies that , and consequently

 ∥f−Q(U∗f+l)∥H⩽∥P∥op∥f−PTf∥H+∥Q∥op∥l∥2.

This finishes the proof.

As a direct consequence of Theorem 2.2 and Lemma 2.4, 3, we obtain the following characterization of the quasi-optimality constant.

###### Corollary 2.5.

If is a bounded and perfect reconstruction operator, then is a bounded oblique projection onto . If denotes the null-space of , then

 μ(Q)=∥QU∗∥op=1cos(φT,W).

In the following we always use the assumption that the angle between the reconstruction and sampling space fulfills . The following lemma shows that this assumption is equivalent to forming a frame for for every frame for . By Theorem 2.2 3 this is necessary for the existence of a quasi-optimal operator. In finite dimensions, for a basis for , can only be a spanning set for if . This means that by the assumption we restrict ourselves to an oversampled regime.

###### Lemma 2.6.

If and are closed subspaces of , then the following are equivalent:

1. .

2. For every frame for with frame bounds and , the projection is a frame for with frame bounds and .

If one of these conditions is satisfied, then the following property holds:

3. , therefore both and are closed subspaces and is pseudo-invertible. Furthermore,

 N(U∗)∩T={0}. (8)
###### Proof.

(i) (ii) Let be a frame for with frame bounds and . The assumption and the definition of imply that

 ∥g∥Hcos(φT,U)⩽∥PUg∥Hfor all g∈T. (9)

In particular, for we obtain with (9)

 A∥g∥2Hcos2(φT,U)⩽A∥PUg∥2H⩽∑j∈N|⟨PUg,uj⟩H|2⩽B∥PUg∥2H⩽B∥g∥2H. (10)

The identity for and now shows that is a frame for with frame bounds and .

(ii) (i) Let be a frame for with upper frame bound and let be a frame for with lower frame bound . Since for , we obtain

 C1∥g∥2H⩽∑j∈N|⟨g,PTuj⟩H|2=∑j∈N|⟨PUg,uj⟩H|2⩽B∥PUg∥2H.

This implies that .

(ii) (iii) Since and are Bessel sequences, both and are bounded, and therefore is also bounded. The entries of are given by

 (U∗T)(j,k)=⟨tk,uj⟩H=⟨tk,PTuj⟩H,

and is a cross-Gramian of two frames for . Let be a dual frame of . Setting we obtain, for ,

 (T∗Uc)k=∑j∈N⟨f,~uj⟩H⟨PTuj,tk⟩H=⟨f,tk⟩H=(T∗f)k.

It follows that

 R(T∗U)=R(T∗). (11)

Since is a frame for , is closed in , and so are and . This implies that both and possess a pseudoinverse (see Appendix 4.1).

To prove (8), let . Then for some and . This means that . Consequently, , and . ∎

## 3. The reconstruction operators

### 3.1. Least squares and the operator Q1

We first consider the reconstruction operator corresponding to the solution of the least squares problem

 Q1d=arg~{}ming∈T∑j∈N∣∣⟨g,uj⟩H−dj∣∣2=arg~{}ming∈T∥U∗g−d∥2. (12)

This approach is analyzed in detail in [9]. Least square approximation is by far the most frequent approximation method in applications and of fundamental importance, since it has the smallest operator norm among all perfect operators.

The following theorem reviews several representations of the operator . The connection of the operator to the oblique projection was already derived in [9, Section 4.1.] for finite dimensional space . Our new contribution is the connection to the canonical dual frame and the systematic discussion of the various representions of a least squares problem. As we will apply the statement several times, we include a streamlined proof. As usual, denotes the Moore-Penrose pseudo-inverse of an operator . For the existence of it suffices to show that the range of is closed (see Appendix 4.1).

###### Theorem 3.1.

Let and be closed subspaces of a Hilbert space such that . Let be a frame for with synthesis operator and frame operator . Let be a frame for with synthesis operator .

Consider the following operators:

1. .

2. The operator is given on by

 A2U∗=PT,S(T)⊥ (13)

and on by

 A2c=0 for c∈R(U∗)⊥. (14)

By (13) is independent of the particular choice of the reconstruction frame for .

3. Let be the canonical dual frame of and be the synthesis operator of .

4. Let and let be the unique minimal norm element of the set

 K:=arg~{}minc∈ℓ2(N) ∥U∗Tc−d∥2. (15)

Let the operator be defined by .

Then all four operators are equal, and provide the unique solution to the least squares problem

 Q1d=arg~{}ming∈T ∑j∈N|⟨g,uj⟩H−dj|2=g∈Targ~{}min∥U∗g−d∥22.
###### Proof.

Step 1. First we check that each , is well defined from to . For this is clear by virtue of Lemma 2.6.

For we need to show that the projection is well defined and bounded on the whole space . According to Lemma 2.4(i) we need to verify that is closed, and that . For this we exploit the frame inequality (10) of , (9), and the fact that , and we obtain

 Acos(φT,U)∥g∥H⩽A∥PUg∥H⩽∥SPUg∥H=∥Sg∥Hfor g∈T. (16)

The lower bound implies that is closed. For the angle we obtain

 (17)

Since and , we continue (17) as follows:

 cos(φT,S(T))⩾infg∈Tg≠0 A∥PUg∥2HB∥g∥H∥PUg∥H=ABcos(φT,U)>0. (18)

It remains to prove that , or, equivalently, that

 (T⊕S(T)⊥)⊥=T⊥∩¯¯¯¯¯¯¯¯¯¯¯¯S(T)=T⊥∩S(T)={0}.

So assume that . Since , we may write every as for some . In particular, there exist and , such that . Then for all , the element satisfies

 0=⟨g,t⟩=⟨STc,Td⟩=⟨UU∗Tc,Td⟩=⟨U∗Tc,U∗Td⟩.

Setting , we obtain . By Lemma 2.6(iii) , and thus , which implies that .

The operator is the synthesis operator with respect to the canonical dual frame of and is therefore bounded by general frame theory.

Now to : By Lemma 2.6 the operator has a closed range and therefore its Moore-Penrose pseudoinverse is well defined. It is well known that is the unique element of of minimal norm. Consequently, is bounded on .

Step 2. We next show that all these operators are equal.

Claim . Since is the unique element of of minimal norm and , we have .

Claim . We define and show that , and . The equality follows from the identity for the Moore-Penrose pseudoinverse applied to . Clearly . To prove the converse inclusion we show that . Using from Lemma 2.6 and we conclude that

 RT=T(U∗T)†U∗T=TPR(T∗U)=TPR(T∗)=TPN(T)⊥=T,

which proves .

Now let , then we have, for all ,

 ⟨T(U∗T)†U∗f,h⟩=⟨(U∗T)†U∗f,T∗h⟩=0.

Since by Lemma 2.6(iii), this means that for all

 0 =⟨(U∗T)†U∗f,T∗Uc⟩=⟨U∗T(U∗T)†U∗f,c⟩ =⟨PR(U∗T)U∗f,c⟩=⟨f,UPR(U∗T)c⟩ =⟨f,UU∗Tc⟩.

In other words, , as claimed.

Claim . We need to show that the operator is the synthesis operator of the canonical dual frame of . The frame operator of can be written in the form

 ~Sf=∑j∈N⟨f,PTuj⟩HPTuj=PT(∑j∈N⟨PTf,uj⟩Huj)=PTUU∗PTf.

By Definition 4.2(iv), the canonical dual frame of is given by with synthesis operator

 A3c=∑j∈Ncj(PTUU∗PT)†PTuj=(PTUU∗PT)†PTUc=(U∗PT)†c, (19)

where we used with for the last equality. Since we have already proved that , we know that the operator is independent of the particular choice of a frame for . We may therefore use the frame with synthesis operator instead of , and as a consequence obtain that , where now we use with . Comparing with (19), we have proved that .

Step 3. Finally we show that each operator provides the unique solution to the least squares fit (12). Since by Lemma 2.6(iii), the solution of the least squares problem

 ~f=arg~{}ming∈T ∥U∗g−d∥22 (20)

is unique. Since , there exists a , such that , and by (20) for every element (cf. (15)). In particular, for the minimal norm element used for the definition of the operator , we obtain