Meshfree Approximation Methods With MATLAB
Meshfree Approximation Methods With MATLAB
e o
587?
I
S
Meshfree Approximation
Methods with MATLAB
1
INTERDISCIPLINARY MATHEMATICAL SCIENCES
Published
V o l . 1: G l o b a l A t t r a c t o r s o f N o n a u t o n o m o u s D i s s i p a t i v e D y n a m i c a l Systems
David N. Cheban
V o l . 2: Stochastic D i f f e r e n t i a l E q u a t i o n s : T h e o r y a n d Applications
A V o l u m e i n H o n o r o f Professor B o r i s L . R o z o v s k i i
eds. Peter H. Baxendale & Sergey V. Lototsky
V o l . 3: A m p l i t u d e E q u a t i o n s f o r Stochastic P a r t i a l D i f f e r e n t i a l E q u a t i o n s
Dirk Blomker
V o l . 5: T h e H i l b e r t - H u a n g T r a n s f o r m and Its A p p l i c a t i o n s
Norden E. Huang & Samuel S. P. Shen
G r e g o r y E . F a s s h a u e r r _ s
W r i f l r l ^ H o n t i f i n
. lb. /
Published by
W o r l d
fif-?, Y~ Scientific Publishing Co. Pte. Ltd.
5 Toh Tuck Link, Singapore 596224
U S A c e : 2 7 W a r r e n
/ t>&^~ ^ Street, Suite 401-402, Hackensack, NJ 07601
/] UK office: 57 Shelton Street, Covent Garden, London W C 2 H 9 H E
M E S H F R E E APPROXIMATION M E T H O D S W I T H MATLAB
(With CD-ROM)
Interdisciplinary Mathematical Sciences Vol. 6
All rights reserved. This book, or parts thereof, may not be reproduced in any form or by any means, electronic or
mechanical, including photocopying, recording or any information storage and retrieval system now known or to
be invented, without written permission from the Publisher.
For photocopying of material in this volume, please pay a copying fee through the Copyright Clearance Center,
Inc., 222 Rosewood Drive, Danvers, M A 01923, U S A . In this case permission to photocopy is not required from
the publisher.
ISBN-13 978-981-270-633-1
ISBN-10 981-270-633-X
ISBN-13 978-981-270-634-8 (pbk)
ISBN-10 981-270-634-8 (pbk)
vii
viii Meshfree Approximation Methods with M A T L A B
as M A T L A B ? H O W d o t h e y c o m p a r e w i t h t r a d i t i o n a l m e t h o d s ? H o w d o t h e v a r i o u s
flavors o f meshfree m e t h o d s differ f r o m one a n o t h e r , a n d h o w are t h e y s i m i l a r t o one
a n o t h e r ? Is t h e r e a g e n e r a l f r a m e w o r k t h a t c a p t u r e s a l l o f these m e t h o d s ? What
s o r t o f a p p l i c a t i o n s are t h e y e s p e c i a l l y w e l l s u i t e d for?
W h i l e w e d o p r e s e n t m u c h o f t h e u n d e r l y i n g t h e o r y for R B F a n d M L S ap-
proximation methods, t h e e m p h a s i s i n t h i s b o o k is n o t o n p r o o f s . For read-
ers w h o are interested i n all the mathematical details and intricacies of the
t h e o r y we r e c o m m e n d t h e t w o e x c e l l e n t recent monographs [ B u h m a n n (2003);
W e n d l a n d (2005a)]. I n s t e a d , o u r o b j e c t i v e is t o m a k e t h e t h e o r y accessible t o a
w i d e a u d i e n c e t h a t i n c l u d e s g r a d u a t e s t u d e n t s a n d p r a c t i t i o n e r s i n a l l s o r t s o f sci-
ence a n d e n g i n e e r i n g fields. W e w a n t t o p u t t h e m a t h e m a t i c a l t h e o r y i n t h e c o n t e x t
o f a p p l i c a t i o n s a n d p r o v i d e M A T L A B i m p l e m e n t a t i o n s w h i c h g i v e t h e r e a d e r a n easy
e n t r y i n t o meshfree a p p r o x i m a t i o n m e t h o d s . T h e s k i l l e d r e a d e r s h o u l d t h e n easily
be able t o m o d i f y t h e p r o g r a m s p r o v i d e d here for h i s / h e r specific p u r p o s e s .
I n a c e r t a i n sense t h e present b o o k was i n s p i r e d b y t h e b e a u t i f u l l i t t l e b o o k [ T r e -
f e t h e n ( 2 0 0 0 ) ] . W h i l e t h e present b o o k is m u c h m o r e e x p a n s i v e ( f i l l i n g m o r e t h a n
1 2
five h u n d r e d pages w i t h f o r t y - s e v e n M A T L A B p r o g r a m s , one M a p l e p r o g r a m , one
h u n d r e d figures, m o r e t h a n f i f t y t a b l e s , a n d m o r e t h a n five h u n d r e d references), i t is
o u r a i m t o p r o v i d e t h e reader w i t h r e l a t i v e l y s i m p l e M A T L A B c o d e t h a t i l l u s t r a t e s
j u s t a b o u t e v e r y aspect discussed i n t h e b o o k .
A l l M A T L A B p r o g r a m s p r i n t e d i n t h e t e x t (as w e l l as a few m o d i f i c a t i o n s dis-
cussed) are also i n c l u d e d o n t h e enclosed C D . T h e f o l d e r MATLAB c o n t a i n s M - f i l e s
a n d d a t a files o f t y p e M A T t h a t h a v e b e e n w r i t t e n a n d t e s t e d w i t h M A T L A B 7. F o r
t h o s e readers w h o d o n o t have access t o M A T L A B 7, t h e folder MATLAB6 c o n t a i n s
versions o f these files t h a t are c o m p a t i b l e w i t h t h e o l d e r M A T L A B release. The
m a i n difference b e t w e e n t h e t w o v e r s i o n s is t h e use o f a n o n y m o u s f u n c t i o n s i n t h e
M A T L A B 7 code as c o m p a r e d t o i n l i n e f u n c t i o n s i n t h e M A T L A B 6 v e r s i o n . T w o
packages f r o m t h e M A T L A B C e n t r a l F i l e E x c h a n g e [ M C F E ] are used t h r o u g h o u t t h e
b o o k : t h e f u n c t i o n h a l t o n s e q w r i t t e n b y D a n i e l D o u g h e r t y a n d used t o g e n e r a t e
sequences o f H a l t o n p o i n t s ; t h e /cd-tree l i b r a r y ( g i v e n as a set o f M A T L A B M E X - f i l e s )
w r i t t e n b y G u y Shechter a n d used t o g e n e r a t e t h e kd-tvee data structure underlying
o u r sparse m a t r i c e s based o n c o m p a c t l y s u p p o r t e d basis f u n c t i o n s . B o t h o f these
packages are discussed i n A p p e n d i x A a n d need t o be d o w n l o a d e d separately. The
folder M a p l e o n t h e C D c o n t a i n s t h e one M a p l e file m e n t i o n e d a b o v e .
T h e m a n u s c r i p t for t h i s b o o k a n d some o f i t s e a r l i e r i n c a r n a t i o n s h a v e b e e n
used i n g r a d u a t e l e v e l courses a n d s e m i n a r s a t N o r t h w e s t e r n U n i v e r s i t y , V a n d e r b i l t
U n i v e r s i t y , a n d t h e I l l i n o i s I n s t i t u t e o f T e c h n o l o g y . S p e c i a l t h a n k s are d u e t o J o n
1
M A T L A B is a trademark of T h e MathWorks, Inc. and is used with permission. T h e Math-
Works does not warrant the accuracy of the text or exercises in this book. T h i s book's use or
discussion of M A T L A B software or related products does not constitute endorsement or sponsor-
ship by T h e MathWorks of a particular pedagogical approach or particular use of the M A T L A B
software.
2
M a p l e is a registered trademark of Waterloo Maple Inc.
Preface ix
C h e r r i e , J o h n E r i c k s o n , P a r i t o s h M o k h a s i , L a r r y S c h u m a k e r , a n d J a c k Z h a n g for
reading various p o r t i o n s of t h e m a n u s c r i p t a n d / o r M A T L A B code a n d p r o v i d i n g
h e l p f u l feedback. F i n a l l y , t h a n k s are d u e t o a l l t h e p e o p l e a t W o r l d Scientific
P u b l i s h i n g C o . w h o h e l p e d m a k e t h i s p r o j e c t a success: R a j e s h B a b u , Y i n g O i
Chiew, L i n d a K w a n , R o k T i n g Tan, and Yubing Zhai.
Greg Fasshauer
Chicago, I L , J a n u a r y 2007
Contents
Preface vii
1. Introduction 1
s
1.1 M o t i v a t i o n : Scattered D a t a I n t e r p o l a t i o n i n K 2
1.1.1 The Scattered D a t a I n t e r p o l a t i o n P r o b l e m 2
1.1.2 Example: I n t e r p o l a t i o n w i t h Distance Matrices 4
1.2 Some H i s t o r i c a l R e m a r k s 13
2. R a d i a l Basis F u n c t i o n I n t e r p o l a t i o n i n M A T L A B 17
2.1 R a d i a l (Basis) F u n c t i o n s 17
2.2 R a d i a l Basis F u n c t i o n I n t e r p o l a t i o n 19
4.1 E x a m p l e 1: G a u s s i a n s 37
4.2 E x a m p l e 2: L a g u e r r e - G a u s s i a n s 38
4.3 E x a m p l e 3: P o i s s o n R a d i a l F u n c t i o n s 39
4.4 E x a m p l e 4: M a t e r n F u n c t i o n s 41
4.5 E x a m p l e 5: G e n e r a l i z e d Inverse M u l t i q u a d r i c s 41
4.6 E x a m p l e 6: T r u n c a t e d P o w e r F u n c t i o n s 42
4.7 E x a m p l e 7: P o t e n t i a l s a n d W h i t t a k e r R a d i a l F u n c t i o n s 43
4.8 E x a m p l e 8: I n t e g r a t i o n A g a i n s t S t r i c t l y P o s i t i v e
Definite Kernels 45
xi
xii Meshfree Approximation Methods with MATLAB
4.9 Summary 45
6. Scattered D a t a I n t e r p o l a t i o n w i t h P o l y n o m i a l Precision 53
8.1 E x a m p l e 1: G e n e r a l i z e d M u l t i q u a d r i c s 67
8.2 E x a m p l e 2: R a d i a l P o w e r s 69
8.3 E x a m p l e 3: T h i n P l a t e Splines 70
10.1 C o n d i t i o n a l l y P o s i t i v e D e f i n i t e F u n c t i o n s a n d p-Norms 79
10.2 Scattered D a t a F i t t i n g o n Manifolds 83
10.3 Remarks 83
11. C o m p a c t l y S u p p o r t e d R a d i a l Basis F u n c t i o n s 85
11.1 O p e r a t o r s for R a d i a l F u n c t i o n s a n d D i m e n s i o n W a l k s 85
11.2 Wendland's Compactly Supported Functions 87
Contents xiii
12.1 A s s e m b l y o f t h e Sparse I n t e r p o l a t i o n M a t r i x 95
12.2 Numerical Experiments w i t h CSRBFs 99
17.1.2 T h e P o w e r F u n c t i o n as I n d i c a t o r for a G o o d S h a p e
Parameter 142
17.1.3 Choosing a G o o d Shape P a r a m e t e r v i a Cross V a l i d a t i o n . 1 4 6
17.1.4 The Contour-Pade A l g o r i t h m 151
17.1.5 Summary 152
17.2 Non-stationary Interpolation 153
17.3 Stationary Interpolation 155
B. l Some I m p o r t a n t C o n c e p t s f r o m M e a s u r e T h e o r y 431
B.2 A Brief S u m m a r y of Integral Transforms 432
xviii Meshfree Approximation Methods with MATLAB
B. 3 T h e S c h w a r t z Space a n d t h e G e n e r a l i z e d F o u r i e r T r a n s f o r m . . . 433
C. l M A T L A B Programs 435
C. 2 Maple Programs 440
Bibliography 451
Index 491
Chapter 1
Introduction
l
2 Meshfree Approximation Methods with MATLAB
Since m a n y o f these a p p l i c a t i o n s e i t h e r c o m e d o w n t o a f u n c t i o n a p p r o x i m a t i o n
p r o b l e m , o r i n c l u d e f u n c t i o n a p p r o x i m a t i o n as a f u n d a m e n t a l c o m p o n e n t , we w i l l
b e g i n o u r discussion w i t h a n d i n fact base a l a r g e p a r t o f t h e c o n t e n t s o f t h i s
book on the multivariate scattered data i n t e r p o l a t i o n problem.
s
1.1 Motivation: Scattered Data Interpolation in M
i . .ft
1. Introduction 3
T h e o r e m 1.1 ( M a i r h u b e r - C u r t i s ) . If Q, C R, s
s > 2 , contains an interior
point, then there exist no Haar spaces of continuous functions except for one-
dimensional ones.
4 Meshfree Approximation Methods with MATLAB
I n o r d e r t o u n d e r s t a n d t h i s t h e o r e m w e need
N o t e t h a t existence o f a H a a r space g u a r a n t e e s i n v e r t i b i l i t y o f t h e i n t e r p o l a t i o n
m a t r i x A, i.e., existence a n d uniqueness o f a n i n t e r p o l a n t o f t h e f o r m (1.1) t o
d a t a specified a t XI,...,XN f r o m t h e space B. A s m e n t i o n e d above, u n i v a r i a t e
p o l y n o m i a l s o f degree N 1 f o r m a n A ^ - d i m e n s i o n a l H a a r space for d a t a g i v e n a t
X i , . . . , X N .
T h e M a i r h u b e r - C u r t i s t h e o r e m t e l l s us t h a t i f w e w a n t t o have a w e l l - p o s e d
m u l t i v a r i a t e scattered d a t a i n t e r p o l a t i o n p r o b l e m we can no longer fix i n advance
t h e set o f basis f u n c t i o n s w e p l a n t o use for i n t e r p o l a t i o n o f a r b i t r a r y s c a t t e r e d d a t a .
For e x a m p l e , i t is n o t possible t o p e r f o r m u n i q u e i n t e r p o l a t i o n w i t h ( m u l t i v a r i a t e )
2
p o l y n o m i a l s o f degree N t o d a t a g i v e n a t a r b i t r a r y l o c a t i o n s i n M . Instead, the
basis s h o u l d d e p e n d o n t h e d a t a l o c a t i o n s . W e w i l l g i v e a s i m p l e e x a m p l e o f s u c h
a n i n t e r p o l a t i o n scheme i n t h e n e x t s u b s e c t i o n .
f (x)
a
s
= 4 ]^[a:d(l - x d ) , X = ( X ! , . . . , x a ) e [o, i ] .s
d=l
1. Introduction 5
S
T h i s f u n c t i o n is zero o n t h e b o u n d a r y o f t h e u n i t c u b e i n R a n d has a m a x i m u m
v a l u e o f one at t h e center o f t h e c u b e . A s i m p l e M A T L A B s c r i p t d e f i n i n g f s is g i v e n
as P r o g r a m C . l i n A p p e n d i x C .
We w i l l use a set o f u n i f o r m l y s c a t t e r e d d a t a sites i n t h e u n i t c u b e a t w h i c h
w e sample o u r t e s t f u n c t i o n f . T h i s w i l l be a c c o m p l i s h e d here ( a n d i n m a n y o t h e r
s
2
Fig. 1.1 289 Halton points in the unit square in R .
A s e x p l a i n e d i n t h e p r e v i o u s s u b s e c t i o n we are i n t e r e s t e d i n c o n s t r u c t i n g a ( c o n -
t i n u o u s ) f u n c t i o n Vf t h a t interpolates t h e samples o b t a i n e d f r o m f s a t t h e set o f
H a l t o n p o i n t s , i.e., s u c h t h a t
Vf(xj) f (xj),
s Xj a Halton point.
A s p o i n t e d o u t above, i f s = 1, t h e n t h i s p r o b l e m is o f t e n s o l v e d u s i n g u n i v a r i a t e
p o l y n o m i a l s or splines. F o r a s m a l l n u m b e r o f d a t a sites p o l y n o m i a l s m a y w o r k
s a t i s f a c t o r i l y . H o w e v e r , i f t h e n u m b e r o f p o i n t s increases, i.e., t h e p o l y n o m i a l degree
g r o w s , t h e n i t is w e l l k n o w n t h a t one s h o u l d use splines ( o r piecewise p o l y n o m i a l s )
t o a v o i d o s c i l l a t i o n s . T h e s i m p l e s t s o l u t i o n is t o use a c o n t i n u o u s piecewise l i n e a r
spline, i.e., t o " c o n n e c t t h e d o t s " . I t is also w e l l k n o w n t h a t one possible basis f o r
t h e space o f piecewise l i n e a r splines i n t e r p o l a t i n g d a t a at a g i v e n set o f p o i n t s i n
[0,1] consists o f t h e shifts o f t h e a b s o l u t e v a l u e f u n c t i o n t o t h e d a t a sites. I n o t h e r
w o r d s , we c a n c o n s t r u c t t h e piecewise l i n e a r s p l i n e i n t e r p o l a n t b y a s s u m i n g Vf is
of the form
N
V (x)
f = ^Tc \x k - x \,
k are [0,1],
fc=i
6 Meshfree Approximation Methods with MATLAB
V {xj)
f = h{ ), Xj j = l,...,N.
(1.3)
\X N - ml \X N - X2\ \x N - XN -CN _
Vf( ) x
= Y2 W Ck x
~ x
k h , xe[0,iy (1.4)
fc=i
1. Introduction 7
Fig. 1.2 A typical basis function for the Euclidean distance matrix centered at the origin in R .
Program 1.1. D i s t a n c e M a t r i x . m
% DM = DistanceMatrix(dsites,ctrs)
% Forms t h e d i s t a n c e m a t r i x of two s e t s o f p o i n t s in R"s,
% i.e., DM(i,j) = II datasite_i - center_j I I_2.
% Input
% dsites: Mxs m a t r i x r e p r e s e n t i n g a set of M data s i t e s in R~s
7, (i.e., each row c o n t a i n s one s-dimensional point)
7o ctrs: Nxs m a t r i x r e p r e s e n t i n g a set of N centers in R~s
7o (one center per row)
7. O u t p u t
7. DM: MxN m a t r i x w h o s e i , j p o s i t i o n c o n t a i n s the Euclidean
7o d i s t a n c e between the i - t h data site and j - t h c e n t e r
8 Meshfree Approximation Methods with MATLAB
1 f u n c t i o n DM = D i s t a n c e M a t r i x ( d s i t e s , c t r s )
2 [M,s] = s i z e ( d s i t e s ) ; [N,s] = size(ctrs);
3 DM = zeros(M,N);
% Accumulate sum of s q u a r e s of c o o r d i n a t e differences
% The n d g r i d command produces two MxN matrices:
'/ . d r , c o n s i s t i n g of N i d e n t i c a l columns ( e a c h containing
% the d-th coordinate of t h e M d a t a sites)
7o c c , c o n s i s t i n g of M i d e n t i c a l rows ( e a c h containing
% the d-th coordinate of t h e N c e n t e r s )
4 f o r d=l:s
5 [[Link]] = n d g r i d ( d s i t e s ( : , d ) , c t r s ( : , d ) ) ;
6 DM = DM + ( d r - c c ) . " 2 ;
7 end
8 DM = sqrt(DM);
N o t e t h a t t h i s s u b r o u t i n e c a n easily b e m o d i f i e d t o p r o d u c e a p - n o r m d i s t a n c e
m a t r i x b y m a k i n g t h e o b v i o u s changes t o lines 6 a n d 8 o f t h e code i n P r o g r a m 1.1.
W e w i l l come back t o t h i s idea i n C h a p t e r 10.
O u r first m a i n s c r i p t is P r o g r a m 1.2. T h i s s c r i p t c a n be u s e d t o c o m p u t e the
distance m a t r i x interpolant to d a t a sampled f r o m the test function / s provided by
P r o g r a m C l . W e use H a l t o n p o i n t s a n d are a b l e t o select t h e space d i m e n s i o n
s and number of points N b y e d i t i n g lines 1 a n d 2 o f t h e code. The subrou-
t i n e MakeSDGrid.m w h i c h w e use t o c o m p u t e t h e e q u a l l y spaced p o i n t s i n t h e s-
dimensional u n i t cube o n line 6 o f D i s t a n c e M a t r i x F i t .mis p r o v i d e d i n A p p e n d i x C.
T h e s e e q u a l l y spaced p o i n t s are used as e v a l u a t i o n p o i n t s a n d t o c o m p u t e errors.
N o t e t h a t since t h e d i s t a n c e m a t r i x i n t e r p o l a n t is o f t h e f o r m ( 1 . 4 ) i t s s i m u l t a n e o u s
e v a l u a t i o n a t t h e e n t i r e set o f e v a l u a t i o n p o i n t s a m o u n t s t o a m a t r i x - v e c t o r p r o d u c t
o f t h e e v a l u a t i o n m a t r i x EM a n d t h e coefficients c. H e r e t h e e v a l u a t i o n m a t r i x has
t h e same s t r u c t u r e as t h e i n t e r p o l a t i o n m a t r i x a n d c a n also be c o m p u t e d u s i n g t h e
subroutine Distancematrix.m ( o n l y u s i n g e v a l u a t i o n p o i n t s i n place o f t h e d a t a
sites, see l i n e 9 o f D i s t a n c e M a t r i x F i t .m). T h e coefficient v e c t o r c is s u p p l i e d d i -
r e c t l y as s o l u t i o n o f t h e l i n e a r s y s t e m Ac = f (see ( 1 . 3 ) a n d t h e M A T L A B e x p r e s s i o n
IM\rhs o n l i n e 1 0 o f t h e p r o g r a m ) . T h e e v a l u a t i o n p o i n t s are s u b s e q u e n t l y used
for t h e e r r o r c o m p u t a t i o n i n lines 1 1 - 1 3 a n d are also u s e d for p l o t t i n g p u r p o s e s i n
t h e last p a r t o f t h e p r o g r a m (lines 1 6 - 3 5 ) . N o t e t h a t for t h i s e x a m p l e w e k n o w t h e
function f s t h a t g e n e r a t e d t h e d a t a , a n d t h e r e f o r e are a b l e t o c o m p u t e t h e e r r o r i n
o u r r e c o n s t r u c t i o n . T h e s u b r o u t i n e s t h a t p r o d u c e t h e 2 D a n d 3 D p l o t s o n lines 2 4
3 2 are p r o v i d e d i n A p p e n d i x C . N o t e t h a t t h e use o f r e s h a p e o n lines 2 2 - 2 3 and
2 7 - 2 9 c o r r e s p o n d s t o t h e use o f meshgrid for p l o t t i n g p u r p o s e s .
P r o g r a m 1.2. D i s t a n c e M a t r i x F i t .m
% DistanceMatrixFit
% S c r i p t that uses E u c l i d e a n d i s t a n c e matrices t o perform
1. Introduction
% s c a t t e r e d d a t a i n t e r p o l a t i o n f o r a r b i t r a r y space dimensions
% C a l l s on: D i s t a n c e M a t r i x , MakeSDGrid, t e s t f u n c t i o n
/ Uses:
0 h a l t o n s e q ( w r i t t e n by D a n i e l Dougherty from MATLAB
% C e n t r a l F i l e Exchange)
1 s = 3;
2 k = 2; N = ( 2 ~ k + l ) ~ s ;
3 n e v a l = 10; M = n e v a l ' s ;
% Use Halton p o i n t s a s d a t a s i t e s and c e n t e r s
4 d s i t e s = haltonseq(N,s);
5 ctrs = dsites;
% Create n e v a l ~ s e q u a l l y spaced e v a l u a t i o n l o c a t i o n s i n t h e
% s-dimensional u n i t cube
6 e p o i n t s = MakeSDGrid(s,neval);
% Create right-hand side vector,
% i . e . , evaluate the t e s t function a t the data s i t e s
7 rhs= testfunction(s,dsites);
% Compute d i s t a n c e m a t r i x f o r t h e d a t a s i t e s and c e n t e r s
8 IM = D i s t a n c e M a t r i x ( d s i t e s , c t r s ) ;
% Compute d i s t a n c e m a t r i x f o r e v a l u a t i o n p o i n t s and c e n t e r s
9 EM = D i s t a n c e M a t r i x ( e p o i n t s , c t r s ) ;
% E v a l u a t e t h e i n t e r p o l a n t on e v a l u a t i o n p o i n t s
/ ( e v a l u a t i o n m a t r i x * s o l u t i o n of i n t e r p o l a t i o n system)
0
10 Pf = EM * ( I M \ r h s ) ;
% Compute exact s o l u t i o n ,
% i . e . , e v a l u a t e t e s t f u n c t i o n on e v a l u a t i o n p o i n t s
11 exact = t e s t f u n c t i o n ( s , e p o i n t s ) ;
% Compute maximum and RMS e r r o r s on e v a l u a t i o n g r i d
12 maxerr = n o r m ( P f - e x a c t , i n f ) ;
13 [Link] = n o r m ( P f - e x a c t ) / s q r t ( M ) ;
J
14 f p r i n t f ( R M S e r r o r : 7 e\n', r m s _ e r r )
0
0
15 f p r i n t f ('Maximum e r r o r : / e\n', maxerr)
0
16 s w i t c h s
17 case 1
18 plot(epoints, Pf)
19 figure; plot(epoints, abs(Pf-exact))
20 case 2
21 fview = [-30,30];
22 xe = r e s h a p e ( e p o i n t s ( : , 2 ) , n e v a l , n e v a l ) ;
23 ye = r e s h a p e ( e p o i n t s ( : , 1 ) , n e v a l , n e v a l ) ;
24 PlotSurf(xe,ye,Pf,neval,[Link],fview);
25 PlotError2D(xe,ye,Pf,exact,maxerr,neval,fview);
26 case 3
10 Meshfree Approximation Methods with M A T L A B
27 xe = reshapeCepoints(:,2),neval,neval,neval);
38 ye = reshape(epoints(:,1),neval,neval,neval);
29 ze = reshape(epoints(:,3),neval,neval,neval);
30 xslice = .25:.25:1; yslice = 1; zslice = [0,0.5];
31 PlotSIices(xe,ye,ze,Pf,neval,xslice,yslice,zslice);
32a PlotErrorSlices(xe,ye,ze,Pf,exact,neval,...
32b xslice,yslice,zslice);
33 otherwise
34 disp('Cannot display plots for s>3')
35 end
RMS-error (1.5)
s
Table 1.1 Distance matrix fit to N Halton points in [0, l ] , s 1, 2, 3.
ID 2D 3D
X X
Fig. 1.3 F i t (left) and absolute error (right) for 5 point distance matrix interpolation in I D .
y 0 0 x y 0 0
Fig. 1.4 F i t s (left) and errors (right) for distance matrix interpolation with 289 points in 2D
(top), and 729 points in 3D (bottom).
s
Table 1.2 Distance matrix fit to N Halton points in [0, l ] , s 4, 5, 6.
4D 5D 6D
metric. Second, as t h e M A T L A B s c r i p t s s h o w , t h e m e t h o d is e x t r e m e l y s i m p l e t o
i m p l e m e n t for a n y space d i m e n s i o n s. For example, no u n d e r l y i n g c o m p u t a t i o n a l
m e s h is r e q u i r e d t o c o m p u t e t h e i n t e r p o l a n t . T h e process o f m e s h g e n e r a t i o n is
a m a j o r f a c t o r w h e n w o r k i n g i n h i g h e r space d i m e n s i o n s w i t h polynomial-based
m e t h o d s such as splines o r finite e l e m e n t s . A l l t h a t is r e q u i r e d for o u r m e t h o d is
t h e p a i r w i s e d i s t a n c e b e t w e e n t h e d a t a sites. T h e r e f o r e , w e h a v e w h a t is k n o w n as
a meshfree (or meshless) method.
T h i r d , t h e a c c u r a c y o f t h e m e t h o d i m p r o v e s i f w e a d d m o r e d a t a sites. I n f a c t ,
i t seems t h a t t h e R M S - e r r o r i n T a b l e s 1.1 a n d 1.2 is r e d u c e d b y a f a c t o r o f a b o u t
s
t w o f r o m one r o w t o t h e n e x t . Since w e use (2* + l ) uniformly distributed random
1. Introduction 13
D o n a l d S h e p a r d , w h o as a n u n d e r g r a d u a t e s t u d e n t a t H a r v a r d U n i v e r s i t y ,
suggested t h e use o f w h a t are n o w c a l l e d Shepard functions in the late
1960s (see C h a p t e r 22). T h e p u b l i c a t i o n [Shepard (1968)] discusses t h e
basic inverse d i s t a n c e w e i g h t e d S h e p a r d m e t h o d a n d some m o d i f i c a t i o n s
thereof. T h e m e t h o d was a t t h e t i m e i n c o r p o r a t e d i n t o a c o m p u t e r p r o -
g r a m , S Y M A P , for m a p m a k i n g .
R o l l a n d H a r d y , w h o was a geodesist at I o w a S t a t e U n i v e r s i t y . H e i n t r o -
d u c e d t h e so-called multiquadrics ( M Q s ) i n t h e e a r l y 1970s (see, e.g., [ H a r d y
(1971)] or C h a p t e r 8 ) . H a r d y ' s w o r k was p r i m a r i l y c o n c e r n e d w i t h a p p l i -
c a t i o n s i n geodesy a n d m a p p i n g .
Meshfree Approximation Methods with MATLAB
R o b e r t L . H a r d e r a n d R o b e r t N . D e s m a r a i s , w h o w e r e aerospace engineers
at M a c N e a l - S c h w e n d l e r C o r p o r a t i o n ( M S C S o f t w a r e ) , a n d N A S A ' s L a n g l e y
Research C e n t e r . T h e y i n t r o d u c e d t h e s o - c a l l e d thin plate splines (TPSs)
i n 1972 (see, e.g., [ H a r d e r a n d D e s m a r a i s ( 1 9 7 2 ) ] o r C h a p t e r 8 ) . T h e i r w o r k
was c o n c e r n e d m o s t l y w i t h a i r c r a f t d e s i g n .
Jean Duchon, a m a t h e m a t i c i a n at the Universite Joseph Fourier i n Greno-
ble, France. D u c h o n suggested a v a r i a t i o n a l a p p r o a c h minimizing the
2 2
integral of V / in R w h i c h also leads t o t h e t h i n p l a t e splines. This
w o r k was d o n e i n t h e m i d 1970s a n d is c o n s i d e r e d t o b e t h e foundation
o f t h e v a r i a t i o n a l a p p r o a c h t o r a d i a l basis f u n c t i o n s (see [ D u c h o n ( 1 9 7 6 ) ;
D u c h o n (1977); D u c h o n (1978); D u c h o n (1980)]) or C h a p t e r 13).
Jean Meinguet, a mathematican a t U n i v e r s i t e C a t h o l i q u e de L o u v a i n i n
L o u v a i n , B e l g i u m . M e i n g u e t i n t r o d u c e d w h a t he c a l l e d surface splines in
t h e l a t e 1970s. Surface splines a n d t h i n p l a t e splines f a l l u n d e r w h a t w e
w i l l refer t o as polyharmonic splines (see, e.g., [ M e i n g u e t ( 1 9 7 9 a ) ; M e i n g u e t
(1979b); M e i n g u e t (1979c); M e i n g u e t (1984)] or C h a p t e r 8).
P e t e r L a n c a s t e r a n d K e s Salkauskas, m a t h e m a t i c i a n s at t h e U n i v e r s i t y o f
Calgary, Canada. T h e y p u b l i s h e d [ L a n c a s t e r a n d Salkauskas ( 1 9 8 1 ) ] i n -
t r o d u c i n g t h e moving least squares method (a g e n e r a l i z a t i o n o f S h e p a r d
functions).
R i c h a r d Franke, a m a t h e m a t i c i a n at the N a v a l Postgraduate School i n M o n -
terey, C a l i f o r n i a . I n [ F r a n k e (1982a)] he c o m p a r e d v a r i o u s s c a t t e r e d d a t a
i n t e r p o l a t i o n methods, a n d concluded M Q s a n d T P S s were the best. Franke
also c o n j e c t u r e d t h a t t h e i n t e r p o l a t i o n m a t r i x for M Q s is i n v e r t i b l e .
W o l o d y m y r ( W a l l y ) M a d y c h , a m a t h e m a t i c i a n at the U n i v e r s i t y of C o n -
necticut, and Stuart A l a n Nelson, a m a t h e m a t i c i a n from Iowa State Univer-
sity. I n 1983 t h e y c o m p l e t e d t h e i r m a n u s c r i p t [ M a d y c h a n d N e l s o n (1983)]
i n w h i c h t h e y p r o v e d Franke's conjecture ( a n d m u c h more) based o n a varia-
t i o n a l approach. However, t h i s m a n u s c r i p t was never published. O t h e r fun-
d a m e n t a l p a p e r s b y these t w o a u t h o r s are, e.g., [ M a d y c h a n d N e l s o n ( 1 9 8 8 ) ;
M a d y c h and Nelson (1990a); M a d y c h and Nelson (1992)].
Charles Micchelli, a m a t h e m a t i c i a n at the I B M W a t s o n Research Center.
M i c c h e l l i p u b l i s h e d t h e p a p e r [ M i c c h e l l i ( 1 9 8 6 ) ] . H e also p r o v e d F r a n k e ' s
conjecture. H i s p r o o f s are r o o t e d i n t h e w o r k o f [ B o c h n e r ( 1 9 3 2 ) ; B o c h n e r
(1933)] a n d [Schoenberg ( 1 9 3 7 ) ; S c h o e n b e r g ( 1 9 3 8 a ) ; S c h o e n b e r g ( 1 9 3 8 b ) ]
on positive definite a n d c o m p l e t e l y m o n o t o n e functions. T h i s is also t h e
approach we w i l l follow t h r o u g h o u t m u c h of this book.
G r a c e W a h b a , a s t a t i s t i c i a n a t t h e U n i v e r s i t y o f W i s c o n s i n . She s t u d i e d t h e
use o f t h i n p l a t e splines for s t a t i s t i c a l p u r p o s e s i n t h e c o n t e x t o f s m o o t h i n g
noisy data a n d d a t a o n spheres, a n d i n t r o d u c e d t h e A N O V A a n d cross
v a l i d a t i o n a p p r o a c h e s t o t h e r a d i a l basis f u n c t i o n s e t t i n g ( s e e , e.g., [Wahba
1. Introduction 15
(1979); W a h b a ( 1 9 8 1 ) ; W a h b a a n d W e n d e l b e r g e r ( 1 9 8 0 ) ] ) . O n e o f t h e first
m o n o g r a p h s o n t h e s u b j e c t is [ W a h b a ( 1 9 9 0 b ) ] .
R o b e r t Schaback, a m a t h e m a t i c i a n a t t h e U n i v e r s i t y o f G o t t i n g e n , Ger-
many. Compactly supported radial basis functions ( C S R B F s ) were i n t r o -
d u c e d i n [Schaback ( 1 9 9 5 a ) ] , a n d a v e r y p o p u l a r f a m i l y o f C S R B F s was
presented b y H o l g e r W e n d l a n d (also a m a t h e m a t i c i a n i n G o t t i n g e n ) i n his
P h . D . thesis (see also [ W e n d l a n d (1995)] a n d C h a p t e r 1 1 ) . B o t h o f these
a u t h o r s have c o n t r i b u t e d e x t e n s i v e l y t o t h e field o f r a d i a l basis f u n c t i o n s .
W e m e n t i o n p a r t i c u l a r l y t h e recent m o n o g r a p h [ W e n d l a n d ( 2 0 0 5 a ) ] .
Chapter 2
A s a first e x a m p l e we p i c k a f u n c t i o n w e l l - r e p r e s e n t e d i n m a n y branches o f m a t h e -
matics, namely the Gaussian
2
^( )= -(^) ,
r e r e i
2
O u r shape parameter e is r e l a t e d t o t h e v a r i a n c e a of the normal distribution
2 2
function by e = l/(2o~ ). I f we c o m p o s e t h e G a u s s i a n w i t h t h e E u c l i d e a n d i s t a n c e
s
f u n c t i o n 11 - j12 we o b t a i n for a n y f i x e d center x k GR a multivariate function
3
$ (sc) = e - " -
fc
e a a ! a ? f c
ll2, x e R .
^fc(x) = (p(\\x - x \\ ).
k 2
S
D e f i n i t i o n 2.1. A f u n c t i o n $ : R > R is c a l l e d radial p r o v i d e d t h e r e exists a
univariate f u n c t i o n <p : [0, 00) > R s u c h t h a t
a n d || || is some n o r m o n M u s u a l l y t h e E u c l i d e a n n o r m .
s
17
18 Meshfree Approximation Methods with MATLAB
u is
2
Fig. 2.1 Gaussian with e = 1 (left) and e = 3 (right) centered at the origin in R .
D e f i n i t i o n 2.1 a n d t h e d i s c u s s i o n l e a d i n g u p t o i t s h o w a g a i n w h y i t m a k e s sense
t o c a l l ip t h e basic f u n c t i o n , a n d &k(\\ H2) (centered at x ) k a r a d i a l basis f u n c t i o n .
O n e single basic f u n c t i o n generates a l l o f t h e basis f u n c t i o n s t h a t are u s e d i n t h e
expansion (1.1).
R a d i a l f u n c t i o n i n t e r p o l a n t s h a v e t h e nice p r o p e r t y o f b e i n g i n v a r i a n t u n d e r a l l
Euclidean transformations (i.e., translations, rotations, a n d reflections). B y this
w e m e a n t h a t i t does n o t m a t t e r w h e t h e r w e first c o m p u t e t h e R B F i n t e r p o l a n t
a n d t h e n a p p l y a E u c l i d e a n t r a n s f o r m a t i o n , o r i f w e first t r a n s f o r m t h e d a t a a n d
t h e n c o m p u t e t h e i n t e r p o l a n t . T h i s is a n i m m e d i a t e consequence o f t h e fact t h a t
E u c l i d e a n t r a n s f o r m a t i o n s are c h a r a c t e r i z e d by orthogonal transformation matri-
ces a n d are t h e r e f o r e 2 - n o r m - i n v a r i a n t . I n v a r i a n c e u n d e r t r a n s l a t i o n , r o t a t i o n a n d
r e f l e c t i o n is o f t e n d e s i r a b l e i n a p p l i c a t i o n s .
Moreover, the application of radial functions to the solution of the scattered data
i n t e r p o l a t i o n p r o b l e m (as w e l l as m a n y o t h e r m u l t i v a r i a t e a p p r o x i m a t i o n p r o b l e m s )
benefits f r o m t h e fact t h a t t h e i n t e r p o l a t i o n p r o b l e m b e c o m e s i n s e n s i t i v e t o the
2. Radial Basis Function Interpolation in M A T L A B 19
d i m e n s i o n s o f t h e space i n w h i c h t h e d a t a sites l i e . I n s t e a d o f h a v i n g t o d e a l w i t h
a multivariate function $ (whose c o m p l e x i t y w i l l increase w i t h i n c r e a s i n g space
d i m e n s i o n s) w e c a n w o r k w i t h t h e same u n i v a r i a t e f u n c t i o n ip for a l l choices o f s.
I n s t e a d o f u s i n g s i m p l e d i s t a n c e m a t r i c e s as w e d i d earlier, w e n o w use a r a d i a l
s
basis f u n c t i o n e x p a n s i o n t o solve t h e s c a t t e r e d d a t a i n t e r p o l a t i o n p r o b l e m i n M. b y
assuming
N
s
V (x)
f = J2ck<p(\\x-x \\ ), k 2 x e R . (2.1)
fc=i
T h e coefficients c k are f o u n d b y e n f o r c i n g t h e i n t e r p o l a t i o n c o n d i t i o n s , a n d t h u s
solving the linear system
A s t h e s o l u t i o n o f t h e s c a t t e r e d d a t a i n t e r p o l a t i o n p r o b l e m hinges e n t i r e l y o n t h e
solution of this system of linear equations we w i l l devote t h e next chapter t o the
q u e s t i o n o f w h e n (i.e., for w h a t t y p e o f basic f u n c t i o n s <p) t h e s y s t e m m a t r i x is
non-singular.
F o r t h e n u m e r i c a l e x a m p l e p r e s e n t e d b e l o w w e r e s t r i c t ourselves t o t h e t w o -
d i m e n s i o n a l case s = 2. A s basic f u n c t i o n ip w e w i l l use b o t h Gaussians a n d t h e
l i n e a r f u n c t i o n <p(r) = r w h i c h gives rise t o t h e E u c l i d e a n d i s t a n c e m a t r i x a p p r o a c h
used earlier.
T h e code o f t h e MATLAB script RBFInterpolation2D.m (see P r o g r a m 2.1)
w e use t o p e r f o r m R B F i n t e r p o l a t i o n i n 2 D is v e r y s i m i l a r t o t h e earlier s c r i p t
D i s t a n c e M a t r i x F i t .m. I t also m a k e s use o f t h e s u b r o u t i n e D i s t a n c e M a t r i x . m .
W h i l e i t is easy t o w r i t e a v e r s i o n o f t h e i n t e r p o l a t i o n s c r i p t t h a t w o r k s for a n y
space d i m e n s i o n s ( j u s t as w e d i d i n D i s t a n c e M a t r i x F i t .m) w e w i l l s t i c k w i t h a
basic 2 D v e r s i o n here.
I n l i n e 1 we define t h e G a u s s i a n R B F as a M A T L A B a n o n y m o u s f u n c t i o n t h a t
accepts a m a t r i x a r g u m e n t ( n a m e l y t h e o u t p u t f r o m D i s t a n c e M a t r i x ) a l o n g w i t h
its shape p a r a m e t e r . N o t e t h a t t h i s f e a t u r e is o n l y a v a i l a b l e since M A T L A B Re-
lease 7. For older M A T L A B versions w e suggest u s i n g a n i n l i n e f u n c t i o n i n s t e a d
(see t h e p r o g r a m s i n t h e f o l d e r M a t l a b 6 o f t h e enclosed C D ) . I f e x e c u t i o n speed
is i m p o r t a n t , t h e n one s h o u l d e x p l i c i t l y p r o v i d e t h e f u n c t i o n ( e i t h e r h a r d c o d e d d i -
r e c t l y w h e r e needed, o r as a n M - f i l e ) . T h i s l a t t e r a p p r o a c h w i l l a l w a y s be m o r e
efficient t h a n t h e i n l i n e o r e v e n a n o n y m o u s f u n c t i o n a p p r o a c h . H o w e v e r , t h e n t h e
i n t e r p o l a t i o n p r o g r a m is n o l o n g e r as generic.
20 Meshfree Approximation Methods with MATLAB
W e c a n replace t h e d e f i n i t i o n o f t h e G a u s s i a n o n l i n e 1 b y t h e d e f i n i t i o n o f t h e
l i n e a r f u n c t i o n <p(r) = r o r a n y o t h e r a d m i s s i b l e R B F w e w i l l e n c o u n t e r l a t e r . In
lines 2 - 6 w e define a t e s t f u n c t i o n t h a t w e w i l l s a m p l e s i m i l a r l y t o t h e function
f
s used i n t h e i n t r o d u c t o r y e x a m p l e . H e r e ( a n d i n m a n y l a t e r e x a m p l e s ) w e use
Franke's function
ffay) = ^ - l / 4 ( ( 9 . r - 2 ) + (9s/-2) )
e
2 2
+ 3 _(l/49)(9a +l) -(l/10)(9 /+l)
e ;
2
!
2
2 2 2 2
+ I - l / 4 ( ( 9 x - 7 ) + (9j/-3) ) _ l _ (
e e 9 a ; -4) -(9y-7)
2 5
w h i c h is a s t a n d a r d t e s t f u n c t i o n f o r 2 D s c a t t e r e d d a t a fitting. N o t e t h a t we used
2
(x,y) t o d e n o t e t h e t w o c o m p o n e n t s o f a; G M . T h e g r a p h o f Franke's function
over t h e u n i t s q u a r e is s h o w n i n F i g u r e 2.2.
f c 2
{(2 + l) } = { 9 , 2 5 , 8 1 , 2 8 9 , 1 0 8 9 , 4 2 2 5 , . . . } . T h e characters u or h ( i n place o f
/ s) are used t o d e n o t e e i t h e r u n i f o r m l y s p a c e d p o i n t s , o r H a l t o n p o i n t s i n t h e u n i t
0
s q u a r e . T h e set o f d a t a p o i n t s is d e f i n e d a n d l o a d e d i n lines 7 a n d 8. A s i n t h e e a r l i e r
e x a m p l e w e consider here o n l y t h e case w h e r e t h e centers f o r t h e R B F s c o i n c i d e
w i t h the d a t a locations (line 9 ) .
A g r i d o f e v a l u a t i o n p o i n t s used t o e v a l u a t e o u r i n t e r p o l a n t for t h e p u r p o s e s
o f r e n d e r i n g a n d e r r o r c o m p u t a t i o n is d e f i n e d i n l i n e s 10 a n d 1 1 . T h e t e s t data
( r i g h t - h a n d side o f t h e i n t e r p o l a t i o n e q u a t i o n s ) are c o m p u t e d o n l i n e 12 w h e r e t h e
t e s t f u n c t i o n is s a m p l e d a t t h e d a t a sites.
T h e m a i n p a r t o f t h e code is g i v e n b y lines 1 3 - 1 7 . N o t e t h a t t h i s p a r t is v e r y
similar t o t h e corresponding segment ( l i n e s 7 - 9 ) i n D i s t a n c e M a t r i x F i t .m. The
o n l y difference is t h a t w e n o w a p p l y t h e basic f u n c t i o n <p t o t h e e n t i r e distance
matrices i n order to o b t a i n the i n t e r p o l a t i o n a n d evaluation matrices.
2. Radial Basis Function Interpolation in MATLAB 21
P r o g r a m 2.1. RBFInterpolation2D.m
% RBFInterpolation2D
S c r i p t t h a t performs b a s i c 2D RBF i n t e r p o l a t i o n
% C a l l s on: D i s t a n c e M a t r i x
/ Define t h e G a u s s i a n RBF and shape parameter
1 r b f = @(e,r) e x p ( - ( e * r ) . ~ 2 ) ; ep = 21.1;
% Define Franke's f u n c t i o n as t e s t f u n c t i o n
2 f l = @(x,y) 0 . 7 5 * e x p ( - ( ( 9 * x - 2 ) . ~ 2 + ( 9 * y - 2 ) . " 2 ) / 4 ) ;
3 f 2 = <3(x,y) 0 . 7 5 * e x p ( - ( ( 9 * x + l ) . ~ 2 / 4 9 + ( 9 * y + l ) . ~ 2 / 1 0 ) ) ;
4 f 3 = @(x,y) 0 . 5 * e x p ( - ( ( 9 * x - 7 ) . ~ 2 + ( 9 * y - 3 ) . ~ 2 ) / 4 ) ;
5 f 4 = @(x,y) 0 . 2 * e x p ( - ( ( 9 * x - 4 ) . ~ 2 + ( 9 * y - 7 ) . " 2 ) ) ;
6 t e s t f u n c t i o n = @(x,y) f l ( x , y ) + f 2 ( x , y ) + f 3 ( x , y ) - f 4 ( x , y ) ;
7 N = 1089; g r i d t y p e = 'h';
% Load d a t a p o i n t s
0
8 name = s p r i n t f ( 'Data2D_y d / s' ,N,gridtype) ; load(name)
o 0
9 ctrs = dsites;
10 n e v a l = 40; g r i d = l i n s p a c e ( 0 , 1 , n e v a l ) ;
11 [xe,ye] = m e s h g r i d ( g r i d ) ; e p o i n t s = [ x e ( : ) y e ( : ) ] ;
% Evaluate the t e s t function a t the data points
12 r h s = t e s t f u n c t i o n ( d s i t e s ( : , 1 ) , d s i t e s ( : , 2 ) ) ;
% Compute d i s t a n c e m a t r i x between t h e d a t a s i t e s and c e n t e r s
13 DM_data = D i s t a n c e M a t r i x ( d s i t e s , c t r s ) ;
% Compute i n t e r p o l a t i o n m a t r i x
14 IM = rbf(ep,DM_data);
% Compute d i s t a n c e m a t r i x between e v a l u a t i o n p o i n t s and c e n t e r s
15 DM_eval = D i s t a n c e M a t r i x ( e p o i n t s , c t r s ) ;
% Compute e v a l u a t i o n m a t r i x
16 EM = rbf(ep,DM_eval);
% Compute RBF i n t e r p o l a n t
% ( e v a l u a t i o n m a t r i x * s o l u t i o n of i n t e r p o l a t i o n system)
17 Pf = EM * ( I M \ r h s ) ;
% Compute exact s o l u t i o n , i . e . ,
% e v a l u a t e t e s t f u n c t i o n on e v a l u a t i o n p o i n t s
18 exact = t e s t f u n c t i o n ( e p o i n t s ( : , 1 ) , e p o i n t s ( : , 2 ) ) ;
0/ Compute e r r o r s on e v a l u a t i o n g r i d
19 maxerr = n o r m ( P f - e x a c t , i n f ) ;
20 [Link] = n o r m ( P f - e x a c t ) / n e v a l ;
21 f p r i n t f ( ' R M S e r r o r : /,e\n', r m s . e r r )
22 f p r i n t f ('Maximum e r r o r : /e\n', maxerr)
0
h = h ,n x s u p m i n ||aj (2.3)
a n d i t i n d i c a t e s h o w w e l l t h e d a t a i n t h e set X fill o u t t h e d o m a i n Q. A g e o m e t r i c
i n t e r p r e t a t i o n o f t h e fill d i s t a n c e is g i v e n b y t h e r a d i u s o f t h e largest possible e m p t y
b a l l t h a t c a n be p l a c e d a m o n g t h e d a t a l o c a t i o n s i n s i d e fl (see F i g u r e 2 . 3 ) . Some-
t i m e s t h e s y n o n y m covering radius is used. I n o u r M A T L A B code w e c a n estimate
t h e fill d i s t a n c e v i a
hX = m a x C m i n C D M - e v a l ' ) ) (2.4)
w h e r e DM_eval is t h e m a t r i x c o n s i s t i n g o f p a i r w i s e d i s t a n c e s b e t w e e n t h e e v a l u a t i o n
p o i n t s ( p l a c e d o n a fine u n i f o r m g r i d i n Q) a n d t h e d a t a sites X (c.f. l i n e 15 o f
P r o g r a m 2 . 1 ) . N o t e t h a t we t r a n s p o s e t h e n o n - s y m m e t r i c e v a l u a t i o n m a t r i x . T h i s
corresponds t o finding for each e v a l u a t i o n p o i n t t h e d i s t a n c e t o t h e c o r r e -
s p o n d i n g closest d a t a site, a n d t h e n s e t t i n g hx,n as t h e w o r s t o f t h o s e d i s t a n c e s .
F i g u r e 2.3 i l l u s t r a t e s t h e fill d i s t a n c e for a set o f 25 H a l t o n p o i n t s . Note that in
t h i s case t h e largest "hole" i n t h e d a t a is n e a r t h e b o u n d a r y .
0.8
0.6
0.4
0.2
Fig. 2.4 Gaussian R B F interpolant with e = 21.1 at N = 289 (left) and at N = 1089 Halton
points (right).
24 Meshfree Approximation Methods with M A T L A B
I f we l o o k a t t h e e n t r i e s i n T a b l e 2 . 1 t h e n w e see t h a t c o n t r a r y t o w h a t w e
a n n o u n c e d e a r l i e r t h e r e s u l t s b a s e d o n t h e d i s t a n c e m a t r i x f i t are m o r e a c c u r a t e
t h a n those o b t a i n e d w i t h Gaussians. T h i s w i l l change, however, i f we t r y t o o p t i m i z e
our choice o f t h e shape p a r a m e t e r e (see t h e r e s u l t s o f t h e n e x t e x p e r i m e n t i n
Table 2.2).
I n a second e x p e r i m e n t w e c o n s i d e r t h e same t e s t f u n c t i o n a n d G a u s s i a n basis
functions. N o w , h o w e v e r , w e w a n t t o s t u d y t h e effects o f t h e s h a p e parameter.
T h e r e f o r e , i n F i g u r e 2.5 w e d i s p l a y b o t h t h e m a x i m u m a n d R M S e r r o r s as a f u n c t i o n
o f t h e s h a p e p a r a m e t e r e f o r four f i x e d d a t a sets ( 8 1 , 289, 1089 a n d 4225 H a l t o n
points). T h e s e curves r e v e a l s o m e o f t h e p r o b l e m s a s s o c i a t e d w i t h r a d i a l basis
function i n t e r p o l a t i o n especially w h e n w o r k i n g w i t h globally s u p p o r t e d basis
f u n c t i o n s , i.e., dense m a t r i c e s . W e see t h a t t h e e r r o r s decrease w i t h d e c r e a s i n g e
(of course, t h e y also decrease w i t h d e c r e a s i n g f i l l d i s t a n c e b u t t h a t is n o t w h a t
we are c o n c e r n e d w i t h n o w ) . H o w e v e r , t h e e r r o r c u r v e s are n o t m o n o t o n i c . W e c a n
i d e n t i f y a n o p t i m a l v a l u e o f e for w h i c h b o t h e r r o r s are m i n i m a l ( t h e m i n i m a o f
t h e t w o e r r o r c u r v e s o c c u r a t a l m o s t t h e s a m e p l a c e ) . M o r e o v e r , t h e r e is a v a l u e o f
e at w h i c h t h e c o m p u t a t i o n a l results become u n p r e d i c t a b l e , a n d t h e error curves
b e c o m e e r r a t i c . T h i s p o i n t is a s s o c i a t e d w i t h severe i l l - c o n d i t i o n i n g o f t h e s y s t e m
m a t r i x . Since M A T L A B issues a w a r n i n g w h e n a t t e m p t i n g t o solve a n i l l - c o n d i t i o n e d
l i n e a r s y s t e m , we refer t o t h e s m a l l e s t v a l u e o f e for w h i c h w e d o n o t see a M A T L A B
w a r n i n g as t h e "safe" v a l u e o f e ( f o r a g i v e n set X o f d a t a sites a n d basic f u n c t i o n
<p>). T h e i n t e r e s t i n g fact a b o u t t h e f o u r p l o t s d i s p l a y e d i n F i g u r e 2.5 is t h a t for t h e
s m a l l e r d a t a sets (N = 81 and N = 289) t h e m i n i m u m e r r o r s are o b t a i n e d f o r a
"safe" e, w h i l e for t h e l a r g e r sets ( i V = 1089 a n d N = 4225) the m i n i m u m errors
are o b t a i n e d i n t h e "unsafe" r a n g e . T h e r e f o r e , w e are c o m p u t i n g i n a c e r t a i n " g r a y
zone". W e are o b t a i n i n g h i g h l y a c c u r a t e s o l u t i o n s f r o m severely i l l - c o n d i t i o n e d
l i n e a r systems. W e w i l l c o m e b a c k l a t e r t o t h i s i n t e r e s t i n g f e a t u r e o f r a d i a l basis
f u n c t i o n i n t e r p o l a t i o n ( c a l l e d uncertainty o r trade-off principle). I t is c o n c e i v a b l e
(and i n fact possible [ F o r n b e r g a n d W r i g h t ( 2 0 0 4 ) ] ) t o o b t a i n even m o r e a c c u r a t e
r e s u l t s b y u s i n g a m o r e s t a b l e w a y t o e v a l u a t e t h e r a d i a l basis f u n c t i o n i n t e r p o l a n t
(see t h e d i s c u s s i o n i n C h a p t e r s 16 a n d 1 7 ) .
10
10 jl
10
10"
10"
10
10 15 20
Fig. 2.5 Maximum (dashed/top curve) and R M S (solid/bottom curve) errors vs. e for 81 (top
left), 289 (top right), 1089 (bottom left), and 4225 Halton points (bottom right).
Ac = y,
N N
Y^CjC A >Q
k jk (3.1)
3= 1 k=l
T N
for c = [ c i , . . .,c ] N G R .
I f t h e q u a d r a t i c f o r m ( 3 . 1 ) is zero o n l y for c = 0, t h e n A is c a l l e d positive
definite.
A n i m p o r t a n t p r o p e r t y o f p o s i t i v e d e f i n i t e m a t r i c e s is t h a t a l l t h e i r eigenvalues
are p o s i t i v e , a n d t h e r e f o r e a p o s i t i v e d e f i n i t e m a t r i x is n o n - s i n g u l a r ( b u t c e r t a i n l y
n o t vice v e r s a ) .
27
28 Meshfree Approximation Methods with MATLAB
D e f i n i t i o n 3.2. A c o m p l e x - v a l u e d c o n t i n u o u s f u n c t i o n <fr : R s
- C is c a l l e d positive
definite on K. i fs
N N
(3.2)
3= 1 k=l
XN &
s N
for a n y N p a i r w i s e different p o i n t s x i , . . . , R , and c = [ c i , . . . , CJV] T
<C .
s
T h e f u n c t i o n <E> is c a l l e d strictly positive definite on M. i f t h e q u a d r a t i c form
(3.2) is zero o n l y for c = 0 .
E x a m p l e 3.1. H e r e , a n d t h r o u g h o u t t h i s b o o k , we w i l l d e n o t e t h e s t a n d a r d i n n e r
t x y
product of x and y i n R s
b y x y. W i t h t h i s n o t a t i o n t h e f u n c t i o n <fr(cc) = e ,
3. Positive Definite Functions 29
for y G R s
fixed, is p o s i t i v e d e f i n i t e o n R s
since t h e q u a d r a t i c f o r m i n D e f i n i t i o n 3.2
becomes
N N N N
3=1 k=l
2
N
> 0.
3= 1
3= 1
is also positive definite. Moreover, if at least one of the <&j is strictly positive
definite and the corresponding c - > 0, then $ is strictly
3 positive definite.
(2) $ ( 0 ) > 0.
(3) $ ( - J T ) - $(aj).
(4) Any positive definite function is bounded. In fact,
for e v e r y c G C. T a k i n g c = 1 a n d c = z ( w h e r e i = y/ 1 ) , r e s p e c t i v e l y , w e c a n see
t h a t b o t h <3?(a;) + <&(ic) a n d z ( $ ( x ) <&(cc)) m u s t be r e a l . T h i s , h o w e v e r , is o n l y
possible i f <E>(x) = <fr(x).
F o r t h e p r o o f o f (4) w e l e t N = 2, X\ 0, x 2 = x, a n d choose c\ = | $ ( a ; ) | a n d
c 2 = Q(x). T h e n t h e q u a d r a t i c f o r m i n D e f i n i t i o n 3.2 is
2 2
^ r ^ T c j - c ^ ^ - x) k = 2<Z>(0)\$(x)\ 2 2
- $ ( - a j ) $ ( x ) | $ ( x ) | - $ (x)\(x)\ > 0.
j=ik=i
2 3
2$(0)\$(x)\ - 2\(x)\ > 0.
2
I f | $ ( x ) | > 0, w e d i v i d e b y |<3?(;E)| a n d t h e s t a t e m e n t f o l l o w s i m m e d i a t e l y . I n case
|<E>(a?)| = 0 t h e s t a t e m e n t h o l d s t r i v i a l l y .
P r o p e r t y (5) f o l l o w s i m m e d i a t e l y f r o m ( 4 ) , a n d P r o p e r t y (6) is a consequence o f a
t h e o r e m b y Schur i n t h e field o f l i n e a r a l g e b r a w h i c h s t a t e s t h a t t h e e l e m e n t w i s e (or
H a d a m a r d ) p r o d u c t o f p o s i t i v e ( s e m i - ) d e f i n i t e m a t r i c e s is p o s i t i v e ( s e m i - ) d e f i n i t e .
For m o r e d e t a i l s w e refer t h e r e a d e r t o [ C h e n e y a n d L i g h t (1999)] o r [ W e n d l a n d
(2005a)].
N N
c c
^2 j k<&(xj - x k ) > 0 (3.4)
3 = 1 fc=l
s T N
for any N pairwise different points Xi,..., x^ G M. , and c = [ c i , . . . , C j v ] G M. .
The function $ is strictly positive definite on R s
if the quadratic form (3.4) is
zero only for c = 0.
3. Positive Definite Functions 31
O n e o f t h e m o s t c e l e b r a t e d r e s u l t s o n p o s i t i v e d e f i n i t e f u n c t i o n s is t h e i r c h a r a c t e r -
i z a t i o n i n t e r m s o f F o u r i e r t r a n s f o r m s e s t a b l i s h e d b y B o c h n e r i n 1932 (for s = 1)
a n d 1933 (for g e n e r a l s).
S
T h e o r e m 3 . 3 ( B o c h n e r ) . A (complex-valued) function <E> e C ( I R ) is positive def-
s
inite on R if and only if it is the Fourier transform of a finite non-negative Borel
s
measure fi on R , i.e.
<2>(cc) = jl(x)
1
J ,ixy x e
N N
1
J E Cj e -zxj-y
E
fc=i
Cite
ix y
k
dLi(y)
Jut*
J=l
J E N
Cje
-txj-y dfj,(y) > 0.
32 Meshfree Approximation Methods with MATLAB
t x y
R e m a r k 3 . 1 . W e c a n see f r o m T h e o r e m 3.3 t h a t t h e f u n c t i o n $(x) = e of
E x a m p l e 3.1 c a n be c o n s i d e r e d as t h e fundamental positive definite function since
a l l o t h e r p o s i t i v e d e f i n i t e f u n c t i o n s are o b t a i n e d as ( i n f i n i t e ) l i n e a r c o m b i n a t i o n s o f
t h i s f u n c t i o n . W h i l e P r o p e r t y (1) o f T h e o r e m 3.1 i m p l i e s t h a t l i n e a r c o m b i n a t i o n s
o f <& w i l l a g a i n b e p o s i t i v e d e f i n i t e , t h e r e m a r k a b l e c o n t e n t o f B o c h n e r ' s Theorem
is t h e fact t h a t i n d e e d all p o s i t i v e d e f i n i t e f u n c t i o n s are g e n e r a t e d b y
X \ \ J { 0 : O is o p e n a n d fi(0) = 0}.
s
T h e o r e m 3 . 4 . Let /J, be a non-negative finite Borel measure on R whose carrier
is a set of nonzero Lebesgue measure. Then the Fourier transform of /j, is strictly
s
positive definite on R .
N o w let v/(2ir) h
AT
c e i X i v
g{y) = Y. i ~ ' '
3= 1
r e m a i n i n g w a y t o m a k e t h e a b o v e i n e q u a l i t y a n e q u a l i t y is i f t h e c a r r i e r o f fi is
c o n t a i n e d i n t h e zero set o f g, i.e., has Lebesgue m e a s u r e zero. T h i s , however, is
ruled out i n the hypothesis o f the theorem.
s
C o r o l l a r y 3 . 1 . Let f be a continuous non-negative function in Li(R ) which is not
s
identically zero. Then the Fourier transform of f is strictly positive definite on M. .
//() = ( f(x)dx.
JB
T h e n t h e c a r r i e r o f fi is e q u a l t o t h e (closed) s u p p o r t o f / . H o w e v e r , since / is
n o n - n e g a t i v e a n d n o t i d e n t i c a l l y e q u a l t o zero, i t s s u p p o r t has p o s i t i v e L e b e s g u e
measure, a n d hence t h e F o u r i e r t r a n s f o r m o f / is s t r i c t l y p o s i t i v e d e f i n i t e b y t h e
preceding theorem.
F i n a l l y , a c r i t e r i o n t o check w h e t h e r a g i v e n f u n c t i o n is s t r i c t l y p o s i t i v e d e f i n i t e
is g i v e n i n [ W e n d l a n d ( 2 0 0 5 a ) ] .
s
T h e o r e m 3.5. Let <& be a continuous function in L i ( R ) . <& is strictly positive
definite if and only if $ is bounded and its Fourier transform is non-negative and
not identically equal to zero.
T h e o r e m 3.5 is o f f u n d a m e n t a l i m p o r t a n c e a n d we w i l l c o m e b a c k t o t h i s t h e o r e m
several t i m e s l a t e r o n . I n f a c t , t h e p r o o f o f T h e o r e m 3.5 i n [ W e n d l a n d (2005a)] shows
t h a t i f <& ^ 0 ( w h i c h i m p l i e s t h a t t h e n also <E ^ 0) w e need t o ensure o n l y
t h a t <l> be n o n - n e g a t i v e i n o r d e r for <3> t o be s t r i c t l y p o s i t i v e d e f i n i t e .
Here
for s = 1,
for s>2,
5
and J ( _ 2 ) / 2 ^ the classical
s Bessel function of the first kind of order (s 2 ) / 2 .
T h e o r e m 3.7. A continuous function <p : [0, oo) * R such that r t> r ~ cp(r) s 1
e
s
L\[0, oo) is strictly positive definite and radial on R if and only if the s-dimensional
Fourier transform
Jo
where fi is a finite non-negative Borel measure on [0, oo) not concentrated at the
origin.
Since t h e e x p o n e n t i a l f u n c t i o n is p o s i t i v e a n d t h e m e a s u r e is n o n - n e g a t i v e , i t f o l l o w s
t h a t /J, m u s t be t h e zero m e a s u r e . H o w e v e r , t h e n cp is i d e n t i c a l l y e q u a l t o zero.
S
T h e r e f o r e , a n o n - t r i v i a l f u n c t i o n cp t h a t is p o s i t i v e d e f i n i t e a n d r a d i a l o n R for a l l
s c a n have n o zeros. T h i s i m p l i e s i n p a r t i c u l a r t h a t
4.1 E x a m p l e 1: Gaussians
W e c a n n o w show t h a t t h e Gaussian
2 x | 1 2
&(x) = e- H , > 0 , (4.1)
is s t r i c t l y p o s i t i v e d e f i n i t e ( a n d r a d i a l ) o n K s
for a n y s. T h i s is d u e t o t h e fact t h a t
t h e F o u r i e r t r a n s f o r m o f a G a u s s i a n is essentially a G a u s s i a n . I n fact,
, 1
a n d t h i s is p o s i t i v e i n d e p e n d e n t o f t h e space d i m e n s i o n s. I n p a r t i c u l a r , for e =
w e have <E> = P l o t s o f G a u s s i a n R B F s were p r e s e n t e d i n F i g . 2 . 1 . C l e a r l y , t h e
Gaussians are i n f i n i t e l y d i f f e r e n t i a b l e . S o m e o f i t s d e r i v a t i v e s (as w e l l as t h o s e o f
m a n y o t h e r R B F s ) are c o l l e c t e d i n A p p e n d i x D .
37
38 Meshfree Approximation Methods with MATLAB
A n o t h e r a r g u m e n t t o s h o w t h a t G a u s s i a n s are s t r i c t l y p o s i t i v e d e f i n i t e a n d r a d i a l
S
on R for a n y s t h a t avoids d e a l i n g w i t h F o u r i e r t r a n s f o r m s w i l l b e c o m e a v a i l a b l e
l a t e r . I t w i l l m a k e use o f c o m p l e t e l y m o n o t o n e f u n c t i o n s .
R e c a l l t h a t P r o p e r t y ( 1 ) o f T h e o r e m 3.1 shows t h a t a n y finite n o n - n e g a t i v e l i n -
ear c o m b i n a t i o n o f ( s t r i c t l y ) p o s i t i v e d e f i n i t e f u n c t i o n s is a g a i n ( s t r i c t l y ) p o s i t i v e
d e f i n i t e . M o r e o v e r , w e j u s t saw t h a t G a u s s i a n s are s t r i c t l y p o s i t i v e d e f i n i t e a n d r a -
S
d i a l o n a l l R . N o w , t h e S c h o e n b e r g c h a r a c t e r i z a t i o n o f f u n c t i o n s t h a t are ( s t r i c t l y )
S
p o s i t i v e d e f i n i t e a n d r a d i a l o n a n y R , T h e o r e m 3.8, s t a t e s t h a t all s u c h f u n c t i o n s
are g i v e n as i n f i n i t e l i n e a r c o m b i n a t i o n s o f G a u s s i a n s . T h e r e f o r e , t h e G a u s s i a n s c a n
be v i e w e d as t h e f u n d a m e n t a l m e m b e r o f t h e f a m i l y o f f u n c t i o n s t h a t are s t r i c t l y
S
positive definite and radial on R for a l l s.
Since Gaussians p l a y a c e n t r a l r o l e i n s t a t i s t i c s t h i s is a g o o d p l a c e t o m e n t i o n
t h a t p o s i t i v e d e f i n i t e f u n c t i o n s are u p t o a n o r m a l i z a t i o n $ ( 0 ) = 1 i d e n t i c a l
w i t h characteristic functions of d i s t r i b u t i o n functions i n statistics.
4.2 E x a m p l e 2: L a g u e r r e - G a u s s i a n s
I n o r d e r t o o b t a i n a g e n e r a l i z a t i o n o f G a u s s i a n s w e s t a r t w i t h t h e generalized La-
2
guerre polynomials L % / o f degree n a n d o r d e r s/2 d e f i n e d b y t h e i r R o d r i g u e s f o r m u l a
(see, e.g., f o r m u l a (6.2.1) i n [ A n d r e w s et al. ( 1 9 9 9 ) ] )
tj. s/2 JTl . .
L / 2 ( i ) = i e < n + / 2
V ^ ( " ' ' ) ' "
A n e x p l i c i t f o r m u l a for t h e g e n e r a l i z e d L a g u e r r e p o l y n o m i a l s is
a n d list t h e i r F o u r i e r t r a n s f o r m s as
_ H"ll 2
n n | | 2 l '
N o t e t h a t t h e d e f i n i t i o n o f t h e L a g u e r r e - G a u s s i a n s d e p e n d s o n t h e space d i m e n s i o n
S
s. T h e r e f o r e t h e y are s t r i c t l y p o s i t i v e d e f i n i t e a n d r a d i a l o n R ( a n d b y L e m m a 3.1
C T
also o n R for a n y a < s).
Laguerre-Gaussian f u n c t i o n s for some s p e c i a l choices o f s a n d n are l i s t e d i n
T a b l e 4 . 1 . F i g u r e 4 . 1 shows a L a g u e r r e - G a u s s i a n for s l , n = 2, a n d for s =
2 , n = 2 d i s p l a y e d w i t h a shape p a r a m e t e r e = 3 a n d scaled so t h a t <&(0) = 1.
M o r e o v e r , t h e L a g u e r r e - G a u s s i a n s are i n f i n i t e l y s m o o t h for a l l choices o f n a n d s.
N o t e t h a t t h e L a g u e r r e - G a u s s i a n s ( w h i l e b e i n g s t r i c t l y p o s i t i v e d e f i n i t e func-
t i o n s ) are n o t p o s i t i v e . Since t h e L a g u e r r e - G a u s s i a n s are o s c i l l a t o r y f u n c t i o n s w e
^ 1
4- Examples of Strictly Positive Definite Radial Functions 39
k n o w f r o m T h e o r e m 3.9 t h a t t h e y c a n n o t be s t r i c t l y p o s i t i v e d e f i n i t e a n d r a d i a l o n
s
R for a l l s. W e w i l l e n c o u n t e r these f u n c t i o n s l a t e r i n t h e c o n t e x t o f a p p r o x i m a t e
m o v i n g least squares a p p r o x i m a t i o n (c.f. Chapter 26).
s n = 1 n = 2
Fig. 4.1 Laguerre-Gaussians with s = l,n = 2 (left) and s = 2, n = 2 (right) centered at the
origin.
4.3 E x a m p l e 3: P o i s s o n R a d i a l Functions
A n o t h e r class o f o s c i l l a t o r y f u n c t i o n s t h a t are s t r i c t l y p o s i t i v e d e f i n i t e a n d r a d i a l
s CT
on R (and all R for a < s) w e r e r e c e n t l y s t u d i e d b y F o r n b e r g a n d c o - w o r k e r s (see
[ F o r n b e r g et al. (2004)] a n d also [ F l y e r ( 2 0 0 6 ) ] ) . T h e s e f u n c t i o n s are o f t h e f o r m
*(*> = ^ P . (*<)
^1. t
40 Meshfree Approximation Methods with MATLAB
T h e o r e m 3.6 a n d t h e r e f o r e c a n be v i e w e d as t h e f u n d a m e n t a l m e m b e r o f t h e f a m i l y
S
o f f u n c t i o n s t h a t are s t r i c t l y p o s i t i v e d e f i n i t e a n d r a d i a l o n R for f i x e d s.
S c h o e n b e r g s h o w e d t h a t t h e f u n c t i o n s f 2 are g i v e n b y
s
1 S
w h e r e co -i
s denotes t h e area o f t h e u n i t s p h e r e S**" i n R , a n d da d e n o t e s t h e
s - 1
usual measure o n S ' .
T h e Poisson f u n c t i o n s are a n o t h e r g e n e r a l i z a t i o n o f G a u s s i a n s ( t h e f u n d a m e n t a l
S
s t r i c t l y positive definite radial function o n R for a l l s) since t h e f o l l o w i n g l i m i t -
i n g r e l a t i o n due t o J o h n v o n N e u m a n n h o l d s (see t h e d i s c u s s i o n i n [Schoenberg
(1938a)]):
- 7 - 2
l i m n (rV2s)
s = e .
s>oo
Since t h e P o i s s o n r a d i a l f u n c t i o n s are d e f i n e d i n t e r m s o f Bessel f u n c t i o n s t h e y are
also band-limited, i.e., t h e i r F o u r i e r t r a n s f o r m has c o m p a c t s u p p o r t . I n fact, t h e
F F
F o u r i e r t r a n s f o r m o f <E> i n R , a < s, is g i v e n b y (see [ F l y e r ( 2 0 0 6 ) ] )
s = 2 s = 3 s = 4 s = 5
2
Fig. 4.2 Poisson functions with s = 2 (left) and s 3 (right) centered at the origin in R .
4- Examples of Strictly Positive Definite Radial Functions 41
4.4 E x a m p l e 4: M a t e r n Functions
A f o u r t h e x a m p l e o f s t r i c t l y p o s i t i v e d e f i n i t e f u n c t i o n s is g i v e n b y t h e class o f
Matern functions w h i c h are q u i t e c o m m o n i n t h e s t a t i s t i c s l i t e r a t u r e (see, e.g.,
[ M a t e r n (1986)] or [Stein ( 1 9 9 9 ) ] )
2
(1+11*11)6-11*11 (3 +311*11+ ||a || )e3
4.5 E x a m p l e 5: G e n e r a l i z e d I n v e r s e Multiquadrics
2
Fig. 4.3 Matern functions with B = i (left) and B = ^5 (right) centered at the origin in R .
f u n c t i o n s w e w i l l be able t o s h o w t h a t i n f a c t w e n e e d t o r e q u i r e o n l y 8 > 0, a n d
S
t h e r e f o r e t h e g e n e r a l i z e d inverse m u l t i q u a d r i c s are s t r i c t l y p o s i t i v e d e f i n i t e o n ]R
for a n y s.
T h e " o r i g i n a l " inverse m u l t i q u a d r i c was i n t r o d u c e d b y H a r d y i n t h e e a r l y 1970s
a n d c o r r e s p o n d s t o t h e value 8 1/2. T h e s p e c i a l choice 8=1 was r e f e r r e d t o as
inverse quadratic i n v a r i o u s p a p e r s o f F o r n b e r g a n d c o - w o r k e r s (see, e.g., [Fornberg
a n d W r i g h t ( 2 0 0 4 ) ] ) . T h e s e t w o f u n c t i o n s are d i s p l a y e d i n F i g u r e 4.4 u s i n g a s h a p e
p a r a m e t e r e = 5.
1-. r , 1
Fig. 4.4 Inverse multiquadric (/3 = ^, left) and inverse quadratic (8=1, right) centered at the
2
origin in R .
4.6 E x a m p l e 6: T r u n c a t e d P o w e r Functions
W e n o w present a n e x a m p l e o f a f a m i l y o f s t r i c t l y p o s i t i v e d e f i n i t e f u n c t i o n s w i t h
compact support. N o t e t h a t d u e t o t h e o b s e r v a t i o n m a d e i n T h e o r e m 3.9 a t t h e e n d
S
of t h e p r e v i o u s c h a p t e r , t h e y c a n n o t be s t r i c t l y p o s i t i v e d e f i n i t e o n R for a l l s.
4- Examples of Strictly Positive Definite Radial Functions 43
ipi(r) = ( 1 - r ) i+ (4.7)
s
give rise t o s t r i c t l y p o s i t i v e d e f i n i t e a n d r a d i a l f u n c t i o n s o n M p r o v i d e d I satisfies
t > L f J + 1- F i n d i n g t h e F o u r i e r t r a n s f o r m o f t h e t r u n c a t e d p o w e r f u n c t i o n is
r a t h e r i n v o l v e d . F o r d e t a i l s we refer t o [ W e n d l a n d ( 2 0 0 5 a ) ] . W e w i l l l a t e r use a
s i m p l e r t e s t based o n m u l t i p l y m o n o t o n e f u n c t i o n s t o e s t a b l i s h t h e s t r i c t p o s i t i v e
definiteness o f t h e t r u n c a t e d p o w e r f u n c t i o n s . I n (4.7) we used t h e cutoff function
() + w h i c h is defined b y
x, for x > 0,
0, for x < 0.
T h e c u t o f f f u n c t i o n c a n be i m p l e m e n t e d c o n v e n i e n t l y i n M A T L A B u s i n g t h e max
f u n c t i o n , i.e., i f f x is a v e c t o r o f f u n c t i o n values o f / for d i f f e r e n t choices o f x, t h e n
m a x ( f x , 0 ) c o m p u t e s (f(x)) . + W e also p o i n t o u t t h a t t h e expressions o f t h e f o r m
( 1 r) e
+ are t o be i n t e r p r e t e d as ( ( 1 r)+)*, i.e., we first a p p l y t h e c u t o f f f u n c t i o n ,
and then the power.
T w o different t r u n c a t e d p o w e r f u n c t i o n s ( w i t h i = 2 , 4 ) are d i s p l a y e d i n F i g -
u r e 4.5. W h i l e n o n e o f t h e t r u n c a t e d p o w e r f u n c t i o n s are d i f f e r e n t i a b l e a t t h e o r i g i n ,
t h e s m o o t h n e s s a t t h e b o u n d a r y o f t h e s u p p o r t increases w i t h I .
Fig. 4.5 Truncated power function with i = 2 (left) and = 4 (right) centered at the origin in
4.7 E x a m p l e 7: P o t e n t i a l s a n d W h i t t a k e r R a d i a l F u n c t i o n s
(4.8)
^1 it.
44 Meshfree Approximation Methods with MATLAB
s
T h e n <& = ip(\\ ||) is s t r i c t l y p o s i t i v e d e f i n i t e a n d r a d i a l o n E p r o v i d e d k > |_|J + 2
(see also T h e o r e m 5.5 b e l o w ) . T h i s c a n be v e r i f i e d b y c o n s i d e r i n g t h e quadratic
form
N N oo N N
- 1
EEw^ *!^ /
3= 1 k=l J
3= 1
c
^2^2 jCkVk-i(t\\xj
k=l
-x \\)f(t)dt
k
w h i c h is n o n - n e g a t i v e since t h e t r u n c a t e d p o w e r f u n c t i o n ^ _ i ( | | | | ) is s t r i c t l y
p o s i t i v e d e f i n i t e b y E x a m p l e 6, a n d / is n o n - n e g a t i v e . Since / is also a s s u m e d t o
be n o t i d e n t i c a l l y e q u a l t o zero, t h e o n l y w a y for t h e q u a d r a t i c f o r m t o e q u a l zero
is i f c = 0 , a n d therefore ip is s t r i c t l y p o s i t i v e d e f i n i t e .
For e x a m p l e , i f we t a k e f(t) = t@, 8 > 0, t h e n we get
- r(fc)r(/3 + 1 )
{ X } ( 4 9 )
~ T{k + 8 + l)\\x\\^- -
W h i l e these f u n c t i o n s are s t r i c t l y p o s i t i v e d e f i n i t e a n d r a d i a l t h e y are also s i n g u l a r
a t t h e o r i g i n a n d therefore n o t useful for o u r p u r p o s e s . H o w e v e r , these f u n c t i o n s
are u p t o s c a l i n g g e n e r a l i z a t i o n s o f t h e Coulomb potential (for 3 = 0 ) , a n d
c a n therefore be g i v e n a p h y s i c a l i n t e r p r e t a t i o n .
a 1
A n o t h e r p o s s i b i l i t y is t o t a k e f(t) = t e~^ , a > 0,/5 > 0. T h e n we get
$ f x = |Ml t f c Q )
- ^r(i+a)r fc) -^ ( e x u 1 Q )
* W /3i+(fc+)/2r(fc+<*+2) e x
i *-- ; L U
(kM ^ ,(k+a+l)/2
{a k)/2 (|Tij) + ( l + o ) M i _ ( f c _ )/
a 2 i ( f c + Q + i)/2 (|NT)) "
Here M M ) I / is t h e W h i t t a k e r - M f u n c t i o n , a c o n f l u e n t h y p e r g e o m e t r i c f u n c t i o n (see,
e.g., Chapter 13 o f [ A b r a m o w i t z a n d S t e g u n (1972)]). W h e n v is a h a l f - i n t e g e r
( w h i c h is, e.g., t h e case for integer k a n d a) f o r m u l a (4.10) s i m p l i f i e s s i g n i f i c a n t l y .
E x a m p l e s for v a r i o u s integer values o f k a n d a are l i s t e d i n T a b l e 4.4. N o t e t h a t
these f u n c t i o n s are n o t defined a t t h e o r i g i n . H o w e v e r , t h e y c a n be m a d e ( o n l y )
c o n t i n u o u s a t t h e o r i g i n . P l o t s o f t w o o f these f u n c t i o n s are p r o v i d e d i n F i g u r e 4.6.
N o t e t h a t o n l y t h e f u n c t i o n s for k > 3 are g u a r a n t e e d t o be s t r i c t l y p o s i t i v e d e f i n i t e
3
and r a d i a l on IR .
Table 4.4 Whittaker radial functions <I> for various choices of k and a.
a fc = 2 fc = 3
2 2
P- 2||a;| + (/? + 2 | | * | | ) e ~ T W 2
P - 4/3j *|| + 6 | x\\
2
- (2/3\\x\\ + 6 | | * | | ) e " T i r 2
1
3 i
P P
Fig. 4.6 Whittaker radial functions for a = 0 and 3 = 1 with k = 2 (left) and k = 3 (right)
2
centered at the origin in R .
4.8 E x a m p l e 8: I n t e g r a t i o n A g a i n s t S t r i c t l y P o s i t i v e
Definite K e r n e l s
<p{r)= / K(t,r)f(t)dt
Jo
S
gives rise t o $ = tp(\\ ||) b e i n g s t r i c t l y p o s i t i v e definite o n R . B y c h o o s i n g / and
K a p p r o p r i a t e l y we can o b t a i n b o t h g l o b a l l y s u p p o r t e d a n d c o m p a c t l y supported
functions.
For example, t h e m u l t i p l y m o n o t o n e f u n c t i o n s i n W i l l i a m s o n ' s c h a r a c t e r i z a t i o n
1
T h e o r e m 5.4 are covered b y t h i s general t h e o r e m b y t a k i n g K(t, r) = (l rt)^T and
/ a n a r b i t r a r y p o s i t i v e f u n c t i o n i n L \ so t h a t dfi(t) = f(t)dt. Also, functions t h a t
S
are s t r i c t l y p o s i t i v e definite a n d r a d i a l o n R for a l l s (or e q u i v a l e n t l y c o m p l e t e l y
rt
m o n o t o n e f u n c t i o n s ) are covered b y c h o o s i n g K(t,r) = e~ .
4.9 Summary
T o s u m m a r i z e t h e t h e o r y s u r v e y e d t h u s far w e c a n say t h a t a n y m u l t i v a r i a t e ( r a d i a l )
f u n c t i o n $ whose F o u r i e r t r a n s f o r m is n o n - n e g a t i v e can be used t o g e n e r a t e a basis
for t h e s c a t t e r e d d a t a i n t e r p o l a t i o n p r o b l e m b y s h i f t i n g i t t o t h e d a t a sites. The
f u n c t i o n $ can be p o s i t i v e , o s c i l l a t o r y , o r have c o m p a c t s u p p o r t . H o w e v e r , i f $ has
S
any zeros t h e n i t c a n n o t be s t r i c t l y p o s i t i v e d e f i n i t e o n R for a l l choices of s.
Chapter 5
W e b e g i n w i t h t h e f o r m e r case. T o t h i s e n d we n o w i n t r o d u c e a class o f f u n c t i o n s
t h a t is v e r y closely r e l a t e d t o p o s i t i v e d e f i n i t e r a d i a l f u n c t i o n s a n d leads t o a s i m p l e
c h a r a c t e r i z a t i o n o f such f u n c t i o n s .
_ e r
E x a m p l e 5 . 2 . T h e f u n c t i o n </?(r) = e , e > 0, is c o m p l e t e l y m o n o t o n e o n [ 0 , o o )
since
W r
(-l)V ( r ) = e e- > 0, t = 0,1,2,....
Some p r o p e r t i e s o f c o m p l e t e l y m o n o t o n e f u n c t i o n s t h a t c a n be f o u n d i n [Cheney
a n d L i g h t (1999); Feller ( 1 9 6 6 ) ; W i d d e r (1941)] are:
(1) A non-negative f i n i t e l i n e a r c o m b i n a t i o n o f c o m p l e t e l y m o n o t o n e f u n c t i o n s is
completely monotone.
(2) T h e p r o d u c t o f t w o c o m p l e t e l y m o n o t o n e f u n c t i o n s is c o m p l e t e l y m o n o t o n e .
47
48 Meshfree Approximation Methods with M A T L A B
T h e o r e m 5.2. A function <p is completely monotone on [0, oo) if and only if <> =
2 s
<p(\\ \\ ) is positive definite and radial on R for all s.
2
w i t h a finite n o n - n e g a t i v e B o r e l m e a s u r e / i . T h e r e f o r e , <&(x) = <^(||cc|j ) has the
representation
/OO
Jo
s
T o see t h a t t h i s f u n c t i o n is p o s i t i v e d e f i n i t e o n a n y R w e consider t h e q u a d r a t i c
form
N N r oo N N
j=l k=l J
j=l k=l
_ e r
E x a m p l e 5.4. Since we showed above t h a t t h e f u n c t i o n s <p(r) = e , s > 0, a n d
p(r) = 1 / ( 1 + r)@, 8 > 0, are c o m p l e t e l y m o n o t o n e o n [0, o o ) , a n d since t h e y are
2
also n o t c o n s t a n t we k n o w f r o m T h e o r e m 5.3 t h a t t h e Gaussians <&(cc) = <>(||cc|| ) =
2 2
e - ^ I M I ^ e > 0, a n d inverse m u l t i q u a d r i c s $(cc) = ^ ( | | a ; | | ) = 1 / ( 1 + | | c c | | ) ^ , 8 > 0,
s
are s t r i c t l y p o s i t i v e d e f i n i t e a n d r a d i a l o n R. for a l l s. N o t o n l y is t h e test for
c o m p l e t e m o n o t o n i c i t y a s i m p l e r one t h a n c a l c u l a t i o n o f t h e F o u r i e r t r a n s f o r m s ,
b u t we also are able t o v e r i f y s t r i c t p o s i t i v e definiteness o f t h e inverse m u l t i q u a d r i c s
w i t h o u t a n y dependence o f s o n 8.
r r 2 r i r
Since c o n v e x i t y o f tp m e a n s t h a t < ^ ( i + ) < ^ ( ) + ^ ( 2 ) ; o r s i p l y p"(r)
m > 0 if
ip" exists, a m u l t i p l y m o n o t o n e f u n c t i o n is i n essence j u s t a c o m p l e t e l y m o n o t o n e
f u n c t i o n whose m o n o t o n i c i t y is " t r u n c a t e d " .
E x a m p l e 5 . 5 . T h e t r u n c a t e d p o w e r f u n c t i o n (c.f. (4.7))
e
<Pi(r) = (1 - r) +
W e saw i n S e c t i o n 4.6 t h a t t h e t r u n c a t e d p o w e r f u n c t i o n s l e a d t o r a d i a l f u n c t i o n s
s
t h a t are s t r i c t l y p o s i t i v e d e f i n i t e o n R p r o v i d e d > [s/2\ + 1.
E x a m p l e 5 . 6 . I f w e define t h e i n t e g r a l o p e r a t o r I b y
/OO
1
p(r)= / (l-r*)*- ^*), (5.2)
( a n d is e q u i v a l e n t p r o v i d e d we e x t e n d W i l l i a m s o n ' s w o r k t o i n c l u d e c o n t i n u i t y a t
the origin).
W e c a n see f r o m Sections 4.6 a n d 4.7 t h a t m u l t i p l y m o n o t o n e f u n c t i o n s give rise
to positive definite r a d i a l functions. S u c h a c o n n e c t i o n was first n o t e d i n [ A s k e y
(1973)] ( a n d i n t h e o n e - d i m e n s i o n a l case b y P o l y a ) u s i n g t h e t r u n c a t e d p o w e r func-
t i o n s o f S e c t i o n 4.6.
I n t h e R B F l i t e r a t u r e t h e f o l l o w i n g t h e o r e m was s t a t e d i n [ M i c c h e l l i ( 1 9 8 6 ) ] ,
a n d t h e n refined i n [ B u h m a n n (1993a)]:
2 e
F^(r) = - 7 = f / (1 - t ) t^J_ (rt)dt
+ 1/2
vi Jo 1
= \/f J
1
(1 - t )2 e
cos(rt)dt
52 Meshfree Approximation Methods with M A T L A B
As we m e n t i o n e d i n t h e i n t r o d u c t i o n i t is n o t a n easy m a t t e r t o use p o l y n o m i a l s
t o p e r f o r m m u l t i v a r i a t e s c a t t e r e d d a t a i n t e r p o l a t i o n . O n l y i f t h e d a t a sites are i n
c e r t a i n special l o c a t i o n s c a n we g u a r a n t e e well-posedness o f m u l t i v a r i a t e p o l y n o m i a l
i n t e r p o l a t i o n . W e n o w address t h i s p r o b l e m .
S
D e f i n i t i o n 6 . 1 . W e c a l l a set o f p o i n t s X = {x\,... ,x^} C R m-unisolvent if
t h e o n l y p o l y n o m i a l o f t o t a l degree at m o s t m i n t e r p o l a t i n g zero d a t a o n X is t h e
zero p o l y n o m i a l .
T h i s d e f i n i t i o n guarantees a u n i q u e s o l u t i o n for i n t e r p o l a t i o n t o g i v e n d a t a at a
subset o f c a r d i n a l i t y M = ( * ) m s
f t n e
p o i n t s x \ , . . . , XN b y a p o l y n o m i a l o f degree
m. Here M is t h e d i m e n s i o n o f t h e linear space o f p o l y n o m i a l s o f t o t a l degree
less t h a n or equal t o m i n s variables.
S
For p o l y n o m i a l i n t e r p o l a t i o n a t N d i s t i n c t d a t a sites i n R t o be a w e l l - p o s e d
p r o b l e m , t h e p o l y n o m i a l degree needs t o be chosen a c c o r d i n g l y , i.e., we need M =
N, a n d t h e d a t a sites need t o f o r m a n m - u n i s o l v e n t set. T h i s is r a t h e r r e s t r i c t i v e .
2
For example, t h i s i m p l i e s t h a t p o l y n o m i a l i n t e r p o l a t i o n at N = 7 points i n R
can n o t be done i n a u n i q u e w a y since we c o u l d either a t t e m p t t o use b i v a r i a t e
q u a d r a t i c p o l y n o m i a l s (for w h i c h M = 6 ) , or b i v a r i a t e c u b i c p o l y n o m i a l s ( w i t h
M = 10). T h e r e exists n o space o f b i v a r i a t e p o l y n o m i a l s for w h i c h M = 7.
W e w i l l see i n t h e n e x t c h a p t e r t h a t m - u n i s o l v e n t sets p l a y a n i m p o r t a n t role i n
t h e c o n t e x t o f c o n d i t i o n a l l y p o s i t i v e definite f u n c t i o n s . T h e r e , however, even t h o u g h
we w i l l be interested i n i n t e r p o l a t i n g N pieces o f d a t a , t h e p o l y n o m i a l degree w i l l
be s m a l l ( u s u a l l y m = 1, 2, 3 ) , a n d t h e r e s t r i c t i o n s i m p o s e d o n t h e l o c a t i o n s o f t h e
d a t a sites b y t h e u n i s o l v e n c y c o n d i t i o n s w i l l be r a t h e r m i l d .
A sufficient c o n d i t i o n ( t o be f o u n d i n [ C h u i (1988)], C h . 9) o n t h e points
X \ , . . . , xpj t o f o r m a n m - u n i s o l v e n t set i n R 2
is
2
T h e o r e m 6 . 1 . Suppose {Lo,..., L m } is a set o / m + 1 distinct lines in R , and that
U = . . . , UM} is a set of M = ( m + l ) ( m + 2 ) / 2 distinct points such that the
53
54 Meshfree Approximation Methods with M A T L A B
first point lies on L Q , the next two points lie on L \ but not on L Q , and so on, so that
the last m + 1 points lie on L m but not on any of the previous lines LQ, . . . , L m _ i .
Then there exists a unique interpolation polynomial of total degree at most m to
arbitrary data given at the points in U. Furthermore, if the data sites {x,..., x^}
2
contain U as a subset then they form an m-unisolvent set on R .
Proof. W e use i n d u c t i o n o n m . F o r m = 0 t h e r e s u l t is t r i v i a l . T a k e R t o be t h e
m a t r i x a r i s i n g f r o m p o l y n o m i a l i n t e r p o l a t i o n a t t h e p o i n t s i n IA, i.e.,
p(t*i)=0, i = l,...,M,
t h e n p is t h e zero p o l y n o m i a l .
For each i = 1 , . . . , m, let t h e e q u a t i o n o f t h e l i n e Li be g i v e n b y
OHX + fay = 7 i ,
2
w h e r e x = (x, y) R .
Suppose n o w t h a t p i n t e r p o l a t e s zero d a t a at a l l t h e p o i n t s U i as s t a t e d above.
Since p reduces t o a u n i v a r i a t e p o l y n o m i a l o f degree m o n L m w h i c h vanishes a t
m + 1 distinct points on L m , i t follows t h a t p vanishes i d e n t i c a l l y o n L m , a n d so
p(x, y) = (a xm + 3y
m - ~/m)q(x, y),
w h e r e q is a p o l y n o m i a l o f degree m 1. B u t n o w q satisfies t h e h y p o t h e s i s o f t h e
x
t h e o r e m w i t h m r e p l a c e d b y m 1 a n d U r e p l a c e d b y U c o n s i s t i n g o f t h e first ( " ^ )
p o i n t s o f U. B y i n d u c t i o n , therefore q = 0, a n d t h u s p = 0. T h i s establishes t h e
uniqueness o f t h e i n t e r p o l a t i o n p o l y n o m i a l . T h e l a s t s t a t e m e n t o f t h e t h e o r e m is
obvious.
A s i m i l a r t h e o r e m was a l r e a d y p r o v e d i n [ C h u n g a n d Y a o ( 1 9 7 7 ) ] . T h e o r e m 6.1
S
c a n be generalized t o R b y using hyperplanes. T h e p r o o f is c o n s t r u c t e d w i t h t h e
h e l p o f a n a d d i t i o n a l i n d u c t i o n o n s. C h u i also gives a n e x p l i c i t expression for t h e
d e t e r m i n a n t o f t h e m a t r i x associated w i t h ( p o l y n o m i a l ) i n t e r p o l a t i o n a t t h e set o f
p o i n t s U.
R e m a r k 6.1. F o r l a t e r reference w e n o t e t h a t ( m l ) - u n i s o l v e n c y o f t h e p o i n t s
x\,..., XN is e q u i v a l e n t t o t h e fact t h a t t h e m a t r i x P w i t h
has f u l l ( c o l u m n - ) r a n k . F o r N = M t h i s is t h e p o l y n o m i a l i n t e r p o l a t i o n m a t r i x .
2
Example 6.1. A s can easily be v e r i f i e d , t h r e e c o l l i n e a r p o i n t s i n R are n o t 1-
u n i s o l v e n t , since a linear i n t e r p o l a n t , i.e., a p l a n e t h r o u g h t h r e e a r b i t r a r y h e i g h t s
a t these t h r e e c o l l i n e a r p o i n t s is n o t u n i q u e l y d e t e r m i n e d . O n t h e o t h e r h a n d , i f a
2
set o f p o i n t s i n R c o n t a i n s t h r e e n o n - c o l l i n e a r p o i n t s , t h e n i t is 1-unisolvent.
x 4-
6. Scattered Data Interpolation with Polynomial Precision 55
6.2 E x a m p l e : R e p r o d u c t i o n of L i n e a r F u n c t i o n s U s i n g
Gaussian RBFs
I f we do insist o n r e p r o d u c t i o n o f l i n e a r f u n c t i o n s t h e n t h e t o p p a r t o f F i g u r e 6.1
shows a Gaussian R B F i n t e r p o l a n t (e = 6) t o t h e b i v a r i a t e l i n e a r f u n c t i o n f(x,y) =
(x + y)/2 based o n 1089 u n i f o r m l y spaced p o i n t s i n t h e u n i t square a l o n g w i t h
t h e absolute error. C l e a r l y t h e i n t e r p o l a n t is n o t c o m p l e t e l y p l a n a r n o t even t o
m a c h i n e precision.
F o r t u n a t e l y , t h e r e is a s i m p l e r e m e d y for t h i s p r o b l e m . A l l we need t o d o
is a d d t h e p o l y n o m i a l f u n c t i o n s x i> 1, x i x, and x t> y t o the- basis
_ e - X l _ e - X J V
{e H' H ,...,e H' H } we have t h u s far been u s i n g t o o b t a i n o u r i n t e r -
polant. However, n o w we have N + 3 u n k n o w n s , n a m e l y t h e coefficients c,
k
k = 1 , . . . , N + 3, i n t h e e x p a n s i o n
N
e 2 x 2 2
V (x)
f = ^ 2 c e - k ^ - ^ + c N + 1 + c N + 2 x + c N + 3 y, x = {x,y) e M ,
fc=i
a n d we have o n l y N c o n d i t i o n s t o d e t e r m i n e t h e m , n a m e l y t h e i n t e r p o l a t i o n con-
ditions
V {x )
S 3 = f( ) Xj = ( Xj + )/2,Vj j = 1 , . . . , N.
W h a t can we do t o o b t a i n a ( n o n - s i n g u l a r ) square system? A s we w i l l see b e l o w ,
we can a d d t h e f o l l o w i n g t h r e e c o n d i t i o n s :
N
^ c f c = 0,
fc=i
56 Meshfree Approximation Methods with M A T L A B
^CkXk = 0,
fc=l
N
y^c ?/fc = o.
fc
fc=i
How do we have t o m o d i f y o u r e x i s t i n g M A T L A B p r o g r a m for s c a t t e r e d data
i n t e r p o l a t i o n t o i n c o r p o r a t e these m o d i f i c a t i o n s ? I f we p r e v i o u s l y d e a l t w i t h t h e
solution of
Ac = y,
T
with A jk = e -e*\\ - \\^
Xj Xk j k = i , . . . c = [ C l ,... ) C 7 v] , and y =
t
[ / ( c c i ) , . . . , / ( C J V ) ] , t h e n we n o w have t o solve t h e augmented system
' A P~ c V
(6.1)
P T
O d 0
P r o g r a m 6 . 1 . R B F I n t e r p o l a t i o n 2 D l i n e a r .m
/ R B F I n t e r p o l a t i o n 2 D l i n e a r
0
% S c r i p t t h a t performs 2D RBF i n t e r p o l a t i o n w i t h r e p r o d u c t i o n of
% l i n e a r functions
% C a l l s on: D i s t a n c e M a t r i x
% Define t h e Gaussian RBF and shape parameter
1 r b f = @(e,r) e x p ( - ( e * r ) . ~ 2 ) ; ep = 6;
% Define l i n e a r t e s t f u n c t i o n
2 t e s t f u n c t i o n = @(x,y) (x+y)/2;
% Number and type of d a t a p o i n t s
5
3 N = 9; g r i d t y p e = ' u ;
% Load d a t a p o i n t s
5
4 name = s p r i n t f ( Data2D_y d%s' ,N, g r i d t y p e ) ; l o a d (name)
o
5 ctrs = dsites;
6 neval = 4 0 ; M = neval~2; g r i d = l i n s p a c e ( 0 , 1 , n e v a l ) ;
7 [[Link]] = m e s h g r i d ( g r i d ) ; e p o i n t s = [ x e ( : ) y e ( : ) ] ;
% Evaluate the t e s t function a t the data p o i n t s .
8 rhs = testfunction(dsites(:,1).dsites(:,2));
6. Scattered Data Interpolation with Polynomial Precision 57
/ Add z e r o s f o r l i n e a r (2D) r e p r o d u c t i o n
0
9 rhs = [rhs; z e r o s ( 3 , l ) ] ;
% Compute d i s t a n c e m a t r i x between t h e d a t a s i t e s and c e n t e r s
10 DM_data = D i s t a n c e M a t r i x ( d s i t e s , c t r s ) ;
% Compute i n t e r p o l a t i o n m a t r i x
11 IM = rbf(ep,DM_data);
% Define 3-column m a t r i x P f o r l i n e a r r e p r o d u c t i o n
12 PM = [ones(N.l) d s i t e s ] ;
% Augment i n t e r p o l a t i o n m a t r i x
J
13 IM = [IM PM; [PM z e r o s ( 3 , 3 ) ] ] ;
% Compute d i s t a n c e m a t r i x between e v a l u a t i o n p o i n t s and c e n t e r s
14 DM_eval = D i s t a n c e M a t r i x ( e p o i n t s , c t r s ) ;
% Compute e v a l u a t i o n m a t r i x
15 EM = rbf(ep,DM_eval);
% Add column f o r constant r e p r o d u c t i o n
16 PM = [ones(M,l) e p o i n t s ] ; EM = [EM PM];
/ Compute RBF i n t e r p o l a n t
0
% ( e v a l u a t i o n matrix * s o l u t i o n of i n t e r p o l a t i o n system)
17 Pf = EM * ( I M \ r h s ) ;
/ Compute maximum e r r o r on e v a l u a t i o n g r i d
0
18 exact = t e s t f u n c t i o n ( e p o i n t s ( : , 1 ) , e p o i n t s ( : , 2 ) ) ;
19 maxerr = n o r m ( P f - e x a c t , i n f ) ;
20 rms_err = norm(Pf-exact)/neval;
21 fprintf('RMS e r r o r : %e\n', rms_err)
5
22 fprintf('Maximum e r r o r : % e \ n , maxerr)
23 fview = [-30,30];
24 plotsurf(xe,ye,Pf,neval,exact,maxerr,fview);
25 ploterror2D(xe,ye,Pf,exact,maxerr,neval,fview);
Fig. 6.1 Top: Gaussian interpolant to bivariate linear function with N = 1089 (left) and as-
sociated abolute error (right). Bottom: Interpolant based on linearly augmented Gaussians to
bivariate linear function with N = 9 (left) and associated abolute error (right).
" A P~ c y
(6.3)
d 0
9 rhs = [rhs; 0 ] ;
12 PM = o n e s ( N , l ) ;
13 IM = [IM PM; [PM' 0 ] ] ;
16 PM = ones(M,l); EM = [EM PM];
a n d for r e p r o d u c t i o n o f b i v a r i a t e q u a d r a t i c p o l y n o m i a l s we c a n use
W e n o w need t o i n v e s t i g a t e w h e t h e r t h e a u g m e n t e d s y s t e m m a t r i x i n (6.3) is n o n -
singular. T h e special case m = 1 ( i n a n y space d i m e n s i o n s), i.e., r e p r o d u c t i o n o f
constants, is covered b y s t a n d a r d results f r o m l i n e a r algebra, a n d w e discuss i t f i r s t .
for a l l c = [ c i , . . . , C N ] T
R N
t h a t satisfy
3= 1
3= 1
T
Proof. A s s u m e [c, d] is a s o l u t i o n o f t h e homogeneous l i n e a r system, i.e., with
T T
y = 0 . W e show t h a t [c, d] = 0 is t h e o n l y possible s o l u t i o n .
T
M u l t i p l i c a t i o n o f t h e t o p b l o c k o f t h e (homogeneous) l i n e a r s y s t e m b y c yields
T T
c Ac + dc P = 0.
T T
F r o m t h e b o t t o m b l o c k o f t h e s y s t e m we k n o w P c = cP = 0, a n d therefore
T
c Ac = 0.
Since t h e m a t r i x A is c o n d i t i o n a l l y p o s i t i v e d e f i n i t e o f o r d e r one b y a s s u m p t i o n we
get t h a t c = 0 . F i n a l l y , t h e t o p b l o c k o f t h e h o m o g e n e o u s linear s y s t e m u n d e r
c o n s i d e r a t i o n states t h a t
Ac + dP = 0,
so t h a t c = 0 a n d t h e fact t h a t P is a v e c t o r o f ones i m p l y d = 0.
6. Scattered Data Interpolation with Polynomial Precision 61
for a n y N p a i r w i s e d i s t i n c t p o i n t s x \ , . . . , XN R , a n d c = [ c i , . . . , c y v ]
S T
& N
satisfying
N
^2cjP(xj) = 0,
3= 1
for any c o m p l e x - v a l u e d p o l y n o m i a l p o f degree at m o s t m 1. T h e f u n c t i o n <E> is
S
called strictly conditionally positive definite of order m on 1R. i f t h e q u a d r a t i c f o r m
(7.1) is zero o n l y for c = 0.
A n i m m e d i a t e o b s e r v a t i o n is
63
64 Meshfree Approximation Methods with M A T L A B
satisfying
N
^r p( ) Cj Xj = 0,
3= 1
for any real-valued polynomial p of degree at most m 1. The function <fr is called
s
s t r i c t l y conditionally positive definite of order m on M if the quadratic form (7.2)
is zero only for c = 0.
T
Proof. T h e p r o o f is a l m o s t i d e n t i c a l t o t h e p r o o f o f T h e o r e m 6.2. A s s u m e [c, d]
is a s o l u t i o n o f t h e homogeneous l i n e a r s y s t e m , i.e., w i t h y = 0. W e show that
T
[c, d] = 0 is t h e o n l y possible s o l u t i o n .
T
Multiplication of the top block by c yields
T T
c Ac + c Pd = 0.
T T T
F r o m t h e b o t t o m b l o c k o f (6.3) we k n o w P c 0. T h i s implies c P = 0 , and
therefore
T
c Ac = 0. (7.3)
7. Conditionally Positive Definite Functions 65
Since t h e f u n c t i o n $ is s t r i c t l y c o n d i t i o n a l l y p o s i t i v e d e f i n i t e o f o r d e r m b y a s s u m p -
T
t i o n we k n o w t h a t t h e q u a d r a t i c f o r m o f A ( w i t h coefficients such t h a t P c = 0)
above is zero o n l y for c = 0. T h e r e f o r e (7.3) tells us t h a t c = 0. T h e u n i s o l v e n c y o f
t h e d a t a sites, i.e., t h e l i n e a r independence o f t h e c o l u m n s o f P (c.f. R e m a r k 6.1),
a n d t h e fact t h a t c = 0 g u a r a n t e e d = 0 f r o m t h e t o p b l o c k
Ac + Pd = 0
of (6.3).
T h e o r e m 7.3 states t h a t s t r i c t l y c o n d i t i o n a l l y p o s i t i v e d e f i n i t e f u n c t i o n s o n W
are characterized b y t h e o r d e r o f t h e s i n g u l a r i t y o f t h e i r generalized F o u r i e r t r a n s -
f o r m at t h e o r i g i n , p r o v i d e d t h a t t h i s generalized F o u r i e r t r a n s f o r m is n o n - n e g a t i v e
a n d non-zero.
Since i n t e g r a l c h a r a c t e r i z a t i o n s s i m i l a r t o Schoenberg's T h e o r e m s 3.6 a n d 3.8
are so c o m p l i c a t e d i n t h e c o n d i t i o n a l l y p o s i t i v e d e f i n i t e case we do n o t p u r s u e t h e
concept o f a c o n d i t i o n a l l y p o s i t i v e definite r a d i a l f u n c t i o n here. The interested
reader is referred t o [ G u o et al. (1993a)] for details. W e w i l l discuss some e x a m -
ples o f r a d i a l f u n c t i o n s v i a t h e F o u r i e r t r a n s f o r m a p p r o a c h i n t h e n e x t c h a p t e r ,
a n d i n C h a p t e r 9 we w i l l e x p l o r e t h e c o n n e c t i o n between c o m p l e t e l y a n d m u l t i p l y
m o n o t o n e f u n c t i o n s a n d c o n d i t i o n a l l y p o s i t i v e definite r a d i a l f u n c t i o n s .
1
Chapter 8
8.1 E x a m p l e 1: G e n e r a l i z e d M u l t i q u a d r i c s
T h e generalized multiquadrics
2 0 s
$(x) = ( l + WxW ) , x e R, 3 e R \ N ,
0 (8.1)
have generalized F o u r i e r t r a n s f o r m s
9I+/3
/ 3 s / 2
<&M = f 7z^ii^ir - ^ / 3 + s /2(ii^ii)> # ,
67
68 Meshfree Approximation Methods with M A T L A B
Fig. 8.1 Hardy's multiquadric with 3 = | (left) and a generalized multiquadric with 3 = |
2
(right) centered at the origin in R .
8.2 E x a m p l e 2: R a d i a l P o w e r s
T h e radial powers
S
{x) = \\x\f, x e 3R , 0 < 3 2 N , (8.2)
have generalized F o u r i e r t r a n s f o r m s
20+S/2Y(S0\
= L_2j.|| ,||-/3-*
a u=AO
{ J 11 11
T(-3/2) ' ^ '
of order m = \3/2~\. Therefore, the functions
$(x) = ( - 1 ) ^ / 2 1 \\xf, 0</32N,
are s t r i c t l y c o n d i t i o n a l l y p o s i t i v e d e f i n i t e o f o r d e r m = \8/2~\ ( a n d h i g h e r ) .
T h i s shows t h a t t h e basic f u n c t i o n $(cc) = ||a?||2 used for t h e d i s t a n c e m a t r i x
fits i n t h e i n t r o d u c t o r y c h a p t e r are s t r i c t l y c o n d i t i o n a l l y p o s i t i v e d e f i n i t e o f o r d e r
one. A c c o r d i n g t o T h e o r e m 7.2 we s h o u l d have used these basic f u n c t i o n s t o g e t h e r
w i t h an a p p e n d e d c o n s t a n t . H o w e v e r , T h e o r e m 9.7 b e l o w p r o v i d e s t h e j u s t i f i c a t i o n
for t h e i r use as a p u r e distance m a t r i x .
I n F i g u r e 8.2 we show r a d i a l cubics {3 = 3, i.e., s t r i c t l y c o n d i t i o n a l l y p o s i t i v e
definite of o r d e r 2) a n d q u i n t i c s (8 = 5, i.e., s t r i c t l y c o n d i t i o n a l l y p o s i t i v e d e f i n i t e
of order 3 ) .
N o t e t h a t we h a d t o exclude even powers i n ( 8 . 2 ) . T h i s is clear since a n even
power c o m b i n e d w i t h t h e square r o o t i n t h e d e f i n i t i o n o f t h e E u c l i d e a n n o r m results
i n a p o l y n o m i a l a n d we have a l r e a d y d e c i d e d t h a t p o l y n o m i a l s c a n n o t be used
for i n t e r p o l a t i o n at a r b i t r a r i l y s c a t t e r e d m u l t i v a r i a t e sites.
N o t e t h a t r a d i a l powers are n o t affected b y a s c a l i n g o f t h e i r a r g u m e n t . I n o t h e r
w o r d s , r a d i a l powers are shape parameter free. T h i s has t h e advantage t h a t t h e
user need n o t w o r r y a b o u t f i n d i n g a " g o o d " value o f e. O n t h e o t h e r h a n d , w e w i l l
see below t h a t r a d i a l powers w i l l n o t be able t o achieve t h e s p e c t r a l convergence
rates t h a t are possible w i t h some o f t h e o t h e r basic f u n c t i o n s such as Gaussians a n d
generalized (inverse) m u l t i q u a d r i c s .
2
Fig. 8.2 Radial cubic (left) and quintic (right) centered at the origin in R .
70 Meshfree Approximation Methods with M A T L A B
8.3 E x a m p l e 3: T h i n P l a t e Splines
have generalized F o u r i e r t r a n s f o r m s
0 + 1 2 l 3 1 + s 2 s 213
= {-l) 2 - / r(B + s/2)8\\\u>\\- -
of o r d e r m = 3 + 1. T h e r e f o r e , t h e f u n c t i o n s
+ 1 2
$(x) = (-l)0 ||x|| 01og||a;||, /?GN,
are s t r i c t l y c o n d i t i o n a l l y p o s i t i v e d e f i n i t e o f o r d e r m = 3 + 1. I n p a r t i c u l a r , we can
use
N
2 2
Vf(x) = "Y^CkHx - c c | | l o g ||cc - x \\
fc k + di + d x 2 + d y, 3 x = (x,y) G R,
fc=i
together w i l l the constraints
N N N
Fig. 8.3 "Classical" thin plate spline (left) and order 3 thin plate spline (right) centered at the
2
origin in R .
2
T h e o r e m 9 . 1 . Let p> G C [ 0 , oo) n C ( 0 , o o ) . Then the function $ = p(\\ | | ) is
S
conditionally positive definite of order m and radial on R for all s if and only if
m m
(l) <>( ) is completely monotone on ( 0 , o o ) .
Proof. T h e fact t h a t c o m p l e t e m o n o t o n i c i t y i m p l i e s c o n d i t i o n a l p o s i t i v e d e f i n i t e -
ness was p r o v e d i n [ M i c c h e l l i (1986)]. M i c c h e l l i also c o n j e c t u r e d t h a t t h e converse
holds a n d gave a s i m p l e p r o o f for t h i s i n t h e case m = 1. For m = 0 t h i s is Schoen-
S
berg's c h a r a c t e r i z a t i o n o f p o s i t i v e definite r a d i a l functions o n R for a l l s i n t e r m s o f
c o m p l e t e l y m o n o t o n e f u n c t i o n s ( T h e o r e m 5.2). T h e r e m a i n i n g p a r t o f t h e t h e o r e m
is s h o w n i n [Guo et al. (1993a)].
73
74 Meshfree Approximation Methods with M A T L A B
E x a m p l e 9.1. T h e f u n c t i o n s
imply
r / 3 1 e
^)(r) = (-l) W - 1) [8 - I + 1 ) ( 1 + rf-
so t h a t
(_i)r/3i^(r/3i) ( r ) ( / 5 _ ^ + 1 ) ( i + r ) /3-r/3i
E x a m p l e 9.2. T h e f u n c t i o n s
2 2
<p(r) = ( - 1 ) ^ / V ^ , 0<8<2N,
imply
^>(r) = ( - 1 ) ^ 1 f - l ) - . . r ^ 2
^
2
so t h a t ( l ) r ^ / l ^ ( r / 3 / 2 l ) j s c o m p l e t e l y m o n o t o n e a n d m = \B/2~\ is t h e smallest
7 7 1
possible m such t h a t ( l ) ^ ) is c o m p l e t e l y m o n o t o n e . Since /? is n o t a n even
integer ip is n o t a p o l y n o m i a l , a n d therefore, t h e r a d i a l powers (c.f. ( 8 . 2 ) )
W 2 ]
*(||*||) = (-l) \\xf, 3>0, 82N,
S
are s t r i c t l y c o n d i t i o n a l l y p o s i t i v e d e f i n i t e o f o r d e r m > \8/2\ and radial on R for
a l l s.
Therefore, we l e t
w h i c h is o b v i o u s l y n o t a p o l y n o m i a l . D i f f e r e n t i a t i n g ip we get
pW(r) +1
= {-lf (3{(3 - l ) . . . ( 8 - e + iy-< logr + p {r),
e 1<<P,
w i t h pe a p o l y n o m i a l o f degree (3 1. T h e r e f o r e , t a k i n g i = (3 we have
<pW(r) = BWogr + C
and
.09+D( ) = ( - l ) / 3 + i [
r >
w h i c h is c o m p l e t e l y m o n o t o n e o n (0, o o ) .
J u s t as we showed earlier t h a t c o m p a c t l y s u p p o r t e d r a d i a l f u n c t i o n s c a n n o t b e
s
s t r i c t l y p o s i t i v e definite o n R for a l l s, i t is i m p o r t a n t t o n o t e t h a t t h e r e are n o
t r u l y c o n d i t i o n a l l y p o s i t i v e definite f u n c t i o n s w i t h c o m p a c t s u p p o r t . M o r e precisely
(see [ W e n d l a n d (2005a)]),
S
T h e o r e m 9 . 4 . Assume that the complex-valued function $ G C ( R ) has compact
support. If & is strictly conditionally positive definite of (minimal) order m, then
m is necessarily zero, i.e., $ is already strictly positive definite.
? r r
Example 9.4. T h e f u n c t i o n v fe( ) = ( 1 ~ )+ is fc-times m o n o t o n e (see E x -
a m p l e 5.5 i n S e c t i o n 5.2). To avoid the i n t e g r a t i o n constant for t h e compactly
s u p p o r t e d t r u n c a t e d p o w e r f u n c t i o n we c o m p u t e t h e a n t i - d e r i v a t i v e v i a t h e i n t e g r a l
o p e r a t o r I of E x a m p l e 5.6 i n S e c t i o n 5.2, i.e.,
oo poo / i \fc
/
Mt)dt = j ( l ~ t ) l d t = ) ~ ^ ( l - r ) k
+
+ 1
.
I f we a p p l y m - f o l d a n t i - d i f f e r e n t i a t i o n we get
(
I-Mr) = / / - W ) = ( , + 1 ) ( f c ; 2 > , , ( + m ) ( l- r ) ^ .
T h e r e f o r e , b y T h e o r e m 9.3, t h e f u n c t i o n
<p(r) = ( 1 - r ) * +
s
is s t r i c t l y c o n d i t i o n a l l y p o s i t i v e d e f i n i t e o f o r d e r m a n d r a d i a l o n R for [s/2\ <
s
k + rn 2, a n d b y T h e o r e m 9.4 i t is even s t r i c t l y p o s i t i v e d e f i n i t e a n d r a d i a l o n M. .
T h i s was also observed i n E x a m p l e 6 o f C h a p t e r 4. I n fact, we saw t h e r e t h a t <p is
s
s t r i c t l y p o s i t i v e d e f i n i t e a n d r a d i a l o n R. for [s/2j < k + m 1.
W e see t h a t we c a n c o n s t r u c t s t r i c t l y p o s i t i v e d e f i n i t e c o m p a c t l y supported
radial functions by anti-differentiating the t r u n c a t e d power f u n c t i o n . T h i s is es-
sentially the approach taken by Wendland to construct his p o p u l a r compactly
s u p p o r t e d r a d i a l basis f u n c t i o n s . W e p r o v i d e m o r e d e t a i l s o f his c o n s t r u c t i o n i n
Chapter 11.
Since a n N x N m a t r i x t h a t is c o n d i t i o n a l l y p o s i t i v e d e f i n i t e o f o r d e r one is p o s i t i v e
d e f i n i t e o n a subspace o f d i m e n s i o n N 1 i t has t h e i n t e r e s t i n g p r o p e r t y that
at least N 1 o f i t s eigenvalues are p o s i t i v e . T h i s follows i m m e d i a t e l y f r o m the
C o u r a n t - F i s c h e r t h e o r e m o f l i n e a r a l g e b r a (see e.g., p . 550 o f [ M e y e r ( 2 0 0 0 ) ] ) :
and
T
Afc = min max x Ax.
dimV=iV-fc+l <*ev
II <= | | = i
9. Conditionally Positive Definite Radial Functions 77
W i t h an a d d i t i o n a l a s s u m p t i o n o n A we c a n m a k e a n even s t r o n g e r s t a t e m e n t .
XN-I = max m i n x Ax T
> min T
c Ac > 0,
dimV=iV-l =>=V c: E--k=
II a> 11 = 1 l|c||=l
the matrix A with entries Ajk = &{xj x)k has N 1 positive and one negative
eigenvalue, and is therefore non-singular.
Proof. S i m p l y consider
N N N N N N
Yl J2 j k[$(xj
c c
- x) k + C] = E ^CjCk&ixj - Xk) + E ^CjCkC.
3=1 k=l 3=1 fc=l j=l fc=l
T h e second t e r m o n t h e r i g h t is zero since <3> is c o n d i t i o n a l l y p o s i t i v e d e f i n i t e o f
c = a n
o r d e r one, i.e., X ^ y L i j 0' d thus the statement follows.
Chapter 10
I n C h a p t e r 1 we used i n t e r p o l a t i o n w i t h d i s t a n c e m a t r i c e s as a m u l t i v a r i a t e g e n e r a l -
i z a t i o n o f t h e piecewise l i n e a r a p p r o a c h . O u r choice o f t h e distance m a t r i x a p p r o a c h
was m o t i v a t e d b y t h e fact t h a t t h e associated basis f u n c t i o n s , &j{x) = \\x X j \ \
w o u l d satisfy t h e dependence o n t h e d a t a sites i m p o s e d o n a m u l t i v a r i a t e i n t e r p o -
lation m e t h o d by the M a i r h u b e r - C u r t i s theorem. We made the (natural?) choice
o f u s i n g t h e E u c l i d e a n ( 2 - n o r m ) d i s t a n c e f u n c t i o n , a n d t h e n showed i n subsequent
chapters t h a t t h e f u n c t i o n &(x) = \\x\\2 is s t r i c t l y c o n d i t i o n a l l y p o s i t i v e d e f i n i t e
s
of o r d e r one a n d r a d i a l o n R , a n d t h u s o u r distance m a t r i x a p p r o a c h was i n d e e d
well-posed v i a M i c c h e l l i ' s T h e o r e m 9.7.
W e n o w b r i e f l y consider s o l v i n g t h e s c a t t e r e d d a t a i n t e r p o l a t i o n p r o b l e m w i t h
r a d i a l f u n c t i o n s based o n o t h e r p - n o r m s . These n o r m s are defined as u s u a l as
T h e c o n t e n t o f t h i s s e c t i o n is m o s t l y t h e s u b j e c t o f t h e p a p e r [ B a x t e r (1991)].
I f we consider o n l y distance m a t r i c e s , i.e., i n t e r p o l a t i o n m a t r i c e s g e n e r a t e d b y
t h e basic f u n c t i o n <&(cc) = ||cc|| , t h e n i t was s h o w n i n [ D y n et al.
p (1989)] t h a t
t h e choice p = 1 leads t o a s i n g u l a r m a t r i x a l r e a d y for v e r y s i m p l e sets o f d i s t i n c t
interpolation points. For e x a m p l e , i f X = { ( 0 , 0 ) , ( 1 , 0 ) , ( 1 , 1 ) , ( 0 , 1 ) } t h e n t h e 1-
n o r m distance m a t r i x is g i v e n b y
"0 1 2 1 "
10 12
2 10 1
12 10
a n d i t is easy t o v e r i f y t h a t t h i s m a t r i x is s i n g u l a r . T h i s r e s u l t has d i s c o u r a g e d
people f r o m u s i n g 1 - n o r m r a d i a l basis f u n c t i o n s .
However, i f we use, e.g., N H a l t o n p o i n t s , t h e n we have never e n c o u n t e r e d a
s i n g u l a r 1-norm distance m a t r i x i n a l l o f o u r n u m e r i c a l e x p e r i m e n t s . I n fact, t h e
79
80 Meshfree Approximation Methods with M A T L A B
6 DM = DM + a b s ( d r - c c ) . " p ;
8 DM = DM."(l/p);
T h e o r e m 1 0 . 1 . Suppose 1 < p < 2 and let A be the p-norm distance matrix with
entries
T h i s t h e o r e m is d e r i v e d f r o m a m u c h earlier t h e o r e m b y Schoenberg r e l a t i n g
c o n d i t i o n a l l y p o s i t i v e definite m a t r i c e s o f o r d e r one a n d E u c l i d e a n d i s t a n c e m a t r i -
ces. W h e n Schoenberg first s t u d i e d c o n d i t i o n a l l y p o s i t i v e d e f i n i t e m a t r i c e s o f o r d e r
one t h i s was i n c o n n e c t i o n w i t h i s o m e t r i c e m b e d d i n g s . Based o n earlier w o r k b y
K a r l M e n g e r [Menger (1928)] Schoenberg d e r i v e d t h e f o l l o w i n g result c h a r a c t e r i z i n g
c e r t a i n c o n d i t i o n a l l y p o s i t i v e definite m a t r i c e s as E u c l i d e a n distance m a t r i c e s (see
[Schoenberg (1937)]).
Ajk = I I ^ -Vk\\l-
N
These points are the vertices of a simplex in ~R .
I n t h e t h i r d r o w o f F i g u r e 10.2 we d i s p l a y t h e i n t e r p o l a n t s t o t h e test f u n c t i o n
2
f(x, y) = (x+y)/2 o n [0, l ] u s i n g distance m a t r i x i n t e r p o l a t i o n based o n 25 e q u a l l y
spaced p o i n t s a n d p - n o r m s w i t h p = 1.001 a n d p = 2. Since we use a p l a i n d i s t a n c e
i n t e r p o l a n t , i.e., $(x) = ||aj|| p i t is r e m a r k a b l e t h a t t h e e r r o r u s i n g t h e p = 1.001-
n o r m is a b o u t t w o orders o f m a g n i t u d e smaller t h a n t h e n e x t best p - n o r m d i s t a n c e
m a t r i x fit a m o n g o u r e x p e r i m e n t s ( w h i c h we o b t a i n e d for p = 100, c.f. F i g u r e 10.2).
T h e use o f different p - n o r m s for different a p p l i c a t i o n s has n o t been s t u d i e d
carefully i n t h e l i t e r a t u r e .
T w o o t h e r results r e g a r d i n g i n t e r p o l a t i o n w i t h p - n o r m r a d i a l basis f u n c t i o n s
can also be f o u n d i n t h e l i t e r a t u r e . I n [ W e n d l a n d (2005a)] we find a reference t o
[ Z a s t a v n y i (1993)] a c c o r d i n g t o w h i c h for space dimensions s > 3 t h e o n l y
s
f u n c t i o n t h a t is p o s i t i v e d e f i n i t e a n d p - n o r m r a d i a l o n M is t h e zero f u n c t i o n .
A g a i n , s o m e w h a t d i s c o u r a g i n g news. H o w e v e r , t h e r e is also g o o d news. T h e f o l l o w -
i n g r a t h e r p o w e r f u l t h e o r e m comes f r o m [ B a x t e r (1991)]. B a x t e r calls t h e m a t r i x
A o f T h e o r e m 10.2 a n almost negative definite m a t r i x (c.f. the remarks following
D e f i n i t i o n 6.2).
r a
Proof. B y Schoenberg's T h e o r e m 10.2 w e c a n w r i t e Ajk = \\yj f P"
N 2
p r o p r i a t e p o i n t s yj ~R B y a s s u m p t i o n p{ ) is c o n d i t i o n a l l y p o s i t i v e d e f i n i t e o f
o r d e r one a n d r a d i a l , a n d therefore B is c o n d i t i o n a l l y p o s i t i v e s e m i - d e f i n i t e o f o r d e r
one. M o r e o v e r , i f Ajk ^ 0 for a l l o f f - d i a g o n a l elements, t h e n j / i , . . . , ? / A T are d i s t i n c t ,
2
a n d therefore B is s t r i c t l y c o n d i t i o n a l l y p o s i t i v e d e f i n i t e o f o r d e r one p r o v i d e d (p( )
is s t r i c t l y c o n d i t i o n a l l y p o s i t i v e d e f i n i t e o f o r d e r one.
2
Fig. 10.1 p-norm Gaussians for p = 1 (left) and p = 10 (right) centered at the origin in R .
A n o t h e r , closely r e l a t e d t h e o r e m b y B a x t e r is
2 2
T h e o r e m 1 0 . 4 . Suppose p( ) and ip( ) are functions that are conditionally pos-
S 2
itive definite of order one and radial on ~R with ip(0) = 0. Then ip o ip( ) is also
S 2
conditionally positive definite of order one and radial on 1R . Indeed, if ip(- ) is
strictly conditionally positive definite of order one and radial and tp vanishes only
2
at zero, then ip o ip( ) is strictly conditionally positive definite of order one and
radial.
T h i s t h e o r e m is a g e n e r a l i z a t i o n o f a classical r e s u l t i n l i n e a r a l g e b r a b y Schur
(see, e.g., [ H o r n a n d J o h n s o n ( 1 9 9 1 ) ; M i c c h e l l i ( 1 9 8 6 ) ] , w h e r e Schur's t h e o r e m was
e x t e n d e d t o cover s t r i c t n e s s ) .
10. Miscellaneous Theory: Other Norms and Scattered Data Fitting on Manifolds 83
T h e r e exists a sizeable b o d y o f l i t e r a t u r e o n t h e t o p i c o f s c a t t e r e d d a t a i n t e r p o l a t i o n
- 1 s
on m a n i f o l d s , especially t h e sphere S ^ i n M . W e w i l l n o t m e n t i o n a n y specific
results here. I n s t e a d w e refer t h e reader t o t h e b o o k [Freeden et al. (1998)], t h e
survey papers [Cheney (1995a); Fasshauer a n d S c h u m a k e r (1998)], as w e l l as m a n y
o r i g i n a l papers such as [ B a x t e r a n d H u b b e r t (2001); B i n g h a m (1973); Fasshauer
(1995a); Fasshauer (1999b); H u b b e r t a n d M o r t o n (2004a); H u b b e r t a n d M o r t o n
(2004b); Levesley et al. (1999); M e n e g a t t o ( 1 9 9 4 b ) ; N a r c o w i c h a n d W a r d ( 2 0 0 2 ) ;
R a g o z i n a n d Levesley (1996); R o n a n d S u n ( 1 9 9 6 ) ; Schoenberg (1942); Schreiner
(1997); W a h b a (1981); W a h b a (1982); X u a n d C h e n e y ( 1 9 9 2 b ) ] .
R a d i a l basis functions o n m o r e general R i e m a n n i a n m a n i f o l d s are s t u d i e d i n ,
e.g., [ D y n et al. (1997); D y n et al. (1999); Levesley a n d R a g o z i n (2002); N a r c o w i c h
(1995); N a r c o w i c h et al. (2003); S c h i m m i n g a n d Belger (1991)].
T h e r e is also a " p o o r m a n ' s s o l u t i o n " t o i n t e r p o l a t i o n o n m a n i f o l d s , especially
the sphere. O n e can use t h e E u c l i d e a n r a d i a l basis f u n c t i o n m e t h o d s discussed t h u s
far, a n d s i m p l y r e s t r i c t t h e i r e v a l u a t i o n t o t h e m a n i f o l d . T h i s a p p r o a c h is o u t l i n e d
i n Section 6 o f [Fasshauer a n d S c h u m a k e r (1998)].
W e w i l l discuss a n o t h e r , r e l a t e d , i n t e r p o l a t i o n p r o b l e m l a t e r . N a m e l y , i n t e r p o -
3
l a t i o n t o p o i n t c l o u d d a t a i n R . I n t h i s case, t h e u n d e r l y i n g m a n i f o l d is u n k n o w n ,
a n d a n o t h e r a p p r o a c h needs t o be t a k e n . See C h a p t e r 30 for details.
10.3 Remarks
M a n y o f t h e results g i v e n i n t h e p r e v i o u s c h a p t e r s c a n be generalized t o v e c t o r -
v a l u e d or even m a t r i x - v a l u e d f u n c t i o n s . Some results a l o n g these lines c a n be f o u n d
i n [ L o w i t z s c h (2002); L o w i t z s c h (2005); M y e r s (1992); N a r c o w i c h a n d W a r d (1994a);
Schaback (1995a)].
We point out that the approach to solving the interpolation problems taken i n
t h e p r e v i o u s chapters a l w a y s assumes t h a t t h e centers, i.e., the points x , k k =
1,...,N, at w h i c h t h e basis f u n c t i o n s are centered, coincide w i t h t h e d a t a sites.
T h i s is a f a i r l y severe r e s t r i c t i o n , a n d i t has been s h o w n i n examples i n t h e c o n t e x t
o f least squares a p p r o x i m a t i o n o f s c a t t e r e d d a t a (see, e.g., [Franke et al. (1994);
F r a n k e et al. (1995)] or [Fasshauer (1995a)]) t h a t b e t t e r r e s u l t s can be achieved
i f t h e centers are chosen different f r o m t h e d a t a sites. Theoretical results i n this
d i r e c t i o n are v e r y l i m i t e d , a n d are r e p o r t e d i n [ Q u a k et al. (1993)] a n d i n [Sun
(1993a)].
84 Meshfree Approximation Methods with M A T L A B
z
0.5
zo.5
2
0.5-
2
Fig. 10.2 p-norm distance matrix fits to f(x, y) = (x + y)/2 on a 5 X 5 grid in [0, l ] unless noted
otherwise. Top: p = 1 (1089 Halton points). 2nd row: p = 10 (left), p = 100 (right). 3rd row:
p = 1.001 (left), p = 2 (right). Bottom: p-norm Gaussian fits for p = 1 (left) and p = 10 (right).
Chapter 11
Compactly Supported
Radial Basis Functions
/>oo
2
*(x) = F,<p{\\x\\) = Ija-H-C/ <p(t)r' J - (t\\x\\)dt,
(a 2)/2
Jo
11.1 O p e r a t o r s for R a d i a l F u n c t i o n s a n d D i m e n s i o n W a l k s
85
86 Meshfree Approximation Methods with M A T L A B
Definition 11.1.
2
(2) For even <p G C(M) we define t h e differential operator V via
N o t e t h a t t h e i n t e g r a l o p e r a t o r X differs f r o m t h e o p e r a t o r / i n t r o d u c e d earlier
(see (5.1)) b y a f a c t o r t i n t h e i n t e g r a n d .
T h e m o s t i m p o r t a n t p r o p e r t i e s o f t h e m o n t e e a n d descente o p e r a t o r s are (see,
e.g., [Schaback a n d W u (1996)] o r [ W e n d l a n d ( 1 9 9 5 ) ] ) :
Theorem 11.1.
(1) Both T> andX preserve compact support, i.e., if if has compact support, then so
do T>p and Xp.
(2) IfpE C(R) and t ^ t(f>(t) G L i [ 0 , o o ) , then VX<p = (p.
2
(3) Ifpe C (R) (<p^l) is even and <p' G L [0, o o ) , then XVp x = ip.
s x
(4) J / t H t ~ p{t) G L i [ 0 , o o ) and s > 3, then F (<p) s = JF _ (X^).
S 2
2 s
(5) Ifpe C (R) is even and 11-+ t <p'(t) G L i [ 0 , o o ) , then F (ip) 8 = T i{pp).
s+
T h e o p e r a t o r s X a n d V a l l o w us t o express s - v a r i a t e F o u r i e r t r a n s f o r m s as (s2)-
o r ( s + 2 ) - v a r i a t e F o u r i e r t r a n s f o r m s , respectively. I n p a r t i c u l a r , a d i r e c t consequence
o f t h e above p r o p e r t i e s a n d t h e c h a r a c t e r i z a t i o n o f s t r i c t l y p o s i t i v e d e f i n i t e r a d i a l
f u n c t i o n s ( T h e o r e m 3.6) is
Theorem 11.2.
s l
(1) Suppose cp G C(R). Ift ^ t ~ p(t) G Z a [ 0 , o o ) and s > 3, then <p is strictly
s
positive definite and radial on R if and only if X<p is strictly positive definite
s 2
and radial on R~.
2 s
(2) If p e C (R) is even and t H- t <p'(t) L i [ 0 , o o ) , then p is strictly positive
s
definite and radial on R if and only if T>p is strictly positive definite and
s + 2
radial on R .
T h i s a l l o w s us t o c o n s t r u c t n e w s t r i c t l y p o s i t i v e d e f i n i t e r a d i a l f u n c t i o n s from
g i v e n ones b y a " d i m e n s i o n - w a l k " t e c h n i q u e t h a t steps t h r o u g h m u l t i v a r i a t e E u -
c l i d e a n space i n even i n c r e m e n t s . T h e e x a m p l e s p r e s e n t e d i n t h e f o l l o w i n g sections
illustrate this technique.
.1,
11. Compactly Supported Radial Basis Functions 87
e
Definition 11.2. W i t h <pe(r) = ( 1 r) + we define
k(
Vs,k = 1 P\_s/2\+k+l-
s
T h e o r e m 1 1 . 3 . The functions (p kSj are strictly positive definite and radial on R
and are of the form
s k K
^ ' ' \0, r > 1,
2 f e
with a univariate polynomial pk 3t of degree [s/2\ + Sk + 1. Moreover, (p k
Sj C (IR)
are unique up to a constant factor, and the polynomial degree is minimal for given
space dimension s and smoothness 2k.
2k
T h i s t h e o r e m states t h a t a n y o t h e r c o m p a c t l y s u p p o r t e d C p o l y n o m i a l func-
t i o n t h a t is s t r i c t l y p o s i t i v e d e f i n i t e a n d r a d i a l o n W w i l l n o t have a smaller p o l y -
n o m i a l degree. O u r o t h e r e x a m p l e s b e l o w ( W u ' s f u n c t i o n s , G n e i t i n g ' s f u n c t i o n s )
i l l u s t r a t e t h i s fact. T h e s t r i c t p o s i t i v e definiteness of Wendland's functions ip k
Sj
+ 1
>,,i(r) = ( l - r [(* + l ) r + l ] ,
e 2 2 2
Vs,2{r) = ( 1 - r) + [( + 4 + 3 ) r + (3 + 6 ) r + 3] ,
3 3 2 3 2 2
ip (r)
at3 = (1 - rY+ [{ + 9 + 23 + 1 5 ) r + {U + 3Q + 4 5 ) r
+ (15^ + 4 5 ) r + 1 5 ] ,
oo
- l
t(l -
t) dt +
e
= f t(l-t) dt
J r
+ 1
= ( , + 1 )(, + 2 ) ( i - r ) ' [ ( ^ i ) r + i ] ,
w h e r e t h e c o m p a c t s u p p o r t o f <f reduces t h e i m p r o p e r i n t e g r a l t o a d e f i n i t e i n t e g r a l
w h i c h c a n be e v a l u a t e d u s i n g i n t e g r a t i o n b y p a r t s . T h e o t h e r t w o cases are o b t a i n e d
similarly b y repeated application o f X.
a n
Table 11.1 Wendland's compactly supported radial functions v'[Link] d *?s,k Vs,fc(l " )
for various choices offcand s 3.
k <P3,k( ) r
&3,k(r) smoothness
r C
0 i
2
1 (1 - r)\ (4r + 1) r\ (5 - 4 r ) C
2 2 4
2 (1 - r)\ ( 3 5 r + 18r + 3) r\ (56 - 88r + 3 5 r ) C
3 2 2 3 6
3 (1 - r)\ ( 3 2 r + 2 5 r + 8r + l ) r \ (66 - 154r + 1 2 1 r - 3 2 r ) C
I n [ W u (1995b)] w e f i n d a n o t h e r w a y t o c o n s t r u c t s t r i c t l y p o s i t i v e d e f i n i t e r a d i a l
functions w i t h compact support. W u starts w i t h t h e f u n c t i o n
2 e
i>(r) = (1 - r ) , + eN,
11. Compactly Supported Radial Basis Functions 89
-oo
2
= f ( l - t Y { l - { 2 r - t Y Y + d t .
T h i s f u n c t i o n is s t r i c t l y p o s i t i v e d e f i n i t e since i t s F o u r i e r t r a n s f o r m is essentially
the square o f t h e F o u r i e r t r a n s f o r m o f ip a n d therefore n o n - n e g a t i v e . J u s t like t h e
W e n d l a n d functions, t h i s f u n c t i o n is a p o l y n o m i a l o n i t s s u p p o r t . I n fact, t h e degree
2e
of t h e p o l y n o m i a l is 4 + 1, a n d ip e e C (R).
N o w , a f a m i l y o f s t r i c t l y p o s i t i v e d e f i n i t e r a d i a l f u n c t i o n s is c o n s t r u c t e d b y a
d i m e n s i o n w a l k u s i n g t h e T> o p e r a t o r .
2 2
D e f i n i t i o n 11.3. W i t h ip {r)e = ( ( 1 - - Y+ * ( 1 - - ) + ) ( 2 r ) we define
k
iP t K = V ip . t
s
T h e f u n c t i o n s ipk,e are s t r i c t l y p o s i t i v e d e f i n i t e a n d r a d i a l o n R for s < 2k + 1,
2 k
are p o l y n o m i a l s o f degree A 2k+l o n t h e i r s u p p o r t a n d i n C ^~ ^ i n the interior
2 k
of t h e s u p p o r t . O n t h e b o u n d a r y t h e smoothness increases t o C ~.
E x a m p l e 1 1 . 2 . For = 3 we c a n c o m p u t e t h e four f u n c t i o n s
k k 2 2
W ) = V ip {r) 3 = V ((l - -f + * (1 - - ) ) ( 2 r ) , + k = 0,1,2,3.
Table 11.2 Wu's compactly supported radial functions ipk,e for various choices of
fc and I = 3.
fc smoothness s
0 2
( l - r ) + ( 5 - r 35r + 101r + 147r + 101r + 3 5 r + 5 r )
3 4 5 6
1
1 2
( 1 - r ) ( 6 + 36r + 8 2 r + 7 2 r + 3 0 r + 5 r )
+
3 4 5
c 4
3
2
2 (1 - r ) + ( 8 + 40r + 4 8 r + 2 5 r + 5 r )
2 3 4
c 5
3 4
(1 - r ) (16 + 29r + 2 0 r + 5 r ) 2 3
c 7
90 Meshfree Approximation Methods with M A T L A B
Table 11.3 Shifted version ipk,e of Wu's compactly supported radial functions tpk,e
for various choices of k and = 3.
r
k ^k,zi ) smoothness s
2 3 4 5 6
0 r ( 4 2 9 - 1287r + 1573r - l O O l r + 3 5 1 r - 6 5 r + 5 r )
+ C 6
1
2 3 4 5
1 r ( 2 3 1 - 561r + 528r - 242r + 5 5 r - 5 r )
+ C 4
3
2 3 4
2 r ( 1 2 6 - 231r + 153r - 4 5 r + 5 r )
+ C 2
5
4 2 3
3 r (70 - 84r + 3 5 r - 5 r ) C 7
Fig. 11.1 Plot of Wendland's functions from Example 11.1 (left) and Wu's functions from E x a m -
ple 11.2 (right).
O t h e r s t r i c t l y p o s i t i v e d e f i n i t e c o m p a c t l y s u p p o r t e d r a d i a l f u n c t i o n s have b e e n p r o -
posed b y G n e i t i n g (see, e.g., [ G n e i t i n g ( 2 0 0 2 ) ] ) . H e s h o w e d t h a t a f a m i l y o f o s c i l l a -
t o r y c o m p a c t l y s u p p o r t e d f u n c t i o n s c a n be c o n s t r u c t e d u s i n g t h e so-called turning
11. Compactly Supported Radial Basis Functions 91
bands operator o f [ M a t h e r o n ( 1 9 7 3 ) ] . S t a r t i n g w i t h a f u n c t i o n ip s t h a t is s t r i c t l y
s
p o s i t i v e definite a n d r a d i a l o n M. for s > 3 t h e t u r n i n g b a n d s o p e r a t o r p r o d u c e s
^ - ( r ) = ^ (r) + ^i^
2 s (11.1)
s _ 2
w h i c h is s t r i c t l y p o s i t i v e d e f i n i t e a n d r a d i a l o n R .
( \ fi V (-\ ^ o ( l + l ) ( l + 2 + s) 2 \
s
w h i c h are s t r i c t l y p o s i t i v e d e f i n i t e a n d r a d i a l o n M p r o v i d e d > (see [ G n e i t i n g
(2002)]). Some specific f u n c t i o n s f r o m t h i s f a m i l y are l i s t e d i n T a b l e 11.4. A l l o f
2
t h e f u n c t i o n s are i n C ( 1 R ) . I f we w a n t s m o o t h e r f u n c t i o n s , t h e n we need t o s t a r t
w i t h a s m o o t h e r W e n d l a n d f a m i l y as d e s c r i b e d b e l o w i n E x a m p l e 11.4.
T r
^ 2,e( ) smoothness
7 i | 5
7/2 (1 - r) ^ (1 + \r - r 2 )
C 2
2 2
5 (1 - r)\ (1 + 5r - 2 7 r ) C
15/2 ( l - r ^ l + f r - _3|i r 2 ) C 2
12 (i - 0 + ( i + 1 2 r
- 104r )
2
C 2
Fig. 11.2 Oscillatory functions of Table 11.4 (left) and Table 11.5 (right).
92 Meshfree Approximation Methods urith M A T L A B
E x a m p l e 1 1 . 4 . A l t e r n a t i v e l y , we c a n o b t a i n a set o f o s c i l l a t o r y f u n c t i o n s t h a t are
3
s t r i c t l y positive definite and r a d i a l o n M b y applying the t u r n i n g bands operator
t o t h e W e n d l a n d f u n c t i o n s </?5,fc t h a t are s t r i c t l y p o s i t i v e d e f i n i t e a n d r a d i a l o n R 5
k Ok(r) smoothness
1 (1 - r)\ (1 + 4 r - 1 5 r ) 2
C 2
2 3 4
2 (1 - r)\ (3 + 18r + 3 r - 1 9 2 r ) C
2 3 4 6
3 (1 - r ) \ (15 + 120r + 2 1 0 r - 8 4 0 r - 3 4 6 5 r ) C
f 27ry 2 f c _ (2r)-r(l-r )
1
2 f c
r. _ -i o q
^ f c + 1 (2r) = K - 1 , 2 , 6 , . . . ,
12(1-r) fc = 0,
for o d d space d i m e n s i o n s s = 2k + 1, a n d as
[ 2 ry 7 2 f c (2r)-TV(l-r2)(l-r ) 2 f c
, _
>2fc+ (2r) = I
2 ^+2 * - 1, A 6,
2 ( a r c c o s r ry/l r) 2
k = 0,
11. Compactly Supported Radial Basis Functions 93
Table 11.6 Euclid's hat functions (defined for 0 < r < 2) for
different values of s.
s V>s(r) smoothness
1 1_ r C
1
2
2
2 ^ ^4arccos (^) r \ / 4 r ^ C
3 1
- sh (( + * ) ~ ) 4 1 6 r r3 c
4 I arccos ( r ) _ _1_^4 - r (20r + r ) 2 3
c
5 1
- 64^ ( t 1 2
+ 8 7 r
+ 3 2 7 r 2
) r 3
~ ( + 2 ? r
) r 3
) c
Fig. 11.3 Euclid's hat functions (left) of Table 11.6 and Buhmann's function of Example 11.6
(right).
A n o t h e r c o n s t r u c t i o n d e s c r i b e d i n [Schaback (1995a)] is t h e r a d i a l i z a t i o n o f t h e
s-fold tensor p r o d u c t o f u n i v a r i a t e 5 - s p l i n e s o f even o r d e r 2m w i t h u n i f o r m k n o t s .
T h e s e f u n c t i o n s d o n o t seem t o have a s i m p l e r e p r e s e n t a t i o n t h a t lends i t s e l f t o
numerical computations. A s c a n be seen f r o m i t s r a d i a l i z e d F o u r i e r t r a n s f o r m , t h e
5
r a d i a l i z e d S - s p l i n e i t s e l f is n o t s t r i c t l y p o s i t i v e d e f i n i t e a n d r a d i a l o n a n y R with
s > 1. For s = 1 o n l y t h e S - s p l i n e s o f even o r d e r are s t r i c t l y p o s i t i v e d e f i n i t e (see,
e.g., [ S c h o l k o p f a n d S m o l a ( 2 0 0 2 ) ] ) .
T h e last f a m i l y o f c o m p a c t l y s u p p o r t e d s t r i c t l y p o s i t i v e d e f i n i t e r a d i a l f u n c t i o n s
we w o u l d like t o m e n t i o n is due t o [ B u h m a n n ( 1 9 9 8 ) ] . B u h m a n n ' s f u n c t i o n s c o n t a i n
94 Meshfree Approximation Methods with M A T L A B
a l o g a r i t h m i c t e r m i n a d d i t i o n t o a p o l y n o m i a l . H i s f u n c t i o n s have t h e g e n e r a l f o r m
/OO
Jo
H e r e 0 < 5 < ^, p > 1, a n d i n o r d e r t o o b t a i n f u n c t i o n s t h a t are s t r i c t l y p o s i t i v e
s
d e f i n i t e a n d r a d i a l o n M for s < 3 t h e c o n s t r a i n t s for t h e r e m a i n i n g p a r a m e t e r s are
A > 0, a n d - 1 < a < =.
W h i l e i t is s t a t e d i n [ B u h m a n n (2000)] t h a t t h e c o n s t r u c t i o n t h e r e encompasses
b o t h W e n d l a n d ' s a n d W u ' s f u n c t i o n s , a n even m o r e g e n e r a l t h e o r e m t h a t shows t h a t
integration of a positive function / L i [ 0 , o o ) against a s t r i c t l y positive definite
k e r n e l K r e s u l t s i n a s t r i c t l y p o s i t i v e d e f i n i t e f u n c t i o n c a n be f o u n d i n [ W e n d l a n d
(2005a)] (see also S e c t i o n 4 . 8 ) . M o r e specifically,
roc
p(r) = / K(t,r)f(t)dt
Jo
is s t r i c t l y p o s i t i v e d e f i n i t e . B u h m a n n ' s c o n s t r u c t i o n t h e n c o r r e s p o n d s t o c h o o s i n g
a s 2
f(t) = t ( l - t )^ a n d K(t, r) = ( 1 - r /t)\.
Chapter 12
We n o w have a n a l t e r n a t i v e w a y t o c o n s t r u c t a n R B F i n t e r p o l a n t t o scattered
s
d a t a i n M. . I f we use t h e c o m p a c t l y s u p p o r t e d r a d i a l f u n c t i o n s o f t h e p r e v i o u s
chapter t h e n t h e m a i n difference t o o u r p r e v i o u s i n t e r p o l a n t s is t h a t n o w t h e i n -
t e r p o l a t i o n m a t r i x c a n be m a d e sparse b y s c a l i n g t h e s u p p o r t o f t h e basic f u n c t i o n
a p p r o p r i a t e l y . T o achieve t h i s w e use as w e d i d earlier t h e basic f u n c t i o n s
(p (r)
e = <p{er). T h u s , a large v a l u e o f e corresponds t o a s m a l l s u p p o r t . I n o t h e r
w o r d s , i f t h e s u p p o r t o f <p is t h e i n t e r v a l [ 0 , 1 ] , t h e n t h e s u p p o r t r a d i u s p o f <p is
We have s t r u c t u r e d t h e s c a t t e r e d d a t a i n t e r p o l a t i o n p r o g r a m i n t h e c o m p a c t l y sup-
p o r t e d case analogous t o t h e code for t h e g l o b a l i n t e r p o l a n t s , i.e., first c o n s t r u c t a
distance m a t r i x , a n d t h e n a p p l y t h e a n o n y m o u s f u n c t i o n r b f t o o b t a i n t h e i n t e r p o -
l a t i o n / e v a l u a t i o n m a t r i x (as o n lines 1 3 - 1 4 a n d 1 5 - 1 6 o f P r o g r a m 2.1). H o w e v e r , i t
t u r n s o u t t h a t i t is easier t o d e a l w i t h t h e c o m p a c t s u p p o r t i f we c o m p u t e t h e "dis-
tance m a t r i x " c o r r e s p o n d i n g t o t h e ( 1 s r ) + t e r m since o t h e r w i s e those entries o f
t h e distance m a t r i x t h a t are zero (since t h e m u t u a l distance between t w o i d e n t i c a l
p o i n t s is zero) w o u l d be "lost" i n t h e sparse r e p r e s e n t a t i o n o f t h e m a t r i x .
T h e M A T L A B code DistanceMatrixCSRBF.m ( P r o g r a m 12.1) c o n t a i n s t w o s i m i -
95
96 Meshfree Approximation Methods with M A T L A B
lar blocks that will be used depending on whether we have more centers than data
sites or vice versa. For example, if there are more data sites than centers (cf. lines 7
16), then we build a kd-txee for the data sites and find for each center x - those 3
data sites within the support of the basis function centered at Xj, i.e., we construct
the (sparse) matrix column by column. In the other case (cf. lines 18-27) we start
with a tree for the centers and build the matrix row by row. This is accomplished by
determining for each data site Xi all centers whose associated^ basis function
covers data site Xi.
The functions kdtree and kdrangequery are provided by the kd-tree library
mentioned above. T h e call in line 7 (respectively 18) of Program 12.1 generates the
kd-txee of all the centers (data sites), and with the call to kdrangequery in line 9
(respectively 20) we find all centers (data sites) that lie within a distance support
of the jth center point (data site). T h e actual distances are returned in the vector
d i s t and the indices into the list of all data sites are provided in idx. T h e distances
for these points only are stored in the matrix DM. For maximum efficiency (in order
to avoid dynamic memory allocation) it is important to have a good estimate of
the number of nonzero entries in the matrix for the allocation statement in lines 4
and 5. The version of the code presented here has the best performance for larger
problems since s p a r s e is only invoked once.
P r o g r a m 12.1. DistanceMatrixCSRBF.m
% DM = D i s t a n c e M a t r i x C S R B F ( d s i t e s , c t r s , e p )
% Forms t h e d i s t a n c e m a t r i x of two s e t s of p o i n t s i n R~s
% f o r compactly supported r a d i a l b a s i s f u n c t i o n s , i . e . ,
7. DM(i,j) = I I d a t a s i t e _ i - c e n t e r _ j I I _2.
7o The CSRBF used with t h i s code must be g i v e n i n s h i f t e d form
7. rbf2(u) = r b f ( r ) , u=l-e*r.
% F o r example, t h e Wendland C2
7o rbf = @(e,r) max(l-e*r , 0 ) . ~ 4 . * ( 4 * e * r + l ) ;
7 becomes
7. r b f 2 = ffl(u) u. "4.*(4*u+5) ;
7o Input
% d s i t e s : Nxs matrix r e p r e s e n t i n g a s e t of N d a t a s i t e s
% i n R"s ( i . e . , each row c o n t a i n s one
7o s-dimensional p o i n t )
% ctrs: Mxs matrix r e p r e s e n t i n g a s e t of M c e n t e r s f o r
7o RBFs i n R~s ( a l s o one c e n t e r p e r row)
7 ep: determines s i z e of support of b a s i s f u n c t i o n .
7o Small ep y i e l d s wide f u n c t i o n ,
7o i . e . , s u p p o r t s i z e = 1/ep
7. Output
7, DM: NxM SPARSE m a t r i x t h a t c o n t a i n s t h e E u c l i d e a n
12. Interpolation with Compactly Supported RBFs in M A T L A B 97
/
0 u - d i s t a n c e ( u = l - e * r ) between t h e i - t h d a t a
% s i t e and t h e j - t h c e n t e r i n t h e i , j p o s i t i o n
% Uses: k-D t r e e package by Guy S h e c h t e r from
'/, MATLAB C e n t r a l F i l e Exchange
1 f u n c t i o n DM = D i s t a n c e M a t r i x C S R B F ( d s i t e s , c t r s , e p )
2 N = size(dsites,1); M = size(ctrs,1);
7 B u i l d k-D t r e e f o r d a t a s i t e s
% F o r each c e n t e r ( b a s i s f u n c t i o n ) , f i n d t h e d a t a s i t e s
% i n i t s support along w i t h u - d i s t a n c e
3 support = 1/ep;
4 nzmax = 25*N; rowidx = zeros(1,nzmax); c o l i d x = zeros(1,nzmax);
5 v a l i d x = zeros(1,nzmax); i s t a r t = 1; i e n d = 0;
6 i f M > N / f a s t e r i f more c e n t e r s than d a t a s i t e s
7 [tmp,tmp,Tree] = k d t r e e ( c t r s , [] ) ;
8 f o r i = 1:N
9 [pts,dist,idx] = kdrangequery(Tree,dsites(i,:),support);
10 newentries = l e n g t h ( i d x ) ;
11 iend = iend + n e w e n t r i e s ;
12 rowidx(istart:iend) = repmat(i,1,newentries);
13 c o l i d x ( i s t a r t : i e n d ) = idx';
14 validx(istart:iend) = l-ep*dist';
15 i s t a r t = i s t a r t + newentries;
16 end
17 e l s e
18 [tmp,tmp,Tree] = k d t r e e ( d s i t e s , [ ] ) ;
19 f o r j = 1:M
20 [[Link],idx] = kdrangequery(Tree,ctrs(j,:),support);
21 newentries = l e n g t h ( i d x ) ;
22 iend = iend + n e w e n t r i e s ;
23 rowidx(istart:iend) = idx';
24 c o l i d x ( i s t a r t : i e n d ) = repmat(j,[Link]);
25 validx(istart:iend) = l-ep*dist';
26 i s t a r t = i s t a r t + newentries;
27 end
f
28 end
29 idx = f i n d ( r o w i d x ) ;
30 DM = s p a r s e ( r o w i d x ( i d x ) , c o l i d x ( i d x ) , v a l i d x ( i d x ) , N , M ) ;
7. Free t h e k-D Tree from memory.
31 kdtree ( [ ] , , Tree) ;
T h e reason for c o d i n g D i s t a n c e M a t r i x C S R B F . m i n t w o different w a y s is so t h a t
we w i l l be able t o speed u p t h e p r o g r a m w h e n d e a l i n g w i t h n o n - s q u a r e ( e v a l u a t i o n )
m a t r i c e s (for e x a m p l e i n t h e c o n t e x t o f M L S a p p r o x i m a t i o n (c.f. Chapter 24).
-3*, *
98 Meshfree Approximation Methods with M A T L A B
1 f u n c t i o n DM = D i s t a n c e M a t r i x C S R B F ( d s i t e s , c t r s , e p )
2 N = size(dsites,l); M = size(ctrs,1);
% B u i l d k-D t r e e f o r d a t a sites
% F o r each c e n t e r ( b a s i s f u n c t i o n ) , f i n d t h e d a t a sites
% i n i t s support along w i t h u - d i s t a n c e
3 support = 1/ep; nzmax = 25*N; DM = spalloc(N,M,nzmax);
4 [tmp, tmp,Tree] = k d t r e e ( d s i t e s , [ ] ) ;
5 f o r j = 1:M
6 Cpts,dist,idx] = kdrangequery(Tree,ctrs(j,:),support);
7 DM(idx.j) = l - e p * d i s t ;
8 end
% F r e e the k-D Tree from memory.
9 kdtree ( [ ] , [ ] , Tree) ;
13 DM_data = D i s t a n c e M a t r i x C S R B F ( d s i t e s , c t r s , e p ) ;
15 DM_eval = D i s t a n c e M a t r i x C S R B F ( e p o i n t s , c t r s , e p ) ;
2
a n d t o define t h e R B F i n s h i f t e d f o r m , i.e., i n s t e a d o f r e p r e s e n t i n g , e.g., t h e C
W e n d l a n d f u n c t i o n c / ^ i o n line 1 b y
1 r b f = @(e,r) m a x ( l - e * r , 0 ) . ~ 4 . * ( 4 * e * r + l ) ; ep=0.7;
we n o w w r i t e
17 c = p c g ( I M , r h s ) ; Pf = EM * c;
e
r a t '"fa--/ *> fc = 2,3,..., (12.1)
ln(/ifc_i/ftfc)
9 1.562729e-001 0.03
25 2.807706e-002 2.4766 0.04
81 4.853006e-003 2.5324 0.12
289 2.006041e-004 4.5965 0.45
1089 1.288000e-005 3.9611 2.75
4225 1.382497e-006 3.2198 47.92
9 1.655969e-001 0.03
25 3.097850e-002 2.4183 0.06
81 4.612941e-003 2.7475 0.20
289 1.305297e-004 5.1432 0.72
1089 4.780575e-006 4.7711 4.06
4225 2.687479e-007 4.1529 55.09
i n t h e sparse s e t t i n g a n d as
S
D e f i n i t i o n 1 3 . 1 . L e t H be a r e a l H i l b e r t space o f f u n c t i o n s / : Q(Q R ) > R w i t h
i n n e r p r o d u c t (-, - ) ^ - A f u n c t i o n K : Q x ft R is called reproducing kernel for Ji
if
\8 f\
x = \f(x)\<--M\\f\\ n
103
104 Meshfree Approximation Methods with M A T L A B
K(x,y) = (K(;y),K(;x)) n
\f(x) - f {x)\
n = \(f - f ,K(-,x)) \
n n < ||/ - f \\ \\K(;x)\\ .
n n n D
N o w i t is i n t e r e s t i n g for us t h a t t h e r e p r o d u c i n g k e r n e l K is k n o w n t o b e p o s i t i v e
definite. Here w e use a s l i g h t g e n e r a l i z a t i o n o f t h e n o t i o n o f a p o s i t i v e d e f i n i t e func-
t i o n t o a p o s i t i v e d e f i n i t e k e r n e l . Essentially, w e replace &(xj X k ) i n D e f i n i t i o n 3.2
b y K(xj, X k ) . A t t h i s p o i n t w e r e m i n d t h e reader t h a t t h e space o f b o u n d e d l i n e a r
f u n c t i o n a l s o n 7i is k n o w n as i t s dual, a n d d e n o t e d b y 7i*.
j=l k=l
N
2
= ||^c -K(-,x )|| ,>0.
J J
j=i
T h u s K is p o s i t i v e d e f i n i t e . T o e s t a b l i s h t h e second c l a i m w e assume K is n o t
s t r i c t l y p o s i t i v e d e f i n i t e a n d show t h a t t h e p o i n t e v a l u a t i o n f u n c t i o n a l s m u s t b e
13. Reproducing Kernel Hilbert Spaces for Strictly Positive Definite Functions 105
linearly dependent. I f K is n o t s t r i c t l y p o s i t i v e d e f i n i t e t h e n t h e r e e x i s t d i s t i n c t
p o i n t s Xi,..., XN a n d n o n z e r o coefficients c - such t h a t 3
N N
c c
Y2Y2 j kK{xj,x ) k = 0.
3=1 k=l
T h e first p a r t o f t h e p r o o f t h e r e f o r e i m p l i e s
3= 1
N o w we t a k e t h e H i l b e r t space i n n e r p r o d u c t w i t h a n a r b i t r a r y f u n c t i o n / 6 H a n d
use t h e r e p r o d u c i n g p r o p e r t y o f K t o o b t a i n
3= 1
N
c K x
= ^2 3(fi (-, 3))n
3= 1
N
3= 1
N
3=1
T h i s , however, i m p l i e s t h e l i n e a r dependence o f t h e p o i n t e v a l u a t i o n f u n c t i o n a l s
S Xj ( / ) f( j)i
x
j 1, ) AT, since t h e coefficients Cj were assumed t o be n o t a l l
zero. A n analogous a r g u m e n t c a n be used t o e s t a b l i s h t h e converse.
T h i s t h e o r e m p r o v i d e s one d i r e c t i o n o f t h e c o n n e c t i o n b e t w e e n s t r i c t l y p o s i t i v e
definite f u n c t i o n s a n d r e p r o d u c i n g kernels. H o w e v e r , we are also i n t e r e s t e d i n t h e
o t h e r d i r e c t i o n . Since t h e R B F s w e have b u i l t o u r i n t e r p o l a t i o n m e t h o d s f r o m are
strictly positive definite functions, we w a n t t o k n o w how t o construct a r e p r o d u c i n g
k e r n e l H i l b e r t space associated w i t h those s t r i c t l y p o s i t i v e d e f i n i t e basic f u n c t i o n s .
11/11* = (f,f)n = ( Y , ^ K { . , X j ) ^ c K { . , x ) )
k k u
3= 1 k=l
N N
S c
= Z2^C3 k{K(-,x ),K(-,x ))n j k
3= 1 k=l
N N
C C K X X
= ^2Yl i k ( 3' k)-
3=1 k=l
Therefore, we define t h e (possibly i n f i n i t e - d i m e n s i o n a l ) space
H (Cl)
K = s p a n { K ( - , y) : y Cl} (13.1)
w i t h a n associated b i l i n e a r f o r m (, )#- g i v e n b y
=
N K N K N K N K
f = ^2cjK(;Xj), XjeCl.
3= 1
Then
N K N K
S s
T h e o r e m 1 3 . 4 . Suppose $ G C ( R ) n L i ( R ) is a real-valued strictly positive def-
inite function. Define
S S 3
g = {fe L (R ) n C(R ):
2 -= e -MR )}
7
ana egwip /ns space with the bilinear form
i ,/ 9 , i / / M I R do;.
Then Q is a real Hilbert space with inner product (, -)g and reproducing kernel
s S
<&( ) . Hence, Q is the native space of & on R , i.e., Q = A / $ ( R ) and both
S
inner products coincide. In particular, every f G A / $ ( R ) can be recovered from its
s S
Fourier transform f G Za(R ) fl L ( R ) . 2
A n o t h e r c h a r a c t e r i z a t i o n o f t h e n a t i v e space is g i v e n i n t e r m s o f t h e eigenfunc-
t i o n s o f a linear o p e r a t o r associated w i t h t h e r e p r o d u c i n g k e r n e l . This operator,
T $ : L 2 ( P ) > L (fl),2 is g i v e n b y
/ $(x,y)(p (y)dy
k = \ <f) (x),
k k A; = 1 , 2 , . . . .
Jn
I n general, M e r c e r ' s t h e o r e m allows us t o c o n s t r u c t a r e p r o d u c i n g k e r n e l H i l b e r t
space 7i b y representing t h e f u n c t i o n s i n 7i as i n f i n i t e linear c o m b i n a t i o n s o f t h e
eigenfunctions, i.e.,
{
oo
/: / = ^ C f c 0 *
108 Meshfree Approximation Methods with M A T L A B
w h e r e we used t h e 7 i - o r t h o g o n a l i t y
(<t>j,4>k)H =
y/Xj^/Xk
of t h e eigenfunctions.
W e n o t e t h a t $ is indeed t h e r e p r o d u c i n g k e r n e l o f Ti since t h e eigenfunction
e x p a n s i o n (13.3) o f <E> a n d t h e o r t h o g o n a l i t y o f t h e eigenfunctions imply
oo oo
(/,$(, x)) n = (Y2cj<pj,^2Xk(f)k<Pk(x))n
j=i k=i
_ CkXk4>k{x)
oo
- ^c <f> ix)
k k = fix).
k=i
a n d t h e n a t i v e space i n n e r p r o d u c t c a n be w r i t t e n as
oo ^
(f,g)rt* = J^T-(/i0fc)L (n)(^,0fc>L (n), a 2 f,g e jV&iQ,).
A k
k=i
13.3 E x a m p l e s o f N a t i v e S p a c e s for P o p u l a r R a d i a l B a s i c
Functions
O n e also f r e q u e n t l y sees t h e d e f i n i t i o n
m a s
W 2 ( f l ) = { / G L ( f i ) n C(Q)
2 : Df G L {Q)
2 for a l l \a\ < m, a GN }, (13.5)
13. Reproducing Kernel Hilbert Spaces for Strictly Positive Definite Functions 109
S
w h i c h applies whenever Q, C R is a b o u n d e d d o m a i n . This interpretation will
make clear t h e c o n n e c t i o n b e t w e e n t h e n a t i v e s spaces o f Sobolev splines a n d t h o s e
3
of p o l y h a r m o n i c splines t o be discussed b e l o w . The n o r m of W ^ R ) is u s u a l l y
given by
I|/IIW7(R-) = { H j D a
/lli (R ) 2
S
||<m
A c c o r d i n g t o ( 1 3 . 4 ) , any s t r i c t l y p o s i t i v e d e f i n i t e f u n c t i o n 3> whose F o u r i e r t r a n s -
f o r m decays o n l y a l g e b r a i c a l l y has a Sobolev space as i t s n a t i v e space. I n p a r t i c u l a r ,
the M a t e r n functions
K -t(\\x\\)\\x\\f>-*
0 8
= ' ^ 2'
of Section 4.4 w i t h F o u r i e r t r a n s f o r m
* () = (i + H | ) -
/9
2 / 9
8
can i m m e d i a t e l y be seen t o have n a t i v e space JV* (R ) P = Wg(R) with B > s/2
( w h i c h is w h y some people refer t o t h e M a t e r n f u n c t i o n s as Sobolev splines).
W e n d l a n d ' s c o m p a c t l y s u p p o r t e d f u n c t i o n s <& fc = <>s,fc(|| l b ) o f C h a p t e r 11 c a n
S]
s 2 + k + 1 2 s
be s h o w n t o have n a t i v e spaces A / $ a k (R ) = W ^ 2 ^ (R ) (where the r e s t r i c t i o n
s > 3 is r e q u i r e d for t h e special case k = 0 ) .
N a t i v e spaces for s t r i c t l y c o n d i t i o n a l l y p o s i t i v e d e f i n i t e f u n c t i o n s c a n also be
c o n s t r u c t e d . H o w e v e r , since t h i s is m o r e t e c h n i c a l , we l i m i t e d t h e discussion a b o v e
t o s t r i c t l y p o s i t i v e d e f i n i t e f u n c t i o n s , a n d refer t h e i n t e r e s t e d reader t o t h e b o o k
[ W e n d l a n d (2005a)] o r t h e p a p e r s [Schaback (1999a); Schaback (2000a)]. W i t h t h e
extension of the theory t o s t r i c t l y c o n d i t i o n a l l y positive definite functions the native
spaces o f t h e r a d i a l powers a n d t h i n p l a t e (or surface) splines o f Sections 8.2 a n d
8.3 c a n be s h o w n t o be t h e so-called Beppo-Levi spaces o f o r d e r k
S S a S s
BL F C ( R ) = { / C{R ) : Df G L (R )
2 for a l l | a | = k, a e N } ,
a
where D denotes a generalized derivative o f o r d e r at (defined i n t h e same s p i r i t as
t h e generalized F o u r i e r t r a n s f o r m , see A p p e n d i x B ) . I n fact, t h e i n t e r s e c t i o n o f a l l
S 3
B e p p o - L e v i spaces B L f c ( R ) o f o r d e r k < m y i e l d s t h e Sobolev space W ^ R ) . In
S
t h e l i t e r a t u r e t h e B e p p o - L e v i spaces B L f c ( R ) are sometimes referred t o as homo-
s
geneous Sobolev spaces of order k. A l t e r n a t i v e l y , t h e B e p p o - L e v i spaces o n R are
defined as
S S S
BL F C ( R ) = { / G C{R ) : /(-)|| | | ? G L (R )},
2
a n d t h e f o r m u l a s g i v e n i n C h a p t e r 8 for t h e F o u r i e r t r a n s f o r m s o f r a d i a l powers a n d
t h i n p l a t e splines show i m m e d i a t e l y t h a t t h e i r n a t i v e spaces are B e p p o - L e v i spaces.
T h e s e m i - n o r m o n BLfc is g i v e n b y
s
Here t h e sine f u n c t i o n is defined for a n y x = (xi,... ,x )
s M as sine a; =
m o r e d e t a i l s o n S h a n n o n ' s s a m p l i n g t h e o r e m see, e.g., Chap-
ter 29 i n t h e b o o k [Cheney a n d L i g h t (1999)] or t h e p a p e r [Unser ( 2 0 0 0 ) ] .
Chapter 14
ill
112 Meshfree Approximation Methods with M A T L A B
order k i f
h) k
\\f-V$ \\ p = 0(h ) for/i-O.
h) k w
M o r e o v e r , i f we c a n also show t h a t | | / - V \\ f p / o(h ), then V has exact
L -approximation
p order k. W e w i l l c o n c e n t r a t e m o s t l y o n t h e case p = oo {i.e.,
p o i n t w i s e estimates), b u t a p p r o x i m a t i o n o r d e r i n o t h e r n o r m s c a n also be s t u d i e d .
I n order t o keep t h e f o l l o w i n g discussion as t r a n s p a r e n t as possible we w i l l r e s t r i c t
ourselves t o s t r i c t l y p o s i t i v e definite f u n c t i o n s . W i t h ( c o n s i d e r a b l y ) m o r e t e c h n i c a l
d e t a i l s t h e f o l l o w i n g can also be f o r m u l a t e d for s t r i c t l y c o n d i t i o n a l l y p o s i t i v e d e f i n i t e
f u n c t i o n s (see [ W e n d l a n d (2005a)] for d e t a i l s ) .
14.2 L a g r a n g e F o r m of t h e I n t e r p o l a n t a n d Cardinal
Basis Functions
Ac = y
T
with Aij = (xi - Xj), i,j = 1,...,N, c = [ci,..., c],
N and y =
T
[f(xi),...,f(xN)] has a u n i q u e s o l u t i o n . I n t h e f o l l o w i n g we w i l l consider the
m o r e general s i t u a t i o n w h e r e $ is a s t r i c t l y p o s i t i v e d e f i n i t e k e r n e l , i.e., t h e entries
of A are g i v e n b y A^ = &(xi,Xj). T h e uniqueness r e s u l t holds i n t h i s case also.
I n order t o o b t a i n t h e c a r d i n a l basis f u n c t i o n s u* , j = 1,... 3 ,N, w i t h the prop-
e r t y u*(xi) = Sij, i.e.,
Uj{Xi) =
we consider t h e linear s y s t e m
S
T h e o r e m 1 4 . 1 . Suppose <fr is a strictly positive definite kernel on I R . Then, for any
distinct points x\,..., XN, there exist functions u* e span{<E>(-, Xj),j = 1,..., N}
such that Uj(xi) = 5ij.
T h e r e f o r e , we can w r i t e t h e i n t e r p o l a n t Vf t o / a t x\,..., XN m t h e c a r d i n a l
form
14- The Power Function and Native Space Error Estimates 113
I t is o f interest t o n o t e t h a t t h e c a r d i n a l f u n c t i o n s d o n o t d e p e n d o n t h e d a t a
values o f t h e i n t e r p o l a t i o n p r o b l e m . O n c e t h e d a t a sites are fixed a n d t h e basic
f u n c t i o n is chosen w i t h a n a p p r o p r i a t e s h a p e p a r a m e t e r (whose o p t i m a l v a l u e w i l l
d e p e n d o n t h e d a t a ) , t h e n t h e c a r d i n a l f u n c t i o n s are d e t e r m i n e d b y t h e l i n e a r s y s t e m
(14.1). W e have p l o t t e d v a r i o u s c a r d i n a l f u n c t i o n s based o n t h e G a u s s i a n basic
f u n c t i o n w i t h shape p a r a m e t e r e = 5 i n F i g u r e s 1 4 . 1 - 1 4 . 3 . T h e dependence o n t h e
d a t a l o c a t i o n s is c l e a r l y a p p a r e n t w h e n c o m p a r i n g t h e different d a t a d i s t r i b u t i o n s
( u n i f o r m l y spaced i n F i g u r e 1 4 . 1 , t e n s o r - p r o d u c t C h e b y s h e v i n F i g u r e 14.2, and
H a l t o n p o i n t s i n F i g u r e 14.3). T h e d a t a sets c a n be seen i n F i g u r e 14.5 b e l o w .
Fig. 14.1 Cardinal functions for Gaussian interpolation (with e = 5) on 81 uniformly spaced
2
points in [0, l ] . Centered at an edge point (left) and at an interior point (right).
Fig. 14.2 Cardinal functions for Gaussian interpolation (with e = 5) on 81 tensor-product Cheby-
2
shev points in [0, l ] . Centered at an edge point (left) and at an interior point (right).
Fig. 14.3 Cardinal functions for Gaussian interpolation (with e = 5) on 81 Halton points in
2
[0, l ] . Centered at an edge point (left) and at an interior point (right).
Fig. 14.4 Cardinal functions for multiquadric interpolation (with e = 5) on 81 Halton points in
2
[0, l ] . Centered at an edge point (left) and at an interior point (right).
P r o g r a m 14.1. RBFCardinalFunction.m
7o R B F C a r d i n a l F u n c t i o n
/ Computes and p l o t s c a r d i n a l f u n c t i o n f o r 2D RBF i n t e r p o l a t i o n
/ C a l l s on: D i s t a n c e M a t r i x
0
1 r b f = @(e,r) e x p ( - ( e * r ) . " 2 ) ; ep = 5 ;
2 N = 81; g r i d t y p e = 'u';
3 neval = 8 0 ; M = neval~2;
14- The Power Function and Native Space Error Estimates 115
/ Load data p o i n t s
0
j=l i=l j = l
Then
s
D e f i n i t i o n 1 4 . 1 . Suppose fiCl a n d $ G C(fl x fl) is s t r i c t l y p o s i t i v e d e f i n i t e
S
on JR . For a n y d i s t i n c t p o i n t s X = {x,..., XN} C fl t h e power function is defined
by
w h e r e u* is t h e v e c t o r o f c a r d i n a l f u n c t i o n s f r o m T h e o r e m 1 4 . 1 .
116 Meshfree Approximation Methods with M A T L A B
U s i n g t h e d e f i n i t i o n o f t h e n a t i v e space n o r m f r o m t h e p r e v i o u s c h a p t e r we c a n
r e w r i t e t h e q u a d r a t i c f o r m Q(u) as
^2 Y2 Yl
N N N
N
2 x $ x
= ($(, x), $ ( - , W n } - )> ("' i))wsi)
N N
+ i j{(-i
u u x
i), H-> x
j))rt*(n)
i=i j=i
N N
3=1 3=1
2
N
(14.2)
3= 1
T T
Q(u) = (x, x) - 2u b(x) 4- u Au. (14.3)
T T
P* (x)
jX = y/Q(u*(x)) = \J${x,x) - 2(u*(x)) b(x) + (u*(x)) Au*(x).
x x x T
P*,x{x) = \/$( , )-(u*( )) b(x)
T
= \J$(x,x) - (u*(x)) Au*(x).
0<P*,*(aO < y / $ ( x , x ) .
P l o t s o f t h e p o w e r f u n c t i o n for t h e G a u s s i a n w i t h e 6 o n t h r e e different p o i n t
sets w i t h N = 81 i n t h e u n i t square are p r o v i d e d i n F i g u r e 14.5. T h e sets o f d a t a
p o i n t s are displayed o n t h e left, w h i l e t h e p l o t s o f t h e power f u n c t i o n are d i s p l a y e d
in the right column. Dependence o f t h e p o w e r f u n c t i o n o n t h e d a t a l o c a t i o n s is
clearly v i s i b l e . I n fact, t h i s c o n n e c t i o n was used i n a recent p a p e r [De M a r c h i et al.
(2005)] t o i t e r a t i v e l y o b t a i n a n o p t i m a l set o f d a t a l o c a t i o n s t h a t are independent
of t h e d a t a values.
A t t h i s p o i n t t h e power f u n c t i o n is m o s t l y a t h e o r e t i c a l t o o l t h a t helps us b e t t e r
u n d e r s t a n d e r r o r estimates since we c a n decouple t h e effects due t o t h e d a t a f u n c t i o n
/ f r o m those due t o t h e basic f u n c t i o n <& a n d t h e d a t a l o c a t i o n s X (see t h e f o l l o w i n g
T h e o r e m 14.2).
T h e power f u n c t i o n is defined i n a n analogous w a y for s t r i c t l y c o n d i t i o n a l l y
p o s i t i v e definite functions.
s s
T h e o r e m 14.2. Let C I , $ e C ( f 2 x f2) be strictly positive definite on R,
and suppose that the points X = {xi,... ,xjy} are distinct. Denote the interpolant
to f G A/$(f2) on X by Vf. Then for every x G O
f(x) = (/,$(", 2 0 ) ^ ( 0 ) .
W e express t h e i n t e r p o l a n t i n i t s c a r d i n a l f o r m a n d a p p l y t h e r e p r o d u c i n g p r o p e r t y
o f 3>. T h i s gives us
N
Vf(x) = J2f(*j)uj(x)
3= 1
N
3= 1
N
118 Meshfree Approximation Methods with M A T L A B
5
z-m-
0 0
1* o o
o o
0.8
0.6
o o
0.4
o o o o
0.2
o o o o o
o o o o o
0*
0.2 0.4 0.6 0.8 0 0
X
0 0
Fig. 14.5 Data sites and power function for Gaussian interpolant with e = 6 based on N = 81
uniformly distributed points (top), tensor-product Chebyshev points (middle), and Halton points
(bottom).
N o w a l l t h a t r e m a i n s t o be done is t o c o m b i n e t h e t w o f o r m u l a s j u s t d e r i v e d a n d
apply the Cauchy-Schwarz inequality. T h u s ,
N
\f(x)-V (x)\
f = (f,&(-,x)-Y2 j( )*(-i j))M-*(n)
u x x
3= 1
14- The Power Function and Native Space Error Estimates 119
<
i=i
= H/II.V(n)-P*,*(aO,
T h i s is analogous t o t h e s t a n d a r d e r r o r e s t i m a t e for p o l y n o m i a l i n t e r p o l a t i o n c i t e d
i n most n u m e r i c a l analysis t e x t s . N o t e , however, t h a t , for a n y g i v e n basic f u n c t i o n
$ , a change o f t h e shape p a r a m e t e r e w i l l have a n effect o n b o t h t e r m s i n t h e e r r o r
b o u n d i n T h e o r e m 14.2 since t h e n a t i v e space n o r m o f / varies w i t h e.
S
T h e o r e m 1 4 . 3 . Let fl C JR , and suppose 3> C(fl x fl) is strictly positive definite
on R . s
Let X = {x\,..., XN} be a set of distinct points in fl, and define the
quadratic form Q(u) as in (14-2). The minimum of Q(u) is given for the vector
u = u*(x) from Theorem 14-1, i-e.,
N
Q(u*(x)) < Q(u) for all u e R .
T T
Q{u) = $(x, x) - 2u b(x) + u Au.
T h e m i n i m u m o f t h i s q u a d r a t i c f o r m is g i v e n b y t h e s o l u t i o n o f t h e linear s y s t e m
Au = b(x).
S
D e f i n i t i o n 1 4 . 2 . A r e g i o n fl C I R satisfies a n interior cone condition i f t h e r e exists
a n angle 6 G ( 0 , 7 r / 2 ) a n d a r a d i u s r > 0 s u c h t h a t for every x fl t h e r e exists a
u n i t v e c t o r (x) such t h a t t h e cone
a T
C={x + Xy: y G E , | | y | | = 1 , y (x)
2 >cosd, AG[0,r]}
is c o n t a i n e d i n fl.
S
Theorem 1 4 . 4 . Suppose fl C I R is bounded and satisfies an interior cone condi-
tion, and let be a non-negative integer. Then there exist positive constants ho, c\,
and C2 such that for all X = {x,..., Ejv} C fl with h ,nx < ho and every x G fl
there exist numbers ui(x),..., ujv(x) with
N
(1) ^^Uj(x)p(xj) = p(x) for all polynomials p GII|,
i=i
N
(2) X>;(a:)| <ci,
3= 1
(3) Uj{x) = 0 if \\x - xj\\ 2 > c h Q.
2 Xj
9 l m
a n d t h e n o t a t i o n D^^iw, ) used b e l o w i n d i c a t e s t h a t t h e o p e r a t o r is a p p l i e d t o
$ ( i y , ) v i e w e d as a f u n c t i o n o f t h e second v a r i a b l e .
14- The Power Function and Native Space Error Estimates 121
T h e m u l t i v a r i a t e T a y l o r e x p a n s i o n o f t h e f u n c t i o n $ ( t o , ) centered a t w is g i v e n
by
$(w,z)= J2 3 i \ z - w f + R(w,z)
\/3\<2k P
'
w i t h remainder
R(w,z)= 2^ /3!
\0\=2k P
'
S
T h e o r e m 14.5. Suppose fl C JR is bounded and satisfies an interior cone condi-
2k
tion. Suppose <& G C (fl x fl) is symmetric and strictly positive definite. Denote
the interpolant to f G A/$(f2) on the set X byVf. Then there exist positive constants
ho and C (independent of x, f and such that
\f(x) - V ( )\
f X < Ch y/C^x)\\f\U^ ,
XtCi a)
provided h ,n
x < ho. Here
Proof. B y T h e o r e m 14.2 we k n o w
\f(x)-Vf(x)\<P*, (*)\\f\U.M- x
Therefore, we n o w derive t h e b o u n d
P*M*) < c h k
x ^ c * { x )
M o r e o v e r , we k n o w f r o m T h e o r e m 14.3 t h a t t h e q u a d r a t i c f o r m Q(u) is m i n i m i z e d
by u u*(x). Therefore, a n y o t h e r coefficient v e c t o r u w i l l y i e l d a n u p p e r b o u n d
on t h e power f u n c t i o n . W e t a k e u = u{x) f r o m T h e o r e m 14.4 so t h a t we are ensured
t o have p o l y n o m i a l p r e c i s i o n o f degree i > 2k 1.
For t h i s specific choice o f coefficients we have
[P^ {x)\
jX
2
< Q(u) = $(x, x) - 2 ^2 Uj$(x : xj) + Y2Y2 u u
i i( iix x
i)i
3 i 3
122 Meshfree Approximation Methods with M A T L A B
Q(u) = ${x,x)-2Y J m ( g
i - x
f R
+ (> x x
i)
\/3\<2k
,/3
D%$(xi,Xi)
W
+ EE * 3
E
\0\<2k
{xj Xi)^ ~\~ R{x>i, Xj)
+* E P
*(* ~ *^ + E E ( i 4 . 5 )
i |/3|<2fc ' i j
Now we c a n a p p l y t h e T a y l o r e x p a n s i o n a g a i n a n d m a k e t h e o b s e r v a t i o n that
D^(xi,Xj)
E
\/3\<2k
0!
(x - Xi) 13
= <E>(xi, x) - R(xi, x). (14.6)
Q(u) = Q(x, x) 23 u
2R(x,
j Xj) ^3 UiR(x{,Xj)
+ 23 U
i [( ii x x
) - R( i, x x
)\ (14.7)
O n e f i n a l T a y l o r e x p a n s i o n we need is ( u s i n g t h e s y m m e t r y o f <E)
$(xi,x) = &(x,Xi)= 2^ ^ ^ ~ ^ - ( x i - x ) p
+ R(x,Xi). (14.8)
\/3\<2k
Q( )u
= _
U
J
R(x, Xj) + R(xj,x) 23 UiR(x Xj)i:
2 k
T h e o r e m 14.5 says t h a t i n t e r p o l a t i o n w i t h a C s m o o t h k e r n e l <E> has approx-
i m a t i o n order k. T h u s , for i n f i n i t e l y s m o o t h s t r i c t l y p o s i t i v e definite functions
such as t h e Gaussians, L a g u e r r e - G a u s s i a n s , P o i s s o n r a d i a l f u n c t i o n s , a n d t h e gen-
eralized inverse m u l t i q u a d r i c s we see t h a t t h e a p p r o x i m a t i o n o r d e r k is a r b i t r a r i l y
h i g h . For s t r i c t l y p o s i t i v e d e f i n i t e f u n c t i o n s w i t h l i m i t e d smoothness such as t h e
M a t e r n f u n c t i o n s , t h e W h i t t a k e r r a d i a l f u n c t i o n s , as w e l l as a l l o f t h e compactly
s u p p o r t e d functions, t h e a p p r o x i m a t i o n o r d e r is l i m i t e d b y t h e smoothness o f t h e
basic f u n c t i o n .
T h e e s t i m a t e i n T h e o r e m 14.5 is s t i l l generic. I t does n o t f u l l y a c c o u n t for t h e
p a r t i c u l a r basic f u n c t i o n <& b e i n g used for t h e i n t e r p o l a t i o n since t h e factor C$(x)
s t i l l depends o n <E. M o r e o v e r , we p o i n t o u t t h a t t h e t e r m C$(x) may include a
h i d d e n dependence o n hx,n- For m o s t basic f u n c t i o n s i t w i l l be possible t o use
C<s>(x) t o "squeeze o u t " a d d i t i o n a l powers o f h. T h i s is t h e reason for s p l i t t i n g t h e
c o n s t a n t i n f r o n t o f t h e fo-power i n t o a generic C a n d a C$(x).
T h e s t a t e m e n t o f T h e o r e m 14.5 c a n be generalized for s t r i c t l y c o n d i t i o n a l l y pos-
i t i v e definite functions a n d also t o cover t h e error for a p p r o x i m a t i n g t h e d e r i v a t i v e s
of / b y derivatives o f Vf. W e state this general theorem w i t h o u t comment (c.f.
[ W e n d l a n d (2005a)] for d e t a i l s ) .
S
Theorem 14.6. Suppose Q C ] R is open and bounded and satisfies an interior
2k
cone condition. Suppose <& C (Q x Q.) is symmetric and strictly conditionally
positive definite of order m on W. Denote the interpolant to f j\f<&(Q) on the
(m 1 )-unisolvent set X by Vf. Fix ex NQ with \ot\ < k. Then there exist positive
constants ho and C (independent of x, f and such that
a a
\D f(x) - D V (x)\ f < C h ^ V C ^ l f U ^ ,
provided h ,n
x < ho- Here
15.1 N a t i v e S p a c e E r r o r B o u n d s for S p e c i f i c B a s i s F u n c t i o n s
\\f--P \\
f Loo(Q) <e*^|/Ufc ( n ) > (15.2)
125
126 Meshfree Approximation Methods with M A T L A B
for all data sites X with sufficiently small fill distance hx,a-
e
Moreover, if ip satisfies even \ipW(r)\ < M, then
e l n l
Wf-PfU^n) <e ^n' ||/|k. n) ( (15-3)
_ e 2 x 2 e T
Example 15.1. For Gaussians &(x) = e H H , e > 0 fixed, we have ip{r) = e~ * ',
e 2 i 2 r 2
so t h a t ipW(r) = (-l) e e- for t > = 0. T h u s , M = e , a n d t h e e r r o r b o u n d
0
(15.3) applies. T h i s k i n d o f e x p o n e n t i a l a p p r o x i m a t i o n o r d e r is u s u a l l y r e f e r r e d t o
as spectral (or even s u p e r - s p e c t r a l ) a p p r o x i m a t i o n o r d e r . W e emphasize t h a t t h i s
nice p r o p e r t y holds o n l y i n t h e n o n - s t a t i o n a r y s e t t i n g a n d for d a t a f u n c t i o n s / t h a t
are i n t h e n a t i v e space of t h e Gaussians such as b a n d - l i m i t e d f u n c t i o n s .
2
Example 15.2. For generalized (inverse) m u l t i q u a d r i c s $(x) = ( 1 + |Ja5|| )^, (3 < 0,
e 6
or 0 < 3 N , we have ip(r) = ( 1 + r ) ^ . I n t h i s case one c a n show t h a t \ip (r)\ < \M
w h e n e v e r ^ > \3~\. H e r e M = 1 + 1/3+11. T h e r e f o r e , t h e e r r o r e s t i m a t e (15.2) applies,
i.e., i n t h e n o n - s t a t i o n a r y s e t t i n g generalized (inverse) m u l t i q u a d r i c s have s p e c t r a l
a p p r o x i m a t i o n order.
/2 2 2 2
E x a m p l e 1 5 . 3 . For L a g u e r r e - G a u s s i a n s $(cc) = L (\\ex\\ )e- ^ , n e > 0 fixed,
2 2 e2r e
w e have ip(r) = Ln (e r)e~ a n d t h e d e r i v a t i v e s ip^ w i l l be b o u n d e d b y ip^ \0) =
2i
Pn{) , where p n is a p o l y n o m i a l o f degree n . T h u s , t h e a p p r o x i m a t i o n p o w e r o f
L a g u e r r e - G a u s s i a n s falls b e t w e e n (15.3) a n d (15.2) a n d L a g u e r r e - G a u s s i a n s have at
least s p e c t r a l a p p r o x i m a t i o n p o w e r .
i n d e p e n d e n t o f x (see [ W e n d l a n d ( 2 0 0 5 a ) ] ) . T h e r e f o r e , t h i s results i n t h e f o l l o w i n g
e r r o r e s t i m a t e s (see, e.g., [ W e n d l a n d ( 2 0 0 5 a ) ] , o r t h e m u c h earlier [ W u a n d Schaback
(1993)] w h e r e o t h e r p r o o f t e c h n i q u e s were used).
K0 > w e e t
E x a m p l e 1 5 . 4 . For t h e M a t e r n f u n c t i o n s <I>(CE) = ~^-^vlt) ' @ '
a
\D~f{x)-D V,{x)\ < Ch -*-\f\^ . x m (15.4)
R a d i a l powers a n d t h i n p l a t e splines c a n be i n t e r p r e t e d as a g e n e r a l i z a t i o n o f u n i -
v a r i a t e n a t u r a l splines. T h e r e f o r e , w e k n o w t h a t t h e a p p r o x i m a t i o n o r d e r e s t i m a t e s
o b t a i n e d v i a t h e n a t i v e space a p p r o a c h are n o t o p t i m a l . For e x a m p l e , for i n t e r -
p o l a t i o n w i t h u n i v a r i a t e piecewise l i n e a r splines ( i . e . , $(cc) = ||cc|| i n x G R ) w e
k n o w t h e a p p r o x i m a t i o n o r d e r t o b e O(h), whereas t h e e s t i m a t e (15.6) y i e l d s o n l y
2 2
a p p r o x i m a t i o n o r d e r 0{h}/ ). S i m i l a r l y , for t h i n p l a t e splines &(x) = ||cc|| l o g ||ic||
2
one w o u l d expect o r d e r 0(h ) i n t h e case o f p u r e f u n c t i o n a p p r o x i m a t i o n . H o w e v e r ,
t h e e s t i m a t e (15.7) yields o n l y 0(h). These t w o e x a m p l e s suggest t h a t i t s h o u l d be
possible t o " d o u b l e " t h e a p p r o x i m a t i o n o r d e r s o b t a i n e d t h u s far.
O n e c a n i m p r o v e t h e e s t i m a t e s for f u n c t i o n s w i t h finite smoothness ( i . e . , M a t e r n
f u n c t i o n s , W e n d l a n d f u n c t i o n s , r a d i a l powers, a n d t h i n p l a t e splines) b y e i t h e r (or
b o t h ) o f t h e f o l l o w i n g t w o ideas:
b y r e q u i r i n g t h e d a t a f u n c t i o n / t o be even s m o o t h e r t h a n w h a t t h e n a t i v e
space prescribes, i.e., b y b u i l d i n g c e r t a i n b o u n d a r y c o n d i t i o n s i n t o t h e n a t i v e
space;
b y u s i n g weaker n o r m s t o m e a s u r e t h e e r r o r .
(1997)])
2k 1+s
11/ - P / I U n ) < Ch + \\f\\ , r s
a ( w k+ + {R3)1 (15.8)
128 Meshfree Approximation Methods with M A T L A B
2 f c + 1 + S S s
w h e r e / is assumed t o lie i n t h e subspace W 2 ( R ) o f Af<j,(R ) = W^"^ . For
r =
e x a m p l e , for t h e p o p u l a r basic f u n c t i o n ^3,1 ( ) 0- ~~ r)+(4r + 1) w e have
6
\\f~V \\ ^ f L <Ch \\f\\ e . w iRs)
N o t e t h a t t h e n u m e r i c a l e x p e r i m e n t s i n T a b l e 12.2 p r o d u c e d R M S - c o n v e r g e n c e r a t e s
o n l y as h i g h as 4.5.
For r a d i a l powers a n d t h i n p l a t e splines o n e o b t a i n s L 2 - e r r o r e s t i m a t e s o f o r -
f3+s 2 k + s
der 0(h ) and 0 ( h ), respectively. T h e s e e s t i m a t e s a r e o p t i m a l , i.e., e x a c t
a p p r o x i m a t i o n orders, as s h o w n i n [ B e j a n c u ( 1 9 9 9 ) ] .
M o r e w o r k o n i m p r o v e d e r r o r b o u n d s c a n b e f o u n d i n , e.g., [ J o h n s o n (2004a)]
or [Schaback ( 1 9 9 9 b ) ] .
T h e e r r o r b o u n d s m e n t i o n e d so f a r were a l l v a l i d u n d e r t h e a s s u m p t i o n t h a t t h e
f u n c t i o n / p r o v i d i n g t h e d a t a c a m e f r o m ( a subspace o f ) t h e n a t i v e space o f t h e
RBF employed i n the interpolation. W e n o w m e n t i o n a f e w recent r e s u l t s that
provide error bounds for i n t e r p o l a t i o n o f functions / not i n t h e n a t i v e space o f
t h e basic f u n c t i o n . I n p a r t i c u l a r , t h e case w h e n / lies i n some S o b o l e v space is o f
great i n t e r e s t . A r a t h e r general t h e o r e m w a s r e c e n t l y g i v e n i n [ N a r c o w i c h et al.
(2005)]. I n t h i s t h e o r e m N a r c o w i c h , W a r d a n d W e n d l a n d p r o v i d e S o b o l e v b o u n d s
for f u n c t i o n s w i t h m a n y zeros. H o w e v e r , since t h e i n t e r p o l a t i o n e r r o r f u n c t i o n is
j u s t such a f u n c t i o n , these b o u n d s have a d i r e c t a p p l i c a t i o n t o o u r s i t u a t i o n . W e
p o i n t o u t t h a t t h i s t h e o r e m a g a i n applies t o t h e n o n - s t a t i o n a r y s e t t i n g .
T h e o r e m 15.2. Let k be a positive integer, 0 < a < l , l < p < 00, 1 < q < 0 0
and let a be a multi-index satisfying k > \ct\ + s/p or, for p = 1, k > \a\ -f- s.
Let X C fl be a discrete set with fill distance h hx,n where fl is a compact set
k+a
with Lipschitz boundary which satisfies an interior cone condition. If u e W (fl)
satisfies u\x = 0, then
k+a a s 1 1
1 , < rh ~\ \- ( /P- /l)+1 | 7/ t ,
k + a
Suppose we have a n i n t e r p o l a t i o n process V : W ' (fl) > V t h a t m a p s S o b o l e v
k + a
f u n c t i o n s t o a f i n i t e - d i m e n s i o n a l subspace V o f W (fl) w i t h t h e a d d i t i o n a l p r o p -
f e + C T
e r t y \pf\ k+a^ w < |/lvK ( f i ) ' t h e n T h e o r e m 15.2 i m m e d i a t e l y y i e l d s t h e e r r o r
estimate
f c + C T s
T h e a d d i t i o n a l p r o p e r t y [Pf\ k+a^ w < |/lw ( n ) * c e r t a i n l y satisfied p r o v i d e d
t h e n a t i v e space o f t h e basic f u n c t i o n is a Sobolev space. T h u s , T h e o r e m 15.2
p r o v i d e s a n a l t e r n a t i v e t o t h e p o w e r f u n c t i o n a p p r o a c h discussed i n t h e p r e v i o u s
15. Refined and Improved Error Bounds 129
ci(l + ||u>||l)- T
< $(w) < c (l + 2 ||a;||l)- , T
- oo, u e R,
s
(15.9)
t h e n t h e above e r r o r e s t i m a t e h o l d s w i t h r = k + a a n d p = 2 p r o v i d e d t h e f i l l
distance is sufficiently s m a l l . E x a m p l e s o f basic f u n c t i o n s w i t h a n a p p r o p r i a t e l y
decaying F o u r i e r t r a n s f o r m are p r o v i d e d b y t h e families o f W e n d l a n d or M a t e r n
functions. I n a d d i t i o n , N a r c o w i c h , W a r d a n d W e n d l a n d s h o w t h a t analogous e r r o r
b o u n d s h o l d for r a d i a l powers a n d t h i n p l a t e splines (whose n a t i v e spaces are B e p p o -
L e v i spaces).
For f u n c t i o n s / outside t h e n a t i v e space o f a basic f u n c t i o n <E> whose F o u r i e r
t r a n s f o r m satisfies (15.9) N a r c o w i c h , W a r d a n d W e n d l a n d p r o v e
T h e o r e m 15.3. Let k and n be integers with 0 < n < k < r and k > s/2, and let
k
f 6 C (fl). Also suppose that X { x i , . . . ,x^} C fl satisfies diam(A') < 1 with
sufficiently small fill distance. Then for any 1 < q < oo we have
We r e m i n d t h e reader t h a t t h e f i l l distance c o r r e s p o n d s t o t h e r a d i u s o f t h e
largest possible e m p t y b a l l t h a t c a n be p l a c e d b e t w e e n t h e p o i n t s i n X. T h e sep-
a r a t i o n distance (c.f. C h a p t e r 16), o n t h e o t h e r h a n d , c a n be i n t e r p r e t e d as t h e
r a d i u s o f t h e largest b a l l t h a t c a n be p l a c e d a r o u n d every p o i n t i n X such t h a t n o
t w o balls o v e r l a p . T h u s , t h e m e s h r a t i o is a measure o f t h e n o n - u n i f o r m i t y o f t h e
d i s t r i b u t i o n o f the points i n X.
S i m i l a r results were o b t a i n e d earlier i n [ B r o w n l e e a n d L i g h t (2004)] (for r a d i a l
powers a n d t h i n p l a t e splines o n l y ) , a n d i n [ Y o o n (2003)] (for shifted surface splines,
see b e l o w ) .
s 2
\f-r \ f L o o <ch^- / \\f\\ c 0 { m
L o w e r b o u n d s o n t h e a p p r o x i m a t i o n o r d e r for a p p r o x i m a t i o n b y p o l y h a r m o n i c
splines were r e c e n t l y p r o v i d e d i n [ M a i o r o v (2005)]. M a i o r o v studies for a n y 1 <
p, q < oo a n d 3/s > (1/p l/q)+ the error E of //^-approximation of functions
b y p o l y h a r m o n i c splines. M o r e precisely,
a s
E(W^{0, l] ),n {p ,3),L {^,lY))>cN-^ ,
N p q
3= 1
where
Ve
h = <p(h-)-
s t u d i e d b y [ B u h m a n n (1989a)] o n i n f i n i t e l a t t i c e s . H o w e v e r , for q u a s i - i n t e r p o l a t i o n
the approximate approximation a p p r o a c h o f M a z ' y a shows t h a t i t is possible t o
choose eo i n such a w a y t h a t t h e level at w h i c h t h e s a t u r a t i o n occurs can be c o n -
t r o l l e d (see, e.g., [ M a z ' y a a n d S c h m i d t ( 1 9 9 6 ) ] ) . T h e r e f o r e , Gaussians m a y v e r y
w e l l be used for s t a t i o n a r y i n t e r p o l a t i o n p r o v i d e d a n a p p r o p r i a t e i n i t i a l shape pa-
rameter is chosen. W e w i l l i l l u s t r a t e t h i s b e h a v i o r i n t h e n e x t c h a p t e r . T h e same
k i n d o f a r g u m e n t also applies t o t h e L a g u e r r e - G a u s s i a n s o f S e c t i o n 4.2.
P<b, ,x{<x)
h < C (e h s^J
h x
E x a m p l e 1 5 . 1 1 . S t a t i o n a r y i n t e r p o l a t i o n w i t h (inverse) m u l t i q u a d r i c s , r a d i a l p o w -
ers a n d t h i n p l a t e splines presents no difficulties. I n fact, [Schaback (1995c)] shows
t h a t t h e n a t i v e space e r r o r b o u n d for t h i n p l a t e splines a n d r a d i a l powers is i n v a r i a n t
under a s t a t i o n a r y scaling. T h e r e f o r e , t h e n o n - s t a t i o n a r y b o u n d o f T h e o r e m 15.3
applies i n t h e s t a t i o n a r y case also. T h e a d v a n t a g e o f scaling t h i n p l a t e splines or
r a d i a l powers comes f r o m t h e a d d e d s t a b i l i t y one can g a i n b y p r e v e n t i n g t h e sepa-
r a t i o n distance f r o m b e c o m i n g t o o s m a l l (see C h a p t e r 16 a n d t h e w o r k o f Iske o n
local p o l y h a r m o n i c spline a p p r o x i m a t i o n , e.g., [Iske ( 2 0 0 4 ) ] ) .
Y o o n p r o v i d e s e r r o r e s t i m a t e s for s t a t i o n a r y a p p r o x i m a t i o n o f r o u g h f u n c t i o n s
(i.e., f u n c t i o n s t h a t are n o t i n t h e n a t i v e space o f t h e basic f u n c t i o n ) b y so-called
shifted surface splines. S h i f t e d surface splines are o f t h e f o r m
I r / 3 S / 2 1 2 7 2
*M = J (- ) " (1 + I N I ) " - , * odd,
1 } 2 s 2 2 1/2
\ ( - l y - ^ + ^ l + \\x\ f~ l log(l + ||a:|| ) > s even,
where s/2 < 6 G N . These f u n c t i o n s i n c l u d e a l l o f t h e (inverse) m u l t i q u a d r i c s ,
r a d i a l powers a n d t h i n p l a t e splines.
Y o o n has t h e f o l l o w i n g t h e o r e m (see [ Y o o n (2003)] for t h e L p case, a n d [ Y o o n
(2001)] for Loo b o u n d s o n l y ) .
have
11/ - V \ \
f L p i a ) < Ch^\f\ , w { R 3 ) , 1 < p < oo,
with
7 p = mm{3,3-s/2 + s/p}.
7 P _ / 3 + Q
I I / - ^ / I I M ) = ^ ) -
15.5 C o n v e r g e n c e w i t h R e s p e c t to the S h a p e P a r a m e t e r
15.6 P o l y n o m i a l I n t e r p o l a t i o n as t h e L i m i t of R B F I n t e r p o l a t i o n
f( )
V x
= Y, M\\ ( - 3)\\)>
C X X x e
W,b]cM,
t o f u n c t i o n values a t N d i s t i n c t d a t a sites t e n d s t o t h e L a g r a n g e i n t e r p o l a t i n g
p o l y n o m i a l o f / as e 0.
T h e m u l t i v a r i a t e case is m o r e c o m p l i c a t e d . However, the l i m i t i n g R B F inter-
p o l a n t ( p r o v i d e d i t exists) is g i v e n b y a low-degree m u l t i v a r i a t e p o l y n o m i a l (see
[Larsson a n d F o r n b e r g ( 2 0 0 5 ) ; Schaback ( 2 0 0 5 ) ; Schaback ( 2 0 0 6 b ) ] ) . F o r e x a m p l e ,
if t h e d a t a sites are l o c a t e d s u c h t h a t t h e y g u a r a n t e e a u n i q u e p o l y n o m i a l i n t e r -
p o l a n t , t h e n t h e l i m i t i n g R B F i n t e r p o l a n t is g i v e n b y t h i s p o l y n o m i a l . I f p o l y n o -
m i a l i n t e r p o l a t i o n is n o t u n i q u e , t h e n t h e R B F l i m i t depends o n t h e choice o f basic
function. H o w e v e r , these s t a t e m e n t s r e q u i r e t h e R B F s t o satisfy a n (unproven)
c o n d i t i o n o n c e r t a i n coefficient m a t r i c e s A j. Pt I n [ L a r s s o n a n d F o r n b e r g (2005)]
t h e a u t h o r s also p r o v i d e a n e x p l a n a t i o n for t h e e r r o r b e h a v i o r for s m a l l values o f
t h e shape p a r a m e t e r , a n d for t h e existence o f a n o p t i m a l ( p o s i t i v e ) value o f e g i v -
ing rise t o a g l o b a l m i n i m u m o f t h e e r r o r f u n c t i o n . F o r t h e special case o f scaled
Gaussians Schaback [Schaback (2005)] shows t h a t t h e R B F i n t e r p o l a n t converges
t o t h e d e B o o r a n d R o n least p o l y n o m i a l i n t e r p o l a n t (see [ d e B o o r a n d R o n ( 1 9 9 0 ) ;
d e B o o r a n d R o n (1992a)] a n d also [ d e B o o r ( 2 0 0 6 ) ] ) as e - > 0.
I n [ F o r n b e r g a n d W r i g h t (2004)] t h e a u t h o r s describe a so-called Contour-Pade
a l g o r i t h m t h a t makes i t possible (for d a t a sets o f r e l a t i v e l y m o d e s t size) t o c o m p u t e
t h e R B F i n t e r p o l a n t for a l l values o f t h e shape p a r a m e t e r e i n c l u d i n g t h e l i m i t i n g
case e -> 0. W e present some n u m e r i c a l r e s u l t o b t a i n e d w i t h G r a d y W r i g h t ' s
M A T L A B toolbox in Chapter 17.
We w i l l later exploit the connection between R B F and p o l y n o m i a l interpolants
t o design n u m e r i c a l solvers for p a r t i a l d i f f e r e n t i a l e q u a t i o n s .
Chapter 16
16.1 S t a b i l i t y a n d C o n d i t i o n i n g of R a d i a l B a s i s Function
Interpolants
A s t a n d a r d c r i t e r i o n for m e a s u r i n g t h e n u m e r i c a l s t a b i l i t y o f a n approximation
m e t h o d is i t s c o n d i t i o n n u m b e r . I n p a r t i c u l a r , for r a d i a l basis f u n c t i o n i n t e r p o l a t i o n
we need t o l o o k at t h e c o n d i t i o n n u m b e r o f t h e i n t e r p o l a t i o n m a t r i x A w i t h entries
Aij = $(xi Xj). For a n y m a t r i x A i t s -condition
2 number is g i v e n b y
1
cond(A) = | | A | | | | A - | |
2 2 = ^ ,
^min
where c r m a x a n d c r i n are t h e largest a n d smallest s i n g u l a r values o f A.
m I f we
c o n c e n t r a t e o n p o s i t i v e d e f i n i t e m a t r i c e s A, t h e n t h e c o n d i t i o n n u m b e r o f A can
also t a k e be c o m p u t e d as t h e r a t i o
Amax
Amin
of t h e largest a n d smallest eigenvalues.
W h a t do we k n o w a b o u t these eigenvalues? F i r s t , G e r s h g o r i n ' s t h e o r e m (see,
e . g . , [Meyer (2000)]) says t h a t
JV
| Amax An | 5: ^ ^ | Aij |
J= I
A m a x < JV$(0)
b y t h e p r o p e r t i e s o f p o s i t i v e d e f i n i t e f u n c t i o n s ( P r o p e r t y (4) i n T h e o r e m 3.1). N o w ,
s
as l o n g as t h e d a t a are n o t t o o w i l d l y d i s t r i b u t e d , N w i l l g r o w as h x Q which
is acceptable. Therefore, t h e m a i n w o r k i n e s t a b l i s h i n g a b o u n d for t h e c o n d i t i o n
n u m b e r o f A lies i n finding lower b o u n d s for A i m n (or c o r r e s p o n d i n g l y u p p e r b o u n d s
135
136 Meshfree Approximation Methods with M A T L A B
_ 1
for t h e n o r m o f t h e inverse | | A | | 2 ) . T h i s is t h e s u b j e c t o f several papers b y B a l l ,
N a r c o w i c h , S i v a k u m a r a n d W a r d [ B a l l et al. (1992); N a r c o w i c h et al. (1994);
N a r c o w i c h a n d W a r d (1991a); N a r c o w i c h a n d W a r d ( 1 9 9 1 b ) ; N a r c o w i c h a n d W a r d
(1992)] w h o m a k e use o f a r e s u l t f r o m [ B a l l (1992)] o n eigenvalues o f distance
m a t r i c e s . B a l l ' s result follows f r o m t h e Rayleigh quotient (or t h e C o u r a n t - F i s c h e r
T h e o r e m 9.5), w h i c h gives t h e smallest eigenvalue o f a s y m m e t r i c p o s i t i v e d e f i n i t e
m a t r i x as
T
c Ac
Amin = mm .
w 1
cK \0 C C
s s
L e m m a 1 6 . 1 . Let X\,..., x^ be distinct points in R and let <3> : R > R be either
strictly positive definite, or strictly conditionally negative definite of order one with
<&(0) < 0. Also, let A be the interpolation matrix with entries A^ = $>(xi Xj). If
the inequality
N N
2
CiCj-Ay > 9\\c\\ 2
i=l j=l
c =
is satisfied whenever the components of c satisfy X^jLi j > then
l
\\A- \\ <6-\ 2
N o t e t h a t for p o s i t i v e d e f i n i t e m a t r i c e s t h e R a y l e i g h q u o t i e n t i m p l i e s 6 = A i n m
qx = \ m i n |jac< - Xj\\ . 2
We c a n p i c t u r e qx as t h e r a d i u s o f t h e largest b a l l t h a t c a n be p l a c e d a r o u n d e v e r y
p o i n t i n X such t h a t no t w o balls o v e r l a p (see F i g u r e 16.1). T h e s e p a r a t i o n d i s t a n c e
is sometimes also referred t o as t h e packing radius. I n our MATLAB code we c a n
c o m p u t e t h e s e p a r a t i o n distance v i a
qX = m i n ( m i n ( D M _ d a t a + e y e ( s i z e ( D M _ d a t a ) ) ) ) / 2
16. Stability and Trade-Off Principles 137
Fig. 16.1 The separation distance for N = 25 Halton points (qx ~ 0.0597).
e 2 2
Example 16.1. For Gaussians &(x) = e~ ^ one o b t a i n s
s 4 0 7 1 s 2
A m i n > C ( v ^ ) - e -
s - / ^ V -
W e see t h a t , for a fixed shape p a r a m e t e r e, t h e lower b o u n d for A i m n goes e x p o -
n e n t i a l l y t o zero as t h e s e p a r a t i o n d i s t a n c e qx decreases. Since we observed above
t h a t t h e c o n d i t i o n n u m b e r o f t h e i n t e r p o l a t i o n m a t r i x A depends o n t h e r a t i o o f
its largest a n d smallest eigenvalues a n d t h e g r o w t h o f A m a x is o f o r d e r N we see
t h a t t h e c o n d i t i o n n u m b e r g r o w s e x p o n e n t i a l l y w i t h decreasing s e p a r a t i o n d i s t a n c e .
T h i s shows t h a t , i f one adds m o r e i n t e r p o l a t i o n p o i n t s i n o r d e r t o i m p r o v e t h e ac-
c u r a c y o f t h e i n t e r p o l a n t ( w i t h i n t h e same d o m a i n fl), t h e n t h e p r o b l e m becomes
138 Meshfree Approximation Methods with M A T L A B
i + i 2 M e
A m i n > C(a, A e ) 4 " e- '/<* )
/ 3 + 1 2 / 3
E x a m p l e 1 6 . 3 . F o r t h i n p l a t e splines $(cc) = ( - l ) ||cc|| l o g \\x\\, 0 G N , one
obtains
s 2
A m i n > C s C / 3 (2M )- - ^^s
2 p
E x a m p l e 1 6 . 4 . F o r t h e r a d i a l powers $ ( J C ) = ( - l ) ^ / ! \\x\\ , 0 < Q 2 N , one
obtains
s / 3
A m i n > C c s / 3 (2M )- - s 4
w i t h a n o t h e r e x p l i c i t l y k n o w n c o n s t a n t c p (different f r o m c p i n E x a m p l e 3 ) . A g a i n ,
t h e decay is o f p o l y n o m i a l order.
E x a m p l e 1 6 . 5 . F o r t h e c o m p a c t l y s u p p o r t e d f u n c t i o n s o f W e n d l a n d $> k(x) Sj
(p k(||x||)
St one o b t a i n s
k+1
A m i n > C(s,k)q%
17.1 I n t e r p o l a t i o n for e 0
T h e e r r o r b o u n d s we r e v i e w e d i n p r e v i o u s c h a p t e r s give us some i n s i g h t i n t o
t h e first issue. I f we k n o w t h a t t h e d a t a come f r o m a v e r y s m o o t h f u n c t i o n , t h e n
a p p l i c a t i o n o f one o f t h e s m o o t h e r basic f u n c t i o n s is called for. O t h e r w i s e , t h e r e
is n o t m u c h t o be g a i n e d f r o m d o i n g so. I n fact, these f u n c t i o n s m a y a d d t o o
m u c h smoothness t o t h e i n t e r p o l a n t . A first a t t e m p t at p r o v i d i n g guidelines for t h e
selection o f a p p r o p r i a t e basic f u n c t i o n s (or kernels) c a n be f o u n d i n [Schaback a n d
W e n d l a n d (2006)]. W e w i l l n o t p u r s u e t h i s issue a n y f u r t h e r .
I n s t e a d we w a n t t o focus o u r a t t e n t i o n o n t h e second q u e s t i o n , i.e., t h e choice o f
t h e shape p a r a m e t e r e. A n u m b e r o f strategies c a n be used t o g u i d e us i n m a k i n g a
decision. W e w i l l assume t h r o u g h o u t t h a t a (fixed) basic f u n c t i o n has been chosen,
a n d t h a t we w i l l use o n l y one value t o scale a l l basis f u n c t i o n s u n i f o r m l y . Clearly,
one can also follow o t h e r strategies such as u s i n g a shape p a r a m e t e r t h a t varies
w i t h j , or even basic f u n c t i o n s t h a t v a r y w i t h j . W h i l e some w o r k has b e e n d o n e
i n these d i r e c t i o n s (see, e.g., [ B o z z i n i et al. ( 2 0 0 2 ) ; K a n s a a n d C a r l s o n (1992);
Schaback a n d W e n d l a n d ( 2 0 0 0 b ) ; F o r n b e r g a n d Z u e v ( 2 0 0 6 ) ] ) , n o t m u c h c o n c r e t e
c a n be said i n these cases.
141
142 Meshfree Approximation Methods with M A T L A B
T h e s i m p l e s t s t r a t e g y is t o p e r f o r m a series o f i n t e r p o l a t i o n e x p e r i m e n t s w i t h v a r y -
i n g shape p a r a m e t e r , a n d t h e n t o p i c k t h e "best" one. T h i s s t r a t e g y c a n be used i f
we k n o w t h e f u n c t i o n / t h a t g e n e r a t e d t h e d a t a , a n d t h e r e f o r e c a n c a l c u l a t e some
s o r t o f e r r o r for t h e i n t e r p o l a n t . O f course, i f we a l r e a d y k n o w / , t h e n t h e exercise
o f f i n d i n g a n i n t e r p o l a n t Vf m a y be m o s t l y p o i n t l e s s . H o w e v e r , t h i s is t h e s t r a t e g y
we used for t h e "academic" examples i n C h a p t e r 2.
I f we do n o t have a n y k n o w l e d g e o f / , t h e n i t becomes v e r y d i f f i c u l t t o decide
w h a t "best" means. O n e ( n o n - o p t i m a l ) c r i t e r i o n w e used i n C h a p t e r 2 was based
o n t h e t r a d e - o f f p r i n c i p l e , i.e., t h e fact t h a t for s m a l l e t h e e r r o r i m p r o v e s w h i l e
t h e c o n d i t i o n n u m b e r g r o w s . W e t h e n defined "best" t o be t h e smallest e for w h i c h
M A T L A B d i d n o t issue a n e a r - s i n g u l a r w a r n i n g .
I n m a n y cases selection o f a n o p t i m a l shape p a r a m e t e r v i a trial and error will
e n d u p b e i n g a r a t h e r s u b j e c t i v e process. H o w e v e r , t h i s m a y p r e s e n t l y be the
approach taken b y most practitioners.
For c o m p a r i s o n w i t h t h e o t h e r s e l e c t i o n m e t h o d s f e a t u r e d b e l o w we present t h r e e
o n e - d i m e n s i o n a l test cases for w h i c h w e k n o w t h e d a t a f u n c t i o n / . W e use
2 2 2 2
JP ( )2 X = | ^ -(9x-2) /4
e + e -(9x+l) /49^ + ^ ~(Qx-7) /4
e _ _Lg-(9x-4) ^
F (x)
3 = ( l - \x - \ \ j (\ + 5\x - \ \ - 27\x - i | 2 ) .
A n o t h e r s t r a t e g y is suggested b y t h e e r r o r a n a l y s i s o f C h a p t e r 14. W e s h o w e d t h e r e
i n T h e o r e m 14.2 t h a t
\f(x)-Vf{x)\<P* (x)\\f\\ ,
tX Wn)
17. Numerical Evidence for Approximation Order Results 143
Table 17.1 Optimal e values based on Gaussian interpolation to various test func-
tions in I D for various choices of N uniform points.
Fi F 2 F 3
Fig. 17.1 Optimal e curves based on Gaussian interpolation in I D for various choices of N uniform
2
points. Data sampled from sine function Fx (left) and C oscillatory function F3 (right).
1
P*,x{x) = ^ ( x , ^ ) - ^ ) ) ^ - ^ ) ,
T
where A is t h e i n t e r p o l a t i o n m a t r i x a n d 6 = [ $ ( - , x i ) , . . . , $ ( - , C j v ) ] . T h i s f o r m u l a
is i m p l e m e n t e d o n lines 1 5 - 1 8 i n t h e M A T L A B p r o g r a m P o w e r f u n c t i o n 2 D .m. W e
c o m p u t e t h e inverse o f A u s i n g t h e f u n c t i o n p i n v w h i c h is based o n t h e s i n g u l a r
value d e c o m p o s i t i o n o f A a n d t h e r e f o r e g u a r a n t e e s m a x i m u m s t a b i l i t y . A l s o , due t o
r o u n d o f f some o f t h e a r g u m e n t s o f t h e s q r t f u n c t i o n o n l i n e 18 come o u t n e g a t i v e .
144 Meshfree Approximation Methods with M A T L A B
This explains the use of the r e a l command. The vectors b(x) are just rows of the
evaluation matrix if x is taken from the grid of evaluation points we used earlier
for error computations and plotting purposes. Except for the loop over the shape
parameter e (lines 12-20) the rest of the program is similar to earlier code.
P r o g r a m 1 7 . 1 . Powerfunction2D.m
% Powerfunction2D
S c r i p t t h a t f i n d s "optimal" shape parameter by computing t h e power
'/, f u n c t i o n f o r the 2D RBF i n t e r p o l a t i o n approach w i t h v a r y i n g e p s i l o n
/, C a l l s on: DistanceMatrix
1 r b f = @(e,r) e x p ( - ( e * r ) . ~ 2 ) ; 7. Define t h e Gaussian RBF
% Parameters f o r shape parameter loop below
2 mine = 0; maxe = 20;
3 ne = 500; ep = linspace(mine,maxe,ne);
'/, Number and type of data p o i n t s
4 N = 81; g r i d t y p e = 'u';
7, R e s o l u t i o n of g r i d f o r power f u n c t i o n norm computation
5 n e v a l = 2 0 ; M = neval"2;
% Load data p o i n t s
, ,
6 name = s p r i n t f ( Data2D_7od7oS ,N,gridtype) ; load(name)
7 ctrs = dsites; % c e n t e r s c o i n c i d e w i t h data s i t e s
8 g r i d = l i n s p a c e ( 0 , 1 , n e v a l ) ; [xe,ye] = m e s h g r i d ( g r i d ) ;
9 epoints = [ x e ( : ) y e ( : ) ] ;
7o Compute d i s t a n c e m a t r i x between e v a l u a t i o n p o i n t s and c e n t e r s
10 DM_eval = D i s t a n c e M a t r i x ( e p o i n t s , c t r s ) ;
7. Compute d i s t a n c e m a t r i x between t h e data s i t e s and c e n t e r s
11 DM_data = D i s t a n c e M a t r i x ( d s i t e s , c t r s ) ;
12 for i=l:length(ep)
7, Compute i n t e r p o l a t i o n m a t r i x
13 IM = r b f ( e p ( i ) , D M _ d a t a ) ;
7o Compute e v a l u a t i o n matrix
14 EM = r b f ( e p ( i ) , D M _ e v a l ) ;
7o Compute power f u n c t i o n a t e v a l u a t i o n p o i n t s
15 invIM = i n v ( I M ) ; phiO = r b f ( e p ( i ) , 0 ) ;
16 f o r j=l:M
17 powfun(j) = r e a l ( s q r t ( p h i O - ( i n v I M * E M ( j , : ) ' ) ' * E M ( j , : ) ' ) ) ;
18 end
7o Compute max. norm of power f u n c t i o n on e v a l u a t i o n g r i d
19 maxPF(i) = max(powfun);
20 end
21 f p r i n t f (' S m a l l e s t maximum norm: 7oe\n' , min(maxPF))
22 f p r i n t f ('at e p s i l o n = 7.f \n' , ep(maxPF==min(maxPF) ) )
&4
17. Numerical Evidence for Approximation Order Results 145
0 5 10 15 20 0 5 10 15 20
e e
Fig. 17.2 Optimal e curves based on power functions for Gaussians in I D (left) and 2D (right)
for various choices of N uniform points.
ID 2D
N optimal e cond(A) N optimal e cond(A)
p}k*) = i>?V(ii*-*jii),
cf
such t h a t
Vfixi) = f it i = 1 , . . . , k - 1, k + 1 , . . . , N,
a n d i f Ek is t h e e r r o r
E k = fk-Vf{x )
k
at t h e one p o i n t Xk n o t used t o d e t e r m i n e t h e i n t e r p o l a n t , t h e n t h e q u a l i t y o f t h e
T
fit is d e t e r m i n e d b y t h e n o r m o f t h e v e c t o r o f errors E \E\,..., Ejy} obtained
by r e m o v i n g i n t u r n one o f t h e d a t a p o i n t s a n d c o m p a r i n g t h e r e s u l t i n g f i t w i t h
t h e ( k n o w n ) value at t h e r e m o v e d p o i n t . I n [ R i p p a (1999)] t h e a u t h o r presented
examples based o n use o f t h e l\ and i 2 n o r m s . W e w i l l m o s t l y use t h e m a x i m u m
n o r m (see line 14 i n t h e code b e l o w ) .
B y a d d i n g a l o o p over e we c a n c o m p a r e t h e e r r o r n o r m s for different values o f
t h e shape p a r a m e t e r , a n d choose t h a t value o f e t h a t y i e l d s t h e m i n i m a l e r r o r n o r m
as t h e o p t i m a l one.
W h i l e a naive i m p l e m e n t a t i o n o f t h e leave-one-out a l g o r i t h m is r a t h e r expensive
4
(on t h e order o f N o p e r a t i o n s ) , R i p p a showed t h a t t h e c o m p u t a t i o n o f t h e e r r o r
c o m p o n e n t s can be s i m p l i f i e d t o a single f o r m u l a
Ek = (17-1)
A
kk
w h e r e Ck is t h e A;th coefficient i n t h e i n t e r p o l a n t Vf based o n t h e full d a t a set, a n d
is t h e kth d i a g o n a l element o f t h e inverse o f t h e c o r r e s p o n d i n g i n t e r p o l a t i o n
17. Numerical Evidence for Approximation Order Results 147
P r o g r a m 1 7 . 2 . L00CV2D.m
% L00CV2D
% S c r i p t t h a t performs leave-one-out c r o s s - v a l i d a t i o n
/ (Rippa's method) t o f i n d a good e p s i l o n f o r 2D RBF i n t e r p o l a t i o n
0
% C a l l s on: D i s t a n c e M a t r i x
1 r b f = @(e,r) e x p ( - ( e * r ) . " 2 ) ; % Gaussian RBF
% Parameters f o r shape parameter loop below
2 mine = 0; maxe = 20; ne = 500;
3 ep = linspace(mine,maxe,ne);
% Number and type of d a t a p o i n t s
4 N = 81; gridtype = 'u';
7, Define t e s t f u n c t i o n
5 t e s t f u n c t i o n = @(x,y) s i n e ( x ) . * s i n c ( y ) ;
% Load data p o i n t s
,
6 name = s p r i n t f ( Data2D_/ d%s',N,gridtype) ; load(name)
0
17 f p r i n t f ( a t e p s i l o n = 7f \n',ep(maxEF==min(maxEF)) )
J
/ P l o t cost f u n c t i o n norm
0
18 f i g u r e ; semilogy(ep,maxEF,'b');
d i m e n s i o n a l e x p e r i m e n t (left) a n d a 2 D e x p e r i m e n t ( r i g h t ) based o n d a t a s a m p l e d
f r o m t h e sine f u n c t i o n F\. E a c h p l o t shows several curves c o r r e s p o n d i n g t o different
choices o f N (set o n l i n e 4 o f L00CV2D.m). T h e o p t i m a l e values are l i s t e d i n
T a b l e 17.3.
0 5 10 15 20 0 5 10 15 20
e e
Fig. 17.3 Optimal e curves based on leave-one-out cross validation for interpolation to the sine
function with Gaussians in I D (left) and 2D (right) for various choices of N uniform points.
ID 2D
N optimal e N optimal
3 0.96 9 0.96
5 1.00 25 1.00
9 0.80 81 1.48
17 0.92 289 1.60
33 1.92
65 1.76
10 10"
W -5 Of "
10
10 5
10" 10"
10 15 20 10 15 20
e
Fig. 17.4 Optimal e curves based on interpolation to the test function with Gaussians for
various choices of N uniform points. Trial and error approach (left), leave-one-out cross validation
(right).
10
10
10 y /
f 4
/ 4*
10" _ // * N=3
* N=5
N=9
10 t N=17
N=33
N=65
10
10 15 20
Fig. 17.5 Optimal e curves based on leave-one-out cross validation for interpolation to I D Franke's
function with Wendland's function 953,1 for various choices of N uniform points (left) and Cheby-
shev points (right).
P r o g r a m 17.3. CostEpsilon.m
7. ceps = C o s t E p s i l o n ( e p , r , r b f , r h s )
% Implements cost f u n c t i o n f o r o p t i m i z a t i o n of shape parameter
7, v i a Rippa's L00CV algorithm
% Example of usage i n L00CV2Dmin.m
1 f u n c t i o n ceps = C o s t E p s i l o n ( e p , r , r b f , r h s )
2 A = rbf(ep,r);
3 invA = p i n v ( A ) ;
4 EF = ( i n v A * r h s ) . / d i a g ( i n v A ) ;
5 ceps = n o r m ( E F ( : ) , i n f ) ;
P r o g r a m 17.4. L00CV2Dmin.m
'/. L00CV2Dmin
% S c r i p t t h a t performs leave-one-out c r o s s - v a l i d a t i o n
7 5
( R i p p a s method) to f i n d a good e p s i l o n f o r 2D RBF i n t e r p o l a t i o n
[Link] the help of MATLAB' s fminbnd
% C a l l s on: D i s t a n c e M a t r i x
7,Requires: CostEpsilon
1 r b f = @(e,r) e x p ( - ( e * r ) . ~ 2 ) ; 7o Gaussian RBF
7o Parameters f o r shape parameter o p t i m i z a t i o n below
2 mine = 0; maxe = 20;
7o Number and type of data p o i n t s
3 N = 81; g r i d t y p e = 'u';
7o Define t e s t f u n c t i o n
4 t e s t f u n c t i o n = @(x,y) s i n e ( x ) . * s i n c ( y ) ;
7o Load data p o i n t s
5 name = s p r i n t f ( 'Data2D_7od7 s ' ,N,gridtype) ; load(name)
0
j=i
for a fixed e v a l u a t i o n p o i n t x as a n a n a l y t i c f u n c t i o n o f e.
T h e key idea is t o represent Vf(x, e) b y a Laurent series i n e, a n d a p p r o x i m a t e
t h e "negative p a r t " o f t h e series b y a P a d e a p p r o x i m a n t , i.e.,
oo
d ek
V (x,e)
f ^r(e) + J2 k >
k=0
where r(e) is t h e r a t i o n a l P a d e a p p r o x i m a n t .
We t h e n r e w r i t e t h e i n t e r p o l a n t i n c a r d i n a l f o r m , i.e., as
Vf{x,e) = ^C V {\\X
3 - Xj\\)
j=i
T
= b (x, e)c
T
= b {x,e)A-\e)f
= (u*(x,e)ff
T
where b(x,e) - 3 = ip (\\x
- xj\\), A(e)ij = <p (\\xi
- x - \ \ ) , c = [a,...,
3 c],
N f =
T
[fi,-- -,fN] , and
_1
u*(x,e) = A ()6(a;,)
0 5 10 15 20 0 5 10 15 20
E E
Fig. 17.6 Optimal e curves based on Contour-Pade for interpolation to the sine function with
Gaussians in I D (left) and 2D (right) for various choices of N uniform points.
Table 17.4 Optimal e values based on Contour-Pade for interpolation to the sine
function with Gaussians in I D and 2 D for various choices of N uniform points.
ID 2D
17.1.5 Summary
uniform Halton
9 3.195983e-001 2.734756e-001
25 5.008591e-002 2.6738 8.831682e-002 2.3004
81 9.029664e-003 2.4717 2.401868e-002 1.7582
289 2.263880e-004 5.3178 1.589117e-003 5.0969
1089 3.323287e-008 12.7339 1.595051e-006 10.8015
4225 1.868286e-008 0.8309 9.510404e-008 4.8203
uniform Halton
9 3.302644e-001 2.823150e-001
25 3.271035e-002 3.3358 1.282572e-001 1.6058
81 1.293184e-002 1.3388 3.407580e-002 1.7898
289 3.786113e-004 5.0941 1.990217e-003 5.3309
1089 3.476835e-008 13.4107 2.286014e-006 10.5905
4225 3.775365e-008 -0.1188 9.868530e-008 5.3724
Fig. 17.7 Maximum errors for non-stationary interpolation to the sine function with Gaussians
(left) and inverse multiquadrics (right) based on N uniformly spaced points in [0,1] and e
1,6,11,16,21,26.
1 / / / l
a n d for inverse m u l t i q u a d r i c s we have s p e c t r a l convergence w i t h h i e ~ . We
c a