Numeric Stability and Speed of B-Splines

Four approaches to the computation of the B-spline \(\beta,\) and discussion of their relative merits in terms of speed and numerical accuracy.

Splines

Splines are piecewise polynomials that are built as the weighted sum of integer-shifted B-splines, themselves being piecewise polynomials, too. The weights are called the spline coefficients. Formally, a spline of degree \(n\in{\mathbb{N}}\) is the mapping

\((1)\)

\[f:{\mathbb{R}}\rightarrow{\mathbb{R}},x\mapsto f(x)=\sum_{k\in{\mathbb{Z}}}\,c[k]\,\beta^{n}(x-k),\]

where \(c\) are the spline coefficients and \(\beta^{n}\) the B-spline, with \(n\) a superscript (not a power). To accurately evaluate the spline \(f\) at the argument \(x,\) it is crucial that the computation of \(\beta^{n}\) be stable numerically. This is easy to achieve for small degrees; when the degree rises, however, one has to face the fact that the involved polynomials have a tendency to be made of delicately balancing terms.

Here, we investigate how the numerical stability depends on the way B-splines are computed. The splinekit library proposes a remarquably stable strategy that results in the very fast and constant-speed evaluation of B-splines.

B-Splines

A B-spline \(\beta^{n}\) can be computed in more than one way. Here are three possibilities.

Classic Formula

The classic, productive representation of the \(m\)-th derivative of a B-spline of degree \(n\in{\mathbb{N}}\) and argument \(x\in{\mathbb{R}}\) is

\((2)\)

\[\frac{{\mathrm{d}}^{m}\beta^{n}(x)}{{\mathrm{d}}x^{m}}=\sum_{k=0}^{n+1}\,{\color{blue}\left(-1\right)^{k}}\,{n+1\choose k}\,\varsigma^{n-m}(x+\frac{n+1}{2}-k).\]

(To compute a non-differentiated B-spline, simply set \(m=0.\)) Unfortunately, the term \({\color{blue}\left(-1\right)^{k}}\) results in additive contributions that tend to cancel each other. This spells numerical trouble, even if the polynomial simple element \(\varsigma^{n}:{\mathbb{R}}\rightarrow{\mathbb{R}},x\mapsto\varsigma^{n}(x)=\frac{1}{2\,n!}\,{\mathrm{sgn}}(x)\,x^{n}\) has the flavor of a well-behaved canonic monomial. Another issue lies with the growth of the binomial coefficients with respect to the degree, which ultimately leads to a delicate balance between large numbers.

De Boor

Other computational recipes have been devised. For instance, the De Boor’s recurrence relation expresses a B-spline of some degree as a weighted combination of B-splines of lesser degree. As defined in \((2),\) it turns out that \(\beta^{0}(x)=\left(\varsigma^{0}(x+\frac{1}{2})-\varsigma^{0}(x-\frac{1}{2})\right),\) which is the base case. Then, the recursive step for \(n\in{\mathbb{N}}+1\) is such that

\((3)\)

\[\beta^{n}(x)=\frac{1}{n}\,\left(\left(x+\frac{n+1}{2}\right)\,\beta^{n-1}(x+\frac{1}{2})-\left(x-\frac{n+1}{2}\right)\,\beta^{n-1}(x-\frac{1}{2})\right).\]

splinekit

In the splinekit library, B-splines of degree \(n=0\) are computed as in the classic or the De Boor’s approach. B-splines of positive degree, however, are computed as the scalar product

\((4)\)

\[\beta^{n}(x)={\mathbf{[\![}}\left|x\right|<\frac{n+1}{2}\,{\mathbf{]\!]}}\,\left(w^{n}[r][\cdot]\right)^{{\mathsf{T}}}\,{\mathbf{v}}^{n}(\chi),\]

where the notation \({\mathbf{[\![}}\cdot\,{\mathbf{]\!]}}\) is that of the Iverson bracket. Moreover, \(r=\left\lceil\xi\right\rceil\in{\mathbb{Z}}\) for \(\xi=\left(\frac{n-1}{2}-x\right),\) and \(\chi=\left(r-\xi\right)\in[0,1).\) Finally, \(\left(w^{n}[r][\cdot]\right)^{{\mathsf{T}}}\) is the \((r+1)\)-th row of the B-spline evaluation matrix \({\mathbf{W}}^{n}\in{\mathbb{Q}}^{\left(n+1\right)\times\left(n+1\right)}\) and \({\mathbf{v}}^{n}(\chi)\in{\mathbb{R}}^{n+1}\) is the Vandermonde vector of argument \(\chi\) and degree \(n.\) In this formulation of a B-spline, the fact that the Vandermonde vector has the domain \([0,1)\) greatly favors numerical stability since the image of each of its components is \([0,1].\)

