




































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  

 

     1 

New Algorithms for L1 Norm Regression 
 
 

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 
In this paper, we propose four algorithms for L1 norm computation of regression parameters, where two of them are more 
efficient for simple and multiple regression models.  However, we start with restricted simple linear regression and corresponding 
derivation and computation of the weighted median problem. In this respect, a computing function is coded.  With discussion 
on the m parameters model, we continue to expand the algorithm to include unrestricted simple linear regression, and two crude 
and efficient algorithms are proposed. The procedures are then generalized to the m parameters model by presenting two new 
algorithms, where the algorithm 4 is selected as more efficient. Various properties of these algorithms are discussed. 

 

Keywords: L1 norm, Regression, Algorithm, Computer program   
 
1. Introduction 
In Bidabad (1989a,b), various aspects of the L1 norm regression were reviewed. We observed that the L1 norm criterion is going 
to find its place in scientific analysis. Since it is not computationally comparable with other criteria such as L2 norm, it needs 
more work to make it a hand tool. The closed form of the solution of the L1 norm estimator has not been derived yet, and 
therefore, makes further inferences of the properties of this estimator difficult. Any attempt to give efficient computational 
algorithms which may introduce significant insight into the different characteristics of the problem is desirable. In this regard, we 
shall try to give a general procedure in this paper to solve the L1 norm linear regression problem. The proposed algorithms are 
based on a special descent method and use a discrete differentiation technique. Primary designs of the algorithms have been 
discussed by Bidabad (1987a,b,88a,b). By manipulating these algorithms, more efficient ones were introduced by Bidabad 
(1989a,b) which has been shown to have better performance than other existing algorithms. 

Consider the following regression model, 
            m 

     yi =  jxji + ui        i=1,...,n                                     (1) 
           j=1 

where j, j=1,...,m are population parameters to be estimated, yi, xji and ui are dependent, independent, and random variables 

respectively. We wish to estimate j's by minimizing the expression,  
             n                 n          m 

     S =   │ui^│ =   │yi -  ßj^xji│                                                                                       (2) 

            i=1             i=1        j=1 

Suppose m=1, then expression (2) reduces to, 
             n                        n                         n 

     S =   │yi - yi^│ =   │yi - ß^xi│ =   Si                                   (3) 

           i=1                    i=1                       i=1 

A typical element, Si=|yi-^xi| can be viewed as a broken line in the Sx^ plane composed of two half-lines. Si attains its 
minimum which is zero at,  

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 

2 
               

     i^ = yi/xi                                                                           (4) 

The slopes of the half-lines to the left and right of i^ are -|xi| and |xi| respectively. So, Si's are all convex, and hence their sum 

S is also convex with a slope at any ^ equal to the sum of the slopes of Si's at that value of ^. Since the slope of each Si 

changes only at the corresponding i^, the minimum of S will lie on one of the i^. Thus, the regression line will pass through 

origin with the slope equal to i^ which minimizes S. On the other hand, to find the L1 norm estimate of  we need to find only 
one observation (see, Karst (1958), Taylor (1974)). Furthermore, Taylor (1974) concludes that: "This implies, of course, that 
the regression line must pass through the observation corresponding to the minimizing i. The regression line, therefore, is 

determined by the point of origin and the observation associated with minimizing i^". But he did not continue this approach, 
that is, minimizing S with respect to subscript i. In this paper as Bidabad (1987a,88a) a first endeavor is to develop this point of 
view. In the next section, after rewriting (3) in a suitable fashion, the value of i is determined by using discrete differentiation. 
     By a similar discussion, it can be concluded that when the number of parameters is m, m observations must lie on the 
regression hyperplane. On the other hand, m equations of the form given by (5) are necessary to specify the regression 
hyperplane (see also Bidabad (1989a,b)). 
            m 

     yi -   j^xji = 0                                                                                                                   (5) 
           j=1 
Now, let us proceed with the simplest regression model. 
 
2. Restricted simple linear regression 

     For model (1) consider the case of one independent variable with no intercept, namely, y i=xi+ui. To find the L1 norm 

estimate of , the following procedure can be suggested. 
            n          n                  n                                        n   

     S =  |ui| =  |yi - xi| =  (yi - xi)sgn(yi - xi) =   |xi|(yi/xi - )sgn(yi/xi - )                 (6) 
          i=1       i=1               i=1                                    i=1 
Let zi=yi/xi and sort zi in descending order. Rename the resultant ordered zi, i=1,...,n to zh, h=1,...,n. Elements of zh should have 
the following property:  
     zh > zl  if h < l  for all h,l=1,...,n. 
Rewrite (6) with ordered observations as, 
             n 

     S =   |xh|(zh - )sgn(zh-)                                                                                                  (7) 
           h=1 
Let us denote the observation which will be on the regression line by (xt+1,yt+1) - that is the (t+1)th observation in the zh array. 
Value of zh is the slope of a ray passing through the origin and the hth observation. Therefore, 
 

     if, h < t+1  then, zh >  and  uh > 0 

     if, h = t+1  then, zh =  and  uh = 0 

     if, h > t+1  then, zh <  and  uh < 0 
We can now rewrite (7) as follows, 
            t                       n 

     S =  |xh|(zh - ) -     |xh|(zh - )                                                                                        (8) 
          h=1                 h=t+1 

To find the minimum of S, we need to differentiate it with respect to  and subscript t, because both  and t are unknowns. 

Note that  has continuous domain and t has discrete domain. Thus, 

     δS         t               n 

     ── = -  │xh│ +  │xh│                                                                                                   (9) 

     δß       h=1          h=t+1 

The differentiation of S with respect to subscript t must be discrete derivative (see, Bender and Orszag (1978), Clarke (1983)): 

     (S)                      S[t+(t)]-S(t) 

     ───  = lim      ──────────  = S(t+1)-S(t) = 



Copyright © CC-BY-NC 2019, BJMSR 

 

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

3 
               

     (t)        (t)─>1         (t) 

 

     t+1                     n                          t                        n 

       │xh│(zh-ß) -    │xh│(zh-ß) -   │xh│(zh-ß) +    │xh│(zh-ß) =  2│xt+1│(zt+1-ß) = 0 

     h=1                 h=t+2                    h=1                 h=t+1 

thus, 

     ß^ = zt+1 = yt+1/xt+1                                                                                                            (10) 

It should be noted that if xt+1=0 then the value of zt+1=l, and when the array z is sorted, infinite values of z locate at 
extremes of two tails of the array and make no problem in (10). Remember that equation (10) is the same as relation (4). 

Equations (9) and (10) are two equations with two unknowns t and . The variable t can be found by rewriting (9) as follows, 
               k              n 

     Dk ≡  │xh│ -    │xh│            k=1,...,n                                                                        (11) 

             h=1        h=k+1 

It is obvious that Dk is an increasing function of k. When k increases from one to n, Dk attains different values, which 
increase from negative to positive. So, initially, k is set equal to one, and Dk is computed accordingly. If Dk is negative, k is 
increased by one, and the procedure is repeated until Dk becomes positive. When Dk reaches the first positive value, then, 
t+1=k. By this procedure, the value of t+1 is found. The observation corresponding to this subscript (t+1=k) is selected 

(xt+1,yt+1). The L1 norm estimate of  is found by substituting the values of xt+1 and yt+1 into (10), (see, Bidabad (1987a,88a)). 
The presented procedure is essentially the weighted median procedure of Laplace (1818). Namely, this solution is the 

steps that Laplace took. The main difference is the mathematical derivation. Laplace found this solution by analyzing the 
algebraic characteristic of the L1 norm objective function; while the above procedure exactly differentiates the objective function. 
Another new contribution of this procedure is the application of the partial discrete differentiation over subscript accompanying 
with conventional differentiation on a regular variable with the continuous domain.  It should be remembered that Boscovich 
(1757) solved this problem by a geometrical procedure and Karst (1958) by an analytical method. 
 
