Subdivision schemes based on weighted local polynomial regression. A new technique for the convergence analysis

The generation of curves and surfaces from given data is a well-known problem in Computer-Aided Design that can be approached using subdivision schemes. They are powerful tools that allow obtaining new data from the initial one by means of simple calculations. However, in some applications, the collected data are given with noise and most of schemes are not adequate to process them. In this paper, we present some new families of binary univariate linear subdivision schemes using weighted local polynomial regression. We study their properties, such as convergence, monotonicity, polynomial reproduction and approximation and denoising capabilities. For the convergence study, we develop some new theoretical results. Finally, some examples are presented to confirm the proven properties.


Introduction
In past years, many techniques have been designed and developed in order to construct curves or surfaces with some properties such as polynomial reproduction or monotonicity-preservation. For example, splines, non-uniform rational B-splines (NURBS) and others (see, e.g.[2,8]).In this context, linear subdivision schemes appears as useful and efficient instruments due to their simple computation (see e.g.[7,17,22]).They consist in obtaining new points from given data using refinement operators and can be classified depending on such operators: if a single operator is used for all the iterations, then the subdivision scheme is called stationary or level-independent (see e.g.[7,14]), otherwise it is denominated non-stationary or level dependent (see e.g.[9,10,16]).They are also classified by the linearity of the operators (see e.g.[12,13]).
There is a vast literature on the generation of subdivision schemes and the study of their properties.An essential property is convergence, which means that the process converges uniformly to a continuous function, for any initial values.Deslauriers and Dubuc, in [11], analysed that the scheme based on centered Lagrange interpolation is convergent using Fourier transform techniques to prove it.
One of the most common studied properties is the reproduction of polynomials, i.e., if the given data are point-values of a polynomial, then the subdivision scheme generates more point-values of such polynomial.This is studied in detail in [15].Its study is interesting since the reproduction is linked with convergence properties and the approximation capability of the scheme.
In some real applications, the given data come from measures that are contaminated by noise and, as a consequence, a suitable subdivision scheme should be used to converge to an appropriate limit function.To this purpose, Dyn et al. in [14] propose a new linear scheme based on least-square methods where the noise is reduced by applying the scheme several times.These schemes are determined by two parameters m and d with d < m: For each m consecutive data values, (y 1 , . . ., y m ), attached to some equidistant knots (x 1 , . . ., x m ), a polynomial regression is performed.The search is constrained to polynomials of degree d and leads to a unique solution, p, that minimizes the regression error concerning the ℓ 2 -norm (least-squares): p = arg min (y l − p(x l )) 2 . ( The subdivision refinement rules can be obtained by evaluating p at a certain point, which, in this work, is assumed to be 0 without loss of generalization.The resulting schemes are linear, which implies some benefits and drawbacks. In [14], the convergence is proved for d = 0, 1, as well as some properties such as polynomials reproduction.In many applied situations, the location of the data is relevant to obtain the approximation, hence a weight function is considered to assign values depending on the distance from the knots x l to 0. These methods, as Shepard's algorithm (see [24]), are called moving least squares (see [20]).In [4,5], the weighted local polynomial regression (WLPR) was used to design a prediction operator for a multiresolution algorithm, leading to good results on image processing when the data was contaminated with some noise.Prediction operators can be considered subdivision operators and their properties can be studied [9].In this paper, we study the family of subdivision schemes based on the prediction operators in [4] and develop a new technique to study their convergence based on some asymptotic behaviour.Also, some properties such as polynomial reproduction, the Gibbs phenomenon in discontinuous data, monotonicity preservation and denoising and approximation capabilities are analysed.We provide some examples to check the theoretical results.
The paper is organized as follows: Firstly, we briefly review the classical components of linear subdivision schemes with the aim to be self-contained in this work.In Section 3, we explain the WLPR and define a general form, leading to new subdivision schemes definitions.Afterward, we study different properties in some particular cases: Starting with d = 0, 1, we analyse the convergence, the smoothness of the limit functions, the monotonicity preservation and the Gibbs phenomenon when the initial data present large gradients.In Section 5, we develop a new technique to study the convergence of a family of schemes and apply it to the case d = 2, 3. We analyse the approximation and noise reduction capabilities of the new schemes in Sections 7 and 8. Finally, some numerical experiments are performed to confirm the theoretical properties, in Section 9, and some conclusions and future work are proposed.

Preliminaries: A brief review of linear subdivision schemes
Let us denote by ℓ ∞ (Z) the set of bounded real sequences with indices in Z.A linear binary univariate subdivision operator S a : ℓ ∞ (Z) → ℓ ∞ (Z) with finitely supported mask a = {a l } l∈Z ⊂ R is defined to refine the data on the level k, f k = {f k j } j∈Z ∈ ℓ ∞ (Z), as: In this work, we only consider level-independent subdivision schemes, meaning that the successive application of a unique operator S a constitutes the subdivision scheme.Hence, we will refer to S a as the subdivision scheme as well.
The binary adjective refers to the two formulas/rules of (2) (corresponding to i = 0 and i = 1) which are characterized by the even mask a 0 = {a 2l } l∈Z and the odd mask a 1 = {a 2l−1 } l∈Z .It is called length of a mask to the number of elements that are between the first and the last non-zero elements, both included.
Remark 2.1.If a linear subdivision scheme is applied to some data g = {G(j) + ǫ j } j∈Z , where G is a smooth function and ǫ = {ǫ j } j∈Z is random data, also called noise, the result is S a g = S a g + S a ǫ, which implies that we can study separately the smooth and the pure noisy cases.
If we apply these rules recursively to some initial data f 0 , it is desirable that the process converges to a continuous function, in the following sense.Definition 1.A subdivision scheme S a is uniformly convergent if for any initial data f 0 ∈ ℓ ∞ (Z), there exists a continuous function F : R → R such that Then, we denote by S ∞ a f 0 = F to the limit function generated from f 0 .We write S a ∈ C d if all the limit functions have such smoothness, S ∞ a f 0 ∈ C d , ∀f 0 ∈ ℓ ∞ (Z).A usual tool for the analysis of linear schemes is the symbol, that we define as follows.
Definition 2. The symbol of a subdivision scheme S a is the Laurent polynomial a(z) = j∈Z a j z −j .
A useful property for a subdivision scheme is the reproduction of polynomials.A necessary condition for convergence is the reproduction of constants.The following lemma determines the relation between the mask, the symbol and the reproduction of the constants.
In such case, the S q scheme is well-defined and called difference scheme.If S q < 1, then S a is convergent.
There exists a direct relationship between the symmetry of S a and the symmetry of its difference scheme, S q .We introduce it in the following result.
Lemma 2.3.If a scheme is odd-symmetric, then its difference scheme is even-symmetric.
Proof.It can be easily checked using the symbols.
Next theorem by Dyn and Levin, [17], links the smoothness of S a and S 2q .
Theorem 2.4.If the scheme based on S 2q is convergent and C m−1 , then S a is convergent and C m .Remark 2.2.We give now a more explicit formula to compute q for the kind of schemes we consider in this paper.We will analyze odd-symmetric subdivision schemes, which implies that the length of the mask is always odd and two possible situations may occur, depending on which sub-mask has the largest support.
Since the sub-masks a 0 and a 1 are finitely supported, from now on we will treat them as vectors containing only their support, which will be important for the theoretical results in Section 5.The first situation is that, for some n ∈ N, the sub-masks are a 0 = {a 0 l } n−1 l=1−n and a 1 = {a 1 l } n l=1−n , while the second one corresponds to a 0 = {a 0 l } n l=−n and a 1 = {a 1 l } n l=1−n (pay attention to the supports).To compute q with a unique formula for both cases, we redefine the mask for the second case, consisting in ā0 l := a 1 l , l = 1 − n, . . ., n, and ā1 l := a 0 l−1 , l = 1 − n, . . ., n + 1.Now the first indices of the supports are 1 − n, in both situations, and the last indices are n − 1 and n (for the first and second sub-mask, respectively) for the first situation and n and n + 1 for the second one.Now, in both cases, the second sub-masks is the largest and we can affirm that there exists some n ∈ N such that with L n = n − 1 or L n = n, so that a 0 = {a 0 l } Ln l=1−n and a 1 = {a 1 l } Ln+1 l=1−n .In any case, now the odd-symmetry is written as Finally, for a subdivision operator written as (3), the difference mask q can be computed as follows: According to Lemma 2.3, S q is an even-symmetric scheme.In particular,

Weighted local polynomial regression (WLPR)
The schemes analysed in the present work has been applied to image processing in a multiresolution context as prediction operator both for point-values as for cell-average discretizations, (see, e.g.[4,5]).They are based on weighted local polynomial regression (WLPR) and they can be defined by inserting a weight function in the minimization problem (1), which emphasizes the points closer to where the new data is attached.In this section, we briefly introduce WLPR and describe some of its properties.For a more detailed description, see [19,21].
Firstly, we fix the space of functions where the regression is performed: Π d , the space of polynomials of degree at most d.Other function spaces could be used as well (see [19]).We can parametrize the polynomials in Π d as where the superscript T is the matrix transposition, A(x) T = (1, x, . . ., x d ) and β ∈ R d+1 .The vectors are considered column vectors in order to perform the matrix multiplication.With this notation, the regression problem (1) can be expressed as The second ingredient is the weight function, ω : R → [0, 1], which assigns a value to the distance between x i and 0, which is the location where p is evaluated in this work.We define ω as and we impose that φ : [0, 1] → [0, 1] is a decreasing function such that φ(0) = 1.With these assumptions it is clear that ω has compact support, [−1, 1], is even, increasing in [−1, 0] and decreasing in [0, 1], and it reaches the maximum at point x = 0.The choice w(0) = 1 assigns the highest weight to the point where p is evaluated.In [21], some functions are proposed, which we compile in Table 1.Observe that many of them have the form φ(x) = (1 − x p ) q with p, q > 0.
The third component is the bandwidth, λ ∈ R + \N.We define Table 1: Weight functions, see [21].
The parameter λ determines how many data values are used in the regression and allows to distribute the weights of the points used in the rank [−λ, λ].By the properties of the function ω, if λ 1 ≤ λ 2 , then w λ1 l ≤ w λ2 l for any l ∈ Z. Finally, we choose a vector norm, typically ℓ 2 is taken for its simplicity, but any ℓ p -norm can be used depending on the characteristics of the problem.The loss function is defined accordingly: With the above elements, we propose these two problems to design the two subdivision rules: Once the fitted polynomial is obtained, it is evaluated at 0 to define the new data: so that only the first coordinate of βi is needed.

2
, the resulting subdivision scheme is the Deslauriers-Dubuc subdivision scheme.
Proof.Let us discuss when this scheme is well defined.Two situations may occur, depending on whether or not d (the polynomial degree) is smaller than the amount of data f k j+l in the minimization problem (6).6) is a least square problem and there is a unique solution [6], otherwise it can be found a polynomial that interpolates the data.Even if the interpolating polynomial is not unique, its evaluation at 0 is exactly f k j .Hence, the even rule is well defined for any λ ∈ R + \N, coinciding with the even rule of the Deslauriers-Dubuc subdivision scheme for d ≥ 2 λ 2 , i.e.
2 and an interpolation problem with unique solution is solved when the equality is reached, coinciding with the Deslaurier-Dubuc odd rule in the last case.However, nor the polynomial neither its value at zero are unique when d + 1 > 2 λ+1 2 , so that the scheme is not well defined in this case.
As conclusion, only if the polynomial degree is , the resulting scheme is the Deslauriers-Dubuc interpolatory subdivision scheme, independently the choice of ω and the loss function L p .
The scheme (7) coincides with the proposed by Dyn et al. in [14] if p = 2 and φ(x) = 1 are used (corresponding to rect in Table 1).Also, the non-linear subdivision scheme presented by Mustafa et al. in [23] can be obtained with the same choice of φ(x) = 1 but with p = 1.
We will analyse the properties of our schemes specifically for the polynomial degrees d = 0, 1, 2, 3, the loss function L 2 and several choices of φ.
We will study how the choice of φ affects the approximation and noise reduction capabilities.We will show that it is not possible to define a φ giving the best approximation and the greatest denoising.In fact, one may decide how much importance to adjudge to each property and find an equilibrium.This decision may be based on the magnitude of the noise and the smoothness of the underlying function.
Observe that, when 2n − 1 < λ < 2n, for some n ∈ N, the even rule (i = 0) support is shorter than the odd (i = 1) one, and just the opposite occurs when 2n < λ < 2n + 1.To simplify, we will discuss in detail the first case, where even and odd masks have lengths 2n − 1 and 2n, respectively, since the second one is analogue and the Remark 2.2 can be taken into account for the consequent analysis.Nevertheless, we deal with both situations along the paper when it can be do it without additional effort.
To give a more explicit definition of the schemes, we solve the quadratic problem posed in (6) with p = 2.In this case, it is a weighted least square problem and its solution is well-known.Let us start with the derivation of the odd sub-mask, a 1 , for 2n − 1 < λ < 2n + 1.For the sake of simplicity, we omit the dependence on d, ω, λ for the following vectors and matrices.If we denote as W 1 the diagonal matrix consisting on the vector we call where the powers (x 1 ) t , t = 0, . . ., d, are computed component-wisely, so that X 1 is a 2n × (d + 1) matrix, and we denote f 1,j,k = (f k j−n+1 , . . ., f k j , f k j+1 , . . ., f k j+n ) T , then the problem of ( 6) can be write as: For the sake of clarity, we write down the above terms: and Since we only need the first coordinate β1 0 , we can use the Cramer's formula instead of solving the full system: Observe that, since the vector w 1 is symmetric, 1) t = 0 for any odd value of t, and t for the even values.Thus, the above expressions can be simplified by placing many zeros and by shorting the range of the remaining sums.Using the linearity of the determinant respect to the first column, , we conclude that the sub-masks coefficients are . By (10), it can also be expressed as , where e 1 is the first element of the canonical basis of R d+1 .Analogously, for 2n − 2 < λ < 2n, we can prove that a 0 = W 0 X 0 ((X 0 ) T W 0 X 0 ) −1 e 1 , so that , where W 0 is the diagonal matrix with diagonal and x 0 = (−2(n − 1), . . ., 2, 0, 2, . . ., 2(n − 1)) T , X 0 = (x 0 ) 0 , (x 0 ) 1 , . . ., (x 0 ) d .
Collecting these developments, for 2n − 1 < λ < 2n, we can define our weighted local polynomial regression-based subdivision as: A direct consequence, by construction, is that the scheme reproduces polynomials up to degree d.
Remark 3.1.Observe that we have considered {1, x, . . ., x d } as basis of Π d , which has led a the linear system with matrix (11).It is possible to consider an orthonormal basis of Π d in a way that the matrix is diagonal, leading to a cleaner mathematical description.However, we preferred the basis {1, x, . . ., x d } because the resulting expression of the subdivision operator is more explicit.A possible benefit of considering an orthonormal basis is that the next results might be more intuitive.Now we prove that (W i ) −1 a i , i = 0, 1, are exactly the evaluations of some polynomial at the grid points x i .
Lemma 3.3.For i = 0, 1, the sub-masks are That is, the vector (W i ) −1 a i coincides with the evaluation of the polynomial A(x) T α i at the points x i , being α i = ((X i ) T W i X i ) −1 e 1 , which expression depends on n, λ, ω.
Proof.By the previous computations, For a 1 (for a 0 is analogous), using (9) we obtain That is, the coordinates of (W 1 ) −1 a 1 are the evaluations of the polynomial A(x) T α 1 at the x 1 grid points.
Moreover, these sub-masks are the only ones that lead to polynomial reproduction and verify that (W i ) −1 a i are polynomial evaluations.This property can be used in practice to easily determine the sub-masks, as we do in Section 6.
Theorem 3.4.The scheme S d,w λ is the unique scheme that reproduces Π d polynomials and its sub-masks have the form a Proof.It is a consequence of Lemma 3.3 together with Proposition 3.2.Suppose that some rule â = W i X i α, for some α ∈ R d+1 , fulfils the reproduction conditions for Π d .Then or, written with matrix multiplications, (X i ) T â = e 1 .Then, The symmetry of the scheme is another consequence of being based on a polynomial regression problem.
Proof.We prove that a 1 j = a 1 1−j for 2n − 1 < λ < 2n + 1 (it can be analogously proven that a 0 j = a 0 −j for 2n − 2 < λ < 2n).Let us consider f j = {δ j,l : l ∈ {−n + 1, . . ., n}}.The coordinates of the sub-mask a 1 can be obtained by applying the rule to f j , for j = −n + 1, . . ., n, and take the first coordinate, where, by ( 6) and ( 7), pj = arg min Then, a 1 j = a 1 1−j provided that (S d,w λ f j ) 1 = (S d,w λ f 1−j ) 1 , or in other words, pj (0) = p1−j (0).Observe that, p1−j = arg min and, performing the change in the summation index l by 1 − l and using Observe the similarity between ( 14) and (15).Since the minimum is unique, it is reached in (15) by pj (−t).Thus p1−j (t) = pj (−t) and, then, By Lemma 3.3, we know that (W 1 ) −1 a i are the evaluations of a polynomial at x i .To take profit of the symmetry, let us write as in ( 13): Since and ω is even, it can be deduced that the polynomials only have even powers.That is A direct consequence is that the subdivision schemes obtained for any weight function of degree d (even number) coincides with the one for d + 1, proven in the following lemma.Proof.The sub-masks of S d+1,w λ can be written in terms of the evaluation of a (d + 1)-degree polynomial, according to Lemma 3.3.Since the odd coefficients are zero, then the leading coefficient is zero, for both rules i = 0, 1.Then, both S d,w λ and S d+1,w λ fulfils the conditions of Theorem 3.4, for the same polynomial degree d, hence they must coincide.
Therefore, we can just study the properties of the subdivision schemes based on the space of polynomials Π d (R) with d an even number.

