The behaviour of splines of different orders is explained in the table.
order | behaviour between knots | continuity at knots |
---|---|---|
0 | piecewise constant | not continuous at all |
1 | piecewise linear | lines are continuous (connected) |
2 | parabolic pieces | 1st derivatives are also continuous |
3 | cubic pieces | 2nd derivatives are also continuous |
n>3 | n-th order polynomials | (n-1)-th derivatives are also continuous |
From here on we work with cubic splines on a random grid unless otherwise indicated. It is quite easy to extend the ideas to lower or higher orders.
Assume we have a set of knots as in
|------------|------|---------|-------------------|--- k0 k1 k2 k3 k4 etc..The set of knots defines the domain where the spline model is valid.
We define a 3rd order polynomial to the first knot segement (k_{0},k_{1}).
f_{0}(x) = p_{0} + p_{1} x_{0} + p_{2} x_{0}^{2} + p_{3} x_{0}^{3} |
For the second segment (k_{1},k_{2}) we extend the function of the previous segment and add
f_{1}(x) = | p_{4} x_{1}^{3} | if x > k_{1} |
0 | if x < k_{1} |
Note that the function f_{1} is continuously differentiable to the 3rd derivative. For x < k_{1} it is all 0; for x > k_{1} it is a 3rd order polynomial and at x = k_{1} the derivatives from the up-side are all 0, equal to those from the low-side.
For the next segments we repeat this procedure.
So in the end we have a spline function f(x) as
f(x) = | p_{0} + p_{1} x_{0} + p_{2} x_{0}^{2} + p_{3} x_{0}^{3} + ∑ p_{i+4} x_{i}^{3} |
This simple algorithm has two disadvantages. Firstly it is not symmetric. We could also have started at the end of the domain and worked our way to the beginning. We would obtain a different spline representation, not better or worse than the first. Secondly the parameters are not localized. Each parameter influences the result over the complete domain. Preferably one parameters affects only one part of the model without influencing the others. In this algorithm a parameter affects all later parameters.
Even with these defects it is a usefull and fast algorithm. The algorithm is implemented in SplinesModel.py.
In figure 1 the model is displayed with all its components.
Figure 1 shows the the splines model and all its components. |
Figure 1 is made with the following piece of python code.
from BayesicFitting import SplinesModel import matplotlib.pyplot as plt x = numpy.linspace( 0, 10, 101, dtype=float ) kn = numpy.linspace( 0, 10, 6, dtype=float ) sm = SplinesModel( knots=kn ) par = [1.989, 0.066, 0.010, 0.016, -0.046, -0.028, 0.210, -0.286] sm.parameters = par part = sm.partial( x, par ) yf = numpy.zeros( 101, dtype=float ) for k in range( sm.npars ) : yf += part[:,k] * par[k] plt.plot( x, yf ) plt.plot( x, sm.result( x, par ), 'k-' ) plt.title( "Splines model and its constituents" ) plt.show()
Carl de Boor (in "A Practical Guide to Splines.", (1978) Springer-Verlag. ISBN 978-3-540-90356-7.) designed a recursive algorithm where each recursion produces a set of base splines of a higher order.
John Foster and Juha Jeronen implemented de Boor's algorithm in Python in the bspline package. We encapsulated the bspline package into our BSplinesModel.py.
In figure 2 a set of basis splines of order 3 is displayed, using parameters all equal to 1.0. The knots are at [0,1,4,5,6,7,8,10]; 8 knots resulting in 10 basis functions and as many parameters. The sum of the components is the black line constant at 1.0.
Figure 2 shows the de Boor splines model and all its components. The position of the knots is indicated with red bars low in the plot. |
Unfortunately the recursive implementation in BsplinesModel is about 20 times slower than the the splines in SplinesModel.
We have a 3rd order polynomial between the first 2 knots.
f(x) = a + b x + c x^{2} + d x^{3}
f_{1} |------------|---- x_{0} x_{1}The property 3-smooth implies that the function and it first and second derivative are 0 at x_{1}, while the function is 1 at x_{0}.
f_{1}(x_{1}) = 0
f_{1}'(x_{1}) = 0
f_{1}"(x_{1}) = 0
f_{1}(x_{0}) = 1
Now we have 4 equations with 4 unknown coefficients (a,b,c,d) which we can solve.
f_{1} f_{2} |------------|------|--- x_{0} x_{1} x_{2}
f_{1}(x) = a_{1} + b_{1} x + c_{1} x^{2} + d_{1} x^{3}
f_{2}(x) = a_{2} + b_{2} x + c_{2} x^{2} + d_{2} x^{3}
The function, f_{1} is 1-smooth at x_{0}. f_{1} and f_{2} are 3-smooth at x_{1} and f_{2} is 3-smooth with 0 at x_{2}. Again we need a normalization at x_{1}.
f_{1}(x_{0}) = 0
f_{1}(x_{1}) - f_{2}(x_{1}) = 0
f_{1}'(x_{1}) - f_{2}'(x_{1}) = 0
f_{1}"(x_{1}) - f_{2}"(x_{1}) = 0
f_{2}(x_{1}) = 0
f_{2}'(x_{1}) = 0
f_{2}"(x_{1}) = 0
f_{1}(x_{1}) = 1
We have 8 equations and 8 coefficients (4 per function) so it can be solved directly.
blob | x_{0} | x_{1} |
---|---|---|
1 | 0-smooth | 3-smooth |
2 | 1-smooth | 2-smooth |
3 | 2-smooth | 1-smooth |
4 | 3-smooth | 0-smooth |
In case we have a cubic spline on a set of 3 knots, x_{0} .. x_{12}, the smoothnes at the knots is given in the table below.
blob | x_{0} | x_{1} | x_{2} |
---|---|---|---|
1 | 0-smooth | 3-smooth | n/a |
2 | 1-smooth | 3-smooth | 3-smooth |
3 | 2-smooth | 3-smooth | 2-smooth |
4 | 3-smooth | 3-smooth | 1-smooth |
5 | n/a | 3-smooth | 0-smooth |
S(x:p) = ∑ p_{k} f_{k}( x )
S(x:p) = ∑ p_{k} ( a_{k} + b_{k} x + c_{k} x^{2} + d_{k} x^{3} )
S(x:p) = ∑ ( p_{k} a_{k} ) + ∑ ( p_{k} b_{k} ) x + ∑ ( p_{k} c_{k} ) x^{2} + ∑ ( p_{k} d_{k} ) x^{3}
The coefficients and the parameters can be multiplied together, forming
one piecewise 3-order polynomial function of the x's.
This algorithm proved to be 10 times faster than the recursive one.
Still slower by a factor 2 than the first one.
The advantages we obtained though, are that the blobs are independent.
when one parameter changed only localized effects are seen. Even when
the position of a knot shifts, adjacent knots are only slightly
affected. This opens the option of a DynamicBasicSplinesModel.