One-dimensional interpolation and differentiation

In this tutorial, we will create a simple one-dimensional interpolation, investigate how to tune the interpolation method, and show how to apply differential operators on the resulting Interpolation object.

Define problem setup and perform interpolation

We start by defining a simple one-dimensional interpolation problem. We will interpolate the oscillatory function

\[f(x) = \exp(\sin(2x^2)) + 0.1(x - \pi/2)^2\]

between $x = -3$ and $x = 3$. For simplicity, we take 25 equidistant points in the interval $[-3, 3]$ as interpolation points.

using KernelInterpolationf(x) = exp(sin(2*x[1]^2)) + 0.1*(x[1] - pi/2)^2x_min = -3x_max = 3N = 25nodeset = NodeSet(LinRange(x_min, x_max, N))values = f.(nodeset)

25-element Vector{Float64}:
 2.561118346612702
 3.6010482803853834
 2.592967848165464
 1.9848324553623153
 3.9645665583549765
 1.957023058914906
 1.3192182268719703
 1.8124194938335236
 3.143477103401214
 3.003796380105595
 ⋮
 2.5151585726832555
 1.0270213304360754
 0.3767404307950323
 0.8574656301584783
 2.7079294969190593
 0.5711157612469084
 1.0221715213705675
 1.8731723209109972
 0.676162754458826

Next, we choose the kernel (radial basis function) for the interpolation. We use the Gaussian kernel with a fixed shape parameter of 0.5 and interpolate the function values.

kernel = GaussKernel{1}(shape_parameter = 0.5)itp = interpolate(nodeset, values, kernel)

Interpolation with 25 nodes, kernel GaussKernel{1}(shape_parameter = 0.5) and polynomial of order 0.

Let's plot the interpolated function and the original function on a finer grid to see how well the interpolation works.

using Plotsmany_nodes = NodeSet(LinRange(x_min, x_max, 200))plot(many_nodes, f, label = "Original function")scatter!(nodeset, f, label = "Original data", markershape = :star)plot!(many_nodes, itp, training_nodes = false, yrange = (0.0, 5.0))

Uhh, that doesn't look too good. What happened?

Finding a well-suited interpolation method

We used the GaussKernel with a rather small shape parameter of 0.5, which leads to an ill-conditioned linear system of equations. We can inspect the condition number of the interpolation matrix to confirm this.

using LinearAlgebraA = system_matrix(itp)cond(A)

9.52474836763431e17

Here, we used the system_matrix function to obtain the interpolation matrix A and calculated the condition number of the matrix. For this specific example the system matrix simply is the kernel_matrix, but for more sophisticated interpolations the system matrix contains additional parts like the polynomial augmentation. The condition number is a measure of how well-conditioned the matrix is. A large condition number indicates that the matrix is ill-conditioned, which usually leads to high numerical errors. To avoid this, we have different options. We can either increase the shape parameter of the kernel or we can use a different kernel. The GaussKernel is known to be rather ill-conditioned and other kernels like the WendlandKernel usually lead to better condition numbers. Here, we choose to increase the shape parameter of the Gaussian kernel to 1.8, which makes the interpolation more localized. Note, however, that you might need to choose another kernel if you increase the number of interpolation points.

kernel = GaussKernel{1}(shape_parameter = 1.8)itp = interpolate(nodeset, values, kernel)plot(many_nodes, f, label = "Original function")plot!(many_nodes, itp, yrange = (0.0, 5.0))

We can see a much better agreement between the original function and the interpolated function. We still observe some undershoots, but this is expected due to the oscillatory nature of the function and the limited number of interpolation points. Let's confirm that increasing the shape parameter improved the condition number of the interpolation matrix.

A = system_matrix(itp)cond(A)

57037.29656744214

Indeed, the condition number is much smaller than before!

Applying differential operators

Sometimes, we are not only interested in interpolating a function, but also in computing its derivatives. Remember that in the simplest case, where no polynomial augmentation is used, the interpolation itp represents a linear combination

\[s(x) = \sum_{j = 1}^N c_j\phi(\|x - x_j\|_2)\]

with $\phi$ given by the radial basis function, in this case the Gaussian. Because we know $\phi$ and its derivatives, we can compute the derivatives of $s$ by differentiating the kernel function. For a general dimension $d$, the partial derivative in the $i$-th direction, $i\in\{1,\ldots,d\}$, of the interpolation is then given by

\[\frac{\partial s}{\partial x_i}(x) = \sum_{j = 1}^N c_j\frac{\partial \phi}{\partial x_i}(\|x - x_j\|_2).\]

KernelInterpolation.jl already provides some common differential operators. For example, we can compute the first derivative of the interpolation itp at a specific point x by using the PartialDerivative operator.

d1 = PartialDerivative(1)x = 0.0itp_dx_at_x = d1(itp, x)

-0.3128030377986403

Let's plot the first derivative of the interpolated function and compare it to the analytical first derivative.

itp_dx = d1(itp)itp_dx_many_nodes = itp_dx.(many_nodes)f_dx(x) = 4*exp(sin(2*x[1]^2))*x[1]*cos(2*x[1]^2) + 0.2*x[1] - pi/10plot(many_nodes, f_dx, label = "Derivative of original function")plot!(many_nodes, itp_dx_many_nodes, label = "Derivative of interpolated function")

Differentiating the interpolation process

So far, we differentiated the interpolation $s$ with respect to its argument $x$, keeping the coefficients $c_j$, i.e. the solution of the linear system, fixed. The whole interpolation process is differentiable as well: the assembly of the system matrix and the linear solve determining the coefficients are ordinary Julia code and can be differentiated with forward-mode automatic differentiation, e.g. with ForwardDiff.jl.

This is useful, for example, to tune the shape parameter of the kernel, which we chose by hand above: the derivative of some error measure with respect to the shape parameter can be obtained by differentiating through interpolate and can then be handed to an optimizer.

using ForwardDifffunction loss(epsilon)    itp = interpolate(nodeset, values, GaussKernel{1}(shape_parameter = epsilon))    return sum(abs2, itp.(many_nodes) .- f.(many_nodes))endForwardDiff.derivative(loss, 1.8)

7.857515631536312

Note that this requires the interpolation system to be well-conditioned. For the shape parameter of 0.5 we started with, the coefficients themselves are dominated by round-off errors, and so is any derivative of them.

The same works for the interpolated values, for the positions of the nodes, and for solving PDEs with solve_stationary. Since the interpolation depends linearly on the values it interpolates, the gradient with respect to them is the vector of the cardinal functions evaluated at the point of interest. For a point close to the interpolation node at $x = 1.0$, the largest entries are the ones belonging to the surrounding nodes.

x_eval = [1.1]ForwardDiff.gradient(v -> interpolate(nodeset, v, kernel)(x_eval), values)

25-element Vector{Float64}:
  0.00033656986595132047
 -0.0012382336868100405
  0.002733563850252959
 -0.004765614901028104
  0.007289621772069951
 -0.01031421408166926
  0.013911550173629028
 -0.018221232974148602
  0.02346298686249623
 -0.02997040898910763
  ⋮
  0.7663323466000272
  0.49071182680768566
 -0.17075354737465986
  0.09319250724820581
 -0.056076538213339445
  0.03337087305492121
 -0.018064485354392097
  0.007887446073231866
 -0.002093611478296018