2.1 Computation of weighted median 

Two series of computations are necessary to compute the weighted median. One sorting algorithm is essential to sort 
the z array and restoring the corresponding subscripts to compute the second part of the calculation to find the optimal value of 
k by using Dk defined by (11). 

Efficient sorting algorithms exist for the first part of the computation. The algorithms 'quicksort' of Hoare (1961,62), 
'quickersort' of Scowen (1965) and 'sort' of Singleton (1969) have desirable performances and efficiencies. For the second part 
of the computation, there is no special purpose procedure, but Bloomfield and Steiger (1980) used the partial sorting of 
Chambers (1971) to give an efficient way to combine the two steps of sorting the z array and finding the value of k. The 
superiority of this procedure is in sorting the smaller segments of z array rather than all its elements. With some modification, 
this procedure is used for our proposed algorithms in the forthcoming sections. This procedure can be stated as the following 
function. 
 
FUNCTION LWMED (n,ys,w,l) 
Step 0) Initialization. 
        Real: ys(n), w(n). 
        Integer: l(n), hi. 
        Set: ii=0, shi=0, slo=0, sz=0, sp=0, sn=0. 
Step 1) Compute left, middle and right sum of weights. 
        Do loop for i=1,n: w(i)=|w(i)|; if ys(i)<0, then sn=sn+w(i), if ys(i)>0, then sp=sp+w(i), if ys(i)=0 then sz=sz+w(i); 

end do. 
        If shi≤slo then go to step 2.b, otherwise go to step 2.a. 
Step 2) Assign subscripts for arrays. 



Copyright © CC-BY-NC 2019, BJMSR 

 

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

4 
               

     a. Let: shi=0. 
        Do loop for i=1,n: if ys(i) ≤0 go to continue, otherwise ii=ii+1, l(ii)=i, continue, end do. 
        Go to step 2.c. 
     b. Let: slo=0. 
        Do loop for i=1,n: if ys(i)>0 go to continue, otherwise ii=ii+1, l(ii)=i, continue, end do. 
     c. Let: lo=1, hi=ii. 
Step 3) Check for solution. 
        If hi>lo+1 then go to step 4, otherwise lwmed=l(lo). 
        If lo=hi return, otherwise if ys(l(lo)) ≤ys(l(hi)) go to step 3.a,  
        otherwise lt=l(lo), l(lo)=l(hi), l(hi)=lt, lwmed=l(lo). 
     a. If shi+w(l(hi))>slo+w(l(lo)) then set lwmed=l(hi), otherwise return. 
Step 4) Divide the string into two halves then sort. 
        Set: mid=(lo+hi)/2, lop=lo+1, lt=l(mid), l(mid)=l(lop), l(lop)=lt. 
     a. If  ys(l(lop)) ≤ys(l(hi))  then  go  to  step  4.b, otherwise lt=l(lop), l(lop)=l(hi), l(hi)=lt. 
     b. If ys(l(lo)) ≤ys(l(hi)) then go to step 4.c, otherwise lt=l(lo), l(lo)=l(hi), l(hi)=lt. 
     c. If ys(l(lop)) ≤ys(l(lo)) then go to step 5, otherwise lt=l(lop), l(lop)=l(lo), l(lo)=lt. 
Step 5) Compute the accumulation of weights. 
        Let: lwmed=l(lo), i=lop, j=hi, xt=ys(lwmed), tlo=slo, thi=shi. 
     a. Set: tlo=tlo+w(l(i)), i=i+1. 
        If ys(l(i))<xt then go to step 5.a, otherwise go to step 5.b. 
     b. Let: thi=thi+w(l(j)), j=j-1. 
        If ys(l(j))>xt then go to step 5.b, otherwise if j≤i then go to step 6,  
        otherwise lt=l(i), l(i)=l(j), l(j)=lt, go to step 5.a. 
Step 6) Test for solution. 
        Let: test=w(lwmed). 
        If i≠j then go to step 6.a, otherwise test=test+w(l(i)), i=i+1, j=j-1. 
     a. If test≥|thi-tlo| then return, otherwise, if tlo>thi then step 6.b, 
        otherwise slo=tlo+test, lo=i, go to step 3. 
     b. Let: shi=thi+test, lo=lop, hi=j. go to step 3. 
END 
 
2.2 Discussion of the m parameters model 

The procedure applied to the simple one-parameter model cannot be simply generalized to the m parameters model. 
This difficulty is due to the fact that we can not reorder observations in such a way that their corresponding residuals are 
increasingly or decreasingly ordered. On the other hand, to apply the discrete differentiation technique, primarily (2) should be 
rewritten as follows: 
             t         m            n            m 

     S =  (yh -  ßjxjh) -    (yh -   ßjxjh)                                                                                (12) 

          h=1      j=1        h=t+1      j=1 

Expression (12), which is free of absolute value sign is the generalization of (8) for m parameters. But our first 
problem is to find a logic which enables us to form (12). On the other hand, we need to reorder observations in such a way that 
when h is less, equal or greater than t+1, uh be greater, equal or less than zero respectively; and as h increases from one to n, the 
corresponding residual uh decreases accordingly. 

If (2) could be written as (12), we could again differentiate it with respect to j for all j=1,...,m and t. Differentiating 

S with respect to j for all j=1,...,m would give, 

     δS         t          n 

     ── = -  xjh +    xjh          j=1,...,m                                                                                 (13) 

     δßj      h=1      h=t+1 

In order to differentiate S with respect to discrete domain variable t we would have, 

     (S)                    S[t+(t)]-S(t) 



Copyright © CC-BY-NC 2019, BJMSR 

 

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

5 
               

     ─── = lim         ──────── = S(t+1)-S(t) = 

     (t)      (t)─>1         (t) 

 

    t+1     m            n         m            t+1     m             n         m                      m 

      (yh- ßjxjh) -    (uh- ßjxjh) -   (yh- ßjxjh) +    (yh- ßjxjh) =2(yt+1- ßjxj,t+1) = 0     (14) 

     h=1   j=1        h=t+2   j=1           h=1    j=1       h=t+1     j=1                  j=1 

m distinct values for t can be computed from (13) by using the following Dkj for each explanatory variable. The procedure to 
compute t is similar to what we did in (11) for the one parameter model. 
               k        n 

     Dkj ≡  xjh -    xjh        j=1,...,m; i=1,...,n                                                                         (15) 

             h=1    h=k+1 

It should be noted that these m values for t, all give a unique value for S in (12), because all of the errors 
corresponding to these t's have zero values. This becomes clear by zero residuals property of the L1 norm regression that there 
exist m points on the regression hyperplane which have zero residuals. In (12), these m points locate sequentially, because we had 
reordered observations in such a way that, when h increases from one to n, the errors decrease from greatest positive to lowest 
negative values. So, these m zero errors will locate one after another nearly in the middle of n elements sequence of (12). 
Corresponding to these m values of t, m observations are recognized. Substitute the values of these m observations into (14), m 

equations are found which can be solved simultaneously for m unknown j's. These last m equations again confirm the above-
cited property of the L1 norm regression. 

Since for the m parameters model, the ordering of residuals before computing optimal values of j's is not possible, 
another device should be adopted to compute the L1 norm estimates of the multiple linear regression parameters. 
3. Unrestricted simple linear regression 

To find the L1 norm estimates of {j,j=1,...,m} for the model (1), we shall try to propose algorithms to search for those 
observations on the optimal regression hyperplane. For m=2 and x1i=0, i=1,...,n; (1) reduces to restricted simple linear 

regression with only one parameter. The L1 norm estimate of the slope 2 for the simple model can be obtained by the algorithm 
given in the previous sections. Let us now consider a simple unrestricted linear model in which m=2 and x1i=1 for all i=1,...,n; 
namely, 

yi = ß1 + ß2x2i + ui                                                                                                      (16) 

