Sparse Grid Integration

Sparse Grid Quadrature

To integrate a sparse grid polynomial approximation we require the sparse grid, function evaluations and the qudarature weights. These are precomputed when the sparse grid is constructed via the compute_quadrature_weights! function.

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

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

f (generic function with 1 method)

and we construct the two-dimension function $[f(x_1), f(x_2)^2]^{\top}$ mapping $[-1,1]^2 \to \mathbb{R}^2$. We seek the value of the integral

\[\int_{\Gamma} [f(x_1), f(x_2)^2]^{\top} \rho(x) \textrm{d} x\]

for a weight function $\rho(x)=0.5^2$.

The constant approximation is $f(x)\approx[1,1]$. The approximate integral using the constant approximation is equal to $[1,1]$.

using SparseGridsKit
n,k = 2,0
mi_set = create_smolyak_miset(n,k)

sg = create_sparsegrid(mi_set)
f_on_grid = [[f(x[1]), f(x[2])^2] for x in get_grid_points(sg)]
# Test integrate_on_sparsegrid
integral_result = integrate_on_sparsegrid(sg, f_on_grid)
result = [1, 1]
println("Quadrature: $integral_result, Integral: $result")

Quadrature: [1.0, 1.0], Integral: [1, 1]

The estimated integral is improved by using a sparse grid with more grid points. The integral can be evaluated explicitly, and it is seen that using an larger sparse grid gets the correct result $[2.0,92/15]$.

n,k = 2,4
mi_set = create_smolyak_miset(n,k)
sg = create_sparsegrid(mi_set)
f_on_grid = [[f(x[1]), f(x[2])^2] for x in get_grid_points(sg)]
# Test integrate_on_sparsegrid
integral_result = integrate_on_sparsegrid(sg, f_on_grid)
result = [2.0,92/15]
println("Quadrature: $integral_result, Integral: $result")

Quadrature: [1.9999999999999998, 6.133333333333334], Integral: [2.0, 6.133333333333334]

This is of course expected –- the polynomial approximation space represented by the larger sparse grid contains the function $[f(x_1), f(x_2)^2]$.

Underlying the integral is the computation of quadrature weights which are computed using compute_quadrature_weights!. This is done automatically when a grid is constructed.

compute_quadrature_weights!(sg)
sg.quadrature_weights

65-element Vector{Float64}:
 -0.043436041083099924
  0.013371390902239579
 -0.04819372348784111
  0.030450582331498402
 -0.30814013872837404
  0.030450582331498402
 -0.04819372348784111
  0.013371390902239579
 -0.043436041083099924
  0.018684351418602807
  ⋮
  0.012184887434668178
  0.018684351418602807
  0.054452776290945464
  0.08158633214085165
  0.09625693230646282
  0.09625693230646282
  0.08158633214085165
  0.054452776290945464
  0.018684351418602807

Sparse Grid Weighted $L_{\rho}^2(\Gamma)$ Norms

Often the weighted $L_{\rho}^2(\Gamma)$ norm is required. Squaring a polynomial approximation results in many higher degree polynomial terms –- it is not as simple as squaring the evaluations at each sparse grid point.

To approximate this integral this efficiently, integrals of the products of interpolation polynomials are precomputed. A maximum level number is selected and for each domain dimension all pairwise $L_{\rho}^2(\Gamma)$ integrals inner products are computed using precompute_lagrange_integrals.

For example, the pairwise integrals are computed, and we can extract the pairwise inner products for the level $2$ and level $3$ polynomials.

using SparseGridsKit, LinearAlgebra
pcl = precompute_lagrange_integrals(7)
level1 = 2
level2 = 3
domaindim = 1
pcl[domaindim][level1,level2,:,:]

65×65 Matrix{Float64}:
  0.0380952    0.151424  0.00952381  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 -0.00952381   0.152381  0.380952       0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0047619   -0.037138  0.00952381     0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0          0.0       0.0            0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0          0.0       0.0            0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0          0.0       0.0         …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0          0.0       0.0            0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0          0.0       0.0            0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0          0.0       0.0            0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0          0.0       0.0            0.0  0.0  0.0  0.0  0.0  0.0  0.0
  ⋮                                  ⋱            ⋮                   
  0.0          0.0       0.0            0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0          0.0       0.0            0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0          0.0       0.0            0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0          0.0       0.0            0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0          0.0       0.0         …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0          0.0       0.0            0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0          0.0       0.0            0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0          0.0       0.0            0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0          0.0       0.0            0.0  0.0  0.0  0.0  0.0  0.0  0.0

