Parallel Newton Methods for Numerical Analysis of Infeasible Linear Programming Problems with Block-Angular Structure of Constraints

Abstract

For the linear programming (LP) problems, perhaps infeasible, with block-angular matrices of constraints, two parallel second order optimization methods are proposed. Being applied to initial LP problem they give its solution if this problem is solvable, and automatically deliver a solution of some relaxed LP problem, otherwise. The methods use penalty and barrier functions as well as Tikhonov regularization of Lagrangian of initial problem. The methods contain a parameter which tends to zero and minimize a residual of constraints. Parallel calculations of matrix operations follow MPI-type paradigm.

Keywords

• linear programming
• infeasibility
• ill-posed problems
• generalized solution
• regularization
• penalty function
• parallel calculations

Introduction

Infeasible linear programs is widely studied class of the improper problems in mathematical programming theory (see [1,2,3,4,5,6] and others). They arise naturally in many applications of linear programming (LP) such as network ows analysis problems, problems of portfolio selection and scheduling problems which re ect complex production processes, e.g. in manufacture and power systems. Their inherent infeasibility may be caused by multiple factors such as errors in the input data, imbalance of its resource constraints and objectives, modelling bias as well as its structural distortions.

Basic approach addressing this issue is to introduce a more general form of solution for unsolvable problem which is in fact an optimum for some repaired or relaxed LP problem [1,2,7,8,9]. Following this approach in [10,11] two numerical methods were proposed which automatically adjust rigth-hand-side vector of initial LP problem if it is infeasible and then nd its generalized solution using penalty and barrier functions as well as Tikhonov regularization of standard Lagrangian [12,13,14]. The methods contain a small parameter which tends to zero and minimize the residual of constraints.

In this paper we discuss how these methods may be t to parallel calculations in MPI-type paradigm under the assumption that constraints matrix of initial LP problem has a block-angular structure.

Copyright ⃝c by the paper's authors. Copying permitted for private and academic purposes. In: A. Kononov et al. (eds.): DOOR 2016, Vladivostok, Russia, published at http://ceur-ws.org

Consider the linear program and its dual one

The Problem State and Preliminaries min f(c; x) : Ax = b; x

0g max f(b; y) : A⊤y

cg ; where A 2 IRm0 n0 , c 2 IRn0 and b 2 IRm0 are given, x 2 IRn0 and y 2 IRm0 are vectors of primal and dual variables respectively, rankA = m0 n0, ( ; ) denotes scalar product.

Assume that dual problem ( 2 ) is feasible, but it isn't known a priory whether primal problem ( 1 ) is feasible or not. If not then it may be transformed (or repaired) to solvable one merely by adjusting of its right-hand-side vector. Note that such problems are known as improper ones of the rst kind [1,2].

If problem ( 1 ) is infeasible, we de ne its generalized solution as a solution x^ of the relaxed problem min f(c; x) : Ax = b + ∆^b; x 0g (=: ^) ; where ∆^b =argminf∥∆b∥: Ax = b + ∆b; x set of feasible solutions for ( 3 ) coincides with 0g, ∥ ∥ is Euclidean norm. Obviously, the M = Arg min ∥Ax x 0 b∥ : If ( 1 ) is solvable, then ∆^b = 0, and generalized solution introduced above coincides with an usual one.

As mentioned above, in [10,11] two numerical methods were proposed which being applied to LP problem ( 1 ) give its solution if such solution exists, or automatically deliver its generalized solution x^, if it is not so. Below we discuss how they may be t to parallel calculations in MPI-type paradigm under the assumption that constraints matrix in ( 1 ) has a block-angular structure.

Two types of infeasible constraint angularity are investigated in this paper: a) block-angular matrix with linking columns and b) block-angular matrix with linking rows.

Both variants exploit parallel algebraic operations with matrices. 2

