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.0f (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-15get_n_grid_points(sg)3Taking 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) == 5truef3(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.0get_n_grid_points(sg)9Multi-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.0011607873639181581Reducing 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.0003590890885740948Function Reference
SparseGridsKit.EvaluationDictionary — Typeevaluations_databaseSparseGridsKit.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 approximateff_on_z: Evaluations offon 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
maxptshas been reached, orproftolattained.
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