GENERALIZATION OF ONE ALGORITHM FOR CONSTRUCTING RECURRENT SPLINES

The task on finding an effective analytical description of a time series is commonly in many application areas. The scope of application of interpolation splines is known to have been limited to processing experimental dependences that do not contain measurement errors or contain errors that can be ignored. In the opposite case, the smoothing splines are employed. There are many methods for smoothing. That relates to the fact that the smoothing conditions are chosen to match a particular applied problem. An expert in the subject often struggles to understand the peculiarities of different smoothing methods and to assign the required parameters of smoothing. Therefore, there is a need to develop algorithms of smoothing whose constraints are interpreted so that a user who may lack a proper mathematical training on this subject can understand them. The task on constructing an optimal smoothing spline implies sorting out all the possible cases of distribution of grid knots among the links of the spline. It is clear that the resource intensity of algorithms that could fully solve such a problem makes them unsuitable for the application in real-time systems [1]. Therefore, it is a relevant task to build algorithms for constructing the smoothing splines, close to optimal, provided there is the possibility to modify a small number of spline links in proportion to the arrival of new data. Such algorithms would make it possible to achieve the approximation of an experimental dependence with acceptable quality for much less time than when solving a given problem in classical statement.


Introduction
The task on finding an effective analytical description of a time series is commonly in many application areas.The scope of application of interpolation splines is known to have been limited to processing experimental dependences that do not contain measurement errors or contain errors that can be ignored.In the opposite case, the smoothing splines are employed.There are many methods for smoothing.That relates to the fact that the smoothing conditions are chosen to match a particular applied problem.An expert in the subject often struggles to understand the peculiarities of different smoothing methods and to assign the required parameters of smoothing.Therefore, there is a need to develop algorithms of smoothing whose constraints are interpreted so that a user who may lack a proper mathematical training on this subject can understand them.
The task on constructing an optimal smoothing spline implies sorting out all the possible cases of distribution of grid knots among the links of the spline.It is clear that the resource intensity of algorithms that could fully solve such a problem makes them unsuitable for the application in real-time systems [1].Therefore, it is a relevant task to build algorithms for constructing the smoothing splines, close to optimal, provided there is the possibility to modify a small number of spline links in proportion to the arrival of new data.
Such algorithms would make it possible to achieve the approximation of an experimental dependence with acceptable quality for much less time than when solving a given problem in classical statement.

Literature review and problem statement
The tasks on constructing smoothing splines in various statements have traditionally attracted attention of experts in computational methods and specialists in different subject areas that employ the apparatus of splines [2,3].However, these studies mainly concern the development of algorithms
One of the important components of the task on constructing smoothing splines is the search for knots, in which spline links are joined.A solution to this auxiliary problem using a genetic algorithm is proposed in paper [4].The same problem for B-splines is solved in article [5] using the author's two-stage algorithm taking into consideration a curvature of the experimental dependence.
In paper [6], in order to establish coordinates of the knots where spline links are joined, it is proposed first to form their redundant sequence followed by chopping using optimization methods.The framework of this approach includes the results obtained by authors of article [7].
The issue of constructing a robust smoothing B-spline is examined in paper [8].The above list shows that most of the results were obtained for B-splines and rely on methods that are time-expensive.
We shall consider two algorithms for constructing smoothing splines using which the problems on robustness and splitting a spline into links, close to optimal, can be solved simultaneously.These algorithms have the common feature that implies a partial imposing of the spline's adjacent links one on another.We note that the traditional approach to the construction of splines implies that adjacent links have one common point.Both algorithms are focused on the post-construction of a new link of the spline in proportion to the arrival of new experimental data without changing all the preceding links of a given spline during approximation of a time series.
The first algorithm is described, in particular, in paper [9].The algorithm allows the use of non-uniform grids.To describe a spline link, a polynomial with a specific structure is used, the properties of which are examined below.The algorithm is based on the synthesis of ideas about Hermitian interpolation and mean-square approximation.The result of its execution is the smoothing spline of order С 0 .
The second algorithm is presented, for example, in paper [10].The existence of a semi-local spline and its convergence were proven for uniform grids [11].To describe a link of the spline, the authors used a polynomial with a representation in a power base.The result of the algorithm execution makes it possible to obtain a spline with a preset order of smoothness [12].
Each of the algorithms has both undeniable advantages and certain disadvantages.Thus, a shortcoming of the first algorithm is the low degree of smoothness of the spline obtained.In addition, local estimation of part of the coefficients of a polynomial from the current spline's link for experimental dependences with significant measurement errors may lead to a significant worsening in the approximating properties of the spline in general.
Effectiveness of applying the second algorithm is limited by the requirement to construct links of the spline strictly equal in length.For the plots of experimental sequences that demonstrate the slow nature of changes, it is advisable to apply longer links of the spline.This, in turn, makes it possible to bring down the cost of resources when maintaining a description of the resulting spline and during its application.Therefore, it is expedient to synthesize these algorithms in order to improve their advantages and to eliminate disadvantages.

