




































Copyright © CC-BY-NC 2019, BJMSR 

 

Bangladesh Journal of Multidisciplinary Scientific Research 
Vol. 1, No. 1 

April-June; 2019 
Published by Centre for Research on Islamic Banking & Finance and Business 

 

     50 

L1 Norm Based Computational Algorithms 
 
 

Bijan Bidabad  
B.A., M.Sc., Ph.D., Post-Doc.  

Professor  
Economics and Chief Islamic Banking Advisor  

Bank Melli, Iran 
E-mail:bijan@bidabad.com 

 
 
Abstract 
This paper gives a rather general review of the L1 norm algorithms. The chronology and historical development of the L1 norm 
estimation theory for the period of 1632-1928 will be surveyed and the algorithms belonging to the after 1928 period will be 
categorized into three main classes of direct descent, simplex type, and other algorithms. 
 
Keywords: L1 norm, Regression, Algorithm, Computer  
 
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,  
                n 

     min :  Σ │yi-ß0-ß1xi1│ 
     ß0,ß1 i=1 
              n                                                                                                                                 (1) 

     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 (1793) provides an algebraic formulation of an algorithm for the L1 norm regression line, which passes through the 
centroid of observations. In Laplace (1799), an extension of L1 norm regression to observations with different weights has also 
been discussed. Prony (1804) gives a geometric interpretation of Laplace's (1799) 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). 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 1x  and 

y be the means of xi1 and yi then,  

     ß0 = y - ß1 1x                                                                                                                       (2) 

Value of ß1 is found by, 
                    n 

mailto:bijan@bidabad.com


Copyright © CC-BY-NC 2019, BJMSR 

 

 www.cribfb.com/journal/index.php/BJMSR       Bangladesh Journal of Multidisciplinary Scientific Research         Vol. 1, No. 1; 2019 

51 
               

     min: S = Σ │yi~ - ß1xi1~│                                                                                                   (3) 
      ß1         i=1 
where xi1~ and yi~ are deviations of xi1 and yi from their 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~│                                                                                (4) 
     i=1                  i=1               i=1 
This procedure to find ß1 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,89a,b) later. Bidabad 
(1987a,88a,89a,b) derives the condition (4) via discrete differentiation method. 
    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 nonunique 
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. 
 
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, 
i)   Direct descent algorithms 
ii)  Simplex type algorithms 
iii) Other algorithms 
which will be discussed in the following sections sequentially. 
 
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 a general 
linear model, which may be summarized as follows (see, Farebrother (1987b)). Suppose, we have n equations with m<n 
unknown parameters. To find L1 norm solution of this overdetermined system of equations; 



Copyright © CC-BY-NC 2019, BJMSR 

 

 www.cribfb.com/journal/index.php/BJMSR       Bangladesh Journal of Multidisciplinary Scientific Research         Vol. 1, No. 1; 2019 

52 
               

Step 1) Select m-1 equations arbitrarily. 
Step 2) Solve these equations for m-1 parameters. 
Step 3) Estimate the remaining mth parameter by Laplace's method. 
Step 4) Recognize the resulting equation in step 3 and add it to the m-1 equations in step 1. 
Step 5) If the set of m equations has recurred m times then stop; otherwise discard the oldest 
            equation and go to step 2. 
Rhodes (1930) explained his algorithm by an example and did not give any proof for convergence. Bruen (1938) reviews the L 1 
norm regression methods presented by earlier authors. He also compares L1, L2, and L∞ norms regressions. 
    Singleton (1940) applied Cauchy's steepest descent method (see, Panik (1976)) for the general linear L1 norm regression. In 
this paper, a geometrical interpretation of gradient on L1 norm polyhedron and some theorems about existence and uniqueness 
of solution and convexity property all were given. Although the paper has not been clearly written, the following steps 
summarize his algorithm. 
Step 1) Select a point ßj(0), j=1,...,m arbitrarily. 
Step 2) Determine the gradient, 
                  n 

        gj
(0)=-Σ sgn(ui^│ßj

(0),j=1,...,m)xij. 
                i=1                  m                         m 

Step 3) Compute, wi
(0)= Σ xijgj

(0), zi
(0)=yi- Σ xijßj

(0). 
                                      j=1                       j=1                     n 

Step 4) Determine the value of t(0) as weighted median of Σ │wit-zi│ by Laplace's method.  
                                                                                           i=1 
            The t value is the length of movement along the direction of the gradient. 
Step 5) Compute ßj

(1)=ßj
(0)+gj

(0)t(0). 
Step 6) Test the optimality condition. Singleton gives a condition to stop, but it is not quite clear. In this step, any other 

criterion relevant to L1 norm function may be displaced. 
Step 7) This step is not well defined by Singleton. In this phase, he tries to choose the best gradient among usable gradients. 

Without this step, algorithm is still operational, because, the steps are all standards of Cauchy steepest descent method 
and instead of choosing the best gradient, we can proceed by going to step 2. 

    Bejar (1956,57) focuses on consideration of residuals rather than on the vector of parameters. He puts forth a procedure with 
the essence of Rhodes (1930). However, he is concerned with two and three parameter linear models. 
    Karst (1958) gives an expository paper for one and two parameter regression models. In his paper, Karst without referring to 
previous literature actually reaches to the Laplace proposition to solve the one parameter restricted linear model, and for the 
two-parameter model, he proposed an algorithm similar to that of Rhodes (1930). His viewpoint is both geometrical and 
algebraic, and no proof of convergence for his iterative method is offered. Sadovski (1974) uses a simple "bubble sort" 
procedure and implements Karst algorithm in Fortran. Sposito (1976) pointed out that the Sadovski's program may not 
converge in general. Sposito and Smith (1976) offered another algorithm to remove this problem. Farebrother (1987c) recodes 
Sadovski's implementation in Pascal language with some improvement such as applying "straight insert sort".  
    Usow (1967b) presents an algorithm for L1 norm approximation for discrete data and proves that it converges in a finite 
number of steps. A similar algorithm on L1 norm approximation for continuous data is given by Usow (1967a). Given the 
function f(x) defined on a finite subset X={x1,...,xn} of an interval on the real line and also linearly independent continuous 

functions Φ1(x),...,Φm(x) where m<n; consider the "polynomial" 
             m 

L(ß,x)=ΣßjΦj(x). In Usow (1967b) the following function is to be minimized: 
           j=1 
                        n                             n    m 

     min: S(ß)= Σ │L(ß,xi)-f(xi)│= Σ │Σ ßjΦj(xi)-f(xi)│                                                          (5) 
      ß              i=1                          i=1 j=1 

where ß is a vector of size m. The above form is general; if we let f(xi)≡yi and Φj(xi)≡xij the function S becomes the standard 
linear regression objective function. Now let the set K be:  

     K={(ß,d)│(ß,d)εEm+1,S≤d}                                                                                                (6) 
Then K is a convex polytope, the vertices of which occur only when L(ß,x)-f(x) is zero at m or more points of X. The Usow's 
algorithm is to descend on K from vertex to vertex along connecting edges of the polytope in such a way that certain 
intermediate vertices are by-passed. This descent continues until the lowest vertex (ß*,d*) is reached. To clarify the algorithm 
assume that we are at the vertex (ßk,dk) on K and the polynomial L(ßk,x) interpolates m points of X denoted by Uk={u1

k,...,um
k}. 

Thus, 
                    m 



Copyright © CC-BY-NC 2019, BJMSR 

 

 www.cribfb.com/journal/index.php/BJMSR       Bangladesh Journal of Multidisciplinary Scientific Research         Vol. 1, No. 1; 2019 

53 
               

     L(Fk,x)= Σ f(ui
k)πi(x)                                                                                                          (7) 

                   i=1 
where, 
     Fk=(f(u1

k),...,f(um
k))                                                                                                            (8) 

                m 

     πi(x)= Σ aj
iΦj(x)         i=1,...,m                                                                                            (9) 

               j=1 

The m coefficients aji are calculated as follows. Form the matrix (Φj(ui
k)) for i,j=1,...,m. Let π(x) and Φ(x) be two mx1 vectors 

for any xεX whose elements are πi(x) and Φi(x) for i=1,...,m respectively. Hence, we can derive, 

     π(x)=[(Φj(ui
k))T]-1Φ(x)                                                                                                    (10) 

aj
i's are the elements of the ith row of the matrix [(Φj(ui

k))T]-1, where the superscript T denotes transposition. Let ei be a zero 

vector of size m where its ith element is equal to one. Then if for some δ, S(Fk-δei)<S(Fk), there is a Tj such that Tjδ>0 and 
S(Fk-Tjei)<S(Fk). Also, 
     S(Fk - Tjei) = min {S(Fk - tei)}                                                                                          (11) 
                             t 

