# Novel Global Harmony Search Algorithm for Least Absolute Deviation

Authors: Longquan Yong
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

The method of least absolute deviation (LAD) finds applications in many areas, due to its robustness compared to the least squares regression (LSR) method. LAD is robust in that it is resistant to outliers in the data. This may be helpful in studies where outliers may be ignored. Since LAD is nonsmooth optimization problem, this paper proposed a metaheuristics algorithm named novel global harmony search (NGHS) for solving. Numerical results show that the NGHS method has good convergence property and effective in solving LAD.


---

## Body

## 1. Introduction

Least squares regression (LSR) method is one of the oldest and the most widely used statistical tools for linear models. Its theoretical properties have been extensively studied and are fully understood. Despite its many superior properties, the LSR estimate can be sensitive to outliers and, therefore, nonrobust [1]. In order to overcome these problems, the researchers have recently investigated the alternative regression method which is least absolute deviation (LAD) method.Least absolute deviation (LAD), also known as least absolute errors (LAE), least absolute value (LAV), or least absolute residual (LAR), or theL
1 norm problem, is a mathematical optimization technique similar to the LSR that attempts to find a function which closely approximates a set of data (
x
i
,
y
i
)
,
i
=
1,2
,
…
,
n. In the simple case, the approximation function is a simple “trend line” in two-dimensional Cartesian coordinates [2]. The method minimizes the sum of absolute errors (SAE) (the sum of the absolute values of the vertical “residuals” between points generated by the function and corresponding points in the data).The method of LAD finds applications in many areas, due to its robustness compared to the LSR method. LAD is robust in that it is resistant to outliers in the data. This may be helpful in studies where outliers may be safely and effectively ignored. Despite its long history and many ground-breaking works, the LAD has not been explored in theory as well as in application to the extent as the LSR [3]. This is largely because computing LAD estimates are more difficult than computing LSR estimates due to the fact that algorithmic method must be employed to calculate LAD estimates. Therefore, LAD estimation method is not popular as LSR method.Over the past few years, a number of approaches have been developed for solving the LAD using classical mathematical programming methods. Since Charnes et al. [4] reduced the LAD method to a linear programming problem, the computational difficulty is now entirely overcome by the availability of computing power and the effectiveness of linear programming. A comprehensive summary of the subject can be found in Koenker and D’Orey [5]. Large sample properties of the LAD estimates are obtained in Koenker and Bassett Jr. [6], Chen et al. [7], and Pollard [8]. Due to these developments in theoretical and computational aspects, the LAD method has become increasingly popular. In particular, it has many applications in econometrics and biomedical studies; see Bassett and Koenker [9], Powell [10], Buchinsky [11], among many others.Since LAD is a class of nondifferentiable optimization problem, this paper proposed a metaheuristics algorithm named novel global harmony search (NGHS) [12] for solving LAD model. The NGHS algorithm includes two important operations: position updating and genetic mutation with a small probability (p
m). The former enables the worst harmony of harmony memory to move to the global best harmony rapidly in each iteration, and the latter can effectively maintain population diversity and prevent the NGHS from trapping into the local optimum.The remaining sections of this chapter are organized as follows. In Section2, LAD model is shown. In Sections 3 and 4, the classical HS and NGHS are described. Numerical results are presented and compared in Section 5. Lastly, Section 6 outlines our conclusions.


## 2. Least Absolute Deviations Model

