Gradient Hermite–Birkhoff Interpolation

• Gradient Hermite–Birkhoff interpolation is defined by constraining approximants with function values and derivative data (e.g., gradients, partial derivatives) at prescribed nodes.
• Recursive Newton-type and barycentric formulations enhance computational efficiency and stability across univariate, multivariate, and spectral implementations.
• The approach is versatile, finding applications in high-dimensional sparse recovery, PDE numerical schemes, and interpolation on manifolds and grids.

Gradient Hermite–Birkhoff interpolation is the interpolation problem in which an approximant is constrained by function values together with derivative data—most commonly gradients, partial derivatives, or linear combinations of derivatives—at prescribed nodes. In the univariate Birkhoff formulation, the data are linear functionals of the form $L_i=\delta_{x_{\beta_i}}\circ D_x^{\alpha_i}$, or more generally $\delta_{x_\beta}\circ P(D_x)$ with $P(D_x)=\sum_k c_kD_x^k$; in several variables, the natural analogue is $$L_{i,\alpha}(p)=\delta_{z_i}\circ\partial^\alpha(p)=\partial^\alpha p(z_i)$$, so gradient components are first-order instances of the same framework (Jiang et al., 12 Nov 2025). Recent work develops this subject in several distinct directions: recursive polynomial constructions, barycentric realizations, Hermite spectral methods on $\mathbb{R}$, sparse high-dimensional recovery from gradient-augmented samples, and multivariate formulations on manifolds, rectilinear grids, and subdivision hierarchies (Wang et al., 21 Jul 2025).

1. Functional formulation and interpolation structure

Classical Hermite interpolation prescribes contiguous derivatives at each node, whereas Birkhoff interpolation permits arbitrary derivative orders and arbitrary nodewise patterns. In the univariate setting, the derivative pattern is encoded by an incidence matrix $E=(e_{i,j})$ with $e_{i,j}\in\{0,1\}$, and the interpolation problem is to find a polynomial subspace $\mathcal F\subset\mathbb F[x]$ such that there exists a unique $p\in\mathcal F$ satisfying

$$p^{(\alpha)}(x_\beta)=y_{\beta,\alpha},\qquad (\beta,\alpha)\in\tilde e.$$

A proper interpolation space is one with $\delta_{x_\beta}\circ P(D_x)$0, where $\delta_{x_\beta}\circ P(D_x)$1 is the number of conditions, and uniqueness holds (Jiang et al., 12 Nov 2025).

The same framework covers gradient data. Each condition is a composition of an evaluation functional and a differential operator,

$\delta_{x_\beta}\circ P(D_x)$2

and the generalized algorithmic setting also allows

$\delta_{x_\beta}\circ P(D_x)$3

In several variables, the formal replacement $\delta_{x_\beta}\circ P(D_x)$4 yields

$\delta_{x_\beta}\circ P(D_x)$5

so a full gradient at $\delta_{x_\beta}\circ P(D_x)$6 appears as the block $\delta_{x_\beta}\circ P(D_x)$7. This identifies gradient Hermite–Birkhoff interpolation as a special case of linear-functional interpolation in which the prescribed functionals are first-order partial derivatives or their linear combinations (Jiang et al., 12 Nov 2025).

A key structural notion is the Newton-type, or triangular, basis. Given ordered functionals $\delta_{x_\beta}\circ P(D_x)$8, a sequence $\delta_{x_\beta}\circ P(D_x)$9 is triangular if $P(D_x)=\sum_k c_kD_x^k$0 when $P(D_x)=\sum_k c_kD_x^k$1 and $P(D_x)=\sum_k c_kD_x^k$2 when $P(D_x)=\sum_k c_kD_x^k$3. The existence of such a sequence is equivalent to linear independence of the functionals. In interpolation language, this is the basis-theoretic form of well-posedness, and it is especially useful when gradient and higher-derivative constraints are added incrementally (Jiang et al., 12 Nov 2025).

