Accumulative Approach in Multistep Diagonal Gradient-Type Method for Large-Scale Unconstrained Optimization

Abstract

This paper focuses on developing diagonal gradient-type methods that employ accumulative approach in multistep diagonal updating to determine a better Hessian approximation in each step. The interpolating curve is used to derive a generalization of the weak secant equation, which will carry the information of the local Hessian. The new parameterization of the interpolating curve in variable space is obtained by utilizing accumulative approach via a norm weighting defined by two positive definite weighting matrices. We also note that the storage needed for all computation of the proposed method is just 𝑂(𝑛). Numerical results show that the proposed algorithm is efficient and superior by comparison with some other gradient-type methods.

1. Introduction

Consider the unconstrained optimization problem: min𝑓(𝑥),𝑥∈𝑅^𝑛,(1.1) where 𝑓∶𝑅^𝑛→𝑅 is twice continuously differentiable function. The gradient-type methods for solving (1.1) can be written as 𝑥_𝑘+1=𝑥_𝑘−𝐵^𝑘−1𝑔_𝑘,(1.2) where 𝑔_𝑘 and 𝐵_𝑘 denote the gradient and the Hessian approximation of 𝑓 at 𝑥_𝑘, respectively. By considering 𝐵_𝑘=𝛼_𝑘𝐼, Barzilai and Borwein (BB) [1] give 𝛼_𝑘+1=𝑠^𝑇𝑘𝑦_𝑘𝑠^𝑇𝑘𝑠_𝑘,(1.3) where it is derived by minimizing ‖𝛼_𝑘+1𝑠_𝑘−𝑦_𝑘‖₂ respect to 𝛼 with 𝑠_𝑘=𝑥_𝑘+1−𝑥_𝑘 and 𝑦_𝑘=𝑔_𝑘+1−𝑔_𝑘. Recently, some improved one-step gradient-type methods [2–5] in the frame of BB algorithm were proposed to solve (1.1). It is proposed to let 𝐵_𝑘 be a diagonal nonsingular approximation to the Hessian and a new approximating matrix 𝐵_𝑘+1 to the Hessian is developed based on weak secant equation of Dennis and Wolkowicz [6] 𝑠^𝑇𝑘𝐵_𝑘+1𝑠_𝑘=𝑠^𝑇𝑘𝑦_𝑘.(1.4) In one-step method, data from one previous step is used to revise the current approximation of Hessian. Later Farid and Leong [7, 8] proposed multistep diagonal gradient methods inspired by the multistep quasi-Newton method of Ford [9, 10]. In this multistep framework, a fixed-point approach for interpolating polynomials was derived from data in previous iterations (not only one previous step) [7–10]. General approach of multistep method is based on the measurement of distances in the variable space where the distance of every iterate is measured from one-selected iterate. In this paper, we are interested to develop multistep diagonal updating based on accumulative approach for defining new parameter value of interpolating curve. From this point, the distance is accumulated between consecutive iterates as they are traversed in the natural sequence. For measuring the distance, we need to parameterize the interpolating polynomial through a norm that is defined by a positive definite weighting matrix, say 𝑀. Therefore, the performance of the multistep method may be significantly improved by carefully defining the weighting matrix. The rest of paper is organized as follows. In Section 2, we discuss a new multistep diagonal updating scheme based on the accumulative approach. In Section 3, we establish the global convergence of our proposed method. Section 4 presents numerical result and comparisons with BB method and one-step diagonal gradient method are reported. Conclusions are given in Section 5.

2. Derivation of the New Diagonal Updating via Accumulative Approach

