Meshless cubature over the disk

by Thin-Plate Splines ?

Alessandro Punzi, Alvise Sommariva a,∗, Marco Vianello a

a Department of Pure and Applied Mathematics, University of Padua (Italy)
and the (symmetric) augmented matrix
 
A B
A= , A = AX,φ = [φj (Pi )] , B = [πk (Pi )] , (6)


BT 0

it is easy to see that the cubature formula can be written in the usual form of
weighted sum of the function values,
n
−1 −1
X
I(s) = hA F, Ii = hA I, Fi = hW, Fi = hw, f i = wj fj , (7)
j=1

h·, ·i denoting the scalar product in Rn . Thus, as usual with interpolatory

formulas, the weights w = {wj } are obtained by solving the linear system
 
w
AW = I , W =   (weights equations) . (8)
z

We give now an error estimate, which takes into account the presence of errors
in the evaluation of the integrals, Ĩ ≈ I, and of possible noise in the sample,
f̃ ≈ f ; the relevant parameters are the fill distance h, and the separation
distance q of the point set X (recall that q ≤ 2h). The cubature error can be
bounded as

|I(f ) − hw̃, f̃ i| ≤ |I(f ) − I(s)| + |hw − w̃, f i| + |hw̃, f − f̃ i|

q
≤ meas(Ω) kf − skL2 (Ω) + kA−1 k2 kf k2 kI − Ĩk2 + kf − f̃ k∞ kw̃k1
n
!
1
q  X
=O F φ (h) + O kI − Ĩk2 + kf − f̃ k∞ |w̃j | , (9)
G φ (q) j=1
where F φ (h) ↓ 0 as h → 0 and G φ (q) ↓ 0 as q → 0, and thus as h → 0.
Technically, such an estimate holds in the “native space” of the RBF, and
requires that Ω satisfies the so-called “interior cone condition”, cf. [2,5,14].
q
The simultaneous appearance of the terms F φ (h) ↓ 0 and 1/G φ (q) ↑ +∞ in
(9) is an occurrence, in the context of numerical cubature, of the well-known
“uncertainty principle” in RBF interpolation, cf. [7]. In practice, however,
there are two distinct situations: for smooth RBF, like Gaussians and (inverse)
Multiquadrics, the rates of F φ (h) and G φ (q) are both exponential , while for
less regular RBF, like TPS, the rates are both algebraic; cf. e.g. [14].

In Table 1, we report the two key parameters related to stability of RBF cuba-
ture, for (inverse) Multiquadrics (MQ and IMQ), Gaussians (G), Wendland’s
compactly supported (W2) and Thin-Plate Splines (TPS). For simplicity, we
consider only unscaled RBF (δ = 1).

3
Table 1
RBF cubature with sets of n = 50 and n = 100 uniform random points in [0, 1] 2 :
spectral norm of the inverses of the collocation matrices and 1-norm of the com-
puted weights vectors (average values on 50 independent trials, rounded to the first
significant digit).
] of rnd pts norms MQ IM Q G W2 TPS
n = 50 kA −1
k2 2E+12 3E+11 5E+15 5E+03 6E+03
kw̃k1 7E+01 9E+01 3E+02 2E+00 1E+00
n = 100 kA−1 k2 2E+16 6E+15 3E+17 5E+04 8E+03
kw̃k1 8E+02 5E+02 1E+03 2E+00 1E+00

The situation seems hopeless concerning the use of smooth RBF like Gaussians
for numerical cubature, in view of the expected exponential magnification of
the integration errors. However, the numerical experiments in [8] have shown
that (9) is by far an overestimate. In any case, it turns out that smooth RBF
are much more sensible to integration errors, increase of data density and per-
turbations in the data, whereas TPS and W2 give reasonably accurate and
more stable cubature formulas. Moreover, the quality of cubature by TPS is
weakly sensitive to the scaling parameter, a phenomenon which can be related
to the well-known scale independence of the condition number of interpolation
by polyharmonic splines (cf. [5, §3.8]). Clearly, this property makes TPS cuba-
ture very attractive for automatic integration, since it avoids the complication
of managing/optimizing the scaling parameter.