3.1 Algorithm 1 
For the model given in (16), the L1 norm objective function S to be minimized will be, 

            n 

     S =  │yi - ß1 - ß2x2i│                                                                                                       (17) 

          i=1 

Let k1 denote a subscript which belongs to the range of one to n, and assume that k1
th observation (yk1,x2k1) is a candidate to be 

on the regression line. If this is the case, then uk1=0 and we can transfer the origin of YxX2 coordinates to the point (yk1,x2k1) 
without any loss. For this, we should rewrite all observations as deviates from the point (yk1,x2k1), 
     yi

k1
 = yi - yk1,             x2i

k1= x2i - x2k1        i=1,...,n                                                             (18) 
Where the terms with superscripts k1 denote deviated values. Rearrange the terms; 
     yi = yi

k1 + yk1,                  x2i =x2i
k1 + x2k1       i=1,...,n                                                             (19) 

              i=1,...,n 
Now, substitute (19) into (17), then, our L1 norm minimization problem can be redefined as, 
                       n 

     min: Sk1 =  |yi
k1-2x2i

k1 + (yk1-1-2x2k1)|                                                                        (20) 

    1,2         i=1 



Copyright © CC-BY-NC 2019, BJMSR 

 

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

6 
               

Since we assumed that the k1th observation is on the regression line, the expression            yk1-1-2x2k1=0. Thus (20) reduces to, 
                       n 

     min: Sk1 =  |yik1 - 2x2i
k1|                                                                                                 (21) 

      2            i=1 
The solution of the optimization problem (21), is the same as the one described before for one parameter linear 

model. Note that when uk1=0, then the minimum value of (21) is equal to the minimum value of (17). 

Let 2 derived by minimizing (21) be denoted by 2
k1. By changing k1 from one to n and minimizing (21), 2

1,...,2
n 

are attained accordingly. Now the question is: what value of k1 minimizes (17)? In other words, which observations are on the 
optimal regression line? Note that in the two parameters L1 norm linear model, there exists at least two observations with zero 
errors. Transferring of YxX2 coordinates to both these points does not change the minimum of (17). Suppose p th and qth are 
those observations on the regression line, thus, 

     up = yp - 1 - 2x2p = 0,         uq = yq - 1 - 2x2q = 0                                                         (22) 
Rewriting (18) and (19) for k1=p,q and substituting them into (21), yields,  

              n 

     Sp =  |yi
p - 2x2i

p|                                                                                                             (23) 
            i=1 
             n 

     Sq =  |yi
q - 2x2i

q|                                                                                                             (24) 
            i=1 
Using (18) and rewriting (23) and (24), it is shown that Sp is equal to Sq,  
              n 

     Sp =  |yi-2x2i - (yp-2x2p)|                                                                                               (25) 
            i=1 
             n 

     Sq =  |yi-2x2i - (yq-2x2q)|                                                                                               (26) 
            i=1 
Sp is equal to Sq if and only if the two parentheses inside the absolute value signs of (25) and (26) are equal. This can be 

concluded by solving the two equations of (22) for 1. That is, 

     1 = yp - 2x2p = yq - 2x2q                                                                                                 (27) 
The equality of Sq and Sp is because the points p and q are on the regression line and therefore up=uq=0. Thus, Sp=Sq, so, the 

values of 2
p and 2

q derived by minimizing either Sp or Sq must be equal. This gives a criterion to find the desired 2 from all 

2
k1 for k1=1,...,n. That is, when 2

p=2
q. The estimated value of 2 is denoted by 2^. Value of 1^ is simply computed by 

(27). Now let us summarize the whole procedure for finding the values of 1^ and 2^ for the model given by (16). 
 
Crude algorithm 1 
Step 0) Set k1=1. 
Step 1) Compute (18). 
Step 2) Minimize (21) by using the weighted median procedure and find ß2

k1. 
Step 3) Check if ß2

k1=ß2
k1-k for 0<k<k1, then set ß2^=ß2

k1 and ß1^=yk1-ß2^x2k1 then stop. 
Step 4) Increase k1 by one and go to step 1. 
 
3.2 Algorithm 2 

The crude algorithm 1 for finding ß1^ and ß2^ is not computationally efficient because it usually requires testing the 
majority of observations. In order to make this algorithm efficient, some elaborations are necessary.  

Instead of setting k1=1, let us set k1 equal to an arbitrary integer value such as a where a is any integer from one to n. 
Suppose now, ua=ya-ß1-ß2x2a=0, and rewrite (21) as, 
                     n 

     min: Sa = Σ │yi
a - ß2x2i

a│                                                                                                 (28) 
      ß2          i=1 
Minimizing (28) gives ß2

a which is equal to  
                yb

a      yb - ya 

     ß2
a = ─── = ─────                                                                                                       (29) 

               x2b
a      x2b - x2a 



Copyright © CC-BY-NC 2019, BJMSR 

 

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

7 
               

On the other hand, by pivoting on the ath observation, another point such as b is found where b refers to that observation which 
has zero error in the minimum solution of (28), so,        ub=yb-ß1-ß2x2b=0. Let us denote the minimum of Sa in (28) as Sa

*, then, 
              n 

     Sa
* = Σ │yi

a - ß2
ax2i

a│                                                                                                       (30) 
             i=1 
Note that, 
     ß1 = ya - ß2

ax2
a = yb - ß2

ax2b                                                                                               (31) 
This comes from multiplying (29) by its denominator and rearranging the terms. Substitute (31) in (17), 
            n                                     n 

     S = Σ │yi-ya-ß2
a(x2i-x2

a)│ = Σ │yi-yb-ß2
a(x2i-x2a)│ 

          i=1                                  i=1 
or, 
            n                           n 

     S = Σ │yi
a-ß2

ax2i
a│ = Σ │yi

b-ß2
ax2i

b│                                                                               (32) 
          i=1                        i=1 
Using (30) the first sum in (32) is Sa

* and the second sum is Sb evaluated at ß2=ß2
a. Thus, it can be concluded that, 

     Sa
* = Sb│                                                                                                                           (33) 

                  │ß2=ß2
a 

Sa
* is at the minimum but Sb│      still can be minimized for other values of ß2.  

                                             │ß2=ß2
a 

Therefore, an important result is derived, that is, 
     Sa* > Sb*                                                                                                                           (34) 
The inequality of (34) guarantees that if we choose an arbitrary point to transfer the origin of coordinates to it and minimizing 
the objective function (21) another point is found, then transferring the origin to this newly found point decreases the total sum 
of absolute errors. Therefore, at each transference we get closer to the minimum of S. By a similar discussion, this conclusion can 
be generalized as,  
     Sa* > Sb* > Sc* > Sd* ...                                                                                                    (35) 
Note that a is an arbitrary starting point. The point b is derived by minimizing Sa, c is derived by minimizing Sb, and d is derived 
by minimizing Sc, and so on. 

Now, the question is, when the minimum value of S is reached? Suppose S*=Sf*; by transferring the origin of 
coordinates to the point f and minimizing Sf, the gth observation is derived. When S*=Sf* by minimizing Sg, fth observation is 
again found, because Sg*=Sf*=S* and the gth and fth observations are both on the L1 norm regression line.  This conclusion can be 
a criterion to stop the procedure. Hence, 
     Sa* > Sb* > Sc* > Sd* ... > Sf* = Sg* = S*                                                                        (36) 

It should be noted that if the minimum solution of S is not unique, that is function S has a horizontal segment, the 
procedure stops when it reaches the first minimum solution. 

Now, let us introduce the whole stages of this algorithm to find the L1 norm estimates of ß1 and ß2 in the simple linear 
model (16). 
 
 
Algorithm 2 
Step 0) Select an arbitrary integer value a between one to n and set k1=a. 
Step 1) Compute (18) with k1=a. 
Step 2) Minimize (21) by using weighted median procedure and find the appropriate observation on the line, namely 

