Bernstein Splines: Efficient Implementation

• Bernstein splines are defined by recasting local polynomial segments into the Bernstein basis, ensuring efficient enforcement of continuity and differentiability.
• The method leverages a three-stage process for coefficient computation and fast geometric evaluation, optimizing performance in isogeometric analysis, optimal control, and PDE discretization.
• Applications include trajectory optimization, multivariate geometric modeling, and advanced PDE solvers, all while maintaining numerical stability through convex recurrences.

A Bernstein splines-based implementation exploits the connection between B-spline basis functions and Bernstein polynomials to achieve computational efficiencies, algorithmic clarity, and, in some cases, improved numerical behavior for a wide array of applications including isogeometric analysis, optimal control, nonlinear PDE discretization, geometric modeling, and data interpolation. This approach systematically expresses splines—whether univariate or multivariate, polynomial or Tchebycheffian, uniform or non-uniform, compactly supported or globally continuous—by recasting all local polynomial pieces on each knot span or mesh cell in the Bernstein (or barycentric Bernstein) basis, then enforcing continuity, derivative, or variational conditions using local-to-global algebraic transforms.

1. Foundations of Bernstein-B-spline Equivalence

B-spline basis functions of degree $p$ on a knot vector $\tau=[\tau_0,\ldots,\tau_{m}]$ can always be written in terms of Bernstein polynomials on each non-empty knot interval $[\tau_j,\tau_{j+1})$ via an affine mapping. For a B-spline with clamped knots, each basis function $N_{i,p}(t)$ reduces to

$$N_{i,p}(t) = \sum_{k=0}^{p} b^{(i,j)}_k B^p_k\left( \frac{t-\tau_j}{\tau_{j+1}-\tau_j} \right)$$

for $t\in[\tau_j,\tau_{j+1})$, with $B^p_k(s) = \binom{p}{k} s^k (1-s)^{p-k}$. These coefficients $b^{(i,j)}_k$ encode all locality and continuity properties and can be constructed efficiently using recurrence relations (Chudy et al., 2024, Chudy et al., 2022).

In the special case of a single interval and uniform knots, $N_{i,p}$ itself is the Bernstein polynomial $B_{i,p}$. For multivariate polynomials over simplicial domains, the barycentric Bernstein basis is uniquely suited. This foundation underpins subsequent algorithmic advantages.

2. Core Algorithms for Bernstein Coefficient Computation

Modern Bernstein-spline implementations leverage a three-stage process to compute all local Bernstein-Bézier coefficients for each relevant span (Chudy et al., 2024, Chudy et al., 2022):

• Phase A: Diagonal Initialization Compute all $b_{p,p}^{(i,j)}$ recursively along the diagonal, exploiting closed-form solutions at interval endpoints.
• Phase B: Boundary Rows Directly assign coefficients at the leftmost ($k=0$) and rightmost ($k=p$) locations via explicit formulas using knot differences and product terms.
• Phase C: Interior Fill by Recurrence Use a main recurrence relation—for example,

$$b_{m,k}^{(i,j)} = \frac{t_j-t_i}{t_{j+1}-t_i}b_{m,k+1}^{(i,j)} + \frac{v_{m,i}}{t_{j+1}-t_i}\left[ (t_{j+1}-t_{m+i+2})b_{m,k}^{(i+1,j)} + (t_{m+i+2}-t_j)b_{m,k+1}^{(i+1,j)} \right],$$

with $v_{m,i} = (t_{m+i+1} - t_i) / (t_{m+i+2}-t_{i+1})$, to fill interior coefficients efficiently.

Total complexity per span is $O(p^2)$. For all $n_e$ non-empty spans, cost is $O(n_e p^2)$, which is asymptotically optimal for coefficient generation (Chudy et al., 2024, Chudy et al., 2022).

3. Geometric Evaluation and Implementation Advantages

Given Bernstein coefficients, fast geometric algorithms—especially de Casteljau-style evaluation—can be used to compute spline values (Chudy et al., 2022). For each evaluation point:

• Identify the containing knot span via binary search.
• Compute the local affine parameter $t$.
• Evaluate the sum $\sum_k b_k B^p_k(t)$ for each affected B-spline basis function via a linear-time de Casteljau process.
• Aggregate the weighted sum across active B-splines.

For tensor-product B-spline surfaces, the approach generalizes naturally. This method enables amortization: if coefficients are shared among multiple curves (e.g., in isogeometric surface evaluation or multi-curve rendering), the Bernstein evaluation approach can be significantly faster than repeated application of the de Boor-Cox algorithm, offering empirical speedups of $30\textrm{–}50$\% for curves and up to $4\textrm{–}6\times$ for surfaces when $M\gg1$ or $N\gg1$ (Chudy et al., 2022).