2. Recursive and barycentric realizations

For polynomial Birkhoff problems, recent recursive algorithms replace direct Gaussian elimination with a Schur-complement and Sylvester-identity construction. If $P(D_x)=\sum_k c_kD_x^k$4 denotes the interpolant satisfying the first $P(D_x)=\sum_k c_kD_x^k$5 conditions and $P(D_x)=\sum_k c_kD_x^k$6 is the next Newton-type basis polynomial, the update has the Newton form

$P(D_x)=\sum_k c_kD_x^k$7

The corresponding basis update is likewise recursive. The decisive algorithmic test is the judgment condition $P(D_x)=\sum_k c_kD_x^k$8; if it fails, the functional order can be pivoted or the candidate degree can be increased. This yields a strongly proper basis, meaning that every leading subsystem remains unisolvent as conditions are added (Jiang et al., 12 Nov 2025).

The same paper emphasizes that this recursion generalizes the Generalized Recursive Polynomial Interpolation Algorithm from Hermite to broader Birkhoff data and computes a lower-degree Newton-type interpolation basis that is also a strongly proper interpolation basis. Compared with existing algorithms that rely on Gaussian elimination to compute the interpolation basis, the recursive approach significantly reduces the computational cost (Jiang et al., 12 Nov 2025).

A complementary implementation route is barycentric Hermite interpolation. For higher-order Hermite–Fejér interpolation at Gauss–Jacobi or Jacobi–Gauss–Lobatto pointsystems, the second barycentric formula admits a fast implementation in which the barycentric weights can be efficiently evaluated and cost linear operations corresponding to the number of grids totally. The crucial numerical point is that the exponentially increasing common factor in the barycentric weights is canceled by the division in the second barycentric form, yielding a superiorly stable method for computing the simplified barycentric weights and leading to a fast implementation with linear operations on the number of grids. Convergence rates are also derived for Hermite–Fejér interpolation at Gauss–Jacobi pointsystems (Xiang et al., 2014).

In the context of gradient Hermite–Birkhoff interpolation, these two realizations serve different regimes. Recursive Newton-type constructions are naturally adapted to arbitrary derivative patterns, while barycentric formulas are most effective when the nodes are structured and the multiplicities are low, as in value-plus-gradient or low-order Hermite conditions.

3. Hermite spectral interpolation and derivative superconvergence on $P(D_x)=\sum_k c_kD_x^k$9

On the unbounded domain $$L_{i,\alpha}(p)=\delta_{z_i}\circ\partial^\alpha(p)=\partial^\alpha p(z_i)$$0, Hermite spectral interpolation uses Hermite functions

$$L_{i,\alpha}(p)=\delta_{z_i}\circ\partial^\alpha(p)=\partial^\alpha p(z_i)$$1

as the approximation space and zeros of $$L_{i,\alpha}(p)=\delta_{z_i}\circ\partial^\alpha(p)=\partial^\alpha p(z_i)$$2 as interpolation nodes. The interpolant $$L_{i,\alpha}(p)=\delta_{z_i}\circ\partial^\alpha(p)=\partial^\alpha p(z_i)$$3 is defined by $$L_{i,\alpha}(p)=\delta_{z_i}\circ\partial^\alpha(p)=\partial^\alpha p(z_i)$$4, where $$L_{i,\alpha}(p)=\delta_{z_i}\circ\partial^\alpha(p)=\partial^\alpha p(z_i)$$5 are the zeros of $$L_{i,\alpha}(p)=\delta_{z_i}\circ\partial^\alpha(p)=\partial^\alpha p(z_i)$$6. For analytic $$L_{i,\alpha}(p)=\delta_{z_i}\circ\partial^\alpha(p)=\partial^\alpha p(z_i)$$7 in a strip $$L_{i,\alpha}(p)=\delta_{z_i}\circ\partial^\alpha(p)=\partial^\alpha p(z_i)$$8, the interpolation error admits a contour representation, and differentiation yields global derivative estimates of the form

