- Email: [email protected]
Approximating surfaces by the moving least absolute deviations method
Applied Mathematics and Computation 219 (2013) 4387–4399
Contents lists available at SciVerse ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Approximating surfaces by the moving least absolute deviations method Ratko Grbic´ a, Klaudija Scitovski b, Kristian Sabo c, Rudolf Scitovski c,⇑ a b c
Faculty of Electrical Engineering, University of Osijek, Kneza Trpimira bb, HR-31 000 Osijek, Croatia Geoinfo AG, Kasernenstrasse 69, 9100 Herisau, Switzerland Department of Mathematics, University of Osijek, Trg Lj. Gaja 6, HR-31 000 Osijek, Croatia
a r t i c l e
i n f o
Keywords: Surface approximating Scattered data fitting Moving least absolute deviations LAD l1 -Norm approximation Weighted median problem
a b s t r a c t In this paper we are going to consider the problem of global data approximation on the basis of data containing outliers. For that purpose a new method entitled the moving least absolute deviations method is proposed. In the region of data in the network of knots weighted least absolute deviations local planes are constructed by means of which a global approximant is defined. The method is tested on the well known Franke’s function. An application in gridding of sonar data is also shown. Ó 2012 Elsevier Inc. All rights reserved.
1. Introduction Problems of global and/or local approximation or interpolation of scattered data obtained by measurements are considered by numerous authors (see e.g. [8,9,11,14,26]). Let us mention just a few of possible applications of such problems in applied research: geodetic observations, plotting geographical maps, determining geological layers (petroleum exploration, geological maps), water surface and water depth information, investigation of heart potentials, etc. [6,17–19,28]. One can find a rather detailed survey of these problems, possible applications and the methods used in [2,11,14]. In this paper we are going to consider the problem of global data approximation on the basis of data containing outliers. Suppose we are given a set of points in space K ¼ fT i ðxi ; yi ; zi Þ 2 R3 : i ¼ 1; . . . ; m; m P 3g X. Thereby zi are the measured values of unknown function z ¼ gðxÞ at points P i ðxi ; yi Þ. Using the given data, one has to approximate an unknown function g by a function g^ (global approximant). ^, the method often used throughout the literature is the moving least squares method In search for the global approximant g (see e.g. [9,4,15,13,31]), where local approximants are determined in terms of ordinary least squares, or the moving total least squares method (see e.g. [24,26,27]), where local approximants are determined in terms of total least squares. But, when among the measurements outliers are expected, it is better to determine local approximants in terms of least absolute deviations (LAD) (see e.g. [20]). For that purpose, we are going to define a numerically efficient scattered data approximation method we entitled the weighted least absolute deviations (WLAD) method. In addition to that, it is also possible to construct real time algorithms for this approach. In the region X we construct networks of knots. In each knot, on the basis of the WLAD method we calculate a local approximant. Since the WLAD method is an iterative procedure, a numerically efficient moving least absolute deviations (MLAD) algorithm is also defined, in which the information on the previously calculated local approximant are used for computation of some local approximant. Next, by using the blending function (see e.g. [11,14]), we can define a smooth approximation g^.
⇑ Corresponding author. E-mail addresses: [email protected] (R. Grbic´), [email protected] (K. Scitovski), [email protected] (K. Sabo), [email protected] (R. Scitovski). 0096-3003/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2012.10.041
R. Grbic´ et al. / Applied Mathematics and Computation 219 (2013) 4387–4399
4388
For a fast (but less accurate) calculation of the sequence of local approximants in the network of knots, a numerically very efficient algorithm is proposed, which is similar to the well-known Shepard’s algorithm (see e.g. [14]). In this case the local approximant in some knot is obtained as a weighted median of the data in some neighborhood of this knot. The proposed method is tested and illustrated on two numerical examples. One practical application to construction of an seafloor topography in the Adriatic is also shown. 2. The best WLAD-plane Let us now define the WLAD-plane problem (see e.g. [7,16,21,25]). Let K ¼ fT i ðxi ; yi ; zi Þ 2 R3 : i 2 Ig be a set of points in space with corresponding data weights wi > 0, where I ¼ f1; . . . ; mg; m P 3 is the set of indices. The best WLAD-plane should be determined, i.e. optimal parameters a ; b ; c 2 R of the function f ðx; y; a; b; cÞ ¼ ax þ by þ c should be determined such that
Gða ; b ; c Þ ¼ min Gða; b; cÞ;
Gða; b; cÞ ¼
ða;b;cÞ2R3
m X wi jzi axi byi cj:
ð1Þ
i¼1
This problem can also be considered as an WLAD-solution to the problem of an overdetermined system of linear equations (see e.g. [1,3,10,30]), whereas in the statistics literature (see e.g. [5]) it is considered as a parameter estimation problem for linear regression. For solving this problem it is necessary to know how to solve the weighted median problem [22,25,29]. Therefore, first we mention the following lemma [22]. Lemma 1. Let ðxi ; yi Þ; i 2 I; I ¼ f1; . . . ; mg; m P 2, be the data, where y1 6 y2 6 6 ym are real numbers, and xi > 0 corresponding data weights. Denote
( J¼
m2I:2
m X
m X
i¼1
i¼1
xi
For J – ;, let us denote
FðaÞ ¼
)
xi 6 0 :
m0 ¼ max J. Furthermore, let F : R ! R be a function defined by the formula
m X
xi jyi aj:
ð2Þ
i¼1
Then P xi ), then the minimum of function F is attained at the point aH ¼ y1 . (i) if J ¼ ; (i.e. 2x1 > m i¼1P Pm0 (ii) if J – ; and 2 i¼1 xi < m x , then the minimum of function F is attained at the point aH ¼ ym0 þ1 . Pm0 Pmi¼1 i (iii) if J – ; and 2 i¼1 xi ¼ i¼1 xi , then the minimum of function F is attained at every point aH from the segment ½ym0 ; ym0 þ1 . The number a 2 R is called the weighted median of the data ðwi ; yi Þ; i 2 I and denoted by medi2I ðwi ; yi Þ. Remark 1. Note that as a consequence of this lemma a pseudo-halving property of data ðwi ; yi Þ; i 2 I, follows directly [25], so it can be said that the weighted median of data is every number a for which it holds
X yi
wi 6
W 2
and
X yi >a
wi 6
W : 2
ð3Þ
The following important theorem gives properties of the best WLAD-plane (see e.g. [7,25]. Theorem 1 (Properties of the best WLAD-plane). Let I ¼ f1; . . . ; mg; m P 3, be a set of indices and K ¼ fT i ðxi ; yi ; zi Þ 2 R3 : i 2 Ig a set of points in space with corresponding data weights wi > 0, such that ðxi ; yi Þ – ðxj ; yj Þ ) zi – zj , whereby not all of the points T i lie on some plane parallel to the third axis. (i) Then there exists the best WLAD-plane which passes through at least three different points from K. (ii) The best weighted WLAD-plane is pseudo-halving.
Proof. Proof of statement (i) can be seen in e.g. [25]. For the purpose of proof (ii), note that due to Lemma 1 the following holds m m X X wi jzi a xi b yi cj P wi jzi a xi b yi c j; i¼1
i¼1
R. Grbic´ et al. / Applied Mathematics and Computation 219 (2013) 4387–4399
4389
where the equality holds if and only if c is the weighted median of data ðwi ; zi a xi b yi Þ; i ¼ 1; . . . ; m. According to ReP mark 1, by denoting W :¼ m i¼1 wi this means
X zi
i þb
wi 6 yi
þc
W 2
X
and zi
>a x
i þb
wi 6 yi
þc
W : 2
ð4Þ
To search for optimal parameters of the best WLAD-plane general methods of minimization without using derivatives can be used (see [12]), but their efficiency in the case of a large number of data is very low. Efficient methods for solving this problem are based on Linear Programming (see e.g. [1]). There are also various specializations of the Gauss–Newton method [10], special methods of combinatorial optimization [16], as well as methods which tend to smooth the LAD-problem by means of M-estimators. h 2.1. Method for searching the best WLAD-plane The methods we proposed try to respond to specific demands of the problem: a large number of data and reaching a solution in real time. First note that (see Theorem 1) the best WLAD-plane must pass through at least three different points from K. Thus we should require that ðxi ; yi Þ – ðxj ; yj Þ ) zi – zj , for all i; j 2 I, whereby not all of the points T i 2 K lie on some plane parallel to the third axis. Algorithm 1. (WLAD) Step 0. (Input) I ¼ f1; . . . ; mg; K ¼ fT i ðxi ; yi ; zi Þ : i 2 Ig; fwi > 0; i 2 Ig; Step 1. Choose two different points T l ; T m 2 K; Step 2. (Searching for the third point) According to Section 2.2, determine the third point T k 2 K n MðT l ; T m Þ and the plane which passes through these three points T l ; T m ; T k ; Step 3. (Preparation for a new loop) If for this plane pseudo-halving property (4) is fulfilled, STOP; Else set: T l ¼ T m ; T m ¼ T k and Go To Step 2. The stopping criterion in the mentioned algorithm is defined by the pseudo-halving property (see Theorem 1). It should be borne in mind that this is a necessary, but not a sufficient condition of the minimum. Other optimality criteria for such problems can be seen in [25]. Since in each step the value of the sum of absolute deviations is decreased and K is a finite set, this procedure in finitely many steps leads towards the best WLAD-plane. One possibility of choosing the third point is described in the next subsection. In the case of a large number of data, calculation of the weighted median of data may require a long computing time. Several fast algorithms can be seen in [10]. In numerical examples at the end of this paper the weighted median of data is calculated by a modification of the algorithm proposed in [10]. 2.2. Searching for the third point (a) Suppose that a set of points K ¼ fT i ðxi ; yi ; zi Þ 2 R3 : i 2 Ig; I ¼ f1; . . . ; mg; m P 3 in space with corresponding data weights wi > 0 is given, such that ðxi ; yi Þ – ðxj ; yj Þ ) zi – zj , whereby not all of the points T i lie on some plane parallel to the third axis. (b) Suppose that in some way two different points T l ðxl ; yl ; zl Þ; T m ðxm ; ym ; zm Þ 2 K are selected. (c) By MðT l ; T m Þ R3 denote the plane which passes through points T l ; T m
xl xm Tðx; yÞ 2 R : yl ym
(
MðT l ; T m Þ ¼
3
) x xm ¼0 y ym
ð5Þ
and by I0 I a set of all indices of points from K which lie in the plane MðT l ; T m Þ
xl x m i2I: yl ym
(
I0 ¼
) xi xm ¼ 0 : yi ym
ð6Þ
Note that l; m 2 I0 . (d1) Suppose that xl – xm . In this case the affine function defining the plane passing through points T l ; T m is of the form
f1 ðx; y; bÞ ¼ zl þ
y0 yl z0 zl ðx xl Þ þ bðy yl Þ; b x0 xl x0 xl
ð7Þ
R. Grbic´ et al. / Applied Mathematics and Computation 219 (2013) 4387–4399
4390
The optimal value of parameter b will be found by minimizing the functional
W1 ðbÞ ¼
X X yl ym zl zm wi jzi f1 ðxi ; yi ; bÞj ¼ wi zi zm ðxi xm Þ b yi ym ðxi xm Þ : x x x x l m l m i2I i2I y y
Note that term yi ym xll xmm ðxi xm Þ for all i 2 I0 vanishes. Therefore, the minimum of functional W1 will be found by minimization of the functional
^ 1 ðbÞ ¼ W
X i2InI0
yl ym zl zm wi zi zm ðxi xm Þ b yi ym ðxi xm Þ : xl xm xl xm
y y
Since yi ym xll xmm ðxi xm Þ – 0 for all i 2 I n I0 , according to Lemma 1, there exist k 2 I n I0 (which also defines the third point T k 2 K n MðT l ; T m Þ, through which this best WLAD-plane passes) and the optimal value z z
b ¼
zk zm xll xmm ðxk xm Þ y y
yk ym xll xmm ðxk xm Þ
ð8Þ
;
^ 1 (and also W1 ) is obtained. on which the minimum of functional W The affine function, which defines the plane passing through points T l ; T m ; T k is given by (2) with parameter b given by (8). (d2) Suppose yl – ym . In this case the affine function passing through points T l ; T m is of the form
! zl zm xl xm f2 ðx; y; aÞ ¼ zm þ aðx xm Þ þ a ðy ym Þ: yl ym yl ym
ð9Þ
x x
Similarly to the previous case, since for all i 2 I n I0 ; xi xm yl ym ðyi ym Þ – 0, the optimal value a of parameter a and the l
m
third point T k 2 K n MðT l ; T m Þ, through which there passes the best WLAD-plane, is obtained by solving the weighted median problem
X i2InI0
! zl zm xl xm wi zi zm ðyi ym Þ a xi xm ðyi ym Þ ! min: a yl ym yl ym
ð10Þ
In that way we obtain z z
a ¼
zk zm yll ymm ðyk ym Þ x x
xk xm yl ym ðyk ym Þ l
:
ð11Þ
m
The affine function defining the plane passing through points T l ; T m ; T k is given by (9) with parameter a given by (11). Note that due to assumption (a), xl ¼ xm and yl ¼ ym cannot occur. 3. Moving least absolute deviations method Let K ¼ fT i ðxi ; yi ; zi Þ 2 R3 : i 2 Ig be a set of points in space, where I ¼ f1; . . . ; mg is a set of indices and P i ðxi ; yi Þ 2 X # R2 are points in some region of interest. Suppose zi ’s are the measured values of the unknown function z ¼ gðx; yÞ at points Pi ðxi ; yi Þ among which outliers appears. We are going to smooth the given data by function g^. Locally, in a neighborhood of the some point Pðx; yÞ 2 X, the surface can be approximated by some basic function such as:
ðiÞ L0P ðn; g; aÞ ¼ a0 ; ðiiÞ L1P ðn; g; aÞ ¼ a0 þ a1 n þ a2 g;
ð12Þ
ðiiiÞ L2P ðn; g; aÞ ¼ a0 þ a1 n þ a2 g þ a3 ng þ a4 n2 þ a5 g2 ; where a ¼ ða0 ; . . . ; an ÞT denotes the vector of parameters.In this paper we are going to approximate function g in a neighborhood of arbitrary point Pðx; yÞ in its domain by a local approximant of the form LP : X ! R,
LP ðn; gÞ ¼ aP n þ bP g þ cP ;
ð13Þ
where parameters aP ; bP ; cP are going to be determined by using the WLAD-method. In order to do this, to each data point Pi ðxi ; yi Þ we will associate a nonnegative weight function wi : X ! ½0; 1 which will attain its maximal value equal to 1 at point P i , and its value will die out with the distance from Pi , becoming zero at all points whose distance P i is less than the given quantity r i > 0. Thereby, number ri should be chosen such that the weighted function wi is not canceled out at least at three data points.
R. Grbic´ et al. / Applied Mathematics and Computation 219 (2013) 4387–4399
4391
Various ways of choosing such weight functions can be found in the literature [9,14,26,27]. For some i 2 I we are going to define the weight function wi at arbitrary point Pðx; yÞ 2 X by
wi ðx; yÞ ¼ W
dðP i ; PÞ ; ri
ð14Þ
where dðPi ; PÞ denotes the distance from Pi ; r i > 0, and the function W : ð0; 1Þ ! ½0; 1 is a nonincreasing function such that
ðiÞ WðuÞ P 0;
ðiiÞ WðuÞ ¼ 0 for u > 1:
Following [26], we define
( WðuÞ ¼
ð1 u3 Þ3 ; 0 6 u 6 1; 0;
ð15Þ
u > 1:
The graph and ContourPlot of such function are depicted in Fig. 1. Supposing that among the data outliers appears, for the given point Pðx; yÞ 2 X, a local approximant of the form (13) is determined in terms of least absolute deviations such that
GðaP ; bP ; cP Þ ¼
m m X X wi ðx; yÞjzi LP ðxi ; yi Þj ¼ wi ðx; yÞjzi aP xi bP yi cP j ! i¼1
i¼1
min
ðaP ;bP ;cP Þ2R3
:
ð16Þ
Algorithm 1 can be used for solving this problem. Finally, we smooth the given data ðxi ; yi ; zi Þ; i ¼ 1; . . . ; m, by the function (global approximant) g^ defined by
g^ðx; yÞ :¼ LP ðx; yÞ: ^ at finitely many points. For simplicity, suppose For purposes, it suffices to calculate the value of the global approximant g that data points Pi are contained in the rectangle X ¼ ½a; b ½c; d. Take an equidistant partition ðni ; gj Þ, where
ba ba þi ; 2 2n1 dc dc þj gj ¼ ; 2 2n2 ni ¼
i ¼ 0; 1; . . . ; n1 ; ð17Þ j ¼ 0; 1; . . . ; n2 ;
where n1 and n2 are chosen such that ba dc . 2n1 2n2 In that way we have defined (see Fig. 2) a network of ð2n1 þ 1Þ ð2n2 þ 1Þ knots
K ¼ fðni ; gj Þ 2 X : i ¼ 0; 1; . . . ; n1 ; j ¼ 0; 1; . . . ; n2 g; and local approximants will be calculated only at these knots. All knots will be arranged in a sequence such that two neighboring members in the sequence are at the same time also neighboring knots in the network. In that way data referring to local approximants at the previous knot will be used for calculation of local approximants at some knot. Construction of a sequence of knots can be done in a number of ways by means of the so-called plane-filling curves [23].
Fig. 1. Weight functions.
R. Grbic´ et al. / Applied Mathematics and Computation 219 (2013) 4387–4399
4392
Fig. 2. Data points and knots.
Fig. 3. Movements at knots of the square region.
Suppose first that the region X ¼ ½a; b ½c; d is square, i.e. that a ¼ c and b ¼ d. Let n1 ¼ n2 . We will construct a sequence of knots such that we start from the central knot denoted by T 0;0 ðn0 ; g0 Þ, and then we extend it toward edges by concentric squares (see also Fig. 3):
T 0;0 ;
T k;km ; m ¼ 1; . . . 2k ; T km;k ; m ¼ 1; . . . 2k ; T k;kþm ; m ¼ 1; . . . 2k ; T kþm;k ; m ¼ 1; . . . 2k ; k ¼ 1; . . . ; n1 : ð18Þ
Members of this sequence of knots will be denoted by T 0 ; T 1 ; . . . ; T n , where n ¼ ð2n1 þ 1Þ ð2n1 þ 1Þ 1. Suppose now (see Fig. 4) that the region X ¼ ½a; b ½c; d is not square, e.g. that b a > d c (analogously in case d c > b a). Chose n1 ; n2 2 N, such that ba dc . First, ð2n2 þ 1Þ ð2n2 þ 1Þ knots T ij ðni ; gj Þ; i; j ¼ 0; 1; . . . ; n2 , where 2n1 2n2
ni ¼
da da ; þi 2 2n2
gj ¼
dc dc þj 2 2n2
placed in the square ½a; d ½c; d will be arranged in a sequence according to the previous scheme (18). After that, to this sequence we add knots that lie on the spiral as follows (see also Fig. 4):
Fig. 4. Movements at knots of the rectangular region.
R. Grbic´ et al. / Applied Mathematics and Computation 219 (2013) 4387–4399
T n2 þl;n2 m ; m ¼ 0; 1; . . . 2n2 ; T n2 þlþ1;n2 m ; m ¼ 0; 1; . . . 2n2 ;
l ¼ 1; 3; . . . ; n1 n2 :
4393
ð19Þ
In that way all knots are arranged in a sequence T 0 ; T 1 ; . . . ; T n , where n ¼ ð2n2 þ 1Þ ð2n2 þ 1Þ þ ðn1 n2 Þ n2 1. Let us denote by J ¼ f0; 1; . . . ; ng a set of all knot indices in the region X. Note that usually the number of data m n. Since the number of data m might be a very large number (in Example 2 that number is 31; 890), if problem (16) is always solved by determining parameters of local approximants, computing time would increase extremely. Therefore, in a way a subset of indices J m I (the corresponding subset of points Km K as well) will be associated to every knot T m , on the basis of which a local approximant will be calculated at that knot. Bearing this in mind, for every point P i ðxi ; yi Þ we calculate its distance to every knot T m ðnm ; gm Þ; m 2 J of the network
d1 ðPi ; T m Þ ¼ kðxi ; yi Þ ðnm ; gm Þk1 ¼ maxfjðxi nm j; jyi gm jg: In that way, for the chosen r m > 0 we define a set of indices J m I of all those points Pi 2 X whose distance to the knot T m is not greater than rm
J m ¼ fi 2 I : d1 ðPi ; T m Þ < r m g:
ð20Þ
Number rm > 0 defines the influence region of the i-th data on the construction of local approximants. Indeed, a local approximant at the knot T m is constructed on the basis of a set of points
Km ¼ fT i 2 K : i 2 Jm g;
ð21Þ
with corresponding data weights wi ¼ wi ðnm ; gm Þ; i 2 J m defined by (14) and (15). S S Note that nm¼0 J m ¼ I, i.e. nm¼0 Km ¼ K, but sets J m , i.e. Km , are not disjoint. By knowing optimal parameters ðam1 ; bm1 ; cm1 Þ of the local approximant Lm1 ðx; yÞ ¼ am1 x þ bm1 y þ cm1 at the knot T m1 , it is easy to find two different points in Km , closest to the plane determined by Lm1 , by which Algorithm 1 will be run for the purpose of searching for the local approximant Lm at the knot T m . Since knots T m1 and T m are neighboring, it can be expected in reality that corresponding local approximants differ a bit. Therefore, implementation of Algorithm 1 will be very fast. Algorithm 2. (MLAD) Step 0. Input: n; J 0 ; K0 ; Step 1. (Searching for the local approximant L0 ): Set wi ¼ wi ðn0 ; g0 Þ for i 2 J 0 and according to Algorithm 1, determine the local approximant L0 at the knot T 0 ; Set m ¼ 1; Step 2. Input: J m ; Km ; Step 3. (Searching for the local approximant Lm ): Choose two different points T l ; T m 2 Km closest to the plane determined by the local approximant Lm1 ; Set wi ¼ wi ðnm ; gm Þ for i 2 J m and according to Algorithm 1, determine the local approximant Lm at the knot T m ; Step 4. (Preparation for a new loop): If [m < n, set m ¼ m þ 1 and Go To Step 2; Else STOP].
Remark 2. Instead of calculating the local approximant at all knots of the lattice (especially in the case of a large number of data), it is meaningful to divide the region of data into disjoint subregions and then to select in a way a representative knot in each of them, at which the local approximant would be calculated. We could take a different approach when choosing the knots at which we calculated local approximants. Using the Dirichlet or Delaunay tessellation (see e.g. [11,27]), we could have subdivided the region X into subregions
X ¼ [li¼1 Xi ;
Xi disjoint
such that each subregion Xi contains roughly the same number of points. Local approximants would then be located at the center of the subregion Xi , and the data used in calculating this local approximant would be taken from a somewhat larger region containing Xi . Since in that way the number of local approximants would decrease, calculation time would speed up. On the other hand, other forms of local approximants given by (12) might be combined here as well. Remark 3. If we were are looking for the local approximant at the knot T m ðnm ; gm Þ in the form
S m ðx; yÞ ¼ cm ; (analogy to the Shepard’s approximant – see e.g. [14]), by minimizing functional
Gðcp Þ ¼
X wi ðnm ; gm Þjzi cP j; i2J m
R. Grbic´ et al. / Applied Mathematics and Computation 219 (2013) 4387–4399
4394
according to Lemma 1 we would obtain the following explicit formula for the local approximant
S m ðx; yÞ ¼ cm ¼ med fðwk ðnm ; gm Þ; zk Þ; k 2 J m g:
ð22Þ
If in Algorithm 2 local approximants are determined in this way instead of using Algorithm 1, the total computing time would decrease significantly, and the influence of present outliers would be eliminated. Hence, in cases when the searched surface does not have considerable oscillations, formula (22) can be applied successfully.
4. Numerical examples and applications The MLAD method for global data approximation on the basis of the data containing outliers will be illustrated on the next example. Example 1. Let us consider a well-known Franke’s test function (see e.g. [8]) f : ½0; 1 ½0; 1 ! R,
f ðx; yÞ ¼
3 ð9x2Þ2 þð9y2Þ2 3 ð9xþ1Þ2 þð9yþ1Þ2 1 ð9x7Þ2 þð9y3Þ2 1 ð9x4Þ2 ð9y7Þ2 4 10 4 e þ e þ e e : 4 4 2 5
Data ðxi ; yi ; zi Þ; i ¼ 1; . . . m; ðm ¼ 1200Þ will be constructed in the following way. First, on the square X ¼ ½0; 1 ½0; 1 we determine random distributed points P i ðxi ; yi Þ; i ¼ 1; . . . m, and after that we define
zi ¼ f ðxi ; yi Þ þ ei ;
ei N ½0; r2 ;
where r2 ¼ :0005, whereby on two spots we add 8 outliers each (see Fig. 5). According to (17), in the region X we introduce a network of knots with n1 ¼ n2 ¼ 7. The knots are arranged according to (18) and denoted by T 0 ; T 1 ; . . . ; T n ; n ¼ 224. According to Algorithm 1 and Algorithm 2, at the knots we calculate local approximants L0 ; L1 ; . . . ; Ln of the form (13). In that way we obtain the LAD-approximating surface which passes through points ðT k ; Lk ðT k ÞÞ; k ¼ 0; 1; . . . ; n (see Fig. 6).
1.0
0.8
0.6
0.4
0.2
0.0 0.0
0.2
0.4
0.6
0.8
1.0
0.4
0.6
0.8
1.0
Fig. 5. Franke’s function and the data.
1.0
0.8
0.6
0.4
0.2
0.0 0.0
0.2
Fig. 6. Approximating surface by WLAD-local approximants of the form (13).
R. Grbic´ et al. / Applied Mathematics and Computation 219 (2013) 4387–4399
4395
Table 1 Computing time, mean squares error and the robust coefficient of determination for Franke’s test function.
Computing time [ms] MSE R2
WLAD-loc. appr. (13)
WLAD-loc. appr. (22)
LS-loc. appr. (13)
43 .00053 .99004
15 .00149 .97511
15 .00719 .98553
1.0
0.8
0.6
0.4
0.2
0.0 0.0
0.2
0.4
0.6
0.8
1.0
Fig. 7. Approximating surface by WLAD-local approximants of the form (22).
For the purpose of comparison, in a similar way (see [14,26]) in each knot we can calculate an LS-local approximant (see Fig. 8). What is obvious is robustness of WLAD-local approximants, which ignore outliers among the data, in contrast to LSlocal approximants, by which there is a significant influence of outliers on results. As an approximation quality measure we can take the Mean Squares Error (MSE) in the data without outliers (see Table 1), or according to [19], i.e. [20], the robust coefficient of determination
R2 ¼ 1
med jzi Li ðT i Þj mad ðzi Þ
2 ;
ð23Þ
where mad ðzi Þ ¼ med jzi medj ðzj Þj is the median absolute deviation. Note that for WLAD-local approximants of the form (13) MSE is acceptable (insignificantly greater than noise ei in data). The robust coefficient of determination R2 for LS-local approximants (13) is relatively good due to a large number of data. It should be noted here that regardless of optimization of Algorithm 1 and Algorithm 2, the computing time necessary for calculation of WLAD-local approximants is considerably greater in comparison to the computing time necessary for calculation of LS-local approximants (see Table 1). If we want to shorten the computing time keeping at the same time robustness to outliers, then in each knot according to a simple formula (22) we can calculate the WLAD-local approximant of the form S m ðx; yÞ ¼ cm . In that way, approximation quality will suffer (see Fig. 7), but the computing time will be significantly shorter (see Table 1). Example 2. Our method was applied to the data obtained by measuring1 the Pakleni Channel seafloor topography in the vicinity of Hvar island in the Adriatic Sea. Measurement data ðxi ; yi ; zi Þ; i ¼ 1; . . . ; 31; 890 are shown in the local coordinate system in Fig. 9 by means of Mathematica-functions ListPlot and ListPlot3D in the region X ¼ ½920; 1570 ½90; 740, and configuration of the data in Fig. 10. The presence of outliers can be seen on two positions. According to (17), in the region X we introduce a network of knots with n1 ¼ n2 ¼ 50. Knots are arranged according to (18) and denoted by T 0 ; T 1 ; . . . ; T n ; n ¼ 10; 200. According to Algorithm 1 and Algorithm 2, at the knots we calculate local approximants L0 ; L1 ; . . . ; Ln of the form (13). In that way we obtain the LAD-approximating surface which passes through points ðT k ; Lk ðT k ÞÞ; k ¼ 0; 1; . . . ; n (see Fig. 11). Also, by means of WLAD-local approximants S 0 ; S 1 ; . . . ; S n of the form (22) we obtain the LAD-approximating surface which passes through points ðT k ; S k ðT k ÞÞ; k ¼ 0; 1; . . . ; n (see Fig. 12), and by means of the LS-local approximant of the form (13) we obtain the LS-approximating surface (see Fig. 13). As can be seen LS-approximating surface is sensitive to outliers. For each of the mentioned approximations in Table 2 the following is given: computing time, mean squares error and the robust coefficient of determination R2 .
1
Measurements were done by the Hydrographic Institute of the Republic of Croatia by means of a multibeam measuring device Kongsberg EM3002D.
4396
R. Grbic´ et al. / Applied Mathematics and Computation 219 (2013) 4387–4399
1.0
0.8
0.6
0.4
0.2
0.0 0.0
0.2
0.4
0.6
0.8
1.0
Fig. 8. Approximating surface by LS-local approximants of the form (13).
Fig. 9. Measurement data of the Pakleni Channel (ListPlot3D).
Fig. 10. Configuration of the data of the Pakleni Channel (ListPlot, ListContourPlot).
5. Global approximant In Section 3 it was shown how the approximating surface value on networks of given knots can be obtained. Next, using ^. the blending function (see e.g. [14,11], we can define a smooth global approximation g Another possibility is to follow [9] and do the convex combination of local approximants. To every local approximant Lm we associate a weight function um , which can be of the same type as the function wi defined by (14) and (15) or they can be defined in some other way (see e.g. [9,14,26]). After that, define a convex combination of local approximants C : X ! R,
R. Grbic´ et al. / Applied Mathematics and Computation 219 (2013) 4387–4399
4397
Fig. 11. Approximating surface of the Pakleni Channel by WLAD-local approximants Lm .
Fig. 12. Approximating surface of the Pakleni Channel by WLAD-local approximants S m .
Fig. 13. Approximating surface of the Pakleni Channel by LS-local approximants.
Table 2 Computing time, mean squares error and the robust coefficient of determination for Pakleni Channel data. WLAD-loc. appr. (13) Computing time [s] MSE R2
12.55 0.618 0.99997
WLAD-loc. appr. (22) 11.31 1.995 0.99915
LS-loc. appr. (13) 11.32 1.311 0.99994
R. Grbic´ et al. / Applied Mathematics and Computation 219 (2013) 4387–4399
4398
Fig. 14. Approximating Franke’s surface by convex combination with the local approximant Lm .
Cðx; yÞ ¼
n X
m¼0
u ðx; yÞ : j¼0 uj ðx; yÞ
xm ðx; yÞLm ðx; yÞ; xm ðx; yÞ ¼ Pn m
ð24Þ
Formula (24) gives an explicit expression for a function, whose surface represents a good approximation of data P ðxi ; yi ; zi Þ; i ¼ 1; . . . ; m. Note also that functions xm satisfy 0 6 xm ðx; yÞ 6 1 and nm¼1 xm ðx; yÞ ¼ 1. Example 3. By means of WLAD-local approximants obtained in illustrative Example 1 we can construct a corresponding convex combination (24) by means of WLAD-local approximants Lm of the form (13) or by means of WLAD-local approximants S m of the form (22). With a good choice of weight functions, convex combination (24) can yield a global approximant of a very good quality, but the computing procedure is very demanding. Graphs of convex combination (24) with the local approximant Lm and the local approximant S m are shown in Figs. 14 and 15, respectively. An approximative value at some point Pðx; yÞ 2 X might be calculated such that first we find the nearest knot T m , and then at point P we calculate the value of the local approximant Lm , i.e. the local approximant S m . The other possibility is to define convex combination (24) with the local approximant Lm , i.e. the local approximant S m , and calculate its value at point P i .
Fig. 15. Approximating Franke’s surface by convex combination with the local approximant S m .
Table 3 Comparison of differnt approximations of Franke’s test function. Point ðxi ; yi Þ 2 X
Nearest knot
Position of nearest knot
Ordinal number of nearest knot
f:672; :381g
{.643,.357}
{2,-2}
12
.38905
f:733; :611g
{.714,.643}
{3,2}
25
.10129
f:239; :845g
{.214,.857}
{-4,5}
111
-.00405
f:099; :079g
{.071,.071}
{-6,-6}
144
.75561
f ðxi ; yi Þ
Lm ðxi ; yi Þ (error)
S m ðxi ; yi Þ (error)
Cðxi ; yi Þ (error)
CS ðxi ; yi Þ (error)
.38705 (.002) .09246 (.009) -.00642 (.002) .74571 (.010)
.32162 (.067) .05458 (.047) .00504 (.009) .71556 (.040)
.42753 (.038) .10113 (.0002) -.00015 (.004) .78971 (.034)
.32148 (.068) .10528 (.004) -.00813 (.004) .72762 (.028)
R. Grbic´ et al. / Applied Mathematics and Computation 219 (2013) 4387–4399
4399
In Table 3, for several randomly selected points ðxi ; yi Þ from various parts of region X the following is compared: value of Franke’s function f, value of the local approximant Lm and the local approximant S m at the nearest knot, convex combination C of local approximants Lm and convex combination CS of local approximants S m . Thereby, the nearest knot is associated to every randomly selected point, the position of this knot is given in accordance with Fig. 3 as well as the ordinal number of that knot in the sequence (18). Also, a corresponding absolute error for every calculated approximation is given in brackets. References [1] A. Barrodale, F.D.K. Roberts, An improved algorithm for discrete l1 linear approximation, SIAM J. Numer. Anal. 10 (1973) 839–848. [2] M.D. Buhmann, Radial Basis Functions, Cambridge Univ. Press, Cambridge, 2003. [3] J.A. Cadzow, Minimum l1 ; l2 and l1 norm approximate solutions to an overdetermined system of linear equations, Digital Signal Process. 12 (2002) 524–560. [4] M. Carley, Moving least squares via orthogonal polynomials, SIAM J. Sci. Comput. 32 (2010) 1310–1322. [5] E. Castillo, R. Mínguez, C. Castillo, A.S. Cofinño, Dealing with the multiplicity of solutions of the l1 and l1 regression models, Eur. J. Oper. Res. 188 (2008) 460–484. [6] R. Chowdhury, B.N. Rao, Probabilistic stability assessment of slopes using high dimensional model representation, Comput. Geotech. 37 (2010) 876– 884. [7] R. Cupec, R. Grbic´, K. Sabo, R. Scitovski, Three points method for searching the best least absolute deviations plane, Appl. Math. Comput. 215 (2009) 983–994. [8] O. Davydov, R. Morandi, A. Sestini, Local hybrid approximation for scattered data fitting with bivariate splines, Comput. Aided Geom. Des. 23 (2006) 703–721. [9] R. Farwig, Multivariate interpolation of scattered data by moving least squares methods, in: J.C. Mason, M.G. Cox (Eds.), Algorithms for Approximation, Clarendon Press, Oxford, 1987, pp. 193–211. [10] C. Gurwitz, Weighted median algorithms for L1 approximation, BIT 30 (1990) 301–310. [11] H. Hagen, Topics in Surface Modelling, SIAM, Philadelphia, 1992. [12] C.T. Kelley, Iterative Methods for Optimization, SIAM, Philadelphia, 1999. [13] Z. Komargodski, D. Levin, Hermite type moving-least-squares approximations, Comput. Math. Appl. 51 (2006) 1223–1232. [14] P. Lancaster, K. Salkauskas, Curve and Surface Fitting: An Introduction, Academic Press, London, 1986. [15] D. Levin, The approximation power of moving least-squares, Math. Comput. 67 (1998) 1517–1531. [16] Y. Li, G.R. Arce, A maximum likelihood approach to least absolute deviation regression, EURASIP J. Appl. Signal Process. 12 (2004) 1762–1769. [17] D. Mirzaei, M. Dehghan, A meshless based method for solution of integral equations, Appl. Numer. Math. 60 (2010) 245–262. [18] M. Oudjene, L. Ben-Ayedb, A. Delaméziére, J.-L. Batoz, Shape optimization of clinching tools using the response surface methodology with moving least-square approximation, J. Mater. Process. Technol. 209 (2009) 289–296. [19] M. Palasean, L. Pearlstine, Estimation of water surface elevations for the Everglades, Florida, Comput. Geosci. 34 (2008) 815–826. [20] P.J. Rousseeuw, A.M. Leroy, Robust Regression and Outlier Detection, Wiley, New York, 2003. [21] K. Sabo, R. Scitovski, I. Vazler, Searching for a best LAD-solution of an overdetermined system of linear equations motivated by searching for a best LAD-hyperplane on the basis of given data, JOTA (J. Optim. Theory Appl.) 149 (2011) 293–314. [22] K. Sabo, R. Scitovski, The best least absolute deviations line – properties and two efficient methods, ANZIAM J. 50 (2008) 185–198. [23] H. Sagan, Space-Filling Curves, Springer-Verlag, New York, 1994. [24] B. Schaffrin, I. Lee, Y. Choi, Y. Felus, Total Least-Squares (TLS) for geodetic straight-line and plane adjustment, Boll. Geod. Sci. Aff. 65 (2006) 141–168. [25] A. Schöbel, Locating Lines and Hyperplanes: Theory and Algorithms, Springer Verlag, Berlin, 1999. [26] R. Scitovski, Š. Ungar, D. Jukic´, Approximating surfaces by moving total least squares method, Appl. Math. Comput. 93 (1998) 219–232. [27] R. Scitovski, Surface generating on the basis of experimental data, Math. Commun. 2 (1997) 59–67. [28] V. Sladek, J. Sladek, Local integral equations implemented by MLS-approximation and analytical integrations, Eng. Anal. Bound. Elem. 34 (2010) 904– 913. [29] I. Vazler, K. Sabo, R. Scitovski, Weighted median of the data in solving least absolute deviations problems, Communications in Statistics – Theory and Methods, 41 (2012) 1455–1465, http://dx.doi.org/10.1080/03610926.2010.539750. [30] Z. Wang, B.S. Peterson, Constrained least absolute deviation neural networks, IEEE Trans. Neural Netw. 19 (2008) 273–283. [31] H. Wendland, Local polynomial reproduction and moving least squares approximation, IMA J. Numer. Anal. 21 (2001) 285–300.
Contents lists available at SciVerse ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Approximating surfaces by the moving least absolute deviations method Ratko Grbic´ a, Klaudija Scitovski b, Kristian Sabo c, Rudolf Scitovski c,⇑ a b c
Faculty of Electrical Engineering, University of Osijek, Kneza Trpimira bb, HR-31 000 Osijek, Croatia Geoinfo AG, Kasernenstrasse 69, 9100 Herisau, Switzerland Department of Mathematics, University of Osijek, Trg Lj. Gaja 6, HR-31 000 Osijek, Croatia
a r t i c l e
i n f o
Keywords: Surface approximating Scattered data fitting Moving least absolute deviations LAD l1 -Norm approximation Weighted median problem
a b s t r a c t In this paper we are going to consider the problem of global data approximation on the basis of data containing outliers. For that purpose a new method entitled the moving least absolute deviations method is proposed. In the region of data in the network of knots weighted least absolute deviations local planes are constructed by means of which a global approximant is defined. The method is tested on the well known Franke’s function. An application in gridding of sonar data is also shown. Ó 2012 Elsevier Inc. All rights reserved.
1. Introduction Problems of global and/or local approximation or interpolation of scattered data obtained by measurements are considered by numerous authors (see e.g. [8,9,11,14,26]). Let us mention just a few of possible applications of such problems in applied research: geodetic observations, plotting geographical maps, determining geological layers (petroleum exploration, geological maps), water surface and water depth information, investigation of heart potentials, etc. [6,17–19,28]. One can find a rather detailed survey of these problems, possible applications and the methods used in [2,11,14]. In this paper we are going to consider the problem of global data approximation on the basis of data containing outliers. Suppose we are given a set of points in space K ¼ fT i ðxi ; yi ; zi Þ 2 R3 : i ¼ 1; . . . ; m; m P 3g X. Thereby zi are the measured values of unknown function z ¼ gðxÞ at points P i ðxi ; yi Þ. Using the given data, one has to approximate an unknown function g by a function g^ (global approximant). ^, the method often used throughout the literature is the moving least squares method In search for the global approximant g (see e.g. [9,4,15,13,31]), where local approximants are determined in terms of ordinary least squares, or the moving total least squares method (see e.g. [24,26,27]), where local approximants are determined in terms of total least squares. But, when among the measurements outliers are expected, it is better to determine local approximants in terms of least absolute deviations (LAD) (see e.g. [20]). For that purpose, we are going to define a numerically efficient scattered data approximation method we entitled the weighted least absolute deviations (WLAD) method. In addition to that, it is also possible to construct real time algorithms for this approach. In the region X we construct networks of knots. In each knot, on the basis of the WLAD method we calculate a local approximant. Since the WLAD method is an iterative procedure, a numerically efficient moving least absolute deviations (MLAD) algorithm is also defined, in which the information on the previously calculated local approximant are used for computation of some local approximant. Next, by using the blending function (see e.g. [11,14]), we can define a smooth approximation g^.
⇑ Corresponding author. E-mail addresses: [email protected] (R. Grbic´), [email protected] (K. Scitovski), [email protected] (K. Sabo), [email protected] (R. Scitovski). 0096-3003/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2012.10.041
R. Grbic´ et al. / Applied Mathematics and Computation 219 (2013) 4387–4399
4388
For a fast (but less accurate) calculation of the sequence of local approximants in the network of knots, a numerically very efficient algorithm is proposed, which is similar to the well-known Shepard’s algorithm (see e.g. [14]). In this case the local approximant in some knot is obtained as a weighted median of the data in some neighborhood of this knot. The proposed method is tested and illustrated on two numerical examples. One practical application to construction of an seafloor topography in the Adriatic is also shown. 2. The best WLAD-plane Let us now define the WLAD-plane problem (see e.g. [7,16,21,25]). Let K ¼ fT i ðxi ; yi ; zi Þ 2 R3 : i 2 Ig be a set of points in space with corresponding data weights wi > 0, where I ¼ f1; . . . ; mg; m P 3 is the set of indices. The best WLAD-plane should be determined, i.e. optimal parameters a ; b ; c 2 R of the function f ðx; y; a; b; cÞ ¼ ax þ by þ c should be determined such that
Gða ; b ; c Þ ¼ min Gða; b; cÞ;
Gða; b; cÞ ¼
ða;b;cÞ2R3
m X wi jzi axi byi cj:
ð1Þ
i¼1
This problem can also be considered as an WLAD-solution to the problem of an overdetermined system of linear equations (see e.g. [1,3,10,30]), whereas in the statistics literature (see e.g. [5]) it is considered as a parameter estimation problem for linear regression. For solving this problem it is necessary to know how to solve the weighted median problem [22,25,29]. Therefore, first we mention the following lemma [22]. Lemma 1. Let ðxi ; yi Þ; i 2 I; I ¼ f1; . . . ; mg; m P 2, be the data, where y1 6 y2 6 6 ym are real numbers, and xi > 0 corresponding data weights. Denote
( J¼
m2I:2
m X
m X
i¼1
i¼1
xi
For J – ;, let us denote
FðaÞ ¼
)
xi 6 0 :
m0 ¼ max J. Furthermore, let F : R ! R be a function defined by the formula
m X
xi jyi aj:
ð2Þ
i¼1
Then P xi ), then the minimum of function F is attained at the point aH ¼ y1 . (i) if J ¼ ; (i.e. 2x1 > m i¼1P Pm0 (ii) if J – ; and 2 i¼1 xi < m x , then the minimum of function F is attained at the point aH ¼ ym0 þ1 . Pm0 Pmi¼1 i (iii) if J – ; and 2 i¼1 xi ¼ i¼1 xi , then the minimum of function F is attained at every point aH from the segment ½ym0 ; ym0 þ1 . The number a 2 R is called the weighted median of the data ðwi ; yi Þ; i 2 I and denoted by medi2I ðwi ; yi Þ. Remark 1. Note that as a consequence of this lemma a pseudo-halving property of data ðwi ; yi Þ; i 2 I, follows directly [25], so it can be said that the weighted median of data is every number a for which it holds
X yi
wi 6
W 2
and
X yi >a
wi 6
W : 2
ð3Þ
The following important theorem gives properties of the best WLAD-plane (see e.g. [7,25]. Theorem 1 (Properties of the best WLAD-plane). Let I ¼ f1; . . . ; mg; m P 3, be a set of indices and K ¼ fT i ðxi ; yi ; zi Þ 2 R3 : i 2 Ig a set of points in space with corresponding data weights wi > 0, such that ðxi ; yi Þ – ðxj ; yj Þ ) zi – zj , whereby not all of the points T i lie on some plane parallel to the third axis. (i) Then there exists the best WLAD-plane which passes through at least three different points from K. (ii) The best weighted WLAD-plane is pseudo-halving.
Proof. Proof of statement (i) can be seen in e.g. [25]. For the purpose of proof (ii), note that due to Lemma 1 the following holds m m X X wi jzi a xi b yi cj P wi jzi a xi b yi c j; i¼1
i¼1
R. Grbic´ et al. / Applied Mathematics and Computation 219 (2013) 4387–4399
4389
where the equality holds if and only if c is the weighted median of data ðwi ; zi a xi b yi Þ; i ¼ 1; . . . ; m. According to ReP mark 1, by denoting W :¼ m i¼1 wi this means
X zi
i þb
wi 6 yi
þc
W 2
X
and zi
>a x
i þb
wi 6 yi
þc
W : 2
ð4Þ
To search for optimal parameters of the best WLAD-plane general methods of minimization without using derivatives can be used (see [12]), but their efficiency in the case of a large number of data is very low. Efficient methods for solving this problem are based on Linear Programming (see e.g. [1]). There are also various specializations of the Gauss–Newton method [10], special methods of combinatorial optimization [16], as well as methods which tend to smooth the LAD-problem by means of M-estimators. h 2.1. Method for searching the best WLAD-plane The methods we proposed try to respond to specific demands of the problem: a large number of data and reaching a solution in real time. First note that (see Theorem 1) the best WLAD-plane must pass through at least three different points from K. Thus we should require that ðxi ; yi Þ – ðxj ; yj Þ ) zi – zj , for all i; j 2 I, whereby not all of the points T i 2 K lie on some plane parallel to the third axis. Algorithm 1. (WLAD) Step 0. (Input) I ¼ f1; . . . ; mg; K ¼ fT i ðxi ; yi ; zi Þ : i 2 Ig; fwi > 0; i 2 Ig; Step 1. Choose two different points T l ; T m 2 K; Step 2. (Searching for the third point) According to Section 2.2, determine the third point T k 2 K n MðT l ; T m Þ and the plane which passes through these three points T l ; T m ; T k ; Step 3. (Preparation for a new loop) If for this plane pseudo-halving property (4) is fulfilled, STOP; Else set: T l ¼ T m ; T m ¼ T k and Go To Step 2. The stopping criterion in the mentioned algorithm is defined by the pseudo-halving property (see Theorem 1). It should be borne in mind that this is a necessary, but not a sufficient condition of the minimum. Other optimality criteria for such problems can be seen in [25]. Since in each step the value of the sum of absolute deviations is decreased and K is a finite set, this procedure in finitely many steps leads towards the best WLAD-plane. One possibility of choosing the third point is described in the next subsection. In the case of a large number of data, calculation of the weighted median of data may require a long computing time. Several fast algorithms can be seen in [10]. In numerical examples at the end of this paper the weighted median of data is calculated by a modification of the algorithm proposed in [10]. 2.2. Searching for the third point (a) Suppose that a set of points K ¼ fT i ðxi ; yi ; zi Þ 2 R3 : i 2 Ig; I ¼ f1; . . . ; mg; m P 3 in space with corresponding data weights wi > 0 is given, such that ðxi ; yi Þ – ðxj ; yj Þ ) zi – zj , whereby not all of the points T i lie on some plane parallel to the third axis. (b) Suppose that in some way two different points T l ðxl ; yl ; zl Þ; T m ðxm ; ym ; zm Þ 2 K are selected. (c) By MðT l ; T m Þ R3 denote the plane which passes through points T l ; T m
xl xm Tðx; yÞ 2 R : yl ym
(
MðT l ; T m Þ ¼
3
) x xm ¼0 y ym
ð5Þ
and by I0 I a set of all indices of points from K which lie in the plane MðT l ; T m Þ
xl x m i2I: yl ym
(
I0 ¼
) xi xm ¼ 0 : yi ym
ð6Þ
Note that l; m 2 I0 . (d1) Suppose that xl – xm . In this case the affine function defining the plane passing through points T l ; T m is of the form
f1 ðx; y; bÞ ¼ zl þ
y0 yl z0 zl ðx xl Þ þ bðy yl Þ; b x0 xl x0 xl
ð7Þ
R. Grbic´ et al. / Applied Mathematics and Computation 219 (2013) 4387–4399
4390
The optimal value of parameter b will be found by minimizing the functional
W1 ðbÞ ¼
X X yl ym zl zm wi jzi f1 ðxi ; yi ; bÞj ¼ wi zi zm ðxi xm Þ b yi ym ðxi xm Þ : x x x x l m l m i2I i2I y y
Note that term yi ym xll xmm ðxi xm Þ for all i 2 I0 vanishes. Therefore, the minimum of functional W1 will be found by minimization of the functional
^ 1 ðbÞ ¼ W
X i2InI0
yl ym zl zm wi zi zm ðxi xm Þ b yi ym ðxi xm Þ : xl xm xl xm
y y
Since yi ym xll xmm ðxi xm Þ – 0 for all i 2 I n I0 , according to Lemma 1, there exist k 2 I n I0 (which also defines the third point T k 2 K n MðT l ; T m Þ, through which this best WLAD-plane passes) and the optimal value z z
b ¼
zk zm xll xmm ðxk xm Þ y y
yk ym xll xmm ðxk xm Þ
ð8Þ
;
^ 1 (and also W1 ) is obtained. on which the minimum of functional W The affine function, which defines the plane passing through points T l ; T m ; T k is given by (2) with parameter b given by (8). (d2) Suppose yl – ym . In this case the affine function passing through points T l ; T m is of the form
! zl zm xl xm f2 ðx; y; aÞ ¼ zm þ aðx xm Þ þ a ðy ym Þ: yl ym yl ym
ð9Þ
x x
Similarly to the previous case, since for all i 2 I n I0 ; xi xm yl ym ðyi ym Þ – 0, the optimal value a of parameter a and the l
m
third point T k 2 K n MðT l ; T m Þ, through which there passes the best WLAD-plane, is obtained by solving the weighted median problem
X i2InI0
! zl zm xl xm wi zi zm ðyi ym Þ a xi xm ðyi ym Þ ! min: a yl ym yl ym
ð10Þ
In that way we obtain z z
a ¼
zk zm yll ymm ðyk ym Þ x x
xk xm yl ym ðyk ym Þ l
:
ð11Þ
m
The affine function defining the plane passing through points T l ; T m ; T k is given by (9) with parameter a given by (11). Note that due to assumption (a), xl ¼ xm and yl ¼ ym cannot occur. 3. Moving least absolute deviations method Let K ¼ fT i ðxi ; yi ; zi Þ 2 R3 : i 2 Ig be a set of points in space, where I ¼ f1; . . . ; mg is a set of indices and P i ðxi ; yi Þ 2 X # R2 are points in some region of interest. Suppose zi ’s are the measured values of the unknown function z ¼ gðx; yÞ at points Pi ðxi ; yi Þ among which outliers appears. We are going to smooth the given data by function g^. Locally, in a neighborhood of the some point Pðx; yÞ 2 X, the surface can be approximated by some basic function such as:
ðiÞ L0P ðn; g; aÞ ¼ a0 ; ðiiÞ L1P ðn; g; aÞ ¼ a0 þ a1 n þ a2 g;
ð12Þ
ðiiiÞ L2P ðn; g; aÞ ¼ a0 þ a1 n þ a2 g þ a3 ng þ a4 n2 þ a5 g2 ; where a ¼ ða0 ; . . . ; an ÞT denotes the vector of parameters.In this paper we are going to approximate function g in a neighborhood of arbitrary point Pðx; yÞ in its domain by a local approximant of the form LP : X ! R,
LP ðn; gÞ ¼ aP n þ bP g þ cP ;
ð13Þ
where parameters aP ; bP ; cP are going to be determined by using the WLAD-method. In order to do this, to each data point Pi ðxi ; yi Þ we will associate a nonnegative weight function wi : X ! ½0; 1 which will attain its maximal value equal to 1 at point P i , and its value will die out with the distance from Pi , becoming zero at all points whose distance P i is less than the given quantity r i > 0. Thereby, number ri should be chosen such that the weighted function wi is not canceled out at least at three data points.
R. Grbic´ et al. / Applied Mathematics and Computation 219 (2013) 4387–4399
4391
Various ways of choosing such weight functions can be found in the literature [9,14,26,27]. For some i 2 I we are going to define the weight function wi at arbitrary point Pðx; yÞ 2 X by
wi ðx; yÞ ¼ W
dðP i ; PÞ ; ri
ð14Þ
where dðPi ; PÞ denotes the distance from Pi ; r i > 0, and the function W : ð0; 1Þ ! ½0; 1 is a nonincreasing function such that
ðiÞ WðuÞ P 0;
ðiiÞ WðuÞ ¼ 0 for u > 1:
Following [26], we define
( WðuÞ ¼
ð1 u3 Þ3 ; 0 6 u 6 1; 0;
ð15Þ
u > 1:
The graph and ContourPlot of such function are depicted in Fig. 1. Supposing that among the data outliers appears, for the given point Pðx; yÞ 2 X, a local approximant of the form (13) is determined in terms of least absolute deviations such that
GðaP ; bP ; cP Þ ¼
m m X X wi ðx; yÞjzi LP ðxi ; yi Þj ¼ wi ðx; yÞjzi aP xi bP yi cP j ! i¼1
i¼1
min
ðaP ;bP ;cP Þ2R3
:
ð16Þ
Algorithm 1 can be used for solving this problem. Finally, we smooth the given data ðxi ; yi ; zi Þ; i ¼ 1; . . . ; m, by the function (global approximant) g^ defined by
g^ðx; yÞ :¼ LP ðx; yÞ: ^ at finitely many points. For simplicity, suppose For purposes, it suffices to calculate the value of the global approximant g that data points Pi are contained in the rectangle X ¼ ½a; b ½c; d. Take an equidistant partition ðni ; gj Þ, where
ba ba þi ; 2 2n1 dc dc þj gj ¼ ; 2 2n2 ni ¼
i ¼ 0; 1; . . . ; n1 ; ð17Þ j ¼ 0; 1; . . . ; n2 ;
where n1 and n2 are chosen such that ba dc . 2n1 2n2 In that way we have defined (see Fig. 2) a network of ð2n1 þ 1Þ ð2n2 þ 1Þ knots
K ¼ fðni ; gj Þ 2 X : i ¼ 0; 1; . . . ; n1 ; j ¼ 0; 1; . . . ; n2 g; and local approximants will be calculated only at these knots. All knots will be arranged in a sequence such that two neighboring members in the sequence are at the same time also neighboring knots in the network. In that way data referring to local approximants at the previous knot will be used for calculation of local approximants at some knot. Construction of a sequence of knots can be done in a number of ways by means of the so-called plane-filling curves [23].
Fig. 1. Weight functions.
R. Grbic´ et al. / Applied Mathematics and Computation 219 (2013) 4387–4399
4392
Fig. 2. Data points and knots.
Fig. 3. Movements at knots of the square region.
Suppose first that the region X ¼ ½a; b ½c; d is square, i.e. that a ¼ c and b ¼ d. Let n1 ¼ n2 . We will construct a sequence of knots such that we start from the central knot denoted by T 0;0 ðn0 ; g0 Þ, and then we extend it toward edges by concentric squares (see also Fig. 3):
T 0;0 ;
T k;km ; m ¼ 1; . . . 2k ; T km;k ; m ¼ 1; . . . 2k ; T k;kþm ; m ¼ 1; . . . 2k ; T kþm;k ; m ¼ 1; . . . 2k ; k ¼ 1; . . . ; n1 : ð18Þ
Members of this sequence of knots will be denoted by T 0 ; T 1 ; . . . ; T n , where n ¼ ð2n1 þ 1Þ ð2n1 þ 1Þ 1. Suppose now (see Fig. 4) that the region X ¼ ½a; b ½c; d is not square, e.g. that b a > d c (analogously in case d c > b a). Chose n1 ; n2 2 N, such that ba dc . First, ð2n2 þ 1Þ ð2n2 þ 1Þ knots T ij ðni ; gj Þ; i; j ¼ 0; 1; . . . ; n2 , where 2n1 2n2
ni ¼
da da ; þi 2 2n2
gj ¼
dc dc þj 2 2n2
placed in the square ½a; d ½c; d will be arranged in a sequence according to the previous scheme (18). After that, to this sequence we add knots that lie on the spiral as follows (see also Fig. 4):
Fig. 4. Movements at knots of the rectangular region.
R. Grbic´ et al. / Applied Mathematics and Computation 219 (2013) 4387–4399
T n2 þl;n2 m ; m ¼ 0; 1; . . . 2n2 ; T n2 þlþ1;n2 m ; m ¼ 0; 1; . . . 2n2 ;
l ¼ 1; 3; . . . ; n1 n2 :
4393
ð19Þ
In that way all knots are arranged in a sequence T 0 ; T 1 ; . . . ; T n , where n ¼ ð2n2 þ 1Þ ð2n2 þ 1Þ þ ðn1 n2 Þ n2 1. Let us denote by J ¼ f0; 1; . . . ; ng a set of all knot indices in the region X. Note that usually the number of data m n. Since the number of data m might be a very large number (in Example 2 that number is 31; 890), if problem (16) is always solved by determining parameters of local approximants, computing time would increase extremely. Therefore, in a way a subset of indices J m I (the corresponding subset of points Km K as well) will be associated to every knot T m , on the basis of which a local approximant will be calculated at that knot. Bearing this in mind, for every point P i ðxi ; yi Þ we calculate its distance to every knot T m ðnm ; gm Þ; m 2 J of the network
d1 ðPi ; T m Þ ¼ kðxi ; yi Þ ðnm ; gm Þk1 ¼ maxfjðxi nm j; jyi gm jg: In that way, for the chosen r m > 0 we define a set of indices J m I of all those points Pi 2 X whose distance to the knot T m is not greater than rm
J m ¼ fi 2 I : d1 ðPi ; T m Þ < r m g:
ð20Þ
Number rm > 0 defines the influence region of the i-th data on the construction of local approximants. Indeed, a local approximant at the knot T m is constructed on the basis of a set of points
Km ¼ fT i 2 K : i 2 Jm g;
ð21Þ
with corresponding data weights wi ¼ wi ðnm ; gm Þ; i 2 J m defined by (14) and (15). S S Note that nm¼0 J m ¼ I, i.e. nm¼0 Km ¼ K, but sets J m , i.e. Km , are not disjoint. By knowing optimal parameters ðam1 ; bm1 ; cm1 Þ of the local approximant Lm1 ðx; yÞ ¼ am1 x þ bm1 y þ cm1 at the knot T m1 , it is easy to find two different points in Km , closest to the plane determined by Lm1 , by which Algorithm 1 will be run for the purpose of searching for the local approximant Lm at the knot T m . Since knots T m1 and T m are neighboring, it can be expected in reality that corresponding local approximants differ a bit. Therefore, implementation of Algorithm 1 will be very fast. Algorithm 2. (MLAD) Step 0. Input: n; J 0 ; K0 ; Step 1. (Searching for the local approximant L0 ): Set wi ¼ wi ðn0 ; g0 Þ for i 2 J 0 and according to Algorithm 1, determine the local approximant L0 at the knot T 0 ; Set m ¼ 1; Step 2. Input: J m ; Km ; Step 3. (Searching for the local approximant Lm ): Choose two different points T l ; T m 2 Km closest to the plane determined by the local approximant Lm1 ; Set wi ¼ wi ðnm ; gm Þ for i 2 J m and according to Algorithm 1, determine the local approximant Lm at the knot T m ; Step 4. (Preparation for a new loop): If [m < n, set m ¼ m þ 1 and Go To Step 2; Else STOP].
Remark 2. Instead of calculating the local approximant at all knots of the lattice (especially in the case of a large number of data), it is meaningful to divide the region of data into disjoint subregions and then to select in a way a representative knot in each of them, at which the local approximant would be calculated. We could take a different approach when choosing the knots at which we calculated local approximants. Using the Dirichlet or Delaunay tessellation (see e.g. [11,27]), we could have subdivided the region X into subregions
X ¼ [li¼1 Xi ;
Xi disjoint
such that each subregion Xi contains roughly the same number of points. Local approximants would then be located at the center of the subregion Xi , and the data used in calculating this local approximant would be taken from a somewhat larger region containing Xi . Since in that way the number of local approximants would decrease, calculation time would speed up. On the other hand, other forms of local approximants given by (12) might be combined here as well. Remark 3. If we were are looking for the local approximant at the knot T m ðnm ; gm Þ in the form
S m ðx; yÞ ¼ cm ; (analogy to the Shepard’s approximant – see e.g. [14]), by minimizing functional
Gðcp Þ ¼
X wi ðnm ; gm Þjzi cP j; i2J m
R. Grbic´ et al. / Applied Mathematics and Computation 219 (2013) 4387–4399
4394
according to Lemma 1 we would obtain the following explicit formula for the local approximant
S m ðx; yÞ ¼ cm ¼ med fðwk ðnm ; gm Þ; zk Þ; k 2 J m g:
ð22Þ
If in Algorithm 2 local approximants are determined in this way instead of using Algorithm 1, the total computing time would decrease significantly, and the influence of present outliers would be eliminated. Hence, in cases when the searched surface does not have considerable oscillations, formula (22) can be applied successfully.
4. Numerical examples and applications The MLAD method for global data approximation on the basis of the data containing outliers will be illustrated on the next example. Example 1. Let us consider a well-known Franke’s test function (see e.g. [8]) f : ½0; 1 ½0; 1 ! R,
f ðx; yÞ ¼
3 ð9x2Þ2 þð9y2Þ2 3 ð9xþ1Þ2 þð9yþ1Þ2 1 ð9x7Þ2 þð9y3Þ2 1 ð9x4Þ2 ð9y7Þ2 4 10 4 e þ e þ e e : 4 4 2 5
Data ðxi ; yi ; zi Þ; i ¼ 1; . . . m; ðm ¼ 1200Þ will be constructed in the following way. First, on the square X ¼ ½0; 1 ½0; 1 we determine random distributed points P i ðxi ; yi Þ; i ¼ 1; . . . m, and after that we define
zi ¼ f ðxi ; yi Þ þ ei ;
ei N ½0; r2 ;
where r2 ¼ :0005, whereby on two spots we add 8 outliers each (see Fig. 5). According to (17), in the region X we introduce a network of knots with n1 ¼ n2 ¼ 7. The knots are arranged according to (18) and denoted by T 0 ; T 1 ; . . . ; T n ; n ¼ 224. According to Algorithm 1 and Algorithm 2, at the knots we calculate local approximants L0 ; L1 ; . . . ; Ln of the form (13). In that way we obtain the LAD-approximating surface which passes through points ðT k ; Lk ðT k ÞÞ; k ¼ 0; 1; . . . ; n (see Fig. 6).
1.0
0.8
0.6
0.4
0.2
0.0 0.0
0.2
0.4
0.6
0.8
1.0
0.4
0.6
0.8
1.0
Fig. 5. Franke’s function and the data.
1.0
0.8
0.6
0.4
0.2
0.0 0.0
0.2
Fig. 6. Approximating surface by WLAD-local approximants of the form (13).
R. Grbic´ et al. / Applied Mathematics and Computation 219 (2013) 4387–4399
4395
Table 1 Computing time, mean squares error and the robust coefficient of determination for Franke’s test function.
Computing time [ms] MSE R2
WLAD-loc. appr. (13)
WLAD-loc. appr. (22)
LS-loc. appr. (13)
43 .00053 .99004
15 .00149 .97511
15 .00719 .98553
1.0
0.8
0.6
0.4
0.2
0.0 0.0
0.2
0.4
0.6
0.8
1.0
Fig. 7. Approximating surface by WLAD-local approximants of the form (22).
For the purpose of comparison, in a similar way (see [14,26]) in each knot we can calculate an LS-local approximant (see Fig. 8). What is obvious is robustness of WLAD-local approximants, which ignore outliers among the data, in contrast to LSlocal approximants, by which there is a significant influence of outliers on results. As an approximation quality measure we can take the Mean Squares Error (MSE) in the data without outliers (see Table 1), or according to [19], i.e. [20], the robust coefficient of determination
R2 ¼ 1
med jzi Li ðT i Þj mad ðzi Þ
2 ;
ð23Þ
where mad ðzi Þ ¼ med jzi medj ðzj Þj is the median absolute deviation. Note that for WLAD-local approximants of the form (13) MSE is acceptable (insignificantly greater than noise ei in data). The robust coefficient of determination R2 for LS-local approximants (13) is relatively good due to a large number of data. It should be noted here that regardless of optimization of Algorithm 1 and Algorithm 2, the computing time necessary for calculation of WLAD-local approximants is considerably greater in comparison to the computing time necessary for calculation of LS-local approximants (see Table 1). If we want to shorten the computing time keeping at the same time robustness to outliers, then in each knot according to a simple formula (22) we can calculate the WLAD-local approximant of the form S m ðx; yÞ ¼ cm . In that way, approximation quality will suffer (see Fig. 7), but the computing time will be significantly shorter (see Table 1). Example 2. Our method was applied to the data obtained by measuring1 the Pakleni Channel seafloor topography in the vicinity of Hvar island in the Adriatic Sea. Measurement data ðxi ; yi ; zi Þ; i ¼ 1; . . . ; 31; 890 are shown in the local coordinate system in Fig. 9 by means of Mathematica-functions ListPlot and ListPlot3D in the region X ¼ ½920; 1570 ½90; 740, and configuration of the data in Fig. 10. The presence of outliers can be seen on two positions. According to (17), in the region X we introduce a network of knots with n1 ¼ n2 ¼ 50. Knots are arranged according to (18) and denoted by T 0 ; T 1 ; . . . ; T n ; n ¼ 10; 200. According to Algorithm 1 and Algorithm 2, at the knots we calculate local approximants L0 ; L1 ; . . . ; Ln of the form (13). In that way we obtain the LAD-approximating surface which passes through points ðT k ; Lk ðT k ÞÞ; k ¼ 0; 1; . . . ; n (see Fig. 11). Also, by means of WLAD-local approximants S 0 ; S 1 ; . . . ; S n of the form (22) we obtain the LAD-approximating surface which passes through points ðT k ; S k ðT k ÞÞ; k ¼ 0; 1; . . . ; n (see Fig. 12), and by means of the LS-local approximant of the form (13) we obtain the LS-approximating surface (see Fig. 13). As can be seen LS-approximating surface is sensitive to outliers. For each of the mentioned approximations in Table 2 the following is given: computing time, mean squares error and the robust coefficient of determination R2 .
1
Measurements were done by the Hydrographic Institute of the Republic of Croatia by means of a multibeam measuring device Kongsberg EM3002D.
4396
R. Grbic´ et al. / Applied Mathematics and Computation 219 (2013) 4387–4399
1.0
0.8
0.6
0.4
0.2
0.0 0.0
0.2
0.4
0.6
0.8
1.0
Fig. 8. Approximating surface by LS-local approximants of the form (13).
Fig. 9. Measurement data of the Pakleni Channel (ListPlot3D).
Fig. 10. Configuration of the data of the Pakleni Channel (ListPlot, ListContourPlot).
5. Global approximant In Section 3 it was shown how the approximating surface value on networks of given knots can be obtained. Next, using ^. the blending function (see e.g. [14,11], we can define a smooth global approximation g Another possibility is to follow [9] and do the convex combination of local approximants. To every local approximant Lm we associate a weight function um , which can be of the same type as the function wi defined by (14) and (15) or they can be defined in some other way (see e.g. [9,14,26]). After that, define a convex combination of local approximants C : X ! R,
R. Grbic´ et al. / Applied Mathematics and Computation 219 (2013) 4387–4399
4397
Fig. 11. Approximating surface of the Pakleni Channel by WLAD-local approximants Lm .
Fig. 12. Approximating surface of the Pakleni Channel by WLAD-local approximants S m .
Fig. 13. Approximating surface of the Pakleni Channel by LS-local approximants.
Table 2 Computing time, mean squares error and the robust coefficient of determination for Pakleni Channel data. WLAD-loc. appr. (13) Computing time [s] MSE R2
12.55 0.618 0.99997
WLAD-loc. appr. (22) 11.31 1.995 0.99915
LS-loc. appr. (13) 11.32 1.311 0.99994
R. Grbic´ et al. / Applied Mathematics and Computation 219 (2013) 4387–4399
4398
Fig. 14. Approximating Franke’s surface by convex combination with the local approximant Lm .
Cðx; yÞ ¼
n X
m¼0
u ðx; yÞ : j¼0 uj ðx; yÞ
xm ðx; yÞLm ðx; yÞ; xm ðx; yÞ ¼ Pn m
ð24Þ
Formula (24) gives an explicit expression for a function, whose surface represents a good approximation of data P ðxi ; yi ; zi Þ; i ¼ 1; . . . ; m. Note also that functions xm satisfy 0 6 xm ðx; yÞ 6 1 and nm¼1 xm ðx; yÞ ¼ 1. Example 3. By means of WLAD-local approximants obtained in illustrative Example 1 we can construct a corresponding convex combination (24) by means of WLAD-local approximants Lm of the form (13) or by means of WLAD-local approximants S m of the form (22). With a good choice of weight functions, convex combination (24) can yield a global approximant of a very good quality, but the computing procedure is very demanding. Graphs of convex combination (24) with the local approximant Lm and the local approximant S m are shown in Figs. 14 and 15, respectively. An approximative value at some point Pðx; yÞ 2 X might be calculated such that first we find the nearest knot T m , and then at point P we calculate the value of the local approximant Lm , i.e. the local approximant S m . The other possibility is to define convex combination (24) with the local approximant Lm , i.e. the local approximant S m , and calculate its value at point P i .
Fig. 15. Approximating Franke’s surface by convex combination with the local approximant S m .
Table 3 Comparison of differnt approximations of Franke’s test function. Point ðxi ; yi Þ 2 X
Nearest knot
Position of nearest knot
Ordinal number of nearest knot
f:672; :381g
{.643,.357}
{2,-2}
12
.38905
f:733; :611g
{.714,.643}
{3,2}
25
.10129
f:239; :845g
{.214,.857}
{-4,5}
111
-.00405
f:099; :079g
{.071,.071}
{-6,-6}
144
.75561
f ðxi ; yi Þ
Lm ðxi ; yi Þ (error)
S m ðxi ; yi Þ (error)
Cðxi ; yi Þ (error)
CS ðxi ; yi Þ (error)
.38705 (.002) .09246 (.009) -.00642 (.002) .74571 (.010)
.32162 (.067) .05458 (.047) .00504 (.009) .71556 (.040)
.42753 (.038) .10113 (.0002) -.00015 (.004) .78971 (.034)
.32148 (.068) .10528 (.004) -.00813 (.004) .72762 (.028)
R. Grbic´ et al. / Applied Mathematics and Computation 219 (2013) 4387–4399
4399
In Table 3, for several randomly selected points ðxi ; yi Þ from various parts of region X the following is compared: value of Franke’s function f, value of the local approximant Lm and the local approximant S m at the nearest knot, convex combination C of local approximants Lm and convex combination CS of local approximants S m . Thereby, the nearest knot is associated to every randomly selected point, the position of this knot is given in accordance with Fig. 3 as well as the ordinal number of that knot in the sequence (18). Also, a corresponding absolute error for every calculated approximation is given in brackets. References [1] A. Barrodale, F.D.K. Roberts, An improved algorithm for discrete l1 linear approximation, SIAM J. Numer. Anal. 10 (1973) 839–848. [2] M.D. Buhmann, Radial Basis Functions, Cambridge Univ. Press, Cambridge, 2003. [3] J.A. Cadzow, Minimum l1 ; l2 and l1 norm approximate solutions to an overdetermined system of linear equations, Digital Signal Process. 12 (2002) 524–560. [4] M. Carley, Moving least squares via orthogonal polynomials, SIAM J. Sci. Comput. 32 (2010) 1310–1322. [5] E. Castillo, R. Mínguez, C. Castillo, A.S. Cofinño, Dealing with the multiplicity of solutions of the l1 and l1 regression models, Eur. J. Oper. Res. 188 (2008) 460–484. [6] R. Chowdhury, B.N. Rao, Probabilistic stability assessment of slopes using high dimensional model representation, Comput. Geotech. 37 (2010) 876– 884. [7] R. Cupec, R. Grbic´, K. Sabo, R. Scitovski, Three points method for searching the best least absolute deviations plane, Appl. Math. Comput. 215 (2009) 983–994. [8] O. Davydov, R. Morandi, A. Sestini, Local hybrid approximation for scattered data fitting with bivariate splines, Comput. Aided Geom. Des. 23 (2006) 703–721. [9] R. Farwig, Multivariate interpolation of scattered data by moving least squares methods, in: J.C. Mason, M.G. Cox (Eds.), Algorithms for Approximation, Clarendon Press, Oxford, 1987, pp. 193–211. [10] C. Gurwitz, Weighted median algorithms for L1 approximation, BIT 30 (1990) 301–310. [11] H. Hagen, Topics in Surface Modelling, SIAM, Philadelphia, 1992. [12] C.T. Kelley, Iterative Methods for Optimization, SIAM, Philadelphia, 1999. [13] Z. Komargodski, D. Levin, Hermite type moving-least-squares approximations, Comput. Math. Appl. 51 (2006) 1223–1232. [14] P. Lancaster, K. Salkauskas, Curve and Surface Fitting: An Introduction, Academic Press, London, 1986. [15] D. Levin, The approximation power of moving least-squares, Math. Comput. 67 (1998) 1517–1531. [16] Y. Li, G.R. Arce, A maximum likelihood approach to least absolute deviation regression, EURASIP J. Appl. Signal Process. 12 (2004) 1762–1769. [17] D. Mirzaei, M. Dehghan, A meshless based method for solution of integral equations, Appl. Numer. Math. 60 (2010) 245–262. [18] M. Oudjene, L. Ben-Ayedb, A. Delaméziére, J.-L. Batoz, Shape optimization of clinching tools using the response surface methodology with moving least-square approximation, J. Mater. Process. Technol. 209 (2009) 289–296. [19] M. Palasean, L. Pearlstine, Estimation of water surface elevations for the Everglades, Florida, Comput. Geosci. 34 (2008) 815–826. [20] P.J. Rousseeuw, A.M. Leroy, Robust Regression and Outlier Detection, Wiley, New York, 2003. [21] K. Sabo, R. Scitovski, I. Vazler, Searching for a best LAD-solution of an overdetermined system of linear equations motivated by searching for a best LAD-hyperplane on the basis of given data, JOTA (J. Optim. Theory Appl.) 149 (2011) 293–314. [22] K. Sabo, R. Scitovski, The best least absolute deviations line – properties and two efficient methods, ANZIAM J. 50 (2008) 185–198. [23] H. Sagan, Space-Filling Curves, Springer-Verlag, New York, 1994. [24] B. Schaffrin, I. Lee, Y. Choi, Y. Felus, Total Least-Squares (TLS) for geodetic straight-line and plane adjustment, Boll. Geod. Sci. Aff. 65 (2006) 141–168. [25] A. Schöbel, Locating Lines and Hyperplanes: Theory and Algorithms, Springer Verlag, Berlin, 1999. [26] R. Scitovski, Š. Ungar, D. Jukic´, Approximating surfaces by moving total least squares method, Appl. Math. Comput. 93 (1998) 219–232. [27] R. Scitovski, Surface generating on the basis of experimental data, Math. Commun. 2 (1997) 59–67. [28] V. Sladek, J. Sladek, Local integral equations implemented by MLS-approximation and analytical integrations, Eng. Anal. Bound. Elem. 34 (2010) 904– 913. [29] I. Vazler, K. Sabo, R. Scitovski, Weighted median of the data in solving least absolute deviations problems, Communications in Statistics – Theory and Methods, 41 (2012) 1455–1465, http://dx.doi.org/10.1080/03610926.2010.539750. [30] Z. Wang, B.S. Peterson, Constrained least absolute deviation neural networks, IEEE Trans. Neural Netw. 19 (2008) 273–283. [31] H. Wendland, Local polynomial reproduction and moving least squares approximation, IMA J. Numer. Anal. 21 (2001) 285–300.