Iterative algorithms to approximate canonical Gabor windows: Computational aspects A.J.E.M. Janssen Philips Research Laboratories WO-02, 5656AA Eindhoven, The Netherlands 6 0 Peter L. Søndergaard 0 ∗ 2 Technical University of Denmark, Department of Mathematics, Building 303, 2800 n a Lyngby, Denmark. J 5 ] A Abstract F . h In this paper we investigate the computational aspects of some recently proposed t a iterative methods for approximating the canonical tight and canonical dual window m of a Gabor frame (g,a,b). The iterations start with the window g while the itera- [ tion steps comprise the window g, the kth iterand γ , the frame operators S and k 1 S corresponding to (g,a,b) and (γ ,a,b), respectively, and a number of scalars. k k v The structure of the iteration step of the method is determined by the envisaged 0 convergenceordermofthemethod.Weconsidertwostrategies forscalingtheterms 9 0 in the iteration step: norm scaling, where in each step the windows are normalized, 1 and initial scaling where we only scale in the very beginning. Norm scaling leads to 0 fast, but conditionally convergent methods, while initial scaling leads to uncondi- 6 0 tionally convergent methods, but with possibly suboptimal convergence constants. / h The iterations, initially formulated for time-continuous Gabor systems, are con- t sidered and tested in a discrete setting in which one passes to the appropriately a m sampled-and-periodized windows and frame operators. Furthermore, they are com- : pared with respect to accuracy and efficiency with other methods to approximate v i canonical windows associated with Gabor frames. X r a Key words: Gabor frame, tight window, dual window, iterative method, scaling, adjoint orbit, Zak transform. MSC: 42C15; 41A25; 47A58; 94A12 Corresponding author. ∗ Email addresses: [email protected](A.J.E.M. Janssen), [email protected](Peter L. Søndergaard). Preprint submitted to Elsevier Science 22nd December 2005 1 Introduction Weconsiderinthispaperiterativeschemesfortheapproximationofthecanon- ical tight and canonical dual windows associated with a Gabor frame. We refer to [6, Ch. 8-10] and [11, Ch. 5-9, 11-13] for recent and comprehensive treatments of the theory of Gabor systems and frames; to fix notations and conventions we briefly give here the main features. We denote for g L2(R) ∈ and a > 0, b > 0 by (g,a,b) the collection of time-frequency shifted windows g , m,n Z, (1.1) na,mb ∈ where for x,y R we denote ∈ g = e2πiytg(t x), t R. (1.2) x,y − ∈ We refer to g as the window and to a and b as the time-shift and frequency- shift parameters, respectively, of the Gabor system (g,a,b). When there are A > 0, B < , called the lower and upper frame bound, respectively, such ∞ that for all f L2(R) there holds ∈ A f 2 ∞ (f,g ) 2 B f 2, (1.3) na,mb k k ≤ | | ≤ k k m,n= X−∞ we call (g,a,b) a Gabor frame. When in (1.3) the second inequality holds for all f L2(R), we have that ∈ f L2(R) Sf := ∞ (f,g )g (1.4) na,mb na,mb ∈ 7→ m,n= X−∞ is well-defined as a bounded, positive, semi-definite linear operator of L2(R). We call S the frame operator of (g,a,b). The frame operator commutes with all relevant shift operators, i.e., we have for all f L2(R) ∈ Sf = (Sf) , m,n Z. (1.5) na,mb na,mb ∈ We shall assume in the remainder of this paper that (g,a,b) is a Gabor frame. Thus the frame operator S is positive definite and therefore boundedly invert- ible.TherearetwowindowscanonicallyassociatedtotheGaborframe(g,a,b). These are the canonical tight window gt and the canonical dual window gd, defined by gt = S 1/2g,gd = S 1g, (1.6) − − 2 respectively. The practical relevance of these windows is that they give rise to Gabor series representations of arbitrary f L2(R) according to ∈ f = ∞ f,gt gt = ∞ f,gd g , (1.7) na,mb na,mb na,mb na,mb m,nX=−∞(cid:16) (cid:17) m,nX=−∞(cid:16) (cid:17) where both series are L2(R)-convergent. Furthermore, the Gabor systems (gt,a,b)and gd,a,b areGaborframesthemselves withframeoperatorsequal to the identit(cid:16)y I and(cid:17)S 1, respectively. − The computation of gt and gd according to (1.6) requires taking the inverse square root and the inverse of the frame operator S, respectively. In the often occurring practical case that ab is a rational number, the frame operator is highly structured which allows relatively efficient methods for computing S 1, − 1 see [27]. The computationof S−2 is much more awkward, even inthe case that ab is rational, and often requires advanced techniques from numerical linear algebra, see for instance [13]. In [20] the calculus of Gabor frame operators was combined with a use of the spectral mapping theorem and Kantorovich’s inequality to analyze an iteration scheme for the approximation of gt that was proposed around 1995 by Feichtinger and Strohmer (independently of one another). In this iteration scheme one sets γ = g and for k = 0,1,... 0 1 1 1 1 I. γ = α γ + β S 1γ ; α = , β = , (1.8) k+1 2 k k 2 k k− k k kγkk k Sk−1γk (cid:13) (cid:13) (cid:13) (cid:13) where S is the frame operator corresponding to (γ ,a(cid:13),b). It(cid:13)was shown in k k [20] that (γ ,a,b) is indeed a Gabor frame, and that γk converges (at least) k γk quadratically to gt . From the numerical results in [2k0]kfor the case that g is gt k k the standard Gaussian window 21/4exp( πt2) and a = b = 1/√2 it appears − that the resulting method compares favorably with other iterative techniques for computing inverse square roots, [24,25,29]. The investigations in [20] were followed by the introduction in [18,19] of two familiesofiterative algorithmsfortheapproximationofgt andgd,inwhich the iteration step involves the initial window g and the frame operator S as well as the current window γ and frame operator S , but frame operator inversion k k as in (1.8) do not occur. The following instances of these two families were analyzed in [18,19]. Again we set γ = g, and for k = 0,1,... 0 II. 3 1 1 1 γ = α γ β S γ ; α = , β = , (1.9) k+1 k k k k k k k 2 − 2 γ S γ k k k k k k k 3 15 5 3 γ = ε γ ε S γ + ε S2γ ; III. k+1 8 k0 k − 4 k1 k k 8 k2 k k 1 1 1 ε = , ε = , ε = , (1.10) k0 γ k1 S γ k2 S2γ k kk k k kk k k kk for the approximation of gt, and IV. 1 1 γ =2α γ β S g; α = , β = , (1.11) k+1 k k k k k k − γ S g k k k k k k γ =3δ γ 3δ S g +δ SS γ ; k+1 k0 k k1 k k2 k k V. − 1 1 1 δ = , δ = , δ = , (1.12) k0 k1 k2 γ S g SS γ k k k k k k k k k k for the approximation of gd. The algorithms II-V are, contrary to algorithm I, conditionally convergent in the sense that the frame bound ratio A of the initial Gabor frame (g,a,b) B should exceed a certain lower bound. Accordingly, in algorithm II and III we have that γk converges to gt quadratically and cubically when A > 1 and γk gt B 2 A > 3,respkecktively. Inalgorikthmk IVandVwe havethat γk converges to gd B 7 γk gd k k k k quadratically and cubically when A > 1 √5 1 and A > 0.513829766..., B 2 − B respectively. A remarkable phenomenon (cid:16)that em(cid:17)erged from the preliminary experiments done with the algorithms around 2002, was the fact that the lower bounds for the algorithms II, III seem far too pessimistic while those for the algorithms IV and V appear to be realistic. In the algorithms just presented, all computed windows are normalized. We shall refer to this as norm scaling. Another possibility that we will investigate, is to replace all scalars α’s, β’s, ε’s and δ’s that occur in (1.9-1.12) by 1, and then initially scale the windows by replacing g by g/Bˆ1/2,S by S/Bˆ. (1.13) ˆ We shall refer to this scaling strategy as initial scaling. If B is (an estimate for) the best upper frame bound maxσ(S) then the algorithms will be un- conditionally convergent, with guaranteed desired convergence order, but with convergence constants that may not be as good as the ones that can be ob- tained by using the norm scaling as described by (1.9-1.12). Matrix versions of algorithms II-III without any scaling have been treated in [5]. In [14,23] the matrix version of algorithm II is considered using norm scal- 4 ing and a scaling method that approximates the optimal scaling. The matrix version of algorithm IV is known as a Schulz iteration, see [28]. The fact that S, and therefore ϕ(S) with ϕ continuous and positive on the spectrum of S, commutes with allrelevant shift operators,allows us to formulatethe iteration steps on the level of the windows themselves. In this paper we investigate the algorithms II-V, using both norm and initial scaling, with more emphasis on computational aspects than in [20,18,19]. Here it is necessary to consider sampled-and-periodized Gabor systems in the style of [16]. This allows for a formulation and analysis of the algorithms I-V in an entirely similar way as was done in [20,18,19]. Thanks to the fact that the involved (canonical) windows and frame operators behave so conveniently un- der theoperationsof sampling andperiodization,the observations doneonthe sampled-and-periodized systems are directly relevant to the time-continuous systems. Wemust restrictheretorationalvaluesofab,andthisgives theframe operator additional structure which can be exploited in the computations as dictated by the recursion steps, also see [2,27,30] for this matter. The notions ”smart”(or, rather, ”smart but risky”) and ”safe”(or, rather, ”safe but conservative”) were introduced in a casual way in [19] to distinguish be- tween cases where the stationary point(s) of the function transforming (frame) operators according to (4.2), (4.9) has a good chance to be well-placed in the middle of and on the ”safe” side of the relevant spectral set, respectively. In the present paper, we choose to refer to the strategies leading to smart and safe modes as ”norm scaling” and ”initial scaling”, respectively, and discard the terms ”smart”and ”safe”altogether. 2 Paper outline and results In Section 3 and 4 we present the basic results of [18,19] on transforming g into γ = ϕ(S)g, where ϕ is a function positive and continuous on σ(S), so as to obtain a window γ whose frame operator S (for approximating gt) or the γ operator (SS )1/2 =: Z (for approximating gd) is closer to (a multiple of) the γ γ identity I than S itself. Here we recall that (gt,a,b) has frame operator I and that gd,a,b has frame operator S 1. Thus we relate the operators S and − Sγ an(cid:16)d their(cid:17)frame bounds, and we present a bound for the distance between (the normed) γ and gt in terms of the frame bounds of S . Similarly, we relate γ the frame bounds of (g,a,b) and the minimum and maximum of σ(Z ), and γ we present a bound for the difference between (the normed) γ and gd in terms of the latter minimum and maximum. Next, in Section 4, the choice of ϕ is specified so as to accommodate the recursions of type II, III and of type IV, V which gives us a means to monitor the frame bound ratio Ak (for gt) and Bk 5 of the ratio between minimum and maximum of σ(Z ) (for gd) during the k iteration process. In Section 5 we present the algorithms using only initial scaling. In Section 6 we elaborate on the observation that all algorithms take place in the closed linear span of the adjoint orbit g j,l Z . Here the g j/b,l/a L ∈ dual lattice representation of frame operators is relnevant as(cid:12)well asoan operator (cid:12) norm to measure the distance of S (for gt) and of (SS )(cid:12)1/2 (for gd) from (a k k multiple of) the identity. The consideration of the algorithms in the space g L reveals a fundamental difference between the algorithms for computing gt and gd that manifests itself in the totally different after-convergence behaviour of the two families of algorithms. In Section 7 we give some considerations in the Zak transform domain, so as to produce examples of Gabor frames for which a specific algorithm diverges. In Section 8 we discuss the discretization and finitization aspects (through sampling and periodization) that have to be taken into account since the algorithms are to be tested numerically. In Section 9 we show how the algorithms can be expressed for discrete, finite Gabor systems, and show that the algorithms are scalar iterations of the sin- gular values of certain matrices. We present an efficient implementation of the iterative algorithms, and we list the window functions we have used to test the algorithms. InSection10wepresentourexperimentalresults,comparethemwithwhatthe theory predicts and with other methods to compute tight and dual windows. We provide examples that show the quadratic and cubic convergence of the algorithms, and the exponential divergence of the dual iterations after the initial convergence. We give an example that breaks the norm scaling schemes for both the tight and dual iterations, and show how various error norms of the iteration step behave. Comparisons with other methods are made: We show that the tight iterations are competitive with respect to computing time and superior with respect to precision. Finally, we show that the number of iterations needed for full convergence of the algorithms are dependent on the frame bound ratio, but independent of the structural properties of the discretization. For initial scaling, we show that it is easy to choose a scaling parameter that gives almost optimal convergence. 6 3 Frame operator calculus and basic inequalities The basic theory to analyze the recursions appears somewhat scattered in [20,18,19]; for the reader’s convenience, we give in Section 3 and 4 a concise yet comprehensive summary of the basic results and ideas. We let (g,a,b) be a GaborframewithframeoperatorS andbestframeboundsA = minσ(S) > 0, B = maxσ(S),whereσ(S)isthespectrumofS.Inthissectionwepresent the basic inequalities expressing the approximation errors in terms of the (frame) boundsontheinvolved(frame)operators.Theseinequalitiesareaconsequence of the calculus of Gabor frame operators, the spectral mapping theorem and Kantorovich’s inequality. Proposition 1 Let ϕ be continuous and positive on [A,B], and set γ := ϕ(S)g. The following holds. (i) (γ,a,b) is a Gabor frame with frame operator S := Sϕ2(S) and best γ frame bounds A := min sϕ2(s), B := max sϕ2(s). (3.1) γ γ s σ(S) s σ(S) ∈ ∈ Furthermore, gt = S 1/2g = S 1/2γ = γt, (3.2) − γ− and γ gt 2 A 1 Q1/4 ; Q = γ. (3.3) (cid:13) γ − gt (cid:13) ≤ − γ s1+Qγ γ Bγ (cid:13)k k k k(cid:13) (cid:16) (cid:17) (cid:13) (cid:13) (ii) Let Z := ((cid:13)SS )1/2 = Sϕ(cid:13)(S), and γ (cid:13) γ (cid:13) E :=minσ(Z ) = min sϕ(s), (3.4) γ γ s σ(S) ∈ F :=maxσ(Z ) = max sϕ(s). (3.5) γ γ s σ(S) ∈ Then gd = Z 1γ,Z γ = S g, (3.6) γ− γ γ and γ gd 2 E 1 R1/2 ; R = γ. (3.7) (cid:13) γ − gd (cid:13) ≤ − γ s1+Rγ γ Fγ (cid:13)k k k k(cid:13) (cid:16) (cid:17) (cid:13) (cid:13) (cid:13) (cid:13) The basic result(cid:13)(i) gives us a(cid:13) clue how to produce a good approximation γ γ of gt : Take ϕ such that sϕ2(s) is flat on σ(S) [A,B] so that the numbkekr gt ⊂ k k Q in (3.3) is close to 1. Similarly, by the basic result (ii), the number R in γ γ (3.7) is close to 1 when ϕ is such that sϕ(s) is flat on σ(S) [A,B], and then we obtain a good approximation of gd . In the next two su⊂bsections, we gd k k 7 use this basic result repeatedly with polynomials ϕ of fixed degree m so as to obtain iterative approximations of gt and gd. 4 Norm scaling 4.1 Iterations for approximating gt We consider iteration schemes γ =g; γ = ϕ (S )γ , k = 0,1,..., (4.1) 0 k+1 k k k for the approximation of gt, where S is the frame operator of (γ ,a,b). We k k use here the basic result (i) repeatedly with g=γ , S = S and γ = γ , S = S = S ϕ2(S ). (4.2) k k k+1 γ k+1 k k k For k = 0,1,... we have that γt = gt and that k γ gt 2 A k 1/4 k 1 Q ; Q = , (4.3) (cid:13) γk − gt (cid:13) ≤ − k s1+Qk k Bk (cid:13)k k k k(cid:13) (cid:16) (cid:17) (cid:13) (cid:13) (cid:13) (cid:13) where A and(cid:13)B are the b(cid:13)est frame bounds of S . The numbers A , B can k k k k k be computed and estimated recursively according to A = A, B = B and 0 0 A = min sϕ2(s) min sϕ2(s), (4.4) k+1 s σ(Sk) k ≥ s [Ak,Bk] k ∈ ∈ B = max sϕ2(s) max sϕ2(s), (4.5) k+1 s σ(Sk) k ≤ s [Ak,Bk] k ∈ ∈ for k = 0,1,.... We should choose ϕ such that sϕ2(s) is flat on σ(S ). To that end there is k k k proposed in [19, Subsec. 5.1] for m = 2,3,... the choice m 1 ϕ (s) = − a α sj, α = Sjγ −1, (4.6) k mj kj kj k k jX=0 (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) where the a are defined by mj m 1 1/2 m 1 − ( 1)l − (1 x)l = − a xj, x > 0. (4.7)   mj − − l=0 l j=0 X X     8 The motivation for this choice is as follows. The left-hand side of (4.7) is the (m 1)th order Taylor approximation of x 1/2 around x = 1, while (when Q − k − is sufficiently close to 1) j S 1 α Sj k , j = 0,...,m 1. (4.8) kj k≈ Sk ! γk − k k k k 2 Hence sϕ2(s) = s1/2ϕ (s) should be expected to be flat on σ(S ), with k k k 1 Q potentia(cid:16)lly of orde(cid:17)r (1 Q )m. k+1 k − − When m = 2,3 we get the iterations II, III in (1.9), (1.10). It is shown in [18, Sec. 4] and [19, Sec. 6] that for m = 2 the quantity Q = Ak increases to 1 and k Bk that γk gt quadratically when k , provided that A > 1. For m = 3 γk → gt → ∞ B 2 it is skhowk n ink[1k8, Sec. 8] that Q increases to 1 and that γk gt cubically k γk → gt when k , provided that A > 3. In [18,19] it was obksekrvedkfork m = 2,3 → ∞ B 7 that the choice of α causes sϕ2(s) to have one or more stationary points in kj k [A ,B ] so that the odds for flatness of sϕ2(s) on [A ,B ] are favourable. k k k k k 4.2 Iterations for approximating gd We consider iteration schemes γ =g; γ = ϕ (Z )γ , k = 0,1,..., (4.9) 0 k+1 k k k for the approximation of gd, where Z = (SS )1/2 with S the frame operator k k k of (γ ,a,b). It is seen from (4.9) by induction that γ = ψ (S)g for some k k k function ψ . Hence by the basic result (i) with ϕ = ψ , we have that k k S = Sψ2(S),Z = Sψ (S), (4.10) k k k k and, by the basic result (ii), that γ gd 2 E k 1/2 k 1 R ; R = , (4.11) (cid:13) γk − gd (cid:13) ≤ − k s1+Rk k Ek (cid:13)k k k k(cid:13) (cid:16) (cid:17) (cid:13) (cid:13) (cid:13) (cid:13) where Ek =(cid:13)minσ(Zk), F(cid:13)k = maxσ(Zk). A further use of the calculus of frame operators as given by the basic result (i), yields Z = Z ϕ (Z ). (4.12) k+1 k k k Consequently, the numbers E , F can be computed and estimated recursively k k according to E = A, F = B and 0 0 9 E = min zϕ (z) min zϕ (z), (4.13) k+1 k k z σ(Zk) ≥ z [Ek,Fk] ∈ ∈ F = max zϕ (z) max zϕ (z), (4.14) k+1 k k z σ(Zk) ≤ z [Ek,Fk] ∈ ∈ for k = 0,1,.... We should choose ϕ such that zϕ (z) is flat on σ(Z ). To that end there is k k k proposed in [19, Subsec. 5.2] for m = 2,3,... the choice m 1 ϕ (z) = − b β zj, β = Zjγ −1, (4.15) k mj kj kj k k jX=0 (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) where the b are defined by mj m 1 m 1 − (1 x)l = − b xj, x > 0. (4.16) mj − l=0 j=0 X X Themotivationfortheproposalissimilartotheoneforthechoiceofϕ in(4.6) k in Subsec. 4.1; we now note that the left-hand side of (4.16) is the (m 1)th − order Taylor approximation of x 1 around x = 1. The implementation of the − resulting recurrence step m 1 Zjγ γ = − b k k , Z = (SS )1/2, (4.17) k+1 mj Zjγ k k jX=0 k k (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) is made feasible by the observation that, thanks to the second item in (3.6), Z γ = S g so that k k k Z2rγ =(SS )rγ , Z2r+1γ = (SS )rS g, r = 0,1,.... (4.18) k k k k k k k k When m = 2,3 we get the recurrences IV, V in (1.11) and (1.12). It is shown in [18, Sec. 5] and [19, Sec. 7], that for m = 2 the quantity R = Ek increases k Fk to 1 and that γk gd quadratically when k , provided that A > γk → gd → ∞ B k k k k 1 √5 1 . For m = 3 it is shown in [18, Sec. 9] that R increases to 1 and 2 − k th(cid:16)at γk (cid:17) gd cubically when k , provided that A > 0.513829766.... γk → gd → ∞ B k k k k In [18,19] it was observed for m = 2,3 that the choice of β causes zϕ (z) to kj k have one or more stationary points in [E ,F ]. k k 10