and ((Fk-Tjei),S(Fk-Tjei)) is a vertex. On the other hand, a point ui
kεUk may be replaced by a point ui

k+1ε{X-Uk} such that the 
polynomial L(ßki,x) interpolating Ui

k={u1
k,...,ui

k+1,...,um
k} and S(ßki)<S(ßk). S(ßKi) is the minimum of all norms obtained if ui

k 
were replaced by the different points of the set {X-Uk}, as indicated by the above relation.  
    In going from vertex (ßk,S(ßk)) to (ßki,S(ßki)), one or more vertices on K might have been by-passed. The nearest vertex to 
(Fk,S(Fk)) and below it on the edge parallel to the ith parameter space coordinate axis, say the vertex ((Fk-trei),S(Fk-trei)), is 
obtained from:  

                           │L(Fk,xs)-f(xs)│ 

     │tr│ = min {──────────}           xsε{X-Uk}                                                           (12) 

                   s             │πi(xs)│ 
The point xr is characterized by, 

     sgn[L(Fk,xl) - f(xl)] = sgn[L(Fki,xl) - f(xl)],  xlε{X-Uk-xr}                                                (13) 

Now, if there is not any δ such that S(Fk-δei)<S(Fk), then S(Fk) could not be reduced by moving on K along the edge parallel to 
the ith parameter space coordinate axis and ui

k should not be replaced by another point from the set {X-Uk}. This iteration is 
repeated m times, once for each point in Uk in succession. The whole cycle is then repeated a finite number of times until the 
solution (ß*,d*) is reached (see also, Abdelmalek (1974)).  
    Relation of this algorithm with the simplex method has been discussed by Abdelmalek (1974). He shows that Usow's 
algorithm is completely equivalent to a dual simplex algorithm applied to a linear programming model with nonnegative 
bounded variables, and one iteration in the former is equivalent to one or more iterations in the latter. Bloomfield and Steiger 
(1980) devise an efficient algorithm based on the proposition of Usow explained above. 
    Sharpe (1971) by applying the L1 norm regression to the portfolio and its rate of return, gives an algorithm for the two 
parameters linear regression model. He argues that for the simple model with the objective function,  
 
            n 

     S = Σ │yi - (ß0 + ß1xi1)│                                                                                                   (14) 
          i=1 
it must be possible to assign half of the points above and half below the regression line. With any given value of ß~

1, we can 
derive ß0 as the median of ß0i=yi-ß~

1xi1. Now segregate the points, such that:  
            nabove                       nbelow 

     S=  Σ   [yi-(ß0+ß1xi1)] -  Σ   [yi-(ß0+ß1xi1)]                                                                        (15) 
            iabove                        ibelow 
Rearranging the terms and note that nabove=nbelow gives, 
     S = k1 + k2ß1                                                                                                                      (16) 
where, 
              nabove      nbelow 

     k1 =  Σ   yi -    Σ    yi                                                                                                          (17) 
              iabove     ibelow 
               nabove       nbelow 

     k2 = -Σ   xi1  +  Σ    xi1                                                                                                       (18) 
               iabove     ibelow 



Copyright © CC-BY-NC 2019, BJMSR 

 

 www.cribfb.com/journal/index.php/BJMSR       Bangladesh Journal of Multidisciplinary Scientific Research         Vol. 1, No. 1; 2019 

54 
               

The overall solution strategy may now be formulated as follows. Any value of ß1 may be chosen at the outset. By computing ß0i, 
segregate the points above and below the line. By using equations (16) through (18), the corresponding segment of the S and ß 1 
relation is calculated. Sign of k2 indicates the appropriate direction for the next iteration. If k2 is positive, only smaller values of 
ß1 need be considered. If k2 is negative, only larger values of ß1 need be considered. If k2 is zero, the initial value of ß1 is a 
solution. 
    When a borderline (median of ß0i) and the direction of search is determined, the nearest intersection of the present borderline 
with another should be found. The calculations can be reduced by comparing slopes (x i1 values) to determine whether or not two 
lines intersect in the region of interest (thus avoiding needless division operations). The new values k1 and k2 must be computed. 
If this alteration causes k2 to change sign, the solution has been obtained, and the algorithm stops.  
    Rao and Srinivasan (1972) interpret Sharpe's procedure as the solution of parametric dual linear programming formulation of 
the problem. They give an alternate and about equally efficient procedure for solving the same problem. Brown (1980) gives a 
distinct but similar approach to those of Edgeworth (1923) and Sharpe (1971). He emphasizes on the median properties of the 
estimator. The similarity comes from the graphical approach of the three authors. Kawara (1979) also develops a graphical 
method for the simple regression model. 
    Bartels and Conn and Sinclair (1978) apply the method of Conn (1976) to the L1 norm solution of an overdetermined linear 
system. Their approach is a minimization technique for piecewise differentiable functions. The algorithm may be reduced as 
follows. 
Step 0) Select an arbitrary point ß. 

Step 1) i) Identify the index set, I = {i1,...,im} = {i│xi
Tß-yi=0}. 

         Let Ith row of X be XI=[xi ,...,xi ] where ijεI and the nullspace N=N(XI
T)={δ│xiδ=0,iεI}. 

                                                 1           m 

         Let the orthogonal projector onto N be denoted by PN. 

         ii) Let Ic be the complement of the set I. Compute the vector h = Σ sgn(xi
Tß-yi)xi. 

                                                                                                               iεIc 
        iii) Compute p=-PNh that is the projection of -h onto the nullspace of XI, as long as this projection is non zero. If p≠0, let 

g=h and go to step 2. 
 
 
                                                                   m 

        iv) compute w according to h=XIw= Σ wjxi ,ijεI. 
                                                                  j=1      j 

        v) if │wj│≤1 for all j=1,...,m stop. In this case, ß is optimal. 

        vi) Find ij0εI such that │wj0│>1. 
        vii) Change I to I-{ij0} and make corresponding changes to XI and N. Compute  
        p=-sgn(wj0)PNxi  and let g=h-sgn(wj0)xi  . 
                                                   j0                                                     j0 

Step 2) Determine A={αl│lεIc,αl=(xl
Tß-yl)/xl

Tp,αl>0} elements and order them such that  

        0<αl <αl <...<αl .  Let τ=1. 
                  1       2               t 

Step 3) If pTg≥2sgn(xl Tß-yl )pTxl , then go to step 5. 
                                    τ            τ            τ 

Step 4) Change g to g-2sgn(xl ß-yl )xl and τ to τ+1 and go to step 3. 

                                                τ          τ      τ 

Step 5) Replace ß by ß+αl p and go to step 1. 

                                          τ 

This algorithm has also been modified for the case of degeneracy (see also, Bartels and Conn and Sinclair (1976)). 
    Bartels and Conn (1977) showed that how L1 norm, restricted L1 norm, L∞ norm regressions, and general linear programming 
can all be easily expressed as a piecewise linear minimization problem. Let U and v be of sizes pxm and px1, respectively. 
Consider the function: 

     Φ(ß) = hTß + Σi│yi-xi
Tß│ + Σjmax(0,vj-uj

Tß)                                                                  (19) 
Where yi-xi

Tß represents the ith element of the residual vector y-Xß and vj-uj
Tß represents the jth element of the residual vector v-

Uß. To minimize Φ with respect to ß the following steps should be taken. 
Step 0) Start with arbitrary ß(k). 

Step 1) Find δ(k) so that Φ(ß(k)+Θδ(k))≤Φ(ß(k)) for all Θ>0 small enough. 

Step 2) Choose Θ(k)≥0 to obtain the largest possible decrease in Φ. 



Copyright © CC-BY-NC 2019, BJMSR 

 

 www.cribfb.com/journal/index.php/BJMSR       Bangladesh Journal of Multidisciplinary Scientific Research         Vol. 1, No. 1; 2019 

55 
               

Step 3) Let ß(k+1)=ß(k)+Θ(k)δ(k). 

When h=0 and the sum on j is vacuous in the Φ(ß) function (19), the resulting simplification of the above steps corresponds 
precisely to the algorithm proposed by Bartels and Conn and Sinclair (1978). Other related problems can be solved by 

modification of Φ(ß) by a quantity µ to obtain a parameterized family of piecewise linear functions Φµ(ß), and take following 
steps to find the minimum solution.  
Step 0) Set µ>0; select any ß=ß(0). 

Step 1) Minimize Φµ(ß) with respect to ß according to the above procedure. 
Step 2) Stop if a prescribed terminating condition on ß is met; otherwise set µ=µ/10 and go to step 1. 
The contribution of this paper is putting a wide class of problems in the mould of two algorithms mentioned above. The 
techniques are easily extended to the models with norm restrictions. 
    Bloomfield and Steiger (1980) proposed a descent method for the L1 norm multiple regression. Their algorithm is also 
