Fast and stable rational RBF-based Partition of Unity interpolation 1 - PowerPoint PPT Presentation

Fast and stable rational RBF-based Partition of Unity interpolation 1 Stefano De Marchi Department of Mathematics “Tullio Levi-Civita” University of Padova SMART 2017 - Gaeta September 18, 2017 1 joint work with A. Martinez (Padova) and E. Perracchione (Padova)

Outline 1 From RBF to Rational RBF (RRBF) Partion of Unity, VSK, Computational issues 2 3 Numerical experiments Future work 4 2 of 27

From RBF to Rational RBF (RRBF) 3 of 27

RBF Approximation Notation 1 Data: Ω ⊂ R d , X ⊂ Ω , test function f X N = { x 1 , . . . , x N } ⊂ Ω f = { f 1 , . . . , f N } , where f i = f ( x i ) 4 of 27

RBF Approximation Notation 1 Data: Ω ⊂ R d , X ⊂ Ω , test function f X N = { x 1 , . . . , x N } ⊂ Ω f = { f 1 , . . . , f N } , where f i = f ( x i ) 2 Approximation setting: kernel K ε , N K (Ω) , N K ( X N ) ⊂ N K (Ω) kernel K = K ε , Strictly Positive Definite (SPD) and radial Examples: globally supported: K ε ( x , y ) = e − ( ε � x − y � ) 2 ( gaussian ), + [ 4 ε 2 � x − y � 2 + 1 ] locally supported: K ε ( x , y ) = ( 1 − ε 2 � x − y � 2 ) 4 ( C 2 ( R 2 ) Wendland ) native space N K (Ω) (where K is the reproducing kernel) finite subspace N K ( X N ) = span { K ( · , x ) : x ∈ X N } ⊂ N K (Ω) 4 of 27

RBF Approximation Notation 1 Data: Ω ⊂ R d , X ⊂ Ω , test function f X N = { x 1 , . . . , x N } ⊂ Ω f = { f 1 , . . . , f N } , where f i = f ( x i ) 2 Approximation setting: kernel K ε , N K (Ω) , N K ( X N ) ⊂ N K (Ω) kernel K = K ε , Strictly Positive Definite (SPD) and radial Examples: globally supported: K ε ( x , y ) = e − ( ε � x − y � ) 2 ( gaussian ), + [ 4 ε 2 � x − y � 2 + 1 ] locally supported: K ε ( x , y ) = ( 1 − ε 2 � x − y � 2 ) 4 ( C 2 ( R 2 ) Wendland ) native space N K (Ω) (where K is the reproducing kernel) finite subspace N K ( X N ) = span { K ( · , x ) : x ∈ X N } ⊂ N K (Ω) Aim Find P f ∈ N K ( X N ) s.t. P f ≈ f 4 of 27

RBF Interpolation Given Ω , X N , f N � Interpolant: P f ( x ) = α k K ( x , x k ) , x ∈ Ω , x k ∈ X N . [Hardy and k = 1 � Gofert 1975] used multiquadrics K ( x , y ) = 1 + ǫ 2 � x − y � 2 . 5 of 27

RBF Interpolation Given Ω , X N , f N � Interpolant: P f ( x ) = α k K ( x , x k ) , x ∈ Ω , x k ∈ X N . [Hardy and k = 1 � Gofert 1975] used multiquadrics K ( x , y ) = 1 + ǫ 2 � x − y � 2 . � N P f ( x ) = P f ( x ) k = 1 α k K ( x , x k ) Rescaled interpolant: ˆ P g ( x ) = where P g � N k = 1 β k K ( x , x k ) is the kernel based interpolant of g ( x ) = 1 , ∀ x ∈ Ω [Deparis at al 2014, DeM et al 2017]. Shepard’s PU method. 5 of 27