$$L_{i,\alpha}(p)=\delta_{z_i}\circ\partial^\alpha(p)=\partial^\alpha p(z_i)$$9

Thus the first- and second-derivative errors are globally $\mathbb{R}$0 and $\mathbb{R}$1, respectively (Wang et al., 21 Jul 2025).

The distinctive contribution is the identification of derivative superconvergence points. For the first derivative, the special points are the zeros $\mathbb{R}$2 of $\mathbb{R}$3; for the second derivative, the special points are the zeros $\mathbb{R}$4 of $\mathbb{R}$5. At these points the dominant term in the differentiated contour integral vanishes. The resulting bounds are smaller than the global derivative bounds, with a factor $\mathbb{R}$6 confirmed numerically as sharp for both first and second derivatives. In the first-derivative case, the analysis distinguishes a general improvement and a sharper central-region improvement near the turning point, but the experiments show that the effective gain in practice is the factor $\mathbb{R}$7 (Wang et al., 21 Jul 2025).

This has direct consequences for gradient Hermite–Birkhoff design on $\mathbb{R}$8. The paper itself treats function-value interpolation, not a Birkhoff problem, but it identifies where gradients are intrinsically most accurate. This suggests the following node geometry for Hermite-function schemes on $\mathbb{R}$9: function-value conditions at zeros of $E=(e_{i,j})$0 or $E=(e_{i,j})$1, first-derivative conditions at zeros of $E=(e_{i,j})$2, and second-derivative conditions at zeros of $E=(e_{i,j})$3. The same pattern reappears in Hermite spectral collocation for second-order ODEs, where the collocation solution superconverges at zeros of $E=(e_{i,j})$4 and its first and second derivatives superconverge at zeros of $E=(e_{i,j})$5 and $E=(e_{i,j})$6, respectively (Wang et al., 21 Jul 2025).

4. Multivariate constructions on manifolds, algebraic sets, and grids

On smooth manifolds, Hermite–Birkhoff interpolation can be formulated in local coordinates by prescribing

$E=(e_{i,j})$7

For gradient interpolation, one takes $E=(e_{i,j})$8, so the local data are $E=(e_{i,j})$9 and the coordinate first derivatives. The interpolant is written as

$e_{i,j}\in\{0,1\}$0

where $e_{i,j}\in\{0,1\}$1 is an incomplete Taylor expansion and $e_{i,j}\in\{0,1\}$2 are cardinal basis functions built from the geodesic distance. The basis functions satisfy a partition of unity, cardinality $e_{i,j}\in\{0,1\}$3, and vanishing derivative conditions at the nodes, so the method imposes Hermite–Birkhoff data without solving linear systems. For localized cardinal functions, the error satisfies an estimate of order $e_{i,j}\in\{0,1\}$4 when complete Taylor data up to order $e_{i,j}\in\{0,1\}$5 are used (Allasia et al., 2017).

A different multivariate route is algebraic. If the nodes are the common zeros of a polynomial map $e_{i,j}\in\{0,1\}$6, then local algebras and Grothendieck residues provide explicit multidimensional Hermite formulas. The interpolation conditions are

$e_{i,j}\in\{0,1\}$7

where $e_{i,j}\in\{0,1\}$8 is the monomial basis of the local algebra at $e_{i,j}\in\{0,1\}$9. Gradient interpolation is the special choice $\mathcal F\subset\mathbb F[x]$0. The explicit residue formula applies when the local ideal is generated by pure powers $\mathcal F\subset\mathbb F[x]$1, so the multiplicity diagram is a parallelepiped. In that class, the paper constructs basis polynomials dual to derivative functionals and gives a concrete two-dimensional example combining value, first derivatives, and a mixed second derivative at a multiple root (Durakov et al., 2022).

