# 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[3]]*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(...).

But, sidelining my "thoughts" about this factors, this doesn't seem to be a useful approximation.

• There is no need to use neural networks for this. The functions Fit and LinearModelFit, NonlinearModelFit, etc. would do the job. – Anton Antonov Dec 29 '17 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.) – Anton Antonov Jan 2 '18 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. – pcalc Jan 3 '18 at 9:10
• Cross-posted: community.wolfram.com/groups/-/m/t/1255676 – Michael E2 Jan 4 '18 at 5:11

## 3 Answers

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]


Radial basis functions:

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]


## Similar discussions/answers

• 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). – pcalc Dec 29 '17 at 18:56
• @pcalc Welcome. I listed similar discussions in my answer. – Anton Antonov Dec 30 '17 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. – pcalc Dec 31 '17 at 13:36

## 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[3]] 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[1], 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[1], 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[1]]}],
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! – Anton Antonov Jan 4 '18 at 15:11
• What @AntonAntonov said. – Daniel Lichtblau Jan 5 '18 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! – pcalc Jan 6 '18 at 14:25
• Looks like I can check off 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.) – J. M.'s technical difficulties Mar 5 '18 at 0:05

You could use ElementwiseLayer to calculate an RBF:

net = NetChain[{
LinearLayer[5],
ElementwiseLayer[Exp[-#^2/2] &],
LinearLayer[1]
}, "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}]