The LAD method is a widely known alternative to the classical LSR method for statistical analysis of linear regression models. Instead of minimizing the sum of squared errors (SSE) in LSR, it minimizes the sum of absolute errors (SAE).Suppose that the data set consists ofn points (data pairs) (
x
i
,
y
i
)
,
i
=
1,2
,
…
,
n, where x
i is an independent variable and y
i is a dependent variable whose value is found by observation. We want to find a function f such that f
(
x
i
)
=
y
i.To attain this goal, we suppose that the functionf is of a particular form containing some parameters which need to be determined. For instance, the simplest form would be linear: f
(
x
,
β
)
=
β
0
+
β
1
x, where β
0 and β
1  are parameters whose values are not known but which we would like to estimate.We now seek estimated values of the unknown parameters that minimize the SAE:(1)
min
⁡
S
=
∑
i
=
1
n
|
f
(
x
i
)
-
y
i
|
.Let(2)
A
T
=
[
1
1
⋯
1
x
1
x
2
⋯
x
n
]
,
β
=
[
β
0
β
1
]
,
b
T
=
[
y
1
,
y
2
,
…
,
y
n
]
.
Then LAD problem is equivalent to following L
1 optimization problem:
(3)
min
⁡
S
=
∑
i
=
1
n
|
f
(
x
i
)
-
y
i
|
=
∥
A
β
-
b
∥
1
.
Meanwhile, then LSR problem is equivalent to following L
2 optimization problem
(4)
min
⁡
S
=
∑
i
=
1
n
|
f
(
x
i
)
-
y
i
|
2
=
∥
A
β
-
b
∥
2
2
.
For a given data set, the method of LAD may produce multiple solutions, whereas the method of LSR always produces only one solution (the regression line is unique).


## 3. Harmony Search Algorithm

Recently, a new class of metaheuristics, named harmony search (HS), has been developed. The HS algorithm proposed in [13] has been developed in an analogy with music improvisation process where musicians in an ensemble continue to polish their pitches in order to obtain better harmony. Jazz improvisation seeks to find musically pleasing harmony similar to the optimum design process which seeks to find optimum solution. The pitch of each musical instrument determines the aesthetic quality, just as the objective function value is determined by the set of values assigned to each decision variable. In addition, HS uses a stochastic random search instead of a gradient search.The steps in the procedure of standard harmony search algorithm (HS) are as follows.Step 1 (initialize the problem and algorithm parameters).
The optimization problem is specified as follows:(5)
Minimize
f
(
x
)
subject
to
x
i
∈
X
i
,
i
=
1,2
,
…
,
N
,
where f
(
x
) is an objective function; x is the set of each decision variable x
i; N is the number of decision variables; X
i is the set of the possible range of values for each decision variable; X
i
:
x
i
L
≤
x
i
≤
x
i
U. The HS algorithm parameters are also specified in this step. These are the harmony memory size (HMS) or the number of solution vectors in the harmony memory; harmony memory considering rate (HMCR); pitch adjusting rate (PAR); and the number of improvisations (Tmax) or stopping criterion.Step 2 (initialize the harmony memory).
The HM matrix is filled with as many randomly generated solution vectors as the HMS(6)
HM
=
[
x
1
f
(
x
1
)
x
2
f
(
x
2
)
⋮
⋮
x
HMS
f
(
x
HMS
)
]
=
[
x
1
1
x
2
1
⋯
x
N
1
f
(
x
1
)
x
1
2
x
2
2
⋯
x
N
2
f
(
x
2
)
⋮
⋮
⋯
⋮
⋮
x
1
HMS
x
2
HMS
⋯
x
N
HMS
f
(
x
HMS
)
]
.Step 3 (improvise a new harmony).
Generate a new harmony is called “improvisation.” A new harmony vector,x
′
=
(
x
1
′
,
x
2
′
,
…
,
x
N
′
), is generated based on three rules: (1) memory consideration, (2) pitch adjustment, and (3) random selection. The procedure works as Algorthim 1.
x
i
′
(
i
=
1,2
,
…
,
N
) is the ith component of x
′, and x
i
j
(
j
=
1,2
,
…
,
HMS
) is the ith component of the jth candidate solution vector in HM. Here r
(
-
1,1
) is uniform distribution random number in the region of (
-
1,1
), b
w is an arbitrary distance bandwidth, and rand
(
) is uniform distribution random number in the region of (
0,1
).Algorithm 1:Generating a new harmony by classical HS algorithm.
For each  i
∈
[
1,2
,
…
,
N
]  do
If  rand() < HMCR
x
i
′
=
x
i
j
(
j
=
1,2
,
…
,
HMS
)  //memory consideration
If rand() < PAR
x
i
′
=
x
i
′
+
r
(
-
1,1
)
×
b
w    // pitch adjustment
x
i
′
=
min
⁡
[
max
⁡
⁡
(
x
i
′
,
x
i
L
)
,
x
i
U
]  //truncation processing
End
Else
x
i
′
=
x
i
L
+
rand
(
)
×
(
x
i
U
-
x
i
L
) //random selection
End
EndStep 4 (update harmony memory).
If the new harmony vectorx
′
=
(
x
1
′
,
x
2
′
,
…
,
x
N
′
) is better than the worst harmony in the HM, judged in terms of the objective function value, the new harmony is included in the HM and the existing worst harmony is excluded from the HM.Step 5 (check stopping criterion).
If the stopping criterion (maximum number of improvisations) is satisfied, computation is terminated. Otherwise, Steps3 and 4 are repeated.