RBF Interpolation Given Ω , X N , f N � Interpolant: P f ( x ) = α k K ( x , x k ) , x ∈ Ω , x k ∈ X N . [Hardy and k = 1 � Gofert 1975] used multiquadrics K ( x , y ) = 1 + ǫ 2 � x − y � 2 . � N P f ( x ) = P f ( x ) k = 1 α k K ( x , x k ) Rescaled interpolant: ˆ P g ( x ) = where P g � N k = 1 β k K ( x , x k ) is the kernel based interpolant of g ( x ) = 1 , ∀ x ∈ Ω [Deparis at al 2014, DeM et al 2017]. Shepard’s PU method. Rational interpolant: R ( x ) = R ( 1 ) ( x ) � N k = 1 α k K ( x , x k ) R ( 2 ) ( x ) = � N k = 1 β k K ( x , x k ) [Jackbsson et al. 2009, Sarra and Bai 2017]. 5 of 27

Rational RBF General framework polynomial case d = 1. a m x m + · · · + a 0 x 0 r ( x ) = p 1 ( x ) p 2 ( x ) = . x n + b n − 1 x n − 1 · · · + b 0 M = m + n + 1 unknowns (Pad´ e approximation). If N > M to get the coefficients we may solve the LS problem  N  � 2  � �    min � f ( x k ) − r ( x k )  .   � �    �  p 1 ∈ Π 1 m , p 2 ∈ Π 1    n k = 1 6 of 27

Rational RBF General framework polynomial case d = 1. a m x m + · · · + a 0 x 0 r ( x ) = p 1 ( x ) p 2 ( x ) = . x n + b n − 1 x n − 1 · · · + b 0 M = m + n + 1 unknowns (Pad´ e approximation). If N > M to get the coefficients we may solve the LS problem  N  � 2  � �    min � f ( x k ) − r ( x k )  .   � �    �  p 1 ∈ Π 1 m , p 2 ∈ Π 1    n k = 1 Let X m = { x k , . . . , x k + m − 1 } , X n = { x j , . . . , x j + n − 1 } ⊂ X N be non empty, such that m + n ≤ N � k + m − 1 R ( x ) = R ( 1 ) ( x ) α i 1 K ( x , x i 1 ) i 1 = k R ( 2 ) ( x ) = , (1) � j + n − 1 β i 2 K ( x , x i 2 ) i 2 = j provided R ( 2 ) ( x ) � 0 , x ∈ Ω . 6 of 27

Rational RBF Posedness In the case m + n = N , asking R ( x i ) = f i , R ( 1 ) ( x i ) − R ( 2 ) ( x i ) f i = 0 = ⇒ B ξ = 0 with ξ = ( α, β ) and B = B 1 + B 2 . 7 of 27

Rational RBF Posedness In the case m + n = N , asking R ( x i ) = f i , R ( 1 ) ( x i ) − R ( 2 ) ( x i ) f i = 0 = ⇒ B ξ = 0 with ξ = ( α, β ) and B = B 1 + B 2 . Obs : depending on X m , X n and on the function values f , we might have cancellations and B can be singular. We must look for non-trivial solutions! 7 of 27

Rational RBF Posedness In the case m + n = N , asking R ( x i ) = f i , R ( 1 ) ( x i ) − R ( 2 ) ( x i ) f i = 0 = ⇒ B ξ = 0 with ξ = ( α, β ) and B = B 1 + B 2 . Obs : depending on X m , X n and on the function values f , we might have cancellations and B can be singular. We must look for non-trivial solutions! Example Given X N and f . Let X m ∩ X n � ∅ , with m + n = N . Suppose that f i = a , i = 1 , . . . , N , a ∈ R , then the system B ξ = 0 admits infinite solutions. 7 of 27

Example Example N = 26 , m = n = 13, Ω = [ 0 , 1 ] , f i = 1 (Left) or f be piecewise constant (Right) . Following [GolubReinsch1975] non-trivial solution by asking || ξ || 2 = 1 i.e. solving constrained problem ξ ∈ R N , || ξ || 2 = 1 || B ξ || 2 . min Figure : The black dots represent the set of scattered data, the red solid line and the blue dotted one are the curves reconstructed via the rational RBF and classical RBF approximation, by the Gaussian kernel. 8 of 27