WLPR-Subdivision schemes for d = 0, 1
In this section we present the WLPR-Subdivision schemes for d = 0, 1 and their properties, by Proposition 3.6, we can just consider d = 0. To simplify the notation, in this section we omit d, w and λ in some variables, such as S := S 0,w λ = S 1,w λ .In this case, the coefficients of the subdivision schemes are easily obtained from ω thanks to Lemma 3.3: If we denote as ||w i || 1 the sum of the components of the vector w i with i = 0, 1, defined in ( 12) and ( 8), and, for 2n < λ < 2n + 1, it can be written in the following way, in agreement with Remark 2.2, Note that if λ ∈ (1, 2) then w 0 = 1 and w 1 = ( 1 2 , 1 2 ), so that the mask for any function ω of the subdivision scheme is a = [1, 2, 1]/2, in other words, the interpolatory Deslauriers-Dubuc scheme [11] (as stated in Proposition 3.1).For λ > 2, if ω(x) = 1 for |x| ≤ 1, then the schemes presented by Dyn et al. in [14] are recovered as we mentioned above.These schemes are for 2n − 1 < λ < 2n: We list some masks for the weight function ω(x) = 1 − |x|, |x| ≤ 1, and for several values of λ: a 0,tria As we can see, all subdivision schemes in this section present a positive mask, since ω is a positive function.Then, the following result on convergence proved in [18], (see also [22,26]) can be applied.
As a direct consequence, the schemes in this section, ( 17) and ( 18), are convergent because the masks are positive.Observe that the condition k ≥ 3 in Proposition 4.1 requires considering λ > 2.
To analyse the smoothness of the limit functions generated by S, we consider the Theorem 2.4.In particular, we will prove that the mask of the difference scheme S q is positive and apply again Proposition 4.1.Thanks to the odd-symmetry of the scheme, the study can be reduced to a half of its coefficients.Lemma 4.3.Let n be a natural number, n ≥ 2, λ ∈ (2n − 1, 2n) and ω a weight function.The coefficients of the difference scheme S q are positive if Proof.By Lemma 2.3, S q is even-symmetric.Since 2n − 1 < λ < 2n, then L n = n − 1 in (5) and we have Then, if q 0 j > 0, for j = 1 − n, . . ., 0, and q 1 j > 0, for j = 1 − n, . . ., −1, the result is proved.First, we check that the coefficients q 0 0 and q 1 1−n are always positive.
Proof.Note that: Consider this basic property: For any a, b, c, d > 0, Firstly, since p ω λ 0 is decreasing, we get by (23): And, again using the monotony of function p ω λ 0 and ( 23): Repeating this process, we get by ( 23): Secondly, we define which is an increasing function since p ω λ 0 is decreasing.We have that Again, using the same strategy, we get: Then, by Lemma 4.3, we conclude that the coefficients of the difference scheme, (4), are positive.
Next lemma allows to easily check the monotonicity of p ω λ 0 .
Proof.By hypothesis, the function is continuous for [0, n − 1] and differentiable in (0, n − 1) (observe that we are considering l a real number here).Hence, it is decreasing provided that its derivative is negative.In addition, and φ ′ /φ is decreasing by hypothesis.
Therefore, we prove the following corollary.
Proof.From Table 2, the function p ω λ 0 is decreasing for any φ(x) = (1 − x p ) q with p ≥ 1 and q > 0.Then, by Lemma 4.3 the coefficients of the difference scheme are positive and by Proposition 4.1 the subdivision scheme S 2q is convergent.The case φ(x) = 1 is studied in [14].
In order to finish this section, we study two properties.Firstly, we analyse if the new family of schemes conserves monotonicity.In our case, the result presented by Yad-Shalom in [25] can be used: Proposition 4.7.Let S a be a convergent subdivision scheme and S q its corresponding difference scheme with a positive mask.If the initial data, f 0 , is non-decreasing then the limit function S ∞ f is non-decreasing.
With this proposition, we can enunciate the following corollary.
Finally, when the initial data presents an isolated discontinuity and a linear subdivision scheme is applied several times some non-desirable effects may appear near the discontinuity, some kind of Gibbs phenomenon (see e.g.[1]).
In [1] it is proved that if the mask of the scheme is non-negative then the Gibbs phenomenon does not appear in the limit function.
In Section 9, we present some examples checking these theoretical results.For d = 0, 1, the resulting mask is positive and we have used classic tools to study its properties.However, for d ≥ 2, the mask are no longer positive.In the next section, we will develop a novel technique based on numerical integration for this goal and we will apply it to prove the convergence of the schemes based on weighted-least squares.

