Primal-dual path-following methods and the trust-region updating strategy for linear programming with noisy data
JJournal of XXX manuscript No. (will be inserted by the editor)
Perturbed Newton Method with Trust-region Time-steppingSchemes for Linear Programming with Uncertain Data
Xin-long Luo ∗ · Yi-yan Yao
Received: date / Accepted: date
Abstract
In this article, we consider the path-following method based on the per-turbed Newton flow with the new trust-region time-stepping scheme for the standardlinear programming problem. For the problem with deficient rank matrix and thenoise right-hand-side vector, we also give the pre-processing method based on thesingular value decomposition. Then, we analyze the global convergence of the newmethod when the initial point is strictly primal-dual feasible. Finally, we test the newmethod for some problems with deficient rank matrices, and compare it with otherpopular interior-point methods such as the path-following method (the subroutinepathfollow.m coded by M. C. Ferris [14, 16]) and Mehrotra’s predictor-corrector al-gorithm (the built-in subroutine linprog.m of the MATLAB environment, which wasimplemented by S. Mehrotra and Y. Zhang [26, 27, 36]). Numerical results show thatthe new method is more robust than those methods for the large-scale deficient-rankproblems without sacrificing its computational efficiency.
Keywords
Continuation Newton method · trust-region technique · linear program-ming · deficient rank · path-following method · uncertain data Mathematics Subject Classification (2010) · · · Xin-long LuoCorresponding author. School of Information and Communication Engineering, Beijing University ofPosts and Telecommunications, P. O. Box 101, Xitucheng Road No. 10, Haidian District, 100876, Bei-jing China, E-mail: [email protected] YaoSchool of Information and Communication Engineering, Beijing University of Posts and Telecommu-nications, P. O. Box 101, Xitucheng Road No. 10, Haidian District, 100876, Beijing China, E-mail:[email protected] a r X i v : . [ m a t h . O C ] J un Luo, Yao
In this article, we are mainly concerned with the linear programming problem withuncertain data as follows:min x ∈ ℜ n c T x , subject to Ax = b , x ≥ , (1)where c and x are vectors in ℜ n , b is a vector in ℜ m , and A is an m × n matrix. For theproblem (1), there are many efficient methods to solve it such as the simplex methodsand the interior-point methods [14, 16, 28, 34, 36]. Those methods are all assumed tohave the consistent system of constraints, i.e. rank ( A , b ) = rank ( A ) .However, for the real-world problem, since it may exist the redundant constraintsand the measurement errors, the rank of matrix A is deficient and the right-hand-side vector b has small noise, which may result in the inconsistent system of con-straints [8]. In order to handle this special case, we give a preprocessing methodbased on the singular value decomposition. Then, according to the first order KKTconditions of the linear programming, we convert the processed problems into theequivalent problem of nonlinear equations with nonnegative constraints. Based onthe system of nonlinear equations with nonnegative constraints, we consider a spe-cial continuous Newton flow with nonnegative constraints, which has the nonnegativesteady-state solution for any nonnegative initial point. In order to improve the robustand efficiency of tracing the trajectory, we construct a perturbed Newton method withthe new time-stepping scheme based on the trust-region technique.The rest of this article is organized as follows. In the next section, we give thecontinuation Newton method based on the perturbed Newton flow with the new trust-region time-stepping scheme to follow the trajectory and obtain its steady-state solu-tion. Then, in section 3, we analyze the global convergence of the new method whenthe initial point is strictly feasible. In section 4, for the problems with deficient rankmatrices and the noise right-hand-side vectors, some promising numerical results ofthe new method also reported, with comparison to other popular interior-point meth-ods such as the path-following method (the subroutine pathfollow.m coded by M. C.Ferris [14, 16]) and Mehrotra’s predictor-corrector algorithm (the built-in subroutinelinprog.m of the MATLAB environment, which was implemented by S. Mehrotra andY. Zhang [26, 27, 36]). Numerical results show that the new method is more robustthan those methods for the large-scale deficient-rank problems without sacrificing itscomputational efficiency. Finally, some discussions are also given in section 5. (cid:107) · (cid:107) denotes the Euclidean vector norm or its induced matrix norm through the paper. x ∗ if and only if it satisfies the following Karush-Kuhn-Tucker conditions (see pp. erturbed Newton Method with Trust-region Time-stepping Schemes for Linear Programming 3 Ax − b = , (2) A T y + s − c = , (3) XSe = , (4) ( x , s ) ≥ , (5)where X = diag ( x ) , S = diag ( s ) , and e = ( , . . . , ) T . (6)For the sake of convenience, we regard the optimality conditions (2)-(5) as a systemof nonlinear equations with nonnegative constraints as follows: F ( z ) = Ax − bA T y + s − cXSe = , ( x , s ) ≥ , and z = ( x , y , s ) . (7)It is not difficult to know that the Jacobian matrix J ( z ) of F ( z ) defined by equation(7) has the following form: J ( z ) = A A T IS X . (8)From the third block XSe = x ∗ i = s ∗ i = z ∗ , where x ∗ i and s ∗ i are the elements of x ∗ and s ∗ , respectively. Thus, theJacobian matrix J ( z ∗ ) defined by (8) may be singular, which results in numericaldifficulties around the solution of the nonlinear system (7) for the Newton’s methodor its variants. Therefore, we consider its small positive perturbed system, which issimilar to the primal-dual central path [14, 28, 34] as follows: F µ ( z ) = F ( z ) − µ e = , ( x , s ) > , µ > z = ( x , y , s ) . (9)We define the strictly feasible region F of the linear programming problem (1)as F = (cid:8) ( x , y , s ) | Ax = b , A T y + s = c , ( x , s ) > (cid:9) . (10)Then, when it exist a strictly feasible interior point ( ¯ x , ¯ y , ¯ s ) ∈ F and the rank ofmatrix A is full, the perturbed system (9) has a unique solution (see Theorem 2.8, p. Luo, Yao
39, [34]). Its existence can be derived by the implicit theorem [13] and its uniquenesscan be proved via considering the strict convexity of the following penalty problem[15]: min c T x − µ n ∑ i = log ( x i ) subject to Ax = b , (11)where µ is a positive parameter.According to the duality theorem of the linear programming (Theorem 13.1, pp.368-369 in [28]), we know c T x ≥ b T y ∗ = c T x ∗ ≥ b T y (12)for any primal-dual feasible solution ( x , y , s ) , where the triple ( x ∗ , y ∗ , s ∗ ) is a primal-dual optimal solution of nonlinear system (7). Thus, when µ > z ∗ µ of perturbed system (9) is an approximation solution of nonlinear system (7).Consequently, x ∗ µ is an approximation of the optimal solution of the original linearprogramming problem (1). It can be illustrated as follows. Since z ∗ µ is the primal-dualfeasible, from inequality (12), we have c T x ∗ µ ≥ b T y ∗ = c T x ∗ ≥ b T y ∗ µ (13)and 0 ≤ ( x ∗ µ ) T s ∗ µ = c T x ∗ µ − b T y ∗ µ . (14)From inequalities (13)-(14), it is not difficult to obtain | c T x ∗ µ − c T x ∗ | ≤ ( x ∗ µ ) T s ∗ µ = n µ . (15)Therefore, x ∗ µ is an approximation of the optimal solution of the original linear pro-gramming problem (1).If we consider the Newton method with a line search strategy for the perturbedsystem (9) [12, 28], we have z k + = z k − ∆ t k J ( z k ) − F µ ( z k ) , (16)where J ( z k ) is the Jacobian matrix of function F µ ( · ) at z k . In equation (16), if weregard z k + = z ( t k + ∆ t k ) , z k = z ( t k ) and let ∆ t k →
0, we obtain the continuous Newtonflow [4, 5, 9, 24, 32] of the perturbed system (9) with the constraints as follows : dz ( t ) dt = − J ( z ) − F µ ( z ) , z = ( x , y , s ) and ( x , s ) > . (17)Actually, if we apply an iteration with the explicit Euler method [31, 35] for the con-tinuous Newton flow (17), we obtain the damped Newton method (16).Since the Jacobian matrix J ( z ) = F (cid:48) µ ( z ) may be singular, we reformulate the con-tinuous Newton flow (17) as the following more general formula [5, 32]: − J ( z ) dz ( t ) dt = F µ ( z ) , z = ( x , y , s ) and ( x , s ) > . (18)The Newton flow (18) has some nice properties. We state them as the following prop-erty 1 [5, 24, 32]. erturbed Newton Method with Trust-region Time-stepping Schemes for Linear Programming 5 Property 1 If z ( t ) is the solution of the continuous Newton flow (18), f ( z ( t )) = (cid:107) F µ ( z ) (cid:107) converges to zero as t tends to infinity. Namely, for every limit point z ∗ of z ( t ) , it is also a solution of nonlinear equations (9). Furthermore, every element F µ , i ( z ) of function F µ ( z ) has the linear convergence rate exp ( − t ) . If the Jacobianmatrix function J ( z ) of function F µ ( z ) is nonsingular, z ( t ) can not converge to itsequilibrium z ∗ on finite interval. Proof.
For the completeness, we restate its proof as follows. Assume that z ( t ) isthe solution of the continuous Newton flow (18), then we have ddt (cid:0) exp ( t ) F µ ( z ) (cid:1) = exp ( t ) J ( z ) dz ( t ) dt + exp ( t ) F µ ( z ) = , which gives F µ ( z ) = F µ ( z ) exp ( − t ) . (19)From equation (19), it is easy to see that every element F µ , i ( z ) of function F µ ( z ) converges to zero with linear convergence rate exp ( − t ) when t tends to infinity. Thus,if the solution z ( t ) of the continuous Newton flow (18) belongs to a compact set, itconverges to a limit point z ∗ when t tends to infinity, and this limit point z ∗ is also asolution of nonlinear equations (9).Furthermore, if the Jacobian matrix J ( z ) of function F µ ( z ) is nonsingular, equa-tion (18) is equivalent to the continuous Newton flow (17). Thus, from the dynamicalproperty of the ordinary differential equation (ODE), we know that the ODE (17) hasa unique solution z ( t ) and it can not converge to its equilibrium z ∗ on finite interval(see pp. 79-82, [30]). (cid:117)(cid:116) The solution x ( t ) of the continuous Newton flow (17) has a nice property. Namely,the inverse Jacobian matrix J ( x ) − can be regarded as the preconditioner of F µ ( x ) such that the solution elements x i ( t ) ( i = , , . . . , n ) have the roughly same con-vergence rates and it mitigates the stiff property of the ODE (ill-conditioned prop-erty) [24]. This property is very useful since it makes us adopt the explicit ODEmethod to trace the trajectory of the Newton flow.2.2 The perturbed Newton methodFrom subsection 2.1, we know that the continuous Newton flow (18) has the niceglobal convergence property. On the other hand, when the Jacobian matrix J ( x ) is sin-gular or nearly singular, the ODE (18) is the system of differential-algebraic equations[3,6,18] and its trajectory can not be efficiently followed by the general ODE methodsuch as the backward differentiation formulas (the built-in subroutine bdf15s.m of theMATLAB environment [26, 31]). Thus, we need to construct the special method tohandle this problem and we expect that the new method has the same global con-vergence as the homotopy continuation methods [1, 25] and the fast convergencenear the solution x ∗ as the traditional optimization methods. In order to attain these Luo, Yao two aims, we consider the continuation Newton method and construct a new time-stepping scheme based on the trust-region technique for problem (18).If we apply the implicit Euler method to the continuous Newton flow (18) [3, 6],we obtain J ( z k + ) z k + − z k ∆ t k = − F µ k + ( z k + ) . (20)The scheme (20) is an implicit formula and it needs to solve a system of nonlinearequations at every iteration. To avoid solving the system of nonlinear equations, wereplace J ( z k + ) with J ( z k ) and replace F µ k + ( z k + ) with F µ k + ( z k ) + J ( z k )( z k + − z k ) in equation (20). Thus, we obtain the perturbed Newton method as follows: z k + = z k − ∆ t k + ∆ t k J ( z k ) − F µ k + ( z k ) . (21)The method (21) is similar to the damped Newton method (16) if we regard ∆ t k / ( + ∆ t k ) in equation (21) as ∆ t k in equation (16). However, from the view ofthe ODE method, they are very different. The damped Newton method (16) is ob-tained by the explicit Euler scheme applied to the continuous Newton flow (18), andits time-stepping ∆ t k is restricted by the numerical stability [18, 31, 35]. Namely, thelarge time-stepping size ∆ t k can not be adopted in the steady-state phase. The method(21) is obtained by the semi-implicit Euler scheme applied to the continuous Newtonflow (18), and its time-stepping size ∆ t k is not restricted by the numerical stability.Therefore, the large time-stepping size can be adopted in the steady-state phase, theperturbed Newton method (21) of which mimics the Newton method and is expectedto have the fast convergence rate. The most of all, the continuation Newton method(21) is favourable to adopt the trust-region technique to accurately follow the trajec-tory of the continuous Newton flow in the transient-state phase and to maintain itsfast convergence rate near the the equilibrium point z ∗ .We select the perturbed factor µ k + as the average of the residual sum [14,28,34]: µ k + = ( x k + ) T s k + n . (22)In equation (21), µ k + is approximated by σ k µ k , where the penalty coefficient σ k issimply selected as follows: σ k = (cid:40) . , when u k > . , u k , when u k ≤ . . (23)Thus, from equations (21)-(23), we have the following iteration scheme: A A T IS k X k ∆ x k ∆ y k ∆ s k = − Ax k − bA T y k + s k − cX k S k e − σ k µ k e = − F σ k µ k ( z k ) (24) erturbed Newton Method with Trust-region Time-stepping Schemes for Linear Programming 7 and ( x k + , y k + , s k + ) = ( x k , y k , s k ) + ∆ t k + ∆ t k ( ∆ x k , ∆ y k , ∆ s k ) , (25)where the perturbed function F µ ( z ) is defined by equation (9) and its Jacobian matrix J ( z ) is computed by equation (8).When the rank of matrix A is full, the linear system (24) can be solved by thefollowing three subsystems: AX k S − k A T ∆ y k = − r kp + AS − k (cid:16) − X k r kd + S k X k e − σ k µ k e (cid:17) , (26) ∆ s k = − r kd − A T ∆ y k , (27) ∆ x k = − S − k ( X k S k e + X k ∆ s k − σ k µ k e ) , (28)where the primal residual r kp and the dual residual r kd are defined by r kp = Ax k − b , (29) r kd = A T y k + s k − c . (30)Since the matrix AX k S − k A T is symmetric positive definite when the rank of matrix A is full and ( X k , S k ) (cid:31)
0, the linear system (26) can be solved efficiently by theCholesky factorization (see pp. 163-165, [17]).2.3 The trust-region time-stepping schemeAnother issue is how to adaptively adjust the time-stepping size ∆ t k at every iter-ation. A popular way to control the time-stepping size is based on the trust-regiontechnique [7, 10, 19, 21, 22, 24]. Its main idea is that the time-stepping size ∆ t k + will be enlarged when the linear model F σµ k ( z k ) + J ( z k )( z k + − z k ) approximates F σµ k ( z k + ) well, and ∆ t k + will be reduced when F σµ k ( z k ) + J ( z k )( z k + − z k ) ap-proximates F σµ k ( z k + ) badly. We adopt the following ratio ρ k as the measurementbetween the linear approximation model and the residual function: ρ k = (cid:107) F σ k µ k ( z k ) (cid:107) − (cid:107) F σ k µ k ( z k + ) (cid:107)(cid:107) F σ k µ k ( z k ) (cid:107) − (cid:107) F σ k µ k ( z k ) + J ( z k )( z k + − z k ) (cid:107) . (31)Thus, based on the measurement model (31), we give the following trust-region time-stepping scheme: ∆ t k + = ∆ t k , if 0 ≤ | − ρ k | ≤ η and ( x k + , s k + ) > , ∆ t k , else if η < | − ρ k | < η and ( x k + , s k + ) > , ∆ t k , others , (32)where the constants are selected as η = . , η = .
75, according to our numericalexperiments.
Luo, Yao A in problem (1) may be deficient andthe constraint system is even inconsistent when the right-hand-side vector b has smallnoise [23]. For the consistent system of constraints with the deficient rank of matrix A ,there are some pre-solving methods to eliminate the redundant constraints [2]. Here,in order to handle the inconsistent system of constraints, we consider the followingthe following best approximation problemmin x ∈ ℜ n (cid:107) Ax − b (cid:107) (33)and obtain the reduced system of constraints of problem (1).Firstly, we factorize the matrix A with its singular value decomposition (see pp.76-80, [17]) as follows: A = U Σ V T , Σ = Σ r
00 0 , Σ r = diag ( λ , λ , . . . , λ r ) (cid:31) , (34)where U ∈ ℜ m × m and V ∈ ℜ n × n are orthogonal matrices, and r is the rank of matrix A . Then, from equation (34), we know that problem (33) equals the following linearleast-squares problem min x ∈ ℜ n (cid:107) Σ V T x − U T b (cid:107) . (35)From equation (35), we obtain its solution V Tr x = b r , (36)where V r = V ( n , r ) , V Tr V r = I , and b r = (( U T b )( r )) ./ diag ( Σ ( r , r )) .Therefore, when the constraints of problem (1) is consistent, problem (1) equalsthe following linear programming problem:min x ∈ ℜ n c T x subject to V Tr x = b r , x ≥ . (37)When the constraints of problem (1) is inconsistent, problem (37) is the best relax-ation approximation of the original problem (1). Consequently, in subsection 2.2, itonly needs to replace matrix A and vector b with matrix V Tr and vector b r respectively,then the perturbed Newton method can handle the deficient rank problem.From the reduced linear programming problem (37), we obtain its KKT condi-tions as follows: V r x − b r = , (38) V Tr y r + s − c = , (39) XSe = , i = , , . . . , n , (40) ( x , s ) ≥ , (41) erturbed Newton Method with Trust-region Time-stepping Schemes for Linear Programming 9 where X = diag ( x ) and S = diag ( s ) . Thus, from the singular decomposition (33) ofmatrix A and equation (36), if ( x r , y r , s r ) is a solution of nonlinear equations (38)-(41), we know that (cid:0) x r , U (cid:0) Σ − r y r ; 0 (cid:1) , s r (cid:1) is a solution of nonlinear equations (2)-(5).Consequently, we recover the solution of the KKT conditions of the original problem(1) from the solution of the KKT conditions of the reduced problem (37).According to the above discussions, we give the detailed descriptions of the per-turbed Newton method with the trust-region time-stepping scheme for linear pro-gramming problems (1) in Algorithm 1. In order to simplify the convergence analysis of Algorithm 1, we assume that (i)the initial point ( x , s ) belongs to the strictly primal-dual feasible region F definedby equation (10), and (ii) the time-stepping size ∆ t k is selected such that ( x k , y k , s k ) satisfies the one-sided ∞ -norm neighborhood N ∞ ( γ ) , which is defined by N − ∞ ( γ ) = { ( x , y , s ) ∈ F | XSe ≥ γ µ e } , (42)where X = diag ( x ) , S = diag ( s ) , e = [ , , . . . , ] T , µ = x T s / n and γ is a small pos-itive constant such as γ = − . Without the loss of generality, we assume that thematrix A ∈ ℜ m × n has full rank. Lemma 1
Assume that the strictly primal-dual feasible point ( x k , y k , s k ) ∈ F sat-isfies the proximity condition (42) , then it exists a sufficiently small positive number α k > such that ( x k ( α k ) , y k ( α k ) , s k ( α k ) = ( x k , y k , s k ) + α k ( ∆ x k , ∆ y k , ∆ s k ) satisfies the proximity condition (42) , where ( ∆ x k , ∆ y k , ∆ s k ) is solved by the linearsystem (24) and α k = ∆ t k + ∆ t k . Proof.
Since ( x k , y k , s k ) is a strictly primal-dual feasible point, from equation (24),we know that A ∆ x k = , A T ∆ y k + ∆ s k = , k = , , . . . . (43)Thus, from equation (43), we have ∆ x Tk ∆ s k = − ∆ x Tk A T ∆ y k = ( A ∆ x k ) T ∆ y k = , k = , , . . . . (44)and Ax k ( α k ) = b , A T y k ( α k ) + s k ( α k ) = c , k = , , . . . . (45) Algorithm 1
Perturbed Newton Method with the trust-region time-stepping schemefor linear programming (The PNMTRLP method)
Input: matrix A ∈ ℜ m × n , vectors b ∈ ℜ m and c ∈ ℜ n for the linear programming problem:min x ∈ ℜ n c T x subject to Ax = b , x ≥ Output: the primal-dual optimal solution: ( solx , soly , sols ) .the maximum error of KKT conditions: KKT Error = max {(cid:107) A ∗ solx − b (cid:107) ∞ , (cid:107) A T ∗ soly + sols − c (cid:107) ∞ , (cid:107) solx . ∗ soly (cid:107) ∞ } .1: Initialize the trust-region parameters: η a = − , η = . , η = . , and set ε = − , ∆ t = . A with the singular value decomposition as follows: A = U Σ V T , Σ = diag ( λ , λ , ..., λ r , , ..., ) , and denote V r = V ( n , r ) , b r = (( U T b )( r )) ./ diag ( Σ ( r , r )) , where r is the rank of matrix A .3: Set bigM = max(max(abs( V Tr ))), bigM = max([norm( b r , inf), norm(c, inf), bigM]).4: Initialize x = bigMfac*bigM*ones(n,1), s = x , y = zeros( r , 1).5: flag success trialstep = 1; % Indicates whether time-stepping size is accepted6: while (itc < maxit) do if (flag success trialstep == 1) then
8: Set itc = itc + 1;9: Compute F k = V Tr x k − b r , F k = V r y k + s k − c , F k = X k s k ; µ k = x Tk s k / n ;10: Compute Resk inf = norm ([ F k ; F k ; F k ] , inf ) .11: if ( Resk inf < ε ) then
12: break;13: end if
14: Set σ k = min ( . , µ k ) .15: Compute F k = F k − σ k µ k ∗ ones ( n , ) , and set F k = [ F k ; F k ; F k ] .16: Use the Cholesky factorization to solve the system of linear equations: V Tr X k S − k V r ∆ y k = − (cid:0) F k + V Tr S − k (cid:0) X k F k − F k (cid:1)(cid:1) .Compute ∆ s k = − F k − V r ∆ y k and ∆ x k = − S − k (cid:0) F k + X k ∆ s k (cid:1) . end if
18: Compute ( x k + , y k + , s k + ) = ( x k , y k , s k ) + ∆ t k + ∆ t k ( ∆ x k , ∆ y k , ∆ s k ) .19: Compute F k + = V Tr x k + − b r , F k + = V r y k + + s k + − c , F k + = X k + s k + − σ k µ k ∗ ones ( n , ) , F k + = [ F k + ; F k + ; F k + ] .20: Compute LinAppF k + = [ F k + ; F k + ; ( F k + X k ( s k + − s k ) + S k ( x k + − x k ))] .21: Compute the ratio ρ k = ( (cid:107) F k (cid:107) − (cid:107) F k + (cid:107) ) / ( (cid:107) F k (cid:107) − (cid:107) LinAppF k + (cid:107) ) .22: if ( ( | ρ k − | ≤ η ) && ( x k + , s k + ) > ) ) then
23: Set ∆ t k + = ∆ t k ;24: else if ( ( η < | ρ k − | ≤ η ) && ( x k + , s k + ) > )) then
25: Set ∆ t k + = ∆ t k ;26: else
27: Set ∆ t k + = . ∆ t k ;28: end if if ( ( ρ k ≥ η a ) && ( x k + , s k + ) > ) ) then
30: Accept the trial step ( x k + , y k + , s k + ) ; flag success trialstep = 1;31: else
32: Set ( x k + , y k + , s k + ) = ( x k , y k , s k ) ; flag success trialstep = 0;33: end if
34: Set k ← k + end while
36: Return ( solx , soly , sols ) = ( x k , y k , s k ) ; KKT Error = max {(cid:107) Ax k − b (cid:107) ∞ , (cid:107) A T y k + s k − c (cid:107) ∞ , (cid:107) x k . ∗ s k (cid:107) ∞ } .erturbed Newton Method with Trust-region Time-stepping Schemes for Linear Programming 11 Since α k = ∆ t k / ( ∆ t k + ) , from equation (25), we have X k ( α k ) S k ( α k ) = ( X k + α k ∆ X k )( S k + α k ∆ S k )= X k S k + α k ( X k ∆ S k + S k ∆ X k ) + α k ∆ X k ∆ S k . (46)Replace X k ∆ S k + S k ∆ X k with equation (28) into equation (46), then we have X k ( α k ) S k ( α k ) e = X k S k e + α k ( σ k µ k − X k S k e ) + α k ∆ X k ∆ S k e = ( − α k ) X k S k e + α k σ k µ k + α k ∆ X k ∆ S k e . (47)From equation (44) and equation (47), we obtain µ ( α k ) = n e T X k ( α k ) S k ( α k ) e = ( − ( − σ k ) α k ) µ k . (48)We denote β kmax = max ≤ i ≤ n {| ∆ x ik | , | ∆ s ik |} . (49)Then, from equation (47), we have X k ( α k ) S k ( α k ) e ≥ ( − α k ) γ µ k + α k σ k µ k − α k β kmax . (50)From equation (48) and inequality (50), we know that the proximity condition X k ( α k ) S k ( α k ) e ≥ γ µ k ( α k ) is satisfied, provided that ( − α k ) γ µ k + α k σ k µ k − α k β kmax ≥ γ ( − α k + α k σ k ) µ k . Reformulate this expression, then we obtain α k ( − γ ) σ k µ k ≥ α k β kmax . (51)Therefore, if we select α k ≤ ( − γ ) σ k µ k β kmax , (52)inequality (51) is true. Consequently, it exists the sufficiently small positive number α k such that ( x k ( α k ) , y k ( α k ) , s k ( α k )) satisfies the proximity condition (42). (cid:117)(cid:116) In the following Lemma 2, we give the lower bounded estimation of ( x k , s k ) whenthe initial point ( x , y , s ) is strictly primal-dual feasible and satisfies the proximitycondition (42). Lemma 2
Assume that the initial point ( x , y , s ) ∈ F is strictly primal-dual fea-sible, and ( x k , y k , s k ) ( k = , , . . . ) generated by Algorithm 1 satisfy the proximitycondition (42) . Furthermore, if it exists a constant C µ such that µ k ≥ C µ > , k = , , . . . , (53) we have the bounded estimations of ( x k , s k ) as follows: < C min ≤ min ≤ i ≤ n { x ik , s ik } ≤ max ≤ i ≤ n { x ik , s ik } ≤ C max , k = , , . . . , (54) where C min and C max are two positive constants. Proof.
Since ( x k , y k , s k ) is generated by Algorithm 1, from equation (48), we have µ k + = x Tk + s k + n = ( − ( − σ k ) α k ) µ k ≤ µ k , k = , , . . . , which gives µ k + = k ∏ i = ( − ( − σ i ) α i ) µ i ≤ µ , k = , , . . . . (55)Furthermore, since the initial point ( x , s , y ) is strictly primal-dual feasible,from Lemma 1 and the mathematical induction, we know that ( x k , y k , s k ) is strictlyprimal-dual feasible. Thus, from equation (45), we have A ( x k − x ) = , and A T ( y k − y ) + ( s k − s ) = , which gives ( x k − x ) T ( s k − s ) = . Rearrange this expression, from equation (55), then we obtain x Tk s + s Tk x = x k s k + x s ≤ n ( µ k + µ ) ≤ n µ , which gives x ik ≤ n µ min ≤ j ≤ n { s j } and s ik ≤ n µ min ≤ j ≤ n { x j } , ≤ i ≤ n , k = , , . . . . Therefore, if we selected C max = max (cid:40) n µ min ≤ j ≤ n { s j } , n µ min ≤ j ≤ n { x j } (cid:41) , we obtain max ≤ i ≤ n { x ik , s ik } ≤ C max , k = , , . . . . (56) erturbed Newton Method with Trust-region Time-stepping Schemes for Linear Programming 13 On the other hand, from the assumption (53) and the proximity condition (42),we have x ik s ik ≥ γ µ k ≥ γ C µ , ≤ i ≤ n , k = , , . . . . Combining the upper estimation (56) of ( x k , s k ) , we obtain x ik ≥ γ C µ max ≤ j ≤ n { s jk } ≥ γ C µ C max , and s ik ≥ γ C µ max ≤ j ≤ n { x jk } ≥ γ C µ C max , k = , , . . . , which concludes that min ≤ i ≤ n { x ik , s ik } ≥ C min , k = , , . . . , if we select C min = γ C µ / C max . (cid:117)(cid:116) Lemma 3
Assume that the initial point ( x , y , s ) ∈ F is strictly primal-dual fea-sible, and ( x k , y k , s k ) ( k = , , . . . ) generated by Algorithm 1 satisfy the proximitycondition (42) . Furthermore, if the assumption (53) is true, it exits positive constantsC ∆ x and C ∆ s such that (cid:107) ∆ s k (cid:107) ≤ C ∆ s and (cid:107) ∆ x k (cid:107) ≤ C ∆ x , k = , , . . . . (57) Proof.
Factorize matrix A with the singular value decomposition (34). Then, fromthe bounded estimation (54) of ( x k , s k ) in Lemma 2, we have vAXS − A T v ≥ C min C max (cid:107) Av (cid:107) ≥ C min λ min C max (cid:107) v (cid:107) , ∀ v ∈ ℜ n , (58)and vAXS − A T v ≤ C max C min (cid:107) Av (cid:107) ≤ C max λ max C min (cid:107) v (cid:107) , ∀ v ∈ ℜ n , (59)where λ min and λ max are the smallest and largest singular values of matrix A , respec-tively. Since ( x , y , s ) is strictly primal-dual feasible, from equations (24)-(25), Al-gorithm 1, and the mathematical induction, we know that ( x k , y k , s k ) ( k = , , . . . ) arestrictly primal-dual feasible. Consequently, from equations (26), (54) and (58)-(59),we obtain C min λ min C max (cid:107) ∆ y k (cid:107) ≤ ∆ y Tk (cid:0) AX k S − k A T (cid:1) ∆ y k = ∆ y Tk AS − k ( X k S k e − σ k µ k e ) ≤ (cid:107) ∆ y k (cid:107)(cid:107) A (cid:107) (cid:13)(cid:13) S − k (cid:13)(cid:13) (cid:107) X k S k e − σ k µ k e (cid:107) ≤ (cid:107) ∆ y k (cid:107) λ max C min ( (cid:107) X k S k e (cid:107) + n σ k µ k ) ≤ (cid:107) ∆ y k (cid:107) λ max C min ( (cid:107) X k S k e (cid:107) + n σ k µ k ) = (cid:107) ∆ y k (cid:107) λ max C min ( + σ k ) n µ k , which gives (cid:107) ∆ y k (cid:107) ≤ C max λ max C min λ min n µ k ≤ C max λ max C min λ min n µ (60) where we use the properties σ k ≤ µ k ≤ µ from inequality(55). Therefore, from equation (27) and inequality (60), we have (cid:107) ∆ s k (cid:107) = (cid:107) − A T ∆ y k (cid:107) ≤ (cid:107) A T (cid:107)(cid:107) ∆ y k (cid:107) ≤ C max λ max C min λ min n µ (cid:44) C ∆ s . Thus, we prove the first part of inequality (57).From equation (28), inequality (54) and the first part of inequality (57), we have (cid:107) ∆ x k (cid:107) = (cid:13)(cid:13) − S − k ( X k S k e + X k ∆ s k − σ µ k e ) (cid:13)(cid:13) ≤ (cid:13)(cid:13) S − k (cid:13)(cid:13) (cid:107) X k S k e + X k ∆ s k − σ k µ k e (cid:107)≤ C min ( (cid:107) X k S k e (cid:107) + (cid:107) X k ∆ s k (cid:107) + (cid:107) σ k µ k e (cid:107) ) ≤ C min ( (cid:107) X k S k e (cid:107) + (cid:107) X k (cid:107)(cid:107) ∆ s k (cid:107) + n σ k µ k (cid:107) ) ≤ C min ( n µ k + C max C ∆ s + n σ k µ k ) ≤ C min ( n µ + C max C ∆ s ) (cid:44) C ∆ x , where we use the properties σ k ≤ µ k ≤ µ from inequality(55). Thus, we also prove the second part of inequality (57). (cid:117)(cid:116) Lemma 4
Assume that the initial point ( x , y , s ) ∈ F is strictly primal-dual fea-sible, and ( x k , y k , s k ) ( k = , , . . . ) generated by Algorithm 1 satisfy the proximitycondition (42) . Furthermore, if the assumption (53) is true, it exists a constant C ∆ t such that ∆ t k ≥ C ∆ t > , k = , , . . . . (61) Proof.
From equation (31) and equations (24)-(25), we have | ρ k − | = (cid:12)(cid:12)(cid:12)(cid:12) (cid:107) F σ k µ k ( z k ) (cid:107) − (cid:107) F σ k µ k ( z k + ) (cid:107)(cid:107) F σ k µ k ( z k ) (cid:107) − (cid:107) F σ k µ k ( z k ) + J ( z k )( z k + − z k ) (cid:107) − (cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12) (cid:107) F σ k µ k ( z k ) (cid:107) − (cid:107) F σ k µ k ( z k + ) (cid:107)(cid:107) F σ k µ k ( z k ) (cid:107) − (cid:107) F σ k µ k ( z k ) − ( ∆ t k ) / ( + ∆ t k ) F σ k µ k ( z k ) (cid:107) − (cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12) (cid:107) F σ k µ k ( z k ) (cid:107) − ( + ∆ t k ) (cid:107) F σ k µ k ( z k + ) (cid:107) (cid:12)(cid:12) ∆ t k (cid:107) F σ k µ k ( z k ) (cid:107)≤ (cid:13)(cid:13) ( + ∆ t k ) (cid:0) F σ k µ k ( z k + ) − F σ k µ k ( z k ) (cid:1) + ∆ t k F σ k µ k ( z k ) (cid:13)(cid:13) ∆ t k (cid:107) F σ k µ k ( z k ) (cid:107) = ∆ t k + ∆ t k (cid:107) ∆ X k ∆ S k e (cid:107)(cid:107) X k s k − σ k µ k e (cid:107) . (62)In the last equality of equation (62), we use the property F σ k µ k ( z k + ) − F σ k µ k ( z k ) = F σ k µ k (cid:18) z k + ∆ t k + ∆ t k ∆ z k (cid:19) − F σ k µ k ( z k )= ∆ t k + ∆ t k J ( z k ) ∆ z k + (cid:18) ∆ t k + ∆ t k (cid:19) ∆ X k ∆ S k e = − ∆ t k + ∆ t k F σ k µ k ( z k ) + (cid:18) ∆ t k + ∆ t k (cid:19) ∆ X k ∆ S k e . erturbed Newton Method with Trust-region Time-stepping Schemes for Linear Programming 15 On the other hand, using the property (cid:107) a (cid:107) ≥ a i ( i = , , . . . , n ) , we have (cid:107) X k s k − σ k µ k e (cid:107) ≥ x ik s ik − σ k µ k , i = , , . . . , n . By summing the n components of the above two sides and the definition µ k = x Tk s k / n ,we obtain (cid:107) X k s k − σ k µ k e (cid:107) ≥ ( − σ k ) µ k ≥ ( − . ) C µ = . C µ , (63)where we use the properties σ k ≤ . µ k ≥ C µ from theassumption (53).Thus, from the bounded estimation (57) of ( ∆ x k , ∆ s k ) and inequalities (62), (63),we obtain | ρ k − | ≤ ∆ t k + ∆ t k (cid:107) ∆ X k ∆ S k e (cid:107) . C µ ≤ ∆ t k + ∆ t k (cid:107) ∆ X k (cid:107)(cid:107) ∆ s k (cid:107) . C µ ≤ ∆ t k + ∆ t k (cid:107) ∆ x k (cid:107)(cid:107) ∆ s k (cid:107) . C µ ≤ ∆ t k + ∆ t k C ∆ x C ∆ s . C µ ≤ η , provided that ∆ t k ≤ . C µ η C ∆ x C ∆ s . (64)Thus, if we assume that K is the first index such that ∆ t K satisfies inequality (64), ac-cording to the trust-region time-stepping scheme (32), ∆ t K + will be enlarged. There-fore, we prove the result (61) if we select C ∆ t as ∆ t K . (cid:117)(cid:116) According to the above discussions, we give the global convergence analysis ofAlgorithm 1.
Theorem 1
Assume that the initial point ( x , y , s ) ∈ F is strictly primal-dual fea-sible, and ( x k , y k , s k ) ( k = , , . . . ) generated by Algorithm 1 satisfy the proximitycondition (42) . Then, we have lim k → ∞ µ k = , (65) where µ k = x Tk s k / n. Proof.
Assume that it exists a positive constant C µ such that µ k ≥ C µ , k = , , . . . . (66)Then, according to the result of Lemma 4, we know that it exists a positive constant C ∆ t such that ∆ t k ≥ C ∆ t , k = , , . . . . Therefore, from equation (48), we have µ k + = µ k ( α k ) = ( − ( − σ k ) α k ) µ k ≤ (cid:18) − ( − σ k ) C ∆ t + C ∆ t (cid:19) µ k ≤ (cid:18) − . C ∆ t + C ∆ t (cid:19) µ k ≤ (cid:18) − . C ∆ t + C ∆ t (cid:19) k + µ , where we use the property σ k ≤ . ( k = , , . . . ) from the definition (23) of σ k ( k = , , . . . ) . Thus, we have µ k →
0, which contradicts the assumption (66). Conse-quently, we obtain lim k → ∞ inf µ k =
0. Since µ k is monotonically decreasing, it is notdifficult to know lim k → ∞ µ k =
0. Furthermore, we obtain lim k → ∞ (cid:107) X k s k (cid:107) = (cid:107) X k s k (cid:107) ≤ (cid:107) X k s k (cid:107) = n µ k and ( x k , s k ) > (cid:117)(cid:116) In this section, we test Algorithm 1 (the PNMTRLP method) for some linear pro-gramming problems with full rank matrix A or deficient rank matrix A , and compareit with the excellent path-following method (pathfollow.m) coded by M. C. Ferrisin [14] and the state-of-the-art Mehrotra’s predictor-corrector algorithm (the built-insubroutine linprog.m of the MATLAB environment, which was implemented by S.mehrotra and Y. Zhang [26, 27, 36]).The tolerance of three methods is set to ε = − by default. We use the maximumabsolute error (KKTError) of the KKT conditions (2)-(5) and the primal-dual gap x T s to measure the error between the numerical optimal solution and the theoreticaloptimal solution.4.1 The problem with the full rank matrixFor the standard linear programming problem with the full rank matrix A , the sparsematrix A of given density 0.2 is randomly generated and we choose feasible x , y , s at random, with x and s each about half-full. The dimension of matrix A varies from10 ×
100 to 300 × KKTError are small. However, the subrou-tine pathfollow.m cannot perform well in some higher-dimensional problems, suchas examples 7 , , , , , , , ,
30, since their solutions do not satisfy the erturbed Newton Method with Trust-region Time-stepping Schemes for Linear Programming 17
Algorithm 2
The standard linear programming problem with the full rank matrix
Input: the number of equality constraints: m ;the number of unknown variables: n . Output: matrix A and vectors b ∈ ℜ m and c ∈ ℜ n .1: density=0.2;2: A = sprandn(m,n,density); % Generate a sparse matrix of give density.3: xfeas = [rand(n/2,1); zeros(n-(n/2),1)];4: sfeas = [zeros(n/2,1); rand(n-(n/2),1)];5: xfeas = xfeas(randperm(n));6: sfeas = sfeas(randperm(n));7: yfeas = (rand(m,1)-0.5)*4;% Choose b and c to make this (x,y,s) feasible.8: b = A*xfeas;9: c = A’*yfeas + sfeas; KKT conditions. From Figure 1, we also find that the subroutine linprog.m performsthe best and the number of its iterations is less than 20, the number of iterations ofthe PNMTRLP method is around 20, and the number of iterations of the subroutinepathfollow.m often reaches the maximum number (i.e. 200 iterations). Therefore, thePNMTRLP method is also an efficient and robust path-following method for linearprogramming problems with full rank matrices.
Exam It e r a t i on s PN M TRLP l inprogpath-following (a) The number of iterations Exam S o l u t i on t i m e PN M TRLP l inprogpath-following (b) The consumed time Fig. 1: The number of iterations and the consumed time.4.2 The problem with the deficient rank matrix and uncertain dataFor some practical problems, since it may exist redundant constraints and measure-ment errors, the rank of matrix A may be deficient and the right-hand-side vec- Table 1: Numerical results of problems with full rank matrices.
Problem ( m × n , r ) PNMTRLP linprog pathfollowKKTError Gap KKTError Gap KKTError GapExam. 1 (10 × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × tor b may have some noise, which may result in the inconsistent system of con-straints. In order to verify the effect of the PNMTRLP method handling those prob-lems, we randomly generate the sparse matrix A and let A ( m − i , : ) = i ∗ A ( m , : ) ( i = , , . . . , m − r ) , where r is the rank of matrix A , then we compare it with the subrou-tine linprog.m for those problems. The dimension of matrix A varies from 10 × × A , we only report the numerical results of the PNMTRLPmethod and the subroutine linprog.m, and put them in Tables 3-8. From Tables 3-8, we can find that the PNMTRLP method can solve all those problems. However, erturbed Newton Method with Trust-region Time-stepping Schemes for Linear Programming 19 Algorithm 3
The standard linear programming problem with the deficient rank ma-trix and the noise right-hand-side vector
Input: the number of equality constraints: m ;the number of unknown variables: n ;the rank of the matrix A : r ;the maximum noise of the right-hand-side vector b : eps. Output: matrix A and vectors b ∈ ℜ m and c ∈ ℜ n .1: density = 0.2;2: A = sprandn(m,n,density); % Generate a sparse matrix of give density.3: for i = , , ..., r do
4: A(m-i, :) = i*A(m, :);5: end for
6: xfeas = [rand(n/2,1); zeros(n-(n/2),1)];7: sfeas = [zeros(n/2,1); rand(n-(n/2),1)];8: xfeas = xfeas(randperm(n));9: sfeas = sfeas(randperm(n));10: yfeas = (rand(m,1) - 0.5)*4;% Choose b and c to make this (x, y, s) feasible.11: b = A*xfeas + 2*(rand(m,1)-0.5)*eps;12: c = A’*yfeas + sfeas;
Table 2: Statistical results of deficient rank problems with the noise data.
PNMTRLP linprognumber of failed problems 0/180 143/180 the subroutine linprog.m can not solve some problems, because those problems donot have feasible solution and the subroutine linprog.m outputs
NaN . Furthermore,although the subroutine linprog.m can obtain the s olutions for some problems, suchas examples A1, A2, A5, A6, etc, the KKT errors and the primal-dual gaps are large,so we conclude that the subroutine linprog.m also fails to solve those problems. Table2 is the statistical results of failed problems solved by the PNMTRLP method andthe subroutine linprog.m. From 2, we can find that the subroutine linprog.m can noteffectively handle the linear programming problem with the deficient rank of matrix A and the noise right-hand-side vector b , because the failed problems of the subroutinelinprog.m are around 80%.We also give the trends of the KKT errors and the primal-dual gaps for the dif-ferent right-hand-side noise from 10 − to 10 − in the right sub-figure and the leftsub-figure in Figure 2, respectively. From Figure 2, we find that the KKT errors andthe primal-dual gaps of the PNMTRLP method are less than 0 .
12 when the right-hand-side noise comes up to 0 .
1. Therefore, the PNMTRLP method can be used inpractical application scenarios, which will greatly reduce the cost of industrial hard-ware facilities.
Exam KK T E rr o r noise = 1e-1noise = 1e-2noise = 1e-3noise = 1e-4noise = 1e-5noise = 1e-6 (a) The maximum KKT error Exam G ap noise = 1e-1noise = 1e-2noise = 1e-3noise = 1e-4noise = 1e-5noise = 1e-6 (b) The primal-dual gap Fig. 2: The KKT error and the primal-dual gap of the PNMTRLP method.
For the standard linear programming problem with the deficient rank of matrix A and the noise right-hand-side vector b , we give a preprocessing method based on thesingular value decomposition method. Then, according to the first order KKT condi-tions of the linear programming, we convert the treated problems into the equivalentproblem of nonlinear equations with nonnegative constraints. Based on the system ofnonlinear equations with nonnegative constraints, we consider a special continuousNewton flow with nonnegative constraints, which has the nonnegative steady-state so-lution for any nonnegative initial point. In order to improve the robust and efficiencyof tracing the trajectory, we construct a perturbed Newton method with the new time-stepping scheme based on the trust-region technique (the PNMTRLP method, Algo-rithm 1). We also prove the global convergence of the PNMTRLP method when theinitial point is strictly primal-dual feasible. Finally, we test some standard linear pro-gramming problems with full rank matrices or deficient rank matrices, and compareit with the path-following method (the subroutine pathfollow.m [14]) and the state-of-the-art Mehrotra’s predictor-corrector algorithm (the built-in subroutine linprog.mof the MATLAB environment [26, 27, 36]). The numerical results show that the PN-MTRLP method performs well and the subroutines pathfollow.m and linprog.m cannot effectively handle the problem with the deficient rank of matrix A and the noiseright-hand-side vector b . Therefore, the PNMTRLP method is worth investigatingfurther as a new path-following method. Acknowledgments
This work was supported in part by Grant 61876199 from National Natural ScienceFoundation of China, Grant YBWL2011085 from Huawei Technologies Co., Ltd., erturbed Newton Method with Trust-region Time-stepping Schemes for Linear Programming 21 and Grant YJCB2011003HI from the Innovation Research Program of Huawei Tech-nologies Co., Ltd..
References
1. E. L. Allgower and K. Georg,
Introduction to Numerical Continuation Methods , SIAM, Philadelphia,PA, 2003.2. E. D. Anderson and K. D. Anderson,
Presolving in linear programming , Mathematical Programming, : 221-245, 1995.3. U. M. Ascher and L. R. Petzold, Computer Methods for Ordinary Differential Equations andDifferential-Algebraic Equations , SIAM, Philadelphia, PA, 1998.4. O. Axelsson and S. Sysala,
Continuation Newton methods , Computer and Mathematics with Applica-tions, : 2621-2637, 2015.5. F. H. Branin, Widely convergent method for finding multiple solutions of simultaneous nonlinear equa-tions , IBM Journal of Research and Development, : 504-521 (1972).6. K. E. Brenan, S. L. Campbell and L. R. Petzold, Numerical solution of initial-value problems indifferential-algebraic equations , SIAM, Philadelphia, PA, 1996.7. A. R. Conn, N. Gould and Ph. L. Toint,
Trust-Region Methods , SIAM, 2000.8. , B. D. Craven and S. M. N. Islam,
Linear programming with Uncertain data: extensions to robustoptimization , Journal of Optimization Theory and Applications, : 673-679, 2012.9. D. F. Davidenko,
On a new method of numerical solution of systems of nonlinear equations (in Rus-sian), Dokl. Akad. Nauk SSSR, : 601-602, 1953.10. P. Deuflhard, Newton Methods for Nonlinear Problems: Affine Invariance and Adaptive Algorithms ,Springer-Verlag, 2004.11. P. Deuflhard, H. J. Pesch and P. Rentrop
A modified continuation method for the numerical solution ofnonlinear two-point boundary value problems by shooting techniques , Numerische Mathematik, :327-343, 1975.12. J. E. Dennis and R. B. Schnabel, Numerical Methods for Unconstrained Optimization and NonlinearEquations , Science Press, 2009.13. E. J. Doedel,
Lecture notes in numerical analysis of nonlinear equations , in B. Krauskopf, H. M.Osinga and J. Gal ´ a n-Vioque (Eds.), Numerical Continuation Methods for Dynamical Systems, 1-50,Springer, 2007.14. M. C. Ferris, O. L. Mangasarian and S. J. Wright, Linear Programming with MATLAB , SIAM,Philadelphia, PA, http://pages.cs.wisc.edu/~ferris/ferris.books.html , 2007.15. A. V. Fiacco and G. P. McCormick,
Nonlinear programming: Sequential Unconstrained MinimizationTechniques , SIAM, 1990.16. C. C. Gonzaga,
Path-following methods for linear programming , SIAM Reivew, (2): 167-224,1992.17. G. H. Golub and C. F. Van Loan, Matrix Computation, fourth editon, The John Hopkins UniversityPress, 2013.18. E. Hairer and G. Wanner, Solving Ordinary Differential Equations II, Stiff and Differential-AlgebraicProblems , 2nd ed., Springer-Verlag, Berlin, 1996.19. D. J. Higham,
Trust region algorithms and timestep selection , SIAM Journal of Numerical Analysis, : 194-210 (1999).20. N. Karmarkar, A new polynomial-time algorithm for linear programming , Combinatorica, ( A second-order pseudo-transient method for steady-state problems , Applied Mathematicsand Computation, : 1752-1762 (2010).22. X.-L. Luo,
A dynamical method of DAEs for the smallest eigenvalue problem , Journal of Computa-tional Science, : 113-119 (2012).23. X.-L. Luo, J.-H. Lv and G. Sun, Continuation method with the trust-region time-stepping schemefor linearly constrained optimization with noise data , http://arxiv.org/abs/2005.05965 , May2020.24. X.-L. Luo, H. Xiao and J.-H. Lv, Continuous Newton method with the trust-region time-steppingscheme , https://arxiv.org/abs/2006.02634 , June 2020.2 Luo, Yao25. J. M. Ortega and W. C. Rheinboldt, Iteration Solution of Nonlinear Equations in Several Variables ,Classics in Applied Mathematics, SIAM, 2000.26. MATLAB 9.6.0 (R2019a), The MathWorks Inc., , 2019.27. S. Mehrotra,
On the implementation of a primal-dual interior point method , SIAM Journal on Opti-mization, : 575-601 (1992).28. J. Nocedal and S. J. Wright, Numerical Optimization , Springer-Verlag, 1999.29. M. J. D. Powell,
Convergence properties of a class of minimization algorithms , in: O.L. Mangasarian,R. R. Meyer and S. M. Robinson, eds.,
Nonlinear Programming 2 , Academic Press, New York, 1975,pp. 1-27.30. L. S. Pontryagin,
The Ordinary differential equations , Higher Education Press (in Chinese Translationfrom Russian), 2006.31. L. F. Shampine, I. Gladwell and S. Thompson,
Solving ODEs with MATLAB , Cambridge UniversityPress, 2003.32. K. Tanabe,
Continuous Newton-Raphson method for solving an underdetermined system of nonlinearequations , Nonlinear Analysis, : 495-503, 1979.33. K. Tanabe, Centered Newton method for mathematical programming , in System Modeling and Opti-mization: Proceedings of the 13th IFIP conference, vol. 113 of Lecture Notes in Control and Informa-tion Systems, Berlin, Springer-Verlag, pp. 197-206, 1988.34. S. J. Wright,
Primal-dual Interior Point Methods , SIAM, Philadelphia, 1997.35. Z.-D. Yuan, J.-G. Fei and D.-G. Liu,
The Numerical Solutions for Initial Value Problems of StiffOrdinary Differential Equations (in Chinese), Science Press, Beijing, 1987.36. Y. Zhang,
Solving large-Scale linear programs by interior-point methods under the MATLAB envi-ronment , Optimization Methods and Software, (1): 1-31, 1998, https://doi.org/10.1080/10556789808805699 .erturbed Newton Method with Trust-region Time-stepping Schemes for Linear Programming 23 A Numerical results of test problems
Table 3: Numerical results of deficient rank problems with noise ε b = − . Problem ( m × n , r ) PNMTRLP linprogCPU (s) KKTError Gap CPU (s) KKTError GapExam. A1 (10 × ,
4) 0.0781 0.0793 1.14E-05 0.1406 15.0250 463.69Exam. A2 (20 × ,
17) 0.063 0.0820 4.78E-05 0.1562 168.7704 17754.71Exam. A3 (30 × ,
23) 0.072 0.0932 0.0006 0.1094 Failed FailedExam. A4 (40 × ,
31) 0.0625 0.0911 0.0002 0.0156 Failed FailedExam. A5 (50 × ,
48) 0.0625 0.0448 0.0005 0.1250 0.9009 159.74Exam. A6 (60 × ,
51) 0.2031 0.0926 0.0008 0.1094 105.0479 24875.35Exam. A7 (70 × ,
68) 0.1718 0.0879 0.0032 0.1406 0.1266 24.59Exam. A8 (80 × ,
79) 0.2343 0.0032 0.0011 0.1250 0.01242 1.11Exam. A9 (90 × ,
89) 0.2031 0.0756 0.0010 0.1094 1.3138 385.91Exam. A10 (100 × ,
97) 0.2031 0.1067 0.0020 0.0781 125.0551 50816.85Exam. A11 (110 × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , Table 4: Numerical results of deficient rank problems with noise ε b = − . Problem ( m × n , r ) PNMTRLP linprogCPU (s) KKTError Gap CPU (s) KKTError GapExam. B1 (10 × ,
5) 0.0469 0.0078 0.0001 0.2344 48.9117 1771.6522Exam. B2 (20 × ,
19) 0.0521 0.0069 0.0002 0.0781 0.0138 0.5149Exam. B3 (30 × ,
21) 0.0625 0.0089 0.0003 0.1094 Failed FailedExam. B4 (40 × ,
35) 0.0625 0.0090 0.0003 0.1094 Failed FailedExam. B5 (50 × ,
47) 0.0625 0.0072 0.0014 0.1094 0.0322 0.3652Exam. B6 (60 × ,
54) 0.0781 0.0099 0.0006 0.0469 0.2956 74.6274Exam. B7 (70 × ,
66) 0.0156 0.0066 0.0005 0.0313 Failed FailedExam. B8 (80 × ,
76) 0.1875 0.0092 0.0013 0.0781 Failed FailedExam. B9 (90 × ,
89) 0.4219 0.0042 0.0012 0.1719 0.0062 0.2883Exam. B10 (100 × ,
98) 0.2813 0.0048 0.0021 0.1094 Failed FailedExam. B11 (110 × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , erturbed Newton Method with Trust-region Time-stepping Schemes for Linear Programming 25 Table 5: Numerical results of deficient rank problems with noise ε b = − . Problem ( m × n , r ) PNMTRLP linprogCPU (s) KKTError Gap CPU (s) KKTError GapExam. C1 (10 × ,
9) 0.0938 0.0005 0.0001 0.2031 22.7266 887.6482Exam. C2 (20 × ,
16) 0.0156 0.0008 0.0001 0.0156 Failed FailedExam. C3 (30 × ,
25) 0.0938 0.0007 0.0001 0.1406 Failed FailedExam. C4 (40 × ,
36) 0.0851 0.0010 0.0003 0.0781 Failed FailedExam. C5 (50 × ,
49) 0.0781 0.0007 0.0007 0.1094 0.0013 0.0005Exam. C6 (60 × ,
59) 0.1406 0.0001 0.0008 0.0938 0.0001 0.0028Exam. C7 (70 × ,
66) 0.0469 0.0008 0.0004 0.1250 Failed FailedExam. C8 (80 × ,
76) 0.1250 0.0006 0.0028 0.0469 Failed FailedExam. C9 (90 × ,
84) 0.4375 0.0007 0.0018 0.0469 Failed FailedExam. C10 (100 × ,
93) 0.2188 0.0008 0.0014 0.0938 Failed FailedExam. C11 (110 × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , Table 6: Numerical results of deficient rank problems with noise ε b = − . Problem ( m × n , r ) PNMTRLP linprogCPU (s) KKTError Gap CPU (s) KKTError GapExam. D1 (10 × ,
6) 0.0469 1.03E-4 0.0001 0.2500 0.0003 1.01E-4Exam. D2 (20 × ,
15) 0.0625 9.36E-5 0.0001 0.1094 Failed FailedExam. D3 (30 × ,
27) 0.1250 8.30E-5 0.0001 0.2656 Failed FailedExam. D4 (40 × ,
38) 0.1563 9.35E-5 0.0004 0.1094 0.0001 0.0005Exam. D5 (50 × ,
49) 0.0313 5.01E-5 0.0004 0.0938 0.3173 53.3101Exam. D6 (60 × ,
53) 0.3281 8.96E-5 0.0008 0.0625 Failed FailedExam. D7 (70 × ,
69) 0.2813 2.92E-5 0.0009 0.0625 128.9256 44395.4123Exam. D8 (80 × ,
79) 0.2813 6.33E-5 0.0020 0.1406 112.1089 40834.5967Exam. D9 (90 × ,
82) 0.0781 9.58E-5 0.0021 0.1094 Failed FailedExam. D10 (100 × ,
95) 0.1250 6.91E-5 0.0032 0.2656 Failed FailedExam. D11 (110 × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , erturbed Newton Method with Trust-region Time-stepping Schemes for Linear Programming 27 Table 7: Numerical results of deficient rank problems with noise ε b = − . Problem ( m × n , r ) PNMTRLP linprogCPU (s) KKTError Gap CPU (s) KKTError GapExam. E1 (10 × ,
5) 0.1250 9.63E-06 2.27E-05 0.2188 Failed FailedExam. E2 (20 × ,
19) 0.0781 3.62E-06 3.26E-05 0.1719 62.0437 5499.26Exam. E3 (30 × ,
28) 0.1719 7.06E-06 0.000299 0.1250 Failed FailedExam. E4 (40 × ,
35) 0.0625 8.93E-06 0.000166 0.0625 Failed FailedExam. E5 (50 × ,
42) 0.1406 9.79E-06 0.000243 0.1094 0.0008 0.0020Exam. E6 (60 × ,
58) 0.2344 1.01E-05 0.000554 0.1250 Failed FailedExam. E7 (70 × ,
69) 0.1875 4.52E-06 0.000787 0.1406 5.8E-06 4.7E-08Exam. E8 (80 × ,
78) 0.2656 1.07E-05 0.001612 0.0469 Failed FailedExam. E9 (90 × ,
82) 0.2813 5.48E-06 0.00198 0.2188 9.97E-06 0.0002Exam. E10 (100 × ,
95) 0.2188 9.59E-06 0.001049 0.1094 Failed FailedExam. E11 (110 × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , Table 8: Numerical results of deficient rank problems with noise ε b = − . Problem ( m × n , r ) PNMTRLP linprogCPU (s) KKTError Gap CPU (s) KKTError GapExam. F1 (10 × ,
2) 0.1094 9.1042E-07 1.95018E-05 0.1406 1.93E+21 701878381Exam. F2 (20 × ,
16) 0.0469 2.4042E-06 0.000196712 0.2500 Failed FailedExam. F3 (30 × ,
22) 0.0469 1.2484E-06 0.000182235 0.2969 7.96E-06 1.5039E-26Exam. F4 (40 × ,
34) 0.1094 3.2483E-06 0.000508886 0.0781 Failed FailedExam. F5 (50 × ,
48) 0.0625 1.1277E-06 0.000227713 0.0938 Failed FailedExam. F6 (60 × ,
58) 0.2031 2.8283E-06 0.000773048 0.1563 8.89E-07 6.8731E-25Exam. F7 (70 × ,
67) 0.2344 1.5834E-05 0.005080376 0.3594 Failed FailedExam. F8 (80 × ,
79) 0.1875 3.4235E-06 0.001094181 0.1563 1.65E-06 5.9364E-08Exam. F9 (90 × ,
82) 0.1250 3.9215E-06 0.001544717 0.1094 3.19E-06 6.4124E-24Exam. F10 (100 × ,
91) 0.2188 1.0383E-05 0.004762294 0.1250 Failed FailedExam. F11 (110 × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × , × ,,