- Email: [email protected]
Scattered noisy Hermite data fitting using an extension of the weighted least squares method
Computers and Mathematics with Applications 65 (2013) 1967–1977
Contents lists available at SciVerse ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
Scattered noisy Hermite data fitting using an extension of the weighted least squares method✩ Tianhe Zhou, Zhong Li ∗ Department of Mathematics, Zhejiang Sci-Tech University, Hangzhou, Zhejiang, 310018, China
article
info
Article history: Received 16 July 2012 Received in revised form 7 January 2013 Accepted 15 April 2013 Keywords: Bivariate spline Scattered Hermite data with noise Extension of the weighted least squares method
abstract Given a set of scattered Hermite data with noise, we use an extension of the weighted least squares method to find the solution based on the bivariate spline. This method can adjust some weights according to different noise sizes to get a better approximation. We show that our method produces a unique spline to fit the data. Also we give the error bound for the method. In addition, we give some probability analysis for our method. Finally, we present some numerical experiments to demonstrate the performance of our method. © 2013 Elsevier Ltd. All rights reserved.
1. Introduction The spline function is an excellent tool for numerical approximation which is made up of piecewise polynomial functions with certain smoothness over given polygonal domains. The problem of reconstruction or approximation of a curve or surface from a scattered data set is one commonly encountered in many numerical applications. There are several methods for dealing with this kind of problem. One way is using the Dm -splines method to solve the problem (see, e.g., [1–5]). The Dm -spline is a particular case of a more general family of splines and based on a criterion of minimization. In [1,2,4] the authors use smoothing splines by choosing the parameter ε to fit data sets. When the data are exact, the convergence of the approximation does not depend too strongly on the choice of ε . It should not tend quickly to infinity, and in particular it might remain bounded (cf. [3]). However, in some applications, the situation is complicated by the great number of data and also because the data are usually not exact. For that case, the error bound of the smoothing Dm -splines for the noisy data has been studied in [3]. In [5], the authors extended the convergence and error estimate results for Dm -splines for noisy data to the variational splines by introducing a new term on the minimized functional weighted by a parameter vector. Another way is using the Bernstein–Bézier splines method to solve the problem (see, e.g., [6–8]). The Bernstein–Bézier splines are different from Dm -splines. They are composed of polynomials and based on a certain triangulation. In [7], the authors use Bernstein–Bézier splines to fit data sets by minimizing the general energy functional when the data are exact. In [8], the authors use the discrete least squares method and the penalized least squares method to fit noisy data. However, in some situations, we can split the data into two groups according to the size of the noise. One is the reliable data group— which means that the noise of the data is equal to zero or far less than 1. The other one is the unreliable data group. In [6], we present a method for fitting the Lagrange scattered data with two noise groups.
✩ This research was supported by the National Natural Science Foundation of China under Grants Nos 11201429 and 51075421 and Science and Technology Project of Zhejiang Province of China (No. 2012C21035). ∗ Corresponding author. Tel.: +86 571 86843540. E-mail address: [email protected] (Z. Li).
0898-1221/$ – see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.camwa.2013.04.011
1968
T. Zhou, Z. Li / Computers and Mathematics with Applications 65 (2013) 1967–1977
In our previous work [9], we introduced an extension of the weighted least squares method for fitting scattered Hermite data which cannot give an error analysis. We just presented some numerical examples there. In this paper, we mainly discuss the error bound for the extension of the weighted least squares method. In addition, we make some probability analysis, going deeper than in [9]. In many numerical applications, we have both Lagrange and Hermite scattered data. In meteorology, for example, the wind velocity value is the derivative of the wind potential function. Let V = {vi = (xi , yi )}Ni=1 be a set of µ,ν points lying in a domain Ω ⊂ R2 with polygonal boundary and {fi , 0 ≤ µ + ν ≤ r , i = 1, . . . , N } be a corresponding set of real numbers where µ, ν and r are integer numbers. Suppose the data are contaminated by noise fi
ν,µ
ν,µ
= Dνx Dµy f (vi ) + ϵi
ν,µ
r where f belongs to the standard Sobolev space W∞ (Ω ) and the ϵi s ∈ Sdr (△) which minimizes
L(s) :=
N
µ,ν
ωi (|△|µ+ν |Dµx Dνy s(vi ) − fi
µ,ν
are noisy terms. We propose to construct a surface
|)2
(1.1)
i=1 0≤µ+ν≤r
µ,ν
where Sdr (△) is a spline space defined in the next section and 0 ≤ ωi ≤ 1 are weight terms. We call this the extension of the weighted least squares method. We use the following method to determine the weight terms: when the data are reliµ,ν µ,ν µ,ν are not reliable, we set the corresponding able (which means |ϵi | = 0 or ≪ 1), we set ωi = 1, and when the data fi µ,ν µ,ν µ,ν weight terms ωi < 1 (e.g. ωi = 0.1, 0.01, . . . depending on the size of |ϵi |). In some circumstances, we only have the µ,ν Lagrange scattered data values; we set weight terms {ωi = 0, µ + ν ≥ 1}. That is the method of [6]. We can see that our method is an extension of the spline data fitting method of [6]. Next we derive the error bound for the method
∥f − s∥∞,Ω ≤ K (1 + C5 )|△|r |f |r ,∞,Ω where |f |r ,∞,Ω denotes the maximum norm of the rth derivatives of f over Ω and ∥f ∥∞,Ω is the standard infinity norm. Here |△| is the maximum of the diameters of the triangles in △, and K and C5 are constants in the error bound which will be introduced in following sections. Finally we give some probability analysis—which constitutes the biggest difference from [6]—and numerical results for our method. The paper is organized as follows. In Section 2 we review some well-known Bernstein–Bézier notation. The extension of the weighted least squares method is introduced in Section 3 together with a discussion on the solution. In Section 4 we prove that the error bound can be improved by averaging the coefficients of several splines which are constructed by fitting different sets of data. The numerical results for our method can be found in the last section. 2. Preliminaries Given a triangulation △ and integers 0 ≤ r < d, we write Sdr (△) := {s ∈ C r (Ω ) : s|T ∈ Pd , for all T ∈ △} for the usual space of splines of degree d and smoothness r, where Pd is the
d+2 2
-dimensional space of bivariate polyno-
mials of degree d. Throughout the paper we shall make extensive use of the well-known Bernstein–Bézier representation of splines. For each triangle T = ⟨v1 , v2 , v3 ⟩ in △ with vertices v1 , v2 , v3 , the corresponding polynomial piece s|T is written in the form s|T =
T T cijk Bijk ,
i+j+k=d
where the BTijk are the Bernstein–Bézier polynomials of degree d associated with T . In particular, if (λ1 , λ2 , λ3 ) are the barycentric coordinates of any point u ∈ R2 relative to the triangle T , then BTijk (u) :=
d! i!j!k!
λi1 λj2 λk3 ,
i + j + k = d·
T T As usual, we associate the Bernstein–Bézier coefficients {cijk }i+j+k=d with the domain points {ξijk := (iv1 +jv2 +kv3 )/d}i+j+k=d .
Definition 1. Let β < ∞. A triangulation △ is said to be β -quasi-uniform provided that |△| ≤ βρ△ , where ρ△ is the minimum of the radii of the incircles of the triangles of △. We say the data locations are evenly distributed over each triangle of △ if there are enough data points for every triangle such that
1/2
∥P ∥∞,T ≤ C
vi ∈T
P (vi )
2
T. Zhou, Z. Li / Computers and Mathematics with Applications 65 (2013) 1967–1977
1969
for any polynomial P and constant C . Given a vertex v , let star0 (v) = v , and starq (v) :=
{T ∈ △ : T ∩ starq−1 (v) = ∅},
q ≥ 1.
Recall that a determining set for a spline space S ⊆ Sd0 (△) is a subset M of the set of domain points such that if s ∈ S and cξ = 0 for all ξ ∈ M , then cξ = 0 for all domain points, i.e., s ≡ 0. The set M is called a minimal determining set (MDS) for S if there is no smaller determining set. It is known that M is a MDS for S if and only if every spline s ∈ S is uniquely determined by its set of B-coefficients {cξ }ξ ∈M . Next we introduce the Markov inequality. Theorem 2.1. Let T := ⟨v1 , v2 , v3 ⟩ be a triangle and fix 1 ≤ q ≤ ∞. Then there exists a constant C depending only on d and q such that, for any polynomial p and any nonnegative integers α and β with 0 ≤ α + β ≤ d,
∥Dαx Dβy p∥q,T ≤
C α+β
ρT
∥p∥q,T
where ρT denotes the radius of largest circle which can be inscribed in T . Recall from [10] that for any given function f ∈ L1 (Ω ), there exists a quasi-interpolatory operator Q mapping f ∈ L1 (Ω ) to Sdr (△) which achieves the optimal approximation order of Sdr (△). That is: Theorem 2.2 (Cf. [10]). Fix d ≥ 3r + 2 and 0 ≤ m ≤ d. Suppose f lies in the Sobolev space Wpm+1 (Ω ) with 1 ≤ p ≤ ∞. Then
∥Dαx Dβy (f − Qf )∥p,Ω ≤ K |△|m+1−α−β |f |m+1,p,Ω for 0 ≤ α + β ≤ m, where |△| is the mesh size of △ (i.e., the diameter of the largest triangle), and |f |m+1,p,Ω is the usual Sobolev semi-norm. If Ω is convex, then the constant K depends only on d, p, m, and on the smallest angle in △. If Ω is nonconvex, it also depends on the Lipschitz constant L∂ Ω associated with the boundary of Ω . 3. Discussion of a solution Let S = Sdr (△) and M = dim S with d ≥ 3r + 2. Clearly, we can see that M ≪ N, where N denotes the number of given data values. Note that using the locally supported basis functions (cf. [11]), any spline function s in the space can be M M represented by s = i=1 ci φi for the set of coefficients {ci }i=1 . The extension of the weighted least squares method is finding Sω ∈ S such that L(Sω ) = min{L(s) : s ∈ S }
(3.1)
where L(s) is defined by (1.1). First, we discuss the existence and uniqueness of the solution Sω ∈ S satisfying (3.1) and an error bound as well. µ,ν
Theorem 3.1. Suppose all the data locations are evenly distributed over each triangle of △. For fixed weight terms {ωi µ + ν ≤ r }Ni=1 , there exists a unique Sω ∈ S satisfying (3.1).
,0 ≤
Proof. The proof of the existence uses an argument similar to that in the proof of Theorem 2 in [6]. We omit the details here. T Next we show the uniqueness of the minimizer Sω . Let c = {cijk , i + j + k = d, T ∈ △} be the coefficient vector associated with s. Then N
L(s) =
µ,ν
ωi (|△|µ+ν |Dµx Dνy s(vi ) − fi
µ,ν
|)2
i=1 0≤µ+ν≤r
=
T ∈△ (xℓ ,yℓ )∈T 0≤µ+ν≤r
µ,ν ωℓ
|△|
µ+ν
2 µ,ν T µ ν T cijk Dx Dy Bijk (xℓ , yℓ ) − fℓ i+j+k=d
= cT Bc − 2bT c + ∥f∥22 µ,ν
µ,ν
1
where f = {(ωℓ ) 2 · |△|µ+ν fℓ
t
B =
0≤µ+ν≤r vℓ ∈t
µ,ν
, ℓ = 1, . . . , N , 0 ≤ µ + ν ≤ r } is a data value vector. Here B = diag(Bt , t ∈ △), with l+m+n=d
ωℓ |△|2µ+2ν Dµx Dνy Btijk (vℓ )Dµx Dνy Btlmn (vℓ ) i+j+k=d
being a matrix of size dˆ × dˆ with dˆ = btijk =
0≤µ+ν≤r vℓ ∈t
µ,ν
d+2 2
µ,ν
and b = (
btijk
ωℓ |△|2µ+2ν fℓ Dµx Dνy Btijk (vℓ ).
, i + j + k = d, t ∈ △) with
1970
T. Zhou, Z. Li / Computers and Mathematics with Applications 65 (2013) 1967–1977
Suppose that we have two solutions Sω and Sˆω . Let c and cˆ be the two coefficient vectors associated with Sω and Sˆω respectively. Since L is a convex functional, we have, for any z ∈ [0, 1], L(zSω + (1 − z )Sˆω ) ≤ zL(Sω ) + (1 − z )L(Sˆω ) = L(Sω )· That is, L(zSω +(1 − z )Sˆω ) is a constant function of z ∈ [0, 1]. It follows that 0 =
d dz
d L dz
(zSω +(1 − z )Sˆω ) = 0 for all z ∈ [0, 1]. That is,
L(zSω + (1 − z )Sˆω )
= 2z (c − cˆ )T B(c − cˆ ) − 2bT (c − cˆ ) for all z ∈ [0, 1]. It follows that
t ∈△ (xℓ ,yℓ )∈t i+j+k=d
µ,ν
t t 2 µ ν t |△|2µ+2ν ωℓ (cijk − cˆijk ) Dx Dy Bijk (vℓ )2 = 0
for any {µ, ν : 0 ≤ µ + ν ≤ r }. That is, c = cˆ , because when the data locations are evenly distributed over each triangle of △, B is invertible. Hence, the minimizer is unique. Next we derive the error bound for the solution. We convert the problem (3.1) into a standard approximation problem in Hilbert space. Let r X := {f ∈ B(Ω ) : f |T ∈ W∞ (T ), all triangles T in △}
where B(Ω ) is the set of all bounded real-valued functions on Ω . For each triangle T in △, let
⟨f , g ⟩XT :=
µ,ν
(xℓ ,yℓ )∈T 0≤µ+ν≤r
ωℓ |△|2µ+2ν Dµx Dνy f (vℓ )Dµx Dνy g (vℓ )·
Then
⟨f , g ⟩X :=
⟨f , g ⟩XT T ∈△
defines a semi-definite inner-product on X . Let ∥f ∥XT and ∥f ∥X be the associated semi-norms. Given f ∈ X , supposing that S ⊆ Sdr (△) is a Hilbert space with respect to ∥ · ∥X , we need to find Sω ∈ S such that L(Sω ) = min L(s) s∈S
where L(s) := ∥f − s∥2X . Define a projection operator P : X → S by Pf = Sω . It is easy to see that the approximation Sω of f is characterized by
⟨f − Sω , s⟩ = 0,
∀s ∈ S .
(3.2)
Moreover, using the Cauchy–Schwarz inequality, we can get
∥Pf ∥X ≤ ∥f ∥X (3.3) for all f ∈ X . More details can be found in [11,12]. Recall from [10] that when d ≥ 3r + 2, Sdr (△) possesses a stable local basis {Bξ }ξ ∈M corresponding to a minimal determining set M . Thus, for any ξ ∈ M , there is a vertex vξ of △ with supp(Bξ ) ⊆ starl (vξ ) and ∥Bξ ∥∞,Ω = 1 for all ξ ∈ M , where l is an integer number related to ξ (see Remark 6.1). Also, there exists a positive constant K1 depending only on d and β such that for any s = ξ ∈M cξ Bξ and any ξ ∈ M, K1 |cξ | ≤ ∥s∥∞,Tξ for some triangle Tξ in △ with Tξ ⊆ supp(Bξ ). Letting µ,ν
µ,ν
ωmin := min{ωi |i = 1, . . . , N , 0 ≤ µ + ν ≤ r , ωi ̸= 0}, µ,ν ω := max ωT |T ∈ △, ωT = ωi , (xi ,yi )∈T 0≤µ+ν≤r
we have Lemma 1. Let {Bξ }ξ ∈M be the basis for S corresponding to the minimal determining set M . Suppose that the data set has a property that for every s ∈ S and every T ∈ △
1/2
K2 ∥s∥∞,T ≤
|△|
2µ+2ν
|△|
2µ+2ν
µ ν
2
µ ν
2
Dx Dy s(vℓ )
,
(3.4)
vℓ ∈T 0≤µ+ν≤r
1/2
K3 ∥s∥∞,T ≥
vℓ ∈T 0≤µ+ν≤r
Dx Dy s(vℓ )
(3.5)
T. Zhou, Z. Li / Computers and Mathematics with Applications 65 (2013) 1967–1977
1971
for some positive constants K2 , K3 . Then there exist positive constants C1 and C2 depending on d, l and β such that
2 ωmin C1 cξ ≤ cξ Bξ ≤ ωC2 cξ2 ξ ∈M ξ ∈M ξ ∈M
2
(3.6)
X
for all {cξ }ξ ∈ M .
Proof. Let s = K12
ξ ∈M
ξ ∈M cξ Bξ .
cξ2 ≤
ξ ∈M
By the definition of the local stable basis we have
∥s∥2∞,Tξ .
Next we use Eq. (3.4) to get
∥s∥∞,T ≤
1/2
1
K2
1/2
1
ω
µ ν
Dx Dy s(vℓ )
ωmin K2
=
|△|
2µ+2ν
2
vℓ ∈T 0≤µ+ν≤r
1
≤
1/2 min K2
|△|
2µ+2ν
vℓ ∈T 0≤µ+ν≤r
µ,ν ωℓ |Dµx Dνy s(vℓ )|2
1/2
∥s∥XT ·
Then we have
ξ ∈M
cξ2 ≤
1
ξ ∈M
ωmin K12 K22
∥s∥2XT . ξ
Note that there may be more than one domain point η ∈ M in Tξ . Let K4 = maxT ∈△ ♯{ξ ∈ M : T ∩ supp(Bξ ) ̸= ∅}, where ♯ denotes the cardinality of set; then
ξ ∈M
cξ2 ≤
K4
ωmin K12 K22
∥s∥2X .
Therefore we have the left hand side of (3.6) with C1 =
∥s∥XT =
(xℓ ,yℓ )∈T 0≤µ+ν≤r
µ,ν
K12 K22 K4
. By using (3.5), we have
µ
ωℓ |△|2µ+2ν |Dx Dνy s(vℓ )|2
≤ ω1/2 K3 ∥s∥∞,T . Let MT := {ξ ∈ M : T ∩ supp(Bξ ) ̸= ∅}; then ♯MT ≤ K4 . Noting that ∥Bξ ∥∞,Ω = 1 for all ξ ∈ M , we have
2 cξ Bξ ∥s ∥ ≤ ∥s ∥ = T ∈△ T ∈△ ξ ∈MT XT 2 ≤ ωK32 cξ Bξ T ∈△ ξ ∈MT ∞,T 2 ≤ ωK32 |cξ | 2 X
2 XT
T ∈△
≤ω
K32
ξ ∈MT
♯MT
≤ K4 ω
|cξ |2
ξ ∈MT
T ∈△
K32
|cξ |2 .
T ∈△ ξ ∈MT
Next let K5 := maxξ ∈M ♯{T ∈ △ : T ∩ supp(Bξ ) ̸= ∅}; then
∥s∥2X ≤ ωK32 K4 K5
|cξ |2 .
ξ ∈M
Finally the proof is complete with C2 = K32 K4 K5 .
1972
T. Zhou, Z. Li / Computers and Mathematics with Applications 65 (2013) 1967–1977
Given a triangle T , let star0 (T ) = T , and starq (T ) :=
{starq (w) : w is a vertex of T },
q ≥ 1. ωC
Lemma 2. Suppose S satisfies the hypotheses of Lemma 1; let C3 = γ (γ + 1)1/2 with γ = ω 2C . Then there exists a constant min 1 0 ≤ σ ≤ 1 such that, for any triangle T ∈ △ and any function g ∈ X with supp(g ) ⊆ T ,
∥Pg ∥XT ≤ C3 σ k ∥g ∥X
(3.7)
whenever T ∈ star2(k+2)l+1 (T ) \ star2(k+1)l+1 (T ) with k ≥ 1. Proof. The proof is similar to that of Theorem 3.1 of [12]. So we put it into Remark 6.3.
Note that for any g ∈ X , we have
∥g ∥XT =
µ,ν
(xℓ ,yℓ )∈T 0≤µ+ν≤r
µ
ωℓ |△|2µ+2ν |Dx Dνy g (vℓ )|2
≤ω
1/2
∥g ∥∞,T +
|△|
µ+ν
µ ν
∥Dx Dy g ∥∞,T .
1≤µ+ν≤r
1/2 n0 + C3 k≥1 σ k nk ; then Lemma 3. Suppose that all the hypotheses of Lemma 2 hold. Let C4 = ω1/2 ωmin K2
|△|µ+ν ∥Dµx Dνy g ∥∞,Ω
∥Pg ∥∞,Ω ≤ C4 ∥g ∥∞,Ω +
(3.8)
1≤µ+ν≤r
where nk := max[star2(k+2)l+1 (T ) \ star2(k+1)l+1 (T )]
n0 = max ♯star4l+1 (T ), T ∈△
T ∈△
and C3 , σ are the constants in (3.7). Proof. Let T be any triangle in △ and
ΩkT := star2(k+2)l+1 (T ) \ star2(k+1)l+1 (T ).
Ω0T := star4l+1 (T ), Note that
∥s∥∞,T ≤
1 1/2 min K2
ω
∥s∥XT .
We can get
∥Pg ∥∞,T ≤
1 1/2
ωmin K2
n0 + C 3
σ nk ∥g ∥XT k
k≥1
as in Theorem 4.1 of [12]. Then
n0 + C 3
∥Pg ∥∞,T ≤
k≥1
1/2
ωmin K2
σ nk k
ω
1/2
∥g ∥∞,T +
|△|
µ+ν
µ ν
∥Dx Dy g ∥∞,T .
1≤µ+ν≤r
1/2 Taking the supremum over all T ∈ △, we get (3.8) with C4 = ω1/2 n0 + C3 k≥1 σ k nk . ωmin K2
Finally we have our main theorem in this paper. Theorem 3.2. Suppose △ is a β -quasi-uniform triangulation and let Sω be the spline minimizing L defined above. Let K be the (r +1)(r +2) constant in Theorem 2.2 and C5 = C4 ; we have 2
∥f − Sω ∥∞,Ω ≤ K (1 + C5 )|△|r |f |r ,∞,Ω r for every function f ∈ W∞ (Ω ).
(3.9)
T. Zhou, Z. Li / Computers and Mathematics with Applications 65 (2013) 1967–1977
1973
Proof. Let sf = Qf be the quasi-interpolation function introduced in Section 2. By the definition of Qf , we can see that Psf = sf and therefore
∥f − Sω ∥∞,Ω = ∥f − Pf ∥∞,Ω ≤ ∥f − sf ∥∞,Ω + ∥sf − Pf ∥∞,Ω ≤ ∥f − sf ∥∞,Ω + ∥Psf − Pf ∥∞,Ω . Next, by Lemma 3 and Theorem 2.2,
∥f − Sω ∥∞,Ω ≤ (1 + C4 )∥f − sf ∥∞,Ω + C4
|△|µ+ν ∥Dµx Dνy (f − sf )∥∞,Ω
1≤µ+ν≤r
≤ K (1 +
(r + 1)(r + 2)
C4 )|△|r |f |r ,∞,Ω 2 = K (1 + C5 )|△|r |f |r ,∞,Ω . Therefore the proof is complete.
µ,ν
From the proof of Theorem 3.2, we can see that if we set all ωi = 0 for µ + ν ≥ 1, we can get the error bound of the µ,ν µ,ν weighted least squares method in [6], and if we set all ωi = 1 for µ + ν = 0 and ωi = 0 for µ + ν ≥ 1, we can get the error bound of the traditional least squares method. 4. Probability analysis In this section, we will prove that the error bound can be improved by averaging the coefficients of several splines which are constructed by fitting different sets of data. µ,ν Let {{fi }j , 0 ≤ µ + ν ≤ 1, i = 1, . . . , N }M j=1 be the M sets of the scattered data with noise. Suppose each of the sets has sufficient number of scattered data and can be used to construct the spline. We use Pf j to denote the spline which is constructed from the jth data set. It is easy to see that Pf 1 , Pf 2 , . . . , Pf M is a sequence of independent random variables; they each have the same distribution with mean µ ¯ and variance σ 2 . Here we suppose that the noisy term has the uniform distribution on the interval [−1, 1]; then it is easy to see that µ ¯ = f and σ 2 = 1/3 (see Remark 6.2). According to the Central Limit Theorem (cf.[13]), we have the distribution of Pf 1 + Pf 2 + · · · + Pf M − Mf
√ σ M
tending to the Standard Normal Distribution as M −→ ∞. That is,
P
Pf 1 + Pf 2 + · · · + Pf M − Mf
√
σ M
≤ c −→ Φ (c ),
as M −→ ∞, where Φ denotes the cumulative area under the Standard Normal Distribution. We use s :=
M 1
M j =1
Pf j
to approximate the original function f . Then given a constant δ > 0, we have
Pf 1 + Pf 2 + · · · + Pf M − Mf P {|s − f | < δ} = P √ σ M
√ √ δ M δ M ∼ ≤ − 1. = 2Φ σ σ
Next we choose δ = 0.1 and δ = 0.01 to check the dependence on the number M and the probability P {|s − f | < δ} (see Tables 1 and 2). In numerical experiments, we find that when M > 12, the distribution tends to the Standard Normal Distribution as well. 5. Numerical experiments The purpose of this section is to demonstrate the performance of our method. The numerical results show that our method is much better than the traditional least squares method. Example 5.1. Suppose △ is a triangulation of the unit square domain Ω := [0, 1] × [0, 1]. Consider 2000 random points (xi , yi ) over △ as shown in Fig. 1. We use the following test functions: f1 (x, y) = x2 + 2y2 , f2 (x, y) = sin(2(x − y)), f3 (x, y) = 0.75 exp(−0.25(9x − 2)2 − 0.25(9y − 2)2 ) + 0.75 exp(−(9x + 1)2 /49 − (9y + 1)/10)
+ 0.5 exp(−0.25(9x − 7)2 − 0.25(9y − 3)2 ) − 0.2 exp(−(9x − 4)2 − (9y − 7)2 ),
1974
T. Zhou, Z. Li / Computers and Mathematics with Applications 65 (2013) 1967–1977 Table 1 The relationship between M and P {|s − f | < δ} with δ = 0.1. The probability P {|s − f | < δ}
The number M
0.99 0.95 0.90 0.85 0.80
225 128 90 69 55
Table 2 The relationship between M and P {|s − f | < δ} with δ = 0.01. The probability P {|s − f | < δ}
The number M
0.99 0.95 0.90 0.85 0.80
22 533 12 805 8 965 6 912 5 461
Table 3 The approximation errors.
f1 (x, y) f2 (x, y) f3 (x, y)
Our method
Traditional least squares method
0.1882 0.0409 0.0180
0.8536 0.6460 0.3017
Fig. 1. The scattered data and the triangulation △.
where f3 is the well-known Franke function. We use these function values and derivative values at the points (xi , yi ), i = µ µ,ν 0, . . . , 2000, to have a set of scattered Hermite data, that is, {(xi , yi , Dx Dνy f (xi , yi ) + ϵi ), 0 ≤ µ + ν ≤ 1, i = 1, . . . , 2000}. µ,ν
µ,ν
We set ϵi = 0 when the corresponding point (xi , yi ) is marked with a cross (which indicates reliable data), and set ϵi to be a random number between −1 to 1 when the corresponding point (xi , yi ) is marked with a dot (which means unreliable data). We approximate these data in spline spaces of S82 (△) using two methods: the traditional least squares method and our proposed extension of the weighted least squares method where the weights are set to be ωi = 1 and ωi = 0.01 for the reliable data and non-reliable data, respectively. Table 3 gives the maximum errors for the two methods when the computation is based on 100 × 100 equally spaced points over Ω . Discussion: From Table 3, we can see that the surface created by the traditional least squares method does not perform well on the test function. However, our method produces much better approximate solutions.
Example 5.2. The spline space, scattered data set and test functions are the same as in Example 5.1. Here, we approximate the sample functions with different weight terms. As in Example 5.1, we set the weights ωi = 1 for reliable data. Furthermore, we choose different weights ωi for unreliable data to check the dependence of the performance of our method on the weights. Some numerical results are given in Table 4.
T. Zhou, Z. Li / Computers and Mathematics with Applications 65 (2013) 1967–1977
1975
Fig. 2. Triangulations for choice N = 2, 32. Table 4 The approximation errors with different weight terms.
f1 (x, y) f2 (x, y) f3 (x, y)
ω = 0.01
ω = 0.001
ω = 0.0001
ω = 0.00001
0.1882 0.0409 0.0180
0.0414 0.0065 0.0173
0.0110 0.0020 0.0167
0.0109 6.51 × 10−4 0.0151
Table 5 Maximum errors for various test functions. N
2
4
8
16
32
f1 f2 f3
0.25221 0.06892 0.04521
0.09072 0.01793 0.01592
0.02974 0.00511 0.00532
0.01114 0.00119 0.00136
3.43 × 10−3 2.83 × 10−4 4.41 × 10−4
Discussion: From Table 4, we see that the error decreases when the weight terms decrease. For the error in approximating f3 (x, y), we observe that it does not decrease quickly when the weight terms decrease. This agrees with the prediction of the error bound (3.9) where the constant C = K (1 + C5 ) has a term which is independent of ωi and at the same time ωmin is too small. Example 5.3. Suppose ♦ is a uniform partition of the unit square domain Ω := [0, 1] × [0, 1] into N 2 subsquares. Let S82 (△♦ ) be a C 2 spline space, where △♦ is the triangulation obtained by inserting the diagonals of each of the subsquares in ♦ (see Fig. 2). We use the same test functions as in Example 5.1. Suppose that there are enough scattered Hermite data in each triangle; this means that the data locations are evenly distributed over each triangle of △♦ , and half of the data are reliable. So we set ωi = 1 and ωi = 0.01 for the reliable data and non-reliable data, respectively. We approximate these functions for choices N = 2, 4, 8, 16, 32, which correspond to repeatedly halving the mesh size. Table 5 gives the maximum error when the computation is based on 100 × 100 equally spaced points over Ω . Table 6 presents the corresponding ratios of errors for successive values of N. Discussion: From Table 6, we can see that the approximation order is close to 3, which is what we expect for order 2 convergence. Example 5.4. In this example, suppose that we have a sufficient number of scattered data. We divide these data into M sets M j and use s := M1 j=1 Pf to approximate the test functions. Here, we use the same test functions as in Example 5.1 to check the dependence of the performance of our method on M. Unlike for Example 5.1, we set all the weight terms as ω = 1 no matter whether the data are reliable or unreliable. Some numerical results are given in Table 7. Discussion: From Table 7, we can see that the error decreases when M increases. That agrees with the Central Limit Theorem. But the error in approximating f2 (x, y) by using M = 250 is bigger than that in using M = 100. This is because the formula
√
P {|s − f | < δ} ∼ = 2Φ δ
3M − 1
√
concerns probability. This means that it has 2Φ δ
3M − 1 probability of happening.
1976
T. Zhou, Z. Li / Computers and Mathematics with Applications 65 (2013) 1967–1977 Table 6 Ratios of convergence for various test functions. Ratio
2/4
4/8
8/16
16/32
f1 f2 f3
2.78 3.85 2.84
3.05 3.51 2.99
2.67 4.28 3.91
3.25 4.21 3.09
Table 7 The approximation errors with different values of M.
f 1 (x , y ) f 2 (x , y ) f 3 (x , y )
M = 20
M = 50
M = 100
M = 250
M = 500
0.0871 0.0904 0.1571
0.0452 0.0521 0.1014
0.0217 0.0212 0.0983
0.0084 0.0109 0.1006
0.0065 0.0081 0.0891
6. Remarks Remark 6.1. It is proved that there exists an integer l such that for each ξ , the support of Bξ is contained in starl (vξ ) for some vertex vξ ∈ △ in [10]. Remark 6.2. In this paper, we assume that the noisy term has a uniform distribution. So we can perform a probability analysis. Of course, we can use the same method to analyze the probability if the noisy term has another distribution with different mean and variance. Remark 6.3. In order to make the paper self-contained, We add the proof of Lemma 2 here. Proof. Let
M0T := {ξ ∈ M : supp(Bξ ) ∩ T ̸= ∅}, MqT := {ξ ∈ M : supp(Bξ ) ∩ starq(2l+1) T ̸= ∅}, N0T := M0T , NqT := MqT \ MqT−1 . Suppose that Pg is given by
Pg =
cξ Bξ ,
ξ ∈M
and let
uq :=
cξ Bξ ,
wq := Pg − uq ,
aq =
ξ ∈MqT
cξ2 ,
ξ ∈NqT
for any q ≥ 0. By (3.6), we know that
aj =
j≥q+1
ξ ̸∈MqT
cξ2 ≤
∥wq ∥2X , ωmin C1
q ≥ 0.
Now, ⟨g − Pg , wq ⟩ = 0 from (3.2), and ⟨g , wq ⟩ = 0 since wq ≡ 0 in T and g ≡ 0 outside of T . Using the fact that supp(wq ) ∩ ξ ∈MT supp(Bξ ) = ∅ for q ≥ 1, it follows that q−1
∥w ∥ = ⟨Pg − uq , wq ⟩ = ⟨g − uq , wq ⟩ = −⟨uq , wq ⟩ ∥wq ∥X , =− cξ Bξ , wq ≤ c B ξ ξ T T ξ ∈Nq ξ ∈Nq 2 q X
X
and thus by (3.6),
2 ∥wq ∥2X ≤ c B ξ ξ ≤ ω C 2 aq . ξ ∈NqT X
T. Zhou, Z. Li / Computers and Mathematics with Applications 65 (2013) 1967–1977
1977
Combining the above, we have
aj ≤
j ≥ q +1
ω C2 aq , ωmin C1
q ≥ 0.
(6.1) ωC
Then applying Lemma 2 in [14] with γ := ω 2C , we see that min 1 aq ≤ (γ + 1)ρ q a0 , where ρ := γ /(γ + 1). Since ∥Pg ∥2 ≤ ∥g ∥2 , we have
a0 ≤
j ≥0
aj =
ξ ∈M
cξ2 ≤
1
ωmin C1
∥Pg ∥2X ≤
1
ωmin C1
∥g ∥2X .
Let τ be a triangle of △ which lies outside of stark (T ). If 1 ≤ k ≤ 2l, then
∥Pg ∥Xτ ≤ ∥Pg ∥X ≤ ∥g ∥X , and (3.7) holds for any 0 ≤ σ ≤ 1 with C3 = 1. Suppose now that k ≥ 2l + 1, and let q := implies ξ ̸∈ MqT , and thus by (3.6),
∥Pg ∥
2 Xτ
k 2l+1
. Then supp(Bξ ) ∩ τ ̸= ∅
2 ω2 C22 c B cξ2 = ωC2 aj ≤ ≤ aq . ξ ξ ≤ ω C2 ωmin C1 ξ ̸∈MqT T j ≥ q + 1 ξ ̸∈Mq X
1
Combining the above estimates yields (3.7) with σ := ρ 4l+2 .
References [1] R. Arcangéli, M.C. Silanes, J.J. Torrens, Multidimensional Minimizing Splines, Kluwer Academic Publishers, 2004. [2] M.C. López Silanes, R. Arcangéli, Sur la convergence des Dm -splines d’ajustement pour des donnĺe¸es exactes ou bruitĺe¸es, Rev. Mat. Univ. Comp. Mad. 4 (2–3) (1991) 279–284. [3] D.L. Ragozin, Error bounds for derivative estimates bases on spline smoothing of exact or noisy data, J. Approx. Theory 37 (1983) 335–355. [4] R. Arcangéli, B. Ycart, Almost sure convergence of smoothing Dm -spline functions, Numer. Math. 66 (1993) 281–294. [5] A. Kouibia, M. Pasadas, Approximation by smoothing variational vector splines for noisy data, J. Comput. Appl. Math. 211 (2008) 213–222. [6] T. Zhou, D. Han, A weighted least squares method for scattered data fitting, J. Comput. Appl. Math. 217 (2008) 56–63. [7] T. Zhou, D. Han, M.J. Lai, Energy minimization method for scattered data Hermite interpolation, Appl. Numer. Math. 58 (2008) 646–659. [8] G. Awanou, M.J. Lai, P. Wenston, The multivariate spline method for scattered data fitting and numerical solutions of partial differential equations, in: G. Chen, M.J. Lai (Eds.), Wavelets and Splines, Nashboro Press, Brentwood, Tennessee, 2006, pp. 24–76. [9] T. Zhou, Z. Li, Scattered noisy data fitting using bivariate splines, Proc. Eng. 15 (2011) 1942–1946. [10] M.J. Lai, L.L. Schumaker, On the approximation power of bivariate splines, Adv. Comput. Math. 9 (1998) 251–279. [11] O. Davydov, L.L. Schumaker, On stable local bases for bivariate polynomial spline spaces, Constr. Approx. 18 (2001) 87–116. [12] M.V. Golitschek, L.L. Schumaker, Bounds on projections onto bivariate polynomial spline spaces with stable local bases, Constr. Approx 18 (2002) 241–254. [13] K. Subrahmaniam, A Primer in Probability, second ed., Mercel Dekker, Inc., 1990. Revised and Expanded. [14] C.D. Boor, A bound on the L∞ norm of L2 approximation by splines in terms of a global mesh ratio, Math. Comp. 30 (1976) 765–771.
Contents lists available at SciVerse ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
Scattered noisy Hermite data fitting using an extension of the weighted least squares method✩ Tianhe Zhou, Zhong Li ∗ Department of Mathematics, Zhejiang Sci-Tech University, Hangzhou, Zhejiang, 310018, China
article
info
Article history: Received 16 July 2012 Received in revised form 7 January 2013 Accepted 15 April 2013 Keywords: Bivariate spline Scattered Hermite data with noise Extension of the weighted least squares method
abstract Given a set of scattered Hermite data with noise, we use an extension of the weighted least squares method to find the solution based on the bivariate spline. This method can adjust some weights according to different noise sizes to get a better approximation. We show that our method produces a unique spline to fit the data. Also we give the error bound for the method. In addition, we give some probability analysis for our method. Finally, we present some numerical experiments to demonstrate the performance of our method. © 2013 Elsevier Ltd. All rights reserved.
1. Introduction The spline function is an excellent tool for numerical approximation which is made up of piecewise polynomial functions with certain smoothness over given polygonal domains. The problem of reconstruction or approximation of a curve or surface from a scattered data set is one commonly encountered in many numerical applications. There are several methods for dealing with this kind of problem. One way is using the Dm -splines method to solve the problem (see, e.g., [1–5]). The Dm -spline is a particular case of a more general family of splines and based on a criterion of minimization. In [1,2,4] the authors use smoothing splines by choosing the parameter ε to fit data sets. When the data are exact, the convergence of the approximation does not depend too strongly on the choice of ε . It should not tend quickly to infinity, and in particular it might remain bounded (cf. [3]). However, in some applications, the situation is complicated by the great number of data and also because the data are usually not exact. For that case, the error bound of the smoothing Dm -splines for the noisy data has been studied in [3]. In [5], the authors extended the convergence and error estimate results for Dm -splines for noisy data to the variational splines by introducing a new term on the minimized functional weighted by a parameter vector. Another way is using the Bernstein–Bézier splines method to solve the problem (see, e.g., [6–8]). The Bernstein–Bézier splines are different from Dm -splines. They are composed of polynomials and based on a certain triangulation. In [7], the authors use Bernstein–Bézier splines to fit data sets by minimizing the general energy functional when the data are exact. In [8], the authors use the discrete least squares method and the penalized least squares method to fit noisy data. However, in some situations, we can split the data into two groups according to the size of the noise. One is the reliable data group— which means that the noise of the data is equal to zero or far less than 1. The other one is the unreliable data group. In [6], we present a method for fitting the Lagrange scattered data with two noise groups.
✩ This research was supported by the National Natural Science Foundation of China under Grants Nos 11201429 and 51075421 and Science and Technology Project of Zhejiang Province of China (No. 2012C21035). ∗ Corresponding author. Tel.: +86 571 86843540. E-mail address: [email protected] (Z. Li).
0898-1221/$ – see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.camwa.2013.04.011
1968
T. Zhou, Z. Li / Computers and Mathematics with Applications 65 (2013) 1967–1977
In our previous work [9], we introduced an extension of the weighted least squares method for fitting scattered Hermite data which cannot give an error analysis. We just presented some numerical examples there. In this paper, we mainly discuss the error bound for the extension of the weighted least squares method. In addition, we make some probability analysis, going deeper than in [9]. In many numerical applications, we have both Lagrange and Hermite scattered data. In meteorology, for example, the wind velocity value is the derivative of the wind potential function. Let V = {vi = (xi , yi )}Ni=1 be a set of µ,ν points lying in a domain Ω ⊂ R2 with polygonal boundary and {fi , 0 ≤ µ + ν ≤ r , i = 1, . . . , N } be a corresponding set of real numbers where µ, ν and r are integer numbers. Suppose the data are contaminated by noise fi
ν,µ
ν,µ
= Dνx Dµy f (vi ) + ϵi
ν,µ
r where f belongs to the standard Sobolev space W∞ (Ω ) and the ϵi s ∈ Sdr (△) which minimizes
L(s) :=
N
µ,ν
ωi (|△|µ+ν |Dµx Dνy s(vi ) − fi
µ,ν
are noisy terms. We propose to construct a surface
|)2
(1.1)
i=1 0≤µ+ν≤r
µ,ν
where Sdr (△) is a spline space defined in the next section and 0 ≤ ωi ≤ 1 are weight terms. We call this the extension of the weighted least squares method. We use the following method to determine the weight terms: when the data are reliµ,ν µ,ν µ,ν are not reliable, we set the corresponding able (which means |ϵi | = 0 or ≪ 1), we set ωi = 1, and when the data fi µ,ν µ,ν µ,ν weight terms ωi < 1 (e.g. ωi = 0.1, 0.01, . . . depending on the size of |ϵi |). In some circumstances, we only have the µ,ν Lagrange scattered data values; we set weight terms {ωi = 0, µ + ν ≥ 1}. That is the method of [6]. We can see that our method is an extension of the spline data fitting method of [6]. Next we derive the error bound for the method
∥f − s∥∞,Ω ≤ K (1 + C5 )|△|r |f |r ,∞,Ω where |f |r ,∞,Ω denotes the maximum norm of the rth derivatives of f over Ω and ∥f ∥∞,Ω is the standard infinity norm. Here |△| is the maximum of the diameters of the triangles in △, and K and C5 are constants in the error bound which will be introduced in following sections. Finally we give some probability analysis—which constitutes the biggest difference from [6]—and numerical results for our method. The paper is organized as follows. In Section 2 we review some well-known Bernstein–Bézier notation. The extension of the weighted least squares method is introduced in Section 3 together with a discussion on the solution. In Section 4 we prove that the error bound can be improved by averaging the coefficients of several splines which are constructed by fitting different sets of data. The numerical results for our method can be found in the last section. 2. Preliminaries Given a triangulation △ and integers 0 ≤ r < d, we write Sdr (△) := {s ∈ C r (Ω ) : s|T ∈ Pd , for all T ∈ △} for the usual space of splines of degree d and smoothness r, where Pd is the
d+2 2
-dimensional space of bivariate polyno-
mials of degree d. Throughout the paper we shall make extensive use of the well-known Bernstein–Bézier representation of splines. For each triangle T = ⟨v1 , v2 , v3 ⟩ in △ with vertices v1 , v2 , v3 , the corresponding polynomial piece s|T is written in the form s|T =
T T cijk Bijk ,
i+j+k=d
where the BTijk are the Bernstein–Bézier polynomials of degree d associated with T . In particular, if (λ1 , λ2 , λ3 ) are the barycentric coordinates of any point u ∈ R2 relative to the triangle T , then BTijk (u) :=
d! i!j!k!
λi1 λj2 λk3 ,
i + j + k = d·
T T As usual, we associate the Bernstein–Bézier coefficients {cijk }i+j+k=d with the domain points {ξijk := (iv1 +jv2 +kv3 )/d}i+j+k=d .
Definition 1. Let β < ∞. A triangulation △ is said to be β -quasi-uniform provided that |△| ≤ βρ△ , where ρ△ is the minimum of the radii of the incircles of the triangles of △. We say the data locations are evenly distributed over each triangle of △ if there are enough data points for every triangle such that
1/2
∥P ∥∞,T ≤ C
vi ∈T
P (vi )
2
T. Zhou, Z. Li / Computers and Mathematics with Applications 65 (2013) 1967–1977
1969
for any polynomial P and constant C . Given a vertex v , let star0 (v) = v , and starq (v) :=
{T ∈ △ : T ∩ starq−1 (v) = ∅},
q ≥ 1.
Recall that a determining set for a spline space S ⊆ Sd0 (△) is a subset M of the set of domain points such that if s ∈ S and cξ = 0 for all ξ ∈ M , then cξ = 0 for all domain points, i.e., s ≡ 0. The set M is called a minimal determining set (MDS) for S if there is no smaller determining set. It is known that M is a MDS for S if and only if every spline s ∈ S is uniquely determined by its set of B-coefficients {cξ }ξ ∈M . Next we introduce the Markov inequality. Theorem 2.1. Let T := ⟨v1 , v2 , v3 ⟩ be a triangle and fix 1 ≤ q ≤ ∞. Then there exists a constant C depending only on d and q such that, for any polynomial p and any nonnegative integers α and β with 0 ≤ α + β ≤ d,
∥Dαx Dβy p∥q,T ≤
C α+β
ρT
∥p∥q,T
where ρT denotes the radius of largest circle which can be inscribed in T . Recall from [10] that for any given function f ∈ L1 (Ω ), there exists a quasi-interpolatory operator Q mapping f ∈ L1 (Ω ) to Sdr (△) which achieves the optimal approximation order of Sdr (△). That is: Theorem 2.2 (Cf. [10]). Fix d ≥ 3r + 2 and 0 ≤ m ≤ d. Suppose f lies in the Sobolev space Wpm+1 (Ω ) with 1 ≤ p ≤ ∞. Then
∥Dαx Dβy (f − Qf )∥p,Ω ≤ K |△|m+1−α−β |f |m+1,p,Ω for 0 ≤ α + β ≤ m, where |△| is the mesh size of △ (i.e., the diameter of the largest triangle), and |f |m+1,p,Ω is the usual Sobolev semi-norm. If Ω is convex, then the constant K depends only on d, p, m, and on the smallest angle in △. If Ω is nonconvex, it also depends on the Lipschitz constant L∂ Ω associated with the boundary of Ω . 3. Discussion of a solution Let S = Sdr (△) and M = dim S with d ≥ 3r + 2. Clearly, we can see that M ≪ N, where N denotes the number of given data values. Note that using the locally supported basis functions (cf. [11]), any spline function s in the space can be M M represented by s = i=1 ci φi for the set of coefficients {ci }i=1 . The extension of the weighted least squares method is finding Sω ∈ S such that L(Sω ) = min{L(s) : s ∈ S }
(3.1)
where L(s) is defined by (1.1). First, we discuss the existence and uniqueness of the solution Sω ∈ S satisfying (3.1) and an error bound as well. µ,ν
Theorem 3.1. Suppose all the data locations are evenly distributed over each triangle of △. For fixed weight terms {ωi µ + ν ≤ r }Ni=1 , there exists a unique Sω ∈ S satisfying (3.1).
,0 ≤
Proof. The proof of the existence uses an argument similar to that in the proof of Theorem 2 in [6]. We omit the details here. T Next we show the uniqueness of the minimizer Sω . Let c = {cijk , i + j + k = d, T ∈ △} be the coefficient vector associated with s. Then N
L(s) =
µ,ν
ωi (|△|µ+ν |Dµx Dνy s(vi ) − fi
µ,ν
|)2
i=1 0≤µ+ν≤r
=
T ∈△ (xℓ ,yℓ )∈T 0≤µ+ν≤r
µ,ν ωℓ
|△|
µ+ν
2 µ,ν T µ ν T cijk Dx Dy Bijk (xℓ , yℓ ) − fℓ i+j+k=d
= cT Bc − 2bT c + ∥f∥22 µ,ν
µ,ν
1
where f = {(ωℓ ) 2 · |△|µ+ν fℓ
t
B =
0≤µ+ν≤r vℓ ∈t
µ,ν
, ℓ = 1, . . . , N , 0 ≤ µ + ν ≤ r } is a data value vector. Here B = diag(Bt , t ∈ △), with l+m+n=d
ωℓ |△|2µ+2ν Dµx Dνy Btijk (vℓ )Dµx Dνy Btlmn (vℓ ) i+j+k=d
being a matrix of size dˆ × dˆ with dˆ = btijk =
0≤µ+ν≤r vℓ ∈t
µ,ν
d+2 2
µ,ν
and b = (
btijk
ωℓ |△|2µ+2ν fℓ Dµx Dνy Btijk (vℓ ).
, i + j + k = d, t ∈ △) with
1970
T. Zhou, Z. Li / Computers and Mathematics with Applications 65 (2013) 1967–1977
Suppose that we have two solutions Sω and Sˆω . Let c and cˆ be the two coefficient vectors associated with Sω and Sˆω respectively. Since L is a convex functional, we have, for any z ∈ [0, 1], L(zSω + (1 − z )Sˆω ) ≤ zL(Sω ) + (1 − z )L(Sˆω ) = L(Sω )· That is, L(zSω +(1 − z )Sˆω ) is a constant function of z ∈ [0, 1]. It follows that 0 =
d dz
d L dz
(zSω +(1 − z )Sˆω ) = 0 for all z ∈ [0, 1]. That is,
L(zSω + (1 − z )Sˆω )
= 2z (c − cˆ )T B(c − cˆ ) − 2bT (c − cˆ ) for all z ∈ [0, 1]. It follows that
t ∈△ (xℓ ,yℓ )∈t i+j+k=d
µ,ν
t t 2 µ ν t |△|2µ+2ν ωℓ (cijk − cˆijk ) Dx Dy Bijk (vℓ )2 = 0
for any {µ, ν : 0 ≤ µ + ν ≤ r }. That is, c = cˆ , because when the data locations are evenly distributed over each triangle of △, B is invertible. Hence, the minimizer is unique. Next we derive the error bound for the solution. We convert the problem (3.1) into a standard approximation problem in Hilbert space. Let r X := {f ∈ B(Ω ) : f |T ∈ W∞ (T ), all triangles T in △}
where B(Ω ) is the set of all bounded real-valued functions on Ω . For each triangle T in △, let
⟨f , g ⟩XT :=
µ,ν
(xℓ ,yℓ )∈T 0≤µ+ν≤r
ωℓ |△|2µ+2ν Dµx Dνy f (vℓ )Dµx Dνy g (vℓ )·
Then
⟨f , g ⟩X :=
⟨f , g ⟩XT T ∈△
defines a semi-definite inner-product on X . Let ∥f ∥XT and ∥f ∥X be the associated semi-norms. Given f ∈ X , supposing that S ⊆ Sdr (△) is a Hilbert space with respect to ∥ · ∥X , we need to find Sω ∈ S such that L(Sω ) = min L(s) s∈S
where L(s) := ∥f − s∥2X . Define a projection operator P : X → S by Pf = Sω . It is easy to see that the approximation Sω of f is characterized by
⟨f − Sω , s⟩ = 0,
∀s ∈ S .
(3.2)
Moreover, using the Cauchy–Schwarz inequality, we can get
∥Pf ∥X ≤ ∥f ∥X (3.3) for all f ∈ X . More details can be found in [11,12]. Recall from [10] that when d ≥ 3r + 2, Sdr (△) possesses a stable local basis {Bξ }ξ ∈M corresponding to a minimal determining set M . Thus, for any ξ ∈ M , there is a vertex vξ of △ with supp(Bξ ) ⊆ starl (vξ ) and ∥Bξ ∥∞,Ω = 1 for all ξ ∈ M , where l is an integer number related to ξ (see Remark 6.1). Also, there exists a positive constant K1 depending only on d and β such that for any s = ξ ∈M cξ Bξ and any ξ ∈ M, K1 |cξ | ≤ ∥s∥∞,Tξ for some triangle Tξ in △ with Tξ ⊆ supp(Bξ ). Letting µ,ν
µ,ν
ωmin := min{ωi |i = 1, . . . , N , 0 ≤ µ + ν ≤ r , ωi ̸= 0}, µ,ν ω := max ωT |T ∈ △, ωT = ωi , (xi ,yi )∈T 0≤µ+ν≤r
we have Lemma 1. Let {Bξ }ξ ∈M be the basis for S corresponding to the minimal determining set M . Suppose that the data set has a property that for every s ∈ S and every T ∈ △
1/2
K2 ∥s∥∞,T ≤
|△|
2µ+2ν
|△|
2µ+2ν
µ ν
2
µ ν
2
Dx Dy s(vℓ )
,
(3.4)
vℓ ∈T 0≤µ+ν≤r
1/2
K3 ∥s∥∞,T ≥
vℓ ∈T 0≤µ+ν≤r
Dx Dy s(vℓ )
(3.5)
T. Zhou, Z. Li / Computers and Mathematics with Applications 65 (2013) 1967–1977
1971
for some positive constants K2 , K3 . Then there exist positive constants C1 and C2 depending on d, l and β such that
2 ωmin C1 cξ ≤ cξ Bξ ≤ ωC2 cξ2 ξ ∈M ξ ∈M ξ ∈M
2
(3.6)
X
for all {cξ }ξ ∈ M .
Proof. Let s = K12
ξ ∈M
ξ ∈M cξ Bξ .
cξ2 ≤
ξ ∈M
By the definition of the local stable basis we have
∥s∥2∞,Tξ .
Next we use Eq. (3.4) to get
∥s∥∞,T ≤
1/2
1
K2
1/2
1
ω
µ ν
Dx Dy s(vℓ )
ωmin K2
=
|△|
2µ+2ν
2
vℓ ∈T 0≤µ+ν≤r
1
≤
1/2 min K2
|△|
2µ+2ν
vℓ ∈T 0≤µ+ν≤r
µ,ν ωℓ |Dµx Dνy s(vℓ )|2
1/2
∥s∥XT ·
Then we have
ξ ∈M
cξ2 ≤
1
ξ ∈M
ωmin K12 K22
∥s∥2XT . ξ
Note that there may be more than one domain point η ∈ M in Tξ . Let K4 = maxT ∈△ ♯{ξ ∈ M : T ∩ supp(Bξ ) ̸= ∅}, where ♯ denotes the cardinality of set; then
ξ ∈M
cξ2 ≤
K4
ωmin K12 K22
∥s∥2X .
Therefore we have the left hand side of (3.6) with C1 =
∥s∥XT =
(xℓ ,yℓ )∈T 0≤µ+ν≤r
µ,ν
K12 K22 K4
. By using (3.5), we have
µ
ωℓ |△|2µ+2ν |Dx Dνy s(vℓ )|2
≤ ω1/2 K3 ∥s∥∞,T . Let MT := {ξ ∈ M : T ∩ supp(Bξ ) ̸= ∅}; then ♯MT ≤ K4 . Noting that ∥Bξ ∥∞,Ω = 1 for all ξ ∈ M , we have
2 cξ Bξ ∥s ∥ ≤ ∥s ∥ = T ∈△ T ∈△ ξ ∈MT XT 2 ≤ ωK32 cξ Bξ T ∈△ ξ ∈MT ∞,T 2 ≤ ωK32 |cξ | 2 X
2 XT
T ∈△
≤ω
K32
ξ ∈MT
♯MT
≤ K4 ω
|cξ |2
ξ ∈MT
T ∈△
K32
|cξ |2 .
T ∈△ ξ ∈MT
Next let K5 := maxξ ∈M ♯{T ∈ △ : T ∩ supp(Bξ ) ̸= ∅}; then
∥s∥2X ≤ ωK32 K4 K5
|cξ |2 .
ξ ∈M
Finally the proof is complete with C2 = K32 K4 K5 .
1972
T. Zhou, Z. Li / Computers and Mathematics with Applications 65 (2013) 1967–1977
Given a triangle T , let star0 (T ) = T , and starq (T ) :=
{starq (w) : w is a vertex of T },
q ≥ 1. ωC
Lemma 2. Suppose S satisfies the hypotheses of Lemma 1; let C3 = γ (γ + 1)1/2 with γ = ω 2C . Then there exists a constant min 1 0 ≤ σ ≤ 1 such that, for any triangle T ∈ △ and any function g ∈ X with supp(g ) ⊆ T ,
∥Pg ∥XT ≤ C3 σ k ∥g ∥X
(3.7)
whenever T ∈ star2(k+2)l+1 (T ) \ star2(k+1)l+1 (T ) with k ≥ 1. Proof. The proof is similar to that of Theorem 3.1 of [12]. So we put it into Remark 6.3.
Note that for any g ∈ X , we have
∥g ∥XT =
µ,ν
(xℓ ,yℓ )∈T 0≤µ+ν≤r
µ
ωℓ |△|2µ+2ν |Dx Dνy g (vℓ )|2
≤ω
1/2
∥g ∥∞,T +
|△|
µ+ν
µ ν
∥Dx Dy g ∥∞,T .
1≤µ+ν≤r
1/2 n0 + C3 k≥1 σ k nk ; then Lemma 3. Suppose that all the hypotheses of Lemma 2 hold. Let C4 = ω1/2 ωmin K2
|△|µ+ν ∥Dµx Dνy g ∥∞,Ω
∥Pg ∥∞,Ω ≤ C4 ∥g ∥∞,Ω +
(3.8)
1≤µ+ν≤r
where nk := max[star2(k+2)l+1 (T ) \ star2(k+1)l+1 (T )]
n0 = max ♯star4l+1 (T ), T ∈△
T ∈△
and C3 , σ are the constants in (3.7). Proof. Let T be any triangle in △ and
ΩkT := star2(k+2)l+1 (T ) \ star2(k+1)l+1 (T ).
Ω0T := star4l+1 (T ), Note that
∥s∥∞,T ≤
1 1/2 min K2
ω
∥s∥XT .
We can get
∥Pg ∥∞,T ≤
1 1/2
ωmin K2
n0 + C 3
σ nk ∥g ∥XT k
k≥1
as in Theorem 4.1 of [12]. Then
n0 + C 3
∥Pg ∥∞,T ≤
k≥1
1/2
ωmin K2
σ nk k
ω
1/2
∥g ∥∞,T +
|△|
µ+ν
µ ν
∥Dx Dy g ∥∞,T .
1≤µ+ν≤r
1/2 Taking the supremum over all T ∈ △, we get (3.8) with C4 = ω1/2 n0 + C3 k≥1 σ k nk . ωmin K2
Finally we have our main theorem in this paper. Theorem 3.2. Suppose △ is a β -quasi-uniform triangulation and let Sω be the spline minimizing L defined above. Let K be the (r +1)(r +2) constant in Theorem 2.2 and C5 = C4 ; we have 2
∥f − Sω ∥∞,Ω ≤ K (1 + C5 )|△|r |f |r ,∞,Ω r for every function f ∈ W∞ (Ω ).
(3.9)
T. Zhou, Z. Li / Computers and Mathematics with Applications 65 (2013) 1967–1977
1973
Proof. Let sf = Qf be the quasi-interpolation function introduced in Section 2. By the definition of Qf , we can see that Psf = sf and therefore
∥f − Sω ∥∞,Ω = ∥f − Pf ∥∞,Ω ≤ ∥f − sf ∥∞,Ω + ∥sf − Pf ∥∞,Ω ≤ ∥f − sf ∥∞,Ω + ∥Psf − Pf ∥∞,Ω . Next, by Lemma 3 and Theorem 2.2,
∥f − Sω ∥∞,Ω ≤ (1 + C4 )∥f − sf ∥∞,Ω + C4
|△|µ+ν ∥Dµx Dνy (f − sf )∥∞,Ω
1≤µ+ν≤r
≤ K (1 +
(r + 1)(r + 2)
C4 )|△|r |f |r ,∞,Ω 2 = K (1 + C5 )|△|r |f |r ,∞,Ω . Therefore the proof is complete.
µ,ν
From the proof of Theorem 3.2, we can see that if we set all ωi = 0 for µ + ν ≥ 1, we can get the error bound of the µ,ν µ,ν weighted least squares method in [6], and if we set all ωi = 1 for µ + ν = 0 and ωi = 0 for µ + ν ≥ 1, we can get the error bound of the traditional least squares method. 4. Probability analysis In this section, we will prove that the error bound can be improved by averaging the coefficients of several splines which are constructed by fitting different sets of data. µ,ν Let {{fi }j , 0 ≤ µ + ν ≤ 1, i = 1, . . . , N }M j=1 be the M sets of the scattered data with noise. Suppose each of the sets has sufficient number of scattered data and can be used to construct the spline. We use Pf j to denote the spline which is constructed from the jth data set. It is easy to see that Pf 1 , Pf 2 , . . . , Pf M is a sequence of independent random variables; they each have the same distribution with mean µ ¯ and variance σ 2 . Here we suppose that the noisy term has the uniform distribution on the interval [−1, 1]; then it is easy to see that µ ¯ = f and σ 2 = 1/3 (see Remark 6.2). According to the Central Limit Theorem (cf.[13]), we have the distribution of Pf 1 + Pf 2 + · · · + Pf M − Mf
√ σ M
tending to the Standard Normal Distribution as M −→ ∞. That is,
P
Pf 1 + Pf 2 + · · · + Pf M − Mf
√
σ M
≤ c −→ Φ (c ),
as M −→ ∞, where Φ denotes the cumulative area under the Standard Normal Distribution. We use s :=
M 1
M j =1
Pf j
to approximate the original function f . Then given a constant δ > 0, we have
Pf 1 + Pf 2 + · · · + Pf M − Mf P {|s − f | < δ} = P √ σ M
√ √ δ M δ M ∼ ≤ − 1. = 2Φ σ σ
Next we choose δ = 0.1 and δ = 0.01 to check the dependence on the number M and the probability P {|s − f | < δ} (see Tables 1 and 2). In numerical experiments, we find that when M > 12, the distribution tends to the Standard Normal Distribution as well. 5. Numerical experiments The purpose of this section is to demonstrate the performance of our method. The numerical results show that our method is much better than the traditional least squares method. Example 5.1. Suppose △ is a triangulation of the unit square domain Ω := [0, 1] × [0, 1]. Consider 2000 random points (xi , yi ) over △ as shown in Fig. 1. We use the following test functions: f1 (x, y) = x2 + 2y2 , f2 (x, y) = sin(2(x − y)), f3 (x, y) = 0.75 exp(−0.25(9x − 2)2 − 0.25(9y − 2)2 ) + 0.75 exp(−(9x + 1)2 /49 − (9y + 1)/10)
+ 0.5 exp(−0.25(9x − 7)2 − 0.25(9y − 3)2 ) − 0.2 exp(−(9x − 4)2 − (9y − 7)2 ),
1974
T. Zhou, Z. Li / Computers and Mathematics with Applications 65 (2013) 1967–1977 Table 1 The relationship between M and P {|s − f | < δ} with δ = 0.1. The probability P {|s − f | < δ}
The number M
0.99 0.95 0.90 0.85 0.80
225 128 90 69 55
Table 2 The relationship between M and P {|s − f | < δ} with δ = 0.01. The probability P {|s − f | < δ}
The number M
0.99 0.95 0.90 0.85 0.80
22 533 12 805 8 965 6 912 5 461
Table 3 The approximation errors.
f1 (x, y) f2 (x, y) f3 (x, y)
Our method
Traditional least squares method
0.1882 0.0409 0.0180
0.8536 0.6460 0.3017
Fig. 1. The scattered data and the triangulation △.
where f3 is the well-known Franke function. We use these function values and derivative values at the points (xi , yi ), i = µ µ,ν 0, . . . , 2000, to have a set of scattered Hermite data, that is, {(xi , yi , Dx Dνy f (xi , yi ) + ϵi ), 0 ≤ µ + ν ≤ 1, i = 1, . . . , 2000}. µ,ν
µ,ν
We set ϵi = 0 when the corresponding point (xi , yi ) is marked with a cross (which indicates reliable data), and set ϵi to be a random number between −1 to 1 when the corresponding point (xi , yi ) is marked with a dot (which means unreliable data). We approximate these data in spline spaces of S82 (△) using two methods: the traditional least squares method and our proposed extension of the weighted least squares method where the weights are set to be ωi = 1 and ωi = 0.01 for the reliable data and non-reliable data, respectively. Table 3 gives the maximum errors for the two methods when the computation is based on 100 × 100 equally spaced points over Ω . Discussion: From Table 3, we can see that the surface created by the traditional least squares method does not perform well on the test function. However, our method produces much better approximate solutions.
Example 5.2. The spline space, scattered data set and test functions are the same as in Example 5.1. Here, we approximate the sample functions with different weight terms. As in Example 5.1, we set the weights ωi = 1 for reliable data. Furthermore, we choose different weights ωi for unreliable data to check the dependence of the performance of our method on the weights. Some numerical results are given in Table 4.
T. Zhou, Z. Li / Computers and Mathematics with Applications 65 (2013) 1967–1977
1975
Fig. 2. Triangulations for choice N = 2, 32. Table 4 The approximation errors with different weight terms.
f1 (x, y) f2 (x, y) f3 (x, y)
ω = 0.01
ω = 0.001
ω = 0.0001
ω = 0.00001
0.1882 0.0409 0.0180
0.0414 0.0065 0.0173
0.0110 0.0020 0.0167
0.0109 6.51 × 10−4 0.0151
Table 5 Maximum errors for various test functions. N
2
4
8
16
32
f1 f2 f3
0.25221 0.06892 0.04521
0.09072 0.01793 0.01592
0.02974 0.00511 0.00532
0.01114 0.00119 0.00136
3.43 × 10−3 2.83 × 10−4 4.41 × 10−4
Discussion: From Table 4, we see that the error decreases when the weight terms decrease. For the error in approximating f3 (x, y), we observe that it does not decrease quickly when the weight terms decrease. This agrees with the prediction of the error bound (3.9) where the constant C = K (1 + C5 ) has a term which is independent of ωi and at the same time ωmin is too small. Example 5.3. Suppose ♦ is a uniform partition of the unit square domain Ω := [0, 1] × [0, 1] into N 2 subsquares. Let S82 (△♦ ) be a C 2 spline space, where △♦ is the triangulation obtained by inserting the diagonals of each of the subsquares in ♦ (see Fig. 2). We use the same test functions as in Example 5.1. Suppose that there are enough scattered Hermite data in each triangle; this means that the data locations are evenly distributed over each triangle of △♦ , and half of the data are reliable. So we set ωi = 1 and ωi = 0.01 for the reliable data and non-reliable data, respectively. We approximate these functions for choices N = 2, 4, 8, 16, 32, which correspond to repeatedly halving the mesh size. Table 5 gives the maximum error when the computation is based on 100 × 100 equally spaced points over Ω . Table 6 presents the corresponding ratios of errors for successive values of N. Discussion: From Table 6, we can see that the approximation order is close to 3, which is what we expect for order 2 convergence. Example 5.4. In this example, suppose that we have a sufficient number of scattered data. We divide these data into M sets M j and use s := M1 j=1 Pf to approximate the test functions. Here, we use the same test functions as in Example 5.1 to check the dependence of the performance of our method on M. Unlike for Example 5.1, we set all the weight terms as ω = 1 no matter whether the data are reliable or unreliable. Some numerical results are given in Table 7. Discussion: From Table 7, we can see that the error decreases when M increases. That agrees with the Central Limit Theorem. But the error in approximating f2 (x, y) by using M = 250 is bigger than that in using M = 100. This is because the formula
√
P {|s − f | < δ} ∼ = 2Φ δ
3M − 1
√
concerns probability. This means that it has 2Φ δ
3M − 1 probability of happening.
1976
T. Zhou, Z. Li / Computers and Mathematics with Applications 65 (2013) 1967–1977 Table 6 Ratios of convergence for various test functions. Ratio
2/4
4/8
8/16
16/32
f1 f2 f3
2.78 3.85 2.84
3.05 3.51 2.99
2.67 4.28 3.91
3.25 4.21 3.09
Table 7 The approximation errors with different values of M.
f 1 (x , y ) f 2 (x , y ) f 3 (x , y )
M = 20
M = 50
M = 100
M = 250
M = 500
0.0871 0.0904 0.1571
0.0452 0.0521 0.1014
0.0217 0.0212 0.0983
0.0084 0.0109 0.1006
0.0065 0.0081 0.0891
6. Remarks Remark 6.1. It is proved that there exists an integer l such that for each ξ , the support of Bξ is contained in starl (vξ ) for some vertex vξ ∈ △ in [10]. Remark 6.2. In this paper, we assume that the noisy term has a uniform distribution. So we can perform a probability analysis. Of course, we can use the same method to analyze the probability if the noisy term has another distribution with different mean and variance. Remark 6.3. In order to make the paper self-contained, We add the proof of Lemma 2 here. Proof. Let
M0T := {ξ ∈ M : supp(Bξ ) ∩ T ̸= ∅}, MqT := {ξ ∈ M : supp(Bξ ) ∩ starq(2l+1) T ̸= ∅}, N0T := M0T , NqT := MqT \ MqT−1 . Suppose that Pg is given by
Pg =
cξ Bξ ,
ξ ∈M
and let
uq :=
cξ Bξ ,
wq := Pg − uq ,
aq =
ξ ∈MqT
cξ2 ,
ξ ∈NqT
for any q ≥ 0. By (3.6), we know that
aj =
j≥q+1
ξ ̸∈MqT
cξ2 ≤
∥wq ∥2X , ωmin C1
q ≥ 0.
Now, ⟨g − Pg , wq ⟩ = 0 from (3.2), and ⟨g , wq ⟩ = 0 since wq ≡ 0 in T and g ≡ 0 outside of T . Using the fact that supp(wq ) ∩ ξ ∈MT supp(Bξ ) = ∅ for q ≥ 1, it follows that q−1
∥w ∥ = ⟨Pg − uq , wq ⟩ = ⟨g − uq , wq ⟩ = −⟨uq , wq ⟩ ∥wq ∥X , =− cξ Bξ , wq ≤ c B ξ ξ T T ξ ∈Nq ξ ∈Nq 2 q X
X
and thus by (3.6),
2 ∥wq ∥2X ≤ c B ξ ξ ≤ ω C 2 aq . ξ ∈NqT X
T. Zhou, Z. Li / Computers and Mathematics with Applications 65 (2013) 1967–1977
1977
Combining the above, we have
aj ≤
j ≥ q +1
ω C2 aq , ωmin C1
q ≥ 0.
(6.1) ωC
Then applying Lemma 2 in [14] with γ := ω 2C , we see that min 1 aq ≤ (γ + 1)ρ q a0 , where ρ := γ /(γ + 1). Since ∥Pg ∥2 ≤ ∥g ∥2 , we have
a0 ≤
j ≥0
aj =
ξ ∈M
cξ2 ≤
1
ωmin C1
∥Pg ∥2X ≤
1
ωmin C1
∥g ∥2X .
Let τ be a triangle of △ which lies outside of stark (T ). If 1 ≤ k ≤ 2l, then
∥Pg ∥Xτ ≤ ∥Pg ∥X ≤ ∥g ∥X , and (3.7) holds for any 0 ≤ σ ≤ 1 with C3 = 1. Suppose now that k ≥ 2l + 1, and let q := implies ξ ̸∈ MqT , and thus by (3.6),
∥Pg ∥
2 Xτ
k 2l+1
. Then supp(Bξ ) ∩ τ ̸= ∅
2 ω2 C22 c B cξ2 = ωC2 aj ≤ ≤ aq . ξ ξ ≤ ω C2 ωmin C1 ξ ̸∈MqT T j ≥ q + 1 ξ ̸∈Mq X
1
Combining the above estimates yields (3.7) with σ := ρ 4l+2 .
References [1] R. Arcangéli, M.C. Silanes, J.J. Torrens, Multidimensional Minimizing Splines, Kluwer Academic Publishers, 2004. [2] M.C. López Silanes, R. Arcangéli, Sur la convergence des Dm -splines d’ajustement pour des donnĺe¸es exactes ou bruitĺe¸es, Rev. Mat. Univ. Comp. Mad. 4 (2–3) (1991) 279–284. [3] D.L. Ragozin, Error bounds for derivative estimates bases on spline smoothing of exact or noisy data, J. Approx. Theory 37 (1983) 335–355. [4] R. Arcangéli, B. Ycart, Almost sure convergence of smoothing Dm -spline functions, Numer. Math. 66 (1993) 281–294. [5] A. Kouibia, M. Pasadas, Approximation by smoothing variational vector splines for noisy data, J. Comput. Appl. Math. 211 (2008) 213–222. [6] T. Zhou, D. Han, A weighted least squares method for scattered data fitting, J. Comput. Appl. Math. 217 (2008) 56–63. [7] T. Zhou, D. Han, M.J. Lai, Energy minimization method for scattered data Hermite interpolation, Appl. Numer. Math. 58 (2008) 646–659. [8] G. Awanou, M.J. Lai, P. Wenston, The multivariate spline method for scattered data fitting and numerical solutions of partial differential equations, in: G. Chen, M.J. Lai (Eds.), Wavelets and Splines, Nashboro Press, Brentwood, Tennessee, 2006, pp. 24–76. [9] T. Zhou, Z. Li, Scattered noisy data fitting using bivariate splines, Proc. Eng. 15 (2011) 1942–1946. [10] M.J. Lai, L.L. Schumaker, On the approximation power of bivariate splines, Adv. Comput. Math. 9 (1998) 251–279. [11] O. Davydov, L.L. Schumaker, On stable local bases for bivariate polynomial spline spaces, Constr. Approx. 18 (2001) 87–116. [12] M.V. Golitschek, L.L. Schumaker, Bounds on projections onto bivariate polynomial spline spaces with stable local bases, Constr. Approx 18 (2002) 241–254. [13] K. Subrahmaniam, A Primer in Probability, second ed., Mercel Dekker, Inc., 1990. Revised and Expanded. [14] C.D. Boor, A bound on the L∞ norm of L2 approximation by splines in terms of a global mesh ratio, Math. Comp. 30 (1976) 765–771.