The aim and objectives of the study
The aim of present study is to examine a possibility to generalize the D. A. Silaev algorithm for irregular grids and to strengthen its computational stability.This would make it possible to reduce resource intensity of the examined algorithm when used in monitoring systems.
To accomplish the aim, the following tasks have been set: -to substantiate the feasibility of transition from representing a polynomial to describe links of the spline in the power base to other types of bases, specifically the bases of Hermite polynomials, the base of N. D. Dicoussar polynomial; -to show a relation between the N. D. Dicoussar polynomial and other known types of polynomials; -to examine the impact of the form of polynomial representation to describe links of the spline on the magnitude of the conditionality numbers of stability matrix and the Gram matrix, which are used in the D. A. Silaev algorithm; -to explore the effectiveness of constraints that are applied to determine the length of the spline's current link in the N. D. Dicoussar algorithm, in order to use them when generalizing on the irregular grids of the D. A. Silaev algorithm; -using test examples, to explore the resistance of the D. A. Silaev algorithm, and modifications of this algorithm, against the errors in experimental data.

1. Relation between the N. D. Dicoussar polynomial and other kinds of polynomials
In cases when the study relates to one link of the spline in order to reduce the number of indexes we shall equate the region for determining this link to section [ ] ; .
a b x x Paper [9] proposed to describe a link of the cubic smoothing spline by a polynomial of the following form: where r=(f a ; f b ; f 0 ) is the vector of weight coefficients at basis elements; d=(d 1 ; d 2 ; d 3 ) is the vector of basis elements in the form of polynomials of the second degree: is the auxiliary "zero" cubic parabola.
The basis elements (in the terminology of paper [9]) (d 1 ; d 2 ; d 3 ) are the basis functions in the traditional sense, associated with knots x a , x 0 and x b .That is, these functions are equal to 1 in their knot and are equal to 0 in all other knots.Function Q is zero in all three knots, as derived from its formula.
We shall show in advance that the set of functions (d 1 ; d 2 ; d 3 ; Q) forms the basis of space of polynomials with degree not higher than the third degree [13].To this end, we shall construct a matrix from coefficients of the specified functions: and compute the determinant: It is obvious that functions { d 1 ; d 2 ; d 3 ; Q Q} are linearly independent and form a basis of the space of polynomials with a degree not higher than the third degree.
Find the matrix that is inverse to matrix D and denote it V D : By using the technology of block matrices for constructing basis functions, which is described in paper [14] (in relation to the construction of basis functions in a finite element method), one can show that polynomial (1) in the matrix form can be recorded as: Author of paper [9] obtianed an estimate of the magnitude of parameter q through the value of function f(x) and its derivatives at the ends of segment [x a ; x b ].We show that it can be achieved in another way.
Let us consider the Hermitian polynomial, which hereafter is referred to as the Hermite polynomial of the first kind, defined by their values and the values of the first derivatives at the ends of current segment [x a ; x b ].Paper [14] shows that it can be represented in the following matrix form: where Because formulae (3) and ( 4) describe the same polynomial in different bases, we can equate the right sides of these formulae: Hence, we obtain With respect to the structure of matrices V and V D , we obtain: Hence, we obtain an expression for parameter q: ( ) which fully coincides with the expression, derived in the N. D. Dicoussar papers.
We note that the value of this parameter does not depend on the localization of knot x 0 .
It will be shown below that it is appropriate to choose the origin of a local coordinate system for the spline's link at point x 0 , that is, in this case, x 0 .=0.The matrix of the N. D. Dicoussar polynomial coefficients (2) then takes the form: This allows us to consider the N. D. Dicoussar polynomial as a hierarchical form of the Hermitian polynomial.