This section motivates to state new implicit updates for diagonal gradient-type method through accumulative approach to determining a better Hessian approximation at each iteration. In multistep diagonal updating methods, weak secant equation (1.4) may be generalized by means of interpolating polynomials, instead of employing data just from one previous iteration like in one-step methods. Our aim is to derive efficient strategies for choosing a suitable set of parameters to construct the interpolating curve and investigate the best norm for measurement of the distances required to parameterize the interpolating polynomials. In general, this method obeys the recursive formula of the form 𝑥_𝑘+1=𝑥_𝑘−𝛼_𝑘𝐵^𝑘−1𝑔_𝑘,(2.1) where 𝑥_𝑘 is the 𝑘th iteration point, 𝛼_𝑘 is step length which is determined by a line search, 𝐵_𝑘 is an approximation to the Hessian in a diagonal form, and 𝑔_𝑘 is the gradient of 𝑓 at 𝑥_𝑘. Consider a differentiable curve 𝑥(𝜏) in 𝑅^𝑛. The derivative of 𝑔(𝑥(𝜏)), at point 𝑥(𝜏^∗), can be obtained by applying the chain rule: 𝑑𝑔𝑑𝜏|||_𝜏=𝜏^∗=𝐺(𝑥(𝜏))𝑑𝑥𝑑𝜏|||_𝜏=𝜏^∗.(2.2) We are interested to derive a relation that will be satisfied by the approximation of Hessian in diagonal form at 𝑥_𝑘+1. If we assume that 𝑥(𝜏) passes through 𝑥_𝑘+1 and choose 𝜏^∗ so that 𝑥𝜏^∗=𝑥_𝑘+1,(2.3) then we have 𝐺𝑥_𝑘+1𝑑𝑥^𝑘+1𝑑𝜏=𝑑𝑔𝑥^𝑘+1𝑑𝜏.(2.4) As in this paper, we use two-step method, therefore; we use information of most recent points 𝑥_𝑘−1, 𝑥_𝑘, 𝑥_𝑘+1 and their associated gradients. Consider 𝑥(𝜏) as the interpolating vector polynomial of degree 2: 𝑥𝜏_𝑗=𝑥_{𝑘+𝑗−1}𝑗=0,1,2.(2.5) The selection of distinct scalar value 𝜏_𝑗 efficiently through the new approach is the main contribution of this paper and will be discussed later in this section. Let ℎ(𝜏) be the interpolation for approximating the gradient vector: ℎ𝜏_𝑗=𝑔𝑥_{𝑘+𝑗−1}𝑗=0,1,2.(2.6) By denoting 𝑥(𝜏₂)=𝑥_𝑘+1 and defining 𝑑𝑥𝜏₂𝑑𝜏=𝑟^𝑘𝑑𝑔𝑥𝜏₂,𝑑𝜏=𝑤^𝑘,(2.7) we can obtain our desired relation that will be satisfied by the Hessian approximation at 𝑥_𝑘+1 in diagonal form. Corresponding to this two-step approach, weak secant equation will be generalized as follows: 𝑟^𝑇𝑘𝐵_𝑘+1𝑟_𝑘=𝑟^𝑇𝑘𝑤_𝑘.(2.8) Then, 𝐵_𝑘+1 can be obtained by using an appropriately modified version of diagonal updating formula in [3] as follows: 𝐵_𝑘+1=𝐵_𝑘+𝑟^𝑇𝑘𝑤_𝑘−𝑟^𝑇𝑘𝐵_𝑘𝑟_𝑘tr𝐹^2𝑘𝐹^𝑘,(2.9) where 𝐹_𝑘=diag((𝑟^𝑘(1))²,(𝑟^𝑘(2))²,…,(𝑟^𝑘(𝑛))²). Now, we attempt to construct an algorithm for finding desired vector 𝑟_𝑘 and 𝑤_𝑘 to improve the Hessian approximation. The proposed method is outlined as follows. First, we seek to derive strategies for choosing a suitable set of values 𝜏₀, 𝜏₁, and 𝜏₂. The choice of {𝜏_𝑗}^2𝑗=0 is such that to reflect distances between iterates 𝑥_𝑘 in 𝑅^𝑛 that are dependent on some metric of the following general form: 𝜙_𝑀𝑧₁,𝑧₂=(𝑧₁−𝑧₂)^𝑇𝑀(𝑧₁−𝑧₂)^1/2.(2.10) The establishment on 𝜏_𝑗 can be made via the so-called accumulative approach where the accumulating distances (measured by the metric 𝜙_𝑀) between consecutive iterates are used to approximate 𝜏_𝑗. This leads to the following definitions (where without loss of generality, we take 𝜏₁ to be origin for value of 𝜏): 𝜏𝜏₁=0,_𝑗=𝜏_𝑗+1−𝜙_𝑀𝑥_𝑘+𝑗,𝑥_{𝑘+𝑗−1}𝑗=0,2.(2.11) Then, we can construct the set {𝜏_𝑗}^2𝑗=0 as follows: 𝜏₀=𝜏₁−𝜙_𝑀𝑥_𝑘,𝑥_𝑘−1_def=−‖‖𝑥_𝑘−𝑥_𝑘−1‖‖_𝑀𝜏₂=−‖‖𝑠_𝑘−1‖‖_𝑀,=𝜙_𝑀𝑥_𝑘+1,𝑥_𝑘_def=‖‖𝑥_𝑘+1−𝑥_𝑘‖‖_𝑀=‖‖𝑠_𝑘‖‖_𝑀,(2.12) where 𝑟_𝑘 and 𝑤_𝑘 are depending on the value of 𝜏. As the set {𝜏_𝑗}^2𝑗=0 measures the distances, therefore they need to be parameterized the interpolating polynomials via a norm defined by a positive definite matrix 𝑀. It is necessary to choose 𝑀 with some care, while improving the approximation of Hessian can be strongly influenced via the choice of 𝑀. Two choices for the weighting matrix 𝑀 are considered in this paper. In first choice, if 𝑀=𝐼, the ‖⋅‖_𝑀 reduces to the Euclidean norm, and then we obtain the following 𝜏_𝑗 values accordingly: 𝜏₂=‖‖𝑠_𝑘‖‖₂,𝜏₁=0,𝜏₀=−‖‖𝑠_𝑘−1‖‖₂.(2.13) The second choice of weighting matrix 𝑀 is to take 𝑀=𝐵_𝑘, where the current 𝐵_𝑘 is diagonal approximation to the Hessian. By these two means, the measurement of the relevant distances is determined by the properties of the current quadratic approximation (based on 𝐵_𝑘) to the objective function: 𝜏₂=𝑠_𝑘𝐵_𝑘s_𝑘^1/2,𝜏₁=0,𝜏₀=−𝑠_𝑘−1𝐵_𝑘𝑠_𝑘−1^1/2.(2.14) Since 𝐵_𝑘 is a diagonal matrix, then it is not expensive to compute {𝜏_𝑗}^2𝑗=0 at each iteration. The quantity 𝛿 is introduced here and defined as follows: 𝛿=^def𝜏²−𝜏₁𝜏₁−𝜏₀,(2.15) and 𝑟_𝑘 and 𝑤_𝑘 are given by the following expressions: 𝑟_𝑘=𝑠_𝑘−𝛿²1+2𝛿𝑠^𝑘−1𝑤,(2.16)_𝑘=𝑦_𝑘−𝛿²1+2𝛿𝑦^𝑘−1.(2.17) To safeguard on the possibility of having very small or very large 𝑟^𝑇𝑘𝑤_𝑘, we require that the condition 𝜀₁‖‖𝑟_𝑘‖‖²²≤𝑟^𝑇𝑘𝑤_𝑘≤𝜀₂‖‖𝑟_𝑘‖‖²²(2.18) is satisfied (we use 𝜀₁=10⁻⁶ and 𝜀₂=10⁶). If not, then we replace 𝑟_𝑘=𝑠_𝑘 and 𝑤_𝑘=𝑦_𝑘. More that the Hessian approximation (𝐵_𝑘+1) might not preserve the positive definiteness in each step. One of the fundamental concepts in this paper is to determine an “improved” version of the Hessian approximation to be used even in computing the metric when 𝑀=𝐵_𝑘 and a weighing matrix as norm should be positive definite. To ensure that the updates remain positive definite, a scaling strategy proposed in [7] is applied. Hence, the new updating formula that incorporates the scaling strategy is given by 𝐵_𝑘+1=𝜂_𝑘𝐵_𝑘+𝑟^𝑇𝑘𝑤_𝑘−𝜂_𝑘𝑟^𝑇𝑘𝐵_𝑘𝑟_𝑘tr𝐹^2𝑘𝐹^𝑘,(2.19) where 𝜂_𝑘=min𝑟^𝑇𝑘𝑤_𝑘𝑟^𝑇𝑘𝐵_𝑘𝑟_𝑘,1.(2.20) This guarantees that the updated Hessian approximation is positive. Finally, the new accumulative MD algorithm is outlined as follows.