A tool for the convergence analysis
The purpose of this section is to provide new theoretical results to analyse the convergence.In Section 4, the convergence was easily proven by the positivity of the mask.However, in Section 6 we will prove the convergence of the scheme based on the regression with polynomials of degrees d = 2, 3, which are no longer positive, so that we cannot follow the same strategy.Nevertheless, as a consequence of Lemma 3.3, the sub-masks can be seen as the evaluation of a second degree polynomial and this fact is advantageous and we will take profit of it in this section.
For any particular value of n, a fixed ω and considering some λ n such that 2n − 1 < λ n < 2n + 1, λ n = 2n, it can be easily computed the difference scheme using the formula (4) and checked if its norm is less than 1, which would imply convergence.Let us call this method the direct inspection.But it serves to prove convergence only for the chosen n, and we wish to prove it for all n ∈ N. Our strategy will consist in proving converge asymptotically, that is, to prove convergence for ∀n > n 0 , for some n 0 ∈ N, and then check the converge for each n ≤ n 0 by direct inspection.
First, we would like to give a general idea about this asymptotic convergence.Thanks to the properties of the space of polynomials Π d , the problem (6) can be formulated using equidistant knots in the interval [−1, 1], such as βi = arg min The last sum is, in fact, a composite integration rule.So that, if n → ∞, then 2n/λ n → 1 and the problem seems (this is not a rigorous argument, but it serves to understand the situation) to converge to arg min for both i = 0, 1.On the one hand, the given data is now a function f (x) which is approximated by a polynomial A(x) T β ∈ Π d in the L p norm with a weight function ω.On the other hand, by Lemma 3.3 the corresponding subdivision sub-masks, say a n,i , fulfils , for some coefficients α n,i ∈ R d+1 .Then, the sub-masks also seem to converge to some continuos function, if some normalization is performed since the sub-masks supports increase with n (see later Remark 5.1 and Section 6 for more details).The results presented in this section exploit this kind of situations.
From now on, we consider a family of subdivision schemes {S a n } ∞ n=1 as in (3).The results in this section allow to prove convergence for n > n 0 , for some n 0 ∈ N, and also provides the value of n 0 , so that it can be checked convergence for n ≤ n 0 by direct inspection.Combining both proofs, we obtain convergence for all n ∈ N. In particular, lim n→∞ S q n ∞ will be computed, which ensures the asymptotic convergence when that limit is less than 1.Here we denote by a n,0 , a n,1 , q n,0 , q n,1 the sub-masks of the masks a n , q n .Theorem 5.1.Let {S a n } ∞ n=1 be a sequence of subdivision schemes that reproduces Π 0 , which odd rules are longer than (or as long as) the even rules, as in for some α > 2, µ > 0, then the first sub-masks of the difference schemes fulfil thus there exists n 0 ∈ N such that q n,0 Moreover, if (25) holds true for α = 3, then where |r ′ (t)|.
Proof.First, we may write q n,0 j in terms of r: Using the composite (backward) rectangle rule, we obtain where θ n j is the integration error, which fulfils With this computation, we will prove (27) first: we use that R(−1) = 0 and the composite (forward) rectangle rule, thus obtaining that where ρ n is the integration error of R(t), If L n = n, we use the composite (backward) rectangle rule and we obtain a similar result: Using all the upper bounds we found, we obtain: From here we deduce that, if α > 2, then the limit when n → ∞ of the right part of (29) is R 1 , which is less than 1.Hence, there exists n 0 ≥ 1 such that q n,0 1 < 1, ∀n > n 0 .In particular, for α = 3, we can find for which value of n 0 the right part of (29) is equal to 1, by solving a second degree equation, arriving to (28).
Remark 5.1.In practice, if the expressions of a n,0 j , a n,1 j are well defined for any j ∈ R (this is the case of S 3,w λ , see (35)), then a practical way to compute r(t) is In Section 6, a complete example of the application of the results of this section will be performed.
On the other hand, S n a reproduces Π 0 if, and only if, S ān does.Hence, the finite difference scheme exists and can be computed with the formula (4).qn,0 Hence, Since R(t) = R(−t) and r(t) = r(−t), we deduce that the formula to compute n 0 , (28), can be used here as well.
The next result is a direct consequence of the previous one.
for some α > 2, µ > 0, then there exists n 0 ∈ N such that In case that α = 3, n 0 can be obtained as in (28).
Proof.By the Theorem 5.2, the flipped version of this scheme fulfils Theorem 5.1 and the claimed inequality is true.
For odd-symmetric subdivision operators, due to Theorem 5.2, the satisfaction of the hypothesis of Theorem 5.1 or Corollary 5.3 is sufficient to ensure convergence.Theorem 5.4.Let {S a n } ∞ n=1 be a set of odd-symmetric subdivision schemes fulfilling the hypothesis of Theorem 5.1.Then, the subdivision scheme S a n is convergent if n > n 0 with n 0 as in (28).

