# Approximation with radial basis functions

I want to approximate some functions with basisfunctions, which can be easy Fourier-transformed. So I got the idea to approximate my function with Gaussian normal-distribution curves. This led me to the idea of approximating my function with radial basis functions (RBF), such as the Gaussian curves are.

Due to this I found the "NeuralNetworks" Mathematica package, which seems to be exactly what I was looking for. I even found this very helpful Example. Poorly this package doesn't seem to be part of Mathematica 11 anymore. The newer built-in Version doesn't seem to include functions, which do similar work.

Because I am quite new at Mathematica I have no idea if and how I can find another way to solve my problem.

Does anybody have an idea how I can approximate my function (maybe similar to the example?) with Gaussian normal-Distribution curves?

I would be very grateful, if anyone could help me!

Thank you very much!

UPDATE:

I applied @Anton Antonov's Suggestion to a part of my problem.

Version: 11.1.1 for Microsoft Windows (64-bit) (April 18, 2017)

f[x_] := AiryAi[x + AiryAiZero]*HeavisideTheta[x]
Ndata = 200;
x = Table[-3 + 12 N[(i + 0.001)/Ndata], {i, 0, Ndata - 1}];
y = f[x];
data = Transpose[{x, y}];
GaussianRadialFunction[k_, s_, r_] := Exp[-(k (r - s))^2];
bFuncs = Flatten[
Table[GaussianRadialFunction[k, s, r], {k, 0.05, 2, 0.15}, {s, 0, 6, 0.3}]];
lm = LinearModelFit[data, bFuncs, r];
Show[ListPlot[data, PlotStyle -> Pink], Plot[lm[x], {x, -5, 15}],Frame -> True] I also tried to use Gauss-functions, which are shifted from the origin (as you can see in the code above).

The approximation looks pretty good, at least in the range, where data was provided. But beside this regions (where I thougt Gaussian curves would start going to zero) the approximation skyrockets downwards (or upwards, depending on the certain approximation).

Normal[lm]


gives something like:

-9.96495*10^8 + 1.05677*10^9 E^(-4. (-6. + r)^2) + 1.14224*10^10 E^(-3.4225 (-6. + r)^2)+ (...)


The constant first term seems to explain this behaviour. I'm a little bit perplexed, why such "huge" Multipliers occur next to each Exp(...).

• There is no need to use neural networks for this. The functions Fit and LinearModelFit, NonlinearModelFit, etc. would do the job. Dec 29, 2017 at 17:39
• I do not observe what you show in your update: see this screenshot. What version of Mathematica are you using? (Mine is 11.2.0.) Jan 2, 2018 at 18:58
• My version is 11.1.1 for Microsoft Windows (64-bit). What happens at your plot for x>10 ? Does it go to zero? The different behaviour for x<0 is, because I multiplied AiryAi with the Heaviside function (this product is the function I want to Fourier-transform, AiryAi without the Heaviside is easy to tranform). But I think this is not essential for my problem with the skyrocketing approximation, because it also happens for positive x. Jan 3, 2018 at 9:10
• Cross-posted: community.wolfram.com/groups/-/m/t/1255676 Jan 4, 2018 at 5:11

## Introduction

Radial basis function approximation has the form $$f(x) \approx \sum_{j=1}^N c_j \phi(\| x-x_j\|)$$ for some function $\phi(r)$ such as a Gaussian. Solving for the coefficients $c_j$ is a linear problem and might be done either as a least-squares fit or as an interpolation problem. As mentioned on the Community posting, http://community.wolfram.com/groups/-/m/t/1255676, there is a basis for RBF interpolation over a regularly spaced grid that has nice properties: it has the "cardinal basis" property, is a numerically accurate approximation of Gaussian RBF, and, for the OP's purpose, has a known Fourier transform.

I am not an expert on RBF approximation, but it is known that RBF has some issues: "error saturation," in which the interpolation error does not vanish as the grid spacing approaches zero; and "boundary layer" at the end points of a truncated, finite interval. Some appreciation of error should be considered in the choice of a method and its parameters.