observation b. 
Step 3) Compute (18) with k1=b. 
Step 4) Minimize (21) and find that observation on the line; observation c. 
Step 5) Check if c=a, then ß2^=yb

c/x2b
c, ß1^=yc-ß2^x2c and stop. 

Step 6) Set a=b and go to step 1. 
 

Operationally, the following steps are taken. 
 
PROGRAM BL1S 
Step 0) Initialization. 
        Parameter: n. 
        Real: y(n), x2(n), z(n), w(n). 



Copyright © CC-BY-NC 2019, BJMSR 

 

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

8 
               

        Integer: l(n). 
        Set: k1=arbitrary, k1r=0, k1s=0, iter=0. 
        Read (y(i), x2(i), i=1,n) 
Step 1) Compute weights and ratios. 
        Do loop for i=1,k1-1: w(i)=x2(i)-x2(k1), z(i)=(y(i)-y(k1))/w(i), end do. 
        Set: w(k1)=0, z(k1)=0. 
        Do loop for i=k1+1,n: w(i)=x2(i)-x2(k1), z(i)=(y(i)-y(k1))/w(i) end do. 
        Set: iter=iter+1. 
Step 2) Compute weighted median. 
        Let: lm=LWMED(n,z,w,l). 
Step 3) Check for optimality. 
        Set: k1s=k1r, k1r=k1. 
        If lm=k1s then go to step 4, otherwise k1=lm. 
        Go to step 1. 
Step 4) Compute the solution. 
        Let b2=z(lm), b1=y(k1)-b2*x2(k1). 
        Print b1, b2, k1, lm, iter. 
        Stop. 
EN 
4. Generalization to m parameters 

Now we extend the above procedure to the restricted two parameters model, namely, 
     yi = ß2x2i + ß3x3i + ui                                                                                                         (37) 
Let, 
            n 

     S = Σ │yi - ß2x2i - ß3x3i│                                                                                                   (38) 
          i=1 
S can be written as, 
            n                                                    n 

     S = Σ │x2i││yi/x2i - ß2 - ß3x3i/x2i│ =  Σ │x2i││yi
s1 - ß2 - ß3x3i

s1│                                  (39) 
          i=1                                                 i=1 
where, 
     yi

s1 = yi/x2i ,         x3i
s1 = x3i/x2i                i=1,...,n                                                              (40) 

Minimization of (39) is similar to that of simple linear model explained in the previous section. An important distinction for 

solving (39) in comparison to (17) is the expression │x2i│ which has been multiplied to │yi
s1-ß2-ß3x3i

s1│. This multiplication 
does not make any problem when (39) is minimized, because if we deviate y i

s1 and x3i
s1 from the point k1 as, 

     yi
s1k1 = yi

s1 - yk1
s1,              x3i

s1k1 = x3i
s1 - x3k1

s1             i=1,...,n                                       (41) 
then we can rewrite (39) similar to (21), so, 
               n 

     Sk1 = Σ │x2i││yi
s1k1 - ß3x3i

s1k1│                                                                                       (42) 
              i=1 
To minimize (42), we should use the following expression, 
               n 

     Sk1 = Σ │x2ix3i
s1k1││yi

s1k1/x3i
s1k1 - ß3│                                                                             (43) 

             i=1 
According to Bidabad (1987a,88a), in applying discrete derivative to (43) the term in the first absolute value sign is used to find 
the subscript of that point which locates on the regression line.  In comparison to the simple unrestricted linear model (16), this 
is the main difference. 

In the case of inclusion of an intercept in the model given by (37), we have, 
            n 

     S = Σ │yi - ß1 - ß2x2i - ß3x3i│                                                                                            (44) 
          i=1 
Let k1 be an arbitrary subscript then, the transference of the origin of coordinates to the k1 th point is done by deviating all 
observations from this point, namely, 
     yi

k1 = yi - yk1,             x2i
k1 = x2i - x2k1,            x3i

k1 = x3i -x3k1          i=1,...,n                                 (45) 
Rearrange the terms of (45) and substitute in (44), then we have, 
            n 



Copyright © CC-BY-NC 2019, BJMSR 

 

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

9 
               

     S = Σ │yi
k1-ß2x2i

k1-ß3x3i
k1+(yk1-ß1-ß2x2k1-ß3x3k1)│                                                           (46) 

          i=1 
If the k1th observation is on the regression plane, then 
     yk1 - ß1 - ß2x2k1 - ß3x3k1 = 0                                                                                               (47) 
So, instead of minimizing (44), the following function is minimized: 
               n 

     Sk1 = Σ │yi
k1 - ß2x2i

k1 - ß3x3i
k1│                                                                                        (48) 

             i=1 
Minimization of (48) is completely similar to that of (38) and can be proceeded as follows, 
               n                                                                 n 

     Sk1 = Σ │x2i
k1││yi

k1/x2i
k1 – ß2 - ß3x3i

k1/x2i
k1│ = Σ │x2i

k1││yi
s1 - ß2 - ß3x3i

s1│                (49) 
             i=1                                                              i=1 
where, 
     yi

s1 = yi
k1/x2i

k1,         x3i
s1 = x3i

k1/x2i
k1         i=1,...,n                                                                 (50) 

Now, again, transfer the origin of the two dimensional Ys1xX3
s1 coordinates to an arbitrary point k2 by deviating yi

s1 and x3i
s1 

from yk2
s1 and x3k2

s1 as follows, 
     yi

s1k2 = yi
s1 - yk2

s1 ,                 x3i
s1k2 = x3i

s1 - x3k2
s1             i=1,...,n                                   (51) 

By rearranging the terms of (51) and substituting them into (49) and assuming that the point k2 is on the regression plane, we 
can rewrite (49) as, 
                 n 

     Sk1k2 = Σ │x2i
k1││yi

s1k2 - ß3x3i
s1k2│                                                                                  (52) 

               i=1 
or, 
                  n 

     Sk1k2 = Σ │x2i
k1x3i

s1k2││yi
s1k2/x3i

s1k2 - ß3│                                                                       (53) 
                i=1 

The objective function (53) can be minimized, as suggested before. Now, the procedure from (49) to (53) can be 
repeated with different values of k2 as in algorithm 2 proposed for the simple linear model. When the last point (M) in the 
process of minimizing (53) is reached, the origin of three dimensional YxX2xX3 coordinates (k1) is transferred to this newly 
found point M and the procedure from (45) to (53) is again repeated with the exception that instead of assigning an arbitrary 
value to k2, we set k2 equal to the previous value of k1. This procedure continues until the point found by minimizing (53) is 
equal to the previous value of k1. The values of ß1^, ß2^ and ß3^ can be computed according to the following formulas,  
     ß3^ = yM

s1k2/x3M
s1k2 ,      2^ = yk2

s1 - ß3^x3k2
s1,           ß1^ = yk1 - ß2^x2k1 - ß3^x3k1

                          (54) 
It should be noted that at each step, it can be proved that we are getting closer to the minimum of S in (44). The 

proof is similar to that of the two parameters model given by (16). 
It can be shown that for any two arbitrary points k1 and k2, the value of Sk1k2 given by (53) for any commutation for 

k1 and k2 remains unchanged. This will be shown in the next sections. However, 
     Sk1k2 = Sk2k1                                                                                                                       (55) 
Now, following the procedure of the proof given by (28) through (36), we can write, 

     Sk1a* = Sk1b│                                                                                                                     (56) 

                        │ß3=ß3
k1a 

where Sk1a* denotes the optimal value of Sk1a with respect to ß3 of (53) for any arbitrary integer a from one to n, namely, 
     Sk1a* = min (Sk1a)                                                                                                              (57) 
                   ß3 