WLPR-Subdivision schemes for d = 2, 3
We consider {λ n } n≥2 such that 2n − 1 < λ n < 2n, then L n = 1 − n.The following computations could be done for 2n < λ n < 2n + 1 as well.First, we compute the coefficients of S 3,w λ (denote it by S n from now on).According to Lemma 3.3, the sub-masks are a n,i = W i X i α i , i = 0, 1, where α i = ((X i ) T W i X i ) −1 e 1 .Then, to compute α we may solve the system We start with i = 1.Using (11) and the symmetry of w 1 and x 1 , Hence, using the Kramer's formula, the three coefficients of α 1 are: Then, by ( 16), the sub-mask coefficients are Similarly, where We first prove convergence ∀n ≥ 2 in the simplest case, φ(x) = 1, in order to be used to the new convergence analysis tools, and later we discuss the general case.In this case, w l = 1, 1 − 2n ≤ l ≤ 2n − 1, so that the mask coefficients can be simplified to It can be easily checked that these operators are odd-symmetric, which for sure we knew by Lemma 3.5.Hence, to prove convergence we can apply Theorem 5.4.
As conclusion, Then, for any n 1 ≥ 3, To compute n 0 , it is also necessary to compute: |r ′ (t)| = 15/2.Now, using formula (28) (case It is desirable to prove convergence for as much values of n as possible, so n 1 should be chosen such that n 0 is as small as possible, but greater or equal than n 1 , due to (37).We computationally found that the compromise is achieved for n 1 = 188, leading to n 0 ≃ 188.506.Hence, according to Theorem 5.4, the subdivision schemes are convergent for n ≥ 189.For smaller values of n, we have computationally checked that This symbolic computation is quick and without rounding errors, so this can be considered a rigorous proof of the convergence.
We can perform some additional computations in order to provide an upper bound of q n,0 1 = q n,1 1 valid for any n ≥ 2. According to (29), We checked that the right side is less than 29/42 for any n ≥ 2236, and we explicitly computed that for n ≤ 2236, q n,0 1 ≤ 29/42.As conclusion, q n,0 and the equality is reached only for n = 4.We tried to prove C 1 regularity with this technique by applying the results to the divided difference schemes, S 2q n , but they do not satisfy (24).

Convergence of the subdivision schemes based on weighted least squares with d = 2, 3 and a general function φ(x)
In this situation, we will study the convergence only for large n values, so that we will not calculate n 0 , because we have not been able to perform the direct inspection without specifying φ.In order to compute r(t) := lim n→∞ (a n,0 tn − a n,1 tn )n 2 , we will define a C 1 function U n j such that a n,i j = U n j (1 − i/2), i = 0, 1, which will allow to write For that purpose, we define σ λn (φ, x, k) and ω(x) = φ(|x|).Observe that the sub-masks ( 33) and (34) can be expressed as Thus, we may define the link function as )) (φ ′ may not exist at 0 or 1).To follow more easily the next computations, we write U n j (x) = φ( j+x−1 λn/2 )U num (x)/U den (x), where U num (x), U den (x) are the numerator and denominator that appear in the last formula.
Taking into account that we proceed to compute the derivative.
Finally, we proceed to compute r(t) = lim n→∞ n 2 2 (U n tn ) ′ (ξ t,n ).To this purpose, we define we observe lim n→∞ 2n/λ n = 1 and we use the following composite integration rule Defining σ λn (φ, x, 0) := n l=1 φ( l−x λn/2 ), so that ∂σ λn ∂x (φ, x, 0) = − 2 λn σ λn (φ ′ , x, 0), we note that and . Taking these comments into account and taking j = tn, we find out that Hence, Clearly, the former expression is valid provided that I 0 (φ)I 4 (φ) − I 2 (φ) 2 = 0. Fortunately, we can use the Schwartz's inequality for the inner product f, g := 1 0 f (x)g(x)φ(x)dx to deduced that We gather in Table 3 the computation of r(t) and R 1 for several choices of φ.Since R 1 < 1 for all of them, we conclude that, for n large enough, any of the corresponding subdivision schemes converge.We realized that the value of R 1 could be greater than one for some extreme choices of φ.An example is φ(x) = 1 + 1000x 2 , but thus kind of functions were discarded in Section 3 due to its practical meaning.
The next two sections are devoted to study the approximation and the noise suppression capability depending on the chosen weight function