Rational RBF Find the coefficients: I [Jackobsson et al 2009] proved the well-posedness of R ( x ) = R ( 1 ) ( x ) � N i = 1 α i K ( x , x i ) R ( 2 ) ( x ) = , � N k = 1 β k K ( x , x k ) B is now a sum of two squared blocks of order N i.e. B is N × 2 N K ( x 1 , x 1 ) · · · K ( x 1 , x N ) − f 1 K ( x 1 , x 1 ) · · · − f 1 K ( x 1 , x N )       . . .   B =  . . .  .    . . .          K ( x N , x 1 ) K ( x N , x N ) − f N K ( x N , x 1 ) − f N K ( x N , x N )  · · · · · ·  The system B ξ = 0 can be written as ( A − DA )( ξ ) = 0 with D = diag ( f 1 , . . . , f N ) , and A i , j = K ( x i , x j ) . But underdetermined! 9 of 27

Rational RBF Find the coefficients: II Since R ( 1 ) ( x i ) = f i R ( 2 ) ( x i ) , find q = ( R ( 2 ) ( x 1 ) , . . . , R ( 2 ) ( x N )) T and then get p = ( R ( 1 ) ( x i ) , ..., R ( 1 ) ( x N )) T as p = D q . 10 of 27

Rational RBF Find the coefficients: II Since R ( 1 ) ( x i ) = f i R ( 2 ) ( x i ) , find q = ( R ( 2 ) ( x 1 ) , . . . , R ( 2 ) ( x N )) T and then get p = ( R ( 1 ) ( x i ) , ..., R ( 1 ) ( x N )) T as p = D q . If p, q are given then the rational interpolant is known by solving A α = p , A β = q . (2) 10 of 27

Rational RBF Find the coefficients: II Since R ( 1 ) ( x i ) = f i R ( 2 ) ( x i ) , find q = ( R ( 2 ) ( x 1 ) , . . . , R ( 2 ) ( x N )) T and then get p = ( R ( 1 ) ( x i ) , ..., R ( 1 ) ( x N )) T as p = D q . If p, q are given then the rational interpolant is known by solving A α = p , A β = q . (2) Existence+Uniqueness of (2) since K is SPD Using the native space norms the above problem is equivalent Problem 1   1 || R ( 1 ) || 2 N K + || R ( 2 ) || 2   min  . (3)      N K  || f || 2 R ( 1 ) , R ( 2 ) ∈N K ,  2 1 / || f || 2 2 || p || 2 2 + || q || 2 2 = 1 , R ( 1 ) ( x k )= f k R ( 2 ) ( x k ) . 10 of 27

Rational RBF Find the coefficients: III || R ( 1 ) || 2 N K = α T A α , || R ( 2 ) || 2 N K = β T A β . and Then, from (2) and symmetry of A || R ( 1 ) || 2 N K = p T A − 1 p , || R ( 2 ) || 2 N K = q T A − 1 q . and Therefore, the Problem 1 reduces to solve Problem 2   1  q T D T A − 1 D q + q T A − 1 q  min  .       || f || 2 q ∈ R N ,  2 1 / || f || 2 2 || D q || 2 2 + || q || 2 2 = 1 . 11 of 27

Rational RBF Find the coefficients: IV [Jackbsson 2009] show that this is equivalent to solve the following generalized eigenvalue problem Problem 3 Σ q = λ Θ q , with 1 1 D T A − 1 D + A − 1 , D T D + I N , Σ = Θ = and || f || 2 || f || 2 2 2 where I N is the identity matrix. 12 of 27

Partion of Unity, VSK, Computational issues 13 of 27

RRBF- PU interpolation Ω = ∪ s j = 1 Ω i (can overlap): that is, each x ∈ Ω belongs to a limited number of subdomains, say s 0 < s . Then we consider, W j non-negative fnct on Ω j , s.t. � s j = 1 W j = 1 is a partition of unity . RRBF-PU interpolant s � I ( x ) = R j ( x ) W j ( x ) , ∀ x ∈ Ω (4) j = 1 where R j is a RRBF interpolant on Ω j i.e. R ( 1 ) � N j i = 1 α j i K ( x , x j ( x ) i ) j R j ( x ) = = , R ( 2 ) � N j k = 1 β j k K ( x , x j ( x ) k ) j N j = | Ω j | , x j k ∈ X N j , k = 1 , . . . , N j . 14 of 27