Radial Basis Function (RBF) methods are important tools for scattered data interpolation and for the solution of Partial Differential Equations in complexly shaped domains. The most straight forward approach used to evaluate the methods involves solving a linear system which is typically poorly conditioned. The Matlab Radial Basis Function toolbox features a regularization method for the ill-conditioned system, extended precision floating point arithmetic, and symmetry exploitation for the purpose of reducing flop counts of the associated numerical linear algebra algorithms.

Radial Basis Function (RBF) methods are important tools for scattered data interpolation and for
the solution of PDEs in complexly shaped domains. The most straight forward approach that is used to
evaluate the method while incorporating the “standard basis functions” involves solving
a linear system which is typically poorly conditioned. Two variations of a method, dubbed the RBF-QR
approach, use a different basis that spans the same space as the standard basis but which in some
cases results in a better conditioned linear system. Software packages which implement the two
RBF-QR approaches are freely available ([

Extended precision floating point arithmetic can be used to accurately evaluate the
ill-conditioned problem in the standard basis. This approach has been used and implemented in
several different software environments that include: Mathematica [

In this work a Matlab [

The

RBF interpolation uses a set of ^{d}

is an infinitely differentiable (compactly supported and global RBFs without a shape parameter
and with less smoothness exist but are not considered here) function of one variable

where

and the inverse quadratic (IQ) RBF

are representative members of the class of global, infinitely differently RBFs containing a shape
parameter that interpolate with exponential accuracy. The two RBFs and their various derivative are
defined respectively in the classes

Many other RBFs exist and may be added by the user to the toolbox by extending rbfx and using gax and iqx as examples.

The entire list of functions associated with an object is available as follows:

First the methods that the class must define are listed and then the static methods inherited from the superclass are listed.

The RBF expansion coefficients are determined by enforcing the interpolation conditions

which result in a

The matrix

is called the system matrix. The solve method which is a static method of the rbfx class can be used to find the expansion coefficients. Static methods can be called in two ways. Either through the class

or from instances of the class

Matlab is optimized for operations involving matrices and vectors. The process of revising loop-based, scalar-oriented code to use Matlab matrix and vector operations is called vectorization. At every opportunity, functions have been vectorized so that they execute as efficiently as possible.

The functions that define RBFs and their derivative matrices take distance matrices as their
arguments. Taking for example 2d where

respectively contain the signed distance between the

stores the distance between centers

and the RBF system matrix is constructed as

The evaluation of the interpolant (2) at _{j}

Continuing the 2d example, consider evaluation points _{1}, _{2}). The distance matrices containing
the distances between the centers and evaluation points

respectively contain the signed distances between the

With the

and the evaluation matrix is constructed via

The interpolant is then evaluated as

Derivatives are approximated by differentiating the RBF interpolant as

where 𝒟 is a linear differential operator. The operator 𝒟 may be a single
differential operator or a linear differential operator such as the Laplacian. Evaluating (13) at
the centers _{𝒟}

That is, 𝒟_{𝒟}

Alternatively, derivatives can be approximated by mutiplying the grid function values _{𝒟}^{–1} since

This is accomplished as

The shape parameter _{k}

and well as the linearly varying shape parameter

Reference [

All three variable shape strategies are available in the MRBFT via

Both equations (6) for the expansion coefficients and (15) for the differentiation matrix assume
that the system matrix is invertible. Both the GA and IQ system matrices are symmetric positive
definite (SPD) if a constant shape parameter is used and therefore they are invertible. While it is
invertible, the system matrix is typically very poorly conditioned and it may cease to be
numerically SPD and the standard algorithm to factorize a SPD matrix may fail. The eigenvalues of
_{min}_{1}
≤ 𝜆_{2} ≤ … ≤ 𝜆_{N}_{max}_{min}_{max}^{18}) in double
precision and 𝒪(10^{36}) in quadruple precision. A regularization technique to
mitigate the effects of the ill-conditioning is discussed in the next section.

Interpolation of a smooth function using both double (F64) and quadruple (F128) precision. The
script

Recent monographs [

Reference [

the regularized system

where

For small ^{–1} is close to ^{–1} and
MDI simply replaces B with (

quantifies how close that (^{–1} and
^{–1} are [_{M}_{M}

All MRBFT functions that involve the factorization of a SPD matrix take two optional arguments: a
regularization parameter

Figure ^{16}) and smaller. When the condition number
of the unregularized method reaches beyond that threshold the regularization keeps the condition
number in the approximately 𝒪(10^{16}) range and the regularized solution is
approximately two decimal places more accurate in this example.

Interpolation with (green solid) and without (blue dashed) regularization. The script

The MRBFT uses the Multiprecision Computing Toolbox for Matlab (MCT) [

The following code that calculates the RBF expansion coefficients gives an example of how to change double precision computations to extended precision:

The only coding difference in the above code between double and extended precision is on lines 2
and 3. Once the number of digits is specified and

The second major class of the MRBFT is

RBF methods place no restriction on the location of centers. However, randomly locating centers
is unlikely to lead to good results. Theory dictates [

Quasirandom sequences [

provide the option to cluster the centers more densely around the boundaries of the domains.
Figure

1000 centers on the unit circle. Left: Hammersley points. Right: clustered Hamersley points.

The third major class of the MRBFT is

A matrix

Center distributions that result in centrosymmetric distance matrices. Left: 2000 centers based
on clustered Hammersley points that have been extended centro symmetrically about the

where _{11}, _{21}, and

are orthogonally similar which allows for efficient eigenvalue calculation. Similar
decompositions are possible for odd

The MRBFT functions for centrosymmetry include the following. Centrosymmetric center distributions in circular domains are produced by

and in domains with either

uses half of the centers given as an argument to return a centrosymmetric center distribution.

Only half of the matrix needs to be formed and stored and then subsequently operated on in centrosymmetric linear algebra routines. The half-sized distance matrices are formed via calling the distanceMatrix functions with the arguments modified as follows:

The half sized distance matrices can be used to form half-sized system and derivative evaluation matrices as

Matrix condition numbers in the 2-norm are calculated via

with a factor of four savings in the dominant term in the flop count compared to the standard algorithm. The RBF expansion coefficients for interpolation problems can be found from the half-sized system matrix as

Two quarter-sized matrices are factorized by a Cholesky factorization and the dominant term in
the flop count is reduced by a factor of four from

with a factor of four savings in computational effort. Centrosymmetric matrix-vector multiplication is asymptotically faster by a factor of two over the standard algorithm and can be accomplished via

Other miscellaneous functions concerning centrosymmetry are

which test matrices for symmetry and

which constructs a full-size centrosymmetric matrix from a half-sized matrix.

The MRBFT distribution includes three directories that contains scripts that provide examples of
using the functions in the toolbox. The

The following example scripts come with the toolbox:

_{xx}_{yy}^{2}
sin(

The problem is solved two ways: 1) standard algorithms, 2) centrosymmetric algorithms. With N = 5000 the accuracy of the two approaches is the same but the centrosymmetric approach is approximately five times faster and requires only half the storage compared to the standard approach.

_{t}_{xx}_{yy}_{2}(1 –

The problem is solved two ways: 1) standard algorithms, 2) centrosymmetric algorithms. With

The following test scripts come with the toolbox:

^{16}).

The following benchmark scripts come with the toolbox:

Additional results demonstrating the benefits of using extending precision with RBF methods and
the efficiency gains from exploiting centrosymmetry can be found in references [

The MRBFT is a freely available collection of Matlab files and scripts that implement RBF methods
for scattered data interpolation and for the numerical solution of PDEs. The toolbox has been
developed over a period of several years while it has been used in the author’s research. The
class

The author uses the MRBFT in his own research and as a result the toolbox is constantly being
improved and new features are being added. Future improvements to the MRBFT include: developing a
class of methods for working with the local RBF method, implementing Mercer’s method for the
GA RBF [

A summary of the implementation of the MRBFT is as follows. The toolbox is implemented in Matlab which is widely used in Mathematics and other scientific disciplines. An object oriented approach is used to organize the functionality of the toolbox.

The MRBFT has been developed over a period of several years as it has been used in the
author’s research. It was used extensively in the preparation of references [

The MRBFT is programmed in Matlab which is available on Windows, OS X, and Linux.

The MRBFT was developed and tested using Matlab versions 2015b through 2016b [

The MRBFT uses the Multiprecision Computing Toolbox for Matlab (MCT) [

The author is the only contributor to the software. However, the contributions of others are welcome.

The potential for others to use the MRBFT is substantial. RBF methods have become very popular in applications. The MRBFT may potentially find use in industrial applications. Research in RBF methods is very active. The MRBF is a potentially useful tool for faculty research, graduate student research, and undergraduate research projects. The object oriented design of the toolbox makes it extremely easy for it to be extended by subclassing rbfx to define new RBFs. The MRBFT can be used to foster the idea of reproducible research in the area of RBF methods. The MRBFT provides core functionality that other research can be built on and for which the code is freely available.

The author has no competing interests to declare.