Evaluation of Splines

How to evaluate a periodic one-dimensional spline at just one argument or at a series of arguments.

Linear-Algebra Formulation

The objects that the class PeriodicSpline1D operates upon are mathematical functions and consist in descriptions of mappings \(f:{\mathbb{R}}\rightarrow{\mathbb{R}},x\mapsto f(x).\) Oftentimes, the class allows one to obtain a new mapping from an existing one, for instance the gradient of a spline (not just a number but the whole function \(\dot{f}\)) out of its original form \(f.\) This is all good, but one is sometimes also interested in the more mundane goal of experimenting with a fixed mapping, typically to ask to what number \(f(x)\in{\mathbb{R}}\) is the argument \(x\in{\mathbb{R}}\) mapped to by \(f:{\mathbb{R}}\rightarrow{\mathbb{R}}.\) The splinekit library offers two possibilities to address this goal.

  • PeriodicSpline1D.at returns the value \(f(x)\) of the spline \(f\) evaluated at \(x.\)

  • PeriodicSpline1D.get_samples returns an array of values at arguments separated by a constant step.

The two possibilities rely on the recipe followed by all one-dimensional uniform polynomial splines of nonnegative integer degree \(n\in{\mathbb{N}},\) according to which an argument \(x\) is mapped to the value

\[f(x)=\sum_{k\in{\mathbb{Z}}}\,c[{k\bmod K}]\,\beta^{n}(x-\delta x-k),\]

where \(K\in{\mathbb{N}}+1\) is a positive integer period, \(c\) is an arbitrary array of \(K\) spline coefficients, \(\beta^{n}\) is a B-spline whose degree \(n\) is typeset in superscript (not a power), and \(\delta x\in{\mathbb{R}}\) is an arbitrary delay. Put simply, a spline is a weighted sum of shifted B-splines.

Since the support of a B-spline is finite, the sum is finite at any given argument \(x.\) Moreover, since B-splines of positive degrees are themselves made of unit-length pieces of polynomials, one can deploy the formalism of linear-algebra to make the expression of the spline recipe take the form

\[\begin{split}\begin{array}{rcl} \forall n\in{\mathbb{N}}+1,\forall x\in{\mathbb{R}}:f(x)&=&\left(\begin{array}{c}c[{\left(-r\right)\bmod K}]\\c[{\left(1-r\right)\bmod K}]\\c[{\left(2-r]\right)\bmod K}\\\vdots\\c[{\left(n-r\right)\bmod K}]\end{array}\right)^{{\mathsf{T}}}\,\left(\begin{array}{ccccc}w_{0,0}^{n}&w_{0,1}^{n}&w_{0,2}^{n}&\cdots&w_{0,n}^{n}\\w_{1,0}^{n}&w_{1,1}^{n}&w_{1,2}^{n}&\cdots&w_{1,n}^{n}\\w_{2,0}^{n}&w_{2,1}^{n}&w_{2,2}^{n}&\cdots&w_{2,n}^{n}\\\vdots&\vdots&\vdots&\ddots&\vdots\\w_{n,0}^{n}&w_{n,1}^{n}&w_{n,2}^{n}&\cdots&w_{n,n}^{n}\end{array}\right)\,\left(\begin{array}{c}1\\v\\v^{2}\\\vdots\\v^{n}\end{array}\right)\\ &=&{\color{blue}{{\mathbf{c}}}}^{{\mathsf{T}}}\,{\color{blue}{{\mathbf{W}}^{n}}}\,{\color{blue}{{\mathbf{v}}^{n}}}, \end{array}\end{split}\]

where \({\color{blue}{{\mathbf{c}}}}\in{\mathbb{R}}^{n+1}\) is a vector whose \(n+1\) components are extracted from the data-dependent array \(c\) at some integer \(r\in{\mathbb{Z}}\) that depends in discrete fashion on \(\left(x-\delta x\right),\) where \({\color{blue}{{\mathbf{W}}^{n}}}\in{\mathbb{R}}^{\left(n+1\right)\times\left(n+1\right)}\) is a matrix that depends on \(n\) only, and where \({\color{blue}{{\mathbf{v}}^{n}}}\in{\mathbb{R}}^{n+1}\) is a Vandermonde vector whose every component belongs to the interval \([0,1]\) and is made to depend continuously on \(\left(x-\delta x\right).\)

Evaluation at a Single Argument

To evaluate a one-dimensional periodic spline \(f\) of positive degree \(n\) at some argument \(x,\) it is enough to compute \({\mathbf{c}}^{{\mathsf{T}}}\,{\mathbf{W}}^{n}\,{\mathbf{v}}^{n}.\) This is precisely what is done by the function call f.at(x) in the piece of code below, which otherwise generates and plots a random spline of some arbitrary degree and delay.