explained in Bloomfield and Steiger (1983). In some steps, this algorithm is related to that of Singleton (1940) and Usow 
(1967b). The basis of this method is to search for a set of m observations which locate on the optimal L1 norm regression. This 
set is found iteratively by successive improvement. In each iteration, one point from the current set is identified as a good 
prospect for deletion. This point is then replaced by the best alternative. The novel features of this method are in an efficient 
procedure for finding the optimal replacement and a heuristic method for identifying the point to be deleted from the pivot. 
    Denote the x1

T,...,xm
T as rows of the independent variables design matrix which correspond to the current set of points 1,...,m; 

and xm
T for replacement. Set, 

     yi = xi
Tß      i=1,...,m-1                                                                                                      (20) 

Redefine ß as, 

     ß = ß0 + tδ                                                                                                                         (21) 

Where ß0 is an arbitrary member of the set and δ vector obeys the following system, 

     xi
Tδ = 0       i=1,...,m-1                                                                                                      (22) 

Given this set of points, the optimum value of S may be found by minimizing the following expression with respect to the scalar 
t. 
      n 

      Σ │yi - xi
T(ß0 + tδ)│                                                                                                        (23) 

     i=1 
Rearranging the terms leads to: 
      n 

      Σ │wi││ri - t│                                                                                                                (24) 
     i=1 

where wi=xi
Tδ and ri=(yi-xi

Tß0)/(xi
Tδ). Value of t may be found by Laplace's weighted median method. Bloomfield and Steiger 

propose a weighted modification of partial quicksort procedure of Chamber (1971) to find the t value efficiently. The data 
point corresponding to the weighted median replaces xk

T. By using the yi from (20) and the additional mth equation, the value of 

ß0 is computed. Vector δ is determined up to scalar multiples of (22). The new set of parameters values are computed by (21). 
Now a point should be deleted. Bloomfield and Steiger do not give an assured way to identify this point. They propose a 
heuristic method based on gradients and use the following quantity: 

              │ Σ  wi - Σ   wi│ - Σ   wi 

              │i:ri<0  i:ri>0  │  i:ri=0 

     ρ   = ──────────────                                                                                         (25) 
                       n 

                       Σ wi 
                      i=1 

Once ρ is calculated for each candidate point for deletion, and delete the point for which ρ is largest.  
    To start the algorithm, any set of m rows of X may be chosen, with the appropriate ß0. Add variables stepwise until a fit for 
ß0 and the corresponding set of m points are derived. At each intermediate step, the fit involves k variables where 0≤k<m and 
the corresponding set of k data points with zero residuals. Improving the fit by entering a new variable, thus increasing k to k+1. 

At each stage, it is the measure ρ that determines whether we set up to a larger model or improve the current one without setting 

up. Suppose at the current stage we are dealing with k variables. If the value of ρ is largest at the variable p, which p is not in the 
set of k variables above, k is increased to k+1 by entering the variable p and improving the fit with k+1 variables. If p is in the 

current set of k variables and ρ is in the largest value, the corresponding point to p is deleted and replaced in a manner described 
before. In this paper relationship of this algorithm to linear programming is also discussed.  



Copyright © CC-BY-NC 2019, BJMSR 

 

 www.cribfb.com/journal/index.php/BJMSR       Bangladesh Journal of Multidisciplinary Scientific Research         Vol. 1, No. 1; 2019 

56 
               

    Seneta and Steiger (1984) proposed an algorithm for the L1 norm solution of a slightly overdetermined system of equations. 
Their proposition is based on the above algorithm of Bloomfield and Steiger. It is more efficient than the former if m is near n. 

Given (xi,yi)εRm+1, i=1,...,n and k=n-m, X=(x1,...,xn)T, their algorithm may be described as follows: 
 
Step 1) Renumber the rows of (X|y) such that XN the bottom m rows of X, is invertible. Solve the k linear equations system 

XN
TN=-XT

T for N, where XT denotes the top k rows of X. 
Step 2) Let D=(A|c) where A=(I|N) is of size kxn and c=Ay. 

Step 3) Let σ =(1,...,k), σc=(k+1,...,n). 

Step 4) Set rσ(i)=bi for i=1,...,k; rσc(i)=0 for i=k+1,...,n. 

Step 5) Let I={i│1≤i≤k,ci=0}. 
Step 6) Do loop for j=1 to m: 

            Let vi=Diσc(j) for j=1,...,k. Let M={i│sgn(ci)≠sgn(vi)} and J={1,...,k}\N\I. 

                         │ΣM│vi│-ΣJ│vi││-1-ΣI│vi│ 

            Let ßj = ────────────────── 
                                                  k 

                                            1+ Σ │vi│ 
                                                i=1 
           End loop. 
Step 7) Set S={1,...,m}. 
Step 8) Choose ßq as max{ßj}. 
                                     s 

Step 9) If ßq>0 go to step 11. 
                  k 

Step 10) If П rσ(i)=0, the problem is degenerate;  and  stop. 
                 i=1                              m 

             Otherwise, solve, yσc(i)= Σ xσc(i)jΘj for Θ, stop. 
                                                   j=1 

Step 11) Let vi=Di,σc(q) for i=1,...,k. 

Step 12) Compute t^=weighted median of c1/v1,...,ck/vk,0 with weights │v1│,...,│vk│,1. 
Step 13) If t^=cp/vp≠0 go to step 16. 
Step 14) Let S=S\{q}; if S=Ø go to step 10. 
Step 15) Go to step 8. 

Step 16) Divide row p of D by Dpσc(q). 

Step 17) For i≠p in D, let (row i)=(row i)-(row p)*Diσc(q). 

Step 18) Commute for pair σ(p) and σc(q). 

Step 19) Let rσ(i)=bi for i=1,...,k; and set rσc(q)=0. 
Step 20) Go to step 5. 
    Seneta (1983) reviews the iterative use of weighted median to estimate the parameters vector in the classical linear model 
when the fitting criterion is L1 norm and also Cauchy criterion.  
    Wesolowsky (1981) presents an algorithm for multiple L1 norm regression based on the notion of edge descent along the 
polyhedron of the objective function. This algorithm is closely related to those of Rhodes (1930) and Bartels and Conn and 
Sinclair (1978) which explained before. Consider the multiple linear regression as before. Select a set of m points (xj1

I,...,xjm
I,yj

I). 
The following system of equations can be solved for a unique set of coefficients. 
             m 

     yj
I -  Σ ßhxjh

I = 0,   j=1,...,m                                           (26) 
           h=1 
An edge is formed of any subset J consisting of m-1 equations. To minimize along an edge, set: 
 

1 11 1 1 11 1

1 1 1,1 1,

, , ,

I I I J J J

m m

I I I J J J

m m mm m m m m

Y x x Y x x

Y x x Y x x  

       
      

         
      

       

I I J JY X Y X               (27) 

 



Copyright © CC-BY-NC 2019, BJMSR 

 

 www.cribfb.com/journal/index.php/BJMSR       Bangladesh Journal of Multidisciplinary Scientific Research         Vol. 1, No. 1; 2019 

57 
               

Let ßp be formed from ß by deleting ßp and let xp be the pth column in XJ and let XpJ be formed by removing xp from XJ.  

Then for a given ßp it can be shown that, 

ßp = XpJ
-1YJ - ßpXpJ

-1xp                                                               (28) 

Let the elements of ßp be 

ßq = rq - sqßp   for q ≠ m                                                                                              (29) 

Substitution in the L1 norm objective function gives, 

1

1
1

2

min :

m

q iq ip

q p

m

i q iqn m
q p

ip q iq m
i q p

ip q iq
qs x x

y r x

x s x
x s x







 



 

  
 



                                                    (30) 

Now the following steps should be taken. 

Step 1) Set k=1, l=0, choose initial values for ß1,...,ßm. Least squares values are one possibility. Let I(1)={j1
(1),...,jm

(1)} be a 

set of m data points chosen in sequence as follows. The point with the smallest squared residual is chosen each time 

subject to the condition that XI
(1) is nonsingular. Find ßq

(1) for q=1,...,m, by solving YI
(1)=XI

(1)ß. Set 

I(k)=(j1
(k),...,jm

(k)); J=(j2
(k),...,jm

(k)). 

Step 2) Set k=k+1. Obtain ßp for the smallest p from (30) by using the weighted median procedure. Set ßp
(k)=ßp and let i be 

the index which defines the lower weighted median ßp in (30) for the lowest possible p. 

Step 3) a) If ßp
(k)-ßp

(k-1) = 0 and if  l>m, go  to  step 4.  Otherwise, set I=(j2
(k-1),...,jm

(k-1),i) and l=l+1, ßq
(k)=ßq

(k-1) 

for all q and go to step 2. 