3 Meshless cubature over the disk.

In [8,9] it has been possible to integrate the basis functions by means of sym-
bolic computation (see [13]), obtaining closed-form formulas for the integrals
{I(φj )}. In fact, concerning integration of RBF over the square, one can re-
sort to radial symmetry (cf. [8]), and more generally integration of TPS over
polygons can be accomplished via Green’s formula, computing a first-level
primitive in one of the variables, and then a second-level primitive along each
side (cf. [9]). In the latter case, the success of the procedure is due to the
peculiarity of the TPS basis and to the fact that the parametrization of each
side is an affine mapping.

Concerning integration of TPS over the disk, the situation is more delicate.
We begin by recalling Green’s formula (cf. [4] and, e.g., [1]), specialized to the
case of the unit disk, i.e.

• Ω = D = {P = (x, y) : x2 + y 2 ≤ 1}

which is not restrictive, in view of the obvious “shift and scale” change of
variables which maps any disk into the unit disk (interpolation and cubature

4
by TPS are weakly sensitive to scaling, as already observed). Now, recall that
we are dealing with

• Duchon’s Thin-Plate Splines (TPS): φ(r) = r 2 log (r).

Setting for j = 1, . . . , n
Z
ψj (P ) = φ(|P − Pj |) dx , (10)

Green’s formula over the disk reads as

Z I Z π
I(φj ) = φ(|P − Pj |) dP = ψj (P ) dy = ψj (cos θ, sin θ) cos θ dθ ,
D ∂D −π
(11)
where the first-level primitive ψj (P ) can be computed explicitly by the pow-
erful symbolic integrator [15],

ψj (P ) = ψj1 (P ) + ψj2 (P ) + ψj3 (P ) , (12)

where
1 2
ψj1 (P ) = − (x − xj )3 − (x − xj )(y − yj )2 ,
9 3
1  
ψj2 (P ) = (x − xj )((x − xj )2 + 3(y − yj )2 ) log (x − xj )2 + (y − yj )2
6
!
2 x − x j
ψj3 (P ) = (y − yj )3 arctan . (13)
3 y − yj

Unfortunately, here symbolic integrators (like Mathematica [13] or the Maple

kernel of Matlab Symbolic Toolbox [6]) fail in computing the second-level
primitive “along the boundary”, i.e. ψj (cos θ, sin θ) cos θ dθ, and thus the last
R

integral in (11) has to be computed numerically. The difficulties are in the

second and third term in (13), whereas the first can be integrated symbolically.
On the contrary, integration on the disk of the polynomial part in (3) (here
m = 1) is trivial in polar coordinates. We stress that, in any case, we have
reduced a 2-dimensional to a 1-dimensional integration problem, with a clear
computational advantage.

Various choices are possible. In view of the periodicity, the simple trapezoidal
rule is appealing, but we have to take into account that both ψj2 and ψj3 have
low regularity. Indeed, ψj2 has a critical point at P = Pj , and ψj3 has a whole
critical line, y = yj . These generate one possible critical point for ψj2 on the
circle (when Pj is on or close to the circle), and two critical points for ψj3 on
the circle (the intersections with the critical line). Then we chose two different
parametrizations in such a way that the critical points are endpoints of the
parameter intervals, and quadrature formulas that automatically cluster nodes
at the endpoints, like Gauss-Legendre or Clenshaw-Curtis.

5
The one-dimensional integrals to be computed are
I Z θj +2π
ψj2 (P ) dy = ψj2 (cos θ, sin θ) cos θ dθ , θj = arg (Pj ) , (14)
∂D θj

and by splitting into two arcs

I Z max {θj+ ,θj− }
ψj3 (P ) dy = ψj3 (cos θ, sin θ) cos θ dθ
∂D min {θj+ ,θj− }