Subscript b in (56) denotes the observation which is found by minimizing Sk1a. Since the ath and bth observations both are on the 
regression hyperplane, by a similar discussion stated by (28) through (32) we can conclude that Sk1a* is equal to Sk1b which 
evaluated at that value of ß3 which is found by minimizing Sk1a. Since the left-hand side of (56) is minimum and its right-hand 
side can be decreased by other values of ß3, it can be concluded that, 
     Sk1a* > Sk1b*                                                                                                                      (58) 
According to (55) we can write, 
     Sk1a* > Sk1b* = Sbk1*                                                                                                         (59) 
Since in each step we discard an observation from the basis and replace it with the newly found one, we have, 
     Sk1a* > Sk1b* = Sbk1* > Sbc* = Scd*                                                                                   (60) 
or generally, 
     Sk1a* > Sk1b* > Sbc* > Scd*                                                                                                (61) 



Copyright © CC-BY-NC 2019, BJMSR 

 

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

10 
               

The minimum solution S* is reached when by entering the f th observation in the basis we can not reduce the objective function 
value. That is the previous observation is again reached, namely, 
     Sk1a* > Sab* > Sbc* > Scd* > ... > Sfg* = Sgf* = S*                                                            (62) 
The relation (62) guarantees that at each step, we are descending down the objective function surface. 
 
 
4.1 Algorithm 3 
     To generalize the above algorithm to the m parameters linear model (1), we should reduce the number of parameters in the 
same fashion as in the case of three parameters model explained above. If the model is restricted, we can make it unrestricted by 
dividing all dependent and independent variables to one of the independent variables as follows. 
            n         m            n                              m 

     S = Σ │yi- Σ ßjxji│= Σ │x2i││yi/x2i -ß2- Σ ßjxji│                                                             (63) 
           i=1      j=2          i=1                          j=3 
If the model is unrestricted, we can make it restricted by deviating all observations from an arbitrary one observation; namely, 
            n              m            n             m 

     S = Σ │yi-ß1- Σ ßjxji│= Σ │yi
k1 - Σ ßjxji

k1│                                                                      (64) 
          i=1           j=2          i=1          j=2 
where, 
     yi

k1 = yi - yk1,             xji
k1 = xji - xjk1      i=1,...,n; j=2,...,m                                               (65)                                                                                                                                              

Therefore, according to the transformations (63), (64) and (65), any m parameters model can be reduced to a simple 
restricted one parameter model and then be estimated as a weighted median problem. To conduct this transformation, if the 
model is unrestricted, we should deviate all observations from an arbitrary one and make the model restricted. Then divide all 
dependent and independent variables to an arbitrary independent variable. This makes the model unrestricted. At this step, we 
have reduced the number of parameters by one. By continuing this process, we can reduce any m parameters model to a one-
parameter model. If the model is restricted, we should start by dividing all of the variables by one of the independent variables 
and make the model unrestricted. Now, we can again reduce one of the parameters by applying the procedures explained above 
in order to transform the unrestricted model to restricted form. By solving the one parameter model, we can then solve for the 
two parameters and then three parameters models and so on. 

The most delicate part of this algorithm is that, at the starting point, once k1,k2,...,k(m-1) are selected arbitrarily, then 
the algorithm assigns the best possible values to the integers k1,k2,...,k(m-1). To explain the procedure, let us deal with the four 
parameters unrestricted linear model. Once, the value of k1 is arbitrarily selected.  Deviate all observations from k1 th observation. 
In this way, the model reduces to a three parameters restricted model. By dividing all of the variables to one of the independent 
variables, our model becomes completely similar to (44). By minimizing (44) according to the algorithm previously explained 
for the three parameters model, subscript M corresponding to Mth point is derived. This is the newly found point which its 
subscript (M) is assigned to k1. The previous value of k1 is assigned to k2, and the previous value of k2 is assigned to k3, and 
the whole procedure is repeated again. The procedure stops when the value of M is equal to k1.  
     The above important assigning technique, which is essential for pivoting on the origins of different size coordinates, can be 
extended to more parameters models as we did before. Again, we should note that at each succeeding step, we get closer to the 
minimum of S in (1). The proof is completely similar to the three parameters model explained before. The whole procedure is as 
follows. 
 
Algorithm 3 
Step 0) Select arbitrary points k1,k2,...,k(m-1). 
Step 1) Set counter=1. 
Step 2) Proceed the transformations (63) through (65) to reduce the model to a one-parameter transformed restricted linear 

model. 
Step 3) Apply weighted median procedure to find the point M. 
Step 4) If M is equal to the previous value of k(counter) go to step 5; otherwise set k(counter)=M and go to step 1. 
Step 5) If counter = m-1 go to step 6; otherwise let counter=counter+1 and go to step 4. 
Step 6) Set k(m)=M. By using the final values of k1,...,k(m) solve the following system of equations for ß j^, for j=1,...,m; 
             m 

             Σ ßj^xjki = yki          i=1,...,m 
            j=1 
           Stop. 
 
 



Copyright © CC-BY-NC 2019, BJMSR 

 

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

11 
               

4.2 Algorithm 4 
In algorithm 3 to solve any m parameters model, the reduced m-1 parameters model should be solved first and 

therefore makes the algorithm rather costly. Because to solve a m parameters model, we should reduce the objective function 
value of the reduced m-1 parameters model down to its minimum value and then going back to the m parameters model. In 
interior steps for m-1,m-2,...,1 this procedure repeats. In order to make this algorithm more efficient, another assigning 
technique may be adopted. 

Once m-1 observations are selected arbitrarily. The S function is reduced to a one-parameter weighted median 
problem as we did in algorithm 3 using (63) and (64). By solving the weighted median problem, a new observation M is found. 
Now, the newly found point is replaced by the most previous observation which entered the basis. The procedure is repeated 
until we can not find any new observation out of the current set of basis points and the deleted point in the previous iteration. 
On the other hand, in this algorithm m-1 point is selected as a basic feasible solution. By pivoting on this solution, a new point 
is found to enter the basis and should be replaced by another one which was previously on the basis. The procedure ends when 
we cannot enter any new point to reduce the objective function value. In this context, it is similar to the simplex procedure of 
linear programming technique. 

Algorithm 4 is also a descent method. Because in each iteration the value of S is decreased. The proof is completely 
similar to that of algorithm 3 stated before. The main difference of this algorithm with algorithm 3 is in choosing the point 
which should be deleted from the basis.  Therefore, in each step, we need not keep some points in the basis to solve smaller size 
models. On the other hand, in each m iterations, all m points in the basis are discarded, and new points enter. In algorithm 3, we 
kept some points in the basis until we reach an optimal point for smaller size reduced model. It should be noted that one point 
may be discarded and reenter the basis through m iterations. This is not contradictory and actually occurs, especially when we 
are near the optimal solution.  

As we will see in the forthcoming section, due to different properties of this algorithm, many manipulations in 
computational steps may be adopted to start calculations. Before making this algorithm more sophisticated, it is suitable to 
summarize its general steps as follows. 
 
Algorithm 4 
Step 0) Select arbitrary points k1,...,k(m-1). 
Step 1) Proceed the transformations (63) through (65) to reduce the model to one-parameter transformed restricted linear 

model. 
Step 2) Apply weighted median procedure to find the point M. 
Step 3) Assign M to one of the k1,...,k(m-1) which is oldest in the basis. 
Step 4) Check if the newly found point M which enters the basis, discard the point which entered in the previous iteration 

repeats m-1 times then go to step 5; otherwise go to step1. 
Step 5) Set k(m)=M. By using the final values of k1,...,k(m), solve the following system of equations for ß j^, for j=1,...,m; 
             m 

            Σ ßj^xjki = yki          i=1,...,m 
           j=1 
        Stop. 
 

After discussing the properties of algorithm 4, we will give more expository steps of computation for this algorithm.  
By using these properties, we can make this algorithm more efficient. 
 
4.3 properties 