We will give code for a straightforward implementation of RBF interpolation over a regularly spaced grid over a truncated, finite interval. We next apply it to the OP's Airy function example. To give a nice interface to the Fourier transform of an RBF interpolation, we give a more sophisticated implementation of the interpolation, as well as its transform, and apply them to the OP's example. Finally, we illustrate the issues with error in RBF interpolation with regard to the OP's example. It will be seen that the RBF approximation is very good overall and worst in the neighborhood of the singularity of the HeavisideTheta factor at x == 0.

## A cardinal basis for Gaussian radial basis function interpolation

A "cardinal basis" for a grid $x_j$ is a set of functions $C_j(x)$ such that $$C_j(x_j)=1,\quad\text{and}\quad C_j(x_k)=0\ \text{for all k\ne j.}$$ The Lagrange interpolation basis is a cardinal basis for a finite grid. The particular advantage of such a basis is that the interpolation formula does not require solving an equation: $$f(x) \approx \sum_j f(x_j)\,C_j(x)\,.\tag{1}$$

Maz'ya & Schmidt, Approximate Approximations, (2007) and Boyd & Wang, "An analytic approximation to the cardinal functions of Gaussian radial basis functions on an infinite lattice", (2009) independently derived the basis for an infinite grid with spacing $h = x_j - x_{j-1}$ and shape parameter $\alpha$: $$C_j(x)=C(x-x_j),\quad\text{where}\quad C(x)= {\alpha ^2\over\pi}\, {\sin \left({\pi x}/{h}\right) \over {\sinh}\left({\alpha^2 x}/{h}\right)}$$ Since for this grid, the formula in (1) is an infinite sum, we will have to truncate it and sum over a finite interval.

A basic implementation for this idea is the following (note that the sum (1) is effected with Dot):

(* cBasis[x_, α_, h_:1] := α^2/Pi Sin[Pi x/h]/Sinh[α^2 x/h]; (* risks divide-by-zero *) *)
cBasis[x_, α_, h_: 1] := Sinc[Pi x/h]/Sinc[I α^2 x/h];    (* uses complex numbers *)