Z min {θj+ ,θj− }+2π  q 

+ ψj3 (cos θ, sin θ) cos θ dθ , θj± = arg ± 1 − yj2 , yj .
max {θj+ ,θj− }
(15)
Since the computation of Gauss-Legendre nodes and weights is expensive at
high degree of exactness, whereas Clenshaw-Curtis formula can be imple-
mented with low cost by the FFT and has in practice a comparable accu-
racy (cf. [12]), we have used the latter. We recall that when several different
functions have to be integrated (as here), the formulation of Clenshaw-Curtis
formula involving explicitly the weights is much more convenient. The weights
can be computed once and for all by the fast algorithm recently proposed in
[13] (based on the IFFT). We have numerical evidence that the Clenshaw-
Curtis formula of degree 32 is sufficient to compute all the three integrals in
(14)-(15) (and then I(φj )) with 10 correct figures. Observe that 2-dimensional
numerical integration of φj would require a much larger number of nodes to
get the same accuracy.

0.8

0.6

0.4

0.2

−0.2

−0.4

−0.6

−0.8

−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8

0.5

0
−3 −2 −1 0 1 2 3

Fig. 1. The nonuniformity effect on a uniform random sample of 1000 points on the
unit disk by mapping in the polar rectangle [−π, π] × [0, 1].

Remark 1 (Polar or cartesian coordinates?) Following [9], we term the ap-

proach described above TPS-Green meshless cubature (in cartesian coordi-
nates). One might think, in principle, that the whole construction is not nec-
essary, since integration over the disk is integration over a rectangle in polar

6
coordinates, and thus one could simply work from the beginning in such co-
ordinates using the cubature method developed in [8]. In such a case, clearly,
the TPS would be themselves TPS in polar coordinates. However, if we start
for example from a random sample generated by the uniform probability dis-
tribution in cartesian coordinates, the mapping from the disk to the polar
rectangle causes a distortion (see Figure 1): the density of the mapped points
decreases approaching the center of the disk, and this deteriorates the quality
of the approximation. The effect of such nonuniformity will be manifest in the
numerical experiments of the next section (compare “TPS-Green polar” with
“TPS-Green cartesian”).

Remark 2 (Possible extensions.) The TPS-Green meshless cubature method

is not restricted to polygons (cf. [9]), or disks. It could be extended to any
bounded domain, even multiply connected, with piecewise regular external
and internal boundaries (∂Ω = Γext ∪ Γint ), namely in computing
I I I
I(φj ) = ψj (P ) dy = ψj (P ) dy − ψj (P ) dy
∂Ω Γext Γint

XZ XZ
= ext
ψj (P ) dy − ψj (P ) dy , (16)
k Γk s Γint
s

provided that the boundary pieces {Γext int

k } and {Γs } are suitably parametrized.
Clearly, apart from the case of a linear piece of the boundary (segment) where
explicit formulas are known (cf. [9]), integration of ψj along the boundary
has to be done numerically. For example, extension to ellipses is immediate
by integrating ψj (a cos θ, b sin θ) cos θ, where a and b are the semiaxes. Again,
one can use Clensahw-Curtis quadrature as above. In such a way we can
work directly with TPS on the ellipse, avoiding the distortion of the sampling
point set occurring via transformation of the ellipse into a square or into a
disk (change of variables in the integration). We have also implemented the
extension in the case of (annular) circular sectors, i.e.

Ω = S(ρ1 , ρ2 , α1 , α2 ) = {P : 0 ≤ ρ1 ≤ |P | ≤ ρ2 , α1 ≤ arg(P ) ≤ α2 } , (17)

whose boundary consists of the two circular arcs |P | = ρi and of the two
segments arg(P ) = αi , i = 1, 2. Here, integration of ψj (P ) along the segments
is explicit, whereas Clenshaw-Curtis quadrature along each arc (the internal
clockwise and the external counterclockwise) is possibly split into two or three
subarcs, to manage as above the critical points.