Before discussing the properties of the proposed algorithms specifically algorithm 4, let us consider the following four 
parameters model L1 norm objective function, 
            n 

     S = Σ │yi-ß0-ß1x1i-ß2x2i-ß3x3i│                                                                                         (66) 
           i=1 
Some steps should be taken to reduce (66) to a weighted median problem by entering three arbitrary integers k1, k2, k3 between 
one and n as follows. Note that function (1) has been slightly redefined to include intercept as a simple term, in equation (66). 
Similar to the transformations (63) through (65), the following steps are taken sequentially. 

Deviate all observations from k1th observation, 
     yi

k1 = yi - yk1         i=1,...,n                                                                                                 (67) 
     xji

k1 = xji - xjk1      i=1,...,n; j=1,2,3                                                                                    (68) 
The objective function becomes, 
               n 



Copyright © CC-BY-NC 2019, BJMSR 

 

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

12 
               

     Sk1 = Σ │yi
k1-ß1-ß2x2i

k1-ß3x3i
k1│                                                                                       (69) 

             i=1 
Divide yik1 and xji

k1 for j=2,3 by x1i
k1 for all i=1,...,n; 

     yi
s1 = yi

k1/x1i
k1       i=1,...,n                                                                                                (70) 

     xji
s1 = xji

k1/x1i
k1     i=1,...,n; j=2,3                                                                                      (71) 

The objective function becomes, 
               n 

     Sk1 = Σ │x1i
k1││yi

s1-ß1-ß2x2i
s1-ß3x3i

s1│                                                                           (72) 
             i=1 
Deviate yi and xji

s1 for j=2,3 from k2th observation, 
     yi

k2 = yi
s1 -yk2

s1       i=1,...,n                                                                                               (73) 
     xji

k2 = xji
s1 - xjk2

s1   i=1,...,n, j=2,3                                                                                     (74) 
The objective function reduces to, 
                  n 

     Sk1k2 = Σ │x1i
k1││yi

k2-ß2x2i
k2-ß3x3i

k2│                                                                             (75) 
               i=1 
Divide yi

k2 and xji
k2 for j=3 by x2i

k2 for all i=1,..,n; 
     yi

s2 = yi
k2/x2i

k2       i=1,...,n                                                                                                (76) 
     xji

s2 = xji
k2/x2i

k2     i=1,...,n; j=3                                                                                         (77) 
The objective function forms as, 
                 n 

     Sk1k2 = Σ │x1i
k1x2i

k2││yi
s2-ß2-ß3x3i

s2│                                                                             (78) 
               i=1 
Deviate yi

s2 and xji
s2 from k3th observation, 

     yi
k3 = yi

s2 -yk3
s2       i=1,...,n                                                                                               (79) 

     xji
k3 = xji

s2 - xjk3
s2   i=1,...,n, j=3                                                                                        (80) 

Rewrite the objective function as, 
                  n 

     Sk1k2 = Σ │x1i
k1x2i

k2x3i
k3││yi

s3-ß3x3i
k3│                                                                           (81) 

                i=1 
Divide yi

k3 by x3i
k3 for all i=1,...,n; 

     yi
s3 = yi

k3/x3i
k3        i=1,...,n                                                                                               (82) 

Finally the objective function becomes, 
                    n 

     Sk1k2k3 = Σ │x1i
k1x2i

k2x3i
k3││yi

s3-ß3│                                                                               (83) 
                  i=1 
Rewrite (83) as, 
                    n 

     Sk1k2k3 = Σ │wi
k1k2k3││ri

s1s2s3-ß3│                                                                                    (84) 
                  i=1 
where, 
     wi

k1k2k3 = x1i
k1x2i

k2x3i
k3                                                                                                      (85) 

     ri
s1s2s3 = yi

s3                                                                                                                        (86) 
The objective function (84) is a weighted median problem which may easily be minimized to find ß3. Equation (84) for the 
general m parameters model  
           m-1 

     yi = Σ xjißj + ui        i=1,...,n                                                                                              (87) 
          j=0 
     x0i=1                    i=1,...,n                                                                                                  (88) 
with the corresponding objective function, 
            n         m-1 

     S = Σ │yi - Σ xjißj│                                                                                                           (89) 
           i=1      j=0 
is,                        n 

     Sk1k2...k(m-1) =  Σ │wi
k1k2...k(m-1)││ri

s1s2...s(m-1)-ßm-1│                                                           (90) 
                          i=1 



Copyright © CC-BY-NC 2019, BJMSR 

 

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

13 
               

where, 

     wi
k1k2...k(m-1) = │x1i

k1x2i
k2...x(m-1)i

k(m-1)│                                                                             (91) 
     ri

s1s2...s(m-1) = yi
s(m-1)                                                                                                            (92) 

Now, we proceed to explain the properties of algorithm 4 by using the formulation described above. 
 
 
Property 1 

If we start the computation with m-1 points which do not belong to the original set of observations, the value of S 
decreases to its global minimum at each iteration. Suppose the points (yn+1,...,xm-1,n+1),...,(yn+m-1,...,xm-1,n+m-1) do not belong to our 
sample points. If n+1,n+2,...,n+m-1 are assigned to k1,k2,...,k(m-1), then minimizing the function in (90) leads to finding a 
new point M which belongs to the sample observations. Replacing k1 by M and deleting the (n+1)th point from the basis and 
minimizing (90) leads to finding another point M' in the sample which will be replaced by the k2 th observation which is out of 
our sample. Thus, in each iteration one of the out of sample points is deleted, and optimization procedure will be done on the 
sample points after m-1 iterations. 

This property has a good performance on ill-conditioned data. That is if for example, rounding error causes an 
observation to get incorrect rounded value, then this point in the next iteration will be replaced by a sample point and thus 
redirect the path of descending to the right points. 
 
Property 2 

Reduction of m-1 parameters model to fewer parameters model can be made by using the computation of m-1 
parameters model and conversely. Thus stepwise and all possible regressions may be computed efficiently.  

To show this property, let us deal with equation (86). Value of ri
s1s2s3 may be computed by substituting (67), (68), 

(70), (71), (73), (74), (76), (77), (79), (80) and (82) into (86). The result is, 

1 2 1 3 1 2 1

1 1 1 1 2 1 1 1 3 1 1 1 2 1 1

2 2 1 2 2 2 1 2 3 2 1 2 2 2 1

1 1 1 1 2 1 1 1 3 1 1 1 2 1 11 2 3

3 3 1 3 2 3 1

1 1 1 1 2 1 1

2 2 1 2

1 1 1

i k k k k k k k

i k k k k k k k

i k k k k k k k

i k k k k k k ks s s

i
i k k k

i k k k

i k

i k

y y y y y y y y

x x x x x x x x

x x x x x x x x

x x x x x x x x
r

x x x x

x x x x

x x x

x x

    
   


   

 
   


 


 






3 3 3 1 3 2 3 1

1 3 1 1 1 2 1 1

2 2 1 2 3 2 1 2 2 2 1

1 2 1 1 1 3 1 1 1 2 1 1

k k k k

k k k k

k k k k k k

k k k k k k

x x x x

x x x x

x x x x x

x x x x x x

 


 


  


  

                                          (93) 

 
By inspection, it is obvious that the contents of rectangle grid lines in (93) are ri

s1, ri
s1s2 and ri

s1s2s3. Generally, ri
s1s2...s(m-1) has the 

following scheme, 

1

1 2

1... ( 1)

1

1 2... ( 1)

s

i

s s

s s m i

s s s m

i

r

r
r

r






                                            (94) 

 
A similar discussion may also be performed for wi

k1k2...k(m-1). Rewrite (91) by substitution of (68), (74) and (80) in (91), 



Copyright © CC-BY-NC 2019, BJMSR 

 

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

14 
               

 

1 2 3

3 3 3 1 3 2 3 13 3 1 3 2 3 1