Approximation capability
To study the approximation capability, we consider the subdivision scheme S d,w λ defined in (7) with d ≥ 0 and λ satisfying the conditions requested in Proposition 3.1.Let F ∈ C d+2 be and consider the initial data f h = {f h j } j∈Z with h > 0 and f h j = F (jh) , j ∈ Z.Let j 0 ∈ Z be any integer, we calculate the approximation error between (S d,w λ f h ) 2j0+i and F ((j 0 + i/2)h), with i = 0, 1, and analyse the largest contribution term.By Taylor's theorem, we have that there exist p i ∈ Π d such that:  Applying the subdivision operator and considering its polynomial reproduction capability, Therefore, the largest contribution to the approximation error is given by We conclude that if two linear schemes are given, with the same approximation order, then the scheme with lesser value of ) provides better approximators, in general.We observe that, if Since the proposed schemes are odd-symmetric, then H(t) = H(−t) and 1 −1 t d+1 H(t)dt = 2I d+1 (H) and the approximation error is given by 2h d+1 n d+1 I d+1 (H) which increases with h, n and I d+1 (H).We will test this formula in Section 9.2.Now, we explore how the selection of φ influences H, with the aim of determining which φ is the best from an approximation point of view.
For d = 0, 1, it is easy to compute H from the expression of a 0 , a 1 in ( 17) and (18).For instance, for 2n−1 < λ < 2n, Hence, 2I 2 (H) = I 2 (φ)/I 0 (φ).In Table 4, we see that the smallest values are reached for φ(x) = e −ξx with large ξ and for φ(x) = (1−x p ) q with large q or small p.We add for comparison H 2 2 , according to Section 8, the smaller it is, the greater is its noise reduction capability.We can see for any scheme that the greater is the approximation capability, the smaller is the noise reduction capability.As conclusion, approximation and noise reduction are incompatible, in this sense, and some equilibrium may be found.This is further discussed in Section 8.1.

