“Good” points for multivariate polynomial interpolation and approximation Marc Van Barel, Matthias Humet, Laurent Sorber Dept. of Computer Science, KU Leuven, Belgium “Cittadella dei Musei”, Cagliari, Sardinia, Italy September 2-5, 2013
Univariate polynomial interpolation and approximation (1) Consider a “difficult” function f defined on a subset Ω of the real line R , e.g., Ω = [ − 1 , 1 ] . “Difficult” means that it is too costly to work directly with the function f . Instead we can use a function g belonging to a vector space V : ◮ approximating f according to a certain criterion, ◮ cheap to determine, ◮ easier to handle, ◮ numerically sound (conditioning, numerical stability). Possible choices for V are: ◮ polynomials up to a certain degree, ◮ rational functions, ◮ trigonometric functions, ◮ . . . Van Barel, Humet, Sorber (KU Leuven) Multivariate polynomial interpol. and approx. Cagliari, September 5, 2013 2 / 23
Univariate polynomial interpolation and approximation (2) In this talk, ◮ the approximating function g is a polynomial function; ◮ the approximation criterion is the ∞ -norm, i.e., � f � ∞ = max x ∈ Ω | f ( x ) | . Problem: not cheap to determine the minmax-approximant. Solution: compute cheaply a nearly-optimal approximant. An example of such a nearly-optimal approximant is the polynomial interpolating the function f in the points { x k ∈ Ω } N 1 with a corresponding small Lebesgue constant. Van Barel, Humet, Sorber (KU Leuven) Multivariate polynomial interpol. and approx. Cagliari, September 5, 2013 3 / 23
Lebesgue constant (univariate case) subset Ω ⊂ R , e.g., Ω = [ − 1 , 1 ] monomials p k ( x ) = x k − 1 , k = 1 , 2 , . . . , N vector space V = P N = { � N k = 1 c k p k } The Lebesgue function Λ and the Lebesgue constant λ associated with a set of points { x k } N 1 is defined by N � Λ( x ) = | l i ( x ) | λ = max x ∈ Ω Λ( x ) i = 1 where l i ( x ) are the Lagrange polynomials � l i ∈ P N � � l i = δ i , j x j For any function f : � f − p � ∞ ≤ ( 1 + λ ) � f − p ∗ � ∞ Hence, nearly-optimal approximant when λ is small. Van Barel, Humet, Sorber (KU Leuven) Multivariate polynomial interpol. and approx. Cagliari, September 5, 2013 4 / 23
Example: Chebyshev points compared to equispaced points in [ − 1 , 1 ] For more details: see the book Nick Trefethen, Approximation Theory and Approximation Practice , SIAM 2013. Lebesgue function for 30 equispaced points 8 5 10 4 6 10 3 4 10 2 2 10 1 0 10 0 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −1 15.2. Lebesgue constants for polynomial interpolation. . . . For Chebyshev points, they satisfy Van Barel, Humet, Sorber (KU Leuven) Multivariate polynomial interpol. and approx. Cagliari, September 5, 2013 5 / 23 2 2
Lebesgue constant (multivariate case) subset Ω ⊂ R n , e.g., Ω = [ − 1 , 1 ] 2 , a triangle, a sphere, . . . monomials p k ( x 1 , . . . , x n ) = x α 1 ( k ) · · · x α n ( k ) = x α ( k ) n 1 in a certain order ( x , α ∈ R n ) vector space P N = { � N k = 1 c k p k } The Lebesgue function Λ and the Lebesgue constant λ associated with a set of points { x k } N 1 is defined by N � Λ( x ) = | l i ( x ) | λ = max x ∈ Ω Λ( x ) i = 1 where l i ( x ) are the Lagrange polynomials � l i ∈ P N � � l i = δ i , j x j For any function f : � f − p � ∞ ≤ ( 1 + λ ) � f − p ∗ � ∞ Hence, nearly-optimal approximant when λ is small. Van Barel, Humet, Sorber (KU Leuven) Multivariate polynomial interpol. and approx. Cagliari, September 5, 2013 6 / 23
Computing the Lebesgue constant Discretize subset Ω ∈ R n with “grid” G ⊂ Ω containing K points N N � � λ = max | l i ( x ) | − → λ ≈ max | l i ( x ) | x ∈ Ω x ∈ G i = 1 i = 1 Example of a mesh G generated by DistMesh for the L-shape consisting of 3475 points. P .-O. Persson, G. Strang, A Simple Mesh Generator in MATLAB . SIAM Review, Volume 46 (2), pp. 329-345, June 2004 Let { φ k } N 1 be a basis for P N Van Barel, Humet, Sorber (KU Leuven) Multivariate polynomial interpol. and approx. Cagliari, September 5, 2013 7 / 23
Minimize Lebesgue constant � � �� { y k } K min λ ≈ � L N � � 1 � { x k } N ∞ 1 � � � � � � { y k } K { y k } K { x k } N s.t. Φ N = L N Φ N 1 1 1 Objective function is not differentiable → need derivative free optimization Many local minima Large problem size, e.g.: n = 2, total degree of 20 → 231 points ∼ 462 variables [Briani, Sommariva, Vianello; 2011] use state of the art Matlab optimization algorithms to compute minima for the simplex, square and disk. → very slow for larger point sets Van Barel, Humet, Sorber (KU Leuven) Multivariate polynomial interpol. and approx. Cagliari, September 5, 2013 8 / 23
Alternatives to obtain low Lebesgue constant Compromise between a small λ and computational effort Special subsets: small λ with little computational cost E.g., direct formula for Padua points on the square [ − 1 , 1 ] 2 Alternative algorithm (“greedy” algorithm) k � � is maximal on � � First phase : add ( k + 1 ) th point where � l ( k ) , i ( x ) 1 the subset Ω , i.e., on the “grid” G i = 1 ⋆ Remember: λ = max x ∈ Ω � k � � � l ( k ) , i ( x ) � i = 1 k + 1 � = 1 ⋆ By adding the point x ∗ , we assure that � � l ( k + 1 ) , i ( x ∗ ) � � i = 1 Second phase : update points one by one 2 N − 1 � � is maximal ⋆ remove point and add it again where � � � l ( N − 1 ) , i ( x ) i = 1 Van Barel, Humet, Sorber (KU Leuven) Multivariate polynomial interpol. and approx. Cagliari, September 5, 2013 9 / 23
Compare points on the square 231 points (total degree 20), result after 20 iterations. Basis of product Chebyshev polynomials: φ k ( x , y ) = T i ( x ) T j ( y ) PADUA: LC = 9.2 OPTIMAL: LC = 8.14 1 1 0.5 0.5 0 0 −0.5 −0.5 −1 −1 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 ALGORITHM: LC = 18.0 1 0.5 0 −0.5 −1 −1 −0.5 0 0.5 1 Van Barel, Humet, Sorber (KU Leuven) Multivariate polynomial interpol. and approx. Cagliari, September 5, 2013 10 / 23
Compute points for the L-shape 231 points (total degree 20), result after 20 iterations. Basis of product Chebyshev polynomials: φ k ( x , y ) = T i ( x ) T j ( y ) Lebesgue constant: 31 . 5 1 0.8 0.6 0.4 0.2 0 y −0.2 −0.4 −0.6 −0.8 −1 −1 −0.5 0 0.5 1 x Van Barel, Humet, Sorber (KU Leuven) Multivariate polynomial interpol. and approx. Cagliari, September 5, 2013 11 / 23
Compute points for the L-shape 231 points (total degree 20), 20 iterations Show Lebesgue constant after each iteration 1000 900 800 700 600 500 LC 400 300 200 100 0 −100 0 5 10 15 20 iteration Van Barel, Humet, Sorber (KU Leuven) Multivariate polynomial interpol. and approx. Cagliari, September 5, 2013 12 / 23
Some problems with L-shape Product Chebyshev polynomials are not a good basis for the L-shape, see figure Alternative basis: orthogonal polynomials with respect to discrete inner product → need to recompute basis when points change too much 14 10 12 10 10 10 8 10 cond( Φ ) 6 10 4 10 2 10 0 10 0 50 100 150 200 250 k Van Barel, Humet, Sorber (KU Leuven) Multivariate polynomial interpol. and approx. Cagliari, September 5, 2013 13 / 23
Alternative optimization algorithm Instead of minimizing the ∞ -norm � �� � { y k } K min λ ≈ � L N � � 1 � { x k } N ∞ 1 � � � � � � { y k } K { y k } K { x k } N s.t. Φ N = L N Φ N 1 1 1 we minimize the Frobenius norm � �� � { y k } K min λ ≈ � L N � � 1 � { x k } N frob 1 � � � � � � { y k } K { y k } K { x k } N Φ N = L N Φ N s.t. 1 1 1 Efficient Matlab implementation of a bound-constrained nonlinear least-squares optimization problem by projected Gauss-Newton dogleg trust-region method Van Barel, Humet, Sorber (KU Leuven) Multivariate polynomial interpol. and approx. Cagliari, September 5, 2013 14 / 23
Alternative optimization algorithm: square 231 points (total degree 20) Lebesgue constant: 7 . 74 (Padua-points: 9 . 20; “best”: 8 . 14 !!) 1 0.8 0.6 0.4 0.2 0 � 0.2 � 0.4 � 0.6 � 0.8 � 1 � 1 � 0.8 � 0.6 � 0.4 � 0.2 0 0.2 0.4 0.6 0.8 1 Van Barel, Humet, Sorber (KU Leuven) Multivariate polynomial interpol. and approx. Cagliari, September 5, 2013 15 / 23
Alternative optimization algorithm: L-shape 231 points (total degree 20) Lebesgue constant: 15 . 34 1 0.8 0.6 0.4 0.2 0 � 0.2 � 0.4 � 0.6 � 0.8 � 1 � 1 � 0.8 � 0.6 � 0.4 � 0.2 0 0.2 0.4 0.6 0.8 1 Van Barel, Humet, Sorber (KU Leuven) Multivariate polynomial interpol. and approx. Cagliari, September 5, 2013 16 / 23
Recommend
More recommend