2 2 1 2 2 2 1 1 3 1 1 1 2 1 11 1 1 1 2 1 1
1 1 1

2 2 1 2 2 2 1 2 3 2 1 2 2 2 11 1 1 1 2 1 1

1 1 1 1 2 1 1 1 3 1 1 1 2

k k k

i

k k k ki k k k

i k k k k k k ki k k k
i k

i k k k k k k ki k k k

i k k k k k k

w

x x x xx x x x

x x x x x x x xx x x x
x x

x x x x x x x xx x x x

x x x x x x x x



  


     
           

    1 1k

 
 
 
 
 
 

 

                                                                                                                                               (95) 
By inspection of (93), (94) and (95) it is obvious that calculation of r i

s1s2...sh and wi
s1s2...sh where 0<h≤m-1 which are necessary for 

h parameters regression may be derived from the m parameters model steps of computation. 
 
Property 3 

By any commutation for kj for j=1,...,m-1, values of ri
s1s2...s(m-1) defined by (92) for all i≠k1,k2,...,k(m-1) remain 

unchanged. This calculation can be made by deleting the suitable terms in (93) and (94) after taking common denominator. 
Because of this property, we can commute for the values of k1,...,k(m-1) in the process of computation without changing the 
value of ri

s1...s(m-1). 
 
Property 4 

By any commutation for kj for j=1,...,m-1, values of wi
k1...k(m-1) for all i=1,...,n defined by (91) remain unchanged. This 

conclusion is made by multiplication of sequential parentheses of (97) and deleting the suitable terms. This property 
accompanying with property 3 guarantees that any commutation for k1,...,k(m-1) does not change the necessary items of (91) 
and (92) in the computation of weighted median problem (90). 
 
Property 5 

Value of rk(m-1)
s1...s(m-1) is equal to zero.  By substituting i=k(m-1) in (93) or (94) this property is concluded. Since, this 

term is equal to zero we should be aware of divide check (dividing by zero) in the process of computation of rk1
s1,rk2

s1s2,...,rk(m-

1)
s1s2...s(m-1). 

 
Property 6 

Values of wh
k1k2...k(m-1)=0 for all h=k1,k2,...,k(m-1). If h is equal to one of the k1,k2,...,k(m-1), then one of the 

parentheses of (95) is equal to zero, and this property is assured. According to this property to avoid divide check, without 
computing wh

k1k2...k(m-1) for h=k1,k2,...,k(m-1) we can set these terms equal to zero. 
 

Property 7 
By any commutation for kj for j=1,...,(m-1), the Mth point derived by minimizing (89) remains unchanged. Because, by 

this commutation values of ri
s1s2...s(m-1) and wi

k1k2...k(m-1) remain unchanged (properties 3,4,5 and 6). 
 
Property 8 

Updating procedure to evaluate different subsets of sample observations can be simply done by assigning zero values to 
their wh

k1k2...k(m-1) in (90) for deleting points.  This makes those unwanted points h's be unaffected in the calculation. 
 
Property 9 

According to properties 6 and 8, since the value of wh
k1k2...k(m-1)=0 for all h=k1,k2,...,k(m-1); the point M which derives 

from solving weighted median problem (90) will not be equal to k1,k2,...,k(m-1). This property guarantees that when we enter 
an observation into the basis, we do not find this point again in the same iteration. 
 
Property 10 

If all kj for j=1,...,m-1 are equal in starting time, no problem occurs in the process of computation, but the execution 
time may be increased. Conversely, if the kj points are selected in such a way that they are near the minimum of S, the algorithm 
converges to the optimal solution with fewer iterations. 
 
Property 11 

Values of wi
k1k2...k(m-1) and ri

s1s2...s(m-1), which are essential to solving the weighted median problem of (90) can be retrieved 
from their values in the previous iteration. That is when the point kh is discarded from the basis and replaced by kj, values of 



Copyright © CC-BY-NC 2019, BJMSR 

 

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

15 
               

wi
k1k2...k(h-1) and ri

s1s2...s(h-1) do not change and need not be computed again. Applying this property to algorithm 4 makes it more 
efficient. 
 
Property 12 

Since at each iteration, values of ri
s1s2...s(m-1) and wi

k1k2...k(m-1) should be computed from the source input data y and X, 
rounding error does not accumulate, and convergence does not perturb. This property makes the algorithm safe from the 
accumulation of rounding error, which exists in simplex-type algorithms. 
 
Property 13 

When m optimal points are reached, the optimal values of ß j^ for j=1,...,m can be computed as a recursive system of 
equations efficiently. Values of ßj^ will be, 

 

( 1)

1

( 2)( 2)
( 1) ( 1)2 ( 1)

2 22
( 1) 3 ( 2) 32 3

1 1 11
( 1) 2 ( 2) 2 2 21 2

( 1) 1 ( 2)0 1

0 0 0 0 0

0 0 0 0

0 0 0

0 0

s m

m M

s ms m
m k mm k m

s ss
m k m kk

s s ss
m k m k kk

m k m kk

y

xy

x xy

x x xy

x xy











 



 
  


 


 


 

   
   
   
   

    
   
   
   

  
  

1

2

2

1

1 2 1 1 1 0 0

m

m

k kx x





















  
  
  
  
  
  
  
  

  
  

  (96) 

The recursive system of equations (96) can be obtained as follows. Solve (67) and (68) for y i and xji, j=1,2,3; and substitute 
them into (66), gives 
               n 

     Sk1 = Σ │yi
k1-ß1x1i

k1-ß2x2i
k1-ß3x3i

k1-(yk1-ß0-ß1x2k1-ß3x3k1)│                                              (97) 
              i=1 
Since we force the regression to pass through the point k1, the expression inside the parentheses of (97) is equal to zero. Hence, 
     yk1-ß0-ß1x2k1-ß3x3k1 = 0                                                                                                     (98) 
We can solve (73) and (74) again for yi

s1 and xji
s1, where, j=2,3 and substitute them into (72) which leads to the following 

objective function, 
             n 

     Sk1k2 = Σ │x1i
k1││yi

k2-ß2x2i
k2-ß3x3i

k2-(yk2
s1-ß1-ß2x2k2

s1-ß3x3k2
s1)│                                   (99) 

            i=1                                                      
By a similar discussion on (97) we have, 
     yk2

s1-ß1-ß2x2k2
s1-ß3x3k2

s1 = 0                                                                                            (100) 
Similarly we can derive, 
               n 

     Sk1k2k3 = Σ │x1i
k1x2i

k2││yi
k3-ß3x3i

k3+(yk3
s2-ß2-ß3x3k3

s2)│                                               (101) 
              i=1                                                      
and based on (101), 
     yk3

s2-ß2-ß3x3k3
s2 = 0                                                                                                         (102) 

Finally, by solving the weighted median problem (83), we find a point M which gives, 
     ß3 = yM

s3                                                                                                                          (103) 
Equations (103), (102), (100) and (98) can be solved sequentially for ß3^,...,ß0^. In a similar fashion, the general solution for m 
parameters model is given by (96). 

Now let us use these properties and give the computational steps of algorithm 4 more expository in the following BL1 
program. 
 
PROGRAM BL1 
Step 0) Initialization. 
        Parameter: n, m, m1=m-1, m2=m-2. 
        Real: y(n), x(n,m1), xsk(m1), yw(n), xkw(n), w(n), ys(n), xs(n,m1), b(m), xw(n,2:m1), ysol(m1), xsol(m1,m1). 
        Integer: l(n), kk(m1). 
        Common: /c1/i1,i2. 
        Read: (y(i),(x(i,j),j=1,m1),i=1,n). 



Copyright © CC-BY-NC 2019, BJMSR 

 

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

16 
               

        Let: iter=0, kr=0, mm=1, (kk(j)=arbitrary,j=1,m1). 
Step 1) Refill working arrays. 
        Do loop for i=1,n: ys(i)=y(i), do loop for j=1,m1: xs(i,j)=x(i,j), end do, end do. 
