Stability Results for Scattered Data Interpolation by Trigonometric - PowerPoint PPT Presentation

Stability Results for Scattered Data Interpolation by Trigonometric Polynomials Daniel Potts Stefan Kunis Institute for Mathematics Institute for Mathematics University of L¨ ubeck University of L¨ ubeck email: potts@math.uni-luebeck.de kunis@math.uni-luebeck.de http://www.math.uni-luebeck.de/potts http://www.math.uni-luebeck.de/kunis

Content • Basics • ’direct’ Problem, matrix vector multiplication, Vandermonde-like, NFFT • ’inverse’ Problem, solving Vandermonde-like systems, INFFT • Interpolation, Stability • Numerical examples, MRI

Basics Geometry torus, sampling set T := R / Z , ( x j ) j =0 ,...,M − 1 =: X ⊂ T separation distance, fill distance h := j =0 ,...,M − 1 dist ( x j , x j +1 ) , min δ := j =0 ,...,M − 1 dist ( x j , x j +1 ) max Ansatz trigonometric polynomials � � e 2 π i k ( · ) : k = − N 2 , . . . , N T N := span 2 − 1 discrete system � e 2 π i kx j � ˆ f ∈ C M , f ∈ C N 2 − 1 , A := j =0 ,...,M − 1; k = − N 2 ,..., N

Matrix vector multiplication - Vandermonde-like matrix - NFFT f ∈ C N given, compute ˆ N 2 − 1 � ˆ f = A ˆ f k e 2 π i kx j , f , f j = j = 0 , . . . , M − 1 k = − N 2 FFT for M = N equispaced nodes, O ( N log N ) operations FFT for non equispaced nodes (Dutt, Rokhlin; Beylkin; P ., Steidl, Tasche), in O ( N log N + M ) operations

Linear system of equations - iNFFT ,,inverse” problem, f ∈ C M given in A ˆ f ≈ f † = A † f fulfills Moore-Penrose pseudo-inverse solution ˆ f � ˆ � f − A ˆ f � 2 → min f � 2 = min . subject to j special case IDFT, Gauß quadrature, M = N, x j = M − 0 . 5 A H W ˆ f = A H W f A = I ⇒ ���� M I 1 direct solver: Reichel, Ammar, Gragg; Faßbender

Approximation problem weighted approximation problem, ω j > 0 , W = diag ( ω j ) M − 1 j =0 , ˆ f � A ˆ f − f � W → min weighted normal equation of first kind ˆ A H W A f = A H W f � �� � Toeplitz dense sampling set δ := j =0 ,...,M − 1 dist ( x j , x j +1 ) max Feichtinger, Gr¨ ochenig, Strohmer (weighted normal equation of first kind ( N < 1 δ )) � 1 + δN � 2 � � A H W A ≤ cond 2 1 − δN

Interpolation problem vanishing residual, i.e. A ˆ f − f = 0 , � N ≥ M (damped) minimisation problem, 0 < ˆ W := diag (ˆ ω k ) k = − N 2 ,..., N 2 − 1 N 2 − 1 � f k | 2 =: � ˆ ˆ k | ˆ f A ˆ ω − 1 f � 2 ˆ → min f = f subject to W − 1 ˆ k = − N 2 example N = 50 , M = 5 nodes, � f � L 2 and � f � L 2 + � f ′ � L 2 minimal, resp. 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 −0.5 0 0.5 −0.5 0 0.5 = 1 + (2 πk ) 2 ω − 1 ω k = 1 ˆ ˆ k normal equation of second kind W A H ˜ W A H ˜ A ˆ f = ˆ ˆ f = f , f

Interpolation with polynomial kernels N 2 − 1 � e − 2 π i kx ˆ ω k e 2 π i ky , K ( x − y ) := k = − N 2 M − 1 � f ( y ) = α l K ( x l − y ) l =0 1 1 1 1 0.5 0.5 0.5 0.5 0 0 0 0 −0.5 0 0.5 −0.5 0 0.5 −0.5 0 0.5 −0.5 0 0.5 Dirichlet Fejer Cesaro Sobolev centers of the kernels and nodes for interpolation are equal � W A H � A ˆ j,l = K ( x j − x l )

Stability aim: find bounds dependent only on N, h for � W A H � � W A H � A ˆ A ˆ λ = λ min , Λ = λ max norm equivalence � f � 2 f ∈ T N ,f ( x j )= f j � f � 2 2 ∼ inf W − 1 ˆ Marcinkiewicz-Zygmund-inequality Λ − 1 � f � 2 f ∈ T N ,f ( x j )= f j � f � 2 W − 1 ≤ λ − 1 � f � 2 2 ≤ inf 2 . ˆ equispaced nodes, circulant interpolation matrix, eigenvalue characterisa- tion � W A H � � A ˆ λ j = M ω k ˆ k = j mod M