Block-Angular Matrices with Linking Columns Consider the case when constraints matrix in infeasible problem ( 1 ) has the following structure ( 1 ) ( 2 ) ( 3 ) ( 4 ) To take the advantage of this structure for parallelization of the calculation process, make use of the mixed barrier-penalty-type function for problem ( 1 ) Consider the unconstrained smooth optimization problem: nd x^ϵ > 0 such that (ϵ; x^ϵ) = min (ϵ; x) :

x>0 As shown in [10], if optimal set of ( 3 ) is bounded, then an unique minimizer x^ϵ exists for every ϵ > 0, even in improper case. This minimizer satis es to nonlinear vector equation ∇x (ϵ; x) = c + ϵ 1AT (Ax b) ϵX 1e = 0 ; where X is diagonal matrix with coordinates xi on its diagonal, i.e. X = diag(x), e = [1; : : : ; 1]T .

Theorem 1 ([10]). Suppose that the feasible set of the dual problem ( 2 ) has an interior point. Then (c; x^ϵ) ! ^; b

Ax^ϵ = ∆^bϵ ! ∆^b as ϵ ! +0.

Corollary 1. Let the optimal vector x^ of the relaxed problem ( 3 ) be unique. Then x^ϵ ! x^ as ϵ ! +0.

To solve equation ( 7 ), one can apply Newton method as follows x(t + 1) = x(t) tw(t);

w(t) = H 1(ϵ; x(t))∇ (ϵ; x(t)); t = 0; 1; : : : : Here x(0) is an appropriate start point, t is a step parameter (usually t 1), H(ϵ; ) is a Hessian of the function (ϵ; ), ∇x (ϵ; ) is its gradient.

Without going deep into all issues of implementation of this algorithm (e.g. see [15] for details), we only consider how to calculate vector w(t) in parallel. In our case Hence, to nd w(t) one has to solve a sparse linear system with

H(ϵ; x) = ϵ 1AT A + ϵX 2 :

H(ϵ; x(t))w = ∇ (ϵ; x(t)) 0 H1

H2 ( 8 ) where Hi = ϵ 1QiT Qi + ϵXi 2 (i = 1; : : : ; n); Mi = ϵ 1PiT Qi (i = 1; : : : ; n) ; n Mm+1 = ϵ 1 ∑ PiT Pi + ϵXn+21 :

i=1 Diagonal matrices Xi = diag(xi) above correspond to the partition of a primal vector x onto sub-vectors xi according to matrix partition ( 4 ).

Sparse structure of the Hessian gives us the opportunity to solve linear system ( 8 ) in parallel mode (see section 4 for details). 3

Block-Angular Matrices with Linking Rows

Now consider the case when