Noise reduction
In this section, we study the application of a subdivision operator to purely noisy data, S a ǫ where all the values ǫ j follows a random distribution E, and are mutually uncorrelated.The results of this study can be applied to any data contaminated with noise due to Remark 2.1.
A direct result is that Since S a ∞ ≥ 1 for any convergent schemes, the best condition is reached for d = 0, 1, for which S d,w λ ∞ = 1, since the mask is positive.Hence, it cannot be concluded from this formula that the noise is reduced.
To reveal the denoising capabilities, a basic statistical analysis can be carried out.If the variance of the refined data is lesser than the variance of the given data, var(E), it indicates a reduction of randomness.Using that var(αX + βY ) = α 2 var(X) + β 2 var(Y ), α, β ∈ R, provided that X, Y are two uncorrelated random distributions, the variance after one subdivision step is Increases with p. Decreases with q.
Decreases with p. Increases with q.

Decreases with ξ
Increases with ξ   Hence, the variance reduction is given by For some schemes studied in this work, this quantity is: The last quantity is less than one owned to the constant reproduction and the positivity of the coefficients.In case that φ(x) = 1, then S 1,rect λ 2 2 = ⌊λ⌋ −1 = (2n − 1) −1 , which is the lowest value that can be obtained with a rule of this length.For d = 2, 3, φ(x) = 1 and 2n − 1 < λ < 2n, which maximum is achieved for n = 2 (i.e. 3 < λ < 4, corresponding to the interpolatory DD4 scheme), which is 1.Two results can be derived: First, if the variance is reduced in each iteration by a factor S d,w λ 2 2 < 1, then the limit function has variance 0. Second, since lim n→∞ S 3,rect λ 2 2 = 0, the noise tends to be completely remove when the mask support tends to ∞.
For any choice of φ(x), an asymptotic result can be given for the noise reduction using an argument similar to Section 7: so that the noise reduction factor behaves asymptotically as Under these assumptions, we observe that the noise is always removed after an iteration when n → ∞.In the Tables 4 and 5 we compute H(t) := lim n→∞ na i tn and the factor H 2 2 for several φ functions, d = 0, 1, 2, 3.

