A LEVINSON-GALERKIN ALGORITHM FOR REGULARIZED ...
A LEVINSON-GALERKIN ALGORITHM FOR REGULARIZED TRIGONOMETRIC APPROXIMATION∗
arXiv:math/9901123v1 [math.NA] 26 Jan 1999
THOMAS STROHMER† Abstract. Trigonometric polynomials are widely used for the approximation of a smooth function f from a set of nonuniformly spaced samples {f (xj )}N−1 j=0 . If the samples are perturbed by noise, controlling the smoothness of the trigonometric approximation becomes an essential issue to avoid overfitting and underfitting of the data. Using the polynomial degree as regularization parameter we derive a multi-level algorithm that iteratively adapts to the least squares solution of optimal smoothness. The proposed algorithm computes the solution in at most O(N M + M 2 ) operations (M being the polynomial degree of the approximation) by solving a family of nested Toeplitz systems. It is shown how the presented method can be extended to multivariate trigonometric approximation. We demonstrate the performance of the algorithm by applying it in echocardiography to the recovery of the boundary of the Left Ventricle. Key words. trigonometric approximation, Toeplitz matrix, Levinson algorithm, regularization. AMS subject classifications. 65T10, 42A10, 65D10, 65F10
1. Introduction. The necessity of recovering a function f from a finite set of −1 nonuniformly spaced measurements {f (xj )}N j=0 arises in areas as diverse as digital signal processing, geophysics, spectroscopy or medical imaging. Often f is a smooth function, in this case trigonometric polynomials may be used to recover f . Usually the measurements f (xj ) are not exactly the samples of a trigonometric polynomial and moreover the data are often distorted by several kinds of error. Hence a complete reconstruction of f from the perturbed data sj = f (xj ) + δj is not possible. In many situations a trigonometric polynomial of low degree (compared to the possibly huge number of samples) provides a good approximation to f . This trigonometric approximation may be found by solving the least squares problem min
p∈PM
N −1 X
|p(xj ) − sj |2 wj
(LSP1)
j=0
where wj > 0 are unspecified weights and PM is the space of trigonometric polynomials of degree less than or equal to M , i.e., ) ( M X 2πikx . (1.1) PM = p : p(x) = ck e k=−M
For given M with 2M + 1 ≤ N we denote the unique polynomial p ∈ PM that solves (LSP1) by p(M) . ∗ The
author has been supported by project S7001-MAT and Schr¨ odinger scholarship J01388-MAT of the Austrian Science foundation FWF. † Department of Mathematics, University of California, Davis, CA-95616, USA, [email protected] 1
2
Thomas Strohmer
Many efficient algorithms have been developed to solve (LPS1), e.g., see the articles [19, 4, 21, 9, 8]. But surprisingly little attention has been paid to the problem of how to control the smoothness of the approximation in order to avoid overfitting and underfitting of the data. An adaptation of the smoothness of the approximation can be achieved of instance by providing a suitable upper bound for the degree M of the space PM in (LSP1). In most of the aforementioned algorithms a necessary requirement to provide useful results in applications is that a good a priori guess of the degree of the trigonometric approximation is available. However a priori it is not clear what is a suitable degree for the solution, in terms of how to choose a reasonable upper bound for M when solving (LSP1). Finding a good choice for M by “trial and error” is certainly not a satisfying alternative. −1 If the data {sj }N j=0 were (i) unperturbed and (ii) stem from sampling a trigonometric polynomial (with degree less than N/2), then the solution of (LSP1) would automatically have the appropriate degree, since the original function could be completely recovered in this case. However the assumptions (i) and (ii) are rarely met in applications and controlling the smoothness of the solution becomes essential to avoid overfitting and underfitting of the data. If we choose the upper bound for the degree in (LSP1) too large, the solution will almost always take on the maximal possible degree, hence being too wiggly and picking up too much noise (overfit), see also Figure 1.1 (a)–(b). In the extreme case 2M + 1 = N, p(M) will become an interpolating polynomial, mostly with strong oscillations and far away from approximating f between the given samples. On the other hand, if we choose M too small, then p(M) will be very smooth but poorly fitting the given data (underfit), Figure 1.1(c) illustrates this behavior.
It is the goal of this paper to derive an efficient algorithm that computes the trigonometric approximation which provides the “optimal” balance between fitting the given data and preserving smoothness of the solution. Here optimality is meant in the sense that the solution has minimal degree among all trigonometric polynomials that satisfy a certain least squares criterion. The algorithm iteratively adapts to the least squares approximation of optimal degree by solving a family of nested Toeplitz systems in at most O(M N +M 2 ) operations, Compare also Figure 1.1(d), which shows the “regularized” trigonometric approximation obtained by the algorithm proposed in this paper, to which we will refer as Levinson-Galerkin algorithm. The paper is organized as follows. In Section 2 we present the main results, including the Levinson-Galerkin algorithm. Some aspects of extending the algorithm to multivariate trigonometric polynomials are discussed in 3. In Section 4 we present some applications in echocardiography. A section with miscellaneous remarks concludes the paper. We use the standard notation from numerical analysis. The inner product is P −1 denoted by h·, ·i, thus for two vectors v, w ∈ CN we have hv, wi = N j=0 vj w j , where the bar denotes complex conjugation. We denote the transpose of a matrix A by AT
3
Regularized trigonometric approximation
0.1 0.08
0.1
original signal noisy samples
0.05
0.06 0.04
0
0.02 0
−0.05
−0.02 −0.04
−0.1
approximation original signal noisy samples
−0.06 −0.15
−0.08 0
200
400
600
800
1000
0
(a) Original signal and perturbed samples
0.1 0.08 0.06
400
600
800
1000
(b) Least squares approximation using a too large upper bound for the degree of the polynomial (overfit)
0.1
approximation original signal noisy samples
approximation original signal noisy samples
0.08 0.06
0.04
0.04
0.02
0.02
0
0
−0.02
−0.02
−0.04
−0.04
−0.06
−0.06
−0.08 0
200
−0.08 200
400
600
800
1000
0
(c) Least squares approximation using a too small upper bound for the degree of the polynomial (underfit)
200
400
600
800
1000
(d) Regularized approximation by proposed Levinson-Galerkin algorithm
Fig. 1.1. Controlling the smoothness of the solution is essential for trigonometric approximation from perturbed data in order to avoid overfitting and underfitting of the data. The proposed Levinson-Galerkin algorithm automatically adapts to the least squares solution of optimal degree.
and its conjugate transpose by A∗ . The norm of p ∈ PM is given by
kpk =
Z1 0
21
|p(x)|2 dx =
M X
k=−M
|ck |2
! 21
.
(1.2)
In some applications it is advantageous to deal with complex-valued polynomials (see also Section 4), hence we do not restrict ourselves to the case of real-valued trigonometric approximation.
4
Thomas Strohmer
We will call a sequence 0 ≤ x0 < · · · < xN −1 < 1 a stable sampling set for PM if it satisfies max(xj+1 − xj ) < j
1 . 2M
(1.3)
2. A Levinson-Galerkin approach to trigonometric approximation. A standard method in numerical analysis to find the optimal balance between fitting the given data and preserving smoothness of the solution is to introduce a regularization parameter. The best value of this regularization parameter is then determined for instance by generalized cross validation [12] or via the L-curve [17]. Here we understand regularization not as a way to stabilize ill-conditioned problems, but in a broader context as a means of finding the best compromise between fitting a given set of data and preserving smoothness of the solution. As we will see, in our case it is not necessary to introduce an additional parameter, since we can regularize the smoothness of the solution by varying the parameter M of the space PM in which we are searching for the solution of (LSP1). Our approach can be seen as a Galerkin-like procedure, since we try to recover f by searching for a solution in a finite-dimensional space spanned by orthogonal polynomials and by increasing the dimension of the space we increase the resolution and add more and more details to the solution. 2.1. Multi-level least squares approximation. The only a priori knowledge we assume about f are the perturbed samples sj = f (tj ) + δj and an estimate of PN −1 the total perturbation δ ≈ j=0 |δj |2 , we do not require further information such as decay properties of the Fourier coefficients of f . We consider following multi-level scheme. For each level M = 0, 1, . . . N2−1 we solve (LSP1) yielding an approximation p(M) . Among all these approximations p(M) we select the “optimal” least squares solution p∗ by defining N −1 X δ |p(M) (xj ) − sj |2 wj ≤ PN −1 p(M) : p∗ = arg min . (2.1) 0≤M≤ N 2−1 |sj |2 j=0
j=0
In other words the regularized trigonometric approximation p∗ has minimal degree among all p(M) satisfying the criterion N −1 X j=0
δ |p(M) (xj ) − sj |2 wj ≤ PN −1 j=0
|sj |2
.
(LSP2)
Equivalently, p∗ could be characterized as follows: among all trigonometric polynomials of minimal degree that satisfy (LSP2), p∗ minimizes (LSP1). A straightforward numerical procedure to compute p∗ would be to solve (LSP1) for each level by one of the fast methods proposed e.g. in [21, 9]. Starting at level M = 0 we compute p(0) by solving (LSP1). If p(0) does not satisfy the stopping criterion (LSP2), we proceed by solving (LSP1) at the next level M = 1 and compute p(1) . We
Regularized trigonometric approximation
5
continue until the stopping criterion (LSP2) is satisfied. The overall computational effort would be O(N M 2 ) operations. In the sequel we show that p∗ can be computed in at most O(N M + M 2 ) operations by solving a family of nested Toeplitz systems. 2.2. A Levinson-Galerkin algorithm. The starting point for the derivation of our algorithm is a theorem of Gr¨ochenig [15, 9]. Theorem 2.1. Given the sampling points 0 ≤ x0 < . . . , xN −1 < 1, samples N −1 {sj }N j=1 , positive weights {wj }j=0 and the degree M with 2M + 1 ≤ N . The polynomial p(M) ∈ PM that solves (LSP1) is given by p(M) (x) =
M X
2πimx c(M) ∈ PM . m e
(2.2)
m=−M (M)
where its coefficients ck
satisfy
c(M) = [T (M) ]−1 b(M) ∈ C(2M+1)
2
(2.3)
where T (M) is a (2M + 1) × (2M + 1) Toeplitz matrix with entries (M)
Tlm =
N −1 X
wj e2πi(l−m)xj
for |l|, |m| ≤ M
(2.4)
for |m| ≤ M .
(2.5)
j=0
and b(M) = m
N −1 X
sj wj e2πimxj
j=0
If in addition {xj }N j=1 is a stable sampling set with max(xj+1 − xj ) := γ and wj = (xj+1 − xj−1 )/2, then the condition number of T (M) is bounded by κ(T
(M)
)≤
1+γ 1−γ
2
.
(2.6)
−1 Remark: T (M) is hermitian, positive semidefinite for all sampling sets {xj }N j=0 , −1 it is invertible if 2M + 1 ≤ N . If the samples {xj }N j=0 are regularly distributed with j xj = N and 2M + 1 ≤ N , then a simple computation yields that T (M) is the identity matrix for wj = N1 .
Following simple observation is crucial for the derivation of the proposed LevinsonGalerkin algorithm. Lemma 2.2. For fixed degree M and M + 1 let T (M) , b(M) and T (M+1) , b(M+1) be the Toeplitz matrices and right hand sides as defined in (2.4) and (2.5), respectively.
6
Thomas Strohmer
Then T (M) and b(M) are embedded in T (M+1) and b(M+1) in the following way:
T (M+1)
=
t0
...
t2(M+1)
.. .
T (M)
.. .
t2(M+1)
...
t0
,
b(M+1)
b−(M+1)
(M) = b bM+1
(2.7)
Proof. (2.7) follows immediately from the definition of T (M) and b(M) in (2.4) and (2.5). Unfortunately the solutions c(M) and c(M+1) of the systems T (M) c(M) = b(M) and T c = b(M+1) are not related is such a simple manner. But we can exploit the nested structure of the family {T (M) } by solving the systems T (M) c(M) = b(M) recursively via a modified Levinson algorithm. The standard Levinson algorithm cannot be applied directly, since it only addresses Toeplitz systems, where the principal leading sub-matrix and the principal leading sub-vector of the right hand side stay unchanged during the recursion, which is not the case here. For T (M+1) it does not matter, if we enlarge T (M) by appending new entries below or above, whwereas the right hand side b(M) cannot be rearranged in such a way, the principal leading subvector of the right hand side will be changed if we switch from b(M) at level M to b(M+1) at level M + 1. (M+1) (M+1)
To adapt Levinson’s algorithm to our situation, we have to split up the change from the system T (M) c(M) = b(M) at level M to the system T (M+1) c(M+1) = b(M+1) at level M + 1 into two separate steps. Instead of indexing the matrix T (M) and the vectors b(M) , c(M) by the degree M , it is therefore advantageous to index them according to their dimension. For clarity of presentation we reserve the subscript (M ) for the degree of the polynomial and its coefficient vector respectively, and use the subscript (ℓ) when we refer to the dimension of the corresponding coefficient vector in Cℓ . Thus for even ℓ, b(ℓ) = [b− ℓ +1 , . . . , b ℓ ]T ∈ Cℓ , and for odd ℓ we set 2 2 b(ℓ) = [b− ℓ−1 , . . . , b ℓ−1 ]T (whence b(1) = b0 ), analogously for c(ℓ) . Further it is useful 2
2
in the sequel to denote t(ℓ) = [t1 , . . . , tℓ ]T . Then the Toeplitz matrix T (ℓ) of size ℓ×ℓ is PN generated by the vector [t0 , (t(ℓ−1) )T ]T with tk = j=1 wj e2πikxj according to (2.4).
Assume we have already solved the system T (M) c(M) = b(M) at level M (with ℓ = 2M + 1) and now we want to switch to the next level M + 1. As we have agreed, we do this in two steps. In the first step (ℓ → ℓ + 1) the Toeplitz system can be written as "
T (ℓ) (t(ℓ) )T E (ℓ)
E (ℓ) t(ℓ) t0
#"
v (ℓ) v ℓ+1 2
#
"
b(ℓ) = b ℓ+1 2
#
,
(2.8)
7
Regularized trigonometric approximation
where E (ℓ) is the rotated identity matrix on Cℓ , 0 . .. E (ℓ) = 1
i.e., 1 . 0
System (2.8) can be solved recursively by the standard Levinson algorithm [18, 13]. To be more detailed, assume that we have already solved the system T (ℓ) c(ℓ) = b(ℓ) for ℓ = 2M + 1 and assume further that the solution of the ℓ-th order Yule-Walker system T (ℓ) y (ℓ) = −t(ℓ) is available. Then the solution of (2.8) can be computed recursively by v ℓ+1 = (b ℓ+1 − [t(ℓ) ]T E (ℓ) c(ℓ) )/βℓ 2
v
(ℓ)
2
=c
(ℓ)
+ v ℓ+1 E (ℓ) y (ℓ) 2
where βℓ = t0 + [t(ℓ) ]T y (ℓ) = (1 − αℓ−1 αℓ−1 )βℓ−1 αℓ = −(tℓ+1 + [t(ℓ) ]T E (ℓ) y (ℓ) )/βℓ z (ℓ) = y (ℓ) + αℓ E (ℓ) y (ℓ) " # z (ℓ) (ℓ+1) y = . αℓ Now we can proceed to the second step (ℓ + 1 → ℓ + 2 = 2(M + 1) + 1), where the Toeplitz system can be expressed as " #" # " # b− ℓ+1 t0 (t(ℓ+1) )∗ v− ℓ+1 2 2 = (ℓ+1) (2.9) t(ℓ+1) T (ℓ+1) v (ℓ+1) b with c(ℓ+2) = [v− ℓ+1 , (v (ℓ+1) )T ]T = c(M+1) . Observe that (2.9) cannot be transformed 2 to a system of the form (2.8) by simple permutations, i.e. just by interchanging rows and columns. Since we have already solved the systems T (ℓ+1) c(ℓ+1) = b(ℓ+1) and T (ℓ+1) y (ℓ+1) = −t(ℓ+1) we can write v (ℓ+1) = (T (ℓ+1) )−1 (b(ℓ+1) − t(ℓ+1) v− ℓ+1 ) = c(ℓ+1) + v− ℓ+1 y (ℓ+1) 2
2
and v− ℓ+1 =(b− ℓ+1 − [t(ℓ+1) ]∗ v (ℓ+1) )/t0 2
2
=(b− ℓ+1 − [t(ℓ+1) ]∗ c(ℓ+1) − [t(ℓ+1) ]∗ v− ℓ+1 y (ℓ+1) )/t0 2
2
=(b− ℓ+1 − [t(ℓ+1) ]∗ c(ℓ+1) )/βℓ+1 , 2
where we have used in the last step that T (ℓ) = [T (ℓ)]∗ which implies that (t(ℓ+1) )∗ y (ℓ+1) is real and therefore t0 + (t(ℓ+1) )∗ y (ℓ+1) = t0 + (t(ℓ+1) )T y (ℓ+1) = βℓ+1 .
8
Thomas Strohmer
Note that at each level M we have to check if the stopping criterion (LSP2) is satisfied. The evaluation of the expression N −1 X
|p(M) (xj ) − sj |2 wj
(2.10)
j=0
can be considerably simplified and by avoiding the evaluation of p(M) at the nonuniformly spaced points xj we can reduce the computational effort from O(M N ) to O(M ) operations. −1 N with the To do this we define the subspace R = {p(xj )}N j=0 : p ∈ P M ⊆ C PN −1 N weighted inner product hy, ziR = j=0 yj z¯j wj for y, z ∈ C . The solution of the −1 N least squares problem (LSP1) is the orthogonal projection of the vector {sj }N j=0 ∈ C onto R and therefore must satisfy h{p(M) (xj )} − {sj }, {p(M) (xj )}iR =
N −1 X
(p(M) (xj ) − sj )p(M) (xj )wj = 0
j=0
which implies h{p(M) (xj )}, {sj }iR = h{p(M) (xj )}, {p(M) (xj )}iR .
(2.11)
Since N −1 X
|sj − p(M) (xj )|2 wj =
N −1 X
|sj |2 wj − 2 Re hs, {p(M) (xj )}iR +
=
N −1 X
|p(M) (xj )|2 wj
j=0
j=0
j=0
N −1 X
|sj |2 wj −
N −1 X
|p(M) (xj )|2 wj
j=0
j=0
by (2.11), and because N −1 X
|p
(M)
2
(xj )| wj =
N −1 X
wj
j=0
j=0
=
M X
X M
2πimxj c(M) m e
m=−M
M X
= hT
c
,c
(M) cn e2πinxj
n=−M
(M)
(M) cm cn
NX −1
wj e2πi(m−n)xj
j=0
m=−M n=−M (M) (M)
X M
(M)
i = hb
(M)
, c(M) i ,
(2.12)
it follows that N −1 X j=0
|sj − p(M) (xj )|2 wj =
N −1 X
|sj |2 wj − hb(M) , c(M) i .
(2.13)
j=0
PN −1 Since j=0 |sj |2 wj has to be computed only once at the beginning of the algorithm, the evaluation of (2.10) can be carried out in O(M ) operations. Summing up we have arrived at following algorithm to compute p∗ .
Regularized trigonometric approximation
9
Algorithm 1 (Levinson-Galerkin algorithm for trigonometric polynomials). Let −1 N −1 the sampling points {xj }N j=0 , sampling values {sj }j=0 , weights wj > 0 and the perturbation estimate δ be given. Then the trigonometric polynomial p∗ that solves (2.1) can be computed in O(N M + M 2 ) operations by following algorithm. PN −1 PN −1 PN −1 PN −1 Initialize: t0 = j=0 wj , t1 = j=0 wj e2πixj , b0 = j=0 sj wj , s = j=0 |sj |2 wj , y (1) = −t1 /t0 , c(1) = b0 /t0 , β0 = t0 , α0 = −t1 /t0 , ε1 = s − b20 /t0 , ℓ = 1.
while εℓ > δ βℓ = (1 − αℓ−1 αℓ−1 )βℓ−1 if ℓ ≡ 1 mod 2 N −1 X sj wj eπi(ℓ+1)xj b ℓ+1 = 2
j=0
b ℓ+1 − hE (ℓ) c(ℓ) , t(ℓ) i
v ℓ+1 =
2
βℓ
2
v
(ℓ)
+ v ℓ+1 E (ℓ) y (ℓ) #2 v (ℓ) = v ℓ+1 " 2 # b(ℓ) = b ℓ+1
=c
c(ℓ+1) b(ℓ+1)
(ℓ)
"
2
elseif ℓ ≡ 0 mod 2 N −1 X b− ℓ = sj wj e−πiℓxj 2
j=0
v− ℓ =
b− ℓ − hc(ℓ) , t(ℓ) i 2
βℓ
2
v
(ℓ)
=c
(ℓ)
c(ℓ+1) = b
(ℓ+1)
=
+ v− ℓ y (ℓ) " #2 v− ℓ 2
v (ℓ)
"
b− ℓ
2
b(ℓ)
#
εℓ+1 = s − hb(ℓ+1) , c(ℓ+1) i end t(ℓ+1) =
N −1 X
wj e2πi(ℓ+1)xj
j=0
αℓ = −
t(ℓ+1) + hE (ℓ) y (ℓ) , t(ℓ) i βℓ
10
Thomas Strohmer
z (ℓ) = y (ℓ) + αℓ E (ℓ) y (ℓ) " # z (ℓ) (ℓ+1) y = αℓ " # t(ℓ) (ℓ+1) t = tℓ+1 ℓ=ℓ+1 end M = ℓ/2 p∗ (t) =
M X
2πimt c(M) m e
m=−M
Remark: Usually one evaluates the final approximation on regularly spaced grid points, hence the last step of the algorithm can be realized by a Fast Fourier transform. The most costly steps are the computation of the entries of t(ℓ) and b(ℓ) . According to Corollary 1 in [9] the entries of T (M) and b(M) can also be computed via FFT by embedding the xj into a regular grid (since the xj can be stored only in finite precision). In this case one automatically gets all entries t0 , . . . , tN −1 at once. However this trick is only useful if the number of points of the regular grid is of the same magnitude as the number of sampling points. Alternatively one may use the numerical attractive formulas of Rokhlin and Dutt [6, 7] or Beylkin [2] for fast evaluation of trigonometric sums at unequally spaced nodes. Algorithm 1 can be simplified for real-valued data, this modification is left to the reader. 3. Multivariate trigonometric approximation. An advantage of the proposed approach, besides its numerical efficiency, is the fact that it can be easily extended to multivariate trigonometric approximation. In this section we briefly discuss some results for the two-dimensional case. 2 We define the space of 2-D trigonometric polynomials PM by M X 2 PM = p : p(x, y) = cj,k e2πi(jx+ky) . (3.1) j,k=−M
To reduce the notational burden, we have assumed in (3.1) that p has equal degree M in each coordinate, the extension to polynomials with different degree in each coordinate is straightforward. −1 2 For an arbitrary sampling set {(xj , yj )}N j=0 ∈ [0, 1) and given degree M the system matrix according to the 2-D version of Theorem 2.1 is [24] (T
(M)
)k,l =
N −1 X j=0
wj e2πi(k−l)(xj +yj ) ,
k, l = 0, . . . 2M .
(3.2)
11
Regularized trigonometric approximation
One can easily verify that T (M) is a block Toeplitz matrix with Toeplitz blocks [24]. In [1] Levinson’s algorithm has been extended to general block Toeplitz systems. With these results at hand, we can easily generalize Algorithm 1 to multivariate trigonometric approximation. Clearly the question of invertibility of T (M) and estimates of the condition number become more difficult. The condition (2M + 1)d ≤ N is necessary in dimension d > 1, but not sufficient any more, since the fundamental theorem of algebra does not hold in dimensions larger than one. In [16] Gr¨ochenig has derived estimates for the condition number of T (M) in higher dimensions. However, these estimates become worse with increasing dimension. Further work remains to be done in that direction. In the following we consider a special case of trigonometric approximation in two dimensions. This case arises when a function is irregularly sampled along lines. A typical example is illustrated in Figure 3.1. Such sampling patterns are encountered for instance in geophysics and medical imaging, see also Section 4.2.
1 0 0 1
1 0 0 1 1 0 1 0
1 0 1 1 0 0 0 1
1 0 1 0
1 0 0 1
1 0 1 0 1 0 0 1
1 0 1 0 1 0 0 1
1 0 0 1
1 0 1 0
1 0 0 1 1 0 1 0 1 0 0 1 1 0 1 0
1 0 1 0
1 0 0 1
1 0 1 0
1 0 0 1
1 0 0 1
1 0 0 1
1 0 1 0
1 0 0 1
1 0 1 0
1 0 0 1
{
0 1 1 0 1 0 0 1
1 0 0 1 1 0 0 1
0 1 1 0 1 0 0 1 0 1 0 1
1 0 0 1
1 0 1 0
1 0 1 0
1 0 0 1 1 0 0 1 1 0 0 1
1 0 0 1
1 0 0 1 1 0 0 1
1 0 0 1 1 0 0 1 11 00 11 00
11 00 11 1 00 0 0 1 1 0 1 0
11 00 00 11
1 0 0 1 1 0 0 1 1 0 1 0 0 1 0 1
δ2
1 0 0 1
1 0 1 0 1 0 0 1 1 0 1 0
1 0 0 1
1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1
1 0 0 1 1 0 0 1
{
y
δ1
1 0 0 1
x
Fig. 3.1. Line-type nonuniform sampling set 2 Corollary 3.1. Let p ∈ PM and let {xj , yjk }, j = 0, . . . , N −1, k = 0, . . . , Nj −1 2 be a sampling set in [0, 1) such that Nj −1
A1 kpk2 ≤
X
|p(yjk )|2 ≤ B1 kpk2
A1 , B1 > 0
(3.3)
k=0
for every p ∈ PM and for all j. Further assume that {xj }j∈Z is a sampling set such that 2
A2 kpk ≤
N −1 X j=0
|p(xj )|2 ≤ B2 kpk2
A2 , B2 > 0
(3.4)
12
Thomas Strohmer
for every p ∈ PM . Then 2
A1 A2 kpk ≤
j −1 N −1 NX X
|p(xj , yjk )|2 ≤ B1 B2 kpk2
(3.5)
j=0 k=0
2 for every p ∈ PM . 1 If {xj } and {yjk } are stable sampling sets with supk (xk+1 − xk ) = δ2 < 2M and 1 2 2 supj,k (yj,k+1 − yjk ) = δ1 < 2M , then Al = (1 − δl ) , Bl = (1 + δl ) , l = 1, 2 and the condition number of the block Toeplitz matrix T (M) is bounded by
κ(T ) ≤
(1 + δ1 )2 (1 + δ2 )2 . (1 − δ1 )2 (1 − δ2 )2
(3.6)
Proof. Let x be fixed. Then y → p(x, y) ∈ PM and hence for all j
A1
Z1
Nj −1 2
|p(xj , y)| dy ≤
X
2
|p(xj , yjk )| ≤ B1
k=0
0
Z1
|p(xj , y)|2 dy
(3.7)
0
by assumption (3.3). It follows that
A1
N −1 Z1 X
2
|p(xj , y)| dy ≤
j=0 0
X
2
|p(xj , yjk )| ≤ B1
N −1 Z1 X
|p(xj , y)|2 dy
(3.8)
j=0 0
j,k
Now let y be fixed. Then x → p(x, y) ∈ PM and
A2
Z1
2
|p(x, y)| dx ≤
N −1 X
2
|p(xj , y)| ≤ B2
j=0
0
Z1
|p(x, y)|2 dx .
(3.9)
0
Since N −1 Z1 X j=0 0
2
|p(xj , y)| dy =
Z1 NX −1 0
|p(xj , y)|2 dy ,
(3.10)
j=0
assertion (3.5) follows by combining (3.8) and (3.9) with (3.10). The estimate of the constants Al , Bl and of the condition number of the block Toeplitz matrix T (M) follow from Theorem 2.1. The proof of Corollary 3.1 is due to Gr¨ochenig [14]. Corollary 3.1 does not only 2 guarantee that p ∈ PM can be recovered from its samples p(xj , yjk ), it provides more. An immediate consequence is, that f can be recovered by an efficient algorithm relying on an successive application of Algorithm 1 and the Gohberg-Semencul representation of the inverse of a Toeplitz matrix. See Section 4.2 for more details and an application in medical imaging.
Regularized trigonometric approximation
13
4. Curve and surface approximation by trigonometric polynomials. Trigonometric polynomials can be used to model the boundary or the surface of smooth objects. Let us consider a two-dimensional object, obtained e.g. by a planar cross-section from a 3-D object and assume that the boundary of this 2-D object is a closed curve in R2 . We denote this curve by f and parametrize it by f (u) = (xu , yu ), where xu and yu are the coordinates of f at “time” u in the x- and y-direction respectively. Obviously we can interpret f as a one-dimensional continuous, complex, and periodic function, where xu represents the real part and yu represents the imaginary part of f (u). It follows from the Theorem of Weierstrass (and from the Theorem of Stone-Weierstrass [22] for higher dimensions) that a continuous periodic function can be approximated uniformly by trigonometric polynomials. If f is smooth, we can fairly assume that trigonometric polynomials of low degree provide an approximation of sufficient precision. Assume that we know only some arbitrary, perturbed points sj = (xuj , yuj ) = f (uj ) + δj , j = 0, . . . , N − 1 of f , and we want to recover f from these points. By a slight abuse of notation we interpret sj as complex number and write sj = xj + iyj .
(4.1)
We relate the curve parameter u to the boundary points sj by computing the distance between two successive points sj−1 , sj via u0 = 0
(4.2)
uj = uj−1 + dj q dj = (xj − xj−1 )2 + (yj − yj−1 )2
(4.3) (4.4)
for j = 1, . . . , N − 1. Via the normalization tj = tj /L with L = uN −1 + dN we force all sampling points to be in [0, 1). Other choices for dj in (4.1) can be found in [5] in conjunction with curve approximation using splines. Having carried out the transformations (4.1)–(4.4), we can solve the problem of recovering the curve f from its perturbed points sj by Algorithm 1. 4.1. Object boundary recovery in Echocardiography. Trigonometric polynomials are certainly not suitable to model the shape of arbitrary objects. However they are often useful in cases where an underlying (stationary) physical process implies smoothness conditions of the object. Typical examples arise in medical imaging, for instance in clinical cardiac studies, where the evaluation of cardiac function using parameters of left ventricular contractibility is an important constituent of an echocardiographic examination [26]. These parameters are derived using boundary tracing of endocardial borders of the Left Ventricle (LV). The extraction of the boundary of the LV comprises two steps, once the ultrasound image of a cross section of the LV is given, see Figure 4.1(a)–(d). First an edge detection is applied to the ultrasound image to detect the boundary of the LV, cf. Figure 4.1(c). However this procedure
14
Thomas Strohmer
may be hampered by the presence of interfering biological structures (such as papillar muscles), the uneveness of boundary contrast, and various kinds of noise [25]. Thus edge detection often provides only a set of nonuniformly spaced, perturbed boundary points rather than a connected boundary. Therefore a second step is required, to recover the original boundary from the detected edge points, cf. Figure 4.1(d). Since the shape of the Left Ventricle is definitely smooth, trigonometric polynomials are particularly well suited to model its boundary. After having transformed the detected boundary points as described in (4.1)–(4.4) we can use Algorithm 1 to recover the boundary. The noise level δ depends on the technical equipment under use, it can be determined from experimental experience. Figure 4.2(a)–(b) demonstrate the importance of determining a proper degree for the approximating polynomial. The approximation displayed in Figure 4.2(a) has been computed by solving (LSP1) where M has been chosen too small, we obviously have underfitted the data. The overfitted approximation obtained by solving (LSP1) using a too large M is shown in Figure 4.2(b). The approximation shown in Figure 4.1(d) has been computed by Algorithm 1, it provides the optimal balance between fitting the data and smoothness of the solution. 4.2. Boundary recovery from a sequence of images. In cardiac clinical studies one is more interested in the behavior of the Left Ventricle over a period of time rather than in a single “snapshot”. Thus for a fixed cross section we are given a sequence of ultrasound images (usually regularly spaced in time) describing the variation of the shape of the LV with time. One cycle from diastole (the state of maximal contraction of the LV), passing systole (the state of maximal expansion) to the next diastole consists typically of about 30 image frames. Since the behavior of the LV is (at least for a short period of time) almost periodic, one can model the varying shape of a fixed cross section of the LV as a distorted two-dimensional torus, which in turn can be interpreted as two-dimensional trigonometric polynomial. Clearly we have to use a different degree for the time coordinate τ and for the spatial coordinate u. Due to interfering biological structures and other distortions it sometimes happens that some of the image frames cannot be used to extract any reliable boundary information. Thus we have to approximate these missing boundaries from the information of the other image frames. To be more precise, assume that an echocardiographic examination provides a sequence of ultrasound images Iτ taken at time points τ = 0, 1, . . . , T − 1, where T is approximately the length of one diastolic cycle (the time points could also be nonuniformly spaced). Assume that some of the images Iτ Nj −1 provide no useful information, so that we can only detect boundary points {sj,k }k=0 −1 from the images Iτj , where {τj }N j=0 is a subset of 0, 1, . . . T − 1. In order to get a complete description of the LV for the time interval [0, T ], we have not only to approximate the boundaries fj from each Ij , but we also have to recover the boundaries corresponding to the missing images. In other words we look for a 2-D trigonometric
Regularized trigonometric approximation
(a) 2-D echocardiography
(b) Cross section of Left Ventricle
(c) Detected boundary points
(d) Recovered boundary of LV computed by Algorithm 1
15
Fig. 4.1. The recovery of the boundary of the Left Ventricle from 2-D ultrasound images is a basic step in echocardiography to extract relevant parameters of cardiac function.
2 polynomial p∗ ∈ PM of appropriate degree M that satisfies p(τj , uj,k ) ≡ (xj,k , yj,k ) where the parameter u is related to sj,k = xj,k + iyj,k by formulas (4.2)–(4.4). This approximation can be computed by the 2-D version of Algorithm 1, as indicated in the beginning of Section 3.
Under certain conditions we can use the 1-D version of Algorithm 1 instead of its 2-D version. As long as the assumptions of Corollary 3.1 are satisfied, we can
16
Thomas Strohmer
(a) Underfitted solution
(b) Overfitted solution
Fig. 4.2. 2 compute p∗ ∈ PM by a successive application of Algorithm 1. We first approximate Nj −1 the boundaries fj for each j separately from its samples {sj,k }k=0 , which yields j (Mj ) different polynomials p ∈ PMj . Having done this, the next step is to recover the missing boundaries at those time points where no information is available. We proceed by approximating successively the missing information “line by line”. We choose u = 0, say, and approximate the missing information from the samples p(Mj ) (u) taken at the time points τj , j = 0, . . . , N − 1. (M) (M) (M) Note that the Toeplitz matrices of the systems Tu cu = bu coincide for all u, since the sampling geometry is constant along the u-coordinate (because we have recovered all samples at each τj ). Thus we have to solve multiple Toeplitz systems with the same system matrix but different right hand side. It is well-known that this can be done efficiently by exploiting the Gohberg-Semencul representation of the inverse of the Toeplitz matrix [11]. In our context this reads as follows. We solve
= b(M) Tu(M) c(M) u u
(4.5)
for one u by Algorithm 1. We can solve now all other systems efficiently by establishing (T (M) )−1 in the Gohberg-Semencul form (4.6) (T (M) )−1 = [L(M) ]∗ L(M) − U (M) [U (M) ]∗ /z0
where L(M) is a lower triangular Toeplitz matrix with z = [z0 , z1 , . . . , z2M ]T as its first column, U (M) is an upper triangular Toeplitz matrix with [0, z1 , . . . , z2M ]T as its last column, z being the first column of (T (M) )−1 . The matrix vector multiplications to (M) (M) compute cu = (T (M) )−1 can now be carried out quickly using the Fast Fourier u bu (M) transform by embedding L and U (M) into circulant matrices.
Regularized trigonometric approximation
17
5. Miscellaneous remarks. For sampling sets with large gaps it can happen that the system T (M) c(M) = b(M) gets ill-conditioned with increasing degree M and therefore Algorithm 1 may become instable [3]. In this case one can use a different, more robust approach, which however comes at higher computational costs [23]. We solve the system T (M) c(M) = b(M) iteratively, e.g. by the conjugate gradient method until a certain stopping criterion is satisfied at iteration k, say, yielding the solution (M) (M+1) ck . We use this solution as initial guess at the next level M + 1 by setting c0 = (M) [0 (ck )T 0]T . The crucial point of such a procedure is to find a suitable stopping criterion at each level to ensure convergence of the iterates, see [23, 20] for more details. The computation of the entries of the Toeplitz matrix in Section 4 involves the nodes uj which in this particular case depend on the (perturbed) samples sj . Therefore not only the right hand side b(M) , but also T (M) is subject to perturbations. Hence in principle one might use the concept of total least squares (see [10]) instead of the standard least squares criterion in (LSP1). A detailed discussion of this modification is beyond the scope of this paper. REFERENCES [1] H. Akaike, Block Toeplitz matrix inversion, SIAM J. Appl. Math., 24 (1973), pp. 234–241. [2] G. Beylkin, On the fast Fourier transform of functions with singularities, Appl. Comp. Harm. Anal., 2 (1995), pp. 363–381. [3] G. Cybenko, The numerical stability of the levinson-durbin algorithm for toeplitz systems of equations, SIAM J. Sci. Statist. Comp., 1 (1980), pp. 303–3190. [4] C. Demeure, Fast QR factorization of Vandermonde matrices, Linear Algebra Appl., 122–124 (1989), pp. 165–194. [5] P. Dierckx, Curve and surface fitting with splines, Monographs on Numerical Analysis, Oxford University Press, 1993. [6] A. Dutt and V. Rokhlin, Fast Fourier transforms for nonequispaced data, SIAM J. Sci. Comp., 14 (1993), pp. 1368–1394. , Fast Fourier transforms for nonequispaced data II, Appl. Comp. Harm. Anal., 2 (1995), [7] pp. 85–10. [8] H. Fassbender, On numerical methods for discrete least-squares approximation by trigonometric polynomials, Math. Comp., 66 (1997), pp. 719–741. ¨ chenig, and T. Strohmer, Efficient numerical methods in non[9] H. G. Feichtinger, K. Gro uniform sampling theory, Numerische Mathematik, 69 (1995), pp. 423–440. [10] R. Fierro, G. Golub, P. Hansen, and D. O’Leary, Regularization by truncated total least squares, SIAM J. Sci. Comp., 18 (1997), pp. 1223–1241. [11] I. Gohberg and A. Semencul, On the inversion of finite Toeplitz matrices and their continuous analogs, Mat. Issled., 2 (1972), pp. 201–233. [12] G. Golub, M. Heath, and G. Wahba, Generalized cross-validation as a method for choosing a good ridge parameter, Technometrics, 21 (1979), pp. 215–223. [13] G. Golub and C. van Loan, Matrix Computations, third ed., Johns Hopkins, London, Baltimore, 1996. ¨ chenig, personal communication. [14] K. Gro [15] , A discrete theory of irregular sampling, Linear Algebra Appl., 193 (1993), pp. 129–150. [16] , Finite and infinite-dimensional models for non-uniform sampling, in SampTA - Sam-
18
Thomas Strohmer
pling Theory and Applications, Aveiro, Portugal, 1997, pp. 285–290. [17] P. Hansen, Analysis of discrete ill-posed problems by means of the L-curve, SIAM Review, 34 (1992), pp. 561–680. [18] N. Levinson, The Wiener RMS (root-mean square) error criterion in filter design and prediction, J. Math. Phys., 25 (1947), pp. 261–278. [19] A. Newbery, Trigonometric interpolation and curve-fitting, Math. Comp., 24 (1970), pp. 869– 876. [20] M. Rauth and T. Strohmer, Smooth approximation of potential fields from noisy scattered data, Geophysics, (to appear). [21] L. Reichel, G. Ammar, and W. Gragg, Discrete least squares approximation by trigonometric polynomials, Math. Comp., 57 (1991), pp. 273–289. [22] W. Rudin, Fourier Analysis on Groups, Wiley Interscience, New York, 1976. [23] O. Scherzer and T. Strohmer, A multi–level algorithm for the solution of moment problems. preprint, 1997. [24] T. Strohmer, Computationally attractive reconstruction of band-limited images from irregular samples, IEEE Trans. Image Proc., 6 (1997), pp. 540–548. [25] M. Suessner, M. Budil, T. Strohmer, M. Greher, G. Porenta, and T. Binder, Contour detection using artifical neuronal network presegmention, in Proc. Computers in Cardiology, Vienna, 1995. [26] D. Wilson, E. Geiser, and J. Li, Feature extraction in 2-dimensional short-axis echocardiographic images, J. Math. Imag. Vision, 3 (1993), pp. 285–298.
arXiv:math/9901123v1 [math.NA] 26 Jan 1999
THOMAS STROHMER† Abstract. Trigonometric polynomials are widely used for the approximation of a smooth function f from a set of nonuniformly spaced samples {f (xj )}N−1 j=0 . If the samples are perturbed by noise, controlling the smoothness of the trigonometric approximation becomes an essential issue to avoid overfitting and underfitting of the data. Using the polynomial degree as regularization parameter we derive a multi-level algorithm that iteratively adapts to the least squares solution of optimal smoothness. The proposed algorithm computes the solution in at most O(N M + M 2 ) operations (M being the polynomial degree of the approximation) by solving a family of nested Toeplitz systems. It is shown how the presented method can be extended to multivariate trigonometric approximation. We demonstrate the performance of the algorithm by applying it in echocardiography to the recovery of the boundary of the Left Ventricle. Key words. trigonometric approximation, Toeplitz matrix, Levinson algorithm, regularization. AMS subject classifications. 65T10, 42A10, 65D10, 65F10
1. Introduction. The necessity of recovering a function f from a finite set of −1 nonuniformly spaced measurements {f (xj )}N j=0 arises in areas as diverse as digital signal processing, geophysics, spectroscopy or medical imaging. Often f is a smooth function, in this case trigonometric polynomials may be used to recover f . Usually the measurements f (xj ) are not exactly the samples of a trigonometric polynomial and moreover the data are often distorted by several kinds of error. Hence a complete reconstruction of f from the perturbed data sj = f (xj ) + δj is not possible. In many situations a trigonometric polynomial of low degree (compared to the possibly huge number of samples) provides a good approximation to f . This trigonometric approximation may be found by solving the least squares problem min
p∈PM
N −1 X
|p(xj ) − sj |2 wj
(LSP1)
j=0
where wj > 0 are unspecified weights and PM is the space of trigonometric polynomials of degree less than or equal to M , i.e., ) ( M X 2πikx . (1.1) PM = p : p(x) = ck e k=−M
For given M with 2M + 1 ≤ N we denote the unique polynomial p ∈ PM that solves (LSP1) by p(M) . ∗ The
author has been supported by project S7001-MAT and Schr¨ odinger scholarship J01388-MAT of the Austrian Science foundation FWF. † Department of Mathematics, University of California, Davis, CA-95616, USA, [email protected] 1
2
Thomas Strohmer
Many efficient algorithms have been developed to solve (LPS1), e.g., see the articles [19, 4, 21, 9, 8]. But surprisingly little attention has been paid to the problem of how to control the smoothness of the approximation in order to avoid overfitting and underfitting of the data. An adaptation of the smoothness of the approximation can be achieved of instance by providing a suitable upper bound for the degree M of the space PM in (LSP1). In most of the aforementioned algorithms a necessary requirement to provide useful results in applications is that a good a priori guess of the degree of the trigonometric approximation is available. However a priori it is not clear what is a suitable degree for the solution, in terms of how to choose a reasonable upper bound for M when solving (LSP1). Finding a good choice for M by “trial and error” is certainly not a satisfying alternative. −1 If the data {sj }N j=0 were (i) unperturbed and (ii) stem from sampling a trigonometric polynomial (with degree less than N/2), then the solution of (LSP1) would automatically have the appropriate degree, since the original function could be completely recovered in this case. However the assumptions (i) and (ii) are rarely met in applications and controlling the smoothness of the solution becomes essential to avoid overfitting and underfitting of the data. If we choose the upper bound for the degree in (LSP1) too large, the solution will almost always take on the maximal possible degree, hence being too wiggly and picking up too much noise (overfit), see also Figure 1.1 (a)–(b). In the extreme case 2M + 1 = N, p(M) will become an interpolating polynomial, mostly with strong oscillations and far away from approximating f between the given samples. On the other hand, if we choose M too small, then p(M) will be very smooth but poorly fitting the given data (underfit), Figure 1.1(c) illustrates this behavior.
It is the goal of this paper to derive an efficient algorithm that computes the trigonometric approximation which provides the “optimal” balance between fitting the given data and preserving smoothness of the solution. Here optimality is meant in the sense that the solution has minimal degree among all trigonometric polynomials that satisfy a certain least squares criterion. The algorithm iteratively adapts to the least squares approximation of optimal degree by solving a family of nested Toeplitz systems in at most O(M N +M 2 ) operations, Compare also Figure 1.1(d), which shows the “regularized” trigonometric approximation obtained by the algorithm proposed in this paper, to which we will refer as Levinson-Galerkin algorithm. The paper is organized as follows. In Section 2 we present the main results, including the Levinson-Galerkin algorithm. Some aspects of extending the algorithm to multivariate trigonometric polynomials are discussed in 3. In Section 4 we present some applications in echocardiography. A section with miscellaneous remarks concludes the paper. We use the standard notation from numerical analysis. The inner product is P −1 denoted by h·, ·i, thus for two vectors v, w ∈ CN we have hv, wi = N j=0 vj w j , where the bar denotes complex conjugation. We denote the transpose of a matrix A by AT
3
Regularized trigonometric approximation
0.1 0.08
0.1
original signal noisy samples
0.05
0.06 0.04
0
0.02 0
−0.05
−0.02 −0.04
−0.1
approximation original signal noisy samples
−0.06 −0.15
−0.08 0
200
400
600
800
1000
0
(a) Original signal and perturbed samples
0.1 0.08 0.06
400
600
800
1000
(b) Least squares approximation using a too large upper bound for the degree of the polynomial (overfit)
0.1
approximation original signal noisy samples
approximation original signal noisy samples
0.08 0.06
0.04
0.04
0.02
0.02
0
0
−0.02
−0.02
−0.04
−0.04
−0.06
−0.06
−0.08 0
200
−0.08 200
400
600
800
1000
0
(c) Least squares approximation using a too small upper bound for the degree of the polynomial (underfit)
200
400
600
800
1000
(d) Regularized approximation by proposed Levinson-Galerkin algorithm
Fig. 1.1. Controlling the smoothness of the solution is essential for trigonometric approximation from perturbed data in order to avoid overfitting and underfitting of the data. The proposed Levinson-Galerkin algorithm automatically adapts to the least squares solution of optimal degree.
and its conjugate transpose by A∗ . The norm of p ∈ PM is given by
kpk =
Z1 0
21
|p(x)|2 dx =
M X
k=−M
|ck |2
! 21
.
(1.2)
In some applications it is advantageous to deal with complex-valued polynomials (see also Section 4), hence we do not restrict ourselves to the case of real-valued trigonometric approximation.
4
Thomas Strohmer
We will call a sequence 0 ≤ x0 < · · · < xN −1 < 1 a stable sampling set for PM if it satisfies max(xj+1 − xj ) < j
1 . 2M
(1.3)
2. A Levinson-Galerkin approach to trigonometric approximation. A standard method in numerical analysis to find the optimal balance between fitting the given data and preserving smoothness of the solution is to introduce a regularization parameter. The best value of this regularization parameter is then determined for instance by generalized cross validation [12] or via the L-curve [17]. Here we understand regularization not as a way to stabilize ill-conditioned problems, but in a broader context as a means of finding the best compromise between fitting a given set of data and preserving smoothness of the solution. As we will see, in our case it is not necessary to introduce an additional parameter, since we can regularize the smoothness of the solution by varying the parameter M of the space PM in which we are searching for the solution of (LSP1). Our approach can be seen as a Galerkin-like procedure, since we try to recover f by searching for a solution in a finite-dimensional space spanned by orthogonal polynomials and by increasing the dimension of the space we increase the resolution and add more and more details to the solution. 2.1. Multi-level least squares approximation. The only a priori knowledge we assume about f are the perturbed samples sj = f (tj ) + δj and an estimate of PN −1 the total perturbation δ ≈ j=0 |δj |2 , we do not require further information such as decay properties of the Fourier coefficients of f . We consider following multi-level scheme. For each level M = 0, 1, . . . N2−1 we solve (LSP1) yielding an approximation p(M) . Among all these approximations p(M) we select the “optimal” least squares solution p∗ by defining N −1 X δ |p(M) (xj ) − sj |2 wj ≤ PN −1 p(M) : p∗ = arg min . (2.1) 0≤M≤ N 2−1 |sj |2 j=0
j=0
In other words the regularized trigonometric approximation p∗ has minimal degree among all p(M) satisfying the criterion N −1 X j=0
δ |p(M) (xj ) − sj |2 wj ≤ PN −1 j=0
|sj |2
.
(LSP2)
Equivalently, p∗ could be characterized as follows: among all trigonometric polynomials of minimal degree that satisfy (LSP2), p∗ minimizes (LSP1). A straightforward numerical procedure to compute p∗ would be to solve (LSP1) for each level by one of the fast methods proposed e.g. in [21, 9]. Starting at level M = 0 we compute p(0) by solving (LSP1). If p(0) does not satisfy the stopping criterion (LSP2), we proceed by solving (LSP1) at the next level M = 1 and compute p(1) . We
Regularized trigonometric approximation
5
continue until the stopping criterion (LSP2) is satisfied. The overall computational effort would be O(N M 2 ) operations. In the sequel we show that p∗ can be computed in at most O(N M + M 2 ) operations by solving a family of nested Toeplitz systems. 2.2. A Levinson-Galerkin algorithm. The starting point for the derivation of our algorithm is a theorem of Gr¨ochenig [15, 9]. Theorem 2.1. Given the sampling points 0 ≤ x0 < . . . , xN −1 < 1, samples N −1 {sj }N j=1 , positive weights {wj }j=0 and the degree M with 2M + 1 ≤ N . The polynomial p(M) ∈ PM that solves (LSP1) is given by p(M) (x) =
M X
2πimx c(M) ∈ PM . m e
(2.2)
m=−M (M)
where its coefficients ck
satisfy
c(M) = [T (M) ]−1 b(M) ∈ C(2M+1)
2
(2.3)
where T (M) is a (2M + 1) × (2M + 1) Toeplitz matrix with entries (M)
Tlm =
N −1 X
wj e2πi(l−m)xj
for |l|, |m| ≤ M
(2.4)
for |m| ≤ M .
(2.5)
j=0
and b(M) = m
N −1 X
sj wj e2πimxj
j=0
If in addition {xj }N j=1 is a stable sampling set with max(xj+1 − xj ) := γ and wj = (xj+1 − xj−1 )/2, then the condition number of T (M) is bounded by κ(T
(M)
)≤
1+γ 1−γ
2
.
(2.6)
−1 Remark: T (M) is hermitian, positive semidefinite for all sampling sets {xj }N j=0 , −1 it is invertible if 2M + 1 ≤ N . If the samples {xj }N j=0 are regularly distributed with j xj = N and 2M + 1 ≤ N , then a simple computation yields that T (M) is the identity matrix for wj = N1 .
Following simple observation is crucial for the derivation of the proposed LevinsonGalerkin algorithm. Lemma 2.2. For fixed degree M and M + 1 let T (M) , b(M) and T (M+1) , b(M+1) be the Toeplitz matrices and right hand sides as defined in (2.4) and (2.5), respectively.
6
Thomas Strohmer
Then T (M) and b(M) are embedded in T (M+1) and b(M+1) in the following way:
T (M+1)
=
t0
...
t2(M+1)
.. .
T (M)
.. .
t2(M+1)
...
t0
,
b(M+1)
b−(M+1)
(M) = b bM+1
(2.7)
Proof. (2.7) follows immediately from the definition of T (M) and b(M) in (2.4) and (2.5). Unfortunately the solutions c(M) and c(M+1) of the systems T (M) c(M) = b(M) and T c = b(M+1) are not related is such a simple manner. But we can exploit the nested structure of the family {T (M) } by solving the systems T (M) c(M) = b(M) recursively via a modified Levinson algorithm. The standard Levinson algorithm cannot be applied directly, since it only addresses Toeplitz systems, where the principal leading sub-matrix and the principal leading sub-vector of the right hand side stay unchanged during the recursion, which is not the case here. For T (M+1) it does not matter, if we enlarge T (M) by appending new entries below or above, whwereas the right hand side b(M) cannot be rearranged in such a way, the principal leading subvector of the right hand side will be changed if we switch from b(M) at level M to b(M+1) at level M + 1. (M+1) (M+1)
To adapt Levinson’s algorithm to our situation, we have to split up the change from the system T (M) c(M) = b(M) at level M to the system T (M+1) c(M+1) = b(M+1) at level M + 1 into two separate steps. Instead of indexing the matrix T (M) and the vectors b(M) , c(M) by the degree M , it is therefore advantageous to index them according to their dimension. For clarity of presentation we reserve the subscript (M ) for the degree of the polynomial and its coefficient vector respectively, and use the subscript (ℓ) when we refer to the dimension of the corresponding coefficient vector in Cℓ . Thus for even ℓ, b(ℓ) = [b− ℓ +1 , . . . , b ℓ ]T ∈ Cℓ , and for odd ℓ we set 2 2 b(ℓ) = [b− ℓ−1 , . . . , b ℓ−1 ]T (whence b(1) = b0 ), analogously for c(ℓ) . Further it is useful 2
2
in the sequel to denote t(ℓ) = [t1 , . . . , tℓ ]T . Then the Toeplitz matrix T (ℓ) of size ℓ×ℓ is PN generated by the vector [t0 , (t(ℓ−1) )T ]T with tk = j=1 wj e2πikxj according to (2.4).
Assume we have already solved the system T (M) c(M) = b(M) at level M (with ℓ = 2M + 1) and now we want to switch to the next level M + 1. As we have agreed, we do this in two steps. In the first step (ℓ → ℓ + 1) the Toeplitz system can be written as "
T (ℓ) (t(ℓ) )T E (ℓ)
E (ℓ) t(ℓ) t0
#"
v (ℓ) v ℓ+1 2
#
"
b(ℓ) = b ℓ+1 2
#
,
(2.8)
7
Regularized trigonometric approximation
where E (ℓ) is the rotated identity matrix on Cℓ , 0 . .. E (ℓ) = 1
i.e., 1 . 0
System (2.8) can be solved recursively by the standard Levinson algorithm [18, 13]. To be more detailed, assume that we have already solved the system T (ℓ) c(ℓ) = b(ℓ) for ℓ = 2M + 1 and assume further that the solution of the ℓ-th order Yule-Walker system T (ℓ) y (ℓ) = −t(ℓ) is available. Then the solution of (2.8) can be computed recursively by v ℓ+1 = (b ℓ+1 − [t(ℓ) ]T E (ℓ) c(ℓ) )/βℓ 2
v
(ℓ)
2
=c
(ℓ)
+ v ℓ+1 E (ℓ) y (ℓ) 2
where βℓ = t0 + [t(ℓ) ]T y (ℓ) = (1 − αℓ−1 αℓ−1 )βℓ−1 αℓ = −(tℓ+1 + [t(ℓ) ]T E (ℓ) y (ℓ) )/βℓ z (ℓ) = y (ℓ) + αℓ E (ℓ) y (ℓ) " # z (ℓ) (ℓ+1) y = . αℓ Now we can proceed to the second step (ℓ + 1 → ℓ + 2 = 2(M + 1) + 1), where the Toeplitz system can be expressed as " #" # " # b− ℓ+1 t0 (t(ℓ+1) )∗ v− ℓ+1 2 2 = (ℓ+1) (2.9) t(ℓ+1) T (ℓ+1) v (ℓ+1) b with c(ℓ+2) = [v− ℓ+1 , (v (ℓ+1) )T ]T = c(M+1) . Observe that (2.9) cannot be transformed 2 to a system of the form (2.8) by simple permutations, i.e. just by interchanging rows and columns. Since we have already solved the systems T (ℓ+1) c(ℓ+1) = b(ℓ+1) and T (ℓ+1) y (ℓ+1) = −t(ℓ+1) we can write v (ℓ+1) = (T (ℓ+1) )−1 (b(ℓ+1) − t(ℓ+1) v− ℓ+1 ) = c(ℓ+1) + v− ℓ+1 y (ℓ+1) 2
2
and v− ℓ+1 =(b− ℓ+1 − [t(ℓ+1) ]∗ v (ℓ+1) )/t0 2
2
=(b− ℓ+1 − [t(ℓ+1) ]∗ c(ℓ+1) − [t(ℓ+1) ]∗ v− ℓ+1 y (ℓ+1) )/t0 2
2
=(b− ℓ+1 − [t(ℓ+1) ]∗ c(ℓ+1) )/βℓ+1 , 2
where we have used in the last step that T (ℓ) = [T (ℓ)]∗ which implies that (t(ℓ+1) )∗ y (ℓ+1) is real and therefore t0 + (t(ℓ+1) )∗ y (ℓ+1) = t0 + (t(ℓ+1) )T y (ℓ+1) = βℓ+1 .
8
Thomas Strohmer
Note that at each level M we have to check if the stopping criterion (LSP2) is satisfied. The evaluation of the expression N −1 X
|p(M) (xj ) − sj |2 wj
(2.10)
j=0
can be considerably simplified and by avoiding the evaluation of p(M) at the nonuniformly spaced points xj we can reduce the computational effort from O(M N ) to O(M ) operations. −1 N with the To do this we define the subspace R = {p(xj )}N j=0 : p ∈ P M ⊆ C PN −1 N weighted inner product hy, ziR = j=0 yj z¯j wj for y, z ∈ C . The solution of the −1 N least squares problem (LSP1) is the orthogonal projection of the vector {sj }N j=0 ∈ C onto R and therefore must satisfy h{p(M) (xj )} − {sj }, {p(M) (xj )}iR =
N −1 X
(p(M) (xj ) − sj )p(M) (xj )wj = 0
j=0
which implies h{p(M) (xj )}, {sj }iR = h{p(M) (xj )}, {p(M) (xj )}iR .
(2.11)
Since N −1 X
|sj − p(M) (xj )|2 wj =
N −1 X
|sj |2 wj − 2 Re hs, {p(M) (xj )}iR +
=
N −1 X
|p(M) (xj )|2 wj
j=0
j=0
j=0
N −1 X
|sj |2 wj −
N −1 X
|p(M) (xj )|2 wj
j=0
j=0
by (2.11), and because N −1 X
|p
(M)
2
(xj )| wj =
N −1 X
wj
j=0
j=0
=
M X
X M
2πimxj c(M) m e
m=−M
M X
= hT
c
,c
(M) cn e2πinxj
n=−M
(M)
(M) cm cn
NX −1
wj e2πi(m−n)xj
j=0
m=−M n=−M (M) (M)
X M
(M)
i = hb
(M)
, c(M) i ,
(2.12)
it follows that N −1 X j=0
|sj − p(M) (xj )|2 wj =
N −1 X
|sj |2 wj − hb(M) , c(M) i .
(2.13)
j=0
PN −1 Since j=0 |sj |2 wj has to be computed only once at the beginning of the algorithm, the evaluation of (2.10) can be carried out in O(M ) operations. Summing up we have arrived at following algorithm to compute p∗ .
Regularized trigonometric approximation
9
Algorithm 1 (Levinson-Galerkin algorithm for trigonometric polynomials). Let −1 N −1 the sampling points {xj }N j=0 , sampling values {sj }j=0 , weights wj > 0 and the perturbation estimate δ be given. Then the trigonometric polynomial p∗ that solves (2.1) can be computed in O(N M + M 2 ) operations by following algorithm. PN −1 PN −1 PN −1 PN −1 Initialize: t0 = j=0 wj , t1 = j=0 wj e2πixj , b0 = j=0 sj wj , s = j=0 |sj |2 wj , y (1) = −t1 /t0 , c(1) = b0 /t0 , β0 = t0 , α0 = −t1 /t0 , ε1 = s − b20 /t0 , ℓ = 1.
while εℓ > δ βℓ = (1 − αℓ−1 αℓ−1 )βℓ−1 if ℓ ≡ 1 mod 2 N −1 X sj wj eπi(ℓ+1)xj b ℓ+1 = 2
j=0
b ℓ+1 − hE (ℓ) c(ℓ) , t(ℓ) i
v ℓ+1 =
2
βℓ
2
v
(ℓ)
+ v ℓ+1 E (ℓ) y (ℓ) #2 v (ℓ) = v ℓ+1 " 2 # b(ℓ) = b ℓ+1
=c
c(ℓ+1) b(ℓ+1)
(ℓ)
"
2
elseif ℓ ≡ 0 mod 2 N −1 X b− ℓ = sj wj e−πiℓxj 2
j=0
v− ℓ =
b− ℓ − hc(ℓ) , t(ℓ) i 2
βℓ
2
v
(ℓ)
=c
(ℓ)
c(ℓ+1) = b
(ℓ+1)
=
+ v− ℓ y (ℓ) " #2 v− ℓ 2
v (ℓ)
"
b− ℓ
2
b(ℓ)
#
εℓ+1 = s − hb(ℓ+1) , c(ℓ+1) i end t(ℓ+1) =
N −1 X
wj e2πi(ℓ+1)xj
j=0
αℓ = −
t(ℓ+1) + hE (ℓ) y (ℓ) , t(ℓ) i βℓ
10
Thomas Strohmer
z (ℓ) = y (ℓ) + αℓ E (ℓ) y (ℓ) " # z (ℓ) (ℓ+1) y = αℓ " # t(ℓ) (ℓ+1) t = tℓ+1 ℓ=ℓ+1 end M = ℓ/2 p∗ (t) =
M X
2πimt c(M) m e
m=−M
Remark: Usually one evaluates the final approximation on regularly spaced grid points, hence the last step of the algorithm can be realized by a Fast Fourier transform. The most costly steps are the computation of the entries of t(ℓ) and b(ℓ) . According to Corollary 1 in [9] the entries of T (M) and b(M) can also be computed via FFT by embedding the xj into a regular grid (since the xj can be stored only in finite precision). In this case one automatically gets all entries t0 , . . . , tN −1 at once. However this trick is only useful if the number of points of the regular grid is of the same magnitude as the number of sampling points. Alternatively one may use the numerical attractive formulas of Rokhlin and Dutt [6, 7] or Beylkin [2] for fast evaluation of trigonometric sums at unequally spaced nodes. Algorithm 1 can be simplified for real-valued data, this modification is left to the reader. 3. Multivariate trigonometric approximation. An advantage of the proposed approach, besides its numerical efficiency, is the fact that it can be easily extended to multivariate trigonometric approximation. In this section we briefly discuss some results for the two-dimensional case. 2 We define the space of 2-D trigonometric polynomials PM by M X 2 PM = p : p(x, y) = cj,k e2πi(jx+ky) . (3.1) j,k=−M
To reduce the notational burden, we have assumed in (3.1) that p has equal degree M in each coordinate, the extension to polynomials with different degree in each coordinate is straightforward. −1 2 For an arbitrary sampling set {(xj , yj )}N j=0 ∈ [0, 1) and given degree M the system matrix according to the 2-D version of Theorem 2.1 is [24] (T
(M)
)k,l =
N −1 X j=0
wj e2πi(k−l)(xj +yj ) ,
k, l = 0, . . . 2M .
(3.2)
11
Regularized trigonometric approximation
One can easily verify that T (M) is a block Toeplitz matrix with Toeplitz blocks [24]. In [1] Levinson’s algorithm has been extended to general block Toeplitz systems. With these results at hand, we can easily generalize Algorithm 1 to multivariate trigonometric approximation. Clearly the question of invertibility of T (M) and estimates of the condition number become more difficult. The condition (2M + 1)d ≤ N is necessary in dimension d > 1, but not sufficient any more, since the fundamental theorem of algebra does not hold in dimensions larger than one. In [16] Gr¨ochenig has derived estimates for the condition number of T (M) in higher dimensions. However, these estimates become worse with increasing dimension. Further work remains to be done in that direction. In the following we consider a special case of trigonometric approximation in two dimensions. This case arises when a function is irregularly sampled along lines. A typical example is illustrated in Figure 3.1. Such sampling patterns are encountered for instance in geophysics and medical imaging, see also Section 4.2.
1 0 0 1
1 0 0 1 1 0 1 0
1 0 1 1 0 0 0 1
1 0 1 0
1 0 0 1
1 0 1 0 1 0 0 1
1 0 1 0 1 0 0 1
1 0 0 1
1 0 1 0
1 0 0 1 1 0 1 0 1 0 0 1 1 0 1 0
1 0 1 0
1 0 0 1
1 0 1 0
1 0 0 1
1 0 0 1
1 0 0 1
1 0 1 0
1 0 0 1
1 0 1 0
1 0 0 1
{
0 1 1 0 1 0 0 1
1 0 0 1 1 0 0 1
0 1 1 0 1 0 0 1 0 1 0 1
1 0 0 1
1 0 1 0
1 0 1 0
1 0 0 1 1 0 0 1 1 0 0 1
1 0 0 1
1 0 0 1 1 0 0 1
1 0 0 1 1 0 0 1 11 00 11 00
11 00 11 1 00 0 0 1 1 0 1 0
11 00 00 11
1 0 0 1 1 0 0 1 1 0 1 0 0 1 0 1
δ2
1 0 0 1
1 0 1 0 1 0 0 1 1 0 1 0
1 0 0 1
1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1
1 0 0 1 1 0 0 1
{
y
δ1
1 0 0 1
x
Fig. 3.1. Line-type nonuniform sampling set 2 Corollary 3.1. Let p ∈ PM and let {xj , yjk }, j = 0, . . . , N −1, k = 0, . . . , Nj −1 2 be a sampling set in [0, 1) such that Nj −1
A1 kpk2 ≤
X
|p(yjk )|2 ≤ B1 kpk2
A1 , B1 > 0
(3.3)
k=0
for every p ∈ PM and for all j. Further assume that {xj }j∈Z is a sampling set such that 2
A2 kpk ≤
N −1 X j=0
|p(xj )|2 ≤ B2 kpk2
A2 , B2 > 0
(3.4)
12
Thomas Strohmer
for every p ∈ PM . Then 2
A1 A2 kpk ≤
j −1 N −1 NX X
|p(xj , yjk )|2 ≤ B1 B2 kpk2
(3.5)
j=0 k=0
2 for every p ∈ PM . 1 If {xj } and {yjk } are stable sampling sets with supk (xk+1 − xk ) = δ2 < 2M and 1 2 2 supj,k (yj,k+1 − yjk ) = δ1 < 2M , then Al = (1 − δl ) , Bl = (1 + δl ) , l = 1, 2 and the condition number of the block Toeplitz matrix T (M) is bounded by
κ(T ) ≤
(1 + δ1 )2 (1 + δ2 )2 . (1 − δ1 )2 (1 − δ2 )2
(3.6)
Proof. Let x be fixed. Then y → p(x, y) ∈ PM and hence for all j
A1
Z1
Nj −1 2
|p(xj , y)| dy ≤
X
2
|p(xj , yjk )| ≤ B1
k=0
0
Z1
|p(xj , y)|2 dy
(3.7)
0
by assumption (3.3). It follows that
A1
N −1 Z1 X
2
|p(xj , y)| dy ≤
j=0 0
X
2
|p(xj , yjk )| ≤ B1
N −1 Z1 X
|p(xj , y)|2 dy
(3.8)
j=0 0
j,k
Now let y be fixed. Then x → p(x, y) ∈ PM and
A2
Z1
2
|p(x, y)| dx ≤
N −1 X
2
|p(xj , y)| ≤ B2
j=0
0
Z1
|p(x, y)|2 dx .
(3.9)
0
Since N −1 Z1 X j=0 0
2
|p(xj , y)| dy =
Z1 NX −1 0
|p(xj , y)|2 dy ,
(3.10)
j=0
assertion (3.5) follows by combining (3.8) and (3.9) with (3.10). The estimate of the constants Al , Bl and of the condition number of the block Toeplitz matrix T (M) follow from Theorem 2.1. The proof of Corollary 3.1 is due to Gr¨ochenig [14]. Corollary 3.1 does not only 2 guarantee that p ∈ PM can be recovered from its samples p(xj , yjk ), it provides more. An immediate consequence is, that f can be recovered by an efficient algorithm relying on an successive application of Algorithm 1 and the Gohberg-Semencul representation of the inverse of a Toeplitz matrix. See Section 4.2 for more details and an application in medical imaging.
Regularized trigonometric approximation
13
4. Curve and surface approximation by trigonometric polynomials. Trigonometric polynomials can be used to model the boundary or the surface of smooth objects. Let us consider a two-dimensional object, obtained e.g. by a planar cross-section from a 3-D object and assume that the boundary of this 2-D object is a closed curve in R2 . We denote this curve by f and parametrize it by f (u) = (xu , yu ), where xu and yu are the coordinates of f at “time” u in the x- and y-direction respectively. Obviously we can interpret f as a one-dimensional continuous, complex, and periodic function, where xu represents the real part and yu represents the imaginary part of f (u). It follows from the Theorem of Weierstrass (and from the Theorem of Stone-Weierstrass [22] for higher dimensions) that a continuous periodic function can be approximated uniformly by trigonometric polynomials. If f is smooth, we can fairly assume that trigonometric polynomials of low degree provide an approximation of sufficient precision. Assume that we know only some arbitrary, perturbed points sj = (xuj , yuj ) = f (uj ) + δj , j = 0, . . . , N − 1 of f , and we want to recover f from these points. By a slight abuse of notation we interpret sj as complex number and write sj = xj + iyj .
(4.1)
We relate the curve parameter u to the boundary points sj by computing the distance between two successive points sj−1 , sj via u0 = 0
(4.2)
uj = uj−1 + dj q dj = (xj − xj−1 )2 + (yj − yj−1 )2
(4.3) (4.4)
for j = 1, . . . , N − 1. Via the normalization tj = tj /L with L = uN −1 + dN we force all sampling points to be in [0, 1). Other choices for dj in (4.1) can be found in [5] in conjunction with curve approximation using splines. Having carried out the transformations (4.1)–(4.4), we can solve the problem of recovering the curve f from its perturbed points sj by Algorithm 1. 4.1. Object boundary recovery in Echocardiography. Trigonometric polynomials are certainly not suitable to model the shape of arbitrary objects. However they are often useful in cases where an underlying (stationary) physical process implies smoothness conditions of the object. Typical examples arise in medical imaging, for instance in clinical cardiac studies, where the evaluation of cardiac function using parameters of left ventricular contractibility is an important constituent of an echocardiographic examination [26]. These parameters are derived using boundary tracing of endocardial borders of the Left Ventricle (LV). The extraction of the boundary of the LV comprises two steps, once the ultrasound image of a cross section of the LV is given, see Figure 4.1(a)–(d). First an edge detection is applied to the ultrasound image to detect the boundary of the LV, cf. Figure 4.1(c). However this procedure
14
Thomas Strohmer
may be hampered by the presence of interfering biological structures (such as papillar muscles), the uneveness of boundary contrast, and various kinds of noise [25]. Thus edge detection often provides only a set of nonuniformly spaced, perturbed boundary points rather than a connected boundary. Therefore a second step is required, to recover the original boundary from the detected edge points, cf. Figure 4.1(d). Since the shape of the Left Ventricle is definitely smooth, trigonometric polynomials are particularly well suited to model its boundary. After having transformed the detected boundary points as described in (4.1)–(4.4) we can use Algorithm 1 to recover the boundary. The noise level δ depends on the technical equipment under use, it can be determined from experimental experience. Figure 4.2(a)–(b) demonstrate the importance of determining a proper degree for the approximating polynomial. The approximation displayed in Figure 4.2(a) has been computed by solving (LSP1) where M has been chosen too small, we obviously have underfitted the data. The overfitted approximation obtained by solving (LSP1) using a too large M is shown in Figure 4.2(b). The approximation shown in Figure 4.1(d) has been computed by Algorithm 1, it provides the optimal balance between fitting the data and smoothness of the solution. 4.2. Boundary recovery from a sequence of images. In cardiac clinical studies one is more interested in the behavior of the Left Ventricle over a period of time rather than in a single “snapshot”. Thus for a fixed cross section we are given a sequence of ultrasound images (usually regularly spaced in time) describing the variation of the shape of the LV with time. One cycle from diastole (the state of maximal contraction of the LV), passing systole (the state of maximal expansion) to the next diastole consists typically of about 30 image frames. Since the behavior of the LV is (at least for a short period of time) almost periodic, one can model the varying shape of a fixed cross section of the LV as a distorted two-dimensional torus, which in turn can be interpreted as two-dimensional trigonometric polynomial. Clearly we have to use a different degree for the time coordinate τ and for the spatial coordinate u. Due to interfering biological structures and other distortions it sometimes happens that some of the image frames cannot be used to extract any reliable boundary information. Thus we have to approximate these missing boundaries from the information of the other image frames. To be more precise, assume that an echocardiographic examination provides a sequence of ultrasound images Iτ taken at time points τ = 0, 1, . . . , T − 1, where T is approximately the length of one diastolic cycle (the time points could also be nonuniformly spaced). Assume that some of the images Iτ Nj −1 provide no useful information, so that we can only detect boundary points {sj,k }k=0 −1 from the images Iτj , where {τj }N j=0 is a subset of 0, 1, . . . T − 1. In order to get a complete description of the LV for the time interval [0, T ], we have not only to approximate the boundaries fj from each Ij , but we also have to recover the boundaries corresponding to the missing images. In other words we look for a 2-D trigonometric
Regularized trigonometric approximation
(a) 2-D echocardiography
(b) Cross section of Left Ventricle
(c) Detected boundary points
(d) Recovered boundary of LV computed by Algorithm 1
15
Fig. 4.1. The recovery of the boundary of the Left Ventricle from 2-D ultrasound images is a basic step in echocardiography to extract relevant parameters of cardiac function.
2 polynomial p∗ ∈ PM of appropriate degree M that satisfies p(τj , uj,k ) ≡ (xj,k , yj,k ) where the parameter u is related to sj,k = xj,k + iyj,k by formulas (4.2)–(4.4). This approximation can be computed by the 2-D version of Algorithm 1, as indicated in the beginning of Section 3.
Under certain conditions we can use the 1-D version of Algorithm 1 instead of its 2-D version. As long as the assumptions of Corollary 3.1 are satisfied, we can
16
Thomas Strohmer
(a) Underfitted solution
(b) Overfitted solution
Fig. 4.2. 2 compute p∗ ∈ PM by a successive application of Algorithm 1. We first approximate Nj −1 the boundaries fj for each j separately from its samples {sj,k }k=0 , which yields j (Mj ) different polynomials p ∈ PMj . Having done this, the next step is to recover the missing boundaries at those time points where no information is available. We proceed by approximating successively the missing information “line by line”. We choose u = 0, say, and approximate the missing information from the samples p(Mj ) (u) taken at the time points τj , j = 0, . . . , N − 1. (M) (M) (M) Note that the Toeplitz matrices of the systems Tu cu = bu coincide for all u, since the sampling geometry is constant along the u-coordinate (because we have recovered all samples at each τj ). Thus we have to solve multiple Toeplitz systems with the same system matrix but different right hand side. It is well-known that this can be done efficiently by exploiting the Gohberg-Semencul representation of the inverse of the Toeplitz matrix [11]. In our context this reads as follows. We solve
= b(M) Tu(M) c(M) u u
(4.5)
for one u by Algorithm 1. We can solve now all other systems efficiently by establishing (T (M) )−1 in the Gohberg-Semencul form (4.6) (T (M) )−1 = [L(M) ]∗ L(M) − U (M) [U (M) ]∗ /z0
where L(M) is a lower triangular Toeplitz matrix with z = [z0 , z1 , . . . , z2M ]T as its first column, U (M) is an upper triangular Toeplitz matrix with [0, z1 , . . . , z2M ]T as its last column, z being the first column of (T (M) )−1 . The matrix vector multiplications to (M) (M) compute cu = (T (M) )−1 can now be carried out quickly using the Fast Fourier u bu (M) transform by embedding L and U (M) into circulant matrices.
Regularized trigonometric approximation
17
5. Miscellaneous remarks. For sampling sets with large gaps it can happen that the system T (M) c(M) = b(M) gets ill-conditioned with increasing degree M and therefore Algorithm 1 may become instable [3]. In this case one can use a different, more robust approach, which however comes at higher computational costs [23]. We solve the system T (M) c(M) = b(M) iteratively, e.g. by the conjugate gradient method until a certain stopping criterion is satisfied at iteration k, say, yielding the solution (M) (M+1) ck . We use this solution as initial guess at the next level M + 1 by setting c0 = (M) [0 (ck )T 0]T . The crucial point of such a procedure is to find a suitable stopping criterion at each level to ensure convergence of the iterates, see [23, 20] for more details. The computation of the entries of the Toeplitz matrix in Section 4 involves the nodes uj which in this particular case depend on the (perturbed) samples sj . Therefore not only the right hand side b(M) , but also T (M) is subject to perturbations. Hence in principle one might use the concept of total least squares (see [10]) instead of the standard least squares criterion in (LSP1). A detailed discussion of this modification is beyond the scope of this paper. REFERENCES [1] H. Akaike, Block Toeplitz matrix inversion, SIAM J. Appl. Math., 24 (1973), pp. 234–241. [2] G. Beylkin, On the fast Fourier transform of functions with singularities, Appl. Comp. Harm. Anal., 2 (1995), pp. 363–381. [3] G. Cybenko, The numerical stability of the levinson-durbin algorithm for toeplitz systems of equations, SIAM J. Sci. Statist. Comp., 1 (1980), pp. 303–3190. [4] C. Demeure, Fast QR factorization of Vandermonde matrices, Linear Algebra Appl., 122–124 (1989), pp. 165–194. [5] P. Dierckx, Curve and surface fitting with splines, Monographs on Numerical Analysis, Oxford University Press, 1993. [6] A. Dutt and V. Rokhlin, Fast Fourier transforms for nonequispaced data, SIAM J. Sci. Comp., 14 (1993), pp. 1368–1394. , Fast Fourier transforms for nonequispaced data II, Appl. Comp. Harm. Anal., 2 (1995), [7] pp. 85–10. [8] H. Fassbender, On numerical methods for discrete least-squares approximation by trigonometric polynomials, Math. Comp., 66 (1997), pp. 719–741. ¨ chenig, and T. Strohmer, Efficient numerical methods in non[9] H. G. Feichtinger, K. Gro uniform sampling theory, Numerische Mathematik, 69 (1995), pp. 423–440. [10] R. Fierro, G. Golub, P. Hansen, and D. O’Leary, Regularization by truncated total least squares, SIAM J. Sci. Comp., 18 (1997), pp. 1223–1241. [11] I. Gohberg and A. Semencul, On the inversion of finite Toeplitz matrices and their continuous analogs, Mat. Issled., 2 (1972), pp. 201–233. [12] G. Golub, M. Heath, and G. Wahba, Generalized cross-validation as a method for choosing a good ridge parameter, Technometrics, 21 (1979), pp. 215–223. [13] G. Golub and C. van Loan, Matrix Computations, third ed., Johns Hopkins, London, Baltimore, 1996. ¨ chenig, personal communication. [14] K. Gro [15] , A discrete theory of irregular sampling, Linear Algebra Appl., 193 (1993), pp. 129–150. [16] , Finite and infinite-dimensional models for non-uniform sampling, in SampTA - Sam-
18
Thomas Strohmer
pling Theory and Applications, Aveiro, Portugal, 1997, pp. 285–290. [17] P. Hansen, Analysis of discrete ill-posed problems by means of the L-curve, SIAM Review, 34 (1992), pp. 561–680. [18] N. Levinson, The Wiener RMS (root-mean square) error criterion in filter design and prediction, J. Math. Phys., 25 (1947), pp. 261–278. [19] A. Newbery, Trigonometric interpolation and curve-fitting, Math. Comp., 24 (1970), pp. 869– 876. [20] M. Rauth and T. Strohmer, Smooth approximation of potential fields from noisy scattered data, Geophysics, (to appear). [21] L. Reichel, G. Ammar, and W. Gragg, Discrete least squares approximation by trigonometric polynomials, Math. Comp., 57 (1991), pp. 273–289. [22] W. Rudin, Fourier Analysis on Groups, Wiley Interscience, New York, 1976. [23] O. Scherzer and T. Strohmer, A multi–level algorithm for the solution of moment problems. preprint, 1997. [24] T. Strohmer, Computationally attractive reconstruction of band-limited images from irregular samples, IEEE Trans. Image Proc., 6 (1997), pp. 540–548. [25] M. Suessner, M. Budil, T. Strohmer, M. Greher, G. Porenta, and T. Binder, Contour detection using artifical neuronal network presegmention, in Proc. Computers in Cardiology, Vienna, 1995. [26] D. Wilson, E. Geiser, and J. Li, Feature extraction in 2-dimensional short-axis echocardiographic images, J. Math. Imag. Vision, 3 (1993), pp. 285–298.