Jupyter Lab notebook

Single argument

Evaluation at Multiple Arguments

We want now to evaluate \(f(x)\) at the \(L\,M\) arguments \(x\in\left\{x_{0}+k/M\right\}_{k=0}^{L\,M-1},\) where \(M\in{\mathbb{N}}+1\) is some positive integer oversampling factor, \(L\in{\mathbb{N}}+1\) is the length of the support over which we want to sample \(f,\) and \(x_{0}\in{\mathbb{R}}\) is the argument of the first sample. To do so, we take advantage of PeriodicSpline1D.get_samples. In the next piece of code, we illustrate how the combination of \(L,\) \(M,\) and \(x_{0}\) can be used to specify the arguments where the spline is evaluated.

Jupyter Lab notebook

Multiple arguments

Experimental Performance

When we evaluate a spline jointly at multiple arguments spaced regularly as above, it is judicious to pay attention on how the terms of \({\mathbf{c}}^{{\mathsf{T}}}\,{\mathbf{W}}^{n}\,{\mathbf{v}}^{n}\) depend on the variables of interest. In particular, when we compare the computation of \(f(x_{0}+q/M)\) to that of \(f(x_{0}+\left(q+M\right)/M),\) a detailed analysis reveals that the vector \({\mathbf{c}}\) changes, while the vector resulting from the product \({\mathbf{W}}^{n}\,{\mathbf{v}}^{n}\) is identical in the two cases. This invariance leads to computational savings.

To ascertain their benefit, we propose to time the determination of regularly spaced spline samples through either repeated calls to PeriodicSpline1D.at or through a combined call to PeriodicSpline1D.get_samples. We report in degree-dependent tables by how many times the global call is faster than the repeated independent calls. For simplicity, we sample a whole period and set \(L=K,\) letting the period \(K\) explore some wide range. We let the oversampling factor \(M\) take values in \(1\) (unit-spaced samples) through \(6\) (samples spaced by \(1/6\)). Finally, to improve the robustness of the reported results, we repeat each timing experiment \(10\) times.

Danger

  • The link below allows you to inspect the notebook. Unfortunately, running it from the browser is meaningless: the timings of the installation-free version are not representative 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 Section correspond to those of the native execution.

Jupyter Lab notebook

Scalar vs array

Hint

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

Results

On a desktop computer of 2021, typical resulting tables are as follows.

Acceleration for Spline Degree n = 1
Period  2  5  10  20  50  100  200  500  1000 
Oversampling 
M = 1  1.00  2.13  3.15  4.41  3.96  5.20  5.36  5.64  4.72 
M = 2  0.98  2.32  3.64  5.13  7.03  7.95  8.39  8.90  8.92 
M = 3  1.43  2.93  4.68  6.73  9.83  11.17  12.68  13.36  13.62 
M = 4  1.50  3.25  5.47  8.24  12.31  14.85  16.44  17.78  17.74 
M = 5  1.74  3.54  6.01  9.39  14.61  17.77  20.19  21.71  22.38 
M = 6  1.72  3.79  6.48  10.35  16.39  20.54  23.96  25.91  26.70 

Acceleration for Spline Degree n = 2
Period  2  5  10  20  50  100  200  500  1000 
Oversampling 
M = 1  1.16  2.09  3.01  3.91  4.78  5.28  5.56  5.66  5.68 
M = 2  1.26  2.35  3.59  5.13  6.96  7.97  8.56  8.97  8.93 
M = 3  1.50  2.94  4.64  6.81  9.83  11.46  12.86  13.34  13.57 
M = 4  1.67  3.24  5.31  8.09  12.45  14.50  16.10  17.47  18.04 
M = 5  1.74  3.57  5.87  9.19  14.50  17.83  20.38  21.80  22.51 
M = 6  1.92  3.92  6.36  10.26  16.41  20.38  23.92  25.73  26.66 

Acceleration for Spline Degree n = 3
Period  2  5  10  20  50  100  200  500  1000 
Oversampling 
M = 1  1.18  1.96  2.86  3.81  4.79  5.25  5.54  5.61  5.67 
M = 2  1.24  2.30  3.49  4.93  6.84  7.89  8.50  8.87  9.11 
M = 3  1.55  2.99  4.58  6.67  9.76  11.47  13.11  13.23  13.66 
M = 4  1.67  3.32  5.30  8.00  12.18  14.72  16.23  17.75  17.86 
M = 5  1.82  3.61  5.92  9.12  14.26  17.32  20.27  21.66  22.68 
M = 6  1.91  3.89  6.39  10.03  16.05  20.47  23.82  25.92  26.19 