4 Numerical tests.

In this section we present several numerical tests of numerical integration

from scattered data over the disk, by TPS-Green meshless cubature. First, we

7
show the behavior of the formulas on three test functions, (already used in [9]
concerning meshless cubature over polygons), dealing with samples of small
size (in the hundreds). Then, we show how to improve efficiency of TPS-Green
cubature on samples of larger size (in the thousands) by splitting the disk into
sectors or annuli. All the tests have been done in Matlab (cf. [6]), with an Intel-
Centrino Duo 2.3Ghz processor and 1Gb RAM. The corresponding Matlab
code can be found at [10]. Our implementation is truly meshless (no kind of
mesh is required or generated), and it needs only the parameters defining the
disk (or more generally the annular sector, i.e. the center, and the radial and
angular intervals) besides the scattered sample.

We have considered the following test functions

q
x2 + y 2 .
f1 (x, y) = exp (x − y) , f2 (x, y) = exp (5(x − y)) , f3 (x, y) =
(18)
∞
Observe that f1 and f2 are C (with f2 varying more rapidly), whereas f3
has a singularity of the gradient in the origin.

In Table 2 we show the parameters relevant to stability of cubature by RBF, i.e.

the spectral norm of the collocation matrix inverse, and the sum of the weights
absolute values (cf. estimate (9)). It is worth observing that the weights are
not all positive, but the sum of the negative ones is, in absolute value, much
smaller than that of the positive ones. Thus, the sum of the weights absolute
values is always of the order of the sum of the weights, which approximates π
(the unit disk area).

In Table 3 we compare the errors of TPS-Green cubature in cartesian and po-

lar coordinates (see Remark 1), and of Monte Carlo formula, on a sequence of
scattered samples in the unit disk obtained by the uniform random distribu-
tion in the enclosing square. The linear systems involved (weights equations)
have been solved by the standard direct solver of Matlab. Observe that TPS-
Green cubature in cartesian coordinates is more accurate than that in polar
coordinates, due to the distortion effect described in Remark 1, and also than
Monte Carlo integration (up to two orders of magnitude). This latter fact is
remarkable thinking to practical applications, where sampling is very costly
or even cannot be refined, and thus a reasonably accurate value of the integral
should be obtained from the available data.
Table 2
Spectral norm of the inverses of the collocation matrices and 1-norm of the computed
weights vectors (cartesian coordinates).
norms 100 pts 200 pts 400 pts 800 pts
kA−1 k2 2E+04 1E+05 2E+05 3E+05
kw̃k1 3.56 3.55 3.44 3.65

These and many other numerical tests (not reported for brevity) show that

8
Table 3
Relative errors of TPS-Green cubature in polar and cartesian coordinates and of
Monte Carlo cubature with n = 100, 200, 400, 800 uniform random points on the
unit disk.
function formula 100 pts 200 pts 400 pts 800 pts
f1 TPS-Green (polar) 3E-03 2E-03 2E-03 3E-04
TPS-Green (cart.) 1E-03 1E-04 1E-05 6E-06
Monte Carlo 4E-02 1E-02 2E-02 1E-03
f2 TPS-Green (polar) 2E-02 1E-02 1E-02 6E-03
TPS-Green (cart.) 3E-02 2E-02 2E-03 6E-04
Monte Carlo 3E-01 1E-01 4E-02 4E-03
f3 TPS-Green (polar) 3E-02 1E-02 1E-02 6E-03
TPS-Green (cart.) 5E-04 4E-04 7E-05 8E-06
Monte Carlo 1E-01 4E-02 2E-02 1E-02

TPS-Green cubature provides a flexible, reasonably accurate and stable mesh-

less method on disks, sectors and annuli, which gives good results even with
relatively small scattered samples (say a size in the hundreds of points). Al-
ready with “moderate” size samples (say a size in the thousands), however, we
face the classical complexity problem of globally supported RBF with direct
solvers. Dealing with interpolation, a number of fast methods have been devel-
oped to accelerate the solution of the collocation equations, typically within
the frame of iterative solvers; cf., e.g., [2,14] and references therein.