Step 2) Store weights and ratios for next iteration. 
        Do loop for i=1,n: w(i)=xkw(i), ys(i)=yw(i), do loop for j=1,m1: xs(i,j)=xw(i,j), end do, end do. 
Step 3) Compute the arguments for weighted median. 
     a. Set: jj=mm, k=kk(jj), ysk=ys(k), i1=1, i2=k-1. 
        Do loop for j=jj,m1: xsk(j)=xs(k,j), end do. 
     b. Do loop for j=jj,m1: call COL1(xsk(j),xs(1,j)) end do. 
        Call COL2(ysk,jj,w,ys,xs(1,jj)). 
        If i2=n go to step 3.c; otherwise set: i1=k+1, i2=n, go to step 3.b. 
     c. Set: w(k)=0. 
        If jj=m1 go to step 4; otherwise i1=1, i2=k-1, go to step 3.d. 
     d. Do loop for j=jj+1,m1: call COL3(xs(1,j),xs(1,jj)), end do. 
        If i=n go to step 3.e; otherwise i1=k+1, i2=n, go to step 3.d. 
     e. If jj≠mm jj=jj+1, go to step 3, otherwise do loop for i=1,n: xkw(i)=w(i), yw(i)=ys(i); do loop for j=jj+1,m1: 

xw(i,j)=xs(i,j), end do; end do. 
        Set: jj=jj+1, go to step 3. 
Step 4) Compute the weighted median. 
        Set: ys(k)=0, iter=iter+1, lm=LWMED(n,ys,w,l). 
Step 5) Test for optimality. 
        If lm=kr go to step 5.b; otherwise iopt=0 go to step 5.a. 
     a. If mm=m1 set mm=1, kr=kk(mm), kk(mm)=lm, go to step 1; otherwise 
        set mm=mm+1, kr=kk(mm), kk(mm)=lm, go to step 2. 
     b. Set: iopt=iopt+1. 
        If iopt≠m1 go to step 5.a, otherwise go to step 6. 
Step 6) Compute the solution. 
        Set: b(m)=ys(lm). 
        Do loop for i=1,m1: ysol(i)=y(kk(i)); do loop for j=1,m1: xsol(i,j)=x(kk(i),j), end do; end do. 
        Set: jj=1. 
     a. Set: ysk=ysol(jj). 
        Do loop for j=jj,m1: xsx(j)=xsol(jj,j), end do. 
        Do loop for i=jj,m1: if i=jj go to continue; otherwise 
        ysol(i)=ysol(i)-ysk, do loop for j=jj,m1: xsol(i,j)=xsol(i,j)-xsk(j), end do; 
        set ysol(i)=ysol(i)/xsol(i,jj), continue, end do. 
     b. Do loop for i=jj,m1: if i=jj go to continue, otherwise, do loop for 
        j=jj+1,m1: xsol(i,j)=xsol(i,j)/xsol(i,jj), end do; continue; end do. 
     c. If jj=m2 go to step 6.d; otherwise go to step 6.a. 
     d. Do loop for i=1,m2: k=m-i, s=ysol(k); do loop for j=k,m1, 
        s=s-b(j+1)*xsol(k,j) end do, b(k)=s, end do. 
        Set: s=y(kk(1)). 
        Do loop for j=1,m1: s=s-b(j+1)*x(kk(1),j), b(1)=s, end do. 
        Print: ((b(j),j=1,m),(kk(j),j=1,m1),lm,iter). 
        Stop. 
END 
 

The major portion of computation in this program is the transformation of two-dimensional arrays. Passing columns 
of these arrays to other subroutines, which involve only one-dimensional arrays saves the time of computation (see, Barrodale 
and Roberts (1974)). Subroutine COL1, COL2, and COL3 have been coded to do this task for subtraction, multiplication and 
division, and for only division respectively. Function LWMED, which is used to compute the weighted median, has been 
introduced in section 2.1. 
SUBROUTINE COL1(v1,v2) 
Step 0) Initialization 
        Real: v2(1). 
        Common /c1/i1,i2. 
Step 1) Subtraction. 



Copyright © CC-BY-NC 2019, BJMSR 

 

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

17 
               

        Do loop for i=i1,i2: v2(i)=v2(i)-v1, end do. 
        Return. 
END 
 
SUBROUTINE COL2(ysk,jj,v1,ys,v2) 
Step 0) Initialization. 
        Real: v1(1),v2(1),ys(1). 
        Common /c1/i1,i2. 
Step 1) Compute weights and ratios. 
        If jj≠1 go to step 1.a,; otherwise do loop for i=i1,i2: 
        v1(i)=v2(i), ys(i)=(ys(i)-ysk)/v2(i). 
        Return. 
     a. Do loop for i=i1,i2: v1(i)=v1(i)*v2(i), ys(i)=(ys(i)-ysk)/v2(i), 
        end do. 
        Return. 
END 
 
SUBROUTINE COL3(v1,v2) 
Step 0) Initialization. 
        Real: v1(1),v2(1),ys(1). 
        Common /c1/i1,i2. 
Step 1) Division. 
        Do loop for i=i1,i2: v1(i)=v1(i)/v2(i), end do. 
        Return. 
END 
 
4.4 Initial value problem 

Selection of the arbitrary points k1,k2,...,k(m-1) at the first stage of computation has an important role in the number 
of iterations necessary to reach the optimal solution. One way to select these m-1 points may be based on applying another 
estimator to guess those points which their residuals are smallest in absolute values among all estimated residuals.  That is  one 
estimator for example least squares is applied on the data, and smallest m-1 residuals in absolute values are selected, and the 
corresponding points to these residuals are nominated for starting algorithms 1, 2, 3 or 4. However, we examined this procedure 
for algorithms 2 and 4 and found that it makes them more efficient. 

In the process of this study, we found some points which might be leading to develop the field of L1 norm 
computation. In the proposed algorithm 4, we could not find a criterion to delete that observation from the current set of 
observations on the basis, which reduces the objective function value more than other points. If anyone finds a criterion like the 
heuristic method of Bloomfield and Steiger (1980), the number of iterations will decrease very much and will reduce 
computation time. 

The main obstacle to applying discrete differentiation technique to operational problems is reordering of the residuals 
in descending order before starting the computation. If this obstacle is removed, a similar algorithm to that of simple restricted 
linear regression for multiple regression may be proposed. 
 
References 
I. Barrodale, F.D.K. Roberts (1974) Algorithm 478: Solution of an overdetermined system of equations in the L 1 norm. 

Comm. ACM, 17, 319-320. 
C.M. Bender, S.A. Orszag (1978) Advanced mathematical methods for scientists and engineers, McGraw-Hill. 
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 

http://www.bidabad.com/doc/lae-I.pdf
http://www.bidabad.com/doc/lae-II.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 

18 
               

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.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. Stat. Comput. 1, 290-301. 
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. 

J. Chambers (1971) Algorithm 410: partial sorting. Comm. ACM, 14, 357-358. 
F.H. Clarke (1983) Optimization and nonsmooth analysis. Wiley, New York C.A.R. Hoare (1961) Algorithm 63 partition; 

64, quicksort; and 65, find., Comm. ACM, 4, July, 321-322. 
C.A.R. Hoare (1962) Quicksort. Comput. J., 5, 10-15. 
O.J. Karst (1958) Linear curve fitting using least deviations. JASA, 53, 118-132. 
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). 
S. Scowen (1965) Algorithm 271 Quickersort. Comm. ACM, 8, 669-670. 
R.S. Singleton (1969) Algorithm 347 Sort. Comm. ACM, 12, 185-186. 
L.D. Taylor (1974) Estimating by minimizing the sum of absolute errors. In 
P. Zarembka (ed.) Frontiers in econometrics. Academic Press. 
 
 
 
 
 
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/). 
 

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
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