2.1. Accumulative MD Algorithm

Step 1. Choose an initial point 𝑥₀∈𝑅^𝑛, and a positive definite matrix 𝐵₀=𝐼.
Let 𝑘∶=0.

Step 2. Compute 𝑔_𝑘. If ‖𝑔_𝑘‖≤𝜖, stop.

Step 3. If 𝑘=0, set 𝑥₁=𝑥₀−(𝑔₀/‖𝑔₀‖). If 𝑘=1 set 𝑟_𝑘=𝑠_𝑘 and 𝑤_𝑘=𝑦_𝑘 go to Step 5.

Step 4. If 𝑘≥2 and 𝑀=𝐼 is considered, compute {𝜏_𝑗}^2𝑗=0 from (2.13).
Else if 𝑀=𝐵_𝑘, compute {𝜏_𝑗}^2𝑗=0 from (2.14).
Compute 𝛿_𝑘, 𝑟_𝑘,𝑤_𝑘 and 𝜂_𝑘, from (2.15), (2.16),  (2.17), and (2.20), respectively.
If 𝑟^𝑇𝑘𝑤_𝑘≤10⁻⁴‖𝑟_𝑘‖₂‖𝑤_𝑘‖₂, set 𝑟_𝑘=𝑠_𝑘 and 𝑤_k=𝑦_𝑘.