b) If ßp
(k)-ßp

(k-1)≠0, set l=0. Calculate ßq
(k) for q≠p from (29). Set I(k)=(j2

(k-1),...,jm
(k-1),i) and go to step 2. 

Step 4) Calculate ßq
(k) for q≠p from (29); set ß*=ß(k) and stop. 

In this paper, Wesolowsky also discusses the problem of multicolinearity and gives an appropriate solution. 

Josvanger and Sposito (1983) modify Wesolowsky's algorithm for the two-parameter simple linear regression model. 

The modification is an alternative way to order observations instead of sorting all of them to find the necessary weighted median 

value. Suppose the problem has been reduced to a weighted median problem. They place smaller values of factors to be sorted 

with corresponding weights below ß1
(k-1) and larger or equal values above it, then recheck the inequalities (4) of the weighted 

median. If the inequalities do not satisfy, then an appropriate adjustment is made. In particular, if the right-hand side is overly 

weighted, then the weight corresponding to the smallest sorting factor is transferred to the left-hand side, and the check is made 

again. A computer program for this algorithm is also given by the authors. 

"Generalized gradient" method introduced by Clarke (see, Clarke (1983)) is a general procedure for nonsmooth 

optimization functions and problems (see, Osborne and Pruess and Womersley (1986)). A subclass of this method is called 

"reduced gradient" explained by Osborne (1985) is a general algorithm which contains linear programming, piecewise linear 

optimization problems, and polyhedral convex function optimization algorithms inside. The reduced gradient algorithm is a 



Copyright © CC-BY-NC 2019, BJMSR 

 

 www.cribfb.com/journal/index.php/BJMSR       Bangladesh Journal of Multidisciplinary Scientific Research         Vol. 1, No. 1; 2019 

58 
               

special case of descent method, which possesses two important characteristics. Identify direction and taking a step in this 

direction to reduce the function value (see also, Anderson and Osborne (1976), Osborne and Watson (1985) Osborne 

(1985,87)). The algorithms of Bartels and Conn and Sinclair (1978), Armstrong and Frome and Kung (1979), Bloomfield and 

Steiger (1980) are all special cases of reduced gradient method. 

Imai and Kato and Yamamoto (1987) present a linear time algorithm for computing the two-parameter L1 norm linear 

regression by applying the pruning technique. Since the optimal solution in the ß0xß1 plane lies at the intersection of data lines, 

so, at each step, a set of data lines which does not determine the optimum solution are discarded. In this paper algebraic 

explanation of the problem is also offered. 

Pilibossian (1987) also gives an algorithm similar to Karst (1958) for the simple two-parameter linear L1 norm 

regression. Bidabad (1987a,b,88a,b) proposed descent methods for the simple and multiple L1 norm regressions. These 

algorithms, with many improvements, will be discussed by Bidabad (1989a,b). Bidabad (1989a,b) introduces four descent 

algorithms which two of them are for simple, and two others are for multiple regression models.His first algorithm is crude and 

tries to check many points to find the optimal solution. The second algorithm is more efficient. Algorithm three is a partial 

descent procedure for the general linear model. Algorithm four is a full descent method which has many proved properties. He 

also proves the convergence of all the above four algorithms. 

Simplex type algorithms 
The essence of linear programming in solving L1 norm problem may be found in the work of Edgeworth (1888). Harris 

(1950) suggested that the L1 norm estimation problem is connected with linear programming. Charnes and Cooper and 

Ferguson (1955) formulated the problem as a linear programming model. This article is the first known to use linear 

programming for this case. Adaptation of linear programming to L1 norm estimation problem is shown below: 

min: 1n
T(w+v) 

ß 

s.to: Xß+In(w-v)=y                                                                                                    (31) 

w,v≥0 

ß unrestricted in sign 

where 1n is a vector of size nx1 of 1's and In is a nth order identity matrix. The vectors v and w are of size nx1 and their 

elements may be interpreted as vertical deviations above and below the fitted regression hyperplane respectively. This problem 

has n equality constraints in m+2n variables. When n is large, this formulation generally requires a large amount of storage and 

computation time. 

Wagner (1959) shows that the formulation of the L1 norm regression may be reduced to m equality constraints linear 

programming problem. Thus, this dual formulation reduces n equations of primal form to m equations of dual form and 

considerably reduces the storage and computation time. 

Fisher (1961) reviews the formulation of the L1 norm estimation in relation to the primal form of linear programming. 

Barrodale and Young (1966) developed a modified simplex algorithm for determining the best fitting function to a set of 

discrete data under the L1 norm criterion. The method is given as Algol codes (for critics see, McCormick and Sposito (1975)). 

Davies (1967) demonstrates the use of the L1 norm regression estimates. Rabinowitz (1968) also discusses the application of 

linear programming in this field. Crocker (1969) cautions against using the L1 norm criterion merely to restrain unwanted 

negative coefficient estimates which occur in the least squares regression.  Multicolinearity is one of the cases which causes this 

result. Robers and Ben-Israel (1969) by using interval linear programming, proposed an algorithm to solve the L1 norm 

estimation problem. Rabinowitz (1970), Shanno and Weil (1970) discuss some connections between linear programming and 

the approximation problem. Barrodale (1970) summarizes the linear and nonlinear L1 norm curve fitting on both continuous 



Copyright © CC-BY-NC 2019, BJMSR 

 

 www.cribfb.com/journal/index.php/BJMSR       Bangladesh Journal of Multidisciplinary Scientific Research         Vol. 1, No. 1; 2019 

59 
               

and discrete data. Spyropoulos and Kiountouzis and Young (1973) suggest two algorithms for fitting general functions and 

particularly fast algorithm with minimum storage requirements for fitting polynomials based on the algebraic properties of linear 

programming formulation. Robers and Robers (1973) have supplied a special version of the general method of Robers and Ben-

Israel (1969), which is designed specifically for the L1 norm problem. A related Fortran code is also provided. 

Barrodale and Roberts (1973) present a modification of the simplex method, which needs a smaller amount of storage, 

and by skipping over simplex vertices is more efficient than the usual simplex procedure. Define the vector ß as a difference of 

two nonnegative vectors c and d; their formulation can be stated as follows, 

min: 1n
T(w+v) 

c,d 

s.to: X(c-d)+In(w-v)=y                                                                                              (32) 

w,v,c,d≥0 

Because of the relationships among variables, the computation can be performed by using only (n+2)x(m+2) amount of array 

storage, including labels for the basic and non-basic vectors.  An initial basis is given by w if all yi are nonnegative. If a yi is 

negative, the sign of the corresponding row is changed, and the unit column from the corresponding element of v is taken as part 

of the basis. The algorithm is implemented in two stages. The first stage restricts the choice of the pivotal column during the 

first m iterations to the elements of the vector cj and dj according to the associated maximum nonnegative marginal costs. The 

vector that leaves the basis causes the maximum decrease in the objective function. Thus the pivot element is not necessarily the 

same as in the usual simplex.  The second stage involves interchanging nonbasic wi or vi with the basic wi or vi. The basic 

vectors corresponding to cj and dj are not allowed to leave the basis. The algorithm terminates when all marginal costs are 

nonpositive (see, Kennedy and Gentle (1980)). Fortran code for this procedure is given by Barrodale and Roberts (1974). Peters 

and Willms (1983) give algorithms accompanying with computer codes for up-and-down dating the solution of the problem 

when a column or row inserted to or deleted from X, or y is changed. These algorithms are all based on Barrodale and Roberts 

(1973,74) procedure. 

Abdelmalek (1974) describes a dual simplex algorithm for the L1 norm problem with no use of artificial variables. For 

this algorithm, the Haar condition (see, Osborne (1985), Moroney (1961)) need not be satisfied anymore. This algorithm 

seemed to be very efficient at the time of publication. An improved dual simplex algorithm for L1 norm approximation is 

proposed by Abdelmalek (1975a). In this algorithm, certain intermediate iterations are skipped, and in the case of ill-

conditioned problems, the basis matrix can lend itself to triangular factorization and thus ensure a stable solution.  Abdelmalek 

(1980a) improves his previous algorithm by using triangular decomposition. A Fortran translation of the algorithm is given by 

Abdelmalek (1980b).  Sposito and McCormick and Kennedy (1975) summarize much of the works on L1 norm estimation 

including problem statement, linear programming formulation, efficient computational algorithms, and properties of the 

estimators. 

Armstrong and Kung (1978) propose an algorithm for a simple two-parameter L1 norm regression. The method is a 

specification of linear programming of Barrodale and Roberts (1973) algorithm. A Fortran code is given too. 

Armstrong and Frome and Kung (1979) use LU (Lower-Upper triangular) decomposition of Bartels and Golub (1969) 

