Radial Basis Approximation

[Rodolphe-Vivant]

Project Radial Basis and Low-Rank

Overview

The main objective of my $1^{st}$ year of Master degree project, under the supervision of Dr. Emmanuel FRANCK, was to implement various Radial kernels in Scimba in order to generate mehsless approximation methods to solve numerically PDE's :

where $\mathcal{L}$ and $\mathcal{L}_{\partial\Omega}$ a linear differential operator, $f$ is a source term and $g$ is a boundary condition.

Radial Functions

Radial basis functions (RBFs) are a class of functions that depend only on the distance between a point $x$ and a center point $\mathbf{\xi}$.

where $\mathbf{\xi}_j$, $j=1,..,N$ are the centers of the RBFs and $\epsilon$ is a parameter that controls the width of the RBF.

We can define the RBF approximant of a certain function $f$ as:

where $\alpha_j$ are the coefficients of the basis functions centered at $\mathbf{\xi}_j$ for $j=1,..,N$.

In other words, $s(x,\epsilon)\in span(K(.,\xi_j)_{j\in[1,N]})$.

In our case, we will need to guarantee a certain regularity of the approximant. The function tested, or the solution of the PDEs we want to approximate, must be in a Sobolev spaces $H^{m}(\mathbb{R}^d)$, where $d$ is the dimension of the domain $\Omega$ we are working on.To do so, we will have to ensure that our kernel are positive (semi) definite.
The choice of the kernel is crucial for the regularity of the approximant. We will test that property trough our examples.

Reformulation of the initial problem

After insterting the RBF approximant into the PDE, we obtain:

where the $x_i\in\mathcal{I}$ are the collocation points inside $\Omega$ and the $x_i\in\mathcal{B}$ are the collocation points on the border $\partial\Omega$.

The system of equations can be written in a matrix form:

where

Finding the approximant

To solve the system of equations, we can use a linear or a non-linear approach.
To find the approximant, we train our model with the following set of parameters $\theta$ :

• Case Linear : $\theta=\lbrace\alpha_j\rbrace$
  • There is only one step of linear regression to find the coefficients $\alpha_j$ of the untrained RBFs.
• Case nonLinear : $\theta=\lbrace\alpha_j,\xi_j,\epsilon_j\rbrace$,
  • First we train the parameters of the RBFs, which are the centers $\xi_j$ and the widths $\epsilon_j$ of the kernels,
  • Then we optimize the coefficients $\alpha_j$ of the approximant through a linear regression step.

Results for $1D$, $2D$ and for the solving of a Poisson equation

I was able to implement :

• Isotropic kernels:
  • Gaussian
  • Powered exponential
  • (Inverse) multiquadric
  • Matern (nu=1/2, 3/2, 5/2)
• Anisotropic kernels:
  • Anisotropic Gaussian

The figures of the results include the linear projection and the non-linear projection of the solution as presented before.

Isotropic Kernel

Gaussian Kernel

• Powered exponential kernel :

$$ \kappa(\epsilon r) = e^{-(\epsilon r)^\beta} $$

• Gaussian kernel : $\beta=2$

• Allows precise approximation of functions $C^\infty(\mathbb{R}^d)$