Arbitrary nodes if K (0) = 1 and C β | K ( x ) | ≤ N β | x | β � � − 1 2 , 1 for x ∈ , then the interpolation matrix ( K ( x j − x l )) j,l =0 ,...,M − 1 has 2 bounded eigenvalues 1 − 2 ζ ( β ) C β ≤ λ ≤ 1 ≤ Λ ≤ 1 + 2 ζ ( β ) C β N β h β N β h β where ζ denotes the Riemann-zeta-function � k � for ˆ ω k ≈ g and under mild assumptions on g , N � � g ( β − 1) � ( ζ ( β ) + 1) � V C β ≤ . (2 π ) β − 2 1 − β ζ ( β ) | g ( β − 1) | V

explicit estimates for Dirichlet’s kernel 1 − (1 − ln 2 h ) ≤ λ ≤ 1 ≤ Λ ≤ 1 + (1 − ln 2 h ) Nh Nh Fejer’s kernel π 2 π 2 1 − 3 N 2 h 2 ≤ λ ≤ 1 ≤ Λ ≤ 1 + 3 N 2 h 2 Jackson’s kernel 16 π 4 16 π 4 1 − 45 N 4 h 4 ≤ λ ≤ 1 ≤ Λ ≤ 1 + 45 N 4 h 4 condition number equispaced arbitrary Nh +(1 − ln 2 h ) Nh +1 Dirichlet Nh − 1 Nh − (1 − ln 2 h ) N 2 h 2 + π 2 N 2 h 2 +1 Fejer 3 N 2 h 2 − π 2 N 2 h 2 − 1 3

Multivariate setting torus, metric T d := R d / Z d , dist ( x , y ) := min j ∈ Z d � ( x + j ) − y � ∞ normal equation, kernels � W A H � A ˆ j,l = K ( x j − x l ) example, Jackson’s kernel, N = 22 C β ∞ for x ∈ T d , then if K ( 0 ) = 1 and | K N ( x ) | ≤ N β � x � β 1 − 2 dζ ( β ) C β N β h β + d − 1 ≤ λ ≤ 1 ≤ Λ ≤ 1 + 2 dζ ( β ) C β N β h β + d − 1

Iterative methods Landweber iteration f l +1 = ˆ ˆ f l + α ˆ W ˆ z l steepest descent z H l ˆ α l = ˆ W ˆ z l v H l W v l conjugated gradient � � W A H W r 0 , ˆ ˆ W A H W A ˆ W A H W r 0 , . . . K l ( A , ˆ r 0 ) := span � r l − r † � W → min CGNR: † � ˆ � ˆ f l − ˆ W − 1 → min CGNE: f residuals r l = f − A ˆ v l = A ˆ z l = A H W r l , f l , W ˆ z l ˆ

approximation problem, N ≤ δ − 1 A H W A ˆ f = A H W f ACT, CGNR (Feichtinger, Gr¨ ochenig, Strohmer) � r l − r † � W ≤ 2 ( Nδ ) l � r 0 − r † � W interpolation problem, N ≥ ch − 1 W A H ˜ W A H ˜ ˆ A ˆ f = ˆ f = f , f CGNE (P ., Kunis) � c β � l † � ˆ † � ˆ � ˆ f l − ˆ � ˆ f 0 − ˆ W − 1 ≤ 2 f f W − 1 N β h β

Examples Franke function, M = 100000 random nodes, N = 512 , L 2 and Sobolev- type, CGNE image reconstruction, M = 30000 random nodes, N = 256 , multiquadric- type, CGNR

Glacier contour data M = 8345 points, N = 256 , multiquadric-type, CGNE

spiral MRI, reconstruction data points INFFT

Software available: NFFT – C subroutine library (Kunis, P . 2002–2004) http://www.math.uni-luebeck.de/potts/nfft Features – Implemented transforms for d dimensions – Arbitrary-size transforms – Works on any platform with a C compiler and the FFTW package – iterative solution of the inverse transform (LANDWEBER, STEEPEST- DESCENT, CGNR, CGNE) NFFT 2.0 manual online Fast Fourier transform at nonequispaced knots, A user’s guide to a C-library (Kunis, P .)

Conclusions • NFFT fast computation of NFFT • iterative method for solving Vandermonde-like systems, i • Applications – MRI – Radon transform (P ., Steidl 2002) – polar FFT, polar IFFT – next step � f − A ˆ f � W + λ � ˆ f � ˆ W → min