To exploit this structure for parallelization of the calculation process, make use of the regularized barrier-type function for problem ( 2 ) (ϵ; x) = (b; y) (Ai; y)); ϵ > 0 ; ( 10 ) ϵ n0 2 ∥y∥2 + ϵ ∑ ln (ci

i=1 (ϵ; y^ϵ) = max (ϵ; y) :

y2Ω where Ai denotes the i-th column of the matrix A. Function ( 10 ) is nite and strongly concave on Ω = fy: AT y < cg. Therefore, for every ϵ > 0 there exists an unique vector y^ϵ such that Note that y^ϵ satis es to nonlinear vector equation ∇y (ϵ; y^ϵ) = b

ϵy^ϵ + ϵAD(y) 1e = 0 ; where diagonal matrix D(y) = diag(c diagonal, e = [1; : : : ; 1]T .

AT y) has elements di = ci (Ai; y) on its Theorem 2 ([11]). Let y^ϵ be the maximizer of (ϵ; ) on y. Then the vector x^ϵ with coordinates x^iϵ = ϵ(ci

(Ai; y^ϵ)) 1; i = 1; : : : ; n0 ; is the minimizer of (ϵ; ) on x > 0. Corollary 2. Suppose that the feasible set of the dual problem ( 2 ) has an interior point. Then ϵy^ϵ ! ∆^b as ϵ ! +0.

To solve equation ( 12 ) one can apply Newton method again. Now its iterations are as follows y(t + 1) = y(t) + tw(t);

w(t) = H 1(ϵ; y(t))∇y (ϵ; y(t)); t = 0; 1; : : : : Here y(0) is an appropriate start point, t is a step parameter (usually t is a Hessian of the function (ϵ; ), ∇y (ϵ; ) is its gradient.

In the case under consideration 1), H(ϵ; )

H(ϵ; y) = ϵI + ϵAD(y) 2AT ; where I is the identity matrix of appropriate order. Therefore, to nd w(t) one has to solve a sparse linear system with where

H(ϵ; y(t))w = ∇y (ϵ; y(t)) C C C ; C C A Hi = ϵIi + ϵQiDi(y) 2QiT ; Mi = ϵPiDi(y) 2QiT ; (i = 1; : : : ; n) ; n Mn+1 = ϵIn+1 + ϵ ∑ PiDi(y) 2PiT :

i=1 Here, diagonal matrices Di(y) = diag(ci QiT yi PiT yn+1) correspond to the partition of a dual vector y onto sub-vectors yi according to matrix partition ( 9 ).

Again sparse structure of the Hessian gives us the opportunity to solve linear system ( 13 ) in parallel mode. 4

Parallel Scheme of Implementation It is easy to see that in the both cases above the structure of the Hessians is just the same. Let us discuss how to solve a sparse linear system of the type

H2 . . . ...

Hn MnT M1 M2 : : : Mn Mn+1

M1T M2T 1

0 f1 C C C C C A ( 13 ) ( 14 ) on massively parallel processor (MPP) or cluster of workstations (COW) with MPIinstructions.

For simplicity assume that a number of processors in our MIMD-computer is equal to n + 1. Let processors with indexes from 1 to n be called the slaves and processor with index n + 1 be called the master. Suppose that input data of LP problem ( 1 ) is distributed between the slaves in such a manner that the slave with index i stores the matrix pair (Qi; Pi) in its local memory. Therefore, each slave with index i can calculate its own part of Hessian (Hi; Mi) and some auxiliary matrices with index i concurrently. The master controls the calculation process and, besides, calculates Hn+1 and other matrices with additional index n + 1. All processors exchange with the messages using MPI-instructions like send, gather and broadcast to coordinate their works.

Parallel resolution of system ( 14 ) may be as follows. At rst, all the slaves calculate its own Cholesky decomposition of Hi = LiLiT concurrently, then form the auxiliary matrices Si = L 1MiT , Ri = SiT Si, and send them to the master. The master gathers this matricesi into a global sum Hn+1 = Hn+1 ∑in=1 Ri and then calculates its own Cholesky decomposition of Hn+1 = Ln+1LTn+1. As a result the global Cholesky decomposition is calculated 0 H1

H2

L2 . . .

Ln S1T S2T : : : SnT Ln+1 Further, all the slaves solve their own angular systems Lizi = fi concurrently, then form the products hi = SiT zi and send them to the master. The master gathers this vectors, forms the global sum fn+1 = fn+1 ∑in=1 hi and solves its own angular system Ln+1zn+1 = fn+1. Thereby, a solution of the following system is obtained

L2 . . .

Ln S1T S2T : : : SnT Ln+1 Finally, the master solves the angular system LTn+1wn+1 = zn+1 and broadcasts wn+1 to all the slaves. The slaves solve their angular systems LiT wi = zi Siwn+1 concurrently. Therefore, the following system is solved

LT 2 . . .

S1 S2 . .

.

LTn Sn

LTn+1 It means that all components of w = [w1; : : : ; wn; wn+1] are found.

Note that similar parallel scheme of calculations was applied by author for sparse linear system in [16], where one can nd the estimates of possible speedup Sn+1

n 1 + Cm 1 ! n (m ! 1) ; ( 15 ) as well as an example of implementation scheme for Parallel Matlab and the encouraging results of numerical experiments with large LP problems. In ( 15 ) C > 0 is some technical constant depending upon concrete computer equipment, m characterizes the average dimension of the blocks of the matrix A in problem ( 1 ). If m is comparatively large then the speedup tends to n. 5

Conclusion

For a linear program with block-angular matrix of constraints, two parallel second order optimization methods are proposed. The methods give usual solution if initial problem is solvable, and automatically deliver a solution of some relaxed LP problem, otherwise. The methods use penalty and barrier functions as well as Tikhonov regularization of Lagrangian of initial LP problem. The methods contain a small parameter which tends to zero and minimize the residual of constraints. Parallel calculations of matrix operations follow MPI-type paradigm.

Acknowledgements. This work was supported by Russian Foundation for Basic Research (project no. 16-07-00266).

References

1. Eremin , I. I. , Mazurov, Vl. D., and Astaf'ev, N. N. : Improper Problems of Linear and Convex Programming . Nauka, Moscow. In Russian ( 1983 )

2. Eremin , I. I. Theory of Linear Optimization. Inverse and Ill-Posed Problems Series . VSP. Utrecht, Boston, Koln, Tokyo ( 2002 )

3. Chinneck , J.W. : Feasibility and viability . in: T. Gal, H.J. Greenberg (Eds.), Advances in Sensitivy Analysis and Parametric Programming , Kluwer Academic Publishers, Boston ( 1997 )

4. Ben-Tal , A. ; Nemirovski , A. : Robust Solutions of Linear Programming Problems Contaminated with Uncertain Data . Math. Progr. 88: 3 , 411 { 424 ( 2000 )

5. Amaral , P. ; Barahona, P. : Connections between the total least squares and the correction of an infeasible system of linear inequalities . Linear Algebra and its Appl . 395 , 191 { 210 ( 2005 )

6. Leon , T ; Liern, V. : A Fuzzy Method to Repair Infeasibility in Linear Constrained Problems . Fuzzy Set and Systems . 122 , 237 { 243 ( 2001 )

7. Dax , A. : The smallest correction of an inconsistent system of linear inequalities . Optimization and Engineering . 2 , 349 { 359 ( 2001 )

8. Eremin , I. I. ; Popov , L.D. : Numerical analysis of improper linear programs using the DELTA-PLAN-ES interactive system . Optimization Methods and Software . Vol. 2 . 69{ 78 ( 1993 )

9. Skarin , V. D.: Regularized Lagrange function and correction methods for improper convex programming problems . Proc. Steklov Inst. Math. Mathematical Programming. Regularization and Approximation, suppl. 1 , S116 { S144 ( 2002 )

10. Popov , L. D.: Use of barrier functions for optimal correction of improper problems of linear programming of the 1st kind . Automation and Remote Control . 73 : 3 , 417 { 424 ( 2012 )

11. Popov , L. D.: Dual approach to the application of barrier functions for the optimal correction of improper linear programming problems of the rst kind . Proceedings of the Steklov Institute of Mathematics. 288:1 , 173 { 179 ( 2015 )

12. Fiacco , A.V. ; McCormick , G.P. : Nonlinear Programming: Sequential Unconstrained Minimization Techniques . John Wiley and Sons, Inc: New York-London-Sidney ( 1968 )

13. Tikhonov , A.N. ; Vasil'ev , F.P. : Methods for solving ill-posed extremal problems , in: Mathematics Models and Numerical Methods . Banach center publ., Warszava . 3 , 291 { 348 ( 1978 )

14. Rockafellar , R.T.: Convex Analysis . Princeton University Press: Princeton, N.J. ( 1970 )

15. Gill , Philip E.: Practical optimization / Philip E. Gill, Walter Murray, Margaret H. Wright . London ; New York : Academic Press ( 1981 )

16. Popov , L. D.: Experience in organizing hybrid parallel calculations in the Evtushenko{ Golikov method for problems with block-angular structure . Automation and Remote Control . 75 : 4 , 622 { 631 ( 2014 )