Robust rational approximation with barycentric representations - PowerPoint PPT Presentation

Robust rational approximation with barycentric representations Silviu Filip , CAIRN team, Univ Rennes, Inria, CNRS, IRISA joint work with Bernhard Beckermann, Yuji Nakatsukasa and Lloyd N. Trefethen March 13, 2018 1 / 17

Rational functions Why are they important? → powerful approximations near singularities or on unbounded domains 2 / 17

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 / 17

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 / 17

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 / 17

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 / 17

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 / 17

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 / 17

A classic example: f ( x ) = | x | , x ∈ [ − 1 , 1] 4 / 17

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 / 17

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 / 17

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 / 17

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 / 17

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 / 17

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 / 17

A classic example: f ( x ) = | x | , x ∈ [ − 1 , 1] 10 0 10 -5 10 -10 10 -15 0 10 20 30 40 50 60 70 80 10 -12 6 4 2 0 -2 -4 -6 -8 10 -12 10 -10 10 -8 10 -6 10 -4 10 -2 10 0 6 / 17

Barycentric representations → many different ways of representing rational functions monomial/monomial, Chebyshev/Chebyshev, partial fraction, etc. 7 / 17

Barycentric representations → many different ways of representing rational functions monomial/monomial, Chebyshev/Chebyshev, partial fraction, etc. → 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 7 / 17

Barycentric representations → many different ways of representing rational functions monomial/monomial, Chebyshev/Chebyshev, partial fraction, etc. → 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 Why use adaptive barycentric formulas? → problem dependent { t k } ⇒ well conditioned representation 7 / 17

The rational Remez algorithm 8 / 17

The rational Remez algorithm Some assumptions: no defect ( d = 0 ) → required diagonal case m = n → for convenience 8 / 17

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 8 / 17

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: 8 / 17

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 8 / 17

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 8 / 17

The rational Remez algorithm Convergence: → usually quadratic [Curtis & Osborne 1966] → guaranteed only if starting reference set is close enough to optimal 9 / 17

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 9 / 17

Step 1 : initial reference set → need suff. good initial guess for { x k } 10 / 17

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) , . . . ) 10 / 17

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 11 / 17

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 ) 11 / 17