2. Approximating properties of the N. D. Dicoussar polynomial
When executing both of the examined algorithms, a system of linear algebraic equations is solved, which is traditional for the least squares method.It is known the matrix of this system coincides with the Gram matrix for basis functions used to represent the approximating ) polynomial.Therefore, an important characteristic of the basis that is chosen to represent a polynomial is the conditionality number of its Gram matrix.Not reducing the generality of results, we shall compute elements of the Gram matrix for the basis functions of polynomial (1) for segment [x a ; x b ]=[-1;1]: where i, jÎ{1, 2, 3, 4} considering that d 4 =Q.
We shall tabulate values for a conditionality number of the the Gram matrix from the position of point 0 .
x Paper [9] proposes fixing the positions of point x 0 as second knot of approximation at segment [x a ; x b ].Our calculations show that this recommendation may be disregared.The dependence of the Gram matrix conditionality number on the position of point x 0 (Fig. 1) shows that a shift of the right end of segment x b leads to a rapid growth in the Gram matrix conditionality number (along with the accumulation of calculation errors).That is why knot x 0 must be chosen in a spline grid knot, which is the closest to the midpoint of segment [x a ; x b ].The results obtained (Table 1) correspond to the point of view on this issue, known from literature, for example, from ref. [15].Based on data from Table 1, we conclude that the application of the N. D. Dicoussar interpolation polynomial (1) in approximate calculations demonstrates considerable advantages in comparison with the Hermitian interpolation polynomial because its basis functions possess better approximating properties.

3. Generalization of the D. A. Silaev algorithm
We shall introduce changes to the algorithm described in paper [10].This relates to the quantity of points, used in the method of least squares, to find coefficients at senior degrees of the polynomial.We also propose changes to the form of a polynomial representation, which describes the current link of the spline.
We shall consistently add knots to the grid at which the problem on the minimization of a functional of least squares is solved until the resulting polynomial deviates from experimental data by less than the specified number d>0.The final number of knots from the current interval at which this condition is fulfilled will be designated (M+1).
The number of points that finally remains in the composition of the current link of the spline will be designated (m+1).
To assess the impact of the form of representation of a polynomial on the stability of the modified algorithm, we shall use such kinds of polynomials.
When building a spline of smoothness С 0 , we shall apply: -polynomials in a power basis; -the N. D. Dicoussar polynomials of form (1).
When building a spline of smoothness С1, we shall use: -polynomials in a power basis; -Hermitian polynomials of the first kind with two knots in which values for the function and its first derivatives are assigned; -Hermitian polynomials of the second kind with three knots in which values for the function are assigned, an in the first knotvalues for the function's derivative.
When building a spline of smoothness С 2 , we shall use: -polynomials in a power basis; -Hermitian polynomials of the second kind with two knots in which values for the function are assigned, an in the first knotvalues for the derivatives of the first two orders.

Constructing a spline of smoothness С 0
The matrices employed in the examined algorithm, when applying a polynomial with a power basis, are described in detail in the paper of its developer [9].To maintain a unified approach, we shall rewrite basis functions of the N. D. Dicoussar polynomial (1) in a local coordinate system with the origin at point x a =0.Then the other two knots will have the coordinates: x b =Mh, x 0 =qh, where q=[M/2], [ ] is a whole part from the number, h is the distance between knots.Then the basis functions will be assigned by formulae: We shall write respective matrices when applying the N. D. Dicoussar polynomial (6).Basis function d 1 , which is associated with knot x a , will be used under conditions of smooth joining, and coefficients at functions d 2 , d 3 and Q will be searched for by the method of least squares.
A condition for joining the l-th and (l+1)-th links of the spline with a smoothness of order С 0 will be consistently written in the form: When implementing a least squares method, we shall use a local numbering of knots for index i, assuming x a =X 0 , X i+1 =X i +h.
A matrix of the least squares method takes a traditional form: or according to denotation applied, for example, in paper [12]: Thus, we obtain a recurrent formula to calculate coefficient +1 l a f from the next link of the spline through the value of coefficient l a f from the current link of the spline, which ensures smoothness of joining the links of the spline of order С 0 : where is the stability matrix.
To examine stability of the algorithm depending on the form of a polynomial representation, we shall tabulate the values: -of the largest module of eigen number of stability matrix U; -of the conditionality number of matrix of the least squares method A 1 , for different ratios of values of M to m.
Table 2 will include ratios between M and m, which lead to the lowest of the largest modules of eigenvalues of stability matrices U at a fixed value of M. We consider h=1 during tabulation.Data from Table 2 testify to significant benefits of polynomials in the N. D. Dicoussar form when applying them for the problem on the construction of a smoothing spline.These benefits are ensured by a significant decrease in the conditionality number of the matrix of least squares method А 1 .

