Sparse Grid Interpolation

Once a sparse grid is constructed, a polynomial approximation is formed by pairing the grid points with function evaluations.

If the multi-index set is admissible, and the sequences of sets of points are nested this approximation is an interpolating polynomial.

One-dimensional Interpolation

Consider a simple example approximating the function $f:[-1,1]\to\mathbb{R}$

\[f(x) = 3x^2 + 2x + 1\]

Ths is expressed as the function f.

f(x) = @. 3.0 * x[1]^2 + 2.0 * x[1] + 1.0

f (generic function with 1 method)

Notice that the function f returns a scalar. It is not essential, but it is advisable to use vector valued functions due to the mutability. In this simple example, the extra allocations due to using scalars will not be an issue. !!!

We can construct a one point approximation (i.e. constant polynomial). The function f is evaluated on the grid points accessed via get_grid_points.

n, k = 1, 0
using SparseGridsKit, LinearAlgebra
mi_set = create_smolyak_miset(n, k)
sg = create_sparsegrid(mi_set)
f_on_grid = [f(x) for x in get_grid_points(sg)]

1-element Vector{Float64}:
 1.0

The sparse grid sg is paired with the point evaluations to interpolate on to target points. Notice that whilst we are in one-dimension, target_points is still a Vector of Vectors

target_points = [[x] for x in -1.0:0.1:1.0]
interpolation_result = interpolate_on_sparsegrid(sg, f_on_grid, target_points)

21-element Vector{Float64}:
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 ⋮
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0

The constant approximation is $f\approx 1$.

It is clear a better approximation can be formed. This is attained with a higher degree polynomial approximation – derived from a larger multi-index set.

mi_set_new = create_smolyak_miset(n, 1)
sg_new = create_sparsegrid(mi_set_new)
f_on_grid_new = [f(x) for x in get_grid_points(sg_new)]

interpolation_result_new = interpolate_on_sparsegrid(sg_new, f_on_grid_new, target_points)

error = norm([interpolation_result_new[i] - f(target_points[i]) for i in eachindex(target_points)])
println("ℓ^2 Error: $error")

ℓ^2 Error: 1.0648885158809838e-15

This is illustrated in the following plot.

plot(SparseGridApproximation(sg,f_on_grid); label="Initial Grid")
plot!(SparseGridApproximation(sg_new,f_on_grid_new); label="Enriched Grid")
plot!(f; label="f")

The number of grid points can be checked using get_n_grid_points and is seen to equal three. It is therefore no surprise that the approximation error is at machine precision - we can exactly represent the degree two polynomial f as a three point interpolant.

get_n_grid_points(sg_new)

3

Multi-dimensional Interpolation

These ideas extend to vector valued functions with multi-dimensional domains. Consider a function

\[f:[-1,1]^4 \to \mathbb{R}^{400}\]

defined via the gaussianpeak Genz test function (see e.g. [2]). A polynomial approximation can again be formed using nested Clenshaw-Curtis points and Smolyak multi-index sets.

using SparseGridsKit, LinearAlgebra, Plots, LaTeXStrings
n, k = 4, 3
mi_set = create_smolyak_miset(n, k)
sg = create_sparsegrid(mi_set)

# Define a complicated function
nparams=100
f(x) = [genz(n, 1.0, w, "quadraticdecay", "gaussianpeak")(x) for (i,w) = enumerate(range(-1,stop=1,length=nparams))]
f_on_grid = [f(x) for x in get_grid_points(sg)]

xplot = range(-1,stop=1,length=100)
plot(xplot, [f([xi, zeros(3)...])[1] for xi in xplot],label=L"f_{1}(x)")
plot!(xplot, [f([xi, zeros(3)...])[25] for xi in xplot],label=L"f_{25}(x)")
plot!(xplot, [f([xi, zeros(3)...])[50] for xi in xplot],label=L"f_{50}(x)")
plot!(xplot, [f([xi, zeros(3)...])[100] for xi in xplot],label=L"f_{100}(x)")
plot!(
    title=L"f(x_1,0,0,...)",
    xlabel=L"x_1"
)

This sparse grid approximation is constructed and evaluated at $100$ test points. The norm of the error vector is computed and is small, but certainly not at machine precision.

nmc = Integer(1e2);
V = [2 * rand(n) .- 1 for ii = 1:nmc]
f_on_v = [interpolate_on_sparsegrid(sg, f_on_grid, [v])[1] for v in V]
error = norm([f_on_v[i] - f(V[i]) for i in eachindex(V)])
println("ℓ^2 Error: $error")

ℓ^2 Error: 0.0046594956954737685

If we add multi-indices to the multi-index set we should have a better, higher degree polynomial approximation. The function $f$ is analytic so we expect the polynomial approximation to converge quickly.

n, k = 4, 6
mi_set = create_smolyak_miset(n, k)
sg = create_sparsegrid(mi_set)
f_on_grid = [f(x) for x in get_grid_points(sg)]
f_on_v = interpolate_on_sparsegrid(sg, f_on_grid, V)
error = norm([f_on_v[i] - f(V[i]) for i in eachindex(V)])
println("ℓ^2 Error: $error")

ℓ^2 Error: 1.9479577654287686e-8

n, k = 4, 7
mi_set = create_smolyak_miset(n, k)
sg = create_sparsegrid(mi_set)
f_on_grid = [f(x) for x in get_grid_points(sg)]
f_on_v = interpolate_on_sparsegrid(sg, f_on_grid, V)
error = norm([f_on_v[i] - f(V[i]) for i in eachindex(V)])
println("ℓ^2 Error: $error")

ℓ^2 Error: 1.65564873190907e-10

Later, the function adaptive_sparsegrid will consider a greedy construction of the multi-index set to reduce the approximation error.

Function Reference

interpolate_on_sparsegrid(sg, f_on_grid, target_points)

interpolate_on_sparsegrid(sga::SparseGridApproximation, target_points)