in maintaining the current basis on the revised simplex procedure. A Fortran translation is also enclosed. Armstrong and Godfrey 

(1979) show that the primal method of Barrodale and Roberts (1973) and the dual method of Abdelmalek (1975a) are 

essentially equivalent. With a given initial basis for the two methods, they show that both algorithms will generate corresponding 

bases at each iteration. The only difference is the choice of initial basis and heuristic rules for breaking ties.  Armstrong  and 

Kung (1982b) present a dual linear programming formulation for the problem. Various basis entry and initialization procedures 

are considered. It has been shown that the dual approach is superior to primal one if a good dual feasible solution is readily 

available (see also, Steiger (1980)). Banks and Taylor (1980) suggest a modification of Barrodale and Roberts (1973) 

algorithm. The objective function is altered to include magnitudes of the elements of both errors and solution vectors. For a 

general discussion on simplex for piecewise linear programming see Fourer (1985a,b) and for a survey of the corresponding 



Copyright © CC-BY-NC 2019, BJMSR 

 

 www.cribfb.com/journal/index.php/BJMSR       Bangladesh Journal of Multidisciplinary Scientific Research         Vol. 1, No. 1; 2019 

60 
               

problem on the L1 norm see Fourer (1986).  Narula and Wellington (1987) propose an efficient linear programming algorithm 

to solve both L1 and L¥ norms linear multiple regressions. The algorithm exploits the special structure and similarities between 

the two problems. 

Brennan and Seiford (1987) develop a geometrical interpretation of linear programming in L1 norm regression.  They 

give a geometric insight into the solving process in the space of observations. McConnell (1987) shows how the method of 

vanishing Jacobians which has been used to optimize quadratic programming problems can also be used to solve the special 

linear programming problem associated with computing linear discrete L1 norm approximation. For the possibility of applying 

other types of linear programming solutions such as Karmarkar solution to L1 norm problem see Meketon (1986). 

 
Other algorithms 

This category consists of algorithms which were not classified in the two last sections. 

Rice (1964c) applies the bisection method to L1 norm regression. In this method at each step, the domain of S is 

broken to two segments, and the appropriate segment is selected for the next iteration. The solution is reached when the last 

segment is less than a predetermined small value. 

Abdelmalek (1971) develops an algorithm for fitting functions to discrete data points and solving the overdetermined 

system of linear equations. The procedure is based on determining L1 norm solution as the limiting case of Lp norm 

approximation when p tends to one from right in the limit. This technique thus obtains a solution to a linear problem by solving 

a sequence of nonlinear problems. 

Schlossmacher (1973) computed the L1 norm estimates of regression parameters by an iterative weighted least squares 

procedure. Instead of minimizing the sum of absolute deviations, he minimized he sum of weighted squared errors with 1/│ui│ 

as weights. Once the least squares is applied to the problem and residuals are computed. The absolute value of the inverse of  the 

residuals are again used as corresponding weights in the next iteration for minimizing the sum of weighted squared errors (see 

also, Holland and Welsh (1977)). Fair (1974) observed that the estimated values of ß did not change after the second or third 

iterations. In cases where any residual is zero, the continuation of the procedure is impossible, because the corresponding weight 

to this residual is infinite. This problem is also discussed by Sposito and Kennedy and Gentle (1977), Soliman and Christensen 

and Rouhi (1988). Absolute convergence of this algorithm has not been proved, but the non-convergent experiment has not 

been reported. 

Soliman and Christensen and Rouhi (1988) used left pseudoinverse (see, Dhrymes (1978) for a description of this 

inverse) to solve the general linear L1 norm regression.  According to this procedure, one should calculate the least squares 

solution using the left pseudoinverse or least squares approximation as ß^=(XTX)-1XTy. Calculate the residuals as u=│y-

Xß^│; where u is nx1 column vector. Select the m observations with the smallest absolute values of the residuals and partition 

the matrices as the selected observations locate on the top, 

, , .

    
    

      
    

    

u X y

u X y

u yX

                                                                                       (33) 

Solve y^=X^ß^ for the top partitions as ß^=X^-1y. Although this procedure is operationally simple, its solution is not the 

same as other exact methods, and no proof is presented to show that the solution is in the neighborhood of the exact solution of 

the L1 norm minimization problem. 

Application of median polish (see, Tukey (1977)) and e-median polish to L1 norm estimation are discussed and 

developed by Bloomfield and Steiger (1983), Kemperman (1984), Sposito (1987a), Bradu (1987a,b). 



Copyright © CC-BY-NC 2019, BJMSR 

 

 www.cribfb.com/journal/index.php/BJMSR       Bangladesh Journal of Multidisciplinary Scientific Research         Vol. 1, No. 1; 2019 

61 
               

Application of Karmarkar's algorithm for linear programming and its relation to L1 norm is given by Sherali and 

Skarpness and Kim (1987). For using homotopy method in L1 norm, see Garcia and Gould (1983), Schellhorn (1987). An 

algorithm for linear L1 norm approximation for the continuous function is given by Watson (1981), (see also, Baboolal and 

Watson (1981)).  

 

References 
N.N. Abdelmalek (1971) Linear approximation for a discrete point set and L1 solutions of overdetermined linear equations. J. 

ACM, 18, 41-47. 
 
N.N. Abdelmalek (1974) On the discrete linear L1 approximation and L1 solutions of overdetermined linear equations. J. of 

Approx. Theory, 11, 38-53. 
 
N.N. Abdelmalek (1975a) An efficient method for the discrete L1 approximation problem. Math. Comput., 29, 844-850. 
 
N.N. Abdelmalek (1980a) L1 solution of overdetermined systems of linear equations. ACM Trans. Math. Soft., 6, 220-227. 
 
N.N. Abdelmalek (1980b) A Fortran subroutine for the L1 solution of overdetermined systems of linear equations. ACM 

Trans. Math. Soft., 6, 228-30. 
 
D.H. Anderson, M.R. Osborne (1976) Discrete linear approximation problems in polyhedral norms. Numer. Math. 26, 179-

189. 
 
R.D. Armstrong, E.L. Frome, D.S. Kung (1979) A revised simplex algorithm for the absolute deviation curve fitting problem. 

Commun. Stat. B8, 175-190. 
 
R.D. Armstrong, J. Godfrey (1979) Two linear programming algorithms for the discrete L1 problem. Math. Comput., 33, 289-

300. 
 
R.D. Armstrong, D.S. Kung (1978) AS132: Least absolute value estimates for a simple linear regression problem. Appl. Stat., 

27, 363-366. 
 
R.D. Armstrong, D.S. Kung (1982b) A dual algorithm to solve linear least absolute value problems. J. Oper. Res. Soc., 33, 

931-936. 
 
T.S. Arthanari, Y. Dodge (1981) Mathematical programming in statistics. John Wiley, Interscience division, New York. 
 
S. Baboolal, G.A. Watson (1981) Computational experience with an algorithm for discrete L1 approximation. Computing, 27, 

245-252. 
 
S.C. Banks, H.L. Taylor (1980) A modification to the discrete L1 linear approximation algorithm of Barrodale and Roberts. 

SIAM J. on Scientific and Stat. Comput. 1, 187-190. 
 
Barrodale (1970) On computing best L1 approximations. In A. Talbot, Approximation theory Academic Press, New York, 

205-215. 
 
Barrodale, F.D.K. Roberts (1973) An improved algorithm for discrete L1 linear approximation. SIAM J. Numer. Anal., 10, 

839-848. 
 
Barrodale, F.D.K. Roberts (1974) Algorithm 478: Solution of an overdetermined system of equations in the L1 norm. 

Commun. ACM, 17, 319-320. 
 
Barrodale, A. Young (1966) Algorithms for best L1 and L∞ linear approximations on a discrete set. Numer. Math., 8, 295-306. 
 



Copyright © CC-BY-NC 2019, BJMSR 

 

 www.cribfb.com/journal/index.php/BJMSR       Bangladesh Journal of Multidisciplinary Scientific Research         Vol. 1, No. 1; 2019 

62 
               

R.H. Bartels, A.R. Conn (1977) LAV regression: A special case of piecewise linear minimization. Commun. Stat., B6, 329-
340. 

 
R.H. Bartels, A.R. Conn, J. Sinclair (1976) The L1 solution to an overdetermined linear system. Proc. 9th Ann. Symp. Interface 

Statist. In C.D.C. Hoaglin (ed.) Boston, Prindle, Weber and Schmidt Inc., 120-7. 
 
R.H. Bartels, A.R. Conn, J. Sinclair (1978) Minimization technique for piecewise differentiable functions: The L1 solution to 

an overdetermined linear system. SIAM J. Numer. Anal. 15, 224-241. 
 
R.H. Bartels, G.H. Golub (1969) The simplex method of linear programming using LU decomposition. Commun. ACM, 12, 