## 4. Novel Global Harmony Search Algorithm

Experiments with the classical HS algorithm over the benchmark problems show that the algorithm suffers from the problem of premature and/or false convergence and slow convergence especially over multimodal fitness landscape. To enrich the searching behavior and to avoid being trapped into local optimum, more improved HS algorithms were presented. Zou et al. [12] proposed an NGHS algorithm and it has been successfully used in reliability problems [14], and 0-1 knapsack problems [15].The NGHS algorithm modifies the improvisation step of the HS such that the new harmony can mimic the global-best harmony in the HM. In Step3 it works as Algorthim 2.Algorithm 2:Generating a new harmony by NGHS algorithm.
For each  i
∈
[
1,2
,
…
,
N
]  do
x
s
=
2
×
x
i
best
-
x
i
worst
x
s
=
min
⁡
[
max
⁡
(
x
s
,
x
i
L
)
,
x
i
U
]     //truncation processing
x
i
′
=
x
i
worst
+
rand
(
)
×
(
x
s
-
x
i
worst
)   //position updating
If  rand
(
)
<
p
m
x
i
′
=
x
i
L
+
rand
(
)
×
(
x
i
U
-
x
i
L
) //genetic mutation
End
EndHere, ‘‘best’’ and ‘‘worst’’ are the indexes of the best and worst harmony in HM, respectively.The NGHS proposed algorithm has strong global search ability in the early stage of optimization and has strong local search ability in the late stage of optimization. In the early stage of optimization, all solution vectors are sporadic in feasible space, so mostx
s is large, which is beneficial to the global search of the NGHS, while in the late stage of optimization, all nonbest solution vectors are inclined to move to the global-best solution vector, so most solution vectors are close to each other. In this case, most x
s is small and most trust regions are narrow, which is beneficial to the local search of the NGHS.In Step4, the NGHS replaces the worst harmony x
worst in HM with the new harmony x
′ even if x
′ is worse than x
worst.The NGHS and the HS are different in following. (i) In Step1, harmony memory considering rate (HMCR), pitch adjusting rate (PAR), and adjusting step (b
w) are excluded from the NGHS, and genetic mutation probability (p
m) is included in the NGHS. (ii) The HS carries out mutation with the probability HMCR × PAR and carries out random selection with the probability 1-HMCR, while the NGHS carries out genetic mutation with the probability p
m. In fact, both operations are exactly the same, and they are used to keep the individual variety better, which can effectively prevent both algorithms from being trapped into the local optimum.


## 5. Computational Results