An equilibrium between approximating and denoising
We have seen that, in order to maximize the approximation and denoising capabilities, the values I 4 (H) and H 2 should be minimized.This is a multi-objective minimization problem, which solutions form a Pareto front that we have estimated using the MATLAB optimization toolbox.Here we will only consider the case d = 2, 3, but a similar analysis can be performed with d = 0, 1.
First, observe Figure 2-left.We find out that φ(x) = (1 − x p ) q is always more convenient than φ(x) = e −ξx , meaning that for each value of ξ there exists some pair (p, q) for which φ(x) = (1 − x p ) q approximates and denoises better than φ(x) = e −ξx .It can also be affirm that φ(x) = 1 is in the Pareto front and it the best for noise reduction and the worst for approximating.In the other extreme would be an interpolatory scheme, with the best approximation capability but the worst denoising power.
In conclusion, we recommend the use of rect to obtain the best denoising.However, with epan the noise increases by 11.11% while the approximation error is reduced by 44.44% compared to rect.If the approximation is desired to be prioritized, φ(x) = (1 − x 4 ) 5 is a good choice, since the noise increases by 31.58% while the approximation error is reduced by 71.43%, compared to rect.The rest of the (p, q) values related to Table 1 are near to be optimal and can be used as well for other approximating-denoising balances.We recommend to never use exp(ξ).
Just to mention that for d = 0, 1 similar conclusions can be obtained.For that polynomial degrees, the weight functions φ(x) = exp(−ξx) are also worse than φ(x) = (1 − x p ) q .The weight function epan is still Pareto optimal, but the pair (p, q) = (4, 5) is not.

Numerical experiments
In this section, we present some numerical examples to show how the new schemes work for the generation of curves.We check that the subdivision schemes are convergent for d = 0, 1, 2, 3 and that the curve present C 2 ) for several choices of φ.Thus, the lower x and y axis values, the better approximation and denoising capabilities, respectively.Blue, φ(x) = (1 − x p ) q for several values (p, q) pairs such that 1 ≤ p ≤ 20, 1  2 ≤ q ≤ 20; red, the Pareto front of the previous pairs; green, φ(x) = exp(−ξx) for 1  2 ≤ ξ ≤ 10.Right, the red line represents the pairs of values (p, q) for which φ(x) = (1 − x p ) q is Pareto-optimal.
(but not G 1 , meaning that kinks can be produced).We analysed the approximating and denoising capabilities to numerically validate the results in Sections 7 and 8.Only for d = 0, 1, we test the conservation of the monotonicity applying the schemes to fit a non-decreasing initial data.Finally, we perform a numerical test using the discretization of a discontinuous function and observe that the proposed methods avoid Gibbs phenomenon in the neighbourhood of an isolated discontinuity for d = 0, 1.