266-268. 
 
J. Bejar (1956) Regression en mediana y la programación lineal, Trabajos de Estadistica 7, 141-58. 
 
J. Bejar (1957) Calculo practico de la regression en mediana Trabajos de Estadistica, 8, 157-173. 
 
Bijan Bidabad (1987a) Least absolute error estimation. The First International Conference on Statistical Data Analysis Based 

on the L1‎‎ norm and Related Methods, Neuchatel, Switzerland. http://www.bidabad.com/doc/lae-I.pdf 
 
Bijan Bidabad (1987b) Least absolute error estimation, part II. Submitted to the First International Conference on Statistical 

Data Analysis Based on the L1‎‎ norm and Related Methods, Neuchatel, Switzerland. 
http://www.bidabad.com/doc/lae-II.pdf 

 
Bijan Bidabad (1988a) A proposed algorithm for least absolute error estimation. Proc. of the Third Seminar of Mathematical 

Analysis. Shiraz Univ., 24-34, Shiraz, Iran.  
 
Bijan Bidabad (1988b) A proposed algorithm for least absolute error estimation, part II. Proc. of the Third Seminar of 

Mathematical Analysis, Shiraz Univ., 35-50, Shiraz, Iran.  
 

Bijan Bidabad (1989a) Discrete and continuous L1‎‎ norm regressions, proposition of discrete approximation algorithms and 
continuous smoothing of concentration surface, Ph.D. thesis, Islamic Azad Univ., Tehran, Iran. 
http://www.bidabad.com/doc/L1-norm-thesis-en.pdf  

 

Bijan Bidabad (1989b) Discrete and continuous L1‎‎ norm regressions, proposition of discrete approximation algorithms and 
continuous smoothing of concentration surface, Ph.D. thesis, Islamic Azad Univ., Tehran, Iran. Farsi translation. 
http://www.bidabad.com/doc/L1-norm-thesis-fa.pdf  

 
Bijan Bidabad (2005). L1 norm based computational algorithms. http://www.bidabad.com/doc/l1-article6.pdf 

Bijan Bidabad (2005). L1 norm solution of overdetermined system of linear equations. http://www.bidabad.com/doc/l1-

article5.pdf 

Bijan Bidabad (2005). L1 norm based data analysis and related methods. http://www.bidabad.com/doc/l1-articl1.pdf  

Bijan Bidabad (2005). New algorithms for the L1 norm regression. http://www.bidabad.com/doc/l1-article2.pdf    

Bijan Bidabad (2005). Comparative study of the L1 norm regression algorithms. http://www.bidabad.com/doc/l1-articl3.pdf   

Bijan Bidabad (2005). Continuous L1 norm estimation of Lorenz curve. http://www.bidabad.com/doc/l1-articl4.pdf  

Bijan Bidabad (1993). Estimating Lorenz curve for Iran by using continuous L1 norm estimation, Economics and Management 

Journal, Islamic Azad University, No. 19, winter 1993, pp. 83-101. http://www.bidabad.com/doc/iraninc-l1.pdf  

Bijan Bidabad (2005). Continuous L1 norm estimation of Lorenz curve when probability density function is known. 

Bijan Bidabad (2005). USA Income distribution counter-business-cyclical trend (Estimating Lorenz curve using continuous L1 

norm estimation). First meeting of the Society for the Study of Economic Inequality (ECINEQ), Palma de Mallorca, 

Spain, July 20-22, 2005. 

http://www.bidabad.com/doc/lae-I.pdf
http://www.bidabad.com/doc/lae-II.pdf
http://www.bidabad.com/doc/L1-norm-thesis-en.pdf
http://www.bidabad.com/doc/L1-norm-thesis-fa.pdf
http://www.bidabad.com/doc/l1-article6.pdf
http://www.bidabad.com/doc/l1-article5.pdf
http://www.bidabad.com/doc/l1-article5.pdf
http://www.bidabad.com/doc/l1-article1.pdf
http://www.bidabad.com/doc/l1-article2.pdf
http://www.bidabad.com/doc/l1-article3.pdf
http://www.bidabad.com/doc/l1-article3.pdf
http://www.bidabad.com/doc/l1-article4.pdf
http://www.bidabad.com/doc/iraninc-l1.pdf


Copyright © CC-BY-NC 2019, BJMSR 

 

 www.cribfb.com/journal/index.php/BJMSR       Bangladesh Journal of Multidisciplinary Scientific Research         Vol. 1, No. 1; 2019 

63 
               

http://www.uib.es/congres/ecopub/ecineq/general.html 

http://www.uib.es/congres/ecopub/ecineq/papers/039Bidabab.pdf  

http://www.bidabad.com/doc/estimating-lorenz-us.pdf   

Bijan Bidabad, Hamid Shahrestani. (2008) An implied inequality index using L1 norm estimation of Lorenz curve. Global 

Conference on Business and Finance Proceedings. Mercedes Jalbert, managing editor, ISSN 1931-0285 CD, ISSN 

1941-9589 Online, Volume 3, Number 2, 2008, The Institute for Business and Finance Research, Ramada Plaza 

Herradura, San Jose, Costa Rica, May 28-31, 2008, pp. 148-163. Global Journal of Business Research, Vol. 4, No. 

1, 2010, pp.29-45. 

http://www.bidabad.com/doc/L1-Implied-inequality-index-4.pdf 

http://www.theibfr.com/archive/ISSN-1941-9589-V3-N2-2008.pdf 

http://www.bidabad.com/doc/SSRN-id1631861.pdf 
 
P. Bloomfield, W. Steiger (1980) Least absolute deviations curve fitting. SIAM J. Sci. Statist. Comput. 1, 290-301. 
 
P. Bloomfield, W. Steiger (1983) Least absolute deviations: theory, applications and algorithms. Birkhauser, Boston. 
 
R.J. Boscovich (1757) De litteraria expeditione per pontificiam ditionem, et synopsis amplioris operis..., 'Bononiensi' 

scientiarum et artum instituto atque academia commetarii, vol.4, 353-396. Reprinted with a Serbo-Croatian 
translation by N. Cubranic, Institute of higher geodesy, University of Zagreb 1961. 

 
R.J. Boscovich (1760) De recentissimis graduum dimensionibus et figura, ac magnitudine terrae inde derivanda. Philosophiae 

recentioris, a benedicto stay in Romano archigynasis publico eloquentare professore, vesibus traditae, libri X, cum 
adnotianibus et supplementas P. Rugerii Joseph Boscovich, S.J., 2, 406-426. 

 
A.L. Bowley (1902) Methods of representing the statistics of wages and other groups not fulfilling the normal law of error, II: 

applications to wage statistics and other groups. J. of the Roy. Stat. Soc., 65, 331-54. 
 
A.L. Bowley (1928) F.Y. Edgeworth's contributions to mathematical statistics. London, Roy. Stat. Soc.. 
 
D. Bradu (1987a) L1 fit, median polish and conjugate gradients. CSIR Tech. Rep. TWISK 509, National Res. Inst. For Math 

Sci. CSIR, Pretoria. 
 

D. Bradu (1987b) An ε-median polish algorithm. CSDA, 5, 327-336. 
 
J.J. Brennan, L.M. Seiford (1987) Linear programming and L1 approximation using the method of vanishing Jacobians. CSDA, 

5, 263-276. 
 
B.M. Brown (1980) Median estimates in a simple linear regression. Australian J. of Stat., 22, 154-165. 
 
C. Bruen (1938) Methods for the combination of observations modal points or most lesser-deviations, mean loci or least 

squares, and mid point of least range or least greatest-deviation. Metron 13, 61-140. 
 
J. Chamber (1971) Algorithm 410: partial sorting. Comm. ACM, 14, 357-358. 
 
J.M. Chambers (1977) Computational methods for data analysis Wiley, New York. 
 
Charnes, W.W. Cooper, R.O. Ferguson (1955) Optimal estimation of executive compensation by linear programming. Manag. 

Sci. 1, 138-151. 
 
E.W. Cheney (1966) Introduction to approximation theory, McGraw-Hill, New York. 
 

http://www.uib.es/congres/ecopub/ecineq/general.htm
http://www.uib.es/congres/ecopub/ecineq/papers/039Bidabab.pdf
http://www.bidabad.com/doc/estimating-lorenz-us.pdf
http://www.bidabad.com/doc/L1-Implied-inequality-index-4.pdf
http://www.theibfr.com/archive/ISSN-1941-9589-V3-N2-2008.pdf
http://www.bidabad.com/doc/SSRN-id1631861.pdf


Copyright © CC-BY-NC 2019, BJMSR 

 

 www.cribfb.com/journal/index.php/BJMSR       Bangladesh Journal of Multidisciplinary Scientific Research         Vol. 1, No. 1; 2019 

