stopping criterion

r

|p1; k (tj ) ’ bj |2 62 ||b ||2 ; (34)

j=1

we take p1; k— as ÿnal approximation. Otherwise we proceed to the next level M =2 and run Algorithm

3:1 again, using p1; k— as initial approximation by setting p2; 0 = p1; k— .

At level M = N the inner level-dependent stopping criterion becomes

r

|pN; k (tj ) ’ bj |2 62 ( ||b || + ||QN f ’ f||)||b ||; (35)

j=1

while the outer stopping criterion does not change since it is level-independent.

Stopping rule (35) guarantees that the iterates of CG do not diverge. It also ensures that CG does

not iterate too long at a certain level, since if M is too small further iterations at this level will

not lead to a signiÿcant improvement. Therefore, we switch to the next level. The outer stopping

criterion (34) controls over-ÿt and under-ÿt of the data, since in presence of noisy data is does not

make sense to ask for a solution pM that satisÿes r |pM (tj ) ’ bj |2 = 0.

j=1

Since the original signal f is not known, the expression ||f ’ QN f|| in (35) cannot be computed.

In [34] the reader can ÿnd an approach to estimate ||f ’ QN f|| recursively.

310 T. Strohmer / Journal of Computational and Applied Mathematics 122 (2000) 297“316

4.2. Solution of ill-conditioned sampling problems

A variety of conditions on the sampling points {tj }j∈Z are known under which the set {sinc(· ’

tj )}j∈Z is a frame for B, which in turn implies (at least theoretically) perfect reconstruction of

a signal f from its samples f(tj ). This does however not guarantee a stable reconstruction from

a numerical viewpoint, since the ratio of the frame bounds B=A can still be extremely large and

therefore the frame operator S can be ill-conditioned. This may happen for instance if in (30)

goes to 1, in which case cond(T ) may become large. The sampling problem may also become

numerically unstable or even ill-posed, if the sampling set has large gaps, which is very common

in astronomy and geophysics. Note that in this case the instability of the system TM aM = yM does

not result from an inadequate discretization of the inÿnite-dimensional problem.

There exists a large number of (circulant) Toeplitz preconditioners that could be applied to the

system TM aM = yM , however it turns out that they do not improve the stability of the problem in

this case. The reason lies in the distribution of the eigenvalues of TM , as we will see below.

Following [38], we call two sequences of real numbers { (n) }n and { (n) }n equally distributed,

k=1 k=1

if

n

1 (n) (n)

lim [F( k) ’ F( k )] =0 (36)

n’∞ n

k=1

for any continuous function F with compact support. 1

Let C be an (n—n) circulant matrix with ÿrst column (c0 ; : : : ; cn’1 ), we write C =circ (c0 ; : : : ; cn’1 ).

√ n’1

The eigenvalues of C are distributed as k = (1= n) l=0 cl e2 ikl=n . Observe that the Toeplitz matrix

An with ÿrst column (a0 ; a1 ; : : : ; an ) can be embedded in the circulant matrix

Cn = circ (a0 ; a1 ; : : : ; an ; an ; : : : ; a1 ): (37)

Theorems 4:1 and 4:2 in [38] state that the eigenvalues of An and Cn are equally distributed as f(x)

where

∞

ak e2 ikx

f(x) = : (38)

k=’∞

The partial sum of the series (38) is

n

ak e2 ikx

fn (x) = : (39)

k=’n

To understand the clustering behavior of the eigenvalues of TM in case of sampling sets with large

gaps, we consider a sampling set in [ ’ M; M ), that consists of one large block of samples and one

large gap, i.e., tj = j=Lm for j = ’mM; : : : mM; for m; L ∈ N. (Recall that we identify the interval

with the torus). Then the entries zk of the Toeplitz matrix TM of (28) (with wj = 1) are

mM

1

e’2 ikj=Lm=(2M +1)

zk = ; k = 0; : : : ; 2M:

2M + 1 j=’mM

(n) (n)

1

In H. Weyl™s deÿnition and are required to belong to a common interval.

k k

T. Strohmer / Journal of Computational and Applied Mathematics 122 (2000) 297“316 311

To investigate the clustering behavior of the eigenvalues of TM for M ’ ∞, we embed TM in a

circulant matrix CM as in (37). Then (39) becomes

mM mM

1

e2 il[k=(4M +1)’j=((2M +1)mL)]

fmM (x) = (40)

