Many data mining tasks can be reformulated as optimization problems, in the solution of which approximation by surfaces plays a key role. The paper proposes a new efficient computational algorithm for spline surfaces over uniform grids. The algorithm is based on a recent result on approximation of a biquartic polynomial by bicubic ones, that ensures C2 continuity of the corresponding four bicubic spline components. As a consequence of this biquartic polynomial based approach to constructing spline surfaces, the classical de Boor's computational task breaks down to a reduced task and a simple remainder one. The comparison of the proposed and classical computational algorithm shows that the former needs less multiplication operations resulting in non negligible speed up.
Recent years a considerable effort has been seen to
develop reliable and efficient data mining tools to discover
hidden knowledge in very large data bases. The
fundamental problem is proposing algorithms to extract some
useful information from very large databases. Fortunately,
many data mining tasks can be reformulated as
optimization problems, where approximation by surfaces plays
a key role [
The goal of the paper is to show that even such standard result as de Boor’s sequential algorithm for construction of interpolating spline surfaces can be improved. It suggests a computational algorithm based on new model equations, in derivation of which biquartic polynomials played an essential role.
The idea of using quartic and biquartic polynomials in
cubic and bicubic spline construction comes from recent
results of Török and Szabó [
The algorithm presented in this paper has a decreased
number of equations and is based on the generalized
results of [
The structure of the article is as follows. Section two is devoted to problem statement. To be self-contained, section three briefly describes de Boor’s algorithm. The next section first provides the definition of bicubic and biquartic polynomials. Then it shortly discusses the interrelation of these bivariate polynomials and their role in the computational schema. Section five contains the proposed sequential computational algorithm based on reduced systems. Section six briefly compares the new and the classical algorithm. The efficiency of the proposed algorithm is shown in the last but one section by computing the theoretical speedup that is approximately 1.33. 2
The section defines the inputs for the spline surface, and the requirements that it should fulfil and based on which it can be constructed.
Consider a uniform grid [u0, u1, . . . , u2m] × [v0, v1, . . . , v2n], (1) where ui = u0 + ihx, i = 1, 2, . . . , 2m, m ∈ N, v j = v0 + jhy,
j = 1, 2, . . . , 2n, n ∈ N.
According to [
j = 0, 1, . . . , 2n y di, j, i = 0, 1, . . . , 2m, j = 0, 2n at boundary horizontals and cross derivatives dix,,jy, i = 0, 2m, j = 0, 2n at the four corners of the grid.
The task is to define a quadruple [zi, j, dix, j, diy, j, dix,,jy] at every grid-point [ui, v j], based on which a uniform bicubic clamped spline surface S of class C2 can be constructed with properties (3) (4) (5)
S(ui, v j) = zi, j, ∂ S(ui, v j) ∂ x = dix, j, ∂ S(ui, v j)
∂ y ∂ 2S(ui, v j) ∂ x∂ y
y
= di, j,
= dix,,jy,
where the adjacent spline segments are twice continuously
differentiable. Our aim is to solve this task with less
equations and less multiplications than the standard spline
construction algorithm [
Paper [
Lemma 1 (de Boor). If the above z values and d derivatives are given, then the values dix, j,
y di, j, dix,,jy, i = 1, . . . , 2m − 1, i = 0, . . . , 2m,
j = 0, . . . , 2n, j = 1, . . . , 2n − 1, i = 1, . . . , 2m − 1, and i = 0, . . . , 2m, j = 1, . . . , 2n − 1 are uniquely determined by the following 2(2m) + (2n) + 5 linear systems of altogether 3(2m)(2n) + (2m) + (2n) − 5 equations: for j = 0, . . . , 2n, di+1, j + 4dix, j + di−1, j = x x (zi+1, j − zi−1, j), dix+,y1, j + 4dix,,jy + dix−,y1, j = (diy+1, j − di−1, j), y where i = 1, . . . , 2m − 1; for i = 0, . . . , 2m, where j = 1, . . . , 2n − 1; for i = 0, . . . , 2m, where j = 1, . . . , 2n − 1.
di, j+1 + 4diy, j + di, j−1 = hy y y
(zi, j+1 − zi, j−1), dix,,jy+1 + 4dix,,jy + dix,,jy−1 = (dix, j+1 − di, j−1),
x 3 hx 3 3 hy 4
The section begins with definition of bicubic and biquartic polynomials. Then it shortly discusses the interrelation of these bivariate polynomials and its role in the computational schema.
The tensor product formulas of bicubic Hermite spline
components, see [
Definition 1. On grid (1) the bicubic Hermite spline components Si, j(x, y) for i = 0, 1, 2, . . . , 2m − 1,
j = 0, 1, 2, . . . , 2n − 1, x ∈ [ui, ui+1], y ∈ [v j, v j+1], are defined as follows
Si, j(x, y) = λ T (x, ui, hx) · ϕ i, j · λ (y, v j, hy), (10) where λ is a vector of basis functions (1+2 t−ht0 )(t−t1)2 T
h2 (t−t0)2(1−2 t−ht1 ) λ (t, t0, h) = (t−t0)h(2t−t1)2 ,
h2 (t−t0)2(t−t1) h2 t1 = t0 + h and ϕ is a matrix of function values and first derivatives
zi, j ϕ i, j = zi+1, j
x di, j x di+1, j For the spline components the following conditions hold For uk, vl defined in (1) the following conditions hold Si, j(uk, vl ) ∂ Si, j(uk, vl )
∂ x ∂ Si, j(uk, vl )
∂ y ∂ 2Si, j(uk, vl ) ∂ x∂ y = = = = zk,l , dkx,l ,
y dk,l , dkx,,ly, k = i, i + 1, l = j, j + 1, k = i, i + 1, l = j, j + 1, k = i, i + 1, l = j, j + 1, k = i, i + 1, l = j, j + 1.
Based on (10) the second derivatives of Si, j(x, y) can be expressed effectively, e.g.
∂ 2Si, j(x, y) ∂ x2 = ∂ 2λ T (x, ui, hx) ∂ x2 · ϕ i, j · λ (y, v j, hy), (11) where ∂ 2λ (t, t0, h) ∂ t2 6(2t−2t0−h) T
h3 6(−2t+2t0+h) = 2(3t−3t0−2h) .
h3 h2 2(3t−3t0−h) h2
The biquartic polynomials are also defined by tensor product.
On grid (1) the biquartic polynomials i = 0, 2, 4, . . . , 2(m − 1), j = 0, 2, 4, . . . , 2(n − 1), x ∈ [ui, ui+2], y ∈ [v j, v j+2],
Fi, j(x, y) for
Fi, j(x, y) = LT(x, ui, hx) · Φi, j · L(y, v j, hy), (12) where L is a vector of basis functions −(1+2 t−ht0 )(t−t1)(t−t2)2 T
4h3 (t−t0)h2(4t−t2)2 L(t, t0, h) = (t−t0)2(t−t1)(1−2 t−ht2 ) , −(t(−t−t0t0))2((t44t−−hh33tt11))((tt−−tt22))2 4h3 t1 = t0 + h, t2 = t0 + 2h and Φ is a matrix of function values and first derivatives zi, j zi+1, j Φi, j = zi+2, j x di, j x di+2, j
The tensor product definition of Fi, j(x, y) by (12) provides a compact way to express first derivatives, e.g. ∂ Fi, j(x, y) ∂ y = LT(x, ui, hx) · Φi, j · ∂ L(y, v j, hy) .
∂ y (13)
Works [
Unlike de Boor’s lemma, we provide only the
interpretation of the main result of [
This interrelation between a biquartic and four bicubic polynomials is illustrated by the schema in Fig. 1. The biquartic polynomial F over [u0, u2] × [v0, v2] is defined by given nine function values z and sixteen derivatives d that set up four quadruples [z, dx, dy, dx,y], two pairs [z, dx], two pairs [z, dy] and a single z. Every bicubic spline component is defined by four quadruples[z, dx, dy, dx,y]. The nine quadruples in the figure are depicted around nine gridpoints. Those eleven directional and cross first derivatives that are computed from the biquartic polynomial F and that are needed to construct the four bicubic spline components S = {S0,0, S0,1, S1,0, S1,1}, are denoted by δ , see Fig. 1.
The below proposed algorithm was developed by generalizing the above described interrelation between biquartic and bicubic polynomials. First biquartic polynomials were handled and then based on them new model equations and formulas for unknown derivatives of the bicubic spline surface were derived.
Figure 2 illustrates the schema of the proposed computational algorithm for (2m + 1) × (2n + 1) bicubic spline surface of class C2. There are 2m + 1 verticals and 2n + 1 horizontals. Rectangles and thick rectangles indicate the boundary of bicubic spline components and biquartic polynomials, respectively. There are two types of objects at every grid-point: known and unknown ones. The given values and derivatives are denoted by z, d and the unknown first derivatives byD, δ . Notice that z is provided at every grid-point and d only along the total grid’s boundary. The most important is where are the unknown D and δ parameters. The D parameters are located only along the thick rectangles, but never in their center. As we shall see later the D parameters will be computed from equations and the δ parameters from explicit formulas. The equations were derived from the equality of second derivatives of spline components and the formulas from the biquartic polynomials.
The derived new model equations for the unknown D
derivatives of the spline surface segments and parts of the
explicite formulas for δ are generalization of model
equations and formulas of the unknown derivatives of spline
curve segments, see [
The below proposed algorithm can be characterized from two aspects • what it computes, • the quality of its outcome.
The algorithm computes D and δ coefficients for bicubic spline surface components from inputs given at equispaced grid-points described in Section 2. The D coefficients are computed from linear systems based on equations (14), (16), (18) – (21). The δ coefficients are gained from explicit formulas (15), (17), (22) – (24). Since the equations for the D parameters were derived from the equality of second derivatives of spline components and the formulas for the δ parameters were gained from biquartic polynomials that as we know grant C2 continuity of their components, the algorithm provides such coefficients that the uniform bicubic spline surface will be of class C2.
Algorithm for computing the unknown first derivatives of the spline surface in three main steps with reduced systems.
Inputs: z and d values, see (2) – (5).
izontals from equation systems.
For each horizontal we construct a system of linear equations to compute the Dx values, located on the inside odd grid-points. Each horizontal represents an independent tridiagonal system of linear equations.
For each horizontal, see Fig. 2, j = 0, 1, . . . , 2n, a tridiagonal system is constructed based on equations D2(i+1), j − 14Dx2i, j + D2(i−1), j = x x
3 12 = hx (z2(i+1), j − z2(i−1), j) − hx (z2i+1, j − z2i−1, j), (14) where i = 1, 2, . . . , m − 1.
formulas.
To finish the computation of all first partial derivatives with respect to x we have to calculate δix, j = 3 1 4hx (zi+1, j − zi−1, j) − 4 (dix+1, j + dix−1, j), (15) where i = 1, 3, . . . , 2m − 1, j = 1, 3, . . . , 2n − 1 .
izontals from equation systems.
For each vertical we construct a system of linear equations to compute the Dy values, located on the inside odd grid-points. Each vertical represents an independent system of linear equations.
For each vertical, i = 0, 1, . . . , 2m, a tridiagonal system is constructed based on equations
Di,2( j+1) − 14Diy,2 j + Di,2( j−1) =
y y where j = 1, 2, . . . , n − 1.
formulas
To finish the computation of all first partial derivatives with respect to y we have to calculate δiy, j = 3 1 y 4hy (zi, j+1 − zi, j−1) − 4 (diy, j+1 + di, j−1), where i = 1, 3, . . . , 2m − 1, j = 1, 3, . . . , 2n − 1.
At this moment all directional derivatives are known: some were provided and the unknown D and δ directional ones were computed in Steps 1 and 2. In the further steps all directional derivatives will be denoted by d and contained in the right hand side of equations and formulas.
tom and top horizontals, and left vertical from equation systems.
We construct systems of linear equations for bottom and
top horizontals and left verticals by [
side grid-points using systems of equations
For the odd verticals, i = 2, 4, 6, . . . , 2m, a tridiagonal system is constructed based on equations
Dix,,jy+2 − 14Dix,,jy + Dx,y
i, j−2 = 1 = 7 (dix−,y2, j+2 + dix−,y2, j−2) − 2dix−,y2, j+ + 73hx (diy−2, j+2 + diy−2, j−2) + 73hy (−dix−2, j+2 + dix−2, j−2)+ 9 9 + 7hx (diy, j+2 + diy, j−2) + 7hxhy (−zi−2, j+2 + zi−2, j−2)+ + 71h2x (−diy−1, j+2 − diy−1, j−2) + 71h2y (dix−2, j+1 − dix−2, j−1)+ + + h3y (dix, j+2 − dix, j−2) + 7h2x7hy (−zi, j+2 + zi, j−2)+
36 7hxhy (zi−1, j+2 − zi−1, j−2 + zi−2, j+1 − zi−2, j−1)− − h6x diy−2, j + 12 hy − 1h8x diy, j + 71h4x4hy (−zi−1, j+1 + zi−1, j−1) + 2h4x diy−1, j, 108 (dix, j+1 + dix, j−1) + 7hxhy (zi, j+1 − zi, j−1)− (21) where j = 4, 6, . . . , 2n − 4.
Mention must be made, that this step was the most critical. At first after computation ofDx,y unknowns along the bottom and top horizontals in Step 3a we got equations with six Dx,y unknowns on the left side. Török suggested to compute the Dxy parameters along the left vertical using de Boors equation in Step 3a and thanks to this three of six Dx,y parameters could be moved to the right side as computed.
formulas
To finish the computation of all first cross derivatives we have to calculate for the even verticals and the even horizontals
1 δix,,jy = 16 (dix+,y1, j+1 + dix+,y1, j−1 + dix−,y1, j+1 + dix−,y1, j−1)− 3 x x − 16hy (dix+1, j+1 − di+1, j−1 + dix−1, j+1 − di−1, j−1)− 3 y y y − 16hx (diy+1, j+1 + di+1, j−1 − di−1, j+1 − di−1, j−1)+ 9
(zi+1, j+1 − zi+1, j−1 − zi−1, j+1zi−1, j−1), +
16hxhy where i = 1, 3, . . . , 2m − 1, j = 1, 3, . . . , 2n − 1; for the even verticals and the odd horizontals
δix,,jy = 43hy (dix, j+1 − dix, j−1) − 41 (dix,,jy+1 + dix,,jy−1), where i = 1, 3, . . . , 2m − 1, j = 2, 4, . . . , 2(n − 1); and for the odd verticals and the even horizontals δix,,jy = 43hy (dix, j+1 − dix, j−1) − 41 (dix,,jy+1 + dix,,jy−1), where i = 2, 4, . . . , 2m − 1, j = 1, 3, . . . , 2n − 1. (22) (23) (24) 6
We sum up the complete spline task and how is it completed by the considered two algorithms. Then we give some details about the role of the biquartic polynomials that are absent in de Boor’s algorithm but played a crucial task in the design of the proposed one.
In case of (2m + 1)(2n + 1) grid-points, to fulfill the complete spline task means to construct (2m)(2n) spline components using (2m + 1)(2n + 1) various quadruples, where one quadruple looks the following way: [z, dx, dy, dx,y]. The input comprises (2m + 1)(2n + 1) z values, 2(2n + 1) dx, 2(2m + 1) dy and 4 dx,y derivatives. The unknown derivatives require computation.
In de Boor’s algorithm every unknown derivative is
computed from equations that are of four types, see (6)–(9)
and [
(2n + 1)(2m − 1) derivatives dx are computed, • from 2 systems with 2m − 1 variables
2(2m − 1) derivatives dx,y are computed, • from 2m + 1 systems with 2n − 1 variables
(2m + 1)(2n − 1) derivatives dy are computed, • from 2m + 1 systems with 2n − 1 variables
(2m + 1)(2n − 1) derivatives dx,y are computed.
The proposed algorithm’s benefit is that lesser number
of unknown derivatives (D parameters) are computed from
systems of equations, see Tab. 3, compared to de Boor’s
algorithm [
The algorithm’s drawback is that it uses more types of relation: six types of equations and five types of explicit formulas. Nevertheless it has to compute approximately 12/5 times less equations within systems – compare Tab. 2 and Tab. 3 from Section 7.
Let us have a closer look at biquartic polynomials and their role in the algorithm’s design. One biquartic polynomial Fi, j(x, y) needs 25 parameters, see Definition 2. We distinguish between four types of F polynomials, see Fig. 2,
they are D parameters, for F from the corners only seven
are unknown and nine are given, these are the d
parameters. After obtaining the D parameters the remaining
derivatives are computed based on F . From every F
eleven delta parameters can be obtained: two pairs of
type [δ x, δ x,y], two pairs of type [δ y, δ x,y] and one triple
[δ x, δ y, δ x,y]. Naturally, δ parameters are functions of D
parameters, see [
After introducing the new algorithm in the previous section and giving a short insight into its design in this one, the next section is devoted to its quantitative characterization. 7
The standard way of solving a tridiagonal linear system b 1 0 | where dK = uyKK , {di = u1i (yi − xi+1)}, i = K − 1, . . . , 1.
The tridiagonal systems of equations for de Boor’s and our algorithm are solved by LU decomposition. All the systems of these algorithms are diagonally dominant with
Dim.
N × N (25) LU γ N (26) Fw N (27) Bw γ N
RHS β N
2γ N + β N + N elements 1, 4, 1 and 1, −14, 1. The LU and backward steps contain a division that is indicated by γ , the ratio between division and multiplication: the performance of a division operation is equivalent to γ multiplications. The proposed and de Boor’s algorithm differ • in number of systems of equations, • in number of equations within systems, • in number of multiplication operations on right hand side (RHS) of equations.
Tab. 1 presents the number of multiplications for solving one general tridiagonal N × N matrix, where β denotes the number of multiplications on the right hand side of an equation.
The second and third columns of Tab. 2 and Tab. 3 provide the count of equations within the given steps (equations) and the count of equations within a system, respectively, for a grid of size (2m + 1) × (2n + 1). The last but one rows contain the total number of equations. The total count of multiplications to solve the tridiagonal systems within de Boor’s and the proposed algorithm are in the last row.
In the proposed algorithm we evaluate in addition to the solution of the tridiagonal systems of equations, see Tab. 3, as well as δ parameters using explicit formulas. Therefore the total number of multiplications in the proposed algorithm based on tables 3 and 4 is 10γ mn + 46mn + 4γ m + 2γ n − 13m + 3n − 10γ − 12.
We can see that the number of multiplication in the proposed algorithm is less. The theoretical speed up for some various grid sizes were computed under assumption that γ = 3.5. Based on Tab. 5 we can conclude that the proposed model is asymptotically 1.33 times faster. We suggested a new efficient sequential algorithm for computation of a spline surface over an equispaced grid. Its theoretically evaluated asymptotic speed up over de Boors algorithm is approximately 1.33. The algorithm has also a very nice property from the view point of parallel computation: the computation of the second half of unknowns based on explicit equations can be parallelized automatically. Naturally, parallel methods, suggested for solving tridiagonal systems of de Boor’s algorithm can be used for solution of the new algorithm’s tridiagonal systems as well. Therefore, the proposed reduced system based algorithm is preferable over de Boor’s algorithm not only for sequential, but also for parallel computation.
The author thanks Cs. Török for the problem statement, his helpful suggestions and support. This work was partially supported by the research grants VEGA 1/0073/15 and VVGS-PF-2015-477.