4. Application Scenarios

4.1 Trajectory Optimization and MPC

In flatness-based MPC, the Bernstein-spline parameterization yields a compact nonlinear program. Differentially-flat outputs are parameterized as

$y(t) = \sum_{i=0}^n c_i N_{i,p}(t),$

with the knot vector clamped to enforce boundary constraints. This reduces the number of variables, yields sparse banded Hessians/Jacobians due to local support, and empirically accelerates MPC solve times by up to $65\%$, reducing average solves from $20\textrm{ms}$ to $7\textrm{ms}$ (Neve et al., 2023).

4.2 Spline-based Discretization in PDEs and Fractional Calculus

For Hilfer fractional differential equations, the Bernstein-spline collocation at local Bernstein nodes allows analytical computation of weighted integrals, efficient assembly of the discretized nonlinear systems, and convergence rates of $O(q^{-\alpha/2}) + O(h^\alpha)$, with high accuracy at singularities and robustness for stiff nonlinearities (Goedegebure et al., 28 Mar 2025).

4.3 Multivariate and Simplicial Meshes

The barycentric Bernstein form is the foundation for explicit $C^r$ B-spline bases on simplicial meshes. By expressing each local polynomial in the Bernstein basis and assembling global continuity conditions via sparse linear systems enforcing value and derivative matching on shared faces, one constructs global basis functions with the desired smoothness (Linger et al., 2014).

4.4 Tchebycheffian and Multi-degree Generalizations

For splines with non-polynomial pieces—Tchebycheffian, mixed, or exponential—the Bernstein-like basis (via generalized powers/interpolants) paired with extraction matrices enable construction of spline spaces with variable local degree and continuity. Internal block-diagonal transforms (B-extraction operators) enforce Sobolev-type smoothness at break points and preserve locality, minimal support, and positivity properties (Hiemstra et al., 2020, Speleers, 2021).

5. Numerical Stability and Practicalities

Stability is achieved through recurrences based on convex combinations; each step involves only knot-difference multipliers and guarantees that error accumulation is at most $O(p\epsilon_\text{machine})$ for floating-point arithmetic. For high-degree splines (up to $p\approx 50$), relative errors remain at machine precision and there is no significant amplification compared to older $O(p^3)$ schemes. Memory requirements are minimal; each span holds an $O(p^2)$ 2D array, suitable for modern vectorized, GPU, or multicore implementation (Chudy et al., 2024, Chudy et al., 2022).

6. Extensions: Surfaces, Manifold Splines, and Non-Polynomial Bases

Subdivision-based Bernstein algorithms extend to surfaces and even to splines on general manifolds, as in the recursive geodesic-average methods for defining Bézier curves on triangle meshes. Subdivision ensures robust, globally defined and smooth limiting curves independent of mesh geometry. Open-uniform Lane–Riesenfeld subdivision and recursive De Casteljau-style bisection methods deliver, respectively, $C^{k-1}$ and $C^1$ continuity for arbitrary control polygons with high robustness and mesh independence (Mancinelli et al., 2021).

For explicit cubic quasi-interpolating splines on regular triangulations, Bernstein-Bézier coefficients are assembled directly from point and gradient values at vertices, with global $O(h^3)$ convergence and $O(h^4)$ superconvergent error at edge midpoints (Barrera et al., 2024).

7. Summary Table: Implementation Steps and Complexity

Step / Component	Algorithmic Cost	Salient Features
Bernstein-Bézier coefficient gen.	$O(p^2)$ per span	Convex recurrence, optimal for any knot sequence
Evaluation (per point/curve)	$O(p^2)$	Amortizable over many curves with shared knots
Multivariate/mesh bases	$O(DN)$, $N\sim d^m$	Linear system assembly for continuity constraints
Spline extraction (Tchebycheffian)	$O(N^2)$	Row and local block null-space operations
PDE/IVP discretization	$O((k+1)q^2)$	Per subinterval, nonlinear fixed-point (collocation)

This approach unifies a wide range of spline technologies: from isogeometric analysis and control theory to spline-PDE solvers and computational geometry, providing a scalable, robust, and explicit computational foundation grounded in the structure of the Bernstein basis (Chudy et al., 2024, Chudy et al., 2022, Neve et al., 2023, Linger et al., 2014, Hiemstra et al., 2020, Speleers, 2021, Goedegebure et al., 28 Mar 2025, Barrera et al., 2024).