Lm(2M + 1) l=’mM j=’mM

whence fmM ’ 1[’1=(2L); 1=(2L)] for M ’ ∞, where 1[’a; a] (x) = 1, if ’a ¡ x ¡ a and 0 else.

Thus the eigenvalues of TM are asymptotically clustered around zero and one. For general nonuni-

form sampling sets with large gaps the clustering at 1 will disappear, but of course the spectral

cluster at 0 will remain. In this case it is known that the preconditioned problem will still have a

spectral cluster at the origin [43] and preconditioning will not be e cient.

Fortunately, there are other possibilities to obtain a stabilized solution of TM aM =yM . The condition

number of TM essentially depends on the ratio of the maximal gap in the sampling set to the

Nyquist rate, which in turn depends on the bandwidth of the signal. We can improve the stability

of the system by adapting the degree M of the approximation accordingly. Thus the parameter

M serves as a regularization parameter that balances stability and accuracy of the solution. This

technique can be seen as a speciÿc realization of regularization by projection, see [12, Chapter 3].

In addition, as described in Section 4.2, we can utilize CG as regularization method for the solution

of the Toeplitz system in order to balance approximation error and propagated error. The multilevel

method introduced in Section 4.1 combines both features. By optimizing the level (bandwidth) and

the number of iterations in each level it provides an e cient and robust regularization technique for

ill-conditioned sampling problems. See Section 5 for numerical examples.

5. Applications

We present two numerical examples to demonstrate the performance of the described methods.

The ÿrst one concerns a 1-D reconstruction problem arising in spectroscopy. In the second example

we approximate the Earth™s magnetic ÿeld from noisy scattered data.

5.1. An example from spectroscopy

The original spectroscopy signal f is known at 1024 regularly spaced points tj . This discrete

sampling sequence will play the role of the original continuous signal. To simulate the situation

of a typical experiment in spectroscopy we consider only 107 randomly chosen sampling values

of the given sampling set. Furthermore, we add noise to the samples with noise level (normalized

by division by 1024 |f(tj )|2 ) of = 0:1. Since the samples are contaminated by noise, we cannot

k=1

expect to recover the (discrete) signal f completely. The bandwidth is approximately = 5 which

translates into a polynomial degree of M ≈ 30. Note that in general and (hence M ) may not be

available. We will also consider this situation, but in the ÿrst experiments we assume that we know

. The error between the original signal f and an approximation fn is measured by computing

||f ’ fn ||2 =||f||2 .

First we apply the truncated frame method with regularized SVD as described in Section 2. We

choose the truncation level for the SVD via formula (21). This is the optimal truncation level in this

case, providing an approximation with least-squares error 0:0944. Fig. 1(a) shows the reconstructed

312 T. Strohmer / Journal of Computational and Applied Mathematics 122 (2000) 297“316

Fig. 1. Example from spectroscopy-comparison of reconstruction methods. (a) Truncated frame method with TSVD,

error=0:0944. (b) Truncated frame method with CG, error=0:1097. (c) Algorithm 3.1 with “correct” bandwidth,

error=0:0876. (d) Using a too small bandwidth, error=0:4645. (e) Using a too large bandwidth, error=0:2412. (f) Mul-

tilevel algorithm, error=0:0959.

T. Strohmer / Journal of Computational and Applied Mathematics 122 (2000) 297“316 313

signal together with the original signal and the noisy samples. Without regularization we get a much

worse “reconstruction” (which is not displayed).

We apply CG to the truncated frame method, as proposed in Section 2.2 with stopping criterion

(22) (for = 1). The algorithm terminates already after 3 iterations. The reconstruction error is with

0:1097 slightly higher than for truncated SVD (see also Fig. 1(b)), but the computational e ort is

much smaller.

Also Algorithm 3:1 (with M = 30) terminates after 3 iterations. The reconstruction is shown in

Fig. 1(c), the least squares error (0:0876) is slightly smaller than for the truncated frame method,

the computational e ort is signiÿcantly smaller.

We also simulate the situation where the bandwidth is not known a priori and demonstrate the

importance of a good estimate of the bandwidth. We apply Algorithm 3:1 using a too small degree

(M = 11) and a too high degree (M = 40). (We get qualitatively the same results using the truncated

frame method when using a too small or too large bandwidth.) The approximations are shown in

Figs. 1(d) and (e). The approximation errors are 0:4648 and 0:2805, respectively. Now we apply