5. Constructing a spline of smoothness С 1
It is obvious from the geometrical content of coefficients of the N. D. Dicoussar polynomial that the form of this polynomial representation leads to significant complication in the algorithm for constructing an approximating spline with a smoothness order larger than С 0 .Given the relation between the N. D. Dicoussar polynomials and the Hermitian polynomials, established in chapter 4, we shall explicitly record the matrices used in the D. A. Silaev algorithm for the Hermitian polynomial (4).
We shall also rewrite the basis functions of polynomial (4) in the local coordinate system with origin at point x a =0.In addition, we shall change the numbering order for basis functions in order to simplify the system of equations describing conditions for the smoothness of the spline.Thus, polynomial (4) takes the following form: where The basis functions of polynomial (8) will then be assigned by formulae: Coefficients at basis functions u 1 and u 2 will be found from the conditions for the smoothness of joining links of the spline, and the coefficients at functions u 3 and u 4 will be searched for by the method of least squares.
Condition for joining the l-th and (l+1)-th links of the spline with a smoothness of order С 1 will be consistently written in the form: A matrix of the least squares method takes a traditional form: or according to the notation used, for example, in paper [12]: Hence, we obtain Recurrent formula (7) takes the following form in this case: ( ) ( ) and ensures smoothness of joining links of the spline of order С 1 .
We shall consider another form of representation of the Hermitian polynomial with three knots in which values for the function are assigned, in the first knotalso the value of the function's derivative: where The basis functions of polynomial (13) will then be assigned by formulae:

qh x M M q h
Coefficients at basis functions n 1 and n 2 will also be found from the conditions for smoothness of joining links of the spline, and coefficients at functions n 3 and n 4 will be searched for by the method of least squares.
Condition for joining the l-th and the (l+1)-th links of the spline with a smoothness of order С 1 will be consistently written in the following form: Other calculation formulae are similar to formulae ( 9)-( 12), which is why we shall not write them.
Results of studying the stability of an algorithm depending on the form of a polynomial representation are given in Table 3.It is obvious that it is appropriate to use the Hermitian polynomials in the form (13), which, all other conditions being equal, ensure the best conditionality of matrix A 1 , which is applied in the implementation of the examined algorithm.

6. Constructing a spline of smoothness С 2
To construct the smoothing spline of order С 2 , we shall use the Hermitian polynomial with two knots, in which values for the function and its two derivatives are assigned in the first knot, and in the second knotonly a value for the function: where The basis functions of polynomial (14) will then be assigned by formulae: Coefficients at basis functions w 1 , w 2 and w 3 will be found from conditions fpr smoothness of joining links of the spline, and coefficient at function w 4 will be searched for by the method of least squares.
Condition for joining the l-th and the (l+1)-th links of the spline with a smoothness of order С 2 will be consistently written in the form: A matrix for the least squares method takes a traditional form: or according to the notation used, for example, in paper [12]: Recurrent formula (7) takes the form in this case: and ensures smoothness for joining the links of spline of order С 2 .
Results of studying the stability of an algorithm, depending on the form of a polynomial representation, will be given in Table 4.As matrix А1 is composed in this case of one element, its conditionality number for all the bases is unity.It is obvious from data in Table 4 that the application of Hermitian polynomials in the form (14) has advantages in comparison with polynomials with a power basis because it ensures better stability of the algorithm to the impact of measurement errors in experimental data.