Step 5. Compute 𝑑_𝑘=−𝐵^𝑘−1𝑔_𝑘 and calculate 𝛼_𝑘>0 such that the Armijo [11], condition holds:
𝑓(𝑥_𝑘+1)≤𝑓(𝑥_𝑘)+𝜎𝛼_𝑘𝑔^𝑇𝑘𝑑_𝑘, where 𝜎∈(0,1) is a given constant.

Step 6. Let 𝑥_𝑘+1=𝑥_𝑘−𝛼_𝑘𝐵^𝑘−1𝑔_𝑘, and update 𝐵_𝑘+1 by (2.19).

Step 7. Set 𝑘∶=𝑘+1, and return to Step 2.

3. Convergence Analysis

This section is devoted to study the convergence of accumulative MD algorithm, when applied to the minimization of a convex function. To begin, we give the following result, which is due to Byrd and Nocedal [12] for the step generated by the Armijo line search algorithm. Here and elsewhere, ‖⋅‖ denotes the Euclidean norm.

Theorem 3.1. Assume that 𝑓 is a strictly convex function. Suppose the Armijo line search algorithm is employed in a way that for any 𝑑_𝑘 with 𝑑^𝑇𝑘𝑔_𝑘<0, the step length, 𝛼_𝑘 satisfies 𝑓𝑥_𝑘+𝛼_𝑘𝑑_𝑘≤𝑓𝑥_𝑘+𝜎𝛼_𝑘𝑔^𝑇𝑘𝑑_𝑘,(3.1) where 𝛼_𝑘∈[𝜏,𝜏^], 0<𝜏<𝜏^ and 𝜎∈(0,1). Then, there exist positive constants 𝜌₁ and 𝜌₂ such that either 𝑓𝑥_𝑘+𝛼_𝑘𝑑_𝑘−𝑓𝑥_𝑘≤−𝜌₁𝑑^𝑇𝑘𝑔_𝑘²‖‖𝑑_𝑘‖‖²(3.2) or 𝑓𝑥_𝑘+𝛼_𝑘𝑑_𝑘−𝑓𝑥_𝑘≤−𝜌₂𝑑^𝑇𝑘𝑔_𝑘(3.3) is satisfied.

We can apply Theorem 3.1 to establish the convergence of some Armijo-type line search methods.

Theorem 3.2. Assume that 𝑓 is a strictly convex function. Suppose that the Armijo line search algorithm in Theorem 3.1 is employed with 𝑑_𝑘 chosen to obey the following conditions: there exist positive constants 𝑐₁ and 𝑐₂ such that −𝑔^𝑇𝑘𝑑_𝑘≥𝑐₁‖‖𝑔_𝑘‖‖²,‖‖𝑑_𝑘‖‖≤𝑐₂‖‖𝑔_𝑘‖‖,(3.4) for all sufficiently large 𝑘. Then, the iterates 𝑥_𝑘 generated by the line search algorithm have the property that liminf_𝑘→∞‖‖𝑔_𝑘‖‖=0.(3.5)

