Pseudospectral Fourier reconstruction with IPRM Karlheinz Grchenig - PowerPoint PPT Presentation

Pseudospectral Fourier reconstruction with IPRM Karlheinz Gröchenig Tomasz Hrycak European Center of Time-Frequency Analysis Faculty of Mathematics University of Vienna http://homepage.univie.ac.at/karlheinz.groechenig/ Dagstuhl, December 2008 Karlheinz Gröchenig (EUCETIFA) IPRM December 2008 1 / 21

Outline What is IPRM? 1 IPRM-Algorithm, Condition Number, Optimality 2 Numerical Simulations 3 Variations 4 logo.eps Karlheinz Gröchenig (EUCETIFA) IPRM December 2008 2 / 21

What is IPRM? IPRM I • Gibbs phenomenon • IPRM = Inverse polynomial reconstruction method (Jung and Shizgal, 2003–07) • Goal: Construct an algebraic polynomial from Fourier coefficients • Find an approximation of a piecewise smooth function from given Fourier coefficients • How many Fourier coefficients are required for accurate construction of algebraic polynomial? • Compression • Relation between Fourier basis and other bases • Gottlieb, Shu; Gelb, Tanner; Tadmore; etc. logo.eps Karlheinz Gröchenig (EUCETIFA) IPRM December 2008 3 / 21

What is IPRM? IPRM II Given function f on [ − 1 , 1 ] and m consecutive Fourier coefficients � 1 1 −⌊ m − 1 ⌋ ≤ k ≤ ⌊ m f ( x ) e − i π kx dx , ˆ f ( k ) = √ 2 ⌋ . 2 2 − 1 Find a polynomial p of degree n − 1 with these Fourier coefficients. Expand p into normalized Legendre polynomials � P k n − 1 � a l � p = P l l = 0 IPRM: solve the system n − 1 � k = −⌊ m − 1 ⌋ , . . . , ⌊ m a l � P l ( k ) = � � f ( k ) 2 ⌋ . 2 l = 0 logo.eps Karlheinz Gröchenig (EUCETIFA) IPRM December 2008 4 / 21

IPRM-Algorithm, Condition Number, Optimality IPRM-Algorithm Input: m Fourier coefficients � f ( k ) , Let A m , n be m × n matrix A m , n with entries � √ a kl = � � 2 ( − ı ) l l + 1 P l ( k ) = 2 j l ( k π ) , (1) k = −⌊ m − 1 2 ⌋ , . . . , ⌊ m 2 ⌋ , l = 0 , . . . , n − 1. Solve overdetermined least squares problem for approximate 1 Legendre coefficients c = [ c 0 , . . . , c n − 1 ] t c ∈ C n � A m , n c − [ � f ( d ) , . . . , � f ( D )] t � 2 , min (2) where d = −⌊ m − 1 2 ⌋ , D = ⌊ m 2 ⌋ . Approximate f by truncated Legendre series 2 n − 1 � c l � f n = P l . (3) l = 0 logo.eps Karlheinz Gröchenig (EUCETIFA) IPRM December 2008 5 / 21

IPRM-Algorithm, Condition Number, Optimality Existence of a Reconstruction Theorem Let d and D be integers such that d ≤ 0 ≤ D, and let p ∈ P M have vanishing D − d + 1 consecutive Fourier coefficients � p ( d ) = � p ( d + 1 ) = . . . = � p ( D − 1 ) = � p ( D ) = 0 . (4) If D − d + 1 ≥ M + 1 , then p = 0 identically. REMARK: A n , n is invertible, and for m > n A m , n has full rank. P M is the space of algebraic polynomials of degree at most M . logo.eps Karlheinz Gröchenig (EUCETIFA) IPRM December 2008 6 / 21

IPRM-Algorithm, Condition Number, Optimality Stability of the Reconstruction Theorem For every α > 1 , every n = 1 , 2 , . . . , and every integer m > α n 2 , the � α condition number of the matrix A m , n does not exceed α − 1 . REMARK: α > 1 can be pushed to α > c for some c ≈ 1 / 2. logo.eps Karlheinz Gröchenig (EUCETIFA) IPRM December 2008 7 / 21

IPRM-Algorithm, Condition Number, Optimality Convergence Rates Theorem Let f = � ∞ l = 0 a l � P l with Legendre coefficients | a l | ≤ ce − β l , (5) where c > 0 and β > 0 , and let f n be the reconstruction by IPRM (3) . If m > n 2 , then � f − f n � ∞ ≤ c ′ ne − β n , (6) for another constant c ′ > 0 . REMARK: Measured by number of Fourier coefficients m = α n 2 , the convergence is root-exponential: � f − f n � ∞ ≤ c ′ √ me − β √ m . logo.eps Karlheinz Gröchenig (EUCETIFA) IPRM December 2008 8 / 21

Numerical Simulations 18 10 16 10 14 10 12 condition number 10 10 10 8 10 6 10 4 10 2 10 0 10 0 50 100 150 number of Legendre polynomials Figure: Condition numbers of the square matrix A n , n for n = 1 , . . . , 150. logo.eps Karlheinz Gröchenig (EUCETIFA) IPRM December 2008 9 / 21

Numerical Simulations Computation versus Proof Experimentally: Smallest singular value λ min ( n ) of A n , n decays exponentially (equivalently: condition number of the square matrix grows exponentially) Current estimate: λ min ≤ 0 . 65 Needed: Behavior of Bessel functions J ν in the non-asymptotic region ν ≤ x ≤ ν 2 . logo.eps Karlheinz Gröchenig (EUCETIFA) IPRM December 2008 10 / 21