rbfInterpolation[f_, {x_, x1_, x2_}, α_, h_] := (* interpolates f over x1 <= x <= x2 *)
With[{fvals = f /. x -> # & /@ Range[x1, x2, h]},
Function[x, fvals.cBasis[x - Range[x1, x2, h], α, h]]
];


The main problem is that cBasis[0, α, h] needs to evaluate to 1, and the first, commented code above yields Indeterminate (from 0/0). The second avoids dividing by zero but does its computations with complex numbers (with a zero imaginary part). One could define cBasis[x_ /; x == 0, α_, h_: 1] := 1 explicitly and make it Listable, but this is slower. In any case, I wanted to start with simpler code. A slightly improved (especially if x == 0 is rare), but somewhat advanced, version will be given below with a more convenient, if more complicated, rbfInterpolation. These improvements will make it easier to handle the Fourier transform of an RBF interpolation.

## Application to OP's Airy function

ff[x_] := AiryAi[x + AiryAiZero] UnitStep[x]; (* HeavsideTheta is undefined at 0 *)
testIF = rbfInterpolation[ff[x], {x, 0, 12}, 1/2, 1/8];


Here's a plot showing the OP's function (blue) covered by the interpolation function (gold), together with the logarithm base-10 of the error (green) below, which we plot beyond the interpolation range:

Plot[{20 ff[x], 20 testIF[x], ff[x] - testIF[x] // RealExponent},
{x, -2, 15},
Frame -> True, FrameTicks -> {
{Automatic, ChartingScaledTicks[{20 # &, #/20 &}]},
{Automatic, Automatic}},
FrameLabel -> {{"Log error", HoldForm@f[x]}, {HoldForm[x], None}}] An understanding of the error in RBF interpolation can help choosing α and h and help explain why the error decreases between the end points. See Error analysis below.

## Improved interpolation code and Fourier transform

One concern is speed, only because this method already appears somewhat slow to me. When possible, the code should preserve packed arrays, and use arrays of real numbers for real values. Also, the Fourier transform of an RBF interpolation is simply $$\hat{f}(t) \approx \sum_j f(x_j) \, \hat{C}_j(t)$$ where $\hat{C}_j(t) = e^{i x_j t} \hat{C}(t)$ and $\hat{C}$ is the Fourier transform of $C$, is given by $$\hat{C}(t) = -\frac{h \, e^{i x_j t} \left(\tanh \left(\frac{\pi (h t-\pi )}{2 \alpha^2}\right)-\tanh \left(\frac{\pi (h t+\pi )}{2 \alpha ^2}\right)\right)}{2 \sqrt{2 \pi }}\,.$$ The only data needed are the grid $x_j$ and the function values $f(x_j)$ on the grid. We can store these in a simple data structure, similar to InterpolatingFunction, instead of in a Function as above.

To fix the divide by zero problem, we can adapt @Carl Woll's answer, which has the property of preserving packed arrays, when the arguments do not lead to 0/0. Even when they do, the unpacking is faster than making cBasis have the attribute Listable.

ClearAll[cBasis];   (* cardinal basis *)

(* for symbolic analysis *)
cBasis /:  Normal[cBasis[x_, α_, h_: 1]] := α^2/Pi Sin[Pi x/h]/Sinh[α^2 x/h];

division[a_, b_] := Quiet@Block[{Indeterminate = 1},  Divide[a, b]];
cBasis[x_?NumericQ | x_?(VectorQ[#, NumericQ] &), α_, h_: 1] :=
division[α^2/Pi Sin[Pi x/h], Sinh[α^2 x/h]];

ClearAll[rbfIF, rbfInterpolation];  (* RBF interpolation *)
rbfInterpolation[f_, {x_, x1_, x2_}, α_, h_] :=
Module[{xdata, fdata},
xdata = Range[x1, x2, h];
fdata = DeveloperToPackedArray[f /. x -> # & /@ xdata];
rbfIF[xdata, fdata, {α, h}]
];
rbfIF[xdata_, fdata_, {α_, h_}][x_] := fdata.cBasis[x - xdata, α, h];
rbfIF /: MakeBoxes[if_rbfIF, form_] := BoxFormArrangeSummaryBox[ (* optional formatting *)
rbfIF, if, None, {{MinMax@First@if}}, {}, form, "Interpretable" -> Automatic];

ClearAll[ft, rbfFT];  (* Fourier transform of the cardinal basis/RBF interpolation *)
Block[{k, x0, α, h, x},
ft[k_, x0_, α_, h_: 1] = (* pre-calculate basic Fourier transform (takes a few secs.) *)
Simplify[
Exp[I k x0] FourierTransform[Normal[cBasis[x, α, h]], x, k,
Assumptions -> α > 0 && h > 0]]
];
rbfFT[rbfIF[args___]] := rbfFT[args];  (* for the transform of a rbfIF, change the head *)
rbfFT[xdata_, fdata_, {α_, h_}][k_] := fdata.ft[k, xdata, α, h];
rbfFT /: MakeBoxes[ft_rbfFT, form_] := BoxFormArrangeSummaryBox[ (* optional formatting *)
rbfFT, ft, None,  {{MinMax@First@ft}}, {}, form, "Interpretable" -> Automatic];


The expressions for cBasis and its Fourier transform are simple elementary expressions. For machine-precision input, they will be vectorized. Note that the symbolic form of cBasis is used to compute the expression for the Fourier transform. One might use different user-interfaces for the symbolic cBasis.

## Fourier transform of the OP's Airy function

testIF = rbfInterpolation[ff[x], {x, 0, 15}, 1/2, 1/8]; (* re-evaluate with new defs *)
testFT = rbfFT@testIF;  (* Fourier transform *)


For comparison here is a numerical integration for the transform (Integrate/FourierTransform failed):

ClearAll[nFT];
Module[{t},
li = NIntegrateLevinIntegrandReduce[ff[x] Exp[I t x] /. UnitStep[_] :> 1, x];
With[{lirules = Sequence @@ (DeleteCases[li["Rules"], _["Variables", _]] /.
HoldPattern["DifferentialMatrices" -> dm_] :> "DifferentialMatrix" -> First[dm])},
nFT[t0_?NumericQ] := Block[{t = t0},
1/Sqrt[2 Pi] NIntegrate[
ff[x] Exp[I t x] /. {UnitStep[_] :> 1} // Evaluate,
{x, 0, Infinity},
Method -> {"LevinRule", lirules}]
]]];


The magnitude and argument of the transform of ff and the interpolation overlay each other:

Plot[{nFT[t], testFT[t]} // Abs // Evaluate, {t, -10, 10}] (* Compare Abs *)
Plot[{nFT[t], testFT[t]} // Arg // Evaluate, {t, -10, 10}, (* Compare Arg *)
MaxRecursion -> 3, Exclusions -> None]  ## Error analysis and the boundary layer phenomenon

There are at least three parameters affecting error, $\alpha$, $h$, and the "round-off plateau," which depends on working precision and is ordinarily around $\epsilon\,\|f\|_\infty$ or a little higher, where $\epsilon$ is the floating-point precision limit (e.g., $MachineEpsilon for machine precision). Connected with$\alpha$and$h$, is the "boundary layer" phenomenon (see below), which affects the error near the end points of the interval. Away from the boundary, the error seems to depend principally on$\alpha$and$h$. According to Boyd, "Error saturation in Gaussian radial basis functions on a finite interval" (2010), the error tends to get worse if$\alpha$gets too small or too large. Just by testing$\alpha = 0.4, 0.5, 0.6$etc., I found$\alpha = 0.5$to be good on the OP's example, although it varies slightly with$h$. The boundary layer effect decreases as ~Exp[-α^2 x/h] according to Boyd, op. cit. (2010). This can be seen for various$h$in the dashed red envelopes in the figure below. The singular point at$x=0$impose by the UnitStep/HeavisideTheta factor means that$x = 0$is a boundary point no matter how far we extend the interpolation. It seems little can be done about the error near$x=0$. On the other hand one can see that at$x \approx 9$, the end point of the OP's interpolation example, the error is quite good for$h = 1/32, 1/64$, because the interpolation range was extended to$x = 12$, far enough out that the boundary layer effect is negligible compared to round-off. (The round-off plateau is around machine precision until the function becomes small and$x$far enough away from the boundary-layer effect of the large peak near$x=4.5$; see the lower grid line at$\log_{10}\epsilon\approx-15.95$in the figure below.) testIF8 = rbfInterpolation[ff[x], {x, 0, 12}, 0.5, 1./8]; testIF16 = rbfInterpolation[ff[x], {x, 0, 12}, 0.5, 1./16]; testIF32 = rbfInterpolation[ff[x], {x, 0, 12}, 0.5, 1./32]; testIF64 = rbfInterpolation[ff[x], {x, 0, 12}, 0.5, 1./64]; errplot = Plot[ (* errors of the approximations *) ff[x] - {testIF8[x], testIF16[x], testIF32[x], testIF64[x]} // RealExponent // Evaluate, {x, -2, 14}, Frame -> True, FrameLabel -> {HoldForm[x], "Log error"}, PlotStyle -> Opacity[0.7], PlotLegends -> Thread[HoldForm[h] == 1/{8, 16, 32, 64}], PlotPoints -> 100] boundarylayer = Show[ (* boundary layer envelopes *) Plot[0.001 Range[5,2,-1] * (* constant factor found by eye *) Exp[-0.5^2 (x)/(1./{8, 16, 32, 64})] // RealExponent // Evaluate, {x, -2, 14}, PlotStyle -> Directive[AbsoluteThickness, Dashing[{0.02, 0.02}], Red]], Plot[3*^-7 Exp[-0.5^2 (12 - x)/(1./{8, 16, 32, 64})] // RealExponent // Evaluate, {x, -2, 14}, PlotStyle -> Directive[AbsoluteThickness, Dashing[{0.02, 0.02}], Red]] ]; Show[errplot, boundarylayer, GridLines -> {None, {-6.6, -9.5, RealExponent[$MachineEpsilon]}},
Frame -> True,
PlotLabel -> Row[{"Boundary layer effects of ", HoldForm[f[x]], " at ",
HoldForm[x == {0, 12}], " (", HoldForm[α == 0.5], ", ",
HoldForm[h == {1/8, 1/16, 1/32, 1/64}], ")"}]] Here is another view of the boundary layer effect at $x=9$ and what the inclusion of the basis functions beyond the boundary add to the interpolation. The basis functions, which in the interpolation are multiplied by the function values, contribute a relatively large amount to the interpolation near $x=9$. We show it for $h = 1/8$ instead of one of the more accurate grid spacings, because it is easier to see what is going on. One can compare the decreasing contribution below with the decreasing error above (see the top two grid lines).

With[{x1 = 9, x2 = 12, α0 = 0.5, h0 = 1./8},
Plot[
Table[cBasis[x - x0, α0, h0] // RealExponent, {x0, x1, x2, h0}] //Evaluate,
{x, x1 - 2, x2 + 2},
PlotStyle -> ReplacePart[
Table[Directive[AbsoluteThickness[1.], Opacity[0.5]], {x0, x1, x2, h0}],
{-1 -> Directive[Black, Opacity]}],
GridLines -> {{x1}, Automatic}, GridLinesStyle -> Lighter@Blue,
Frame -> True,
FrameLabel -> {HoldForm[x], HoldForm[Log10@HoldForm[C][x - x0]]},
PlotLabel ->
Row[{"Boundary layer effects at ", HoldForm[x == x1], " (",
HoldForm[α == α0], ", ", HoldForm[h == h0], ")"}]]
] The following show the magnitude of the error in the Fourier transform of ff and the approximation as well as the magnitude of the phase error, for the four approximations.

Plot[nFT[t] -
Through[(rbfFT @@@ {testIF8, testIF16, testIF32, testIF64})[t]] //
RealExponent // Evaluate,
{t, -10, 10},
MaxRecursion -> 3, Exclusions -> None,
PlotLegends -> Thread[HoldForm[h] == 1/{8, 16, 32, 64}],
PlotLabel -> "Log error"]

Plot[Arg@nFT[t] -
Arg@Through[(rbfFT @@@ {testIF8, testIF16, testIF32, testIF64})[t]] //
RealExponent // Evaluate,
{t, -10, 10},
MaxRecursion -> 3, Exclusions -> None,
PlotLegends -> Thread[HoldForm[h] == 1/{8, 16, 32, 64}],
PlotLabel -> "Log error in Arg"]  • Impressive post! Jan 4, 2018 at 15:11
• What @AntonAntonov said. Jan 5, 2018 at 0:23
• I'm speechless - this post is incredible! Thanks a lot! I'll work through it and will get back after. Thank you very much! Jan 6, 2018 at 14:25
• Looks like I can check the implementation of a Gaussian RBF off my list. ;) (IIRC, the ill-conditioning for $\alpha$ can be ameliorated by using a complex interpolation, but I haven't had time to explore that avenue.) Mar 5, 2018 at 0:05

Here is an answer with LinearModelFit using the data in the example referenced in the question.

Data:

Ndata = 20;
x = Table[10 N[i/Ndata], {i, 0, Ndata - 1}];
y = Sin[0.1 x^2];
data = Transpose[{x, y}]
ListPlot[data] GaussianRadialFunction[k_, r_] := Exp[-(k r)^2];
bFuncs = Table[GaussianRadialFunction[k, r], {k, 0.01, 1, 0.04}]
Plot[bFuncs, {r, 0, 10}, PlotRange -> All] Linear model fit with the basis functions:

lm = LinearModelFit[data, bFuncs, r];
Show[ListPlot[data, PlotStyle -> Pink], Plot[lm[x], {x, 0, 10}], Frame -> True] • Thank you very much- this seems to be exactly what I was looking for! I'm really surprised, that even unshifted Gauss-curves give such a great approximation, I thougt, that I would have to vary k and µ,the shift from the origin, as well. Is there a way to find out the coefficients of the approximation? CoefficientList only works for polynomials, as far as I know (and tried). Dec 29, 2017 at 18:56
• @pcalc Welcome. I listed similar discussions in my answer. Dec 30, 2017 at 17:55
• I tried to follow these posts, but (probably because of my lack of Mathematica-knowledge) I coudln't find a solution of my problem on my own. Your suggested solution seems to work, but there is still one issue, I can't clear. So I edited my post above and added my proceeding. I would be glad, if you still have some hints for me. Dec 31, 2017 at 13:36

You could use ElementwiseLayer to calculate an RBF:

net = NetChain[{
LinearLayer,
ElementwiseLayer[Exp[-#^2/2] &],
LinearLayer
}, "Input" -> 2]

trained =
NetTrain[net, {{0, 0} -> {0}, {1, 0} -> {1}, {0, 1} -> {1}, {1,
1} -> {0}}]

ContourPlot[trained[{x, y}], {x, 0, 1}, {y, 0, 1}] • There is a problem about this method, do you know how to solve for different sigma and miu for the different neurons. Dec 17, 2022 at 23:27