Application to noisy geometric data
We start with one of the experiments presented in [14] which consists of a star-shaped curve given by: with samples taken at t 0 j = jπ/25 with j ∈ Z.That is, we consider f 0 := F | t 0 , t 0 = {t 0 j } j∈Z , i.e. f 0 j = F (t 0 j ).Because of the periodicity of the function, we can focus on j = 0, . . ., 49.We add Gaussian noise in each component, defining f 0 = f 0 + ǫ σ with ǫ σ = {(ε σ,1 j , ε σ,2 j )} 49 j=0 , being ε σ,l j ∼ N (0, σ), l = 1, 2, j = 0, . . ., 49 and σ ∈ {0.5, 1}.In Figure 3, we illustrate the results only for two interesting choices of φ, according to the conclusions in Section 8.1.Nevertheless, the results obtained with the rest of weight functions are graphically similar and they are shown in detail in Table 6.
Without noise, the smaller is λ the more accurate are the results for any φ and d.Measuring the approximation error as where t 5 j := 2 −5 t 0 j , j ∈ Z.As expected, we can see in Table 6 that the approximation error is always smaller for d = 2, 3 than d = 0, 1.But also, if we sort the weight function by the approximation order, for any 0 ≤ d ≤ 3, they would be exactly in the same order as if we sort them by the theoretical approximation power in Tables 4 and 5. Nevertheless, it has to be taken into account that the results in Tables 4 and 5 have an asymptotic nature, when h → 0. These behaviours are also visible in the first row of Figure 3.
We measure the noise reduction capability of the schemes with the quantities S 5 ǫ 0.5 ∞ and S 5 ǫ 1 ∞ , in Table 6.The ones that show results closer to zero are the schemes with higher denoising capacity.In general, the noise is more reduced when λ is larger or d is smaller.Comparing the sorting of the weight functions by its theoretical denoising capability, according to Tables 4 and 5, and by the numbers in Table 6, we see that both orderings are the same, in general.Only in some particular cases this ordering is slightly changed.The reason may be that the results on Tables 4 and 5 are asymptotical, for λ → ∞.But also, the reduction is in terms of the variance of the statistical distribution.Hence, the same experiment should be repeated many times and the results averaged in order to obtain a more consistent comparison.S 1,rect 9.5 S 3,rect 15.5  S 1,p4q5 9.5 S 3,p4q5 Figure 3: Several subdivision schemes (by columns) applied to the star-shaped data in (39).In the first row, they are applied to the original data.In the second and third row, the data is contaminated by normal noise with σ = 0.5 and σ = 1, respectively.
In Figure 3, we can see how important the choice of the weight function is to increase the approximation capability (and only losing a bit of denoising capability).In turn, taking d = 2, 3 gives better approximations and λ can also be increased to reduce noise.
Of course, during our study we generated much more graphics than the ones here presented.In some of them, specially in presence of noise, artefacts may appear, such as auto-intersections or kinks, proving that it does not provide G 1 curves, even if the scheme is C 1 .By taking λ larger, the artefacts usually disappear and curves become softer.
This threshold is not a real constrain in practice, since the noise is usually greater than the approximation error (see first row of Table 7).If an approximation error tending to zero is needed, n ∝ h − 1 2 can be chosen, for instance.Table 7: The approximation error at x = 0 after five iterations of S 3,rect λ k applied to data with and without noise.

Avoiding Gibbs phenomenon
In this section we confirm that the subdivision schemes based on weighted-least squares with d = 0, 1 avoid Gibbs phenomenon, as stated in Corollary 4.9.To study it, we propose the following experiment.We discretize the function: x ∈ [0, 0.5]; − sin(πx), x ∈ (0.5, 1], in the interval [0, 1] with 33 equidistant points, x i = i • h, i = 0, . . ., 32 and h = 1 32 and apply the subdivision schemes.We show the results in Figure 5.It is clearly visualize that the Gibbs phenomenon does not appear around the discontinuity, but there is diffusion, instead.The larger is λ, the more diffusion, specially when rect is used.

Monotonicity
Finally, we introduce the last example in order to see numerically that the new family of the schemes conserves the monotonicity of the data, for d = 0, 1, proved in Corollary 4.8.We apply S 1,rect and S 1,trwt to the data collected in Table 8 (see [3]) and obtain Figure 6.

Conclusions and future work
In this work, a family of subdivision schemes based on weighted local polynomial regression has been analysed.We introduced the general form of this type of schemes and prove that the schemes corresponding to the polynomial degrees d = 2k and d = 2k + 1 coincide, for k = 0, 1, 2 . . .In particular, we analysed in detail the cases d = 0, 1, 2, 3 with positive weight functions, ω, with compact support.
In the first part of the paper, for d = 0, 1, we took advantage of the positivity of the mask to prove the convergence.Also, under some conditions of the ω functions, the C 1 regularity of the limit function was demonstrated.Afterward, some properties were proved as monotonicity and elimination of the Gibbs phenomenon effect.In the second part, we developed a general technique to analyse the convergence of a family of linear schemes and used it in the case d = 2, 3.
The last sections have been dedicated to discussing noise removal and approximation capabilities.We showed how the weight function φ determines these properties and that it is not possible to find a φ maximizing both capabilities approximation and noise reduction.This led to a multi-objective optimization problem in which optimal solutions were found along a Pareto front.Some numerical tests were presented to confirm the theoretical results.
For future works, we can consider the following ideas: The C 1 regularity of the cases d = 2, 3 were not proven.New theoretical tools such as those presented in Section 5 and their application to these schemes can be done.
We considered several weight functions φ from the literature.Now that we know the influence of φ in the approximation and denoising capabilities, it could be designed φ trying to improve them.Taking into account that the noise contribution is usually greater than the approximation error on the final curve, the use of an optimized weight function can be even more interesting than augmenting the polynomial degree, since some properties related to the monotonicity and the Gibbs phenomenon are only available for d = 0, 1.
If the data present some outliers, a different loss function can provide better results.Mustafa et al. in [23] proposed a variation of Dyn's schemes changing the ℓ 2 -norm by the ℓ 1 -norm in the polynomial regression but they do not prove their properties.The theoretical study of this scheme, as well as the use of different weight functions, can be considered in the future.

Figure 6 :
Figure 6: Limit curves for monotone data using subdivision schemes with rect (red line) and trwt (black line) weight functions, d = 0, 1.

Table 3 :
The function r(t) and the value R 1 of Theorem 5.1 for S 3,w λ several choices of φ and 2n − 1 < λn < 2n .