Discussion of results on a possibility of the generalization of the examined algorithm for irregular grids
It is clear from the results of theoretical research that when constructing the spline of smoothness С 0 it is appropriate to apply polynomials in the N. D. Dicoussar form; that of smoothness С 1 -the Hermitian polynomials of second kind (13).To construct the spline of the examined kind of smoothness С 2 , it is expedient to apply polynomials of higher power, since when applying the third-power polynomials, in order to determine by the method of least squares, there is only one coefficient, the model turns out to be too rigid and the stability of the algorithm is significantly worse than in the two previous cases.
To generate an experimental sequence, we shall construct a function that is a combination of the three Lorentz functions with parameters that are given in Table 5.
where A j is the amplitude; w j is the semi-width; x j0 is the frequency position of a maximum.We shall add to the values of function y(x i ) an error that has a normal distribution law with parameters a=0; s=0.075.
We shall generate a sequence of N=120 in the interval [−3 ; 5].Using rule 2s, we shall add points to each spline's link until the polynomial starts to deviate from experimental points by larger than 2s.The final number of points at which this requirement is fulfilled will be designated (М+1).
In the first interval, we shall solve an auxiliary task on approximation with additional constraints that take into consideration the physical essence of the problem.Specifically: the approximating polynomial at the beginning of the segment accepts a positive value and has a positive first derivative at the same point.Under such conditions, the algorithm selected 12 links of the spline with smoothness С 0 with the polynomials from the current link of the spline represented in the N. D. Dicoussar form (1) (Fig. 2).As follows from the conditions for the construction of such a spline, its derivative at the points of connection has gaps (Fig. 3).
Computational experiments to construct the spline of smoothness С 1 with the representation of polynomials from current links in the second Hermitian form (13) revealed the following.The best approximation results, in terms of their robustness, are obtained when the imposing depth of the spline's links increases by one point relative to the number of points that are determined based on the estimates of eigenvalues modules of stability matrices.In this case, the algorithm selects 14 links (Fig. 4).It is clear that the spline derivative chart in this case is continuous.
Fig. 2-4 show that at segments with a slow change in the examined parameter the algorithm selects longer spline links than at sections with a fast ascending or descending of experimental dependence.As expected, this makes it possible to reduce the total number of the spline's links compared to their number at uniform splitting.Thus, while ensuring an acceptable quality of approximation, the resulting spline requires less memory for its storing and using.
The disadvantages that remain inherent to the algorithm after its generalization include subjective determining of the acceptable magnitude of deviation in the chart of a polynomial from the current link of the spline from experimental data.Note that in the case of an unknown probabilistic law of distribution (or its parameters), by which measurement errors abide, such a disadvantage is common to all the algorithms for constructing the smoothing splines of such kind.
Prospects for the further research include the extension of algorithm to use, in its composition, polynomials of higher degrees to ensure the second order of smoothness of joining links of the spline with the simultaneous avoidance of the emergence of parasitic oscillations at an increase in the length of experimental dependence.It is also expedient to devise a formal criterion for determining the length of the current link of the spline in order to ensure the property of robustness for a smoothing spline.Further research will also relate to studying the question of alignment of the specified criterion with the criterion of algorithm stability.

Conclusions
1.It is shown that in a transition from representing a polynomial in a power basis to the representation in the basis of the N. D. Dicoussar polynomial (when constructing a spline of smoothness С 0 ) or to the representation in the bases of the appropriate Hermitian polynomials (in the construction of a spline of smoothness С 1 and С 2 ), there is a significant reduction in the conditionality number of the Gram matrix for the basis functions of the recommended forms of a polynomial representation.This ensures better approximating properties of the examined spline.
2. It is shown that the N. D. Dicoussar polynomial is a special case of the Hermitian polynomial, which is represented in a hierarchical form.
3. It was established that the form of a polynomial representation affects the magnitude of eigenvalues of the stability matrix, starting at the second order of smoothness while joining links of the spline.Such an impact is missing at the lower orders of smoothness of joining the links of the spline.
4. In all the test examples, investigated by Authors, the use of constraints, suggested by N. D. Dicoussar for preliminary determining the length of a link of the spline, made it possible to obtain a smoothing spline, based on the generalized D. A. Silaev algorithm, of acceptable quality, with links of different lengths.
5. An analysis of the quality of the obtained solutions to test problems revealed that in the processing of experimental dependences with significant measurement errors, there may occur the need to reduce the final lengths of spline's links compared to the lengths that are determined by the D. A. Silaev rule.Therefore, in the presence of significant errors in experimental data, it is appropriate to run a preliminary analysis employing the methods of mathematical statistics.This would make it possible, when executing the algorithm, to pay special attention to ensuring its stability.