A new search direction for full-Newton step infeasible interior-point method in linear optimization
aa r X i v : . [ m a t h . O C ] F e b A new search direction for full-Newton stepinfeasible interior-point method in linearoptimization
B. KheirfamDepartment of Applied MathematicsAzarbaijan Shahid Madani UniversityTabriz, Iran [email protected]
Abstract
In this paper, we study an infeasible interior-point method for linear op-timization with full-Newton step. The introduced method uses an algebraicequivalent transformation on the centering equation of the system which de-fines the central path. We prove that the method finds an ε -optimal solutionof the underlying problem in polynomial time. keywords : Linear optimization, infeasible interior-point methods, new search direc-tions, polynomial complexity. AMS : 90C05, 90C51.
Interior-point methods (IPMs) for linear optimization (LO) began when Karmarkar[6] published his exceptional paper in 1984. After that, several variants of this al-gorithm were presented. In the meantime, we can talk about feasible and infeasibleIPMs. In feasible IPMs we presume that a strictly feasible point is at hand which thealgorithm can be immediately beginning. Usually find such a starting point is not sim-ple. In that case an infeasible IPM (IIPM) should be used. These methods begin froman arbitrary positive point and try to reach both feasibility and optimality. IIPMswere first introduced by Lustig [18] and Tanabe [22]. The first feasible IPM withfull-Newton step for LO was presented by Roos et al. [21]. Determining the searchdirections plays a very important role in IPMs. In 2003, Darvay [2] utilizes the AETtechnique on the centering equation of the system defining the central path for LO.He uses the square root function in the AET strategy and then applies the Newtonmethod for obtain the search directions. This method is extended in [1, 23, 24, 25], re-spectively, to convex quadratic optimization (CQO), second-order cone optimization1SOCO), symmetric optimization (SO) and the Cartesian P ∗ ( κ ) linear complementar-ity problem (LCP). Kheirfam and Haghighi [16] have proposed an IPM for P ∗ ( κ )-LCPwhich uses the function ψ ( t ) = √ t/ √ t ) in the AET technique. An infeasibleversion of the method proposed in [21] has presented by Roos in [19] which needs afeasibility step and three centering steps in each main iteration. Some generalizationsand versions of the method can be seen in [15, 8, 9, 17, 27, 10]. The author is improvedthis algorithm so that the algorithm performs only one feasibility step in each itera-tion and does not need centering steps [20]. Kheirfam [11, 12, 13, 14] extended thealgorithm proposed in [20] to HLCP, the Cartesian P ∗ ( κ )-LCP, the convex quadraticsymmetric cone optimization (CQSCO) and SO. By considering the AET techniquebased on the function ψ ( t ) = t − √ t , Darvay et al. [3] have introduced a full-Newtonstep IPM for LO. Kheirfam [7] has presented an infeasible version of this algorithmfor SDLCP. Darvay et al. [5] published a corrector-predictor IPM (CP-IPM) for LOusing the function ψ ( t ) = t − √ t for AET. Darvay and Tak´acs [4] proposed an IPMfor LO based on a new type of AET on the centering equation of the central path.Motivated by the aforementioned works, in this paper we aim to present a full-Newton step IIPM for LO using the AET ψ ( t ) = ψ ( √ t ) for the centering equationof the central path. The method uses the function ψ ( t ) = t in order to determinethe new search directions and performs only one feasibility step in a main iteration.In fact, our method is an infeasible version of the method proposed in [4]. We provethat the proposed algorithm enjoys the best-known iteration complexity for IIPMs.The paper is organized in the following way. In the next section, we remember theproblem pair (P) and (D). We state the perturbed problems corresponding to (P) and(D) and then provided the central path. In Sect. 3, the new search directions basedon the new type of AET using the function ψ ( t ) = t is discussed, and finally thealgorithm is presented. Section 4 consists of the complexity analysis of the introducedIIPM with the new search directions. In Section 5, some concluding remarks arefollowed. Let us consider the LO problem in the standard form( P ) min { c T x : Ax = b, x ≥ } , where A ∈ R m × n with rank ( A ) = m, b ∈ R m and c ∈ R n . The dual of this problemcan be written in the following standard form:( D ) max { b T y : A T y + s = c, s ≥ } . In accordance with the routine of IIPMs, we consider the starting point ( x , y , s ) = ξ ( e, , e ) such that k ( x ∗ ; s ∗ ) k ∞ ≤ ξ for some primal-dual optimal solution ( x ∗ , y ∗ , s ∗ ),where e is the all-one vector and ξ is a positive scalar. It should be noted that forthe optimal solution ( x ∗ , y ∗ , s ∗ ) the inequality k ( x ∗ ; s ∗ ) k ∞ ≤ ξ is true if and only if0 ≤ x ∗ ≤ ξe, ≤ s ∗ ≤ ξe. (1)2or an IIPM, a triple ( x, y, s ) is called an ε -solution of (P) and (D) ifmax (cid:8) x T s, k b − Ax k , k c − A T y − s k (cid:9) ≤ ǫ, where ε is a accuracy parameter. Following [19], for any 0 < ν ≤ ν ) and (D ν ) as follows:( P ν ) min { ( c − νr c ) T x : b − Ax = νr b , x ≥ } , ( D ν ) max { ( b − νr b ) T y : c − A T y − s = νr c , s ≥ } , where r b := b − Aξe and r c := c − ξe. It is simply seen that ( x , y , s ) = ξ ( e, , e ) isa feasible solution of the problem pair (P ν ) and (D ν ) if ν = 1. We conclude that if ν = 1, then (P ν ) and (D ν ) satisfy the interior point condition (IPC). We recall thefollowing lemma. Lemma 1. (Theorem 5.13 in [26])
The original problems, (P) and (D) are feasibleif and only if for each ν satisfying < ν ≤ the perturbed problems (P ν ) and (D ν ) satisfy the IPC. In the view of Lemma 1, we assume that the original problem pair (P) and (D) isfeasible and ν ∈ (0 , ν ) and (D ν ) exists; b − Ax = νr b , x ≥ ,c − A T y − s = νr c , s ≥ ,xs = µe, (2)has a unique solution ( x ( µ, ν ) , y ( µ, ν ) , s ( µ, ν )) for every µ >
0. This solution consistsof the µ -centers of the perturbed problems (P ν ) and (D ν ). Note that for x, s > µ > xs = µe ⇔ xsµ = e ⇔ r xsµ = e ⇔ xsµ = r xsµ . (3)Now the perturbed central path can be equivalently stated as follows: b − Ax = νr b , x ≥ ,c − A T y − s = νr c , s ≥ ,xsµ = r xsµ . (4)In the sequel, the parameters µ and ν always satisfy the relation µ = νµ = νξ . In accordance with the Darvay’s idea, we consider the function ψ defined and contin-uously differentiable on the interval ( k , ∞ ), where 0 ≤ k <
1, such that 2 tψ ′ ( t ) − ′ ( t ) > , ∀ t > k . Now, if we apply the AET method to (4), then we get b − Ax = νr b , x ≥ , (5) c − A T y − s = νr c , s ≥ , (6) ψ (cid:16) xsµ (cid:17) = ψ (cid:16)r xsµ (cid:17) . (7)Let ( x, y, s ) be a feasible solution of the perturbed pair (P ν ) and (D ν ). We considerthe notation f ( x, y, s ) = ν + r b − b − Axν + r c − c + A T y + sψ (cid:16) xsµ (cid:17) − ψ (cid:16)r xsµ (cid:17) = 0 , where ν + = (1 − θ ) ν and θ ∈ (0 , J f ( x, y, s ) ∆ x ∆ y ∆ s = − f ( x, y, s ) , where J f ( x, y, s ) denotes the Jacobian matrix of f at ( x, y, s ). After some computa-tions, we obtain the following system: A ∆ x = θνr b ,A T ∆ y + ∆ s = θνr c , µ (cid:0) s ∆ x + x ∆ s (cid:1) = − ψ ( xsµ ) + ψ ( q xsµ ) ψ ′ ( xsµ ) − √ xsµ ψ ′ ( q xsµ ) . (8)Defining the scaled search directions d x := v ∆ xx , d s := v ∆ ss , where v = r xsµ , (9)we can give the scaled form of system (8):¯ Ad x = θνr b , ¯ A T ∆ yµ + d s = θνvs − r c ,d x + d s = p v , (10)where p v := 2 ψ ( v ) − ψ ( v )2 vψ ′ ( v ) − ψ ′ ( v ) , and ¯ A := A diag (cid:0) xv (cid:1) . If we use the function ψ : ( √ , ∞ ) → R , ψ ( t ) = t introduced in [4], then we obtain p v = v − v v − e . (11)4fter a full-Newton step, the new iterate is given by x + := x + ∆ x, y + := y + ∆ y, s + := s + ∆ s. (12)Furthermore, in each iteration of the algorithm, a quantity is needed to measure howfar an iterate is from the central path. We consider the proximity measure defined by δ ( v ) := δ ( x, s ; µ ) = k p v k (cid:13)(cid:13)(cid:13) v − v v − e (cid:13)(cid:13)(cid:13) , (13)which was first suggested for a feasible IPM in [4].Let q v = d x − d s . Then d x = p v + q v , d s = p v − q v , d x d s = p v − q v , (14)and k q v k k d x − d s k k d x + d s k − d Tx d s = k p v k − d Tx d s . (15)Suppose that for some µ ∈ (0 , µ ], our algorithm begins from a feasible solution( x, y, s ) of the problem pair (P ν ) and (D ν ) with ν = µµ , and such that δ ( x, s ; µ ) ≤ τ, τ ∈ (0 , x + , y + , s + ) of (P + ν ) and(D + ν ), where ν + = (1 − θ ) ν, θ ∈ (0 , µ is decreased to µ + = (1 − θ ) µ and such that δ ( x + , s + ; µ + ) ≤ τ . This procedure is repeated until an ε -solution isfound. We are now in a position to state the theoretical framework of the infeasibleinterior-point algorithm as follows: Algorithm1 : an infeasible interior − point algorithm Input :Accuracy parameter ε > θ, < θ < τ > . begin x := ξe ; y := 0; s := ξe ; µ := νξ ; ν = 1; while max( x T s, k r b k , k r c k ) > ε dobegin solve the system (10) and use (9) to obtain (∆ x, ∆ y, ∆ s );( x, y, s ) := ( x, y, s ) + (∆ x, ∆ y, ∆ s );update of µ and ν : µ := (1 − θ ) µ ; ν := (1 − θ ) ν ; endend . Analysis of the algorithm
Here, we will prove that Algorithm 1 is well-defined. The main goal of our analysis isto find some values for the parameters τ and θ such that x + > s + >
0, and wehave δ ( x + , s + ; µ + ) ≤ τ. In the following section, we obtain an upper bound for theproximity measure after an iteration of the algorithm. δ ( v + ) In the next lemma, we give a condition on the proximity measure which ensuresthe feasibility of a full-Newton step. In what follows, we use the notation ω = (cid:0) k d x k + k d s k (cid:1) . Lemma 2.
The iterate ( x + , y + , s + ) with v > √ e is strictly feasible if δ ( v ) + ω < .Proof. Let 0 ≤ α ≤
1. We define x ( α ) := x + α ∆ x and s ( α ) := s + α ∆ s. Using (9),the third equation of (10) and (14) one can find x ( α ) s ( α ) µ = xsv ( v + αd x )( v + αd s ) = v + αv ( d x + d s ) + α d x d s = (1 − α ) v + α ( v + vp v ) + α (cid:16) p v − q v (cid:17) (16) ≥ (1 − α ) v + α e + α p v − α q v , where the inequality is due to α ≥ α and the following inequality: v + vp v − e = v + v − v v − e − e = v v − e − e = ( v − e ) v − e ≥ . (17)The inequality x ( α ) s ( α ) > (cid:13)(cid:13)(cid:13) − p v q v (cid:13)(cid:13)(cid:13) ∞ ≤ (cid:13)(cid:13)(cid:13) p v (cid:13)(cid:13)(cid:13) ∞ + (cid:13)(cid:13)(cid:13) q v (cid:13)(cid:13)(cid:13) ∞ ≤ k p v k k q v k δ ( v ) − d Tx d s ≤ δ ( v ) + k d x kk d s k ≤ δ ( v ) + ω < , where the equality is due to (15), the third inequality uses from the Cauchy-Schwarzinequality and the last inequality holds due to the assumption of the lemma. Thus, x ( α ) s ( α ) >
0, for 0 ≤ α ≤ x ( α ) and s ( α ) do not change sign on the interval[0 , x (0) = x > s (0) = s > x (1) = x + > s (1) = s + >
0. Thus, the proof is completed.In correspondence to the definition (13), we have δ ( v + ) = δ ( x + , s + ; µ + ) = 12 (cid:13)(cid:13)(cid:13) v + − v v − e (cid:13)(cid:13)(cid:13) , where v + = r x + s + µ + . emma 3. Let δ ( v ) + ω < (1 − θ ) and v > √ e . Then, v + > √ e and δ ( v + ) ≤ p − δ ( v ) − ω (cid:0) θ √ n + 10 δ ( v ) + ω (cid:1) √ − θ (2(1 − δ ( v ) − ω ) − (1 − θ )) . Proof.
Let α = 1. Then from (16) it follows that v = x + s + µ + = v + vp v + p v − q v − θ = e + ( v − e ) v − e + p v − q v − θ = e + (cid:0) v − ev (cid:1) p v − q v − θ ≥ e − q v − θ , where the second equality is due to (17) and the inequality follows from the fact that9 v − e ≥ . e >
0. Consequently, we havemin( v + ) ≥ s − k q v k ∞ − θ ≥ s − k q v k − θ ≥ r − δ ( v ) − ω − θ , (18)where the last inequality follows from (15) and the Cauchy-Schwarz inequality.From δ ( v ) + ω < (1 − θ ) it follows that min( v + ) > √ , hence v + > √ e . Now,we have δ ( v + ) = 12 (cid:13)(cid:13)(cid:13) v + − v v − e (cid:13)(cid:13)(cid:13) = 12 (cid:13)(cid:13)(cid:13) v + v − e (cid:0) e − v (cid:1)(cid:13)(cid:13)(cid:13) ≤ min( v + )2(2 min( v + ) − (cid:13)(cid:13) e − v (cid:13)(cid:13) ≤ p (1 − θ )(1 − δ ( v ) − ω )2(2(1 − δ ( v ) − ω ) − (1 − θ )) (cid:13)(cid:13) e − v (cid:13)(cid:13) . (19)On the other hand, one has (cid:13)(cid:13) e − v (cid:13)(cid:13) = (cid:13)(cid:13)(cid:13) e + (cid:0) v − ev (cid:1) p v − q v − θ − e (cid:13)(cid:13)(cid:13) ≤ − θ (cid:16) k θe k + (cid:13)(cid:13)(cid:13)(cid:16) v − ev (cid:17) p v − q v (cid:13)(cid:13)(cid:13)(cid:17) ≤ − θ (cid:16) θ √ n + 9 k p v k k q v k (cid:17) = 11 − θ (cid:0) θ √ n + 10 δ ( v ) + ω (cid:1) . Substituting this bound into (19) gives us exactly the desired result. Thus, the proofis completed. 7 .2 Upper bound for ω Following [20], let N := { ζ : ¯ Aζ = 0 } denote the null space of the matrix ¯ A . Then, the { ζ : ¯ Aζ = θνr b } affine space equals N + d x . Since the row space of ¯ A is the orthogonalcomplement N ⊥ of N , thus d s ∈ θνvs − r c + N ⊥ . Also note that
N ∩ N ⊥ = { } , andthe affine spaces N + d x and N ⊥ + d s meet in a unique point q . Applying a similarargument to Lemma 3.4 in [20], we can find2 ω ≤ k q k + (cid:18) k q k + (cid:13)(cid:13)(cid:13) v − v v − e (cid:13)(cid:13)(cid:13)(cid:19) = k q k + (cid:0) k q k + 2 δ ( v ) (cid:1) . (20)Again from [20], we have k q k ≤ θ (cid:0) n + k v k (cid:1) min( v ) . (21)By definition (13) , we have2 δ ( v ) = (cid:13)(cid:13)(cid:13) v − v v − e (cid:13)(cid:13)(cid:13) = (cid:13)(cid:13)(cid:13) v + v v − e ( e − v ) (cid:13)(cid:13)(cid:13) ≥ (cid:13)(cid:13) e − v (cid:13)(cid:13) ≥ (cid:0) k v k − k e k (cid:1) , which implies k v k ≤ k e k + 4 δ ( v ) = √ n + 4 δ ( v ) . Furthermore, we have 4 δ ( v ) ≥ k e − v k ≥ | − v i | , i = 1 , . . . , n. This gives min( v ) ≥ − δ ( v ). Combining these two inequalities with (21), we willget k q k ≤ θ (cid:16) n + (cid:0) √ n + 4 δ ( v ) (cid:1) (cid:17) − δ ( v ) . (22) θ and τ In this section, we require finding values θ and τ such that if δ ( v ) ≤ τ holds, then δ ( v + ) ≤ τ . From Lemma 3, it suffices to have p − δ ( v ) − ω (cid:0) θ √ n + 10 δ ( v ) + ω (cid:1) √ − θ (2(1 − δ ( v ) − ω ) − (1 − θ )) ≤ τ, (23)provided that δ ( v ) + ω < (1 − θ ). One can easily see the right-hand-side of (22)is monotonically increasing with respect to δ ( v ) <
1. Hence, invoking δ ( v ) ≤ τ , wehave k q k ≤ θ (cid:16) n + (cid:0) √ n + 4 τ (cid:1) (cid:17) − τ .
8y substituting the above result into (20) and using again δ ( v ) ≤ τ , we obtain ω ≤ (cid:20)(cid:18) θ (cid:0) n + (cid:0) √ n + 4 τ (cid:1) (cid:1) − τ (cid:19) + (cid:18) θ (cid:0) n + (cid:0) √ n + 4 τ (cid:1) (cid:1) − τ + 2 τ (cid:19) (cid:21) =: f ( τ ) . We claim that χ ( t ) := √ − t − t ) − (1 − θ ) , ≤ t ≤
12 (1 − θ ) , (24)is increasing. Hence, 0 ≤ δ ( v ) + ω ≤ τ + f ( τ ) implies χ ( δ ( v ) + ω ) ≤ χ ( τ + f ( τ )).Therefore, δ ( v ) + ω ≤ (1 − θ ) and (23) will certainly hold if τ + f ( τ ) ≤
12 (1 − θ ) , y ( τ ) := χ ( τ + f ( τ ))( θ √ n + 10 τ + f ( τ ))2 √ − θ ≤ τ. If we take τ = and θ = n , n ≥
4, then τ + f ( τ ) = 0 . < . ≤ (1 − θ )and y ( τ ) ≤ . < . Hence, we may state the following result. Lemma 4. If τ = and θ = n , n ≥ , then δ ( v ) ≤ τ implies δ ( v + ) ≤ τ. Lemma 4 establishes the proposed algorithm is well-defined, in the sense that theproperty δ ( x, s ; µ ) := δ ( v ) ≤ τ is maintained in all iterations.In each main iteration, both the barrier parameter µ and the norms of the residualvectors are reduced by the factor 1 − θ . Hence, the total number of main iterationsis bounded above by 1 θ log max { nξ , k r b k , k r c k} ε . Now, we state our main result.
Theorem 5. If (P) and (D) are feasible and ξ > such that k ( x ∗ ; s ∗ ) k ∞ ≤ ξ forsome optimal solutions x ∗ of (P) and ( y ∗ , s ∗ ) of (D) , then after at most n log max { nξ , k r b k , k r c k} ǫ iterations, the algorithm finds an ǫ -optimal solution of (P) and (D) . The method presented in this paper is a full-Newton step IIPM for LO based on theAET proposed in [4]. The method is used in each iteration only one feasibility step.Our method analysis is different from the existing IIPMs based on the AET becauseit uses a different AET. The obtained complexity bound coincides with the currentbest-known theoretical iteration bound for IIPMs.9 eferences [1] Achache, M.,
A new primal-dual path-following method for convex quadratic pro-gramming , Comp. Appl. Math., 25(1), 97-110 (2006).[2] Darvay, Zs.,
New interior-point algorithms in linear programming , Adv. Model.Optim., 5(1), 51-92 (2003).[3] Darvay, Zs., Papp, I.M., Tak´acs, P.R.,
Complexity analysis of a full-Newton stepinterior-point method for linear optimization , Period. Math. Hung., 73, 27-42(2016).[4] Darvay, Zs., Tak´acs, P.R.,
New method determining search directions for interior-point algorithms in linear optimization , Optim. Lett., 12, 1099-1116 (2018).[5] Zs. Darvay, Zs., Ill´es, T., Kheirfam, B., Rig´o, P.R.,
A corrector-predictor interior-point method with new search direction for linear optimization , Cent. Eur. J.Oper. Res., 28(3), 1123-1140 (2020).[6] Karmarkar, N.K,
A new polynomial time algorithm for linear programming , Com-binatorica, 4, 375-395 (1984).[7] Kheirfam, B.,
An infeasible interior point method for the monotone SDLCP basedon a transformation of the central path , J. Appl. Math. Comput., 57(1), 685-702(2018).[8] Kheirfam, B.,
A new complexity analysis for full-Newton step infeasible interior-point algorithm for P ∗ ( κ ) -horizontal linear complementarity problems , J. Optim.Theory Appl., 161(3), 853-869 (2014).[9] Kheirfam, B., A full Nesterov-Todd step infeasible interior-point algorithm forsymmetric optimization based on a specific kernel function , NACO, 3(4), 601-614 (2013).[10] Kheirfam, B.,
A new infeasible interior-point method based on Darvay’s techniquefor symmetric optimization , Ann. Oper. Res., 211(1), 209-224 (2013).[11] Kheirfam, B.,
An improved full-Newton step O ( n ) infeasible interior-pointmethod for horizontal linear complementarity problem , Numer. Algorithms,71(3), 491-503 (2016).[12] Kheirfam, B., A full step infeasible interior-point method for Cartesian P ∗ ( κ ) -SCLCP , Optim. Lett., 10(3), 591-603 (2016).[13] Kheirfam, B., An improved and modified infeasible interior-point method for sym-metric optimization , Asian-Eur. J. Math., 9(2), 1650059 (13 pages) (2016).[14] Kheirfam, B.,
An infeasible full-NT step interior point algorithm for
CQSCO,Numer. Algorithms, 74(1), 93-109 (2017).1015] Kheirfam, B., Mahdavi-Amiri, N.,
A full Nesterov-Todd step infeasible interior-point algorithm for symmetric cone linear complementarity problem , Bull. IranianMath. Soc., 40(3), 541-564 (2014).[16] Kheirfam, B., Haghighi, M.,
A full-Newton step feasible interior-point algorithmfor P ∗ ( κ ) -LCP based on a new search direction , Croat. Oper. Res. Rev., 7(2),277-290 (2016).[17] Liu, Z., Sun, W., Tian, F., A full-Newton step O ( n ) infeasible interior-pointalgorithm for linear programming based on kernel function , Appl. Math. Optim.,60, 237-251 (2009).[18] Lustig, I.J., Feasibility issues in a primal-dual interior-point method for linearprogramming , Math. Program., 49(1-3), 145-162 (1991).[19] Roos, C.,
A full-newton step O ( n ) infeasible interior- point algorithm for linearoptimization , SIAM J. Optim., 16(4), 1110-1136 (2006).[20] Roos, C., An improved and simplified full-Newton step O ( n ) infeasible interior-point method for linear optimization , SIAM J. Optim., 25(1), 102-114 (2015).[21] Roos, C., Terlaky, T., Vial, J-Ph., Theory and Algorithms for Linear Optimiza-tion. An Interior-Point Approach , John Wiley and Sons, Chichester, UK, (1997).[22] Tanabe, K.,
Centered Newton method for linear programming: Interior and ’ex-terior’ point method (in Janpanese). In: New Methods for Linear Programming.K. Tone (Ed.) 3, pages 98-100, 1990.[23] Wang, G.Q., Bai, Y.Q.,
A primal-dual path-following interior-point algorithmfor second-order cone optimization with full Nesterov-Todd step , Appl. Math.Comput., 215(3), 1047-1061 (2009).[24] Wang, G.Q., Bai, Y.Q.,
A new full Nesterov-Todd step primal-dual path-followinginterior-point algorithm for symmetric optimization , J. Optim. Theory Appl.,154(3), 966-985 (2012).[25] Wang, G.Q., Fan, X.J., Zhu, D.T., Wang, D.Z.,
New complexity analysis of afull-Newton step feasible interior-point algorithm for P ∗ ( κ ) -LCP , Optim. Letters,9(6), 1105-1119 (2015).[26] Ye, Y., Interior point algorithms . Wiley-Interscience Series in Discrete Mathe-matics and Optimization. John Wiley & Sons Inc., New York, 1997. Theory andanalysis, A Wiley-Interscience Publication.[27] Zhang, L., Sun, L., Xu, Y.,