Copyright © CC-BY-NC 2019, CRIBFB | AFBR Australian Finance & Banking Review; Vol. 3, No. 1; 2019 ISSN 2576-1196 E-ISSN 2576-120X Published by Centre for Research on Islamic Banking & Finance and Business, USA 43 L1‎‎ Norm Based Data Analysis and Related Methods (1632-1989) Bijan Bidabad 1 Abstract This paper gives a rather general view on the L1‎‎ norm criterion on the area of data analysis and related topics. We tried to cover all aspects of mathematical properties, historical development, computational algorithms, simultaneous equations estimation, statistical modeling, and application of the L1‎‎ norm in different fields of sciences. Keywords: L1 norm, Regression, Algorithm, Computer 1. Introduction Although the L1‎‎ norm is an old topic in science, lack of a general book or paper on this subject induced me to gather a relatively complete list of references in this paper. The methods related to L1‎‎ norm are very broad and summarizing them is very difficult. However, it has been tried to have a glance at almost all related areas. The sections are designed as separate modules, so that the reader may skip some of the sections without loss of continuity of the subject. While the least squares method of estimation of the regression parameters is the most commonly used procedure, some alternative techniques have received widespread attention in recent years. Conventionally, interest in other methods of estimation has been generated by the unsatisfactory performance of least squares estimators in certain situations when some model assumptions fail to hold or when large correlations exist among the regressors. However, the least squares regression is very far from optimal in many non-Gaussian situations, especially when the errors follow distributions with longer tails. In particular, when the variance of the error is infinite. While intuition may dispel consideration of errors with infinite variance, in many cases, studies have shown that, in fact, certain distributions with infinite variances may be quite appropriate models. An infinite variance means thick tail error distribution with lots of outliers. Of course, observed distributions of economic variables will never display infinite variances. However, the important issue is not that the second moment is actually infinite, but the interdecile range in relation to the interquartile range is sufficiently large that one is justified in acting as though the variance is infinite. Even when the majority of the errors in the model follow a normal distribution, it often occurs that a small number of observations are from a different distribution. That is, the sample is contaminated with outliers. Since least squares gives a lot of weight to outliers, it becomes extremely sample dependent and it is well known that the performance of this estimator is markedly degraded in this situation. It has been stated that even when errors follow a normal distribution, an alternative to least squares may be required; especially if the form of the model is not exactly known or any other specification error exists. Further, least squares is not very satisfactory if the quadratic loss function is not a satisfactory measure of loss. Loss denotes the seriousness of the nonzero prediction error to the investigator, where prediction error is the difference between the predicted and the observed values of the response variable. It has been shown that for certain economic problems least absolute errors gives more satisfactory results than least squares, because the former is less sensitive than the latter to extreme errors, and consequently is resistant to outliers. It should be noted that the least absolute errors estimates have maximum likelihood properties and hence are asymptotically efficient when the errors follow the Laplace distribution. Although least absolute errors estimator is very old, it has emerged in the literature again and has attracted attention in the last two decades because of unsatisfactory properties of least squares. Now, this method is discussed in econometrics textbooks such as Kmenta (1986) and Maddala (1977). Many Master's and Ph.D. dissertations have been written on this subject in different departments such as Lawson (1961), Burgoyne (1965), Gentleman (1965), Barrodale (1967), Oveson (1968), Lewis (1969), Cline (1970), Hunt (1970), Groucher (1971), Henriksson (1972), 1 (B.A., M.Sc., Ph.D., Post-Doc.) Professor of Economics and Chief Economic Advisor to Bank Melli Iran. http://www.bidabad.com bijan@bidabad.com bidabad@yahoo.com http://www.bidabad.com/ mailto:bijan@bidabad.com mailto:bidabad@yahoo.com Copyright © CC-BY-NC 2019, CRIBFB | AFBR www.cribfb.com/journal/index.php/afbr Australian Finance & Banking Review Vol. 3, No. 1; 2019 44 Bassett (1973), Forth (1974), Anderson (1975), Ronner (1977), Nyquist (1980), Clarke (1981), Kotiuga (1981), Gonin (1983), Busovaca (1985), Kim ( ), Bidabad (1989a,b) which are more recent (see bibliography for the corresponding departments and universities). Robust property of this estimator is its advantage to deal with large variance error distributions. Since many economic phenomena such as distribution of personal income, security returns, speculative prices, stock and commodity prices, employment, asset size of business firms, demand equations, interest rate, treasury cash flows, insurance, price expectations, and many other economic variables fall within the category of infinite variance (see, Ganger and Orr (1972), Nyquist and Westlund (1977), Fama (1965), Goldfeld and Quandt (1981), Sharpe (1971)) it is necessary to turn the economists attention to this estimator. There are many other works which confirm the superiority of least absolute to least squares estimator such as interindustry demand analysis of Arrow and Hoffenberg (1959), investment models of Meyer and Glauber (1964), security and portfolios analysis of Sharpe (1971), Danish investment analysis of Kaergard (1987) and so forth. Many new economic theories weaken the assumption of rationality of human behavior. This relative irrationality is a major source of large variances and outliers in economic data. Therefore the least absolute errors estimator becomes a relevant estimator in the cases that rationality is a strong assumption. Another major application of this estimator is on data with measurement errors. This type of errors makes variances large and forces the observations to locate far from reality, which obviously causes outliers. Existence of two important types of measurement errors, sampling, and nonsampling errors, specifically in countries with poor statistics such as developing countries make this estimator a basic tool of analysis. Unknown specification errors in regression models because of the complexity of human behavior always occurs in the mathematical formulation of human-related problems. Specification error occurs whenever the formulation of the regression equation or one of the underlying assumptions is incorrect. In this context when any assumption of the underlying theory or the formulation of the model does not hold, a relevant explanatory variable is omitted, or an irrelevant one is included, qualitative change of the explanatory variable is disregarded, incorrect mathematical form of the regression is adopted, or incorrect specification of the way in which the disturbance enters the regression equation is used and so on; specification error exits (see also, Kmenta (1986)). Since specification errors are not always clear to researcher, least squares is a poor estimator, and other alternatives as least absolute errors estimators become attractive. Although the least absolute errors estimator benefits from optimal properties in many econometric problems, it is not a commonly used tool. This is to some extent due to difficulties of calculus with absolute value functions. When the model is enlarged, and equations enter simultaneously, difficulties of computation increase. Another problem with this estimator is that the properties of the solution space is not completely clear, and the corresponding closed form of the solution have not been derived yet. Thus these three important problems of algebraic closed form, computational difficulties, and solution space properties are the main obstacles that prevent the regular use of L1‎‎ norm estimator. Any attempt to remove these obstacles are worthy. 2. Lp norm and regression analysis Given a point u=(u1,...,un) in Rn, Minkowski norm or Lp norm can be written as the following expression, n ││u││p = dp(u,0) =[ │ui│ p]1/p (1) i=1 When p=2, we are confronted with Euclidian or L2 norm. Thus, Euclidian distance is a special case of Lp distance (see, Ralston and Rabinowitz (1985)). The following overdetermined system of equations is given, y = Xß + u (2) where, y is a nx1 vector of dependent variables, X, a nxm matrix of independent or explanatory variables with n>m, ß, a mx1 vector of unknown parameters and u is a nx1 vector of random errors. The problem is to find the unknown vector ß such that the estimated value of y be close to its observed value. A class of procedures which obtains these http://www.cribfb.com/journal/index.php/afbr Copyright © CC-BY-NC 2019, CRIBFB | AFBR www.cribfb.com/journal/index.php/afbr Australian Finance & Banking Review Vol. 3, No. 1; 2019 45 estimated values is Lp norm minimization criterion (see, Narula (1982)). In this class, ││u││p is minimized to find the ß vector, min S = min ││u││p = min ││y-Xß││p = ß ß ß n n m min [ │yi-xiß│p ]1/p = min [ │yi- ßjxij│ p ]1/p ==> ß i=1 ß i=1 j=1 n n m min │yi- xiß│p = min │yi- ßjxij│ p (3) ß i=1 ß i=1 j=1 where yi is the ith element of y and xi is the ith row of the matrix X. Any value of p [1, ] may be used to find ß in (3) (see, Money et al. (1978a), Rice (1933)), but each value of p is relevant for special types of error distributions. Many authors have investigated this problem (see, Barrodale (1968), Barr et al (1980a,b,c,81a,b), Money et al (1978b,82), Gonin and Money (1985a,b), Sposito and Hand (1980), Sposito and Hand and Skarpness (1983), Sposito (1987b)). However, justification of p comes from the following theorem (see, Kiountouzis (1971), Rice and White (1964), Hogan (1976), Taguchi (1974,78)). Theorem: If in model (2), X is nonstochastic and E(u)=0, E(uuT)= ²I, and u distributed with f(u)=h.exp(-k│u│p), where h and k are constants and p [1, ]; then the "best" ß with maximum likelihood properties is a vector which comes from minimization of (3). Certain values of p have particular importance (see, Box and Tiao (1962), Theil (1965), Anscombe (1967), Zeckhauser and Thompson (1970), Blattberg and Sargent (1971), Kadiyala (1972), Maddala (1977)). L norm minimization of (3) is called Tchebyshev or uniform norm minimization or minimum maximum deviations and has the maximum likelihood properties when u has a uniform probability distribution function. When p=2, we are confronted with the least squares method. In this case, if the errors’ distribution is normal, it is the best unbiased estimator (see, Anderson (1962), Theil (1971)). When p=1, we have L1‎‎ norm or Gershgorin norm minimization problem. It is also called least or minimum sum of absolute errors (MSAE, LSAE), minimum or least absolute deviations, errors, residuals, or values (MAD, MAE, MAR, MAV, LAD, LAE, LAR, LAV), L1‎‎ norm fit, approximation, regression or estimation. Harter (1974a,b,75a,b,c,76) monumental papers provide a chronology of works on nearly all the estimators which include L1‎‎ norm estimation too. A concise review of data analysis based on the L1‎‎ norm is presented by Dodge (1987), and a brief discussion is given by Gentle (1977) too. Narula and Wellington (1982) and Narula (1987) give a brief and concise presentation of L1‎‎ norm regression. Blattberg and Sargent (1971) show that if the errors of the regression follow the second law of Laplace (two-tailed exponential distribution) with probability density function f(u)=(1/2 ).exp(-│u│/ ) (4) where var(u)=2 ², then L1‎‎ norm minimization leads to maximum likelihood estimator. 3. Properties of the L1‎‎ norm estimation Similar to other criteria, the L1‎‎ norm estimation has its own properties, which are essential in computational and statistical viewpoints. The more important properties are as follows. 3.1 Invariance property An estimator ß^(y,X) of population parameter ß is invariant if, http://www.cribfb.com/journal/index.php/afbr Copyright © CC-BY-NC 2019, CRIBFB | AFBR www.cribfb.com/journal/index.php/afbr Australian Finance & Banking Review Vol. 3, No. 1; 2019 46 ß^( y,X) = ß^(y,X), [0, ) (5) Gentle and Sposito (1976), Koenker and Bassett (1978) have proved that the Lp norm estimator of ß is invariant when the regression model is linear. The Lp norm estimator is not invariant for general nonlinear models. The invariance property is the homogeneity of degree one of the ß^ solution function. 3.2 Transformation of variables If Rm, by transforming y to y+X the optimal value of ß^ will increase by , (see, Koenker and Bassett (1978)); ß^(y+X ,X) = ß^(y,X) + (6) If A is a mxm nonsingular matrix, the transformation of X to XA premultiplies optimal ß^ by the inverse of A (see, Taylor (1974), Koenker and Bassett (1978), Bassett and Koenker (1978)). ß^(y,XA) = A-1ß^(y,X) (7) 3.3 convexity of the objective function To show the convexity of S in (3), suppose m=1; the objective function (3) reduces to n n S = │yi - ß1xi1│ = Si (8) i=1 i=1 where Si=│yi-ß1xi1│. If we plot Si as a function of ß1, then we will have a broken line in Sxß1 plane, and its function value is zero at ßi1=yi/xi1. The slope of the half-lines to the left and right of ßi1 are -│xi1│ and │xi1│, respectively. So, Si's are all convex, and hence their sum S is also convex with slope at any ß1 equal to the sum of the slopes of the Si's at that value of ß1 (see, Karst (1958), Taylor (1974)). Consider now (3) when m=2, n n S = │yi - ß1xi1 – ß2xi2│ = Si (9) i=1 i=1 Where Si=│yi-ß1xi1-ß2xi2│. We may plot Si as a function of ß1 and ß2. Every Si is composed of two half-planes in Sxß1xß2 space that intersect in the ß1xß2 plane. Thus Si is convex downward which its minimum locates on the intersection of the two half-planes. Since Si's are all convex, their sum S surface is convex too. Extension to m independent variables is straightforward. In this case, each Si consists of two m dimensional half- hyperplanes in Sxß1x...xßm space intersecting in the ß1x...xßm hyperplane, and as before is convex in the opposite direction of the S axis. S, which is the sum of all these half-hyperplanes forms a polyhedron hypersurface which is convex too. 3.4 Zero residuals in the optimal solution L1‎‎ norm regression hyperplane always passes through r of thy n data points, where r is rank of the X matrix. Usually, X is of full rank, and thus r is equal to m. So, for the number of parameters, there exist zero residuals for the minimal solution of (3). This implies that L1‎‎ norm regression hyperplane must pass through m observation points (see, Karst (1959), Taylor (1974), Money et al. (1978), Appa and Smith (1973), Gentle and Sposito and Kennedy (1977)). This phenomenon is because of the polyhedron shape of the S. It is obvious that the minimum solution occurs on at least one of the corners of S, and the corners of S are the loci of changes in slopes of the polygonal http://www.cribfb.com/journal/index.php/afbr Copyright © CC-BY-NC 2019, CRIBFB | AFBR www.cribfb.com/journal/index.php/afbr Australian Finance & Banking Review Vol. 3, No. 1; 2019 47 hypersurface. Note that these corners and also edges of S will be above the intersections of m subset of the following n hyperplanes. m yi - ßjxij = 0 i {1,...,n} (10) j=1 Since each of these hyperplanes corresponds to a particular m subset of observations, there will be m observations that lie on the regression hyperplane (see, Taylor (1974)). 3.5 Optimality condition This condition is derived from the Kuhn-Tucker necessary condition of nonlinear programming and proved by Gonin and Monpy (1987b) and Charalambous (1979). Define A={i│yi-xiß *=0} and I={i│yi-xiß *╪0}; in linear L1‎‎ norm regression, a necessary and sufficient condition for ß* to be a global L1‎‎ norm solution is the existence of multipliers i [-1,1] such that: ixi + sgn(yi-xiß *)xi = 0 (11) i A i I (see also, El-Attar and Vidyasagar and Dutta (1976). Appa and Smith (1973) showed that this solution is a hyperplane such that: │n+ - n-│ m (12) where n+ and n- are the number of observations above and below the regression hyperplane, respectively. 3.6 Unique and non-unique solutions Since S is a convex polyhedron hypersurface, it always has a minimum. This solution is often unique. Sometimes the shape of S is such that a line or a closed polygon or polyhedron or hyperpolyhedron segment of S is parallel to ß1x...xßm hyperplane. On this case the L1‎‎ norm regression parameters are not unique and infinite points of the mentioned hyperpolyhedron are all solutions (see, Moroney (1961), Sielken and Hartley (1973), Taylor (1974), Farebrother (1985), Sposito (1982), Harter (1977)). 3.7 Interior and sensitivity analysis Narula and Wellington (1985) Showed that the L1‎‎ norm estimates might not be affected by certain data points. Thus deleting those points does not change the estimated values of the regression parameters. In another discussion, they called sensitivity of L1‎‎ norm estimates, determined the amounts by which the value of response variable yi can be changed before the parameters estimates are affected. Specifically, if the value of yi increases or decreases without changing the sign of ui, the solution of the parameters will not change (see, Gauss (1809), Farebrother (1987b)). For the topology of L1‎‎ norm approximation and its properties see Kripke and Rivlin (1965), Vajda (1987), Hromadka II et al. (1987). Other properties of the L1‎‎ norm regression are discussed by Gentle and Kennedy and Sposito (1976,77), Assouad (1977), Sposito and Kennedy and Gentle (1980), Bassett (1987,88a,b). 4. Chronology and historical development (1632-1928) The origin of L1‎‎ norm estimation may be traced back to Galilei (1632). In determining the position of a newly discovered star, he proposed the least possible correction in order to obtain a reliable result (see, Ronchetti (1987) for some direct quotations). Boscovich (1757) for the first time, formulated and applied the minimum sum of absolute errors for obtaining the best fitting line given three or more pairs of observations for a simple two-variable regression model. He also restricts the line to pass through the means of the observation points. That is, http://www.cribfb.com/journal/index.php/afbr Copyright © CC-BY-NC 2019, CRIBFB | AFBR www.cribfb.com/journal/index.php/afbr Australian Finance & Banking Review Vol. 3, No. 1; 2019 48 n min : │yi-ß0-ß1xi1│ ß0,ß1 i=1 n (13) s.to: (yi-ß0-ß1xi1)=0 i=1 Boscovich (1760) gives a simple geometrical solution to his previous suggestion. This paper has been discussed by Eisenhart (1961) and Sheynin (1973). In a manuscript, Boscovich poses the problem to Simpson and Simpson gives an analytical solution to the problem (see, Stigler (1984)). Laplace (1773) provides an algebraic formulation of an algorithm for the L1‎‎ norm regression line, which passes through the centroid of observations. In Laplace (1779), the extension of L1‎‎ norm regression to observations with different weights has also been discussed. Prony (1804) gives a geometric interpretation of Laplace's (1779) method and compares it with other methods through an example. Svanberg (1805) applies Laplace's method in determining a meridian arc, and Von Lindenau (1806) uses this method in determination of the elliptic meridian. Gauss (1809) suggests the minimization of the sum of absolute errors without constraint. He concludes that this criterion necessarily sets m of the residuals equal to zero, where m is the number of parameters, and further, the solution obtained by this method is not changed if the value of the dependent variable is increased or decreased without changing the sign of the residual. This conclusion is recently discussed by Narula and Wellington (1985) which explained in the previous section under the subject of interior and sensitivity analysis. He also noted that Boscovich or Laplace estimators which minimize the sum of absolute residuals with zero-sum of residuals constraint, necessarily set m-1 of the residuals equal to zero (see, Stigler (1981), Farebrother (1987b)). Mathieu (1816) used Laplace's method to compute the eccentricity of the earth. Van Beeck-Calkoen (1816) advocates the using of the least absolute values criterion in fitting curvilinear equation obtained by using powers of the independent variable. Laplace (1818) adapted Boscovich's criterion again and gave an algebraic procedure (see, Farebrother (1987b)). Let x1* and y* be the means of xi1 and yi then, ß0 = y* - ß1x1* (14) Value of ß1 is found by, n min: S = │yi~ - ß1xi1~│ (15) ß1 i=1 where, xi1~ and yi~ are deviations of xi1 and yi from these means respectively. By rearranging the observations in descending order of yi~/xi1~ values, Laplace notes that S is infinite when ß1 is infinite and decreases as ß1 is reduced. ß1 reaches the critical value yt~/xt1~ when it again begins to increase. This critical value of ß1 is determined when, t-1 n t │xi1~│ < ½ │xi1~│ │xi1~│ (16) i=1 i=1 i=1 This procedure to find a1 is called weighted median and has been used in many other algorithms such as Rhodes (1930), Singleton (1940), Karst (1958), Bloomfield and Steiger (1980), Bidabad (1987a,b,88a,b) later. Bidabad (1987a,b,88a,b) derives the condition (16) via discrete differentiation method. http://www.cribfb.com/journal/index.php/afbr Copyright © CC-BY-NC 2019, CRIBFB | AFBR www.cribfb.com/journal/index.php/afbr Australian Finance & Banking Review Vol. 3, No. 1; 2019 49 Fourier (1824) formulates least absolute residuals regression as what we would now call linear programming; that is the minimization of a linear objective function subject to linear inequality constraints. Edgeworth (1883) presents a philosophical discussion on differences between minimizing mean square errors and mean absolute errors. Edgeworth (1887a,b) proposed a simple method for choosing the regression parameters. By fixing m-1 of the parameters, he used Laplace's procedure to determine the optimal value of the remaining parameter. Repeating this operation for a range of values for m-1 fixed parameters, he obtained a set of results for each of m possible choices of the free parameters. Edgeworth drops the restriction of passing through the centroid of data. Turner (1887) discusses the problem of non-unique solutions under the least absolute error criterion as a graphical variant of Edgeworth (1887a) as a possible drawback to the method. Edgeworth (1888) replies to Turner's criticism by proposing a second method for choosing the two parameters of least absolute error regression of a simple linear model which makes no use of the median loci of his first method. Edgeworth, in this paper, followed Turner's suggestion for graphical analysis of steps to reach the minimum solution. Before referring to double median method of Edgeworth (1923), it should be noted that Bowley (1902) completes the Edgeworth's (1902) paper by a variant of double median method which presented after him by Edgeworth (1923). This variant ignores the weights attached to errors. Edgeworth (1923) discussed the more general problem of estimating the simple linear regression parameters by minimizing the weighted sum of the absolute residuals. He restates the rationale for the method and illustrates its usage through several examples. He also considers the nonunique solution problem. His contribution is called double median method. Estienne (1926-28) proposes replacing the classical theory of errors of data based on least squares with what he calls a rational theory based on the least absolute residual procedure. Bowley (1928) summarizes the Edgeworth's contributions to mathematical statistics, which includes his work on L1‎‎ norm regression. Dufton (1928) also gives a graphical method of fitting a regression line. Farebrother (1987b) summarizes the important contributions to L1‎‎ norm regression for the period of 1793- 1930. For more references see also Crocker (1969), Harter (1974a,b,75a,b,c,76), Dielman (1984). Up to 1928, all algorithms had been proposed for simple linear regression. Though some of them use algebraic propositions, are not so organized to handle multiple L1‎‎ norm regression problem. In the next section, we will discuss the more elaborated computational methods for simple and multiple L1‎‎ norm regressions not in a chronological sense; because many digressions have occurred. We may denote the period of after 1928 the time of modern algorithms in the subject of L1‎‎ norm regression. 5. Computational algorithms Although a closed form of the solution of L1‎‎ norm regression has not been derived yet, many algorithms have been proposed to minimize its objective function (see, Cheney (1966), Chambers (1977), Dielman and Pfaffenberger (1982,84)). Generally, we can classify all L1‎‎ norm algorithms in three major categories as, direct descent algorithms, simplex type algorithms, and other algorithms which will be discussed in the following sections sequentially. 5.1 Direct descent algorithms The essence of the algorithms which fall within this category is finding a steep path to descend down the polyhedron of the L1‎‎ norm regression objective function. Although the Laplace's method (explained hereinbefore) is a special type of direct descent algorithms; the origin of this procedure in the area of L1‎‎ norm can be traced back to the algorithms of Edgeworth which were explained in the previous section. Rhodes (1930) found Edgeworth's graphical solution laborious; therefore, he suggested an alternative method for the general linear model, which may be summarized as follows (see, Farebrother (1987b)). Suppose, we have n equations with m