the multilevel algorithm of Section 4.1 which does not require any initial choice of the degree M .

The algorithm terminates at “level” M = 22, the approximation is displayed in Fig. 1(f), the error

is 0:0959, thus within the error bound , as desired. Hence without requiring explicit information

about the bandwidth, we are able to obtain the same accuracy as for the methods above.

5.2. Approximation of geophysical potential ÿelds

Exploration geophysics relies on surveys of the Earth™s magnetic ÿeld for the detection of anoma-

lies which reveal underlying geological features. Geophysical potential ÿeld data are generally ob-

served at scattered sampling points. Geoscientists, used to looking at their measurements on maps or

proÿles and aiming at further processing, therefore need a representation of the originally irregularly

spaced data at a regular grid.

The reconstruction of a 2-D signal from its scattered data is thus one of the ÿrst and crucial steps

in geophysical data analysis, and a number of practical constraints such as measurement errors and

the huge amount of data make the development of reliable reconstruction methods a di cult task.

ˆ

It is known that the Fourier transform of a geophysical potential ÿeld f has decay |f(!)| =

O(e’|!| ). This rapid decay implies that f can be very well approximated by band-limited functions

[30]. Since, in general, we may not know the (essential) bandwidth of f, we can use the multilevel

algorithm proposed in Section 4.1 to reconstruct f.

The multilevel algorithm also takes care of following problem. Geophysical sampling sets are often

highly anisotropic and large gaps in the sampling geometry are very common. The large gaps in the

sampling set can make the reconstruction problem ill-conditioned or even ill-posed. As outlined in

Section 4.2 the multilevel algorithm iteratively determines the optimal bandwidth that balances the

stability and accuracy of the solution.

Fig. 2(a) shows a synthetic gravitational anomaly f. The spectrum of f decays exponentially,

thus the anomaly can be well represented by a band-limited function, using a “cut-o -level” of

|f(!)|60:01 for the essential bandwidth of f.

We have sampled the signal at 1000 points (uj ; vj ) and added 5% random noise to the sampling

values f(uj ; vj ). The sampling geometry “ shown in Fig. 2 as black dots “ exhibits several features

one encounters frequently in exploration geophysics [30]. The essential bandwidth of f would imply

314 T. Strohmer / Journal of Computational and Applied Mathematics 122 (2000) 297“316

Fig. 2. Approximation of synthetic gravity anomaly from 1000 nonuniformly spaced noisy samples by the multilevel

algorithm of Section 4.1. The algorithm iteratively determines the optimal bandwidth (i.e. level) for the approximation.

(a) Contour map of synthetic gravity anomaly, gravity is in mGal. (b) Sampling set and synthetic gravity anomaly. (c)

Approximation by multi-level algorithm. (d) Error between approximation and actual anomaly.

to choose a polynomial degree of M = 12 (i.e., (2M + 1)2 = 625 spectral coe cients). With this

choice of M the corresponding block Toeplitz matrix TM would become ill-conditioned, making

the reconstruction problem unstable. As mentioned above, in practice we usually do not know the

essential bandwidth of f. Hence we will not make use of this knowledge in order to approximate f.

We apply the multilevel method to reconstruct the signal, using only the sampling points {(uj ; vj )},

the samples {f (uj ; vj )} and the noise level =0:05 as a priori information. The algorithm terminates

at level M = 7. The reconstruction is displayed in Fig. 2(c), the error between the true signal and

the approximation is shown in Fig. 2(d). The reconstruction error is 0:0517 (or 0:193 mGal), thus

of the same order as the data error, as desired.

T. Strohmer / Journal of Computational and Applied Mathematics 122 (2000) 297“316 315

References

[1] J. Benedetto, Irregular sampling and frames, in: C.K. Chui (Ed.), Wavelets: A Tutorial in Theory and Applications,

Academic Press, New York, 1992, pp. 445“507.

[2] J. Benedetto, W. Heller, Irregular sampling and the theory of frames, I, Math. Notes 10 (1990) 103“125.

[3] A. Beurling, P. Malliavin, On the closure of characters and the zeros of entire functions, Acta Math. 118 (1967)

79“93.

[4] G. Beylkin, On the fast Fourier transform of functions with singularities, Appl. Comp. Harm. Anal. 2 (4) (1995)

363“381.

[5] A.A. Bjork, V. Pereyra, Solution of Vandermonde systems of equations, Math. Comp. 24 (1970) 893“903.

