Adaptive Sparse Grids
Sparse grids can be adaptively constructed using a Gerstner–Griebel [4] type construction. This is implemented in the adaptive_sparsegrid
function.
One dimensional example
To illustrate the adaptive algorithm consider functions that can be exactly approximated using a sparse grid approximation. Consider the polynomial $f:[-1,1]\mapsto\R$
\[f(x) = 3x^2 + 2x +1.\]
using SparseGridsKit
nparams = 1
f(x) = @. 3.0*x[1]^2 + 2.0*x[1] + 1.0
f (generic function with 1 method)
The adaptive_sparsegrid
function can be called with a function $f$ and the dimension of the function domain nparams
.
(sg, f_on_Z) = adaptive_sparsegrid(f, nparams)
SparseGridApproximation(sg, f_on_Z)
SparseGridApproximation(SparseGrid(1, [[-1.0, 1.0]], [1 2], [0, 1], Real[1.0; 0.0; -1.0;;], [0.16666666666666669, 0.6666666666666667, 0.16666666666666669], Points[CCPoints{Float64}([-1.0, 1.0])], Level[Doubling()], SparseGridsKit.sparse_grid_data{1}(SparseGridsKit.SparseTerm{1}[SparseGridsKit.SparseTerm{1}(([2, 1],)), SparseGridsKit.SparseTerm{1}(([2, 2],)), SparseGridsKit.SparseTerm{1}(([2, 3],))], [(2,), (1,), (3,)], [1, 2, 3], [1, 1, 1], [Dict([2, 3] => 3, [2, 1] => 2, [1, 1] => 1, [2, 2] => 1)], [[[0.0], [1.0, 0.0, -1.0]]], [[0.0, 1.0, -1.0]], [[[1.0], [0.16666666666666669, 0.6666666666666667, 0.16666666666666669]]])), [6.0, 1.0, 2.0], [[-1.0, 1.0]])
The function f
can be exactly represented by a three point interpolant, which is identified in the adaptive algorithm.
using LinearAlgebra
test_points = range(-1,stop=1,length=100)
error = norm([[f(x)] - interpolate_on_sparsegrid(sg,f_on_Z,x) for x in test_points])
println("Error: $error")
Error: 3.6502539264068535e-15
get_n_grid_points(sg)
3
Taking powers $k$ of the polynomial $f$ gives a polynomial $f^k$ of polynomial degree $2^k$. Consequently, this is represented exactly using $2^k+1$ interpolation points. The adaptive algorithm identifies this.
f2(x) = f(x).^2
(sg, f2_on_Z) = adaptive_sparsegrid(f2, nparams)
# Expect three point approx cubic (1 iteration to suffice)
error = norm([f2(x)] ≈ interpolate_on_sparsegrid(sg,f2_on_Z,x) for x in test_points),
println("Error: $error")
(10.0, nothing)
get_n_grid_points(sg) == 5
true
f3(x) = f(x).^3
(sg, f3_on_Z) = adaptive_sparsegrid(f3, nparams)
# Expect three point approx cubic (1 iteration to suffice)
error = norm([f3(x)] ≈ interpolate_on_sparsegrid(sg,f3_on_Z,x) for x in test_points)
println("Error: $error")
Error: 10.0
get_n_grid_points(sg)
9
Multi-dimensional example
The adaptive algorithm is also applicable for higher dimensional functions. Consider the genz
gaussianpeak
example.
using SparseGridsKit, LinearAlgebra
n = 8
C = 1.0
W = 0.0
T = "exponentialdecay"
N = "gaussianpeak"
f = genz(n::Int, C::Float64, W::Float64, T::String, N::String)
#5 (generic function with 1 method)
An adaptive approximation is formed. This produces a (relatively) accurate approximation of the high-dimensional function $f$. The iterations halt when the stopping criterion of max profit p_max < proftol
.
# Approximate
(sg, f_on_Z) = adaptive_sparsegrid(f, n)
nmc = 100
test_points = [2 * rand(8) .- 1 for ii = 1:nmc]
f_approx_on_test = interpolate_on_sparsegrid(sg,f_on_Z,test_points)
norm([f(x) - f_approx_on_test[i] for (i,x) in enumerate(test_points)])
0.0015471677772621146
Reducing proftol
allows further iterations and a more accurate approximation.
(sg, f_on_Z) = adaptive_sparsegrid(f, n, proftol=1e-5)
f_approx_on_test = interpolate_on_sparsegrid(sg,f_on_Z,test_points)
norm([f(x) - f_approx_on_test[i] for (i,x) in enumerate(test_points)])
0.00041235982884738363
Function Reference
SparseGridsKit.EvaluationDictionary
— Typeevaluations_database
SparseGridsKit.EvaluationDictionary
— MethodEvaluationDictionary(sg,f_on_z)
Creates an EvaluationDictionary from a sparse grid and function evaluations
Arguments
sg
: Sparse gridf_on_z
: Function evaluations on the sparse grid
Returns
evaluations
: EvaluationDictionary
SparseGridsKit.EvaluationDictionary
— MethodEvaluationDictionary()
Creates an empty EvaluationDictionary
Returns
evaluations
: EvaluationDictionary
SparseGridsKit.adaptive_estimate
— Methodadaptive_estimate(sg, datastore, pcl, type=:deltaint)
Estimates the profit of adding multi-indices {α} in reduced margin to the sparse grid.
Arguments
sg
: Sparse grid.datastore
: EvaluationDictionary.rule
: Level function(s) for sparse grid.knots
: Knot function(s) for sparse grid.type
: Type of profit computation. Default is:deltaint
.
Returns
- Vector of profits for each multi-index α.
SparseGridsKit.adaptive_mark
— Functionadaptive_mark(α, p_α, θ=1e-4)
Dorfler style marking of multi-indices based on profits
Arguments
α
: Multi-indices in reduced margin.p_α
: Vector of profits.θ
: Threshold for marking. Default is1e-4
.
SparseGridsKit.adaptive_refine
— Methodadaptive_refine(sg, α_marked, rule, knots)
Refines the sparse grid based on marked multi-indices
Arguments
sg
: Sparse grid.α_marked
: Marked multi-indices.rule
: Level function(s) for sparse grid.knots
: Knot function(s) for sparse grid.
Returns
sg
: Refined sparse grid.f_on_z
: Function evaluations on the refined sparse grid.
SparseGridsKit.adaptive_solve!
— Methodadaptive_solve!(f, datastore, sg)
Adds function evaluations on sg to datastore
Arguments
f
: Functiondatastore
: EvaluationDictionarysg
: Sparse grid
Returns
datastore
: Updated EvaluationDictionary
SparseGridsKit.adaptive_sparsegrid
— Methodadaptive_sparsegrid(f, nparams; maxpts = 100, proftol = 1e-4, rule=Doubling(), knots=CCPoints(), θ=1e-4, type=:deltaint)
Constructs an adaptive sparse grid for approximating the function f
in nparams
dimensions.
Arguments
f
: Function to be approximated taking arguments x with length(x)=nparams.nparams
: Dimension of domain.maxpts
: (Optional) Maximum number of points to include in the sparse grid. Default is100
.proftol
: (Optional) Tolerance for profits. Default is1e-4
.rule
: (Optional) Level function(s) for sparse grid. Default isDoubling()
.knots
: (Optional) Knot function(s) for sparse grid. Default isCCPoints()
.θ
: (Optional) Threshold for marking. Default is1e-4
.type
: (Optional) Type of profit computation. Default is:deltaint
.costfunction
(Optional) Cost function for fidelity related multi-indices.
Returns
sg
: Final sparse grid used to approximatef
f_on_z
: Evaluations off
on grid points in sparse gridsg
SparseGridsKit.adaptive_terminate
— Methodadaptiveterminate(sg, pα, maxpts, proftol)
Test profits and sparse grid to determine loop termination
Arguments
sg
: Sparse grid.p_α
: Vector of profits.maxpts
: Maximum number of grid points allowed.proftol
: Profit tolerance.
Returns
- Flag to indicate loop termination if
maxpts
has been reached, orproftol
attained.
SparseGridsKit.add_evaluation!
— Methodadd_evaluation!(evaluations, z, f_at_z)
Adds the evaluation fatz into the evaluations_database
Arguments
evaluations
: EvaluationDictionaryz
: Pointf_at_z
: Function evaluation at z
Returns
evaluations
: Updated EvaluationDictionary
SparseGridsKit.add_evaluations!
— Methodadd_evaluations!(evaluations, sg::SparseGrid, f_on_z)
Adds the evaluations fonz into the evaluations_database
Arguments
evaluations
: EvaluationDictionarysg
: Sparse gridf_on_z
: Function evaluations on the sparse grid
Returns
evaluations
: Updated EvaluationDictionary
SparseGridsKit.compute_profit
— Methodcompute_profit(sg_α, f_α, f, cost, pcl; type=:deltaint)
Computes the "profit" of a sparse grid supplemented by a multi-index α
Arguments
sg_α
: Enhanced sparse grid with αsg
: Sparse gridf_α
: Function evaluations at grid points in αf
: Function evaluations on the sparse grid.cost
: Cost of adding α to the sparse grid.type
: Type of profit computation. Default is:deltaint
. Options are:deltaint
,:deltaintcost
,:Linf
,:Linfcost
.
Returns
- Computed profit as Expected change in approximation.
SparseGridsKit.retrieve_evaluations
— Methodretrieve_evaluations(evaluations, sg)
Retrieves the evaluations from the evaluations_database
Arguments
evaluations
: EvaluationDictionarysg
: Sparse grid
Returns
f_on_z
: Function evaluations on the sparse grid