The matrix \({\mathbf{W}}^{n}\) depends on the degree only and can be precomputed and cached. Its component at the \((r+1)\)-th row and \((c+1)\)-th column is rational and is defined as

\((5)\)

\[w_{r+1,c+1}^{n}=w^{n}[r][c]=\frac{1}{c!}\,\left(\left.\frac{{\mathrm{d}}^{c}\beta^{n}(x)}{{\mathrm{d}}x^{c}}\right|_{x=\frac{n-1}{2}-r}+\frac{1}{2}\,{\mathbf{[\![}}c=n\,{\mathbf{]\!]}}\,\left(-1\right)^{n-r}\,{n+1\choose r+1}\right).\]

The B-spline derivatives that appear in \((5)\) are computed as rational numbers, through \((2).\) This leads to a rational value for \(w^{n}[r][c],\) which is then cached as its float approximation. With this strategy, the balancing act that was governed in \((2)\) by the term \({\color{blue}\left(-1\right)^{k}}\) is now applied to rational values and can be achieved exactly in the rational domain.

Experiments

We want to compare four computational approaches.

  1. In the first (ground-truth) approach, we rely on the observation that \(\beta^{n}:{\mathbb{R}}\rightarrow{\mathbb{R}}\) is defined in \((2)\) in such a way that \(\beta^{n}\) would map a rational number to a number that is rational, too, so that \(\beta^{n}:{\mathbb{Q}}\rightarrow{\mathbb{Q}}.\) This allows us to compute exact, ground-truth values by relying on Python’s built-in Fraction type.

  2. In the second (classic) approach, we compute \(\beta^{n}\) as in \((2),\) this time relying on Python’s built-in float type.

  3. In the third (De Boor) approach, we compute \(\beta^{n}\) as in \((3),\) again relying on Python’s built-in float type.

  4. In the fourth (splinekit) approach, we compute \(\beta^{n}\) as in \((4),\) relying on Python’s built-in float type.

What we measure is by how far the methods depart from the ground truth, in a mean-square sense expressed in dB units.

Danger

  • The link below allows you to access and run the installation-free notebook. While the reported SNRs are representative, the timings are not because the Python kernel is WebAssembly-based and does not run natively.

  • If you want to test for realistic timings on your own computer, then you will have to first install in full the splinekit library. Only after that will you be able to launch the notebook either as a regular, full-fledged Jupyter Lab or as a module executed by the native Python kernel.

  • The timings reported in the Results Sections correspond to those of the native execution.

Jupyter Lab notebook

Numeric stability and speed of B-splines

Hint

The notebook is available for download in compressed form from here. Decompression is achieved from the terminal with gunzip bspline_numeric_stability.ipynb.gz.

Critical Regime

When expressed as a signal-to-noise ratio (SNR), the error is most noticeable near the end of the left and right tails of the B-spline. In this critical regime, the error is computed over \(400\) samples taken at random, consistently between methods to allow for their fair comparison. We also report the computational time.

Results

On a desktop computer of 2021, a typical resulting table is as follows.

Accuracy and time in the critical regime
 Ground-Truth  Classic  De Boor  splinekit 
 Degree  SNR[dB]   Time[s]  SNR[dB]   Time[s]  SNR[dB]   Time[s]  SNR[dB]   Time[s] 
0 ∞  3.1e-03 ∞  2.6e-04 ∞  1.6e-04 ∞  1.2e-04 
1 ∞  7.2e-03 ∞  5.5e-04 ∞  3.2e-04 ∞  1.4e-03 
2 ∞  1.0e-02 309.4  6.7e-04 ∞  6.3e-04 327.2  1.5e-03 
3 ∞  1.4e-02 297.7  7.5e-04 320.9  1.3e-03 318.6  1.5e-03 
4 ∞  1.8e-02 282.3  8.7e-04 319.3  2.5e-03 313.0  1.3e-03 
5 ∞  2.2e-02 262.1  9.9e-04 317.5  5.0e-03 310.6  1.3e-03 
6 ∞  2.7e-02 243.6  1.2e-03 317.2  9.9e-03 303.3  1.3e-03 
7 ∞  3.1e-02 224.9  1.4e-03 314.9  2.0e-02 302.2  1.3e-03 
8 ∞  3.7e-02 199.3  1.4e-03 315.2  4.0e-02 295.9  1.3e-03 
9 ∞  4.3e-02 171.9  1.6e-03 315.3  7.9e-02 294.5  1.3e-03 
10 ∞  5.0e-02 143.8  1.7e-03 313.5  1.6e-01 290.8  1.3e-03 
11 ∞  5.7e-02 120.2  1.8e-03 314.2  3.2e-01 287.6  1.3e-03 
12 ∞  6.5e-02 94.4  2.0e-03 312.2  6.4e-01 281.7  1.3e-03 
13 ∞  7.4e-02 72.6  2.3e-03 310.9  1.3e+00 284.1  1.3e-03 
14 ∞  8.3e-02 45.4  2.4e-03 310.3  2.5e+00 277.3  1.3e-03 
15 ∞  9.2e-02 15.9  2.6e-03 310.0  5.1e+00 273.7  1.3e-03 
16 ∞  1.0e-01 -13.4  2.7e-03 312.6  1.0e+01 270.1  1.3e-03 