Proof. By (3.4), we have that either (3.2) or (3.6) becomes 𝑓𝑥_𝑘+𝛼_𝑘𝑑_𝑘−𝑓𝑥_𝑘≤−𝑐‖‖𝑔_𝑘‖‖²,(3.6) for some positive constants. Since 𝑓 is strictly convex, it is also bounded below. Then, (3.1) implies that 𝑓(𝑥_𝑘+𝛼_𝑘𝑑_𝑘)−𝑓(𝑥_𝑘)→0 as 𝑘→∞. This also implies that ‖𝑔_𝑘‖→0 as 𝑘→∞ or at least liminf_𝑘→∞‖‖𝑔_𝑘‖‖=0.(3.7)

To prove that the accumulative MD algorithm is globally convergent when applied to the minimization of a convex function, it is sufficient to show that the sequence {‖𝐵_𝑘‖} generated by (2.19)-(2.20) is bounded both above and below, for all finite 𝑘 so that its associated search direction satisfies condition (3.4). Since 𝐵_𝑘 is diagonal, it is enough to show that each element of 𝐵_𝑘 says 𝐵^𝑘(𝑖);  𝑖=1,…,𝑛 is bounded above and below by some positive constants. The following theorem gives the boundedness of {‖𝐵_𝑘‖}.

Theorem 3.3. Assume that 𝑓 is strictly convex function where there exists positive constants 𝑚 and 𝑀 such that 𝑚‖𝑧‖²≤𝑧^𝑇∇²𝑓(𝑥)𝑧≤𝑀‖𝑧‖²,(3.8) for all 𝑥,𝑧∈𝑅^𝑛. Let {‖𝐵_𝑘‖} be a sequence generated by the accumulative MD method. Then, ‖𝐵_𝑘‖ is bounded above and below for all finite 𝑘, by some positive constants.

Proof. Let 𝐵^𝑘(𝑖) be the 𝑖th element of 𝐵_𝑘. Suppose 𝐵₀ is chosen such that 𝜔₁≤𝐵^0(𝑖)≤𝜔₂;𝑖=1,…,𝑛, where 𝜔₁ and 𝜔₂ are some positive constants.
Case 1. If (2.18) is satisfied, we have 𝐵₁=⎧⎪⎨⎪⎩𝜂⁰𝐵₀;if𝑟^𝑇0𝑤₀<𝑟^𝑇0𝐵₀𝑟₀𝐵₀+𝑟^𝑇0𝑤₀−𝑟^𝑇0𝐵₀𝑟_𝑘tr𝐹²⁰𝐹⁰;if𝑟^𝑇0𝑤₀≥𝑟^𝑇0𝐵₀𝑟₀.(3.9) By (2.18) and the definition of 𝜂_𝑘, one can obtain 𝜀₁𝜔₂≤𝜂₀≤1.(3.10) Thus, if 𝑟^𝑇0𝑤₀<𝑟^𝑇0𝐵₀𝑟₀, then 𝐵₁=𝜂₀𝐵₀ satisfies 𝜀₁𝜔₁𝜔₂≤𝐵^1(𝑖)≤𝜔₂.(3.11) On the other hand, if 𝑟^𝑇0𝑤₀≥𝑟^𝑇0𝐵₀𝑟₀, then 𝐵^1(𝑖)=𝐵^0(𝑖)+𝑟^𝑇0𝑤₀−𝑟^𝑇0𝐵₀𝑟₀tr𝐹²⁰𝑟^0(𝑖)²,(3.12) where 𝑟^0(𝑖) is the 𝑖th component of 𝑟₀. Letting (𝑟^0(𝑀)) be the largest component (in magnitude) of 𝑟₀, that is, (𝑟^0(𝑖))²≤(𝑟^0(𝑀))²; for all 𝑖, then it follows that ‖𝑟₀‖²≤𝑛(𝑟^0(𝑀))², and the property of 𝑟^𝑇0𝑤₀≥𝑟^𝑇0𝐵₀𝑟₀, (3.12) becomes 𝜔₁≤𝐵^0(𝑖)≤𝐵^1(𝑖)≤𝜔₂+𝑛𝜀²−𝜔₁tr𝐹²⁰𝑟^0(𝑀)⁴≤𝜔₂+𝑛𝜀₂−𝜔₁.(3.13) Hence, 𝐵^1(𝑖) is bounded above and below, for all 𝑖 in both occasions.
Case 2. If (2.18) is violated, then the updating formula for 𝐵₁ becomes 𝐵^1(𝑖)=𝜂₀𝐵^0(𝑖)+𝑠^𝑇0𝑦₀−𝜂₀𝑠^𝑇0𝐵₀𝑠₀tr𝐸²⁰𝑠^0(𝑖)²,(3.14) where 𝑠^0(𝑖) is the 𝑖th component of 𝑠₀, 𝐸₀=diag((𝑠⁰⁽¹⁾)²,(𝑠⁰⁽²⁾)²,…,(𝑠^0(𝑛))²), and 𝜂₀=min(1,𝑠^𝑇0𝑦₀/𝑠^𝑇0𝐵₀𝑠₀).
Because 𝜂₀≤1 also implies that 𝑠^𝑇0𝑦₀−𝑠^𝑇0𝐵₀𝑠₀≥0, then this fact, together with the convexity property (3.8), and the definition of 𝜂 give min1,𝑚𝜔₂𝜔¹≤𝜂₀𝐵^0(𝑖)≤𝐵^1(𝑖)≤𝐵^0(𝑖)+𝑀−𝜔¹‖𝑠₀‖²tr𝐸²⁰𝑠^0(𝑖)².(3.15) Using the similar argument as above, that is, by letting 𝑠^0(𝑀) be the largest component (in magnitude) of 𝑠₀, then it follows that min1,𝑚𝜔₂𝜔¹≤𝐵^1(𝑖)≤𝜔₂+𝑛𝑀−𝜔₁.(3.16)
Hence, in both cases, 𝐵^1(𝑖) is bounded above and below, by some positive constants. Since the upper and lower bound for 𝐵^1(𝑖) is, respectively, independent to 𝑘, we can proceed by using induction to show that 𝐵^𝑘(𝑖) is bounded, for all finite 𝑘.

