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

import numpy as np
import plotly.graph_objects as go
from scipy.interpolate import RBFInterpolator

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

Best Fit Radial Basis Functions¶

Since RBFs were also solved using a linear system, the extension to least square regression is straightforward. Recall that since for most kernels $\varphi_i(\|x-x_i\| \rightarrow \infty)\rightarrow 0$ , the approximated function will go to zero far from the data.

The modern implementation of RBFs accounts for the global trend of the data through a polynomial least squares fit alongside normal RBFs for local features.

\begin{aligned} y(x) &= \sum_i^N \omega_i \varphi_i(\|x-x_i\|) + \sum_i^N P_i(x_i) b_i \end{aligned}

(1)

Where $P_i$ is an order $n<m$ polynomial. The Numpy RBFInterpolator object fits this equation to:

\begin{aligned} [\Phi(x_i, x_j) -\lambda I]\omega +P(x_i) b &= y_i \\ P(x_i)^T \omega &=0 \end{aligned}

(2)

where $\lambda = 0$ recovers an exact fit and $\lambda > 0$ effectively shifts the fitting of the $x_i=x_j$ terms to the best fit polynomial.

Example: Toy Gaussian over a Quadratic¶

Let’s sample a Gaussian placed over a quadratic function, and use a scipy.interpolate.RBFInterpolator over the sampled points to see how it recovers the true function.

# Define the function
def f(x):
  return np.exp(-(x/2)**2) + 0.1*x**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 an RBFInterpolator object
rbf = RBFInterpolator(np.array([x_d]).T, y_d.T, kernel='gaussian', epsilon=1, degree=2)

# Interpolate at the x_toy values
y_rbf = rbf(np.array([x_toy]).T)

# Plot the results using Plotly
fig = go.Figure()
fig.add_trace(go.Scatter(x=x_toy, y=y_toy, mode='lines', name='True Function'))
fig.add_trace(go.Scatter(x=x_d, y=y_d, mode='markers', marker=dict(size=8, color='red'), name='Data Points'))
fig.add_trace(go.Scatter(x=x_toy, y=y_rbf, mode='lines', line=dict(dash='dash', color='green'), name='RBF Interpolation'))

fig.update_layout(title='1D RBF Interpolation', xaxis_title='x', yaxis_title='y')
fig.show()

Example: RBF Interpolation in 2D¶

We can also naturally extend RBFs to 2D surfaces, which is where they truly shine!

# 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)

# Prepare sample points for RBFInterpolator
sample_points = np.column_stack((x_samples, y_samples))

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

# Fit RBF interpolator (using Gaussian kernel, can adjust epsilon and degree as needed)
rbf = RBFInterpolator(sample_points, z_samples, kernel='gaussian', epsilon=1, degree=1)

# Evaluate on grid
grid_points = np.column_stack((X.ravel(), Y.ravel()))
Z_fit = rbf(grid_points).reshape(X.shape)

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 = go.Scatter3d(
    x=x_samples, y=y_samples, z=z_samples,
    mode='markers', marker=dict(color='red', size=3), name='Samples'
 )

from plotly.subplots import make_subplots

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

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

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

fig.update_layout(
    title='2D RBF Interpolation Comparison',
    scene=dict(
        xaxis_title='x', yaxis_title='y', zaxis_title='z',
        zaxis=dict(range=[-.7, .7]),
        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]),
        aspectratio=dict(x=1, y=1, z=0.5)
    ),
    width=1000, height=500
)

fig.show()