On non-equally spaced rectilinear grids, classical multivariate coordinate Hermite interpolation provides a complete-data baseline. For grid $\mathcal F\subset\mathbb F[x]$2 with multiplicity vector $\mathcal F\subset\mathbb F[x]$3, there exists a unique interpolating polynomial in

$\mathcal F\subset\mathbb F[x]$4

matching all partial derivatives $\mathcal F\subset\mathbb F[x]$5 for $\mathcal F\subset\mathbb F[x]$6. The interpolant has a compact single-summation closed form using local lower-unitriangular matrices $\mathcal F\subset\mathbb F[x]$7, and the interpolation ideal is $\mathcal F\subset\mathbb F[x]$8, where $\mathcal F\subset\mathbb F[x]$9. The paper also proves continuity of Hermite polynomials defined on adjacent $p\in\mathcal F$0-dimensional grids, establishing spline behavior (Kechriniotis et al., 2023). Gradient Hermite–Birkhoff interpolation differs from this classical setting by deleting some of the derivative conditions inside the full rectangular block, but the coordinate Hermite theory supplies the algebraic and spline-theoretic reference model.

5. Sparse, high-dimensional, and multiresolution formulations

In high dimensions, one influential extension is compressive Hermite interpolation. Here the data at each sampling point $p\in\mathcal F$1 consist of the function value and all first-order partial derivatives, and the goal is sparse polynomial approximation in a tensor-product basis via weighted $p\in\mathcal F$2-minimization. The main theoretical result is that, for the same asymptotic sample complexity as the function-only setting, gradient-augmented measurements yield an approximation error bound in a stronger Sobolev norm rather than merely the $p\in\mathcal F$3-norm. For Chebyshev and Legendre polynomial approximations, the sample complexity estimate is algebraic in the sparsity $p\in\mathcal F$4 and at most logarithmic in the dimension $p\in\mathcal F$5, mitigating the curse of dimensionality to a substantial extent (Adcock et al., 2017).

This framework is not exact interpolation in the classical finite-dimensional sense; it is randomized recovery in a sparse polynomial model. Nevertheless, it is Hermite–Birkhoff in its data model: each sample contributes function and gradient conditions, and the linear measurement matrix contains rows built from basis values and basis partial derivatives. In noiseless regimes and on exact supports, the recovered polynomial behaves as an interpolant; with truncation or noise, it is more precisely a gradient-augmented approximation (Adcock et al., 2017).

A multiresolution counterpart is provided by multivariate generalized Hermite subdivision schemes. One chooses an ordered multiset $p\in\mathcal F$6, where each $p\in\mathcal F$7 specifies which derivative is carried in the $p\in\mathcal F$8-th data channel. Gradient data correspond to including the first-order multi-indices $p\in\mathcal F$9. The refinement uses a matrix mask $$p^{(\alpha)}(x_\beta)=y_{\beta,\alpha},\qquad (\beta,\alpha)\in\tilde e.$$0 and the derivative-scaling matrix

$$p^{(\alpha)}(x_\beta)=y_{\beta,\alpha},\qquad (\beta,\alpha)\in\tilde e.$$1

so the iterates encode dyadically refined values of the selected derivatives. The interpolation property of the associated basis vector function $$p^{(\alpha)}(x_\beta)=y_{\beta,\alpha},\qquad (\beta,\alpha)\in\tilde e.$$2 is

$$p^{(\alpha)}(x_\beta)=y_{\beta,\alpha},\qquad (\beta,\alpha)\in\tilde e.$$3

which is the subdivision analogue of a Hermite–Birkhoff cardinal basis. The paper also introduces linear-phase moments as the condition guaranteeing polynomial interpolation and proves constructively that, for any given positive integer $$p^{(\alpha)}(x_\beta)=y_{\beta,\alpha},\qquad (\beta,\alpha)\in\tilde e.$$4, there always exist convergent smooth generalized Hermite subdivision schemes with linear-phase moments such that their basis vector functions are spline functions in $$p^{(\alpha)}(x_\beta)=y_{\beta,\alpha},\qquad (\beta,\alpha)\in\tilde e.$$5 and have linearly independent integer shifts (Han, 2021).

