- Email: [email protected]
Locally optimal knots and tension parameters for exponential splines
Journal of Computational and Applied Mathematics 196 (2006) 94 – 114 www.elsevier.com/locate/cam
Locally optimal knots and tension parameters for exponential splines Karl O. Riedela, b,∗ a Mathematical Institute, University of Cologne, Weyertal 86-90, 50931 Cologne, Germany bWeißensteinerstraße 74, 73525 Schwäbisch Gmünd, Germany
Received 24 August 2004; received in revised form 11 August 2005
Abstract Tension can be applied to cubic splines in order to avoid undesired spurious oscillations. This leads to the well-known (exponential) spline in tension. It is crucial but unfortunately difficult to find suitable tension parameters of interpolating splines in tension. Instead of heuristics, we propose a simultaneous knot placing and tension setting algorithm for least-squares splines in tension which includes interpolating splines in tension as a special case. Moreover, the splines presented here are the foundation of exponential surface splines on fairly arbitrary meshes [K.O. Riedel, Two-dimensional splines on fairly arbitrary meshes, ZAMM—Z. Angew. Math. Mech. 85(3) (2005) 176–188]. © 2005 Elsevier B.V. All rights reserved. MSC: 65D07; 65D10; 65D05 Keywords: Spline; Tension; Exponential spline; Interpolation; Approximation
1. Introduction Cubic splines, i.e. twice continuously differentiable functions consisting of piecewise cubic polynomials, are widespread. They are used to interpolate or approximate data. Interpolating cubic splines definitely show less oscillations than interpolating polynomials. Nevertheless, in some cases spurious oscillations of the splines occur which are not indicated by the data. These oscillations can be reduced by shortening the length of the spline. To this purpose, an additional tension parameter is introduced. This leads to the (exponential) spline in tension, i.e., a twice continuously differentiable function consisting piecewise of a linear combination of the functions 1, x, exp(i x), and exp(−i x), where i > 0 is the interval tension parameter. Some authors call this kind of spline in tension just exponential spline, whereas others understand by an exponential spline a generalization or a different spline with exponential terms. Bsplines for splines in tension have been introduced by [16], see also [17,18]. Explicit formulas also can be found in [26]. We will use them in order to calculate smoothing (least-squares) splines in tension. Often convexity constraints are given with or imposed to interpolating or approximating splines. Especially, the tension parameters of tension splines can be suitably adapted in order to fulfill these shape constraints. Yet, in this ∗ Corresponding address: Weißensteinerstraße 74, 73525 Schwäbisch Gmünd, Germany. Tel.:+49(0)1706837882; fax: +49(0)6990553820.
E-mail address: [email protected]. 0377-0427/$ - see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2005.08.012
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114
95
paper, we do not deal with convexity constraints, but simultaneously choose knot positions and tension parameters which minimize the least-squared error. This yields visually pleasing curves. Smooth curves as well as sharp corners are captured by a minimal number of knots, which implies a highly reduced tendency to spurious oscillations, i.e., a strong tendency to fulfill convexity constraints. This is an alternative at least to imposing convexity constraints. Finding suitable knot positions of a least-squares spline is important but unfortunately difficult. Dierckx [7] has proposed a successive knot placing algorithm which minimizes the least-squared error as a function of the knot positions. We generalize this approach considering the dependence on the interval tension parameters of exponential splines simultaneously. Further, we introduce a local minimal distance condition in order to avoid overfitting. This, automatically, leads to the interpolating spline if the number of knots and data points is identical. Thus, the same algorithm can be used in order to find interpolating as well as least-squares splines. Considering interpolating splines, this includes the generation of suitable tension parameters which have been determined by heuristic methods in former approaches. This paper is organized as follows: In Section 2, we discuss the problem of spurious oscillations. Splines in tension, and for data without sharp corners also cubic splines with suitable knots are a remedy. The simultaneous knot placing and tension setting algorithm is described in Section 3. Finally, we demonstrate its efficiency in finding suitable leastsquares splines for various data in Section 4. Interpolating splines are automatically included. 1.1. Related work The spline in tension was first analytically modeled by Schweikert [31]. He further defines an extraneous inflection point as a zero crossing of the second derivative which is not motivated by the data and shows that the exponential spline is free of extraneous inflection points for a tension large enough. Späth [32] simplifies this proof. Thus, finding tension parameters for a spline without spurious oscillations in principle is possible. One choice is infinite tension leading to a polygonal line. However, the task of finding tension parameters for a smooth, visually pleasant spline without extraneous inflection points is neither uniquely solvable nor trivial. Späth [32] and Rentrop [25] propose heuristic methods in order to determine tension parameters for which the exponential spline is free of extraneous inflection points. Heß [13] and Renka [23] also describe methods for finding suitable tension parameters of interpolating splines. Heidemann [12] uses the tension parameters proposed in [25] for least-squares splines. Lynch [19] presents an iterative scheme for constructing shape preserving exponential splines under the stringent constraint of uniform tension. Sapidis et al. [28] propose a generalized Newton Raphson method for finding suitable tension parameters. Pruess [22] shows that the approximation of sufficiently smooth functions with exponential splines is of fourth order, as for cubic splines. For large tension parameters the exponential spline is, essentially, locally linear. He exploits this fact in order to construct convex or monotone spline approximants. Nielson [20,21] proposed a polynomial alternative to the exponential spline, which he called the -spline. For a rational spline with interval and point tension see [10,29,30]. For a derivation of cubic and exponential splines from a variational principle and a comparison see [1]. Related work for monotony preserving interpolation is found, e.g., in [8,9,5,6,14]. For a software package on exponential splines see [24]. Jüttler [15] proposes the use of a reference curve in order to impose convexity constraints, i.e., given inflection points, to the approximating spline. Costantini et al. [3,4] use variable degree polynomial splines instead of exponential splines in order to get shape-preserving approximants for 2D and 3D parametric curves.
2. Suitable knots and tension parameters: an example We consider an example of Späth [32]. Only the interpolating cubic spline shows spurious oscillations which disappear for the exponential spline with uniform tension 0 = 7 as shown in Fig. 1. Using suitable knots, there are also cubic splines without spurious oscillations—see Fig. 2. We found these knots by trial and error by hand. Unfortunately, it is not easy to find optimal knot positions automatically. Moreover, in order to form sharp corners with cubic splines a good placement of the knots is difficult and could be impossible. For an alternative, we refer to Section 3.3, where we use least-squares exponential splines with as few as possible knots and optimize the interval tension parameters and the knot positions simultaneously.
96
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114 10 9 8 7 6 5 4 3
0
1
2
3
4
5
6
7
8
9
10
0
1
2
3
4
5
6
7
8
9
10
0
1
2
3
4
5
6
7
8
9
10
1 0.5 0
1 0.5 0
Fig. 1. (Top) Späth data (stars) are interpolated by a cubic spline (dashed) and an exponential spline (full line) with uniform tension 0 = 7. Natural boundary conditions are used for both splines. The corresponding exponential B-splines (bottom) are intermediate between cubic B-splines (middle) and hat functions, the basis functions of piecewise linear splines.
10
10
8
8
6
6
4
4 0
2
4
6
8
10
0
2
4
6
8
10
Fig. 2. Cubic splines with natural boundary conditions interpolating Späth data (stars). The knots of the splines are depicted by squares. Changing the position of just one knot marked by arrows leads to strongly increased oscillations.
3. Data approximation by splines in tension 3.1. Least-squares splines Dealing with noisy data it does not make sense to use an interpolating spline. Like for interpolation, an approximating spline is preferable to a polynomial, because it shows less oscillations. In this section, we state how to find the smoothing spline which minimizes the least-squared error for given knots and tension parameters. For finding suitable knots and tension parameters see Section 3.3. Given the data points (tr , xr ),
r = 1, . . . , m,
we use at most the same number of B-splines with natural boundary conditions Bi (t),
i = 0, . . . , n + 1, n + 2 m,
(3.1)
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114
97
with the knots 0 , . . . , n+1 which are placed such that there exist n data points tri ∈ {t2 , . . . , tm−1 }, i = 1, . . . , n, which fulfill 0 = t1 ,
i−1 tri i , i = 1, . . . , n,
n+1 = tm .
(3.2)
In order to find the B-spline representation of the approximating spline, we search the coefficients ci , i = 0, . . . , n + 1, which minimize 2 m n+1 2 wr xr − ci Bi (tr ) = min (3.3) r=1
i=0
for given weighting factors wr > 0, r =1, . . . , m. For the numerical examples presented here we use wr =1, r =1, . . . , m. Denoting ⎞ ⎞ ⎛ ⎛ w1 x1 c0 . . c = ⎝ .. ⎠ , xw = ⎝ .. ⎠ , E = (wr Bi (tr ))(r,i) , cn+1 w m xm for r = 1, . . . , m, i = 0, . . . , n + 1, we search the coefficient vector c of the least-squares spline which minimizes min Ec − xw 22 .
(3.4)
c
Lemma 1. Given data points and knots fulfilling Eq. (3.2), the least-squares spline which minimizes (3.3) is unique. Proof. The least-squares spline can be defined by its coefficient vector c. In order to show that c is unique we have to show that E has full rank (see, e.g., [34]). Thus, it is sufficient to show Ec = 0 implies ci = 0, i = 0, . . . , n + 1. Having r > 0, r = 1, . . . , m, we can set without loss of generality r = 1, r = 1, . . . , m. Then Ec = 0 can be written as n+1
ci Bi (tr ) = 0,
r = 1, . . . , m.
i=0
Denoting s(t) :=
n+1
ci Bi (t)
i=0
this reads s(tr ) = 0,
r = 1, . . . , m.
Using the notations tr0 := t1 and trn+1 := tm this implies s(tri ) = 0,
i = 0, . . . , n + 1 where i−1 tri i , i = 1, . . . , n.
The spline s(t) is twice continuously differentiable and the theorem of Rolle yields n + 1 points i with s (i ) = 0,
tri−1 < i < tri , i = 1, . . . , n + 1.
Applying Rolle again we have n points i fulfilling s (i ) = 0,
i < i < i+1 , i = 1, . . . , n.
(3.5)
98
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114
Setting 0 := tr0 and n+1 := trn+1 and considering natural boundary conditions we have s (i ) = 0,
i = 0, . . . , n + 1,
i.e., n + 2 zeros i within n + 1 intervals [i , i+1 ], i = 0, . . . , n. Hence, there exists an interval J := [j , j +1 ] which contains two zeros k and k+1 . Considering Eq. (3.5) we have k < k+1 < k+1 , i.e., k+1 ∈ J . Further, Eq. (3.2) guarantees the existence of a data point til ∈ J . Summarizing, we have s(til ) = 0,
til ∈ J ,
s (k+1 ) = 0, s (k ) = 0,
(3.6)
k+1 ∈ J ,
(3.7)
k ∈ J ,
s (k+1 ) = 0,
(3.8)
k+1 ∈ J .
(3.9)
Within interval J for a given tension parameter > 0 the exponential spline s(t) can be written as s(t) = a + bt + cet + de−t , s (t) = b + cet − de−t , s (t) = c2 et + d2 e−t ,
t ∈ J, t ∈ J,
t ∈ J.
Hence, Eq. (3.8) reads −k c 2
e k +d 2 e
=0, >0
(3.10)
>0
which yields cd 0.
(3.11)
Subtraction of Eq. (3.8) from (3.9) gives c 2 (ek+1 − ek ) +d 2 (e−k+1 − e−k ) =0,
>0
<0
which implies cd 0.
(3.12)
Considering Eqs. (3.11) and (3.12), we have cd = 0. Hence, Eq. (3.10) implies c = d = 0. Hence, (3.7) reads b=0 and (3.6) reads a = 0. Thus, we have s(t) ≡ 0,
t ∈ [j , j +1 ].
With the tension parameter ˜ > 0 on the next interval [j +1 , j +2 ] we write s(t) as ˜ ˜ − j +1 ) + ce ˜ −˜ (t−j +1 ) . s(t) = a˜ + b(t ˜ (t−j +1 ) + de
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114
99
Considering that s is twice continuously differentiable yields 2 ˜ = 0, s (j +1 ) = ˜ (c˜ + d)
i.e., ˜ c˜ = −d. Hence, s(j +1 ) = a˜ = 0 and s (j +1 ) = b˜ + 2c˜˜ = 0, i.e., ˜ b˜ = −2c˜. Eq. (3.2) guarantees another zero trj +2 ∈ (j +1 , j +2 ], i.e., ˜ r − j +1 ) + ce ˜ s(trj +2 ) = −2c˜(t j +2
˜ (trj +2 −j +1 )
− ce ˜
−˜ (trj +2 −j +1 )
= 0,
which reads ˜ r − j +1 )] =0. ˜ r − j +1 )) − (t c˜ [sinh((t j +2 j +2
>0
This implies c˜ = 0
and b˜ = 0.
Thus, we have s(t) ≡ 0,
t ∈ [j +1 , j +2 ].
Iteration also to the left-hand side yields s(t) =
n+1
ci Bi (t) ≡ 0,
t ∈ [t1 , tm ].
(3.13)
i=0
As the B-splines Bi (t), i = 0, . . . , n + 1 are linearly independent on [t1 , tm ], Eq. (3.13) implies ci = 0, i = 0, . . . , n + 1. For calculating the coefficients c of the least-squares spline by minimizing (3.4) see, e.g., [33]. 3.2. Pseudo-arclength parametrization So far, in the functional setting the abscissae and ordinates of the data points were not treated equally. Only functions x(t) could be represented. For many applications this is appropriate. Yet, in order to be able to include also parametric curves like, e.g., the exponential spline representation of contours of camera images [27], we use the parametric setting (x(t), y(t)). Moreover, generalizing rectangular meshes, exponential surface splines are introduced in [27] on fairly arbitrary meshes, which are represented by the more general parametric exponential splines presented here. Using parametric splines, Eq. (3.4) is to be solved twice, whereas the computational effort for calculating the B-splines is not increased. If computational time does not matter, also data with a functional dependency can be represented by parametric splines without changing the code, provided that the function x = x(t) = sx (t) is invertible. Usually, the
100
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114
monotonicity of the data (ti , xi ) leads to a monotonous, and hence invertible, spline sx (t). In such a situation, parametric splines look much the same as their functional counterparts. This is why we do not explicitly investigate functional splines, although only the data presented in Figs. 10 and 11 demand for the parametric setting. If the functional setting is required explicitly the code can be adapted by minor changes, i.e., using the given abscissae instead of pseudo-arclength parametrization and leaving out the calculation of the second least-squares spline sy (t). Parametric spline curves suggest the use of pseudo-arclength parametrization. Thus, starting from the data points r = 1, . . . , m,
(xr , yr ),
(3.14)
we introduce t1 := 0, r tr := (xk − xk−1 )2 + (yk − yk−1 )2 ,
r = 2, . . . , m.
(3.15)
k=2
For given tension parameters i > 0, i = 1, . . . , n + 1, and joint knots i , i = 0, . . . , n + 1, in the interval [t1 , tm ] we calculate joint B-splines in tension [16,18,26] Bi (t),
t ∈ [t1 , tm ], i = 0, . . . , n + 1.
Now, the problem splits into two parts: for each partial set of data points (tr , xr ), (tr , yr ),
r = 1, . . . , m, r = 1, . . . , m,
we find the coefficients of the least-squares spline according to Section 3.1. Denote them by cx,i and cy,i , respectively. Hence, the least-squares spline with pseudo-arclength parametrization is given by n+1 n+1 sx (t) := cx,i Bi (t), sy (t) := cy,i Bi (t) , t ∈ [t1 , tm ]. (3.16) i=0
i=0
This treatment includes interpolating splines, if there are as many knots as data points. 3.3. Finding optimal knot positions and tension parameters As demonstrated in Section 2 finding suitable knots and tension parameters is difficult. We use the knot placing algorithm of Dierckx [7] for successively placing knots. However, we do not only adapt the knots of (cubic) splines, but as well the tension parameters of exponential splines simultaneously. We further introduce a local minimal distance condition, which insures that the distance between the knots is not smaller than the distance between the nearest data points. Thus, we avoid oscillations caused by overfitting. We always use the two outermost data points t1 and tm as knots 0 and n+1 . For given interior knots Tn := {1 , . . . , n } and interval tension parameters Ln := {1 , . . . , n+1 } we consider equadr (Tn , Ln ) :=
m
wr2 {(xr − sx (tr ))2 + (yr − sy (tr ))2 }
(3.17)
r=1
as a measure for the discrepancy between the data points (xr , yr ), r = 1, . . . , m, and the least-squares spline (sn,x (t), sn,y (t)), t ∈ [t1 , tm ], with n interior knots. In order to get tr , r = 1, . . . , m, we use pseudo-arclength parametrization according to Section 3.2.
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114
101
3.3.1. Initial values for the knots and tension parameters The success of a free knot least-squares algorithm for minimizing equadr (Tn , Ln ) will strongly depend on the choice of the initial values, i.e., the initial interior knots T0n and the initial tension parameters L0n . Indeed, due to the nonconvexity of equadr (Tn , Ln ) we can only expect to find a local minimum in the vicinity of T0n , L0n . Therefore, in order to find the global minimum, the determination of good initial knots and interval tension parameters is an important but unfortunately also a non-trivial task. A straightforward equidistant distribution of knots, for example, will rarely yield the global minimum of equadr (Tn , Ln ). Furthermore, we do not know in advance the required number of knots for an acceptable least-squares fit. Therefore, we use a scheme which will hopefully return the successive least-squares splines (sn,x (t), sn,y (t)), n = 1, 2, . . . , and which uses the locally optimal knot/tension set T∗n , L∗n in order to find suitable initial values T0n+1 , L0n+1 for determining (sn+1,x (t), sn+1,y (t)): Algorithm 1 (Successive knot placing strategy). • Initialization: Choose a scaled (see Section 3.4) tension parameter, e.g., med = 6, leading to a spline far apart from a cubic spline as well as from a polygonal line: 01 =
t1 + t m , 2
01 =
med , 01 − t1
02 =
med . tm − 01
• For n = 1, 2, . . . until the approximating spline is good enough: ◦ Choose the allowed interval [min , max ] for the scaled tension parameters. ◦ Starting with {Tn , Ln } determine the optimal knot positions and tension parameters T∗n , L∗n and the corresponding least-squares spline (sn,x (t), sn,y (t)). For this, we choose a minimization algorithm which limits the interval tension parameters to i ∈ [min /(i − i−1 ), max /(i − i−1 )], i = 1, . . . , n + 1, and forces the distance between the knots to be at least as large as the distance of the neighboring data points. ◦ Compute the relative error measure numbers i , i = 1, . . . , n + 1:
qi +mi 2 2 2 r=qi +1 wr ([xr − sn,x (tr )] + [yr − sn,y (tr )] ) i = , m − mi tqi < ∗i−1 tqi +1 < tqi +2 < · · · < tqi +mi ∗i < tqi +mi +1 , ∗0 := t0 ,
∗n+1 := tm .
◦ Find suitable initial knot positions and tension parameters in order to determine (sn+1,x (t), sn+1,y (t)): Let l denote the interval with the largest relative error measure l = max{1 , . . . , n+1 } and define for this interval the new tension parameter ∗new := max{∗l , 2min /(∗l − ∗l−1 )}. This leads to 0i
=
0i
=
i = 1, . . . , l − 1, ∗i , (∗l−1 + ∗l )/2, i = l, ∗i−1 , i = l + 1, . . . , n + 1, i = 1, . . . , l − 1, ∗i , ∗new , i = l, l + 1, ∗i−1 , i = l + 2, . . . , n + 2.
◦ Stop when the approximation is good enough. So, this algorithm successively locates the initial position of an additional knot in the middle of the knot interval [l−1 , l ], i.e., in a part of [t1 , tm ] where on the one hand (sn,x (t), sn,y (t)) has a large weighted sum of squared residuals, but where on the other hand the concentration of knots is not too high (m − ml not too large). Within the frame of Algorithm 1, we proceed with describing scaled tension parameters (Section 3.4) and a minimization algorithm for finding optimal knots which are not closer than the data points (Section 3.5). In this way, we avoid overfitting. Finally, we address the task of fixing the number of the knots (Section 3.6).
102
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114 1 λscal=20 λscal=6 λscal=2
0.8 0.6
cubic
0.4 0.2 0
0
1
2
3
4
5
Fig. 3. Exponential B-splines with different tension ranging from cubic splines to hat functions. Scaled tension parameters ranging from 2 to 20 include a quite large portion of possible exponential splines.
tli o
o
*o τi-1
o
tli+1 x o
o
τi-1 + τi 2
*o
o
τi
Fig. 4. Fixing locally the minimal knot distance.
3.4. Scaled tension parameters The shape of an exponential spline depends on the interval tension parameters as well as on the distance between the data points. This means that interpolating or least-squares splines are not size invariant. Linear compression of the data points leads to larger squared second derivatives and hence to a spline with a more pronounced cubic behavior. This can be overcome by using scaled interval tension parameters i,scal = i (i − i−1 ),
i = 1, . . . , n + 1,
(3.18)
for the knots 0 , . . . , n+1 . Scaled tension parameters yield size invariant splines. In other words, stretching the data points by a factor 2 and dividing the tension parameters by 2 leads to the exponential spline stretched by a factor 2. For a proof see [26]. Fig. 3 shows exponential B-splines with different scaled tension parameters. Consider that scaled tension parameters of 2 and 20 lead to a spline quite close to a cubic spline and a hat function, respectively. A scaled tension parameter of, e.g., scal = 6 lies in between these two extremes and is therefore suitable for initializing the minimization algorithm. 3.5. A minimization algorithm avoiding overfitting Consider the data points (xi , yi ), i = 1, . . . , m, with pseudo-arclength parametrization ti , i = 1, . . . , m. Within the frame of Algorithm 1, we minimize the error function equadr (Tn , Ln ) :=
m
wr2 {(xr − sx (tr ))2 + (yr − sy (tr ))2 }
(3.19)
r=1
with suitable inner knots 1 , . . . , n and interval tension parameters 1 , . . . , n+1 . In order to avoid overfitting, we force the distance between two consecutive knots i−1 and i to be at least as large as the distance of the two data points tli and tli +1 neighboring the middle of the knot interval [i−1 , i ]—see Fig. 4. Consider, without loss of generality, knots in ascending order. Defining li such that i−1 + i tli = max tl : tl , i = 1, . . . , n + 1, 2
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114
103
we get the nonlinear inequality constraints i − i−1 tli +1 − tli ,
i = 1, . . . , n + 1.
(3.20)
With this constraint, Eq. (3.2), which guarantees the uniqueness of the least-squares spline, is fulfilled. Moreover, the knots coincide with the data points if there are as many knots as data points—see Fig. 9. In shallow regions of the error function a large relative change in the scaled tension parameter causes only small changes of the exponential spline and, hence, of the error function. Therefore, in order to keep away from such shallow regions, we introduce the linear inequality constraint min i (i − i−1 ) max ,
i = 1, . . . , n + 1.
(3.21)
Finally, we keep the knots within the data support by the inequalities t0 i tm ,
i = 1, . . . , n.
(3.22)
3.5.1. Relative alteration rate of the tension parameters So far, we did not take into account that the knots and tension parameters play different roles and need not be of the same order of magnitude. One could predominantly focus either on finding good knot positions or tension parameters. In order to permit this choice we introduce a weighting factor c > 0 and define ˆ i := c i ,
i = 1, . . . , n + 1.
We minimize ˆ 1 ˆ n+1 ˆ ˆ eˆquadr (Tn ; 1 , . . . , n+1 ) := equadr Tn ; , . . . , c c
(3.23)
under the constraints (3.20), (3.22), and min
ˆ i (i − i−1 ) max . c
(3.24)
ˆ Thus, having, e.g., c > 1 the error term eˆquadr is flattened in -direction. This results in a reduced alteration rate of the tension parameters compared to the alteration rate of the knots, but this is true only for the same initial values. For the minimization problem the dynamics is more complicated. There is a kind of dynamical equilibrium between the optimization of the knot positions and of the tension parameters. Fairly good optimized knots lead to a reduced alteration rate of the knots and, hence, to a relatively large alteration rate of the tension parameters and vice versa. This is the reason why the choice of c is not crucial. We get similar results using c = 0.01 and 10 for adapting an exponential spline to titanium heat data, which were considered by de Boor [2]—see Fig. 6. In order to minimize Eq. (3.23) with constraints (3.20), (3.22), and (3.24) we use the function ‘fmincon’ for constrained minimization in Matlab’s optimization toolbox. 3.6. Fixing the number of the knots Starting with one interior knot, Algorithm 1 yields successive least-squares splines (sn,x (t), sn,y (t)), n = 1, 2, . . . , approximating m data points, but we still have to decide which is the optimal number of the knots. This problem is analogous to the choice of the optimal degree g of a least-squares polynomial Pg (x) [11,7]. If g is too small the corresponding fit is not accurate enough; if it is too large the results are overly affected by the statistical errors in the data values. The usual way to estimate the optimal degree g of a least-squares polynomial is to examine the root-mean-square residuals g , computed as the weighted sum of squared residuals divided by m − g − 1, i.e., the number of degrees of freedom m − dim(Pg ). In a satisfactory case, these g will decrease steadily as g increases and then settle down to a
104
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114
fairly constant value. Likewise, for our problem we examine the numbers n :=
equadr (Tn , Ln ) m − (n + 2)
and accept the number n of interior knots which first gives an approximately constant value for n . Alternatively, we set a threshold for the maximal distance 2 2 emax := max (xr − sx (tr )) + (yr − sy (tr )) r
(3.25)
(3.26)
between the data points (xr , yr ), r = 1, . . . , m, and the smoothing spline. Then, we stop the knot placing Algorithm 1 when emax is below the threshold. It depends on the application which of these two conditions is preferable. 4. Numerical examples Exponential splines are a generalization of cubic splines—the limiting case of zero tension leads to cubic splines. Some data allow the representation with cubic, as well as exponential splines. Fig. 5 shows exponential splines as an alternative to cubic splines. The simultaneous adaptation of knots and interval tension parameters is demonstrated in Fig. 6. Compared to cubic splines, the greater flexibility of exponential splines is combined with a lower tendency to overfitting. A comparison is shown for a representative example in Figs. 7 and 8. The limiting case for finding tension parameters for interpolating exponential splines is shown in Fig. 9. Finally, Figs. 10 and 11 show an example for which the parametric setting is necessary. Here, we use Algorithm 1 in order to find exponential least-squares splines for various data. Having exact data points, this approach leads to finding interpolating exponential splines as a special case. This includes an alternative to heuristics for finding the tension parameters, as proposed, e.g., in [32,25], or [12]. In Fig. 5, we show the optimized least-squares exponential splines with n = 1, . . . , 6 interior knots for titanium heat data [2], which are available, e.g., in Matlab. When adding the fourth knot n is only slightly decreased. This indicates that the fourth knot does not improve the result. On the other hand, adding a further knot at the peak of the data (lower left plot) reduces n , equadr , and emax . However, having only few data points at the data peak and no error estimate of the data, it is not possible to discriminate whether this leads to a more accurate representation of the data or to overfitting. The least-squares spline with three knots and four tension parameters (middle left) is comparable to the cubic least-squares splines for the same data with five and six knots shown in [7]. With the same number of parameters, cubic least-squares splines are somewhat more accurate but contain a higher tendency to undesirable oscillations. For this example we also investigate the influence of the weighting factor c —see Section 3.5.1. For the same initial values a small c leads to a large alteration rate of the tension parameters compared to the knot positions. In Fig. 6, we show the initial values and the first four iterations of the minimization algorithm for c = 10 and 0.01, respectively. In the first iteration, the scaled tension parameters are much more changed for c = 0.01 than for c = 10, but already in the second iteration the opposite is true for the third scaled tension parameter. In the fourth iteration the first, the second, and the third tension parameters show larger changing for c = 0.01, whereas the fourth one shows larger changing for c = 10. This demonstrates that, locally, there is a kind of dynamical equilibrium between the adaptation rate of the knot positions and the tension parameters. This explains why there is a large range of weighting factors c leading to the same result with comparable accuracy in comparable time—see Table 1. Fig. 7 shows another example found in [25]. Considering the values of n , equad , emax or the graphical appearance of the splines in correspondence to the data shows that at least six interior knots are desirable. Moreover, having 6 < 7 indicates that the seventh knot is not needed. Restricting the scaled tension parameters for all plots, apart from the bottom right one, to the interval [2, 20], we use the larger range [1, 40] for the bottom right plot. Using also stronger tension reduces the quadratic approximation error and leads to a somewhat different spline having six interior knots, too. For a comparison we also show the results of least-squares approximation for the same data with cubic splines in Fig. 8. Clearly, cubic splines are not suitable to form sharp corners, which leads to undesired spurious oscillations. Compared to exponential splines more knots have to be allocated, which further increases undesired oscillations in one part, whereas in other parts of the data knots are lacking—compare Figs. 7 and 8 (third row right), respectively.
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114 σ1= 0.086; equadr= 3.96; emax= 0.992
σ2= 0.0501; equadr= 2.25; emax= 0.667
2.5
2.5
2
2
1.5
1.5
1
1
0.5 600
700
800
900
1000
σ3= 0.000138; equadr= 0.00608; emax= 0.0405 2.5
0.5 600
2
1.5
1.5
1
1
700
800
900
1000
σ5= 7.57e-005; equadr= 0.00318; emax= 0.0288 2.5
0.5 600
1.5
1.5
1
1
800
900
1000
900
1000
700
800
900
1000
σ6= 7.74e-005; equadr= 0.00317; emax= 0.00289
2
700
800
2.5
2
0.5 600
700
σ4= 0.000137; equadr= 0.00591; emax= 0.0409 2.5
2
0.5 600
105
0.5 600
700
800
900
1000
Fig. 5. Finding successively optimal knots and tension parameters for titanium heat data (dots). The optimized least-squares splines are shown by solid lines and its knots are marked by squares. We start with one interior knot and two tension parameters for the two corresponding intervals in the upper left plot. Above each plot we print the summarized squared error equadr and the maximal error emax . Alternatively, the values n may be used for fixing the number of the knots.
106
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114 equadr= 0.614 (black), 0.509 (grey)
equadr= 0.79539
Iteration: 1
20
20
10
10
0
0 equadr= 0.553 (black), 0.508 (grey)
equadr= 0.538 (black), 0.416 (grey)
Iteration: 2
Iteration: 3
20
20
10
10
0
0 equadr= 0.5 (black), 0.389 (grey)
equadr= 0.0060803 (black/grey)
Iteration: 4
Iteration: 66
20
20
10
10
0
0
Fig. 6. We start finding the optimal positions and tension parameters for three interior knots as shown in the upper left plot. The interior knots are marked by squares, the scaled tension parameters, which are restricted to the range [2, 20], are marked by bars below the corresponding interval. The ticks in y-direction exclusively mark the strength of the scaled tension. The successive plots show the first four iterations of the minimization algorithm. For c = 0.01 (gray curve and bars) the minimizer more drastically changes the tension parameters than for c = 10 (black curve and bars). However, finally one reaches the same minimum (lower right plot).
In order to demonstrate the transition from least squares to interpolating splines, Fig. 9 shows the data proposed by Späth [32] as representative example. It is remarkable that already four interior knots lead to a very accurate approximating spline. With eight interior knots we have as many knots as data points and get an interpolating exponential spline in tension. The nonlinear local distance condition described in Section 3.5 successfully avoids overfitting and leads to an interpolating exponential spline with coinciding knots and data points when having as many knots as data points. Thus, we use the same algorithm for finding suitable tension parameters of least squares, as well as interpolating splines. Finally, we consider data presented by Jüttler [15]. Jüttler uses a reference curve in order to impose convexity constraints. This leads to shape preserving least-squares splines. Optimizing 42 degrees of freedom, he arrives at a shape preserving approximation with a least-squared error equadr = 9.5. Costantini et al. [3,4] use variable degree
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114 σ2= 0.229; equadr= 1.83; emax= 0.576
σ1= 0.231; equadr= 2.08; emax= 0.662 5
5
4
4
3
3
2
2
1
1
0
5
-10
-5
0
5
10
σ3= 0.176; equadr= 1.23; emax= 0.635
0
5
4
4
3
3
2
2
1
1
0
-10
-5
0
5
10
σ4= 0.176; equadr= 0.272; emax= 0.382
0 -10
-5
0
5
10
-10
5
4
4
3
3
2
2
1
1
-10
-5
0
5
10
-5
0
5
10
σ6= 0.00068; equadr= 0.00272; emax= 0.0263
σ5= 0.0247; equadr= 0.124; emax= 0.251 5
0
0
-10
-5
0
5
10
5
σ6= 0.000441; equadr= 0.00176; emax= 0.0263 5
4
4
3
3
2
2
1
1
σ7= 0.000888; equadr= 0.00266; emax= 0.029
0
107
-10
-5
0
5
10
0
-10
-5
0
5
10
Fig. 7. Least-squares exponential spline for Rentrop data (dots). The knots are marked by squares. For this example the larger range [1, 40] for scaled tension parameters gives a slightly better result (bottom right) than for the range [2, 20] (third row, right). All plots apart from the bottom right restrict the scaled tension parameters to [2, 20].
polynomial splines, which have similar characteristics as exponential splines. Considering the Jüttler data, they get a graphically comparable least-squares approximation—the least-squared error is not mentioned—with 15 intervals. This means 14 inner knots and 15 variable polynomial degrees for the 15 intervals, i.e., totally 29 degrees of freedom.
108
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114 σ2= 0.234; equadr= 1.88; emax= 0.559
σ1= 0.285; equadr= 2.57; emax= 0.694 5
5
4
4
3
3
2
2
1
1
0
0 -10
-5
0
5
10
-10
0
5
10
σ4= 0.0814; equadr= 0.488; emax= 0.352
σ3= 0.128; equadr= 0.899; emax= 0.42 5
5
4
4
3
3
2
2
1
1
0
-5
0 -10
-5
0
5
10
-10
σ5= 0.097; equadr= 0.485; emax= 0.35
0
5
10
σ6= 0.0678; equadr= 0.271; emax= 0.373
5
5
4
4
3
3
2
2
1
1
0
-5
0 -10
-5
0
5
10
-10
-5
0
5
10
σ7= 0.00564; equadr= 0.0169; emax= 0.0876 5
σ8= 0.00545; equadr= 0.0109; emax= 0.0636 5
4
4
3
3
2
2
1
1
0
0 -10
-5
0
5
10
-10
-5
0
5
10
Fig. 8. Least-squares cubic spline for the same data as in Fig. 7. The knots are marked by squares. Cubic splines show much more undesired oscillations than exponential splines—compare Fig. 7.
In Fig. 10, we use Algorithm 1 in order to find a least-squares exponential spline approximating the Jüttler data. Using eight interior knots and nine interval tension parameters, i.e., 17 degrees of freedom, we arrive at an approximation with a least-squared error equadr = 4.61—bottom right. Restricting the scaled tension parameter scal to 0.01, we use the same algorithm in order to find an approximating cubic spline. With eight interior knots, i.e., a total amount
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114 σ2= 0.209; equadr= 1.25; emax= 0.73
σ1= 0.426; equadr= 2.98; emax= 0.918 10
10
8
8
6
6
4
4
2
2
0
0
2
4
6
8
10
σ3= 0.0211; equadr= 0.106; emax= 0.253
0
10
8
8
6
6
4
4
2
2 0
2
4
6
8
10
σ5= 0.00103; equadr= 0.00308; emax= 0.0417
0
10
8
8
6
6
4
4
2
2 0
2
4
6
8
10
0
0
0
σ7= 7.53e-010; equadr= 7.53e-010; emax= 2.14e-005 10
8
8
6
6
4
4
2
2 0
2
4
6
8
10
4
6
8
10
2
4
6
8
10
2
4
6
8
10
equadr= 3.55e-030; emax= 1.26e-015
10
0
2
σ6= 2.18e-006; equadr= 4.36e-006; emax= 0.00318
10
0
0
σ4= 0.00199; equadr= 0.00794; emax= 0.0467
10
0
109
0
0
2
4
6
8
10
Fig. 9. Exponential least-squares spline for Späth data (dots). The knots are marked by squares. The local minimal distance condition prevents them from becoming closer than the data points. As many knots as data points lead to an interpolating exponential spline with suitable tension parameters.
110
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114 σ1= 18; equadr= 1080; emax= 7.28
σ2= 6.36; equadr= 375; emax= 4.86
10
10
5
5
0
0
-5
-5
-10 -5
10
0
5
10
15
σ3= 1.47; equadr= 85.1; emax= 2.17
-10 -5
10
5
5
0
0
-5
-5
-10 -5
0
5
10
15
-10 -5
10
5
5
0
0
-5
-5
0
5
10
15
-10 -5 10
5
5
0
0
-5
-5
0
5
10
10
15
σ4= 1.26; equadr= 71.5; emax= 2.26
0
5
10
15
0
5
10
15
σ8= 0.087; equadr= 4.61; emax= 0.667
σ7= 0.104; equadr= 5.6; emax= 0.683 10
-10 -5
5
σ6= 0.146; equadr= 8.02; emax= 0.709
σ5= 0.774; equadr= 43.4; emax= 2.05 10
-10 -5
0
15
-10 -5
0
5
10
15
Fig. 10. Least-squares exponential spline for Jüttler data (dots). The knots are marked by squares. Clearly, this example demands for the parametric setting.
of 8 degrees of freedom, we arrive at an approximating cubic spline with a least-squared error equadr = 6.43—see Fig. 11. Finally, Fig. 12 shows the curvature of the approximating exponential spline (left) with eight interior knots and the corresponding cubic spline (right). Negative curvature is dashed. Thus, there are exactly two inflection points for the exponential, as well as the cubic spline approximation like for the shape preserving approximation of Jüttler. The
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114 σ1= 18.5; equadr= 1111; emax= 7.24
σ2=10.9; equadr= 642; emax= 5.93
10
10
5
5
0
0
-5
-5
-10
-5
0
5
10
15
σ3= 2.05; equadr= 119; emax= 2.51
10
-10
5
0
0
-5
-5
-5
0
5
10
-5
15
-10
-5
10
5
5
0
0
-5
-5
-5
0
5
10
5
10
15
0
5
10
15
σ6= 0.17; equadr= 9.38; emax= 0.763
σ5= 1.32; equadr= 73.9; emax= 2.55 10
-10
0
σ4= 1.54; equadr= 88; emax= 2.39
10
5
-10
111
15
-10 -5
0
5
10
15
σ8= 0.121; equadr= 6.43; emax= 0.655
σ7= 0.123; equadr= 6.65; emax= 0.783 10
10
5
5
0
0
5
5
-10
-10 -5
0
5
10
15
-5
0
5
10
15
Fig. 11. Least-squares cubic spline for the same data as in Fig. 10. The knots are marked by squares. Due to the lack of sharp corners, the cubic spline is suitable.
approximating splines with less than eight interior knots shown in Figs. 10 and 11 also have not more than two inflection points. This means that at least for this example Algorithm 1 finds approximations without extraneous inflection points as desired. Note that the curvature of a parametric cubic spline (sx (t), sy (t)) need not be linear, actually it is in the 2D case sx (t)sy (t) − sy (t)sx (t). As expected, variations in the curvature are more pronounced for the exponential spline. Summarizing, we find approximations closer to the data which further demand for fewer degrees of freedom than the mentioned approaches which impose convexity constraints. Few degrees of freedom are preferable, because unnecessary
112
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114
Table 1 Performance of the minimization algorithm for three interior knots in dependence of the position-tension weighting factor c c
equadr
Iterations
100 50 10 1 0.1 0.01 0.001 0.0001 0.00001
0.00631891 0.00609905 0.00608029 0.00608025 0.00608025 0.00608024 0.00608028 0.00609962 No result
72 60 68 96 50 59 58 40 —
For the range c ∈ [0.001, 10] we get accurate minimization results. The number of necessary iterations is comparable. Outside this range we find inaccurate results or even worse numerical difficulties caused by very shallow surfaces used for minimization.
0.6 0.3 0.4 0.2 0.2 0.1 Curvature
Curvature
0
-0.2
-0.4
0
-0.1
-0.6 -0.2 -0.8 -0.3 -1 0
10
20
30
40 t
50
60
70
0
10
20
30
40
50
60
70
t
Fig. 12. (Left) Curvature of the exponential spline in Fig. 10 (bottom right). (Right) Curvature of the cubic spline in Fig. 11 (bottom right). Negative curvature is dashed. In both cases there are exactly two inflection points.
degrees of freedom are time consuming and can lead to undesired effects like spurious oscillations—unless avoided by convexity constraints. As demonstrated in Fig. 2 the knot positions are crucial. The curvature of the Jüttler data is moderate. Therefore, as demonstrated, cubic least-squares splines with locally optimal knots lead to an approximation close to the data with only 8 degrees of freedom. As undesired oscillations are related to overfitting, which is avoided by Algorithm 3.3, we can deal without convexity constraints.
5. Conclusion Least-squares splines in tension can be optimized by a simultaneous knot placing and tension setting algorithm. The dynamical equilibrium between the alteration rate of the knot positions and the tension parameters makes it easy to find
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114
113
a suitable relative alteration rate of the tension parameters compared to the knot positions. Scaled tension parameters are useful in order to keep the minimization away from shallow regions. Using a local minimal distance condition avoids overfitting and leads to interpolating splines with coincident knots and data points as a special case. This includes finding suitable tension parameters of interpolating exponential splines. The algorithm presented here usually does not violate natural convexity requirements inherited by the data. If such a case should occur, enlarging the tension of the critical interval(s) is a remedy. The simultaneous knot placing and tension setting approach presented here in combination with variable degree polynomial splines [4] is possible if the minimizer can deal with discrete polynomial degrees instead of smooth tension parameters. The generalization of our approach to 3D curves is straightforward.
Acknowledgements I wish to thank my academic advisor Prof. Dr. R. Seydel for proposing and supporting this work on exponential splines. Further, my thanks go to Dr. K.-D. Reinsch for motivating and profound discussions. Prof. Dr. B. Jüttler kindly provided the data used in Figs. 10 and 11. Last but not least, I wish to thank the referees for their supplements and improvements. References [1] B.A. Barsky, Exponential and polynomial methods for applying tension to an interpolating spline curve, Computer Vision, Graphics, Image Process. 27 (1984) 1–18. [2] C. de Boor, A Practical Guide to Splines, Springer, Berlin, Heidelberg, New York, 1978. [3] P. Costantini, F. Pelosi, Shape-preserving approximation by space curves, Numer. Algorithms 27 (2001) 237–264. [4] P. Costantini, F. Pelosi, Shape-preserving approximation of spatial data, Adv. Comput. Math. 20 (2004) 25–51. [5] R. Delbourgo, J.A. Gregory, C 2 rational quadratic spline interpolation to monotonic data, IMA J. Numer. Anal. 3 (1983) 141–152. [6] R. Delbourgo, J.A. Gregory, The determination of the derivative parameters for a monotonic rational quadratic interpolant, IMA J. Numer. Anal. 3 (1983) 141–152. [7] P. Dierckx, Curve and Surface Fitting with Splines, Oxford Science Publications, Clarendon Press, Oxford, 1993. [8] F.N. Fritsch, R.E. Carlson, Monotone piecewise cubic interpolation, SIAM, J. Numer. Anal. 17 (1980) 238–246. [9] J.A. Gregory, R. Delbourgo, Piecewise rational quadratic interpolation to monotonic data, IMA J. Numer. Anal. 2 (1982) 123–130. [10] A. Gregory, M. Sarfraz, A rational cubic spline with tension, Comput. Aided Geom. Design 7 (1990) 1–13. [11] J.G. Hayes, Curve fitting by polynomials in one variable, in: J.G. Hayes (Ed.), Numerical Approximation to Functions and Data, The Athlone Press, London, 1970, pp. 43–64. [12] U. Heidemann, Linearer Ausgleich mit Exponentialsplines bei automatischer Bestimmung der Intervallteilungspunkte, Computing 36 (1986) 217–227. [13] W. Heß, J.W. Schmidt, Convexity preserving interpolation with exponential splines, Computing 36 (1986) 335–342. [14] D.J. Higham, Monotonic piecewise cubic interpolation, with applications to ODE plotting, TR 229/90, Department of Computer Science, University of Toronto, 1990. [15] B. Jüttler, Shape preserving least-squares approximation by polynomial parametric spline curves, Comput. Aided Geom. Design 14 (1997) 731–747. [16] P.E. Koch, T. Lyche, Exponential B-splines in tension, in: C.K. Chui, L.L. Schumaker, J.D. Ward (Eds.), Approximation Theory VI, vol. 2, ISBN 0-12-174587-2, Academic Press, London, 1989, pp. 361–364. [17] P.E. Koch, T. Lyche, Construction of exponential tension B-splines of arbitrary order, in: P.J. Laurent, A. Le Méhauté, L.L. Schumaker (Eds.), Curves and Surfaces, ISBN 0-12-438660-1, Academic Press, London, 1991, pp. 255–258. [18] P.E. Koch, T. Lyche, Interpolation with exponential B-splines in tension, Computing Supplementum 8, in: G. Farin, H. Hagen, H. Noltemeier (Eds.), Geometric Modelling, ISBN 3-211-82399-9, Springer, Wien, New York, 1993, pp. 173–190. [19] R.W. Lynch, A method for choosing a tension factor for spline under tension interpolation, M.S. Thesis, University of Texas at Austin, 1982. [20] G.M. Nielson, Computation of -splines, TR 044-433-11, Department of Mathematics, Arizona State University, 1974. [21] G.M. Nielson, Some piecewise polynomial alternatives to splines under tension, in: R.E. Barnhill, R.F. Riesenfeld (Eds.), Computer Aided Geometric Design, 1974, pp. 209–235. [22] S. Pruess, Properties of splines in tension, J. Approx. Theory 17 (1976) 86–96. [23] R.J. Renka, Interpolatory tension splines with automatic selection of tension factors, SIAM J. Sci. Statist. Comput. 8 (3) (1987) 393–415. [24] R.J. Renka, Tspack: tension spline curve-fitting package, ACM Trans. Math. Software 19 (1993) 81–94. [25] P. Rentrop, An algorithm for the computation of the exponential spline, Numer. Math. 35 (1980) 81–93. [26] K.O. Riedel, Aspects of image processing—splines, anisotropic diffusion, and biological models, Ph.D. Thesis, ISBN 3-8322-1236-1, Shaker Verlag, Aachen, 2003. [27] K.O. Riedel, Two-dimensional splines on fairly arbitrary meshes, ZAMM—Z. Angew. Math. Mech. 85 (3) (2005) 176–188.
114
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114
[28] N.S. Sapidis, P.D. Kaklis, T.A. Loukakis, A method for computing the tension parameters in convexity-preserving spline-in-tension interpolation, Numer. Math. 54 (1988) 179–192. [29] M. Sarfraz, A geometrical rational spline with tension controls: an alternative to the weighted nu-spline, Punjab Univ. J. Math. 26 (1993) 27–40. [30] M. Sarfraz, Freeform rational bicubic spline surfaces with tension control, Facta Univ. Ser. Math. Inform. 9 (1994) 83–93. [31] D.G. Schweikert, An interpolating curve using a spline in tension, J. Math. Phys. 45 (1966) 312–317. [32] H. Späth, Exponential spline interpolation, Computing 4 (1968) 225–233. [33] J. Stoer, Numerische Mathematik 1, Springer, Berlin, Heidelberg, New York, 1999. [34] J. Stoer, R. Bulirsch, Introduction to Numerical Analysis, Springer, Berlin, Heidelberg, New York, 1980.
Locally optimal knots and tension parameters for exponential splines Karl O. Riedela, b,∗ a Mathematical Institute, University of Cologne, Weyertal 86-90, 50931 Cologne, Germany bWeißensteinerstraße 74, 73525 Schwäbisch Gmünd, Germany
Received 24 August 2004; received in revised form 11 August 2005
Abstract Tension can be applied to cubic splines in order to avoid undesired spurious oscillations. This leads to the well-known (exponential) spline in tension. It is crucial but unfortunately difficult to find suitable tension parameters of interpolating splines in tension. Instead of heuristics, we propose a simultaneous knot placing and tension setting algorithm for least-squares splines in tension which includes interpolating splines in tension as a special case. Moreover, the splines presented here are the foundation of exponential surface splines on fairly arbitrary meshes [K.O. Riedel, Two-dimensional splines on fairly arbitrary meshes, ZAMM—Z. Angew. Math. Mech. 85(3) (2005) 176–188]. © 2005 Elsevier B.V. All rights reserved. MSC: 65D07; 65D10; 65D05 Keywords: Spline; Tension; Exponential spline; Interpolation; Approximation
1. Introduction Cubic splines, i.e. twice continuously differentiable functions consisting of piecewise cubic polynomials, are widespread. They are used to interpolate or approximate data. Interpolating cubic splines definitely show less oscillations than interpolating polynomials. Nevertheless, in some cases spurious oscillations of the splines occur which are not indicated by the data. These oscillations can be reduced by shortening the length of the spline. To this purpose, an additional tension parameter is introduced. This leads to the (exponential) spline in tension, i.e., a twice continuously differentiable function consisting piecewise of a linear combination of the functions 1, x, exp(i x), and exp(−i x), where i > 0 is the interval tension parameter. Some authors call this kind of spline in tension just exponential spline, whereas others understand by an exponential spline a generalization or a different spline with exponential terms. Bsplines for splines in tension have been introduced by [16], see also [17,18]. Explicit formulas also can be found in [26]. We will use them in order to calculate smoothing (least-squares) splines in tension. Often convexity constraints are given with or imposed to interpolating or approximating splines. Especially, the tension parameters of tension splines can be suitably adapted in order to fulfill these shape constraints. Yet, in this ∗ Corresponding address: Weißensteinerstraße 74, 73525 Schwäbisch Gmünd, Germany. Tel.:+49(0)1706837882; fax: +49(0)6990553820.
E-mail address: [email protected]. 0377-0427/$ - see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2005.08.012
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114
95
paper, we do not deal with convexity constraints, but simultaneously choose knot positions and tension parameters which minimize the least-squared error. This yields visually pleasing curves. Smooth curves as well as sharp corners are captured by a minimal number of knots, which implies a highly reduced tendency to spurious oscillations, i.e., a strong tendency to fulfill convexity constraints. This is an alternative at least to imposing convexity constraints. Finding suitable knot positions of a least-squares spline is important but unfortunately difficult. Dierckx [7] has proposed a successive knot placing algorithm which minimizes the least-squared error as a function of the knot positions. We generalize this approach considering the dependence on the interval tension parameters of exponential splines simultaneously. Further, we introduce a local minimal distance condition in order to avoid overfitting. This, automatically, leads to the interpolating spline if the number of knots and data points is identical. Thus, the same algorithm can be used in order to find interpolating as well as least-squares splines. Considering interpolating splines, this includes the generation of suitable tension parameters which have been determined by heuristic methods in former approaches. This paper is organized as follows: In Section 2, we discuss the problem of spurious oscillations. Splines in tension, and for data without sharp corners also cubic splines with suitable knots are a remedy. The simultaneous knot placing and tension setting algorithm is described in Section 3. Finally, we demonstrate its efficiency in finding suitable leastsquares splines for various data in Section 4. Interpolating splines are automatically included. 1.1. Related work The spline in tension was first analytically modeled by Schweikert [31]. He further defines an extraneous inflection point as a zero crossing of the second derivative which is not motivated by the data and shows that the exponential spline is free of extraneous inflection points for a tension large enough. Späth [32] simplifies this proof. Thus, finding tension parameters for a spline without spurious oscillations in principle is possible. One choice is infinite tension leading to a polygonal line. However, the task of finding tension parameters for a smooth, visually pleasant spline without extraneous inflection points is neither uniquely solvable nor trivial. Späth [32] and Rentrop [25] propose heuristic methods in order to determine tension parameters for which the exponential spline is free of extraneous inflection points. Heß [13] and Renka [23] also describe methods for finding suitable tension parameters of interpolating splines. Heidemann [12] uses the tension parameters proposed in [25] for least-squares splines. Lynch [19] presents an iterative scheme for constructing shape preserving exponential splines under the stringent constraint of uniform tension. Sapidis et al. [28] propose a generalized Newton Raphson method for finding suitable tension parameters. Pruess [22] shows that the approximation of sufficiently smooth functions with exponential splines is of fourth order, as for cubic splines. For large tension parameters the exponential spline is, essentially, locally linear. He exploits this fact in order to construct convex or monotone spline approximants. Nielson [20,21] proposed a polynomial alternative to the exponential spline, which he called the -spline. For a rational spline with interval and point tension see [10,29,30]. For a derivation of cubic and exponential splines from a variational principle and a comparison see [1]. Related work for monotony preserving interpolation is found, e.g., in [8,9,5,6,14]. For a software package on exponential splines see [24]. Jüttler [15] proposes the use of a reference curve in order to impose convexity constraints, i.e., given inflection points, to the approximating spline. Costantini et al. [3,4] use variable degree polynomial splines instead of exponential splines in order to get shape-preserving approximants for 2D and 3D parametric curves.
2. Suitable knots and tension parameters: an example We consider an example of Späth [32]. Only the interpolating cubic spline shows spurious oscillations which disappear for the exponential spline with uniform tension 0 = 7 as shown in Fig. 1. Using suitable knots, there are also cubic splines without spurious oscillations—see Fig. 2. We found these knots by trial and error by hand. Unfortunately, it is not easy to find optimal knot positions automatically. Moreover, in order to form sharp corners with cubic splines a good placement of the knots is difficult and could be impossible. For an alternative, we refer to Section 3.3, where we use least-squares exponential splines with as few as possible knots and optimize the interval tension parameters and the knot positions simultaneously.
96
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114 10 9 8 7 6 5 4 3
0
1
2
3
4
5
6
7
8
9
10
0
1
2
3
4
5
6
7
8
9
10
0
1
2
3
4
5
6
7
8
9
10
1 0.5 0
1 0.5 0
Fig. 1. (Top) Späth data (stars) are interpolated by a cubic spline (dashed) and an exponential spline (full line) with uniform tension 0 = 7. Natural boundary conditions are used for both splines. The corresponding exponential B-splines (bottom) are intermediate between cubic B-splines (middle) and hat functions, the basis functions of piecewise linear splines.
10
10
8
8
6
6
4
4 0
2
4
6
8
10
0
2
4
6
8
10
Fig. 2. Cubic splines with natural boundary conditions interpolating Späth data (stars). The knots of the splines are depicted by squares. Changing the position of just one knot marked by arrows leads to strongly increased oscillations.
3. Data approximation by splines in tension 3.1. Least-squares splines Dealing with noisy data it does not make sense to use an interpolating spline. Like for interpolation, an approximating spline is preferable to a polynomial, because it shows less oscillations. In this section, we state how to find the smoothing spline which minimizes the least-squared error for given knots and tension parameters. For finding suitable knots and tension parameters see Section 3.3. Given the data points (tr , xr ),
r = 1, . . . , m,
we use at most the same number of B-splines with natural boundary conditions Bi (t),
i = 0, . . . , n + 1, n + 2 m,
(3.1)
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114
97
with the knots 0 , . . . , n+1 which are placed such that there exist n data points tri ∈ {t2 , . . . , tm−1 }, i = 1, . . . , n, which fulfill 0 = t1 ,
i−1 tri i , i = 1, . . . , n,
n+1 = tm .
(3.2)
In order to find the B-spline representation of the approximating spline, we search the coefficients ci , i = 0, . . . , n + 1, which minimize 2 m n+1 2 wr xr − ci Bi (tr ) = min (3.3) r=1
i=0
for given weighting factors wr > 0, r =1, . . . , m. For the numerical examples presented here we use wr =1, r =1, . . . , m. Denoting ⎞ ⎞ ⎛ ⎛ w1 x1 c0 . . c = ⎝ .. ⎠ , xw = ⎝ .. ⎠ , E = (wr Bi (tr ))(r,i) , cn+1 w m xm for r = 1, . . . , m, i = 0, . . . , n + 1, we search the coefficient vector c of the least-squares spline which minimizes min Ec − xw 22 .
(3.4)
c
Lemma 1. Given data points and knots fulfilling Eq. (3.2), the least-squares spline which minimizes (3.3) is unique. Proof. The least-squares spline can be defined by its coefficient vector c. In order to show that c is unique we have to show that E has full rank (see, e.g., [34]). Thus, it is sufficient to show Ec = 0 implies ci = 0, i = 0, . . . , n + 1. Having r > 0, r = 1, . . . , m, we can set without loss of generality r = 1, r = 1, . . . , m. Then Ec = 0 can be written as n+1
ci Bi (tr ) = 0,
r = 1, . . . , m.
i=0
Denoting s(t) :=
n+1
ci Bi (t)
i=0
this reads s(tr ) = 0,
r = 1, . . . , m.
Using the notations tr0 := t1 and trn+1 := tm this implies s(tri ) = 0,
i = 0, . . . , n + 1 where i−1 tri i , i = 1, . . . , n.
The spline s(t) is twice continuously differentiable and the theorem of Rolle yields n + 1 points i with s (i ) = 0,
tri−1 < i < tri , i = 1, . . . , n + 1.
Applying Rolle again we have n points i fulfilling s (i ) = 0,
i < i < i+1 , i = 1, . . . , n.
(3.5)
98
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114
Setting 0 := tr0 and n+1 := trn+1 and considering natural boundary conditions we have s (i ) = 0,
i = 0, . . . , n + 1,
i.e., n + 2 zeros i within n + 1 intervals [i , i+1 ], i = 0, . . . , n. Hence, there exists an interval J := [j , j +1 ] which contains two zeros k and k+1 . Considering Eq. (3.5) we have k < k+1 < k+1 , i.e., k+1 ∈ J . Further, Eq. (3.2) guarantees the existence of a data point til ∈ J . Summarizing, we have s(til ) = 0,
til ∈ J ,
s (k+1 ) = 0, s (k ) = 0,
(3.6)
k+1 ∈ J ,
(3.7)
k ∈ J ,
s (k+1 ) = 0,
(3.8)
k+1 ∈ J .
(3.9)
Within interval J for a given tension parameter > 0 the exponential spline s(t) can be written as s(t) = a + bt + cet + de−t , s (t) = b + cet − de−t , s (t) = c2 et + d2 e−t ,
t ∈ J, t ∈ J,
t ∈ J.
Hence, Eq. (3.8) reads −k c 2
e k +d 2 e
=0, >0
(3.10)
>0
which yields cd 0.
(3.11)
Subtraction of Eq. (3.8) from (3.9) gives c 2 (ek+1 − ek ) +d 2 (e−k+1 − e−k ) =0,
>0
<0
which implies cd 0.
(3.12)
Considering Eqs. (3.11) and (3.12), we have cd = 0. Hence, Eq. (3.10) implies c = d = 0. Hence, (3.7) reads b=0 and (3.6) reads a = 0. Thus, we have s(t) ≡ 0,
t ∈ [j , j +1 ].
With the tension parameter ˜ > 0 on the next interval [j +1 , j +2 ] we write s(t) as ˜ ˜ − j +1 ) + ce ˜ −˜ (t−j +1 ) . s(t) = a˜ + b(t ˜ (t−j +1 ) + de
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114
99
Considering that s is twice continuously differentiable yields 2 ˜ = 0, s (j +1 ) = ˜ (c˜ + d)
i.e., ˜ c˜ = −d. Hence, s(j +1 ) = a˜ = 0 and s (j +1 ) = b˜ + 2c˜˜ = 0, i.e., ˜ b˜ = −2c˜. Eq. (3.2) guarantees another zero trj +2 ∈ (j +1 , j +2 ], i.e., ˜ r − j +1 ) + ce ˜ s(trj +2 ) = −2c˜(t j +2
˜ (trj +2 −j +1 )
− ce ˜
−˜ (trj +2 −j +1 )
= 0,
which reads ˜ r − j +1 )] =0. ˜ r − j +1 )) − (t c˜ [sinh((t j +2 j +2
>0
This implies c˜ = 0
and b˜ = 0.
Thus, we have s(t) ≡ 0,
t ∈ [j +1 , j +2 ].
Iteration also to the left-hand side yields s(t) =
n+1
ci Bi (t) ≡ 0,
t ∈ [t1 , tm ].
(3.13)
i=0
As the B-splines Bi (t), i = 0, . . . , n + 1 are linearly independent on [t1 , tm ], Eq. (3.13) implies ci = 0, i = 0, . . . , n + 1. For calculating the coefficients c of the least-squares spline by minimizing (3.4) see, e.g., [33]. 3.2. Pseudo-arclength parametrization So far, in the functional setting the abscissae and ordinates of the data points were not treated equally. Only functions x(t) could be represented. For many applications this is appropriate. Yet, in order to be able to include also parametric curves like, e.g., the exponential spline representation of contours of camera images [27], we use the parametric setting (x(t), y(t)). Moreover, generalizing rectangular meshes, exponential surface splines are introduced in [27] on fairly arbitrary meshes, which are represented by the more general parametric exponential splines presented here. Using parametric splines, Eq. (3.4) is to be solved twice, whereas the computational effort for calculating the B-splines is not increased. If computational time does not matter, also data with a functional dependency can be represented by parametric splines without changing the code, provided that the function x = x(t) = sx (t) is invertible. Usually, the
100
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114
monotonicity of the data (ti , xi ) leads to a monotonous, and hence invertible, spline sx (t). In such a situation, parametric splines look much the same as their functional counterparts. This is why we do not explicitly investigate functional splines, although only the data presented in Figs. 10 and 11 demand for the parametric setting. If the functional setting is required explicitly the code can be adapted by minor changes, i.e., using the given abscissae instead of pseudo-arclength parametrization and leaving out the calculation of the second least-squares spline sy (t). Parametric spline curves suggest the use of pseudo-arclength parametrization. Thus, starting from the data points r = 1, . . . , m,
(xr , yr ),
(3.14)
we introduce t1 := 0, r tr := (xk − xk−1 )2 + (yk − yk−1 )2 ,
r = 2, . . . , m.
(3.15)
k=2
For given tension parameters i > 0, i = 1, . . . , n + 1, and joint knots i , i = 0, . . . , n + 1, in the interval [t1 , tm ] we calculate joint B-splines in tension [16,18,26] Bi (t),
t ∈ [t1 , tm ], i = 0, . . . , n + 1.
Now, the problem splits into two parts: for each partial set of data points (tr , xr ), (tr , yr ),
r = 1, . . . , m, r = 1, . . . , m,
we find the coefficients of the least-squares spline according to Section 3.1. Denote them by cx,i and cy,i , respectively. Hence, the least-squares spline with pseudo-arclength parametrization is given by n+1 n+1 sx (t) := cx,i Bi (t), sy (t) := cy,i Bi (t) , t ∈ [t1 , tm ]. (3.16) i=0
i=0
This treatment includes interpolating splines, if there are as many knots as data points. 3.3. Finding optimal knot positions and tension parameters As demonstrated in Section 2 finding suitable knots and tension parameters is difficult. We use the knot placing algorithm of Dierckx [7] for successively placing knots. However, we do not only adapt the knots of (cubic) splines, but as well the tension parameters of exponential splines simultaneously. We further introduce a local minimal distance condition, which insures that the distance between the knots is not smaller than the distance between the nearest data points. Thus, we avoid oscillations caused by overfitting. We always use the two outermost data points t1 and tm as knots 0 and n+1 . For given interior knots Tn := {1 , . . . , n } and interval tension parameters Ln := {1 , . . . , n+1 } we consider equadr (Tn , Ln ) :=
m
wr2 {(xr − sx (tr ))2 + (yr − sy (tr ))2 }
(3.17)
r=1
as a measure for the discrepancy between the data points (xr , yr ), r = 1, . . . , m, and the least-squares spline (sn,x (t), sn,y (t)), t ∈ [t1 , tm ], with n interior knots. In order to get tr , r = 1, . . . , m, we use pseudo-arclength parametrization according to Section 3.2.
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114
101
3.3.1. Initial values for the knots and tension parameters The success of a free knot least-squares algorithm for minimizing equadr (Tn , Ln ) will strongly depend on the choice of the initial values, i.e., the initial interior knots T0n and the initial tension parameters L0n . Indeed, due to the nonconvexity of equadr (Tn , Ln ) we can only expect to find a local minimum in the vicinity of T0n , L0n . Therefore, in order to find the global minimum, the determination of good initial knots and interval tension parameters is an important but unfortunately also a non-trivial task. A straightforward equidistant distribution of knots, for example, will rarely yield the global minimum of equadr (Tn , Ln ). Furthermore, we do not know in advance the required number of knots for an acceptable least-squares fit. Therefore, we use a scheme which will hopefully return the successive least-squares splines (sn,x (t), sn,y (t)), n = 1, 2, . . . , and which uses the locally optimal knot/tension set T∗n , L∗n in order to find suitable initial values T0n+1 , L0n+1 for determining (sn+1,x (t), sn+1,y (t)): Algorithm 1 (Successive knot placing strategy). • Initialization: Choose a scaled (see Section 3.4) tension parameter, e.g., med = 6, leading to a spline far apart from a cubic spline as well as from a polygonal line: 01 =
t1 + t m , 2
01 =
med , 01 − t1
02 =
med . tm − 01
• For n = 1, 2, . . . until the approximating spline is good enough: ◦ Choose the allowed interval [min , max ] for the scaled tension parameters. ◦ Starting with {Tn , Ln } determine the optimal knot positions and tension parameters T∗n , L∗n and the corresponding least-squares spline (sn,x (t), sn,y (t)). For this, we choose a minimization algorithm which limits the interval tension parameters to i ∈ [min /(i − i−1 ), max /(i − i−1 )], i = 1, . . . , n + 1, and forces the distance between the knots to be at least as large as the distance of the neighboring data points. ◦ Compute the relative error measure numbers i , i = 1, . . . , n + 1:
qi +mi 2 2 2 r=qi +1 wr ([xr − sn,x (tr )] + [yr − sn,y (tr )] ) i = , m − mi tqi < ∗i−1 tqi +1 < tqi +2 < · · · < tqi +mi ∗i < tqi +mi +1 , ∗0 := t0 ,
∗n+1 := tm .
◦ Find suitable initial knot positions and tension parameters in order to determine (sn+1,x (t), sn+1,y (t)): Let l denote the interval with the largest relative error measure l = max{1 , . . . , n+1 } and define for this interval the new tension parameter ∗new := max{∗l , 2min /(∗l − ∗l−1 )}. This leads to 0i
=
0i
=
i = 1, . . . , l − 1, ∗i , (∗l−1 + ∗l )/2, i = l, ∗i−1 , i = l + 1, . . . , n + 1, i = 1, . . . , l − 1, ∗i , ∗new , i = l, l + 1, ∗i−1 , i = l + 2, . . . , n + 2.
◦ Stop when the approximation is good enough. So, this algorithm successively locates the initial position of an additional knot in the middle of the knot interval [l−1 , l ], i.e., in a part of [t1 , tm ] where on the one hand (sn,x (t), sn,y (t)) has a large weighted sum of squared residuals, but where on the other hand the concentration of knots is not too high (m − ml not too large). Within the frame of Algorithm 1, we proceed with describing scaled tension parameters (Section 3.4) and a minimization algorithm for finding optimal knots which are not closer than the data points (Section 3.5). In this way, we avoid overfitting. Finally, we address the task of fixing the number of the knots (Section 3.6).
102
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114 1 λscal=20 λscal=6 λscal=2
0.8 0.6
cubic
0.4 0.2 0
0
1
2
3
4
5
Fig. 3. Exponential B-splines with different tension ranging from cubic splines to hat functions. Scaled tension parameters ranging from 2 to 20 include a quite large portion of possible exponential splines.
tli o
o
*o τi-1
o
tli+1 x o
o
τi-1 + τi 2
*o
o
τi
Fig. 4. Fixing locally the minimal knot distance.
3.4. Scaled tension parameters The shape of an exponential spline depends on the interval tension parameters as well as on the distance between the data points. This means that interpolating or least-squares splines are not size invariant. Linear compression of the data points leads to larger squared second derivatives and hence to a spline with a more pronounced cubic behavior. This can be overcome by using scaled interval tension parameters i,scal = i (i − i−1 ),
i = 1, . . . , n + 1,
(3.18)
for the knots 0 , . . . , n+1 . Scaled tension parameters yield size invariant splines. In other words, stretching the data points by a factor 2 and dividing the tension parameters by 2 leads to the exponential spline stretched by a factor 2. For a proof see [26]. Fig. 3 shows exponential B-splines with different scaled tension parameters. Consider that scaled tension parameters of 2 and 20 lead to a spline quite close to a cubic spline and a hat function, respectively. A scaled tension parameter of, e.g., scal = 6 lies in between these two extremes and is therefore suitable for initializing the minimization algorithm. 3.5. A minimization algorithm avoiding overfitting Consider the data points (xi , yi ), i = 1, . . . , m, with pseudo-arclength parametrization ti , i = 1, . . . , m. Within the frame of Algorithm 1, we minimize the error function equadr (Tn , Ln ) :=
m
wr2 {(xr − sx (tr ))2 + (yr − sy (tr ))2 }
(3.19)
r=1
with suitable inner knots 1 , . . . , n and interval tension parameters 1 , . . . , n+1 . In order to avoid overfitting, we force the distance between two consecutive knots i−1 and i to be at least as large as the distance of the two data points tli and tli +1 neighboring the middle of the knot interval [i−1 , i ]—see Fig. 4. Consider, without loss of generality, knots in ascending order. Defining li such that i−1 + i tli = max tl : tl , i = 1, . . . , n + 1, 2
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114
103
we get the nonlinear inequality constraints i − i−1 tli +1 − tli ,
i = 1, . . . , n + 1.
(3.20)
With this constraint, Eq. (3.2), which guarantees the uniqueness of the least-squares spline, is fulfilled. Moreover, the knots coincide with the data points if there are as many knots as data points—see Fig. 9. In shallow regions of the error function a large relative change in the scaled tension parameter causes only small changes of the exponential spline and, hence, of the error function. Therefore, in order to keep away from such shallow regions, we introduce the linear inequality constraint min i (i − i−1 ) max ,
i = 1, . . . , n + 1.
(3.21)
Finally, we keep the knots within the data support by the inequalities t0 i tm ,
i = 1, . . . , n.
(3.22)
3.5.1. Relative alteration rate of the tension parameters So far, we did not take into account that the knots and tension parameters play different roles and need not be of the same order of magnitude. One could predominantly focus either on finding good knot positions or tension parameters. In order to permit this choice we introduce a weighting factor c > 0 and define ˆ i := c i ,
i = 1, . . . , n + 1.
We minimize ˆ 1 ˆ n+1 ˆ ˆ eˆquadr (Tn ; 1 , . . . , n+1 ) := equadr Tn ; , . . . , c c
(3.23)
under the constraints (3.20), (3.22), and min
ˆ i (i − i−1 ) max . c
(3.24)
ˆ Thus, having, e.g., c > 1 the error term eˆquadr is flattened in -direction. This results in a reduced alteration rate of the tension parameters compared to the alteration rate of the knots, but this is true only for the same initial values. For the minimization problem the dynamics is more complicated. There is a kind of dynamical equilibrium between the optimization of the knot positions and of the tension parameters. Fairly good optimized knots lead to a reduced alteration rate of the knots and, hence, to a relatively large alteration rate of the tension parameters and vice versa. This is the reason why the choice of c is not crucial. We get similar results using c = 0.01 and 10 for adapting an exponential spline to titanium heat data, which were considered by de Boor [2]—see Fig. 6. In order to minimize Eq. (3.23) with constraints (3.20), (3.22), and (3.24) we use the function ‘fmincon’ for constrained minimization in Matlab’s optimization toolbox. 3.6. Fixing the number of the knots Starting with one interior knot, Algorithm 1 yields successive least-squares splines (sn,x (t), sn,y (t)), n = 1, 2, . . . , approximating m data points, but we still have to decide which is the optimal number of the knots. This problem is analogous to the choice of the optimal degree g of a least-squares polynomial Pg (x) [11,7]. If g is too small the corresponding fit is not accurate enough; if it is too large the results are overly affected by the statistical errors in the data values. The usual way to estimate the optimal degree g of a least-squares polynomial is to examine the root-mean-square residuals g , computed as the weighted sum of squared residuals divided by m − g − 1, i.e., the number of degrees of freedom m − dim(Pg ). In a satisfactory case, these g will decrease steadily as g increases and then settle down to a
104
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114
fairly constant value. Likewise, for our problem we examine the numbers n :=
equadr (Tn , Ln ) m − (n + 2)
and accept the number n of interior knots which first gives an approximately constant value for n . Alternatively, we set a threshold for the maximal distance 2 2 emax := max (xr − sx (tr )) + (yr − sy (tr )) r
(3.25)
(3.26)
between the data points (xr , yr ), r = 1, . . . , m, and the smoothing spline. Then, we stop the knot placing Algorithm 1 when emax is below the threshold. It depends on the application which of these two conditions is preferable. 4. Numerical examples Exponential splines are a generalization of cubic splines—the limiting case of zero tension leads to cubic splines. Some data allow the representation with cubic, as well as exponential splines. Fig. 5 shows exponential splines as an alternative to cubic splines. The simultaneous adaptation of knots and interval tension parameters is demonstrated in Fig. 6. Compared to cubic splines, the greater flexibility of exponential splines is combined with a lower tendency to overfitting. A comparison is shown for a representative example in Figs. 7 and 8. The limiting case for finding tension parameters for interpolating exponential splines is shown in Fig. 9. Finally, Figs. 10 and 11 show an example for which the parametric setting is necessary. Here, we use Algorithm 1 in order to find exponential least-squares splines for various data. Having exact data points, this approach leads to finding interpolating exponential splines as a special case. This includes an alternative to heuristics for finding the tension parameters, as proposed, e.g., in [32,25], or [12]. In Fig. 5, we show the optimized least-squares exponential splines with n = 1, . . . , 6 interior knots for titanium heat data [2], which are available, e.g., in Matlab. When adding the fourth knot n is only slightly decreased. This indicates that the fourth knot does not improve the result. On the other hand, adding a further knot at the peak of the data (lower left plot) reduces n , equadr , and emax . However, having only few data points at the data peak and no error estimate of the data, it is not possible to discriminate whether this leads to a more accurate representation of the data or to overfitting. The least-squares spline with three knots and four tension parameters (middle left) is comparable to the cubic least-squares splines for the same data with five and six knots shown in [7]. With the same number of parameters, cubic least-squares splines are somewhat more accurate but contain a higher tendency to undesirable oscillations. For this example we also investigate the influence of the weighting factor c —see Section 3.5.1. For the same initial values a small c leads to a large alteration rate of the tension parameters compared to the knot positions. In Fig. 6, we show the initial values and the first four iterations of the minimization algorithm for c = 10 and 0.01, respectively. In the first iteration, the scaled tension parameters are much more changed for c = 0.01 than for c = 10, but already in the second iteration the opposite is true for the third scaled tension parameter. In the fourth iteration the first, the second, and the third tension parameters show larger changing for c = 0.01, whereas the fourth one shows larger changing for c = 10. This demonstrates that, locally, there is a kind of dynamical equilibrium between the adaptation rate of the knot positions and the tension parameters. This explains why there is a large range of weighting factors c leading to the same result with comparable accuracy in comparable time—see Table 1. Fig. 7 shows another example found in [25]. Considering the values of n , equad , emax or the graphical appearance of the splines in correspondence to the data shows that at least six interior knots are desirable. Moreover, having 6 < 7 indicates that the seventh knot is not needed. Restricting the scaled tension parameters for all plots, apart from the bottom right one, to the interval [2, 20], we use the larger range [1, 40] for the bottom right plot. Using also stronger tension reduces the quadratic approximation error and leads to a somewhat different spline having six interior knots, too. For a comparison we also show the results of least-squares approximation for the same data with cubic splines in Fig. 8. Clearly, cubic splines are not suitable to form sharp corners, which leads to undesired spurious oscillations. Compared to exponential splines more knots have to be allocated, which further increases undesired oscillations in one part, whereas in other parts of the data knots are lacking—compare Figs. 7 and 8 (third row right), respectively.
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114 σ1= 0.086; equadr= 3.96; emax= 0.992
σ2= 0.0501; equadr= 2.25; emax= 0.667
2.5
2.5
2
2
1.5
1.5
1
1
0.5 600
700
800
900
1000
σ3= 0.000138; equadr= 0.00608; emax= 0.0405 2.5
0.5 600
2
1.5
1.5
1
1
700
800
900
1000
σ5= 7.57e-005; equadr= 0.00318; emax= 0.0288 2.5
0.5 600
1.5
1.5
1
1
800
900
1000
900
1000
700
800
900
1000
σ6= 7.74e-005; equadr= 0.00317; emax= 0.00289
2
700
800
2.5
2
0.5 600
700
σ4= 0.000137; equadr= 0.00591; emax= 0.0409 2.5
2
0.5 600
105
0.5 600
700
800
900
1000
Fig. 5. Finding successively optimal knots and tension parameters for titanium heat data (dots). The optimized least-squares splines are shown by solid lines and its knots are marked by squares. We start with one interior knot and two tension parameters for the two corresponding intervals in the upper left plot. Above each plot we print the summarized squared error equadr and the maximal error emax . Alternatively, the values n may be used for fixing the number of the knots.
106
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114 equadr= 0.614 (black), 0.509 (grey)
equadr= 0.79539
Iteration: 1
20
20
10
10
0
0 equadr= 0.553 (black), 0.508 (grey)
equadr= 0.538 (black), 0.416 (grey)
Iteration: 2
Iteration: 3
20
20
10
10
0
0 equadr= 0.5 (black), 0.389 (grey)
equadr= 0.0060803 (black/grey)
Iteration: 4
Iteration: 66
20
20
10
10
0
0
Fig. 6. We start finding the optimal positions and tension parameters for three interior knots as shown in the upper left plot. The interior knots are marked by squares, the scaled tension parameters, which are restricted to the range [2, 20], are marked by bars below the corresponding interval. The ticks in y-direction exclusively mark the strength of the scaled tension. The successive plots show the first four iterations of the minimization algorithm. For c = 0.01 (gray curve and bars) the minimizer more drastically changes the tension parameters than for c = 10 (black curve and bars). However, finally one reaches the same minimum (lower right plot).
In order to demonstrate the transition from least squares to interpolating splines, Fig. 9 shows the data proposed by Späth [32] as representative example. It is remarkable that already four interior knots lead to a very accurate approximating spline. With eight interior knots we have as many knots as data points and get an interpolating exponential spline in tension. The nonlinear local distance condition described in Section 3.5 successfully avoids overfitting and leads to an interpolating exponential spline with coinciding knots and data points when having as many knots as data points. Thus, we use the same algorithm for finding suitable tension parameters of least squares, as well as interpolating splines. Finally, we consider data presented by Jüttler [15]. Jüttler uses a reference curve in order to impose convexity constraints. This leads to shape preserving least-squares splines. Optimizing 42 degrees of freedom, he arrives at a shape preserving approximation with a least-squared error equadr = 9.5. Costantini et al. [3,4] use variable degree
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114 σ2= 0.229; equadr= 1.83; emax= 0.576
σ1= 0.231; equadr= 2.08; emax= 0.662 5
5
4
4
3
3
2
2
1
1
0
5
-10
-5
0
5
10
σ3= 0.176; equadr= 1.23; emax= 0.635
0
5
4
4
3
3
2
2
1
1
0
-10
-5
0
5
10
σ4= 0.176; equadr= 0.272; emax= 0.382
0 -10
-5
0
5
10
-10
5
4
4
3
3
2
2
1
1
-10
-5
0
5
10
-5
0
5
10
σ6= 0.00068; equadr= 0.00272; emax= 0.0263
σ5= 0.0247; equadr= 0.124; emax= 0.251 5
0
0
-10
-5
0
5
10
5
σ6= 0.000441; equadr= 0.00176; emax= 0.0263 5
4
4
3
3
2
2
1
1
σ7= 0.000888; equadr= 0.00266; emax= 0.029
0
107
-10
-5
0
5
10
0
-10
-5
0
5
10
Fig. 7. Least-squares exponential spline for Rentrop data (dots). The knots are marked by squares. For this example the larger range [1, 40] for scaled tension parameters gives a slightly better result (bottom right) than for the range [2, 20] (third row, right). All plots apart from the bottom right restrict the scaled tension parameters to [2, 20].
polynomial splines, which have similar characteristics as exponential splines. Considering the Jüttler data, they get a graphically comparable least-squares approximation—the least-squared error is not mentioned—with 15 intervals. This means 14 inner knots and 15 variable polynomial degrees for the 15 intervals, i.e., totally 29 degrees of freedom.
108
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114 σ2= 0.234; equadr= 1.88; emax= 0.559
σ1= 0.285; equadr= 2.57; emax= 0.694 5
5
4
4
3
3
2
2
1
1
0
0 -10
-5
0
5
10
-10
0
5
10
σ4= 0.0814; equadr= 0.488; emax= 0.352
σ3= 0.128; equadr= 0.899; emax= 0.42 5
5
4
4
3
3
2
2
1
1
0
-5
0 -10
-5
0
5
10
-10
σ5= 0.097; equadr= 0.485; emax= 0.35
0
5
10
σ6= 0.0678; equadr= 0.271; emax= 0.373
5
5
4
4
3
3
2
2
1
1
0
-5
0 -10
-5
0
5
10
-10
-5
0
5
10
σ7= 0.00564; equadr= 0.0169; emax= 0.0876 5
σ8= 0.00545; equadr= 0.0109; emax= 0.0636 5
4
4
3
3
2
2
1
1
0
0 -10
-5
0
5
10
-10
-5
0
5
10
Fig. 8. Least-squares cubic spline for the same data as in Fig. 7. The knots are marked by squares. Cubic splines show much more undesired oscillations than exponential splines—compare Fig. 7.
In Fig. 10, we use Algorithm 1 in order to find a least-squares exponential spline approximating the Jüttler data. Using eight interior knots and nine interval tension parameters, i.e., 17 degrees of freedom, we arrive at an approximation with a least-squared error equadr = 4.61—bottom right. Restricting the scaled tension parameter scal to 0.01, we use the same algorithm in order to find an approximating cubic spline. With eight interior knots, i.e., a total amount
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114 σ2= 0.209; equadr= 1.25; emax= 0.73
σ1= 0.426; equadr= 2.98; emax= 0.918 10
10
8
8
6
6
4
4
2
2
0
0
2
4
6
8
10
σ3= 0.0211; equadr= 0.106; emax= 0.253
0
10
8
8
6
6
4
4
2
2 0
2
4
6
8
10
σ5= 0.00103; equadr= 0.00308; emax= 0.0417
0
10
8
8
6
6
4
4
2
2 0
2
4
6
8
10
0
0
0
σ7= 7.53e-010; equadr= 7.53e-010; emax= 2.14e-005 10
8
8
6
6
4
4
2
2 0
2
4
6
8
10
4
6
8
10
2
4
6
8
10
2
4
6
8
10
equadr= 3.55e-030; emax= 1.26e-015
10
0
2
σ6= 2.18e-006; equadr= 4.36e-006; emax= 0.00318
10
0
0
σ4= 0.00199; equadr= 0.00794; emax= 0.0467
10
0
109
0
0
2
4
6
8
10
Fig. 9. Exponential least-squares spline for Späth data (dots). The knots are marked by squares. The local minimal distance condition prevents them from becoming closer than the data points. As many knots as data points lead to an interpolating exponential spline with suitable tension parameters.
110
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114 σ1= 18; equadr= 1080; emax= 7.28
σ2= 6.36; equadr= 375; emax= 4.86
10
10
5
5
0
0
-5
-5
-10 -5
10
0
5
10
15
σ3= 1.47; equadr= 85.1; emax= 2.17
-10 -5
10
5
5
0
0
-5
-5
-10 -5
0
5
10
15
-10 -5
10
5
5
0
0
-5
-5
0
5
10
15
-10 -5 10
5
5
0
0
-5
-5
0
5
10
10
15
σ4= 1.26; equadr= 71.5; emax= 2.26
0
5
10
15
0
5
10
15
σ8= 0.087; equadr= 4.61; emax= 0.667
σ7= 0.104; equadr= 5.6; emax= 0.683 10
-10 -5
5
σ6= 0.146; equadr= 8.02; emax= 0.709
σ5= 0.774; equadr= 43.4; emax= 2.05 10
-10 -5
0
15
-10 -5
0
5
10
15
Fig. 10. Least-squares exponential spline for Jüttler data (dots). The knots are marked by squares. Clearly, this example demands for the parametric setting.
of 8 degrees of freedom, we arrive at an approximating cubic spline with a least-squared error equadr = 6.43—see Fig. 11. Finally, Fig. 12 shows the curvature of the approximating exponential spline (left) with eight interior knots and the corresponding cubic spline (right). Negative curvature is dashed. Thus, there are exactly two inflection points for the exponential, as well as the cubic spline approximation like for the shape preserving approximation of Jüttler. The
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114 σ1= 18.5; equadr= 1111; emax= 7.24
σ2=10.9; equadr= 642; emax= 5.93
10
10
5
5
0
0
-5
-5
-10
-5
0
5
10
15
σ3= 2.05; equadr= 119; emax= 2.51
10
-10
5
0
0
-5
-5
-5
0
5
10
-5
15
-10
-5
10
5
5
0
0
-5
-5
-5
0
5
10
5
10
15
0
5
10
15
σ6= 0.17; equadr= 9.38; emax= 0.763
σ5= 1.32; equadr= 73.9; emax= 2.55 10
-10
0
σ4= 1.54; equadr= 88; emax= 2.39
10
5
-10
111
15
-10 -5
0
5
10
15
σ8= 0.121; equadr= 6.43; emax= 0.655
σ7= 0.123; equadr= 6.65; emax= 0.783 10
10
5
5
0
0
5
5
-10
-10 -5
0
5
10
15
-5
0
5
10
15
Fig. 11. Least-squares cubic spline for the same data as in Fig. 10. The knots are marked by squares. Due to the lack of sharp corners, the cubic spline is suitable.
approximating splines with less than eight interior knots shown in Figs. 10 and 11 also have not more than two inflection points. This means that at least for this example Algorithm 1 finds approximations without extraneous inflection points as desired. Note that the curvature of a parametric cubic spline (sx (t), sy (t)) need not be linear, actually it is in the 2D case sx (t)sy (t) − sy (t)sx (t). As expected, variations in the curvature are more pronounced for the exponential spline. Summarizing, we find approximations closer to the data which further demand for fewer degrees of freedom than the mentioned approaches which impose convexity constraints. Few degrees of freedom are preferable, because unnecessary
112
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114
Table 1 Performance of the minimization algorithm for three interior knots in dependence of the position-tension weighting factor c c
equadr
Iterations
100 50 10 1 0.1 0.01 0.001 0.0001 0.00001
0.00631891 0.00609905 0.00608029 0.00608025 0.00608025 0.00608024 0.00608028 0.00609962 No result
72 60 68 96 50 59 58 40 —
For the range c ∈ [0.001, 10] we get accurate minimization results. The number of necessary iterations is comparable. Outside this range we find inaccurate results or even worse numerical difficulties caused by very shallow surfaces used for minimization.
0.6 0.3 0.4 0.2 0.2 0.1 Curvature
Curvature
0
-0.2
-0.4
0
-0.1
-0.6 -0.2 -0.8 -0.3 -1 0
10
20
30
40 t
50
60
70
0
10
20
30
40
50
60
70
t
Fig. 12. (Left) Curvature of the exponential spline in Fig. 10 (bottom right). (Right) Curvature of the cubic spline in Fig. 11 (bottom right). Negative curvature is dashed. In both cases there are exactly two inflection points.
degrees of freedom are time consuming and can lead to undesired effects like spurious oscillations—unless avoided by convexity constraints. As demonstrated in Fig. 2 the knot positions are crucial. The curvature of the Jüttler data is moderate. Therefore, as demonstrated, cubic least-squares splines with locally optimal knots lead to an approximation close to the data with only 8 degrees of freedom. As undesired oscillations are related to overfitting, which is avoided by Algorithm 3.3, we can deal without convexity constraints.
5. Conclusion Least-squares splines in tension can be optimized by a simultaneous knot placing and tension setting algorithm. The dynamical equilibrium between the alteration rate of the knot positions and the tension parameters makes it easy to find
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114
113
a suitable relative alteration rate of the tension parameters compared to the knot positions. Scaled tension parameters are useful in order to keep the minimization away from shallow regions. Using a local minimal distance condition avoids overfitting and leads to interpolating splines with coincident knots and data points as a special case. This includes finding suitable tension parameters of interpolating exponential splines. The algorithm presented here usually does not violate natural convexity requirements inherited by the data. If such a case should occur, enlarging the tension of the critical interval(s) is a remedy. The simultaneous knot placing and tension setting approach presented here in combination with variable degree polynomial splines [4] is possible if the minimizer can deal with discrete polynomial degrees instead of smooth tension parameters. The generalization of our approach to 3D curves is straightforward.
Acknowledgements I wish to thank my academic advisor Prof. Dr. R. Seydel for proposing and supporting this work on exponential splines. Further, my thanks go to Dr. K.-D. Reinsch for motivating and profound discussions. Prof. Dr. B. Jüttler kindly provided the data used in Figs. 10 and 11. Last but not least, I wish to thank the referees for their supplements and improvements. References [1] B.A. Barsky, Exponential and polynomial methods for applying tension to an interpolating spline curve, Computer Vision, Graphics, Image Process. 27 (1984) 1–18. [2] C. de Boor, A Practical Guide to Splines, Springer, Berlin, Heidelberg, New York, 1978. [3] P. Costantini, F. Pelosi, Shape-preserving approximation by space curves, Numer. Algorithms 27 (2001) 237–264. [4] P. Costantini, F. Pelosi, Shape-preserving approximation of spatial data, Adv. Comput. Math. 20 (2004) 25–51. [5] R. Delbourgo, J.A. Gregory, C 2 rational quadratic spline interpolation to monotonic data, IMA J. Numer. Anal. 3 (1983) 141–152. [6] R. Delbourgo, J.A. Gregory, The determination of the derivative parameters for a monotonic rational quadratic interpolant, IMA J. Numer. Anal. 3 (1983) 141–152. [7] P. Dierckx, Curve and Surface Fitting with Splines, Oxford Science Publications, Clarendon Press, Oxford, 1993. [8] F.N. Fritsch, R.E. Carlson, Monotone piecewise cubic interpolation, SIAM, J. Numer. Anal. 17 (1980) 238–246. [9] J.A. Gregory, R. Delbourgo, Piecewise rational quadratic interpolation to monotonic data, IMA J. Numer. Anal. 2 (1982) 123–130. [10] A. Gregory, M. Sarfraz, A rational cubic spline with tension, Comput. Aided Geom. Design 7 (1990) 1–13. [11] J.G. Hayes, Curve fitting by polynomials in one variable, in: J.G. Hayes (Ed.), Numerical Approximation to Functions and Data, The Athlone Press, London, 1970, pp. 43–64. [12] U. Heidemann, Linearer Ausgleich mit Exponentialsplines bei automatischer Bestimmung der Intervallteilungspunkte, Computing 36 (1986) 217–227. [13] W. Heß, J.W. Schmidt, Convexity preserving interpolation with exponential splines, Computing 36 (1986) 335–342. [14] D.J. Higham, Monotonic piecewise cubic interpolation, with applications to ODE plotting, TR 229/90, Department of Computer Science, University of Toronto, 1990. [15] B. Jüttler, Shape preserving least-squares approximation by polynomial parametric spline curves, Comput. Aided Geom. Design 14 (1997) 731–747. [16] P.E. Koch, T. Lyche, Exponential B-splines in tension, in: C.K. Chui, L.L. Schumaker, J.D. Ward (Eds.), Approximation Theory VI, vol. 2, ISBN 0-12-174587-2, Academic Press, London, 1989, pp. 361–364. [17] P.E. Koch, T. Lyche, Construction of exponential tension B-splines of arbitrary order, in: P.J. Laurent, A. Le Méhauté, L.L. Schumaker (Eds.), Curves and Surfaces, ISBN 0-12-438660-1, Academic Press, London, 1991, pp. 255–258. [18] P.E. Koch, T. Lyche, Interpolation with exponential B-splines in tension, Computing Supplementum 8, in: G. Farin, H. Hagen, H. Noltemeier (Eds.), Geometric Modelling, ISBN 3-211-82399-9, Springer, Wien, New York, 1993, pp. 173–190. [19] R.W. Lynch, A method for choosing a tension factor for spline under tension interpolation, M.S. Thesis, University of Texas at Austin, 1982. [20] G.M. Nielson, Computation of -splines, TR 044-433-11, Department of Mathematics, Arizona State University, 1974. [21] G.M. Nielson, Some piecewise polynomial alternatives to splines under tension, in: R.E. Barnhill, R.F. Riesenfeld (Eds.), Computer Aided Geometric Design, 1974, pp. 209–235. [22] S. Pruess, Properties of splines in tension, J. Approx. Theory 17 (1976) 86–96. [23] R.J. Renka, Interpolatory tension splines with automatic selection of tension factors, SIAM J. Sci. Statist. Comput. 8 (3) (1987) 393–415. [24] R.J. Renka, Tspack: tension spline curve-fitting package, ACM Trans. Math. Software 19 (1993) 81–94. [25] P. Rentrop, An algorithm for the computation of the exponential spline, Numer. Math. 35 (1980) 81–93. [26] K.O. Riedel, Aspects of image processing—splines, anisotropic diffusion, and biological models, Ph.D. Thesis, ISBN 3-8322-1236-1, Shaker Verlag, Aachen, 2003. [27] K.O. Riedel, Two-dimensional splines on fairly arbitrary meshes, ZAMM—Z. Angew. Math. Mech. 85 (3) (2005) 176–188.
114
K.O. Riedel / Journal of Computational and Applied Mathematics 196 (2006) 94 – 114
[28] N.S. Sapidis, P.D. Kaklis, T.A. Loukakis, A method for computing the tension parameters in convexity-preserving spline-in-tension interpolation, Numer. Math. 54 (1988) 179–192. [29] M. Sarfraz, A geometrical rational spline with tension controls: an alternative to the weighted nu-spline, Punjab Univ. J. Math. 26 (1993) 27–40. [30] M. Sarfraz, Freeform rational bicubic spline surfaces with tension control, Facta Univ. Ser. Math. Inform. 9 (1994) 83–93. [31] D.G. Schweikert, An interpolating curve using a spline in tension, J. Math. Phys. 45 (1966) 312–317. [32] H. Späth, Exponential spline interpolation, Computing 4 (1968) 225–233. [33] J. Stoer, Numerische Mathematik 1, Springer, Berlin, Heidelberg, New York, 1999. [34] J. Stoer, R. Bulirsch, Introduction to Numerical Analysis, Springer, Berlin, Heidelberg, New York, 1980.