Radial Basis Functions - Executable Engineering: Interactive Numerical Methods in Python

Radial basis functions are an n-dimensional interpolation technique that doesn’t rely on polynomials. Rather, we define a radial basis function, called a kernel, applied to each data point:

\varphi_i(||x-x_i||)

(1)

Commonly, we say $\varphi_i(x=x_i)\equiv 1$ .

The kernel only depends on the Euclidian distance between the associated data point, $x_i$ and the evaluation point $x$ (and are therefore radial).

The interpolation function $y(x)$ is the weighted sum of the $N$ kernels:

y(x) = \sum_{j=1}^N \omega_j \varphi(||x-x_j||)

(2)

To determine the weights $\omega_j$ , we evaluate the function at each data point $x_i$ :

y(x_i) = \sum_{j=1}^N \omega_j \varphi(||x_i-x_j||) = y_i

(3)

Applying this to all $N$ data points generates a linear system to solve for $\omega_j$ : $$

\begin{bmatrix} \varphi(\|x_1 - x_1\|) & \varphi(\|x_1 - x_2\|) & \cdots & \varphi(\|x_1 - x_N\|) \\ \varphi(\|x_2 - x_1\|) & \varphi(\|x_2 - x_2\|) & \cdots & \varphi(\|x_2 - x_N\|) \\ \vdots & \vdots & \ddots & \vdots \\ \varphi(\|x_N - x_1\|) & \varphi(\|x_N - x_2\|) & \cdots & \varphi(\|x_N - x_N\|) \end{bmatrix}

(4)

\omega_1 \ \omega_2 \ \vdots \ \omega_N \end{bmatrix} = \begin{bmatrix} y_1 \ y_2 \ \vdots \ y_N \end{bmatrix} $$

which we know how to solve!

RBF Kernels¶

Kernels are defined with $r = \| x-x_i\|$ and a tuning parameter $\epsilon$ . Some common simple kernels are:

Kernel	Formula
Gaussian	$e^{-\epsilon^2 r^2}$
Inverse quadratic	$\frac{1}{1+[\epsilon r ]^2}$
Inverse multiquadric	$\frac{1}{\sqrt{1+[\epsilon r ]^2}}$

Determination of optimal $\epsilon$ is a nuanced question, but a good rule of thumb is to use the average distance between samples.

$\epsilon = \frac{1}{avg \|x_i-x_j\|}$

We will see that its choice has a large effect on the condition number of the matrix, and therefore the numerical robustness!

Let’s see the kernels:

import numpy as np
import plotly.graph_objects as go

try:
    import executable_engineering as exe
except ImportError:
    %pip install -q executable_engineering
    import executable_engineering as exe

# Define the radial basis functions
def gaussian(r, epsilon):
  return np.exp(-(epsilon * r)**2)

def inverse_quadratic(r, epsilon):
  return 1 / (1 + (epsilon * r)**2)

def inverse_multiquadric(r, epsilon):
  return 1 / np.sqrt(1 + (epsilon * r)**2)

epsilon = 1

# Create a range of r values
r_values = np.linspace(0, 10, 100)

# Calculate the function values for each kernel
gaussian_values = gaussian(r_values, epsilon)
inverse_quadratic_values = inverse_quadratic(r_values, epsilon)
inverse_multiquadric_values = inverse_multiquadric(r_values, epsilon)

# Plot the radial basis functions
fig = go.Figure()
fig.add_trace(go.Scatter(x=r_values, y=gaussian_values, name='Gaussian'))
fig.add_trace(go.Scatter(x=r_values, y=inverse_quadratic_values, name='Inverse Quadratic'))
fig.add_trace(go.Scatter(x=r_values, y=inverse_multiquadric_values, name='Inverse Multiquadric'))
fig.update_layout(title='Radial Basis Functions (epsilon = 1)', xaxis_title='r', yaxis_title='φ(r)')
fig.show()

1D RBF Interpolation¶

For these kernels, the influence of a data point diminishes the further from them you get.

In general, $\varphi_i(r=0)$ is not necessarily 1, and $\varphi(r \rightarrow \infty) \ne 0$ , but this requires one more key factor to implement robustly.

Example - Our Toy problem from last lecture (Gaussian sampled at 10 points, equally spaced)¶

# Define the function
def f(x):
  return np.exp(-(x/2)**2)

def gaussian(r, epsilon):
  return np.exp(-(epsilon * r)**2)

# Create x values for plotting
x_toy = np.linspace(-6, 6, 100)
y_toy = f(x_toy)

# Sample 11 times at 1-interval intervals
x_d = np.arange(-5, 6, 1)
y_d = f(x_d)

