Iterative rational least squares fitting

1 Introduction

In the last decades, many papers have been devoted to Shepard-type operators, which are rational, positive operators interesting in classical approximation theory and in scattered data interpolation problems (see, e.g., [1, 3, 8, 10, 11, 12, 18]). Recently, in [4], Shepard-type operators have been studied to construct Shepard-type curves useful in CAGD. Such curves overcome some of the original Shepard operator’s drawbacks and have some advantages compared to the Bézier case. A weighted progressive iterative approximation (WPIA in short) technique for Shepard-type curves was also developed in [4]. The key idea was to iteratively change the control points of the active curve to deform towards the target shape represented by the points data. So, by adjusting the control points and by a proper weight, the WPIA process generates sequences of Shepard-type curves converging to the global Shepard-type interpolating curve at the original control points. Similar techniques for Bézier and B-spline curves were examined, for example, in [2, 9, 14, 15, 17]. However, in the WPIA method, the number of control points is equal to that of the data points, which is a drawback for a very large number of data points.

On the other hand, the least squares fitting (LSF in short) is one of the most commonly used mathematical tools in practice, particularly useful to fit a large number of data points by a polynomial or an approximation tool of lower degree. However, the corresponding system of linear equations can be ill-conditioned or even singular (see, e.g., [13, 16]).

The aim of the present paper is to overcome the above drawback of the WPIA technique by a progressive iterative approximation process for Shepard-type least squares fitting (SLSPIA in short) (Section 2). In each iteration, every control point is adjusted by a weighted sum of some difference vectors between the data points and their corresponding points on the fitting curve. By adjusting the control points iteratively, the SLSPIA technique generates a sequence of Shepard-type curves, and the limit curve is the Shepard-type least squares curve fitting the given data points. The convergence result for such a procedure, an explicit expression for the error and the corresponding decay rate are shown in Theorem 2.2. An optimal value of the weight parameter giving the fastest convergence rate is presented in Theorem 2.4. Moreover, a simple method to compute a practical value of the weight is discussed giving the convergence rate comparable to the optimal one. The SLSPIA technique is efficient since it avoids the computational cost of solving a large system of linear equations corresponding to the Shepard-type LSF method. We remark that such system is well-conditioned (Lemma 4.3). The above technique is interesting in CAGD because, by playing on a shape parameter – the iteration level –, the designer can draw several curves modeling the given data. Thanks to the above properties, our technique improves the algorithms in [4, 6, 9].

Then, in Section 3, starting from a rational operator of Shepard-type, we deduce a numerical procedure for the roots finding improving the chord method. Section 4 contains the proofs of the main results, which are based on some results on the eigenstructure of Shepard-type operators and direct estimates for our operators. Finally, numerical experiments confirming our theoretical results are shown in Section 5.

2 The progressive iterative approximation for Shepard-type least squares fitting

In this section, we introduce the SLSPIA technique. First we recall the definition of parametric Shepard-type curves [4]. Let An⁢(t)=[An,0⁢(t),An,1⁢(t),…,An,n⁢(t)]T, where

for 0≤i≤n, n∈N, t∈[0,1], ti=in, i=0,1,…,n, 0<ns⁢λ≤C with C any fixed positive constant and s>2 even.

Given the control vector P=[P0,P1,…,Pn]T, Pi∈Rd, i=0,1,…,n, d≥2, the parametric Shepard-type curve Sn⁢[P,t] is defined by

The properties of such curves interesting in CAGD were studied in [4]. In particular, we recall that Sn⁢[P] is a rational curve, preserving points, lying in the convex hull of the control polygon defined by P and satisfying the pseudo-local control property, i.e., each function An,j⁢(t), 0≤j≤n, attains its maximum value close to 1 at t=tj; in other words, the point Pj influences strongly the shape of the curve in a neighborhood of t=tj.

Generalizations of such curves were discussed in [6, 5, 7]. An efficient WPIA technique based on Sn curves for data fitting was also introduced and studied in [4]. However, this requires that the number of control points be equal to that of data points.