[6] P.L. Butzer, W. Splettsto er, R.L. Stens, The sampling theorem and linear prediction in signal analysis. Jahresbericht

der DMV 90, 1988, pp. 1“70.

[7] P.G. Casazza, O. Christensen, Approximation of the inverse frame operator and applications to Weyl-Heisenberg

frames, J. Approx. Theory, accepted for publication.

[8] R. Chan, M. Ng, Conjugate gradient methods for Toeplitz systems, SIAM Rev. 38 (3) (1996) 427“482.

[9] O. Christensen, Frames containing Riesz bases and approximation of the frame coe cients using ÿnite dimensional

methods, J. Math. Anal. Appl. 199 (1996) 256“270.

[10] O. Christensen, Moment problems and stability results for frames with applications to irregular sampling and Gabor

frames, Appl. Comp. Harm. Anal. 3 (1) (1996) 82“86.

[11] R. Du n, A. Schae er, A class of nonharmonic Fourier series, Trans. Amer. Math. Soc. 72 (1952) 341“366.

[12] H.W. Engl, M. Hanke, A. Neubauer, Regularization of Inverse Problems, Kluwer Academic Publishers, Dordrecht,

1996.

[13] H.G. Feichtinger, Coherent frames and irregular sampling, Proceedings of the NATO Conference on Recent Advances

in Fourier Analysis and its Applications, Pisa, NATO ASI Series C, Vol. 315, 1989, pp. 427“ 440.

[14] H.G. Feichtinger, K. Grochenig, T. Strohmer, E cient numerical methods in non-uniform sampling theory, Numer.

Math. 69 (1995) 423“440.

[15] H.G. Feichtinger, K.H. Grochenig, Theory and practice of irregular sampling, in: J. Benedetto, M. Frazier (Eds.),

Wavelets: Mathematics and Applications, CRC Press, Boca Raton, FL, 1994, pp. 305“363.

[16] K.M. Flornes, Y.I. Lyubarskii, K. Seip, A direct interpolation method for irregular sampling, Appl. Comp. Harm.

Anal. 7 (3) (1999) 305“314.

[17] I.C. Gohberg, I.A. Fel™dman, Convolution Equations and Projection Methods for Their Solution, American

Mathematical Society, Providence, RI, 1974 (Translated from the Russian by F.M. Goldware, Translations of

Mathematical Monographs, Vol. 41.).

[18] G.H. Golub, C.F. van Loan, Matrix Computations, 3rd Edition, Johns Hopkins, London, Baltimore, 1996.

[19] K. Grochenig, A discrete theory of irregular sampling, Linear Algebra Appl. 193 (1993) 129“150.

[20] K. Grochenig, Irregular sampling, Toeplitz matrices, and the approximation of entire functions of exponential type,

Math. Comp. 68 (1999) 749“765.

[21] K. Grochenig, Non-uniform sampling in higher dimensions: from trigonometric polynomials to band-limited functions,

in: J.J. Benedetto, P.J.S.G Ferreira (Eds.), Modern Sampling Theory: Mathematics and Applications, Birkhauser,

Boston, to appear.

[22] M. Hanke, Conjugate Gradient Type Methods for Ill-Posed Problems, Longman Scientiÿc & Technical, Harlow,

1995.

[23] M.L. Harrison, Frames and irregular sampling from a computational perspective, Ph.D. Thesis, University of

Maryland, College Park, 1998.

[24] J.R. Higgins, Sampling Theory in Fourier and Signal Analysis: Foundations, Oxford University Press, Oxford, 1996.

[25] H. Landau, Necessary density conditions for sampling and interpolation of certain entire functions, Acta Math. 117

(1967) 37“52.

[26] F.A. Marvasti, Nonuniform sampling, in: R.J. Marks II (Ed.), Advanced Topics in Shannon Sampling and

Interpolation Theory, Springer, Berlin, 1993, pp. 121“156.

[27] A.S. Nemirovski , Regularizing properties of the conjugate gradient method in ill-posed problems, Zh. Vychisl. Mat.

i Mat. Fiz. 26 (3) (1986) 332“347, 477.

316 T. Strohmer / Journal of Computational and Applied Mathematics 122 (2000) 297“316

[28] A.S. Nemirovski , B.T. Polyak, Iterative methods for solving linear ill-posed problems under precise information I,

Eng. Cybernet. 22 (1984) 1“11.