Sparse grid
A sparse grid is a multidimensional discretization constructed from a restricted combination of tensor-product grids. It is used for numerical integration, interpolation, and the approximation of solutions to partial differential equations. Sparse grids retain much of the approximation order of a full tensor-product construction while requiring substantially fewer sampling points when the dimension is moderate and the target function has bounded mixed derivatives.
The standard construction is based on a hierarchy of one-dimensional approximation operators. Their tensor products are selected according to a constraint on the sum of their refinement levels, rather than allowing every coordinate to reach the maximum level simultaneously. This selection changes the asymptotic point count from exponential growth in the dimension to a product of an essentially one-dimensional term and a dimension-dependent power of the refinement level. The reduction does not remove the curse of dimensionality, but it modifies its effect for function classes possessing mixed regularity.
Sparse grids are closely associated with the combination formula introduced by Sergey A. Smolyak in 1963. Later work connected the same construction to hierarchical finite-element spaces, adaptive interpolation, and stochastic collocation. The term also refers more broadly to anisotropic and dimension-adaptive variants derived from the original level-selection principle.
Mathematical construction
Let (U_i) denote a one-dimensional approximation operator at level (i), and define the hierarchical increment
[ \Delta_i = U_i-U_{i-1}, ]
with (U_0=0). For a (d)-dimensional function, the full tensor-product operator at uniform level (n) is
[ U_n^{\otimes d}
\sum_{1\leq i_1,\ldots,i_d\leq n} \Delta_{i_1}\otimes\cdots\otimes\Delta_{i_d}. ]
This decomposition contains contributions for every multi-index in the corresponding (d)-dimensional box. Many of those contributions represent simultaneous high resolution in all coordinates and therefore account for most of the degrees of freedom.
The isotropic sparse-grid operator retains only increments satisfying a total-level condition:
[ A_{n,d}
\sum_{\lvert\mathbf{i}\rvert_1\leq n+d-1} \Delta_{i_1}\otimes\cdots\otimes\Delta_{i_d}, ]
where (\mathbf{i}=(i_1,\ldots,i_d)) and (\lvert\mathbf{i}\rvert_1=i_1+\cdots+i_d). Indexing conventions vary by the initial level assigned to the one-dimensional hierarchy, but the selected multi-indices form the same type of discrete simplex.
The construction can also be written as a signed sum of ordinary tensor-product operators:
[ A_{n,d}
\sum_{n\leq \lvert\mathbf{i}\rvert_1\leq n+d-1} (-1)^{n+d-1-\lvert\mathbf{i}\rvert_1} \binom{d-1}{n+d-1-\lvert\mathbf{i}\rvert_1} U_{i_1}\otimes\cdots\otimes U_{i_d}. ]
This expression is the Smolyak algorithm. It permits the constituent tensor problems to be evaluated separately before their results are combined with integer coefficients. Repeated nodes can be identified when the one-dimensional rules are nested, whereas non-nested rules generally produce a larger union of points.
Complexity and approximation error
Suppose a one-dimensional grid at level (i) contains (O(2^i)) points. A full (d)-dimensional tensor grid at level (n) then contains
[ O(2^{nd}) ]
points. Under the same scaling, a conventional isotropic sparse grid contains
[ O!\left(2^n n^{d-1}\right) ]
points. The power (n^{d-1}) arises from the number of multi-indices near the boundary of the total-level set.
For piecewise multilinear interpolation of a function with bounded mixed second derivatives, the error commonly has the form
[ O!\left(2^{-2n}n^{d-1}\right). ]
In terms of the total number (N) of grid points, this becomes a one-dimensional algebraic rate accompanied by powers of (\log N). The precise exponents depend on the selected basis, the norm used to measure error, and the mixed smoothness of the approximated function. A function that is smooth in every coordinate separately does not automatically satisfy the required mixed-derivative conditions.
The relevant function spaces are often formulated through Sobolev spaces with dominating mixed smoothness. Their norms control derivatives involving combinations of coordinate directions. Sparse grids are therefore distinguished from unstructured point reduction: their complexity estimates follow from a particular relation between hierarchical basis coefficients and mixed regularity.
Hierarchical representation
Sparse-grid interpolation is commonly expressed in a hierarchical basis. At each level, basis functions are associated with points not already represented by coarser levels. The coefficient of such a function is its hierarchical surplus, which measures the difference between the sampled value and the approximation inherited from preceding levels.
For a multi-index (\mathbf{i}) and a point index (\mathbf{j}), the approximation may be written as
[ f_{n,d}(\mathbf{x})
\sum_{\lvert\mathbf{i}\rvert_1\leq n+d-1} \sum_{\mathbf{j}\in B_{\mathbf{i}}} \alpha_{\mathbf{i},\mathbf{j}} \phi_{\mathbf{i},\mathbf{j}}(\mathbf{x}), ]
where (B_{\mathbf{i}}) identifies the hierarchical points, (\phi_{\mathbf{i},\mathbf{j}}) is a tensor-product basis function, and (\alpha_{\mathbf{i},\mathbf{j}}) is the associated surplus. For piecewise linear bases, each one-dimensional factor has compact support. Polynomial constructions instead use rules related to Lagrange interpolation, often on nested abscissae such as Clenshaw–Curtis quadrature points.
Hierarchical surpluses also provide local information about unresolved variation. Adaptive constructions enlarge the index set or refine selected spatial regions according to the contribution of these coefficients. The resulting grid no longer has the exact symmetry of the isotropic Smolyak set, although it retains downward closure: an included fine-level index is accompanied by the coarser indices required for its hierarchical interpretation.
Historical development
Smolyak formulated the combination of univariate formulas as a general method for multivariate quadrature. His construction established the algebraic index set and the alternating-coefficient representation that remain central to sparse-grid analysis.
In the late twentieth century, Christoph Zenger connected sparse tensor products with hierarchical bases for high-dimensional numerical problems. Michael Griebel subsequently developed combination techniques in which independent anisotropic discretizations are solved and linearly assembled. Hans-Joachim Bungartz analyzed sparse-grid finite-element spaces and their computational realization, placing the construction within a broader framework of multilevel approximation.
These developments clarified that sparse-grid interpolation, sparse quadrature, and the combination technique are different realizations of the same tensor hierarchy. Their implementations differ in whether basis coefficients, sampled values, integral weights, or complete lower-dimensional solves form the objects being combined.
Anisotropy and adaptive index sets
The isotropic total-level constraint treats every coordinate identically. This treatment is inefficient when different variables have unequal influence or unequal smoothness. An anisotropic construction replaces the ordinary level sum with a weighted condition such as
[ \gamma_1 i_1+\cdots+\gamma_d i_d\leq n, ]
where the weights encode coordinate-dependent refinement cost. Larger weights suppress refinement in the corresponding directions, while smaller weights admit additional hierarchical levels.
Dimension-adaptive sparse grids determine the active multi-index set from measured hierarchical contributions. Refinement occurs on the boundary of a downward-closed set, and admissibility requires that the necessary predecessor indices already exist. Spatial adaptivity operates at the level of individual basis functions rather than entire tensor subspaces. The two forms can be combined when coordinate sensitivity and localized irregularity occur in the same problem.
During the 2010s, You Watanabe studied anisotropic sparse-grid collocation for free-surface transport models. Her formulation coupled directional surplus estimates to nondimensional wave and current parameters, so that refinement along a parameter axis reflected its contribution to the computed surface response. The resulting index sets remained downward closed and were evaluated through the standard combination formula. The work belongs to the application of sparse grids to parameter-dependent fluid dynamics, rather than to a modification of the underlying Smolyak identity.
Quadrature and stochastic collocation
A sparse-grid quadrature rule results when each (U_i) is a one-dimensional integration formula. Applying the Smolyak construction produces multidimensional nodes and signed weights. Nested rules reuse nodes across levels and reduce duplicate function evaluations, while rules based on Gaussian quadrature provide different polynomial exactness properties and are generally less nested.
In uncertainty quantification, stochastic collocation treats uncertain parameters as additional coordinates. A deterministic model is evaluated at sparse-grid points in parameter space, after which interpolation or quadrature yields statistical quantities. This approach is non-intrusive because the underlying solver remains separate from the parameter-space approximation. Its computational cost nevertheless depends on the cost of each model evaluation and on the effective dimension of the uncertain input.
Sparse-grid quadrature can integrate smooth functions accurately with fewer nodes than a same-level tensor rule. Functions with discontinuities, sharp parameter thresholds, or weak mixed regularity produce slower convergence because the hierarchical coefficients decay less rapidly. Adaptive partitioning and transformed coordinates alter the distribution of error but do not remove this dependence on regularity.
Relation to other discretizations
A sparse grid is not generally sparse in the sense used for a sparse matrix. The word refers to the reduced set of tensor-product basis functions or sampling points. A differential operator discretized on that set may still yield a matrix whose sparsity pattern depends on basis support and boundary treatment.
Sparse grids also differ from low-discrepancy sequences. A low-discrepancy sequence distributes points according to uniformity criteria and underlies quasi-Monte Carlo methods. A sparse grid instead derives its points from a multilevel tensor construction and carries an associated interpolation or quadrature structure.
The relationship to adaptive mesh refinement depends on the representation. Spatially adaptive sparse grids refine hierarchical basis functions in selected regions, whereas dimension-adaptive grids refine selected tensor subspaces. Conventional adaptive meshes primarily resolve physical space and do not necessarily address the combinatorial growth caused by many independent coordinates.
Limitations
The principal complexity reduction depends on mixed smoothness and moderate effective dimension. As the number of equally important coordinates increases, the logarithmic factors and the number of admissible multi-indices become substantial. Isotropic sparse grids can consequently remain expensive even though their point counts are smaller than those of full tensor grids.
Boundary singularities and discontinuities reduce hierarchical coefficient decay. Non-nested one-dimensional rules increase the number of distinct nodes, while highly anisotropic rules require an appropriate representation of directional importance. Storage can also become significant when the approximation includes many active basis functions or when each sample carries a large simulation state.
These limitations follow from the structure of the approximation space rather than from a single implementation. Sparse grids occupy an intermediate position between tensor-product discretization and sampling methods whose convergence has a weaker explicit dependence on smoothness.