64 
               

F.H. Clarke (1983) Optimization and nonsmooth analysis. Wiley, New York 
 
A.R. Conn (1976) linear programming via a nondifferentiable penalty function. SIAM J. Numer. Anal., 13, 145-154. 
 
D.C. Crocker (1969) Linear programming technique in regression analysis, the hidden danger. A.I.E.E. Trans., 1, 112-126. 
 
M. Davies (1976) Linear approximation using the criterion of least total deviations. J. Roy. Stat. Soc. B29, 101-109. 
 
T.E. Dielman (1984) Least absolute value estimation in regression models: An annotated bibliography. Comm. Stat. 13, 513-

41. 
 
T. Dielman, R. Pfaffenberger (1982) LAV (Least Absolute Value) estimation in linear regression: A review, TIMS studies in 

the Manag. Sci.,19, 31-52. 
 
T. Dielman, R. Pfaffenberger (1984) Computational algorithms for calculating least absolute value and Chebyshev estimates 

for multiple regression. Amer. J. Math. Manag. Sci., 4, 169-197. 
 
P.J. Dhrymes (1978) Mathematics for econometrics. Springer-Verlag, New York. 
 
Y. Dodge (1987) An introduction to statistical data analysis L1-norm based. In Y. Dodge (ed.) Statistical data analysis based 

on the L1 norm and related methods. North-Holland. Reprinted in CSDA, 5, 239-254. 
 
A.F. Dufton (1928) Correlation. Nature, 121, 866. 
 
F.Y. Edgeworth (1883) The method of least squares. Philosophical Magazine, 16, 360-375. 
 
F.Y. Edgeworth (1887a) On observations relating to several quantities. Hermathena, 6, 279-285. 
 
F.Y. Edgeworth (1887b) A new method of reducing observations relating to several quantities. Philosophical Magazine, 24, 

222-223. 
 
F.Y. Edgeworth (1888) On a new method of reducing observation relating to several quantities. Philosophical Magazine, 25, 

184-191. 
 
F.Y. Edgeworth (1902) Method of representing statistics of wage and other groups not fulfilling the normal law of error, I: 

mathematical considerations. J. Roy. Stat. Soc., 65, 325-331. 
 
F.Y. Edgeworth (1923) On the use of medians for reducing observations relating to several quantities. Philosophical Magazine, 

6th series, 46, 1074-1088. 
 
Eisenhart (1961) Boscovich and the combination of observations. Ch. 9 of Whyte (1961, 200-212) reprinted in Kendall and 

Plackett (1977) studies in the history of statistics and probability, vol.II, Charles Griffin and Co. Ltd., High 
Wycombe 88-100. 

 
J.E. Estienne (1926-28) Introduction a une theorie rationnelle des erreurs d'observation. Revue d'artillerie 97(1926), 421-441; 

98(1928), 542-562; 100(1927), 471-487. 
 
R.C. Fair (1974) On the robust estimation of econometric models. Ann. Econ. Soc. Measurement, 3, 667-77. 
 
R.W. Farebrother (1987b) The historical development of the L1 and L∞ estimation procedures. In Y. Dodge (ed.) Statistical 

data analysis based on the L1 norm and related methods. North-Holland. 37-64. 
 
R.W. Farebrother (1987c) A simple recursive procedure for the L1 norm fitting of a straight line. Work. Pap., Dept. of 

Econometrics and Social Stat. University of Manchester, Manchester, M13 9PL, UK. 
 
W.D. Fisher (1961) A note on curve fitting with minimum deviations by linear programming. JASA, 11, 359-362. 



Copyright © CC-BY-NC 2019, BJMSR 

 

 www.cribfb.com/journal/index.php/BJMSR       Bangladesh Journal of Multidisciplinary Scientific Research         Vol. 1, No. 1; 2019 

65 
               

 
R. Fourer (1985a) A simplex algorithm for piecewise-linear programming I: derivation and proof. Math. Prog., 33, 204-233. 
 
R. Fourer (1985b) A simplex algorithm for piecewise-linear programming II: Finiteness, feasibility and degeneracy. Tech. Rep., 

85-03 (revised), Dept. of Ind. Engin. and Manag. Sci.,The Tech. Inst., Northwestern University, Evanston, Illinois. 
 
R. Fourer (1986) A simplex algorithm for piecewise-linear programming III: Computational analysis and applications. Tech. 

Rep., 86-03, Dept. of Ind. Engin. and Manag. Sci., The Tech. Inst., Northwestern University, Evanston, Illinois. 
 
J.B.I. Fourier (1824) Solution d'une question particuliere au calcul des inegalites, second extrait. Histoire de l'academie des 

sciences pour 1824, 47-55. Reprinted in oeuvres de Fourier, 2. Paris, 1980, Gauthier-Villars, 325-328. 
 
G. Galilei (1632) Dialogo dei massimi sistemi. 
 
C.B. Garcia, F.G. Gould (1983) An application of homotopy to solving linear programs. Math. Prog. 27, 263-282. 
 
C.F. Gauss (1809) Theoria motus corporum coelestium. In F. Perthes, I.H. Besser, Sectionbus conicis solem ambientium, 

Hamburg. Reprinted in his werke, vol. 7, F. Pethes, Gotha 1871. English translation by C.H. Davis, Little, Brown 
and Co., Boston, 1857. Reprinted by Dover Pub. New York, 1963. 

 
T.E. Harris (1950) Regression using minimum absolute deviations. Am. Statist., 4, 14-15. 
 
H.L. Harter (1974a) The method of least squares and some alternative, I. Int. Stat. Rev., 42, 147-174. 
 
H.L. Harter (1974b) The method of least squares and some alternative, II. Int. Stat. Rev., 42, 235-264. 
 
H.L. Harter (1975a) The method of least squares and some alternative, III. Int. Stat. Rev., 43, 1-44. 
 
H.L. Harter (1975b) The method of least squares and some alternative, IV Int. Stat. Rev., 43, 125-190, 273-278. 
 
H.L. Harter (1975c) The method of least squares and some alternative, V. Int. Stat. Rev., 43, 269-272. 
 
H.L. Harter (1976) The method of least squares and some alternative, VI. Int. Stat. Rev., 44, 113-159. 
 
P.W. Holland, R.E. Welsch (1977) Robust regression using iteratively reweighted least-squares. Comm. Stat., A6, 813-827. 
 
L. Horvath (1987) Asymptotic normality of Lp-norms of density estimators. Tech. Rep. series of Lab Res. Stat. Prob., no.3, 

Carleton University, Ottawa, Canada. 
 
H. Imai, K. Kato, P. Yamamoto (1987) A linear-time algorithm for linear L1 approximation of points. Tech. Rep. CSCE-87-

C30. Dept. of Comp. Sci. and Commun. Engin., Kyushu University 36, Fukuoka 812, Japan. 
 
L.A. Josvanger, V.A. Sposito (1983) L1-norm estimates for the simple regression problem. Comm. Stat. B12, 215-21. 
 
O.J. Karst (1958) Linear curve fitting using least deviations. JASA, 53, 118-132. 
 
Y. Kawara (1979) Straight line fitting by minimizing the sum of absolute deviations. J. of the Japan Stat. Soc., 9, 47-64. 
 
J.H.B. Kemperman (1984) Least absolute value and median polish. In Y.L. Tong (ed.), Inequalities in statistics and probability 

(IMS Lecture notes monograph series, vol.5), Inst. of Math. Stat., Hayward, CA, 84-113. 
 
W.J. Kennedy, J.E. Gentle (1980) Statistical computing. New York, Marcel Dekker. 
 
P.S. Laplace (1793) Sur quelques points du system du monde. Memoires de l'Academie Royale des Science de Paris. Annee 

1789, 1-87. Reprinted in Oeuvres completes de Laplace II. Paris, Gauthier-Villars, 1985, 477-558. 
 



Copyright © CC-BY-NC 2019, BJMSR 

 

 www.cribfb.com/journal/index.php/BJMSR       Bangladesh Journal of Multidisciplinary Scientific Research         Vol. 1, No. 1; 2019 

66 
               

P.S. Laplace (1799) Traite des mecanique celeste, 2. Paris; J.B.M. Depart. Reprinted as oeuvres completes de Laplace, 2. Paris; 
Gauthier-Villars 1878, 116-165. 

 
P.S. Laplace (1812) Theorie analytique des probabilites, Mme courcier Paris 1820 Reprinted in his oeuvres, vol.7, Imprimerie 

Royale, Paris, 1847, and Gauthier-Villars et fils, Paris 1886. 
 
P.S. Laplace (1818) Duexieme supplement to Laplace (1812). 
 
C.L. Mathieu (1816) Sur les experiences du pendule, faites par les navigateurs espagnol, en differens points du globe. 