# Create a matrix of the radial basis functions
phi_matrix = np.zeros((len(x_d), len(x_d)))

epsilon = 1

for i in range(len(x_d)):
  for j in range(len(x_d)):
    phi_matrix[i, j] = gaussian(np.abs(x_d[i] - x_d[j]), epsilon)

#~~ How do we solve for w_i?
# Take a look at the matrix!

np.set_printoptions(precision=2, suppress=True)
print(phi_matrix)
print(f'The matrix condition number is {np.linalg.cond(phi_matrix):.2e}')
weights = np.linalg.solve(phi_matrix, y_d)

# Define the interpolation function
def interpolation_function(x, weights, x_d, epsilon):
  y = 0
  for i in range(len(x_d)):
    y += weights[i] * gaussian(np.abs(x - x_d[i]), epsilon)
  return y

# Interpolate y_fit
y_fit = [interpolation_function(x, weights, x_d, epsilon) for x in x_toy]

# Plot the results
fig = go.Figure()
fig.add_trace(go.Scatter(x=x_toy, y=y_toy, name='Original Function'))
fig.add_trace(go.Scatter(x=x_d, y=y_d, mode='markers', name='Data Points', marker=dict(color='red')))
fig.add_trace(go.Scatter(x=x_toy, y=y_fit, name='Interpolation', line=dict(dash='dash')))
fig.update_layout(title='Radial Basis Function Interpolation (Gaussian Kernel)', xaxis_title='x', yaxis_title='y')
fig.show()

2D RBF Interpolation¶

Note this is a great result, but it works because the true function tends to zero outside of the data samples.

Example - 2D gaussian¶

# Define the target function
def f(x, y):
    return np.exp(-.5*x**2 - .2*y**2) * np.sin(x) * np.cos(.5*y)

x_grid_vals = np.linspace(-3, 3, 100)
y_grid_vals = np.linspace(-3, 3, 100)
X, Y = np.meshgrid(x_grid_vals, y_grid_vals)
Z_true = f(X, Y)

num_samples = 100
x_samples = np.random.uniform(-3, 3, num_samples)
y_samples = np.random.uniform(-3, 3, num_samples)
z_samples = f(x_samples, y_samples)

def gaussian_rbf(r, epsilon):
    return np.exp(-(epsilon * r)**2)

dist_matrix = np.sqrt((x_samples[:, np.newaxis] - x_samples)**2 + (y_samples[:, np.newaxis] - y_samples)**2)
average_distance = np.mean(dist_matrix)
eps = 1 / average_distance
print(f"Using epsilon = {eps:.4f}")

phi_matrix = gaussian_rbf(dist_matrix, eps)

print(f'The matrix condition number is, {np.linalg.cond(phi_matrix):.2e}')

weights = np.linalg.solve(phi_matrix, z_samples)

Z_fit = np.zeros_like(X)
for i in range(num_samples):
    r = np.sqrt((X - x_samples[i])**2 + (Y - y_samples[i])**2)
    Z_fit += weights[i] * gaussian_rbf(r, eps)

from plotly.subplots import make_subplots

surface_true = go.Surface(x=X, y=Y, z=Z_true, colorscale='Blues', showscale=False, name='Original')
surface_fit = go.Surface(x=X, y=Y, z=Z_fit, colorscale='Viridis', showscale=False, name='RBF Fit')
scatter_samples_1 = go.Scatter3d(
    x=x_samples, y=y_samples, z=z_samples,
    mode='markers', marker=dict(color='red', size=3), name='Samples'
)
scatter_samples_2 = go.Scatter3d(
    x=x_samples, y=y_samples, z=z_samples,
    mode='markers', marker=dict(color='red', size=3), name='Samples'
)

fig = make_subplots(
    rows=1, cols=2,
    specs=[[{'type': 'surface'}, {'type': 'surface'}]],
    subplot_titles=('Original Function', 'RBF Fit')
)

fig.add_trace(surface_true, row=1, col=1)
fig.add_trace(scatter_samples_1, row=1, col=1)

fig.add_trace(surface_fit, row=1, col=2)
fig.add_trace(scatter_samples_2, row=1, col=2)

fig.update_layout(
    title_text='2D RBF Interpolation Comparison',
    scene=dict(
        xaxis_title='x', yaxis_title='y', zaxis_title='z',
        zaxis=dict(range=[-.7, .7]),
        aspectmode='manual',
        aspectratio=dict(x=1, y=1, z=0.5)
    ),
    scene2=dict(
        xaxis_title='x', yaxis_title='y', zaxis_title='z',
        zaxis=dict(range=[-.7, .7]),
        aspectmode='manual',
        aspectratio=dict(x=1, y=1, z=0.5)
    )
)
fig.show()