We present an algorithm to approximate large dataset by Radial Basis Function
(RBF) techniques. The method couples a fast domain decomposition procedure with a
localized stabilization method. The resulting algorithm can efficiently deal with large
problems and it is robust with respect to the typical instability of kernel methods

Cavoretto, R.

DE MARCHI, Stefano

De Rossi, A

Perracchione, E

Santin, Gabriele

English

Cavoretto, Roberto

DE ROSSI, Alessandra

Perracchione, Emma

Institutional Research Information System University of Turin

Proceedings of the 15th International Conferenceon Computational and Mathematical Methodsin Science and Engineering, CMMSE 20156–10 July, 2015.RBF approximation of large datasets by partition of unityand local stabilizationR. Cavoretto1, S. De Marchi2, A. De Rossi1, E. Perracchione1 and G.Santin21 Department of Mathematics “G. Peano”, University of Torino2 Department of Mathematics, University of Padovaemails: roberto.cavoretto@unito.it, demarchi@math.unipd.it,alessandra.derossi@unito.it, emma.perracchione@unito.it,gsantin@math.unipd.itAbstractWe present an algorithm to approximate large datasets by Radial Basis Function(RBF) techniques. The method couples a fast domain decomposition procedure with alocalized stabilization method. The resulting algorithm can efficiently deal with largeproblems and it is robust with respect to the typical instability of kernel methods.Key words: Approximation, Radial Basis Functions, partition of unity, stabilization1 IntroductionRadial Basis Function (RBF) methods provide a flexible way to recover multivariate func-tions from scattered data. In the recent years, such techniques have been successfully appliedin various fields of applied mathematics, computer science and engineering to solve severalproblems, such as approximation of functions and numerical solutions of PDEs.Although effective, these methods require special attention when applied to large prob-lems. Namely, difficulties often arise on the side of the computational cost, due to the use ofusually dense matrices, and on the side of stability of the algorithms, and these limitationsare sometimes an obstacle to their application on concrete and real-world problems.c©CMMSE ISBN: 978-84-617-2230-3Page 317 of  1268PU-WSVDAs regards stability, it is known that the cause of ill-conditioning is the particularalgorithm used to construct the solution. In particular, it has been proven that the approx-imation operator itself is stable in the function space associated with the RBF kernel (see[6]), while the instability comes from the use of the unstable standard basis of translates. Tothis end, in recent years a lot of efforts have been played to study and construct better con-ditioned bases. A general theoretical framework for change of basis in kernel-based spacesis developed in [12]. For specific kernels, namely the Gaussian and inverse multiquadrickernels, almost perfectly conditioned bases have been constructed ([8, 9, 10, 7]), and somemethod to produce better bases for a generic kernel have been developed too (see [4, 5, 11]).To overcome the aforementioned limitations, we propose a new method based on thecoupling of a fast domain decomposition technique with a local stabilization method relyingon the change of basis presented in [4]. This stabilization method have been chosen thanksto its applicability, without restriction, to any kernel. The algorithm presented here isintended as an exploratory attempt to couple such methods, and we plan to improve it byreplacing the domain decomposition/points search technique used in this paper with therecently developed block search algorithms (see [2, 3]).2 Kernel based approximationAssuming to have a finite set of n distinct data points Xn = {xi, 1 ≤ i ≤ n} in a boundeddomain Ω ⊂ Rd and a set of data values Fn = {fi, 1 ≤ i ≤ n}, our goal is to recover acontinuous function Rn : Ω→ R that interpolates the data. To construct Rn one considersa positive definite and symmetric kernel K : Ω×Ω→ R and looks for a solution in the formRn =n∑i=1ciK(·, xi). (1)The coefficient vector c = [c1, . . . , cn]T is determined by imposing the interpolation conditionRn(xi) = fi, 1 ≤ i ≤ n, i.e., by solving the linear system Ac = f , where A = [K(xi, xj)]ni,j=1is the kernel matrix and f = [f1, . . . , fn]T . Since the kernel is positive definite, so it is thematrix A, hence the solution c exists and it is unique for any given data. The kernels wewill deal with are always radial, meaning that there exist a univariate function φ : R≥0 → Rand a shape parameter ε > 0 such that K(x, y) = φ(ε‖x− y‖2) for all x, y ∈ Ω. The shapeparameter has the effect of localize or flatten the kernel, and it has a crucial role in thebehavior of the approximant.We recall that associated to a positive definite kernel K on Ω there is a uniquely definedreproducing kernel Hilbert space NK(Ω), called the native space of K on Ω. When the datacomes from the sampling of a function f ∈ NK(Ω), i.e., fi = f(xi), the operator f → Rnis the projector from NK(Ω) into the subspace Vn = span{K(·, xi) : 1 ≤ i ≤ n}. Forf ∈ NK(Ω) we have Rn → f as Xn gets dense into Ω, and explicit convergence estimatesc©CMMSE ISBN: 978-84-617-2230-3Page 318 of  1268R. Cavoretto, S. De Marchi, A. De Rossi, E. Perracchione, G. Santinare known (see e.g. [16]), depending on the smoothness of K and Ω and on the distributionof Xn in Ω.Nevertheless, the actual computation of Rn is often numerically intractable using thestandard method described above. In particular, the kernel matrix is generally dense andof the same size n of the data, hence the construction and the solution of the linear systemis computationally expensive. Moreover, it is in general required to tune the parameter εto fit the data, but reducing it is a source of instability, especially when K is smooth.We will see in Section 3 how to deal with large data sets by splitting the problem insmaller subproblems, and in Section 4 how to locally employ the change of basis to controlthe instability.3 Domain decompositionWe want to avoid to solve a global problem involving a large and possibly dense matrix. Todo so, we focus on the partition of unity (PU) method (see [15]).The idea is to decompose the problem, solve (many) small local approximation sub-problems, and then blend together the results in a global approximant. We give here ageneral review of this method, while further implementation details are provided in Section5.We consider an open and finite covering {Ω(k)}Nk=1 of the bounded set Ω, that isΩ ⊂N⋃k=1Ω(k),where we require mild overlapping of the patches. Associated to the covering one considersa partition of unity, i.e., a set of nonnegative and continuous functions {W (k)}Nk=1 withsupp(W (k)) ⊂ Ω(k),N∑k=1W (k)(x) = 1 for all x ∈ Ω.With these tools in hand, if we are able to produce a local approximant R(k) for any subsetΩ(k), it is immediate to recover a global approximant Rn on Ω asRn(x) =N∑k=1W (k)(x)R(k)(x), x ∈ Ω.But the kernel K is still positive definite when restricted to each Ω(k), so it is possibleto uniquely solve the approximation problem restricted to the local domains. Namely,for each Ω(k), 1 ≤ k ≤ N , we consider the restricted set of n(k) data locations defined asc©CMMSE ISBN: 978-84-617-2230-3Page 319 of  1268PU-WSVDX(k)n = Xn∩Ω(k), n(k) = card(X(k)n ), and the corresponding values F (k)n = {f(x) : x ∈ X(k)n }.The local approximant R(k) is then computed as in (1), i.e,R(k) =n(k)∑i=1c(k)i K(·, x(k)i ),where now it suffices to solve a much smaller linear system to determine the coefficients inthe above expansion.It is important to recall that, under further assumptions on the covering and on thepartition of unity, it is possible to prove that the global approximant keeps the same ap-proximation order of the local ones (see [15]), but we will not go into the details here.Instead we emphasize that, besides the theoretical framework, the success of the PUmethod relies on a good point searching technique. Although we will employ in Section5 a standard algorithm, we plan to improve the method by updating the current domaindecomposition method with the recently developed block-based algorithms (see [2, 3]). Westress here that this update will presumably lead to a reduction of the computational cost,while nothing will change in terms of accuracy. In this view, the numerical experimentspresented in Section 5 should be indicative, as regards the error, of the behavior of theenhanced method.4 Local stabilizationThe partition of unity method allows to deal with large datasets in an efficient way, andpartially reduces the ill-conditioning found in many cases when dealing with RBF problems,since it reduces the size of the problem. Nevertheless, especially in the flat limit ε → 0,using a domain decomposition technique can be not enough. To this end, we want toemploy locally a stabilization method that permits to face in an even more stable way largeapproximation problems. The coupling of the PU algorithm with the local stabilization willlead to the full algorithm of the next Section.For any subdomain Ω(k) there is an associated space of functions NK(Ω(k)) given bythe restriction of NK(Ω) to Ω(k). As in the global setting, the interpolation operator isthe projection into the finite subspace V(k)n = span{K(·, x(k)j ) : x(k)j ∈ X(k)n } generated bythe restricted set of data X(k)n . To this subspace we apply the change of basis described in[4, 5]. Namely, we substitute the standard basis {K(·, x(k)i ) : x(k)i ∈ X(k)n } with the WeightedSingular Value Decomposition (WSVD) basis {u(k)j }n(k)j=1 ,u(k)j =n(k)∑i=1d(k)ij K(·, x(k)i ).c©CMMSE ISBN: 978-84-617-2230-3Page 320 of  1268R. Cavoretto, S. De Marchi, A. De Rossi, E. Perracchione, G. SantinThe invertible coefficient matrix D(k) = [d(k)ij ]n(k)i,j=1 is computed through a singular valuedecomposition of the local kernel matrix A(k) = [K(x(k)i , x(k)j )]n(k)i,j=1, i.e.,A(k) = U (k)Σ(k)(U (k))T , D(k) = U (k)(Σ(k))−1/2, Σ(k) = diag(σ(k)1 , . . . , σ(k)n(k)).This basis has been constructed to mimic in a discrete way the eigenbasis defined throughthe Mercer’s Theorem (see e.g. [13]), where the L2(Ω) inner product is replaced by itsdiscrete pointwise version ℓ2(Xn). In particular, as proven in [4], the basis enjoys thefollowing properties.Proposition 4.1 The WSVD basis {u(k)j }n(k)j=1 has the following properties:(i) it is NK(Ω(k))-orthonormal,(ii) it is ℓ2(X(k)n )-orthogonal with norms {σ(k)j }j,(iii) (u(k)j , f)ℓ2(X(k)n )= σ(k)j (u(k)j , f)NK(Ω(k)), ∀f ∈ NK(Ω(k)), 1 ≤ j ≤ n(k),(iv) σ(k)n(k)≥ . . . ≥ σ(k)1 > 0.(v)∑n(k)j=1 σ(k)j = n(k) φ(0).We want now to approximate the restriction of f ∈ NK(Ω) to Ω(k), say f (k), using this basis.Thanks to the orthogonality property (i), the local interpolant of f (k) can be representedasR(k) =n(k)∑j=1(f (k), u(k)j )NK(Ω(k))u(k)j .Moreover, each element of this new basis is associated to a singular value of the matrix A.Since we expect the singular values to decay very fast to zero, we can drop the last elementsof the basis to remove from the interpolant the terms responsible for the ill-conditioning.To be more precise, for any m(k) ≤ n(k), the truncated interpolantR(k)m(k)=m(k)∑j=1(f (k), u(k)j )NK(Ω(k))u(k)kis the solution of the least-squares problemming∈span{u(k)1 ,...,u(k)m(k)}‖f (k) − g‖ℓ2(X(k)n ).c©CMMSE ISBN: 978-84-617-2230-3Page 321 of  1268PU-WSVDIndeed, thanks to property (iii), we have for any m(k) ≤ n(k)R(k)m(k)=m(k)∑j=1(f (k), u(k)j )NK(Ω(k))u(k)j =m(k)∑j=1(σ(k)j )−1(f (k), u(k)j )ℓ2(X(k)n )u(k)j , (2)and the rightmost term in the above expression is the ℓ2(X(k)n ) projection of f (k) into thesubspace span{u(k)1 , . . . , u(k)m(k)}. Obviously, for m(k) = n(k) the latter approximant coincideswith the interpolant, while we can look for a suitable truncation index m(k) < n(k) toremove the elements of the basis corresponding to the singular values below some prescribedtolerance.The use of this basis allows to reduce the instability, but it has some drawbacks. Namely,in case of severe ill-conditioning the truncation approach may be over regularizing. More-over, the computation of the SVD is expensive, especially since we want to drop part of itsfactors. To overcome the second problem we will approximate the singular value decomposi-tion using the Lanczos method (see e.g [14]), as done in [5]. This approximation introducesa perturbation in the basis that is thoroughly analyzed in the cited paper, nevertheless,the form (2) of the approximant is still valid for the approximated basis. As regards theregularizing effect of the truncation, it has been observed that the method, when appliedglobally, is capable to deal with only a limited instability. In fact, when ε→ 0 or when thekernel is too smooth, the original algorithm needs to leave out too much elements of thebasis, and the resulting approximant may be meaningless. Instead, when repeatedly appliedto the small subdomains Ω(k), the regularizing effect simply improves the approximation,as shown in the next Section.5 Numerical experimentsWe present here two experiments to test our algorithm on two small problems in dimensiond = 2. For both examples we compare the behavior of two different RBF kernels K(x, y) =φ(ε‖x− y‖2) with different smoothness, namely the Gaussian φ(r) = exp(−ε2r2) (which isa C∞ kernel) and the C2 Mate´rn kernel φ(r) = exp(−εr)(1+ εr). In both examples the testfunction is the bivariate Franke’s test functionf(x, y) = 34e− 14((9x−2)2+(9y−2)2) +34e−( 149(9x+1)2+ 110(9y+1))+12e− 14((9x−7)2+(9y−3)2) −15e−((9x−4)2+(9y−7)2),approximated on the unit square Ω = [0, 1]2. The shape parameter ε varies on a log-spacedgrid in [10−4, 102]. The error is measured by means of the root mean squared error (RMSE)on a equally spaced grid of 40d points in Ω.c©CMMSE ISBN: 978-84-617-2230-3Page 322 of  1268R. Cavoretto, S. De Marchi, A. De Rossi, E. Perracchione, G. SantinIn both examples the PU method uses a covering of Ω based on balls {Ω(k)}Nk=1, wherethe centers {y(k)}Nk=1 of the balls lie on an equally spaced 2D grid in Ω and the radiiare the same for any k = 1, . . . , N . The points are organized through a kd-tree, as in[1]. We use n Halton points, where n = 65d, 129d, 257d, which correspond respectivelyto N = 16d, 32d, 64d. The common radius of the balls is chosen to be 1/√N , so thatthe boundary of every patch intersects the centers of the nearby patches. As a partitionof unity we use a superposition of C2 Wendland’s functions W (k)(x) = ψ(ν‖x − y(k)‖),ψ(r) = (1−νr)4+(4νr+1), where the shape parameter ν > 0 is scaled to control the support ofthe partition. We remark that these choices of Ω(k) and W (k) fulfills the requirement on thecovering and the partition in order to preserve globally the order of the local approximation.As a stopping rule for the locally employed Lanczos algorithm we use an extension ofthe rule proposed in [5]. Namely, since we know that the trace of the local kernel matrix A(k)is n(k) φ(0) (see property (v) in Proposition 4.1), we stop the iteration when the distancebetween the trace of the Lanczos matrix and n(k) φ(0) is below a fixed tolerance τ , hereτ = 10−14.The first example compares the approximation obtained with the new method with theone obtained with a global application of the change of basis described in Section 4. Sincewe use here a global approach, the experiment is limited to n = 65d. The results are shownin Figure 1. For a not too small shape parameter, for both kernels, the global method andthe new algorithm behaves in the same way. As ε → 0 the global method loses accuracy,while the local one is still able to compute an accurate and stable approximant.In the second example we compare the behavior of the PU method without local regu-larization with the new algorithm. Figure 2 presents the results for this test. As expected,there is a first phase where no instability is present and the two methods behaves exactlyin the same way. Instead, in the flat limit case the standard PU method starts to becomeunstable, while the stabilized method retain its accuracy and it is able to effectively com-pute the approximants. The effect is much more evident for the smoother Gaussian (whichis smoother) than for the Mate´rn kernel.Other tests have been performed for d > 2, with similar results.AcknowledgementsThe first, third and fourth authors have been partially supported by the funds of the Univer-sity of Torino, project ex 60% “Metodi numerici nelle scienze applicate”. The second andfifth authors have been partially supported by the funds of the University of Padova, projectCPDA124755 “Multivariate approximation with application to image reconstruction”.c©CMMSE ISBN: 978-84-617-2230-3Page 323 of  1268PU-WSVD10−410−210010210−410−310−210−1100εRMSE10−410−210010210−810−610−410−2100εRMSEFigure 1: Comparison between the global approximant computed with the WSVD basis(dotted lines) and the new method (solid lines), for the C2 Mate´rn kernel (left) and theGaussian kernel (right), for 10−4 ≤ ε ≤ 102 and n = 652 Halton points in the unit square.References[1] R. Cavoretto, A numerical algorithm for multidimensional modeling of scattereddata points, Comput. Appl. Math. 34 (2015) 65–80.[2] R. Cavoretto, A. De Rossi, A meshless interpolation algorithm using a cell-basedsearching procedure, Comput. Math. Appl. 67 (2014) 1024–1038.[3] R. Cavoretto and A. De Rossi, A trivariate interpolation algorithm using a cube-partition searching procedure, submitted (2015)[4] S. De Marchi, G. Santin, A new stable basis for radial basis function interpolation,J. Comput. Appl. Math. 253 (2014) 1–13.[5] S. De Marchi, G. Santin, Fast computation of orthonormal basis for RBF spacesthrough Krylov space methods, to appear in BIT DOI: 10.1007/s10543-014-0537-6(2015).[6] S. De Marchi, R. Schaback, Stability of kernel-based interpolation, Adv. Comput.Math. 32 (2010) 155–161.[7] G. E. Fasshauer, M. J. McCourt, Stable evaluation of Gaussian radial basis func-tion interpolants, SIAM J. Sci. Comput. 34 (2012) A737–A762.c©CMMSE ISBN: 978-84-617-2230-3Page 324 of  1268R. Cavoretto, S. De Marchi, A. De Rossi, E. Perracchione, G. Santin10−410−210010210−610−410−2100102εRMSE10−410−210010210−1010−51001051010εRMSEFigure 2: Comparison between the standard PU approximant (dotted lines) and the newmethod (solid lines), for the C2 Mate´rn kernel (left) and the Gaussian kernel (right), for10−4 ≤ ε ≤ 102 and n = 652 (black line), n = 1292 (light gray), n = 2572 (dark gray)Halton points in the unit square.[8] B. Fornberg, E. Larsson, N. Flyer, Stable computations with Gaussian radialbasis functions, SIAM J. Sci. Comput. 33 (2011) 869–892.[9] B. Fornberg, C. Piret, A stable algorithm for at radial basis functions on a sphere,SIAM J. Sci. Comput. 30 (2007/08) 60–80.[10] B. Fornberg, G. Wright, Stable computation of multiquadrics interpolants for allvalues of the shape parameter, Comput. Math. Appl. 47 (2004) 497–523.[11] S. Mu¨ller and R. Schaback, A Newton basis for kernel spaces, J. Approx. Theory,161 (2009) 645–655.[12] M. Pazouki, R. Schaback, Bases for kernel-based spaces, J. Comput. Appl. Math.236 (2011) 575–588.[13] W. Pogorzelski, Integral Equations and Their Applications. Vol. I, InternationalSeries of Monographs in Pure and Applied Mathematics, vol. 88, Pergamon Press,Oxford, 1966.[14] H. D. Simon and H. Zha, Low-rank matrix approximation using the Lanczos bidiag-onalization process with applications, SIAM J. Sci. Comput. 21 (2000) 2257–2274.c©CMMSE ISBN: 978-84-617-2230-3Page 325 of  1268PU-WSVD[15] H. Wendland, Fast evaluation of radial basis functions: methods based on partitionof unity, Approximation theory, X (St. Louis, MO, 2001), Vanderbilt Univ. Press,Nashville, TN, (2002) 473–483.[16] H. Wendland, Scattered Data Approximation, Cambridge Monogr. Appl. Comput.Math., vol. 17, Cambridge Univ. Press, Cambridge, 2005.c©CMMSE ISBN: 978-84-617-2230-3Page 326 of  1268

RBF approximation of large datasets by partition of unity and local stabilization

We present an algorithm to approximate large dataset by Radial Basis Function
(RBF) techniques. The method couples a fast domain decomposition procedure with a
localized stabilization method. The resulting algorithm can efficiently deal with large
problems and it is robust with respect to the typical instability of kernel methods

RBF approximation of large datasets by partition of unity and local stabilization

Abstract

Similar works

Full text

Available Versions

Institutional Research Information System University of Turin

Archivio istituzionale della ricerca - Università di Padova

Padua Research Archive