Connaissance des tems, 314-332. 
 
C.R. McConnell (1987) On computing a best discrete L1 approximation using the method of vanishing Jacobians. CSDA, 5, 

277-288. 
 
G.F. McCormick, V.A. Sposito (1975) A note on L1 estimation based on the median positive quotient. Appl. Stat., 24, 347-

350. 
 
M.S. Meketon (1986) Least absolute value regression. Work. Pap., AT&T Bell Laboratories, Holmdel, N.J. 
 
R.M. Moroney (1961) The Haar problem in L1. Proc. Amer. Math. Soc., 12, 793-795. 
 
S.C. Narula, J.F. Wellington (1985) Interior analysis for the minimum sum of absolute errors regression. Technometrics, 27, 

181-188. 
 
S.C. Narula, J.F. Wellington (1987) An efficient algorithm for the MSAE and MMAE regression problems. Work. Pap., 

Virginia CommonWealth University, Richmond, VA 23284. 
 
M.R. Osborne (1987) The reduced gradient algorithm. In Y. Dodge (ed.) Statistical data analysis based on the L1 norm and 

related methods. North-Holland. 95-108. 
 
M.R. Osborne, S.A. Pruess, R.S. Womersley (1986) Concise representation of generalized gradients. J. of Austra. Math. Soc., 

Ser. B, 28, 57-74. 
 
M.R. Osborne, G.A. Watson (1985) An analysis of the total approximation problem in separable norms, and an algorithm for 

the total L1 problem. SIAM J. Sci. Stat. Comp., 6, 410-424. 
 
M.J. Panik (1976) Classical optimization: foundation and extensions. North-Holland, Amsterdam. 
 
U. Peters, C. Willms (1983) Up- and down-dating procedures for linear L1 regression. OR Spektrum 5, 229-239. 
 
P. Pilibossian (1987) A direct solving algorithm for a linear regression according to L1-norm criteria. Work. Pap., L.S.T.A. 

Universite, Paris VI 
 
P. Rabinowitz (1968) Application of linear programming to numerical analysis. SIAM Rev., 10, 121-159. 
 
P. Rabinowitz (1970) Mathematical programming and approximation. In A. Talbot (ed.) Approximation Theory. Academic 

Press, 217-231. 
 
M.R. Rao, V. Srinivasan (1972) A note on Sharpe's algorithm for minimum sum of absolute deviations in a simple regression 

problem. Manag. Sci., 19, 222-225. 
 
E.C. Rhodes (1930) Reducing observations by the method of minimum deviations. Philo. Mag., 7th series, 9, 974-92. 
 
J.R. Rice (1964c) The approximation of functions, vol. I, linear theory. Reading Mass:, Addison-Wesley. 
 



Copyright © CC-BY-NC 2019, BJMSR 

 

 www.cribfb.com/journal/index.php/BJMSR       Bangladesh Journal of Multidisciplinary Scientific Research         Vol. 1, No. 1; 2019 

67 
               

P.D. Robers, A. Ben-Israel (1969) An interval programming algorithm for discrete linear L1 approximation problem. J. Approx. 
Theory, 2, 323-336. 

 
P.D. Robers, S.S. Robers (1973) Algorithm 458: discrete linear L1 approximation by interval linear programming. Comm. 

ACM, 16, 629-633. 
 
E. Ronchetti (1987) Bounded influence in regression: a review. In Y. Dodge (ed.) Statistical data analysis based on the L 1 norm 

and related methods. North-Holland, 65-80. 
 
A.N. Sadovski (1974) AS74: L1-norm fit of a straight line. Appl. Stat. 23, 244-248. 
 
J.P. Schellhorn (1987) Fitting data through homotopy methods In Y. Dodge (ed.) Statistical data analysis based on the L 1 

norm and related methods. North-Holland. 131-138. 
 
E.J. Schlossmacher (1973) An iterative technique for absolute deviations curve fitting. JASA 68, 857-865. 
 
E. Seneta (1983) The weighted median and multiple regression. Austral. J. Stat., 25(2), 370-377. 
 
E. Seneta, W.L. Steiger (1984) A new LAD curve-fitting algorithm: slightly overdetermined equation system in L1. Discrete 

Applied Math., 7, 79-91. 
 
Shanno, R.L. Weil (1970) Linear programming with absolute value functionals. Oper. Res., 19, 120-124. 
 
W.F. Sharpe (1971) Mean-absolute deviation characteristic lines for securities and portfolios. Manag. Sci., 18, B1-B13. 
 
H.D. Sherali, B.O. Skarpness, B. Kim (1987) An assumption-free convergence analysis for a perturbation of the scaling 

algorithm for linear programs, with application to the L1 estimation problem. Dept. of Ind. Engin. and OR, Virginia 
Polytechnic Inst. and State University, Blacksburg, Virginia. 

 
O.B. Sheynin (1973) R.J. Boscovich's work on probability. Archive for history of exact sciences, vol. 9, 306-324, and vol. 28, 

173. 
 
R.R. Singleton (1940) A method for minimizing the sum of absolute values of deviations. Annals of math. Stat., 11, 301-310. 
 
S.A. Soliman, G.S. Christensen, A. Rouhi (1988) A new technique for curve fitting based on minimum absolute deviations. 

CSDA, 6(4), 341-352. 
 
V.A. Sposito (1976) A remark on algorithm AS74, L1 norm fit of a straight line. Appl. Stat., 25, 96-97. 
 
V.A. Sposito (1987a) On median polish and L1 estimators.CSDA, 5, 155-162. 
 
V.A. Sposito, W.J. Kennedy, J.E. Gentle (1977) AS110: Lp norm fit of a straight line. Appl. Stat., 26, 114-118. 
 
V.A. Sposito, G.F. McCormick, W.J. Kennedy (1975) L1 estimation strategies based on the simplex algorithm. In Proc. of the 

eighth symposium on the interface, J.W. France (ed.) Health science computing facility. UCLA, Los Angeles. 
 
V.A. Sposito, W.C. Smith (1976) On a sufficient and necessary condition for L1 estimation. Appl. Stat., 25, 154-157. 
 
K. Spyropoulos, E. Kiountouzis, A. Young (1973) Discrete approximation in the L1 norm. Comp. J., 16, 180-186. 
 
W.L. Steiger (1980) Linear programming via L1 curve fitting beats simplex. Abstracts, AMS, 80T-C26, 385-386. 
 
S.M. Stigler (1981) Gauss and invention of least squares. Annals of Stat., 9, 465-474. 
 
S.M. Stigler (1984) Studies in the history of probability and statistics XL, Boscovich, Simpson and a 1760 manuscript note on 

fitting a linear relation. Biometrica, 71, 3, 615-620. 



Copyright © CC-BY-NC 2019, BJMSR 

 

 www.cribfb.com/journal/index.php/BJMSR       Bangladesh Journal of Multidisciplinary Scientific Research         Vol. 1, No. 1; 2019 

68 
               

 
J. Svanberg (1805) Exposition des operations faites en lappnie pour la determination d'un arc du meridien en 1801, 1802 et 

1803,... Stockholm. 
 
J.W. Tukey (1977) Exploratory data analysis. Reading, Mass. Addison-Wesley. 
 
H.H. Turner (1887) On Mr. Edgeworth's method of reducing observations relating to several quantities. Phil. Mag. (5 th series), 

24, 466-470. 
 
K.H. Usow (1967a) On L1 approximation: computation for continuous functions and continuous dependence. SIAM J. of 

Numer. Anal., 4, 70-88. 
 
K.H. Usow (1967b) On L1 approximation: computation for discrete functions and discretization effect. SIAM J. Numer. 

Anal., 4, 233-244. 
 
J.F. Van Beeck-Calkoen (1816) Ver de theoric der Gemiddelde Waardij. Verhandlingen der K. Nederlandandsch Instituut Can 

Wetenschappen, 2, 1-19. 
 
B.A. Von Lindenau (1806) Uber den Gebrauch der Gradmessungen zur bestimmung der gestalt der erde. Monatliche 

correspondenz zur befar derung der Erd-und Himmels-kunde, 14, 113-158. 
 
H.M. Wagner (1959) Linear programming technique for regression analysis. JASA, 54, 202-212. 
 
G.A. Watson (1981) An algorithm for linear L1 approximation of continuous functions. IMA J. Num. Anal., 1, 157-167. 
 
G.O Wesolowsky (1981) A new descent algorithm for least absolute value regression problem. Comm. Stat., B10, 479-491. 
 
 
 
 
Copyrights  
Copyright for this article is retained by the author(s), with first publication rights granted to the journal. This is an open-access 
article distributed under the terms and conditions of the Creative Commons Attribution license 
(http://creativecommons.org/licenses/by/4.0/). 