In the case of numerical cubature, however, there is a simple alternative for

improving efficiency, even using a direct solver, at least until conditioning re-
mains acceptable. Due to additivity of the integral, we can split the integration
domain, and hence the data, into subdomains, like (sub)sectors or (sub)annuli.
The integral on each subdomain can be computed with TPS-Green cubature,
where each local system of weights equations can be much smaller than the
original global one, depending on the number of subdomains and on the dis-
tribution of the sampling points.

In Table 4 we show an example of the behavior of TPS-Green cubature with

the data splitting technique just described, by subdividing the unit disk into
16 equal area annuli. Observe that with no splitting, the bulk of the compu-
tational process is clearly given by the assembling and solution of the weights
equations. When the splitting is applied, the cost of the weights equations is
pulled down, whereas the cost of integration of the TPS is roughly doubled
(each ψj is integrated along two circles on the average, namely the boundaries
of an annulus). The resulting effect is a remarkable reduction of the global
CPU time (speed-up ratio about 9 in this example).

9
Table 4
Relative errors and CPU times (seconds): TPS-Green cubature of f (x, y) =
exp (5(x2 + y 2 )) by 3000 uniform random points on the unit disk, directly and by
subdividing into 16 equal area annuli.
subdivisions error CPU weights CPU basis cub CPU tot speed-up
none 5E-4 17.3 1.1 18.4
16 1E-4 0.3 1.8 2.1 8.8

References

[1] T.M. Apostol, Calculus, vol. II, 2nd edition, Blaisdell, 1969.

[2] M.D. Buhmann, Radial Basis Functions: Theory and Implementation,

Cambridge Monographs on Applied and Computational Mathematics, vol. 12,
Cambridge University Press, 2003.

[3] R. Cools and D. Laurie (Eds.), Numerical evaluation of integrals, J. Comput.

Appl. Math. 112 (1999), no. 1-2.

[4] G. Green, An Essay on the Application of Mathematical Analysis to the Theories

of Electricity and Magnetism, Nottingham, 1828.

[5] A. Iske, Multiresolution Methods in Scattered Data Modelling, Lecture Notes in

Computational Science and Engineering, vol. 37, Springer, 2004.

[6] The MathWorks, MATLAB Documentation Set, 2006 version (available online
at http://www.mathworks.com).

[7] R. Schaback, Error estimates and condition numbers for radial basis function
interpolation, Adv. Comput. Math. 3 (1995), 251–264.

[8] A. Sommariva and M. Vianello, Numerical cubature on scattered data by radial

basis functions, Computing 76 (2005), 295–310.

[9] A. Sommariva and M. Vianello, Meshless cubature by Green’s formula, Appl.

Math. Comput. 183 (2006), 1098–1107.

[10] A. Sommariva and M. Vianello, GreenDisk: Matlab code for meshless

cubature over disks by Thin-Plate Splines, November 2006 (downloadable at
http://www.math.unipd.it/∼marcov/software.html).

[11] A. Sommariva and R. Womersley, Integration by RBF over

the sphere, 2006, submitted (preprint UNSW AMR05/17, available online at
http://www.maths.unsw.edu.au).

[12] L.N. Trefethen, Is Gauss quadrature better than Clenshaw-Curtis? , SIAM Rev.,
to appear.

[13] J. Waldvogel, Fast construction of the Fejér and Clenshaw-Curtis quadrature

rules, BIT 46 (2006), 195–202.

10
[14] H. Wendland, Scattered Data Approximation, Cambridge Monographs on
Applied and Computational Mathematics, vol. 17, Cambridge University Press,
2005.

[15] Wolfram Research, Inc., The Wolfram Integrator , 2006 version (available online
at http://integrals.wolfram.com).

11