On the other hand, especially for points sets of large size, LSF is one the most commonly used mathematical tools. Therefore, a natural answer to the above restriction of the WPIA method for a large number of data is the construction of a weighted progressive iterative approximation technique for Shepard-type least squares fitting (SLSPIA) curves. Similarly to the classical WPIA, the SLSPIA process starts with an initial Shepard-type curve and constructs a sequence of fitting curves by adjusting the control points iteratively. The limiting curve is the Shepard-type least squares fitting curve to the data set. More precisely, assume that {Qi}i=0m is an ordered points set to be fitted. At the beginning of the iteration, we select {Pi}i=0n, n≤m, from {Qi}i=0m as the control points set and construct a piece of the blending curve 𝒫(0)⁢(t), i.e.,

with Pi(0)=Pi, i=0,1,…,n. The choice of the initial points will be discussed later. Letting

we construct the first adjusting vector for the i-th control point as

with μ a constant satisfying the condition

Here λn is the largest eigenvalue of AT⁢A, AT being the transpose of A, with A the collocation matrix of the blending basis An,i on {0=t¯0,t¯1,…,t¯m=1} defined by (cf. [4])

with t¯i=i/m, i=0,…,m (cf. Lemma 4.2). Then we can get the new control points by

and the new curve

Similarly, if 𝒫(k)⁢(t) is the k-th curve after the k-th iteration, we put

and the (k+1)-th curve is defined as

In this way, we have a curve sequence {𝒫(k)⁢(t)}k=0∞.

Now we prove that the limit of 𝒫(k)⁢(t), for k tending to ∞, is just the Shepard-type least squares fitting to the data points {Qi}i=0m. Indeed, denoting by Bk the k-th iteration of any matrix B, we have the following theorem.

The SLSPIA process makes it possible to construct a sequence of control polygons converging to the control polygon of a Shepard-type LSF curve. Therefore, the parameter k can be used as a shape parameter in order to model different shapes, obtaining, as extreme cases, the Shepard-type curve at some starting points Pi(0), i=0,1,…,n, and the Shepard-type LSF curve fitting the original control points.

For brevity, only the SLSPIA method for Shepard-type curves is handled, and the results for LSPIA process in the Shepard-type tensor product surfaces case can be deduced similarly [4].

In equation (2.4), a weight is required to generate the adjusting vector. The following theorem shows how to select an appropriate weight to improve the convergence rate of the SLSPIA method.

Although μ=2/(λ0+λn) provides the fastest convergence rate in theory (Theorem 2.4), it needs a large amount of computation to calculate the greatest and smallest eigenvalues. To avoid the computation of the eigenvalues, we follow [13] and propose a simple method to determine the weight μ. Indeed, we set

with C:=maxici, i=0,1,…,n, ci being the sum of the i-th row elements of AT⁢A (see Section 5 for the discussion on μ¯).

3 An improved chord length method

In this section, we discuss a numerical method for solving a nonlinear equation by a Shepard-type operator.

Let f∈C2⁢([a,b]), with a,b∈R, such that

and

From this assumption, we deduce that there is a unique root x¯∈(a,b) for f, i.e., f⁢(x¯)=0. Under assumptions (3.1) and (3.2), only one of the following two cases may occur:

1. f⁢(a)⁢f′′⁢(a)>0,

2. f⁢(b)⁢f′′⁢(b)>0.

Then consider the Shepard-type operator defined by

with s odd. Obviously,

In other words, T⁢(f) can be considered as a rational approximation of f. Hence, letting a∼b if and only if |a/b|±1≤C, it is reasonable that T⁢(f)∼f in a suitable interval and that we can use the root of T⁢(f) (easily computable) to approximate the root of f.

Then, in case (a), consider the iterative procedure

with

with x0=b and s odd.

In case (b), consider the iterative procedure

with

with x0=a and s odd.

Geometrically, xn+1 is the abscissa of the intersection of the graph of T⁢(f) with the x-axis.

We remark that if s=1, then the iteration procedure (3.3)–(3.6) gives back the classical chord method. It is well known that the chord method is converging linearly to x¯ under the above assumptions, while if f∉C′⁢([a,b]) or f′⁢(x¯)=0, the convergence is not guaranteed or the order is no longer linear.

Now we prove that the procedure (3.3)–(3.6) overcomes such drawbacks.

4 Proofs of the main results

5 Implementation and examples

5.1 Numerical experiments on Section 3