6. Regularity, singularity, and PDE-oriented usage

A persistent issue in multivariate Hermite and Hermite–Birkhoff interpolation is regularity. Dimension counting alone does not guarantee unisolvence. For multivariate Hermite interpolation of type total degree, the singularity theory is explicit: all Hermite interpolation of type total degree on $$p^{(\alpha)}(x_\beta)=y_{\beta,\alpha},\qquad (\beta,\alpha)\in\tilde e.$$6 points in $$p^{(\alpha)}(x_\beta)=y_{\beta,\alpha},\qquad (\beta,\alpha)\in\tilde e.$$7 is singular if $$p^{(\alpha)}(x_\beta)=y_{\beta,\alpha},\qquad (\beta,\alpha)\in\tilde e.$$8, and all Hermite interpolations of type total degree on $$p^{(\alpha)}(x_\beta)=y_{\beta,\alpha},\qquad (\beta,\alpha)\in\tilde e.$$9 points with $\delta_{x_\beta}\circ P(D_x)$00 are singular. For $\delta_{x_\beta}\circ P(D_x)$01 and $\delta_{x_\beta}\circ P(D_x)$02, only three cases and one case can produce regular interpolation schemes, respectively (Meng et al., 2013).

These results are especially relevant to gradient data. In the paper’s notation, prescribing function values and all first derivatives at node $\delta_{x_\beta}\circ P(D_x)$03 is the case $\delta_{x_\beta}\circ P(D_x)$04. Among small-node regimes, the five-node, three-dimensional case $\delta_{x_\beta}\circ P(D_x)$05, $\delta_{x_\beta}\circ P(D_x)$06, $\delta_{x_\beta}\circ P(D_x)$07, $\delta_{x_\beta}\circ P(D_x)$08 is identified as an almost regular scheme, while many nearby total-degree configurations are singular (Meng et al., 2013). The conclusion is that gradient Hermite–Birkhoff interpolation in several variables must be designed with explicit attention to the interpolation space, not merely to the number of scalar conditions.

A PDE-oriented use of gradient Hermite interpolation appears in numerical schemes for heat equations. A collocation polynomial is first used to approximate spatial derivatives, producing numerical gradients such as

$\delta_{x_\beta}\circ P(D_x)$09

and then cubic Hermite interpolation reconstructs midpoint values from nodal values and gradients. The resulting numerical gradient scheme has convergence order $\delta_{x_\beta}\circ P(D_x)$10 under the discrete maximum norm; moreover, the same order is retained when the space step size is just twice the one of the high-order compact difference method, which accelerates the computational process and makes the result much smoother to some extent (Li et al., 2013). In two dimensions, analogous formulas use nodal values together with $\delta_{x_\beta}\circ P(D_x)$11- and $\delta_{x_\beta}\circ P(D_x)$12-derivative approximants to reconstruct edge and cell-center values (Li et al., 2013).

Taken together, these developments show that gradient Hermite–Birkhoff interpolation is less a single algorithm than a family of interpolation principles. Its common core is the use of derivative-valued linear functionals—exact gradients, partial derivatives, or operator-derived surrogates—to constrain an approximant. The main design variables are the interpolation space, the node geometry, and the admissible derivative pattern. On $\delta_{x_\beta}\circ P(D_x)$13, Hermite spectral theory singles out superconvergent derivative nodes; on manifolds and grids, Taylor–cardinal and coordinate-Hermite constructions provide explicit bases; in sparse and subdivision settings, gradient data are embedded into recovery and refinement operators; and in multivariate polynomial spaces, singularity theory sharply limits which gradient data configurations can be unisolvent.