In this section we solve an LAD problem in order to illustrate the implementation and efficiency of the NGHS method. All the experiments were performed on MatlabR2009a system with Intel(R) Core(TM)4
×
3.3 GHz and 2 GB RAM.The data pairs(
x
i
,
y
i
)
,
i
=
1,2
,
…
,
23, are listed in Tables 1 and 2, where x
i is an independent variable and y
i is a dependent variable whose value is found by observation.Table 1
Given data.
i
1
2
3
4
5
6
7
8
9
10
11
12
x
i
20
19.6
19.6
19.4
18.4
19
19
18.3
18.2
18.6
19.2
18.2
y
i
1.0
1.2
1.1
1.4
2.3
1.7
1.7
2.4
2.1
2.1
1.2
2.3
i
13
14
15
16
17
18
19
20
21
22
23
x
i
18.7
18.5
18
17.4
16.5
17.2
17.3
17.8
17.3
18.4
16.9
y
i
1.9
2.4
2.6
2.9
4.0
3.3
3.0
3.4
2.9
1.9
3.9Table 2
Given data with outlier.
i
1
2
3
4
5
6
7
8
9
10
11
12
x
i
20
19.6
19.6
19.4
18.4
19
19
18.3
18.2
18.6
19.2
18.2
y
i
1.0
1.2
1.1
1.4
2.3
1.7
1.7
2.4
2.1
2.1
1.2
2.3
i
13
14
15
16
17
18
19
20
21
22
23
x
i
18.7
18.5
10
17.4
16.5
17.2
17.3
17.8
17.3
18.4
16.9
y
i
1.9
2.4
2.6
2.9
4.0
3.3
3.0
3.4
2.9
1.9
3.9For data pairs listed in Tables1 and 2, solved by LSR model [16], respectively, the results of data fitting are shown in Figure 1.Figure 1
Data fitting results by LSR method. Heref
1
(
x
)
=
-
0.8940
x
+
18.6746 is fitting results of data in Table 1, and f
2
(
x
)
=
-
0.2288
x
+
6.4055 is fitting results of data in Table 2.From Figure1, LSR method is sensitive to outliers and, therefore, nonrobust. Its performance in terms of accuracy and statistical inferences may be compromised when the errors are large and heterogeneous.Following we use LAD model to fit given data with outlier in Table2. Since LAD is a kind of nondifferentiable optimization problem, we solved it by NGHS. In order to show the optimization capabilities of the NGHS, the results of classical HS, HSCH [17], and HSWB [18] are listed too. To make the comparison fair, the populations for all the competitor algorithms were initialized using the same random seeds. The HS-variants algorithm parameters were set the same parameters: harmony memory size HMS = 15, harmony memory consideration rate HMCR = 0.85, pitch adjusting rate PAR = 0.35, and Tmax = 400. In NGHS, we set p
m
=
0.005. If using a random number, we set the random-number generator to the state of 0 so that the same data can be regenerated.To judge the accuracy of different algorithms, 10 independent runs of each of the four algorithms were carried out and the best, the mean, and the worst fitness values and the standard deviation (Std) were recorded in Table3.Table 3
The Statistical results for 10 runs on given data in Table2.
Algorithm
Best
Mean
Worst
Standard
Mean_Time (s)
LAD
HS
1.0830
E
+
01
2.2618
E
+
01
5.9845
E
+
01
1.7983
E
+
01
2.9411
E
-
03
HSCH
1.0723
E
+
01
2.9073
E
+
01
1.4305
E
+
02
4.0622
E
+
01
3.4392
E
-
03
HSWB
1.1861
E
+
01
1.9807
E
+
02
7.5315
E
+
02
2.3679
E
+
02
3.5885
E
-
03
NGHS
9.9041
E
+
00
1.2058
E
+
01
1.4555
E
+
01
1.7163
E
+
00
3.5500
E
-
03Figure2 shows the convergence and its boxplot figure of the best fitness in the population for the different algorithms. The values plotted for every generation are averaged over 10 independent runs. As can be seen, the NGHS algorithm is the best.The convergence and its boxplot of the best fitness for LAD.
(a)
(b)Moreover, the results of data fitting by NGHS algorithm are shown in Figure3.Figure 3
Data fitting results by LSR and LAD. Heref
1
(
x
)
=
-
0.8940
x
+
18.6746 is fitting results of data in Table 1 by LSR, f
2
(
x
)
=
-
0.2288
x
+
6.4055 is fitting results of data in Table 2 by LSR, and f
3
(
x
)
=
-
0.7306
x
+
15.5965 is fitting results of data in Table 2 by LAD.Unlike the LSR method, the LAD method is not sensitive to outliers and produces robust estimates.


## 6. Conclusion

The LAD method is a widely known alternative to the classical LSR method for statistical analysis of linear regression models. Since LAD is a class of nondifferentiable optimization problem, this paper given NGHS algorithm for solving. Unlike the LSR, the LAD is not sensitive to outliers and produces robust estimates. Moreover, NGHS algorithm is effective for solving LAD model. The future works include investigating the application of LAD or other algorithm for solving LAD.



---

*Source: 632975-2014-07-21.xml*