4. Numerical Results

In this section, we examine the practical performance of our proposed algorithm in comparison with the BB method and standard one-step diagonal gradient-type method (MD). The new algorithms are referred to as AMD1 and AMD2 when 𝑀=𝐼 and 𝑀=𝐵_𝑘 are used, respectively. For all methods we employ Armijo line search [11] where 𝜎=0.9. All experiments in this paper are implemented on a PC with Core Duo CPU using Matlab 7.0. For each run, the termination condition is that ‖𝑔_𝑘‖≤10⁻⁴. All attempts to solve the test problems were limited to a maximum of 1000 iterations. The test problems are chosen from Andrei [13] and Moré et al. [14] collections. The detailed test problem is summarized in Table 1. Our experiments are performed on a set of 36 nonlinear unconstrained problems, and the problems vary in size from 𝑛=10 to 10000 variables. Figures 1, 2, and 3 present the Dolan and Moré [15] performance profile for all algorithms subject to the iteration, function call, and CPU time.

Table 1 

Test problem and its dimension.

Figure 1 

Performance profile based on Iteration for all problems.

Figure 2 

Performance profile based on function call.

Figure 3 

Performance profile based on CPU time per iteration.

From Figure 1, we see that AMD2 method is the top performer, being more successful than other methods in the number of iteration. Figure 2 shows that AMD2 method requires the fewest function calls. From Figure 3, we observe that the AMD2 method is faster than MD and AMD1 methods and needs reasonable time to solve large-scale problems when compared to the BB method. At each iteration, the proposed method does not require more storage than classic diagonal updating methods. Moreover, a higher-order accuracy in approximating the Hessian matrix of the objective function makes AMD method need less iterations and less function evaluation. The numerical results by the tests reported in Figures 1, 2, and 3 demonstrate clearly the new method AMD2 shows significant improvements, when compared with BB, MD, and AMD1. Generally, 𝑀=𝐵_𝑘 performs better than 𝑀=𝐼. It is most probably due to the fact that 𝐵_𝑘 is a better Hessian approximation than the identity matrix 𝐼.

5. Conclusion

In this paper, we propose a new two-step diagonal gradient method as view of accumulative approach for unconstrained optimization. The new parameterization for multistep diagonal gradient-type method is developed via employing accumulative approach. The new technique is devised for interpolating curves which are the basis of multistep approach. Numerical results show that the proposed method is suitable to solve large-scale unconstrained optimization problems and more stable than other similar methods in practical computation. The improvement that our proposed methods bring does come at a complexity cost of 𝑂(𝑛) while others are about 𝑂(𝑛²) [9, 10].