This precomputation step is generally not too expensive but uses a moderate amount of memory.

@elapsed precompute_lagrange_integrals(7)

0.005006233

sizeof(pcl)

8

Additional arguments must be provided to precompute_lagrange_integrals if specific knots and rules are use.

The precomputed $L_{\rho}^2(\Gamma)$ integrals from precompute_lagrange_integrals are used with the sparse grid approximation structure to exactly compute the $L_{\rho}^2(\Gamma)$ norm of a sparse grid polynomial function. These are paired with appropriate pairwise products of the function evaluation which are computed with precompute_pairwise_norms. For example, for a scalar valued function we may require

\[f(x_1) f(x_2) \forall x_1,x_2 \in Z\]

where $Z$ is the set of sparse grid points.

For example, the function $f(x)=x$ has the weighted $L_{\rho}^2([-1,1])$ norm equal to

\[\Vert x^2 \Vert_{L_{\rho}^2([-1,1])} = \frac{1}{\sqrt{3}}\]

for weight function $\rho=0.5$.

n,k = 1,1
mi_set = create_smolyak_miset(n,k)
pcl = precompute_lagrange_integrals(7)

sg = create_sparsegrid(mi_set)
f_on_grid = [x[1] for x in get_grid_points(sg)]
pairwise_norms = precompute_pairwise_norms(f_on_grid, product=(x,y)->dot(x,y))
l2_integral_result = integrate_L2_on_sparsegrid(sg, f_on_grid, pcl)
result = sqrt(1/3)
println("Quadrature: $l2_integral_result, Results: $result")

Quadrature: 0.5773502691896258, Results: 0.5773502691896257

Similarly, for $x^2$ we get

\[\Vert x^2 \Vert_{L_{\rho}^2([-1,1])}=\frac{1}{\sqrt{5}}.\]

n,k = 1,3
mi_set = create_smolyak_miset(n,k)
pcl = precompute_lagrange_integrals(7)

sg = create_sparsegrid(mi_set)
f_on_grid = [x[1]^2 for x in get_grid_points(sg)]
l2_integral_result = integrate_L2_on_sparsegrid(sg, f_on_grid, pcl)
result = sqrt(1/5)
println("Quadrature: $l2_integral_result, Results: $result")

Quadrature: 0.447213595499958, Results: 0.4472135954999579

Finally we consider the $L^2_{\rho}([-1,1]^2)$ norm of $f(x_1)$. This can be explicitly computed to be $\sqrt{92/15}$.

n,k = 2,4
mi_set = create_smolyak_miset(n,k)
pcl = precompute_lagrange_integrals(7)

sg = create_sparsegrid(mi_set)
f_on_grid = [f(x[1]) for x in get_grid_points(sg)]
pairwise_norms = precompute_pairwise_norms(f_on_grid, product=(x,y)->x.*y)
l2_integral_result = integrate_L2_on_sparsegrid(sg, f_on_grid, pcl)
result = sqrt(92/15)
println("Quadrature: $l2_integral_result, Results: $result")

Quadrature: 2.476556749467556, Results: 2.4765567494675613

If the evaluation of the function $f:\Gamma\maps\to\R^n$ represents a discrete approximation of a function $F$, we may use a mass matrix $Q$ to compute

\[(x_1^{\top} Q x_2)^{1/2} \quad \forall x_1,x_2 \in Z.\]

This is useful for integrals of the type $L^2_\rho(\Gamma; L^2(D))$ where $D$ represents the domain of the approximated function $F$ and is often seen in the approximation of parametric PDE solutions. This is achieved by passing product=(x,y)->dot(x,Q,y) as an argument in precompute_pairwise_norms.

Function Reference

compute_quadrature_weights!(sg)

integrate_on_sparsegrid(sg, f_on_grid, precomputed_lagrange_integrals)

integrate_L2_on_sparsegrid(sg, f_on_grid, precomputed_lagrange_integrals; product=dot, precomputed_pairwise_norms=nothing)