Discussion for the Critical Regime

  • The numerical stability of the classic approach is consistently the worst at all degrees. Moreover, it even collapses when the degree rises: for degree \(n=16\) already, there is more noise than signal, at least over the domain being examined, which is \((-8.5,-6.5)\cup(6.5,8.5)\) for \(n=16.\) The computational cost exhibits a linear increase over the degrees explored here. Yet, the classic approach is still the fastest over the goldilocks range of degrees \(n\in\{3,4,5,6\}\), a range where its accuracy can also be considered sufficient for most purposes.

  • The numerical stability of the De Boor’s approach is always the best (310dB corresponds to about 51 significant bits) but comes at a computational cost that increases combinatorially with the degree. Yet, the trivial code proposed here comes the fastest over the small range of degrees \(n\in\{0,1,2\}\). For such low degrees, however, one has to realize that the computational cost depends more on language idiosyncrasies (e.g., namespace lookup, recursivity vs loops, checks on type and validity of the parameters) than on algorithmic considerations.

  • The numerical stability of the splinekit library is consistently much higher than that of the classic approach and degrades more slowly with the degree (270dB corresponds to about \(45\) significant bits). Indeed, one has to reach the very high degree \(n=94\) before the \(0\) dB threshold is crossed, at least over the domain being examined, which is \((-47.5,-45.5)\cup(45.5,47.5)\) for \(n=94.\)

  • For degree \(n=16\), splinekit is more than seven thousand times faster than De Boor while, for degree \(n=7\) and higher, the splinekit strategy is always faster than the other three. The computation time of splinekit remains constant through all degrees, a property that holds up to \(n=94\) and beyond.

Global Regime, Low Degrees

In the global regime, the error is computed over \(400\) samples taken at random according to a uniform distribution that is not restricted to the tails of B-splines but that covers their whole support, consistently between methods to allow for their fair comparison.

Here, we consider only low degrees. This allows us to include results for the (slow) De Boor’s method.

Results

On a desktop computer of 2021, a typical resulting table is as follows.

Accuracy and time in the global regime, low degrees
 Ground-Truth  Classic  De Boor  splinekit 
 Degree  SNR[dB]   Time[s]  SNR[dB]   Time[s]  SNR[dB]   Time[s]  SNR[dB]   Time[s] 
0 ∞  4.2e-03 ∞  3.8e-04 ∞  1.6e-04 ∞  1.3e-04 
1 ∞  9.0e-03 ∞  5.4e-04 ∞  3.1e-04 ∞  1.4e-03 
2 ∞  1.1e-02 308.6  6.6e-04 322.4  6.3e-04 319.8  1.3e-03 
3 ∞  1.4e-02 299.6  7.8e-04 319.6  1.3e-03 319.1  1.3e-03 
4 ∞  1.8e-02 287.2  9.1e-04 316.7  2.5e-03 316.5  1.3e-03 
5 ∞  2.2e-02 276.6  1.0e-03 316.6  5.0e-03 314.7  1.3e-03 
6 ∞  2.7e-02 264.3  1.2e-03 316.0  9.9e-03 317.3  1.3e-03 
7 ∞  3.2e-02 258.0  1.3e-03 317.1  2.0e-02 320.7  1.3e-03 
8 ∞  3.8e-02 243.2  1.5e-03 317.1  3.9e-02 313.7  1.3e-03 
9 ∞  4.4e-02 230.4  1.6e-03 316.4  7.9e-02 314.8  1.3e-03 
10 ∞  5.1e-02 216.5  1.8e-03 315.9  1.6e-01 314.3  1.3e-03 
11 ∞  5.8e-02 212.0  1.9e-03 314.8  3.2e-01 313.3  1.2e-03 
12 ∞  6.6e-02 198.8  2.1e-03 315.5  6.3e-01 313.4  1.2e-03 
13 ∞  7.5e-02 192.8  2.3e-03 313.5  1.3e+00 313.4  1.3e-03 
14 ∞  8.4e-02 178.4  2.5e-03 314.4  2.5e+00 311.4  1.3e-03 
15 ∞  9.4e-02 172.7  2.6e-03 314.0  5.1e+00 317.0  1.3e-03 
16 ∞  1.0e-01 161.0  2.8e-03 312.9  1.0e+01 310.1  1.3e-03 