Though the SLSPIA method can be started with arbitrary initial control points, an appropriate selection of initial control points makes the LSPIA convergent quickly. In an implementation, following [13], we select the initial control points {Pi}i=0n as

(5.1)P0=Q0,Pn=Qm,Pi=Qf⁢(i),i=1,2,…,n-1,f⁢(i)=[(n+1)⁢in].

In this subsection, we present an example to show the efficiency and validity of our SLSPIA-type techniques.

5.1.1 Example 1

Consider a helix of radius 5 given by (see, e.g., [4, 15])

(x⁢(t),y⁢(t),z⁢(t))=(5⁢cos⁡t,5⁢sin⁡t,t),t∈[0,6⁢π].

A sequence of 301 control points Pi, i=0,…,300, is sampled from the helix as

(x⁢(si),y⁢(si),z⁢(si)),si=π50⁢i,i=0,1,…,300.

The above point set is fitted using SLSPIA curves given by (2.4) and (2.5) with 31 initial control points chosen as in (5.1), where s=4, λ=1.08×10-6 (see (2.1)) and μ¯ is given by (2.8). The original point set and initial Shepard-type curves are shown in Figure 1; the Shepard-type curves after the first and second iteration steps are shown in Figures 2 and 3.

Figure 1

SLSPIA curve at iteration level 0.

Figure 2

SLSPIA curve at iteration level 1.

Figure 3

SLSPIA curve at iteration level 2.

Figure 4

Limit SLSPIA fitting curve.

Moreover, denoting by ∥⋅∥ the usual Euclidean norm, we consider the fitting errors

Ek=∑j=0m∥Qj-∑i=0nAn,i⁢(t)⁢Pi(k)∥2

(scaled so that E0=1); if |Ek+1-Ek|<10-4, we stop the iteration process and let the last Ek be E∞. The limit fitting curve is presented in Figure 4. From Figures 1–4, we can see that LSPIA constructs a series of curves approximating the given data points progressively.

We also plot Ek, k=0,1,…,20, in Figure 5. From Figure 5, we can see that all Ek decrease very fast in some starting steps and {Ek} is a monotone decreasing sequence.

Figure 5

Plot of Ek at different iterations.

Figure 6

Plot of Dk (dotted line, left axis) and dk (continuous line, right axis).

Furthermore, let Dk=∑i=0n∥Pi(k)-Pi(∞)∥ (scaled so that D0=1). We plot Dk (k=1,2,…) in Figure 6, when we choose the weight μ¯=2/C. We use the same stop criterion as above. Figure 6 also shows the difference dk between the above Dk and Dk computed when the weight μ=2/(λ0+λn) (see (2.7)). Note that the scale for Dk and dk is different. We can see that the convergence rates corresponding to the above two weights are similar, as announced in Section 2.

5.1.2 Example 2

Consider the curve (x⁢(t),y⁢(t),z⁢(t)), t∈[0,6⁢π], where

x⁢(t)=(10.5⁢2⁢π⁢(exp⁡(-(t-2⁢π)2)+exp⁡(-(t-4⁢π)2))),

y⁢(t)=(10.5⁢2⁢π⁢(exp⁡(-(t-2⁢π-1)2)+exp⁡(-(t-4⁢π)2-1))),

z⁢(t)=t.

Figure 7

SLSPIA fitting curve at iteration level 0.

Figure 8

SLSPIA fitting curve at iteration level 1.

Figure 9

SLSPIA fitting curve at iteration level 2.

Figure 10

Limit SLSPIA fitting curve.

Figure 11

Plot of Ek at different iterations.

A sequence of 301 control points Pi, i=0,…,300, is sampled from the curve as

(x⁢(si),y⁢(si),z⁢(si)),si=π50⁢i,i=0,1,…,300.

The above point set is fitted using SLSPIA curves given by (2.4) and (2.5) with 31 initial control points chosen as in (5.1), where s=4, λ=1.08×10-5 (see (2.1)) and μ¯ is given by (2.8). The original point set and initial Shepard-type curves are shown in Figure 7, the Shepard-type curves after the first and second iteration steps are shown in Figures 8 and 9, while the limiting fitting curve is shown in Figure 10. Figure 11 shows the fast monotone decrease of Ek analogously to the helix case.