Numerical Simulations 8 7 6 condition number 5 4 3 2 1 0 10 20 30 40 50 60 70 80 90 100 number of Legendre polynomials Figure: Condition numbers of the matrix A ⌈ n 2 ⌉ , n for n = 1 , . . . , 100. 3 logo.eps Karlheinz Gröchenig (EUCETIFA) IPRM December 2008 11 / 21

Numerical Simulations 10 10 8 10 condition number 6 10 4 10 2 10 0 10 0 10 20 30 40 50 60 70 80 90 100 number of Legendre polynomials 20 , 1 1 40 , 1 Figure: Condition numbers of the matrix A ⌈ α n 2 ⌉ , n for α = 60 . logo.eps Karlheinz Gröchenig (EUCETIFA) IPRM December 2008 12 / 21

Numerical Simulations 1.25 1.2 condition number 1.15 1.1 1.05 1 0 10 20 30 40 50 60 70 80 90 100 number of Legendre polynomials Figure: Condition numbers of the matrix A n 2 , n for n = 1 , . . . , 100. logo.eps Karlheinz Gröchenig (EUCETIFA) IPRM December 2008 13 / 21

Numerical Simulations 3 10 condition number 2 10 1 10 0 10 0 10 20 30 40 50 60 70 80 90 100 number of Fourier coefficients Figure: Condition numbers of the matrix A m , 20 with m Fourier coefficients for m = 1 , . . . , 100. logo.eps Karlheinz Gröchenig (EUCETIFA) IPRM December 2008 14 / 21

Numerical Simulations 2 10 m = n 2 m = n 1 10 0 10 relative maximum error −1 10 −2 10 −3 10 −4 10 −5 10 −6 10 0 5 10 15 20 25 30 35 40 45 50 number of Legendre polynomials 1 Figure: Relative maximum reconstruction errors for the function x − 0 . 3 ı on the interval [ − 1 , 1 ] with the two versions of IPRM. logo.eps Karlheinz Gröchenig (EUCETIFA) IPRM December 2008 15 / 21

Numerical Simulations 0 10 m = n 2 m = n relative maximum error −5 10 −10 10 −15 10 0 5 10 15 20 25 30 35 40 45 50 number of Legendre polynomials 1 Figure: Relative maximum reconstruction errors for the function x − 1 . 0 ı on the interval [ − 1 , 1 ] with the two versions of IPRM. logo.eps Karlheinz Gröchenig (EUCETIFA) IPRM December 2008 16 / 21

Numerical Simulations 0 10 m = n m = n 2 −2 10 relative maximum error −4 10 −6 10 −8 10 −10 10 0 5 10 15 20 25 30 35 40 45 50 number of Fourier coefficients 1 Figure: Relative maximum reconstruction errors for the function x − 1 . 0 ı on the interval [ − 1 , 1 ] with the two versions of IPRM. logo.eps Karlheinz Gröchenig (EUCETIFA) IPRM December 2008 17 / 21

Variations Piecewise Polynomials from Fourier Coefficients Fix nodes − 1 = a 0 < a 1 < . . . < a L − 1 < a L = 1 , (7) and consider P M , a = { f : f | ( a j − 1 , a j ) is polynomial of degree M } dim P a , M = L ( M + 1 ) Theorem Let d and D be integers such that d ≤ 0 ≤ D, and let p ∈ P M , a have vanishing D − d + 1 consecutive Fourier coefficients � p ( d ) = � p ( d + 1 ) = . . . = � p ( D − 1 ) = � p ( D ) = 0 . (8) If D − d + 1 ≥ L ( M + 1 ) , then p = 0 identically. logo.eps Karlheinz Gröchenig (EUCETIFA) IPRM December 2008 18 / 21

Variations Piecewise Constant Functions with free nodes L � p = p j χ ( t j − 1 , t j ) . (9) j = 1 Theorem Let p a step function on [ − 1 , 1 ] with at most L − 1 points of discontinuity, and let d, and D ∈ Z be such that d ≤ 0 ≤ D. If D − d + 1 ≥ 2 L − 1 , then p is uniquely determined by its D − d + 1 consecutive Fourier coefficients � p ( d ) , � p ( d + 1 ) , . . . , � p ( D − 1 ) , � p ( D ) . Reconstruction by Prony’s spectral estimator, Used in compressed sensing by M. Vetterli as “Occam’s razor” logo.eps Karlheinz Gröchenig (EUCETIFA) IPRM December 2008 19 / 21

Variations To Do List • Optimality of order of condition number • Open question: Is n →∞ κ ( A α n 2 , n ) = e β/α lim • Condition numbers for piecewise polynomials with fixed nodes • Variable degrees for piecewise polynomials with fixed nodes • Method for piecewise polynomials with free nodes • Reconstruction from arbitrary frequencies, from random frequencies logo.eps Karlheinz Gröchenig (EUCETIFA) IPRM December 2008 20 / 21

Variations Summary • Rigorous convergence analysis of IPRM • First proof of existence of the square IPRM (invertibility of A n , n ) • n × n IPRM is acceptable for entire functions n 2 × n IPRM is reliable for meromorphic functions • n 2 × n IPRM useful in applications because it handles noisy signals • and uses all available Fourier coefficients Thank you! Further questions also to tomasz.hrycak@univie.ac.at logo.eps Karlheinz Gröchenig (EUCETIFA) IPRM December 2008 21 / 21