Discussion for the Global Regime, Low Degrees

  • The timings of all methods are consistent with the timings discussed in the critical regime. One concludes that the computational costs do not depend on the value of the arguments of B-splines.

  • The numerical stability of the classic approach is again consistently the worst at all degrees. Yet, contrarily to the critical regime, it remains globally serviceable for the degrees considered here.

  • The numerical stability of the De Boor’s approach is approximately the same in the critical regime or in the global regime.

  • The numerical stability of the splinekit library is on par with that of the De Boor’s approach and consistently much higher than that of the classic approach. It does not degrade with the degree.

Global Regime

In the global regime, the error is computed over \(400\) samples taken at random according to a uniform distribution that is not restricted to the tails of B-splines but that covers their whole support, consistently between methods to allow for their fair comparison.

Here, we consider degrees too high to allow the (slow) De Boor’s method to terminate computations in a reasonable time.

Results

On a desktop computer of 2021, a typical resulting table is as follows.

Accuracy and time in the global regime
 Ground-Truth  Classic  splinekit 
 Degree  SNR[dB]   Time[s]  SNR[dB]   Time[s]  SNR[dB]   Time[s] 
0 ∞  3.3e-03 ∞  3.8e-04 ∞  1.1e-04 
1 ∞  7.9e-03 ∞  5.6e-04 ∞  1.3e-03 
2 ∞  1.2e-02 307.2  6.8e-04 319.3  1.3e-03 
3 ∞  1.4e-02 300.6  8.0e-04 318.6  1.3e-03 
4 ∞  1.8e-02 285.4  9.1e-04 316.9  1.2e-03 
5 ∞  2.2e-02 275.8  1.1e-03 315.4  1.3e-03 
6 ∞  2.8e-02 267.0  1.2e-03 316.3  1.4e-03 
7 ∞  3.2e-02 257.1  1.4e-03 319.5  1.3e-03 
8 ∞  3.7e-02 243.9  1.5e-03 314.2  1.2e-03 
9 ∞  4.8e-02 231.7  1.6e-03 314.2  1.3e-03 
10 ∞  5.3e-02 219.7  1.8e-03 313.8  1.2e-03 
11 ∞  5.9e-02 209.2  1.9e-03 313.4  1.3e-03 
12 ∞  6.6e-02 200.4  2.0e-03 313.2  1.3e-03 
13 ∞  7.4e-02 189.9  2.3e-03 311.7  1.3e-03 
14 ∞  8.4e-02 182.4  2.5e-03 313.1  1.3e-03 
15 ∞  9.3e-02 169.6  2.6e-03 317.8  1.3e-03 
16 ∞  1.0e-01 160.3  2.7e-03 309.6  1.3e-03 
17 ∞  1.1e-01 147.7  2.9e-03 310.3  1.3e-03 
18 ∞  1.2e-01 135.4  3.0e-03 311.6  1.3e-03 
19 ∞  1.4e-01 115.7  3.1e-03 311.7  1.3e-03 
20 ∞  1.5e-01 105.2  3.3e-03 310.6  1.3e-03 
21 ∞  1.6e-01 96.4  4.4e-03 310.1  1.3e-03 
22 ∞  1.8e-01 84.6  4.7e-03 310.1  1.3e-03 
23 ∞  2.0e-01 81.5  4.8e-03 310.1  1.3e-03 
24 ∞  2.1e-01 65.2  5.2e-03 309.5  1.3e-03 
25 ∞  2.2e-01 59.4  5.3e-03 310.9  1.3e-03 
26 ∞  2.4e-01 49.3  5.6e-03 310.7  1.3e-03 
27 ∞  2.6e-01 41.6  5.9e-03 312.4  1.3e-03 
28 ∞  2.9e-01 29.3  6.3e-03 309.7  1.3e-03 
29 ∞  3.0e-01 20.9  6.5e-03 309.9  1.3e-03 
30 ∞  3.2e-01 8.8  6.8e-03 309.7  1.3e-03 
31 ∞  3.4e-01 0.6  6.9e-03 321.3  1.3e-03 

Discussion for the Global Regime

  • The numerical stability of the classic approach degrades sharply near the end of the range of degrees that we considered here, to the point of uselessness.

  • The numerical stability of the splinekit library is as good as the floating-point representation of numbers will allow, at all degrees shown here.