Rational minimax approximation via adaptive barycentric representations Silviu Filip , CAIRN team, Univ Rennes, Inria, CNRS, IRISA joint work with Bernhard Beckermann, Yuji Nakatsukasa and Lloyd N. Trefethen January 22, 2018 1 / 18
Rational functions Why are they important? → powerful approximations near singularities or on unbounded domains 2 / 18
Rational functions Why are they important? → powerful approximations near singularities or on unbounded domains Some applications: elementary + special functions recursive filter design matrix exponentials & stiff PDEs optimal control problems ... 2 / 18
Rational minimax approximation Input: f ∈ C ([ a, b ]) , target type ( m, n ) ∈ N 2 � p � Output: r ∗ ∈ R m,n = q , p ∈ R m [ x ] , q ∈ R n [ x ] s.t. � f − r ∗ � ∞ is minimal . → denote this minimax error with E m,n ( f ) 3 / 18
Rational minimax approximation Input: f ∈ C ([ a, b ]) , target type ( m, n ) ∈ N 2 � p � Output: r ∗ ∈ R m,n = q , p ∈ R m [ x ] , q ∈ R n [ x ] s.t. � f − r ∗ � ∞ is minimal . → denote this minimax error with E m,n ( f ) → theoretical results: 3 / 18
Rational minimax approximation Input: f ∈ C ([ a, b ]) , target type ( m, n ) ∈ N 2 � p � Output: r ∗ ∈ R m,n = q , p ∈ R m [ x ] , q ∈ R n [ x ] s.t. � f − r ∗ � ∞ is minimal . → denote this minimax error with E m,n ( f ) → theoretical results: existence & unicity of r ∗ [de la Vallée Poussin, Walsh] 3 / 18
Rational minimax approximation Input: f ∈ C ([ a, b ]) , target type ( m, n ) ∈ N 2 � p � Output: r ∗ ∈ R m,n = q , p ∈ R m [ x ] , q ∈ R n [ x ] s.t. � f − r ∗ � ∞ is minimal . → denote this minimax error with E m,n ( f ) → theoretical results: existence & unicity of r ∗ [de la Vallée Poussin, Walsh] presence of the defect : → if r ∗ = p ∗ /q ∗ in irreducible form, then its defect is d = min { m − deg p ∗ , n − deg q ∗ } 3 / 18
Rational minimax approximation Input: f ∈ C ([ a, b ]) , target type ( m, n ) ∈ N 2 � p � Output: r ∗ ∈ R m,n = q , p ∈ R m [ x ] , q ∈ R n [ x ] s.t. � f − r ∗ � ∞ is minimal . → denote this minimax error with E m,n ( f ) → theoretical results: existence & unicity of r ∗ [de la Vallée Poussin, Walsh] presence of the defect : → if r ∗ = p ∗ /q ∗ in irreducible form, then its defect is d = min { m − deg p ∗ , n − deg q ∗ } Alternation Theorem [Achieser 1930]: → f − r ∗ equioscillates at least m + n + 2 − d times 3 / 18
A classic example: f ( x ) = | x | , x ∈ [ − 1 , 1] 4 / 18
A classic example: f ( x ) = | x | , x ∈ [ − 1 , 1] → consider best uniform approximations: degree 8 polynomial Error curve 0.04 0.03 0.02 0.01 0 -0.01 -0.02 -0.03 -0.04 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 4 / 18
A classic example: f ( x ) = | x | , x ∈ [ − 1 , 1] → consider best uniform approximations: degree 8 polynomial type (4 , 4) rational function Error curve 0.04 0.03 0.02 0.01 0 -0.01 -0.02 -0.03 -0.04 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 4 / 18
A classic example: f ( x ) = | x | , x ∈ [ − 1 , 1] → asymptotic behavior E n, 0 ( f ) ∼ β/n, β = 0 . 2801 ... [Varga & Carpenter 1985] E n,n ( f ) ∼ 8 e −√ n , [Newman 1964, Stahl 1993] 5 / 18
A classic example: f ( x ) = | x | , x ∈ [ − 1 , 1] → asymptotic behavior E n, 0 ( f ) ∼ β/n, β = 0 . 2801 ... [Varga & Carpenter 1985] E n,n ( f ) ∼ 8 e −√ n , [Newman 1964, Stahl 1993] → rational minimax approximations can be difficult to compute e.g. [Varga, Ruttan & Carpenter 1991] conjecture Stahl’s result using 200-digit arithmetic for n � 80 5 / 18
A classic example: f ( x ) = | x | , x ∈ [ − 1 , 1] → asymptotic behavior E n, 0 ( f ) ∼ β/n, β = 0 . 2801 ... [Varga & Carpenter 1985] E n,n ( f ) ∼ 8 e −√ n , [Newman 1964, Stahl 1993] → rational minimax approximations can be difficult to compute e.g. [Varga, Ruttan & Carpenter 1991] conjecture Stahl’s result using 200-digit arithmetic for n � 80 → codes (the Remez algorithm): Maple: numapprox[minimax] Mathematica: MinimaxApproximation 5 / 18
A classic example: f ( x ) = | x | , x ∈ [ − 1 , 1] → asymptotic behavior E n, 0 ( f ) ∼ β/n, β = 0 . 2801 ... [Varga & Carpenter 1985] E n,n ( f ) ∼ 8 e −√ n , [Newman 1964, Stahl 1993] → rational minimax approximations can be difficult to compute e.g. [Varga, Ruttan & Carpenter 1991] conjecture Stahl’s result using 200-digit arithmetic for n � 80 → codes (the Remez algorithm): Maple: numapprox[minimax] Mathematica: MinimaxApproximation Chebfun (Matlab): minimax 5 / 18
Barycentric representations → many different ways of representing rational functions 6 / 18
Barycentric representations → many different ways of representing rational functions → barycentric form for type ( n, n ) rational functions n n � r ( z ) = N ( z ) α k β k � � D ( z ) = z − t k z − t k k =0 k =0 Notation: { α k } , { β k } barycentric coefficients { t k } support points 6 / 18
Barycentric representations Why use adaptive barycentric formulas? → problem dependent { t k } ⇒ well conditioned representation 7 / 18
Barycentric representations Why use adaptive barycentric formulas? → problem dependent { t k } ⇒ well conditioned representation Example: → the adaptive Antoulas-Anderson (AAA) algorithm [Nakatsukasa, Sète & Trefethen 2018]: greedy least squares approximation 7 / 18
Example : f ( x ) = | x | , x ∈ [ − 1 , 1] , type (20 , 20) p/q vs N/D 2 1 | D | q 0 -1 -0.5 0 0.5 1 10 0 | D | 10 -5 10 -10 q 10 -15 -1 -0.5 0 0.5 1 8 / 18
The rational Remez algorithm 9 / 18
The rational Remez algorithm Some assumptions: no defect ( d = 0 ) → required diagonal case m = n → for convenience 9 / 18
The rational Remez algorithm Some assumptions: no defect ( d = 0 ) → required diagonal case m = n → for convenience Step 1: choose a reference set a � x 0 < · · · < x 2 n +1 � b 9 / 18
The rational Remez algorithm Some assumptions: no defect ( d = 0 ) → required diagonal case m = n → for convenience Step 1: choose a reference set a � x 0 < · · · < x 2 n +1 � b → iterate the following steps until convergence: 9 / 18
The rational Remez algorithm Some assumptions: no defect ( d = 0 ) → required diagonal case m = n → for convenience Step 1: choose a reference set a � x 0 < · · · < x 2 n +1 � b → iterate the following steps until convergence: Step 2: find r ∈ R n,n and λ ∈ R s.t. f ( x k ) − r ( x k ) = ( − 1) k +1 λ, k = 0 , . . . , 2 n + 1 9 / 18
The rational Remez algorithm Some assumptions: no defect ( d = 0 ) → required diagonal case m = n → for convenience Step 1: choose a reference set a � x 0 < · · · < x 2 n +1 � b → iterate the following steps until convergence: Step 2: find r ∈ R n,n and λ ∈ R s.t. f ( x k ) − r ( x k ) = ( − 1) k +1 λ, k = 0 , . . . , 2 n + 1 Step 3: among local extrema of f − r , take 2 n + 2 new points a � x ′ 0 < · · · < x ′ 2 n +1 � b, f − r alternates in sign + at least one global extrema over [ a, b ] and | f ( x ′ k ) − r ( x ′ k ) | � | λ | , k = 0 , . . . , 2 n + 1 9 / 18
The rational Remez algorithm Convergence: → usually quadratic [Curtis & Osborne 1966] → guaranteed only if starting reference set is close enough to optimal 10 / 18
The rational Remez algorithm Convergence: → usually quadratic [Curtis & Osborne 1966] → guaranteed only if starting reference set is close enough to optimal What can go wrong? → no pole-free solution in Step 2 10 / 18
Step 1 : initial reference set → need suff. good initial guess for { x k } 11 / 18
Step 1 : initial reference set → need suff. good initial guess for { x k } Our strategy: use Carathéodory-Fejér (CF) approximation [Trefethen & Gutknecht 1983] 11 / 18
Step 1 : initial reference set → need suff. good initial guess for { x k } Our strategy: use Carathéodory-Fejér (CF) approximation [Trefethen & Gutknecht 1983] AAA-Lawson approx. (adaptively re-weighted least squares AAA variant) 11 / 18
Step 1 : initial reference set → need suff. good initial guess for { x k } Our strategy: use Carathéodory-Fejér (CF) approximation [Trefethen & Gutknecht 1983] AAA-Lawson approx. (adaptively re-weighted least squares AAA variant) extrapolation from lower degree approx. ( (2 , 2) , (3 , 3) , (4 , 4) , . . . ) 11 / 18
Step 2 : find r → find r = N/D ∈ R n,n s.t. N ( x k ) = D ( x k )( f ( x k ) − ( − 1) k +1 λ ) , k = 0 , . . . , 2 n + 1 12 / 18
Step 2 : find r → find r = N/D ∈ R n,n s.t. N ( x k ) = D ( x k )( f ( x k ) − ( − 1) k +1 λ ) , k = 0 , . . . , 2 n + 1 → matrix form f ( x 0 ) − 1 f ( x 1 ) 1 Cα = − λ Cβ, ... − 1 ... f ( x 2 n +1 ) C ∈ R (2 n +2) × ( n +1) Cauchy matrix, C k,j = 1 / ( x k − t j ) 12 / 18
Recommend
More recommend