• Approximation of the $\mathcal{C}^\infty(\mathbb{R}^2$ function $f(x,y)=\sin(\cos(x)\times \sin(y))/(2\times 0.4^2)$ with Gaussian Kernels' approximated space.

Gaussian RBF's approximation of f

Exponential Kernel

• Powered exponential kernel with $\beta=1$ :

$$ \kappa(\epsilon r) = e^{-(\epsilon r)^\beta} $$

• Allows precise approximation of functions $C^0(\mathbb{R}^d)$

Figure: Exponential RBF's approximation of f

• Less precise than the Gaussian Kernel for this $C^\infty(\mathbb{R}^2)$ function.

(Inverse) Multiquadratic Kernel

• (Inverse) Multiquadric kernel : $\kappa(\epsilon r) = (1+\epsilon^2 r^2)^\beta$, $\beta \in \mathbb{R} \setminus \mathbb{N}_0$
• $\beta < 0$ : Inverse MultiQuadric (IMQ) kernel
• $\beta > 0$ : MultiQuadratic (MQ) kernel

Figure: Multiquadric RBF's approximation of f, $\beta = -2$

Good result for the projection of this $C^\infty(\mathbb{R}^2)$ function. It was expected as the (I)MQ is a simplistic form of the Gaussian Kernel. This why we get comparable result with a lesser precision.

Matern Kernel

• Matern kernel:

$$ \kappa(x_i, x_j) = \frac{1}{\Gamma(\nu),2^{\nu-1}} \left( \frac{\sqrt{2\nu}}{l} d(x_i, x_j) \right)^\nu K_\nu \left( \frac{\sqrt{2\nu}}{l} d(x_i, x_j) \right) $$

where $d(\cdot,\cdot)$ is the Euclidean distance, $\Gamma$ is the gamma function, and $K_\nu$ is the modified Bessel function of the second kind of order $\nu$.

• Generalization of the Radial Basis Function (RBF) kernel:
  • Used in statistics and approximation theory.
  • $\nu = \frac{1}{2}$ : Exponential kernel
  • $\nu \to \infty$ : Gaussian kernel
  • $\nu = \frac{3}{2}$ : Good approximation for once differentiable functions
  • $\nu = \frac{5}{2}$ : Good approximation for twice differentiable functions
  • Generates classical Sobolev spaces $H^{\nu + \frac{d}{2}}(\mathbb{R}^d)$

For $\mathcal{C}^\infty$ function $f(x) = \cos(2\pi x)$ Figure: Gaussian, Matern ($\nu = 0.5$), and Exponential RBF's approximation of $f$

For $\mathcal{C}^0$ function $f(x) = 1 - |x|$ Figure: Gaussian, Matern ($\nu = 0.5$), and Exponential RBF's approximation of $f$

As expected for $\mathcal{C}^\infty$, the Gaussian Kernel works better whereas for the $\mathcal{C}^0$, it's the Matern Kernel and the Exponential that get better results.

Anisotropic Gaussian Kernel

The kernel is given by:

$$ K(\mathbf{x}, \mathbf{z}) = e^{-(\mathbf{x} - \mathbf{z})^T E (\mathbf{x} - \mathbf{z})} $$

where $E$ is a symmetric positive definite matrix. If $E = \epsilon^2 I_d$, we obtain the isotropic Gaussian kernel.

• Test on an anisotropic rotated gaussian function

Anisotropic Gaussian RBF Result

The projection is good.

• Test on our $\mathcal{C}^\infty$ function

Anisotropic Gaussian RBF Result

Again, we obtain a very good projection of our function. The results are logically close to the Gaussian Kernel's results.

Result for the Poisson equation solving

• Solving the following Poisson equation in a 2D domain $\Omega = [-1,1]^2$ with the RBFs:

$$ \begin{cases} -\Delta u(x,y) = f(x,y), & (x,y) \in \Omega \\ u(x,y) = g(x,y), & (x,y) \in \partial \Omega \end{cases} $$

where $f(x,y) = -4$ and $g(x,y) = x^2 + y^2$.

Gaussian RBF's approximation

Exponential RBF's

Multiquadratic RBF's ($\beta=-2$)

Matern RBF's ($\nu=0.5$) approximation

Matern RBF's ($\nu=2.5$) approximation

The results are quite precise for the Gaussian, the IMQ, and the Matern kernel ($\nu=2.5$) as theorically expected. Indeed, those kernels allow precise approximation for functions that are at least twice differentiable as the Strong Form solving of the PDE impones.
The Exponential and the Matern ($\nu=0.5$) can't approach the solution as precisely as the others (again as theoritically expected).

Conclusion

In the future, I believe it would be interesting to start by implementing a more complete Matern Kernel with other orders than $\nu=0.5$, $\nu=1.5$ and $\nu=2.5$.The current Mattern Kernel is imported from gpytorch and their version only allows those three values.
Then, implementing other Radial Kernels like the Oscillatory kernel, learned kernels and maybe the most important one, the Compactly supported kernels. Indeed, the current Kernels are not compactly supported and therefore, the Matrices generated are often dense. Compactly supported Kernels would allow the use of sparce matrices and so enable far better performances.
And last but not least, testing the performance of the RBFs with different types of PDEs and ODEs with different boundary conditions.