. 72
( 83 .)


for some ÿxed ¿ 1. Denote this approximation (at iteration k— ) by p1; k— . If p1; k— satisÿes the outer
stopping criterion
|p1; k (tj ) ’ bj |2 62 ||b ||2 ; (34)

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
|pN; k (tj ) ’ bj |2 62 ( ||b || + ||QN f ’ f||)||b ||; (35)

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.
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
1 (n) (n)
lim [F( k) ’ F( k )] =0 (36)
n’∞ n

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)

ak e2 ikx
f(x) = : (38)

The partial sum of the series (38) is
ak e2 ikx
fn (x) = : (39)

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
e’2 ikj=Lm=(2M +1)
zk = ; k = 0; : : : ; 2M:
2M + 1 j=’mM

(n) (n)
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
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
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


[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)
[4] G. Beylkin, On the fast Fourier transform of functions with singularities, Appl. Comp. Harm. Anal. 2 (4) (1995)
[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,
[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,
[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.


. 72
( 83 .)