Acceleration for Spline Degree n = 4
Period  2  5  10  20  50  100  200  500  1000 
Oversampling 
M = 1  1.19  1.93  2.77  3.58  4.63  5.14  5.49  5.61  5.68 
M = 2  1.26  2.27  3.50  4.90  6.78  7.80  8.60  8.87  9.20 
M = 3  1.52  2.96  4.54  6.53  9.67  11.36  12.44  13.28  13.43 
M = 4  1.74  3.33  5.27  7.86  11.82  14.51  16.59  17.88  17.94 
M = 5  1.77  3.78  5.92  9.04  14.23  17.09  19.86  21.93  22.42 
M = 6  1.85  4.00  6.36  9.87  15.98  20.05  23.92  25.98  26.25 

Acceleration for Spline Degree n = 5
Period  2  5  10  20  50  100  200  500  1000 
Oversampling 
M = 1  1.21  1.86  2.61  3.50  4.55  5.12  5.53  5.56  5.69 
M = 2  1.23  2.40  3.33  4.65  6.58  7.74  8.51  8.85  9.05 
M = 3  1.49  3.02  4.32  6.41  9.57  10.90  12.31  13.21  13.43 
M = 4  1.71  3.51  5.14  7.66  11.86  14.16  16.05  17.43  17.88 
M = 5  1.82  3.79  5.79  8.83  13.87  17.24  19.97  21.59  22.20 
M = 6  1.93  4.19  6.39  9.85  15.83  19.82  23.03  25.60  26.53 

Acceleration for Spline Degree n = 6
Period  2  5  10  20  50  100  200  500  1000 
Oversampling 
M = 1  1.19  1.84  2.47  3.35  4.41  4.98  5.46  5.58  5.62 
M = 2  1.26  2.36  3.27  4.59  6.49  7.61  8.38  8.91  8.96 
M = 3  1.54  3.05  4.27  6.21  9.22  10.12  12.25  13.07  13.37 
M = 4  1.74  3.56  5.21  7.55  11.82  13.53  16.15  17.50  17.77 
M = 5  1.84  3.93  5.81  8.70  13.51  16.72  19.62  21.53  22.30 
M = 6  1.90  4.22  6.35  9.71  15.67  19.47  23.18  25.44  26.61 

Acceleration for Spline Degree n = 7
Period  2  5  10  20  50  100  200  500  1000 
Oversampling 
M = 1  1.15  1.82  2.32  3.18  4.29  4.92  5.33  5.60  5.58 
M = 2  1.23  2.37  3.17  4.42  6.32  7.39  8.20  8.71  8.97 
M = 3  1.59  3.05  4.27  6.03  8.79  11.06  12.15  13.11  13.36 
M = 4  1.75  3.52  5.04  7.37  9.16  13.77  15.80  17.45  18.00 
M = 5  1.83  3.89  5.71  8.53  13.43  16.77  19.49  21.47  22.55 
M = 6  1.98  4.22  6.30  9.55  15.05  19.28  23.05  25.61  26.28 

Acceleration for Spline Degree n = 8
Period  2  5  10  20  50  100  200  500  1000 
Oversampling 
M = 1  1.16  1.79  2.23  3.03  4.12  4.82  5.30  5.53  5.58 
M = 2  1.29  2.39  3.13  4.34  6.22  7.26  8.22  8.74  8.91 
M = 3  1.59  3.05  4.13  5.88  9.05  10.61  12.07  12.87  13.47 
M = 4  1.69  3.56  4.98  7.24  9.13  13.60  15.73  17.38  17.58 
M = 5  1.92  3.97  5.68  8.35  13.16  16.59  19.47  21.11  22.08 
M = 6  2.01  4.27  6.32  9.42  15.10  19.25  22.72  25.02  25.68 

Acceleration for Spline Degree n = 9
Period  2  5  10  20  50  100  200  500  1000 
Oversampling 
M = 1  1.15  1.75  2.16  2.92  4.02  4.70  5.23  5.50  5.60 
M = 2  1.30  2.37  3.06  4.25  6.03  7.17  8.12  8.78  8.92 
M = 3  1.59  3.06  4.13  5.76  8.96  10.33  12.16  12.98  13.42 
M = 4  1.79  3.56  4.95  7.10  10.77  13.69  14.68  17.40  17.55 
M = 5  1.85  3.91  5.66  8.20  12.90  16.29  19.33  21.35  22.24 
M = 6  1.96  4.28  6.33  9.45  14.22  19.29  22.71  25.16  25.79 

Discussion

The gain in computational efficiency is substantial. It increases with the oversampling factor and with the period of the spline. The degree of the spline does not play much of a role.