=Paper=
{{Paper
|id=Vol-1623/papermp11
|storemode=property
|title=Parallel Newton Methods for Numerical Analysis of Infeasible Linear Programming Problems with Block-Angular Structure of Constraints
|pdfUrl=https://ceur-ws.org/Vol-1623/papermp11.pdf
|volume=Vol-1623
|authors=Leonid D. Popov
|dblpUrl=https://dblp.org/rec/conf/door/Popov16
}}
==Parallel Newton Methods for Numerical Analysis of Infeasible Linear Programming Problems with Block-Angular Structure of Constraints==
https://ceur-ws.org/Vol-1623/papermp11.pdf
Parallel Newton Methods for Numerical Analysis
of Infeasible Linear Programming Problems
with Block-Angular Structure of Constraints
Leonid D. Popov
1
Krasovskii Institute of Mathematics and Mechanics UB RAS, Ekaterinburg, Russia
2
Ural Federal University, Ekaterinburg, Russia
popld@imm.uran.ru
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 mathe-
matical programming theory (see [1,2,3,4,5,6] and others). They arise naturally in many
applications of linear programming (LP) such as network flows analysis problems, prob-
lems of portfolio selection and scheduling problems which reflect 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 find 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 fit 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
254 L. D. Popov
1 The Problem State and Preliminaries
Consider the linear program
min {(c, x) : Ax = b, x ≥ 0} (1)
and its dual one
max {(b, y) : A⊤ y ≤ c} ; (2)
where A ∈ IRm0 ×n0 , c ∈ IRn0 and b ∈ IRm0 are given, x ∈ IRn0 and y ∈ 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 first kind [1,2].
If problem (1) is infeasible, we define its generalized solution as a solution x̂ of the
relaxed problem
min {(c, x) : Ax = b + ∆b̂, x ≥ 0} (=: υ̂) , (3)
where ∆b̂ =argmin{∥∆b∥: Ax = b + ∆b, x ≥ 0}, ∥ · ∥ is Euclidean norm. Obviously, the
set of feasible solutions for (3) coincides with
M = Arg min ∥Ax − b∥ .
x≥0
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 fit
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
Q1 P1
Q2 P2
. .
A= .. .. . (4)
Qn−1 Pn−1
Qn Pn
Parallel Newton Methods for Numerical Analysis of Infeasible LP 255
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)
1 ∑n0
Φ(ϵ; x) = (c, x) + ∥Ax − b∥2 − ϵ ln xi , ϵ>0 . (5)
2ϵ i=1
Consider the unconstrained smooth optimization problem: find x̂ϵ > 0 such that
Φ(ϵ; x̂ϵ ) = min Φ(ϵ; x) . (6)
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 satisfies to nonlinear vector
equation
∇x Φ(ϵ; x) = c + ϵ−1 AT (Ax − b) − ϵX −1 e = 0 , (7)
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) − αt w(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
H(ϵ; x) = ϵ−1 AT A + ϵX −2 .
Hence, to find w(t) one has to solve a sparse linear system
H(ϵ; x(t))w = ∇Φ(ϵ; x(t)) (8)
with
H1 M1T
H2 M2T
.. ..
H(ϵ; x) = . . ,
Hn MnT
M1 M2 . . . Mn Mn+1
256 L. D. Popov
where
Hi = ϵ−1 QTi Qi + ϵXi−2 (i = 1, . . . , n), Mi = ϵ−1 PiT Qi (i = 1, . . . , n) ,
∑
n
−1 −2
Mm+1 = ϵ PiT Pi + ϵXn+1 .
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
Q1
Q2
..
.
A= . (9)
Qn−1
Qn
P1 P2 ... Pn−1 Pn
To exploit this structure for parallelization of the calculation process, make use of the
regularized barrier-type function for problem (2)
ϵ ∑n0
Ψ (ϵ; x) = (b, y) − ∥y∥2 + ϵ ln (ci − (Ai , y)), ϵ>0 , (10)
2 i=1
where Ai denotes the i-th column of the matrix A. Function (10) is finite and strongly
concave on Ω = {y: AT y < c}. Therefore, for every ϵ > 0 there exists an unique vector
ŷϵ such that
Ψ (ϵ; ŷϵ ) = max Ψ (ϵ; y) . (11)
y∈Ω
Note that ŷϵ satisfies to nonlinear vector equation
∇y Ψ (ϵ; ŷϵ ) = b − ϵŷϵ + ϵAD(y)−1 e = 0 , (12)
where diagonal matrix D(y) = diag(c − AT y) has elements di = ci − (Ai , y) on its
diagonal, e = [1, . . . , 1]T .
Theorem 2 ([11]). Let ŷϵ be the maximizer of Ψ (ϵ; ·) on y. Then the vector x̂ϵ with
coordinates
x̂iϵ = ϵ(ci − (Ai , ŷϵ ))−1 , i = 1, . . . , n0 ,
is the minimizer of Φ(ϵ; ·) on x > 0.
Parallel Newton Methods for Numerical Analysis of Infeasible LP 257
Corollary 2. Suppose that the feasible set of the dual problem (2) has an interior
point. Then ϵŷϵ → ∆b̂ as ϵ → +0.
To solve equation (12) one can apply Newton method again. Now its iterations are
as follows
y(t + 1) = y(t) + αt w(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 ≡ 1), H(ϵ; ·)
is a Hessian of the function Ψ (ϵ; ·), ∇y Ψ (ϵ; ·) is its gradient.
In the case under consideration
H(ϵ; y) = ϵI + ϵAD(y)−2 AT ,
where I is the identity matrix of appropriate order. Therefore, to find w(t) one has to
solve a sparse linear system
H(ϵ; y(t))w = ∇y Ψ (ϵ; y(t)) (13)
with
H1 M1T
H2 M2T
.. ..
H(ϵ; x) = . . ,
Hn MnT
M1 M2 . . . Mn Mn+1
where
Hi = ϵIi + ϵQi Di (y)−2 QTi , Mi = ϵPi Di (y)−2 QTi , (i = 1, . . . , n) ,
∑
n
Mn+1 = ϵIn+1 + ϵ Pi Di (y)−2 PiT .
i=1
Here, diagonal matrices Di (y) = diag(ci − QTi 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
H1 M1T w1 f1
H2 M2T
w2 f2
. . . .
.. .. .. = .. (14)
Hn Mn T wn fn
M1 M2 . . . Mn Mn+1 wn+1 fn+1
258 L. D. Popov
on massively parallel processor (MPP) or cluster of workstations (COW) with MPI-
instructions.
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 first, all the slaves calculate
its own Cholesky decomposition of Hi = Li LTi concurrently, then form the auxiliary
matrices Si = L−1 T T
i Mi , Ri = Si Si , and send them to ∑the master. The master gath-
n
ers this matrices into a global sum H̄n+1 = Hn+1 − i=1 Ri and then calculates its
own Cholesky decomposition of H̄n+1 = Ln+1 LTn+1 . As a result the global Cholesky
decomposition is calculated
T
H1 M1T L1 L1 S1
H2 M2T L2 LT2 S2
.. .. .. .. ..
. . = . . . .
Hn Mn T Ln Ln Sn
T
M1 M2 . . . Mn Mn+1 S1T S2T . . . SnT Ln+1 LTn+1
Further, all the slaves solve their own angular systems Li zi = fi concurrently, then
form the products hi = SiT zi and send them ∑to the master. The master gathers this
n
vectors, forms the global sum f¯n+1 = fn+1 − i=1 hi and solves its own angular system
Ln+1 zn+1 = f¯n+1 . Thereby, a solution of the following system is obtained
L1 z1 f1
L2 z2 f2
.. .. ..
. . = . .
Ln zn fn
S1T S2T . . . SnT Ln+1 zn+1 fn+1
Finally, the master solves the angular system LTn+1 wn+1 = zn+1 and broadcasts wn+1 to
all the slaves. The slaves solve their angular systems LTi wi = zi − Si wn+1 concurrently.
Therefore, the following system is solved
T
L1 S1 w1 z1
LT2 S2
w2 z2
.. . . .
. .. .. = .. .
LTn Sn wn zn
LTn+1 wn+1 zn+1
It means that all components of w = [w1 , . . . , wn , wn+1 ] are found.
Parallel Newton Methods for Numerical Analysis of Infeasible LP 259
Note that similar parallel scheme of calculations was applied by author for sparse
linear system in [16], where one can find the estimates of possible speedup
n
Sn+1 ≈ → n (m → ∞) , (15)
1 + Cm−1
as well as an example of implementation scheme for Parallel Matlab and the encour-
aging 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 regular-
ization 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 Contam-
inated 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 Prob-
lems. Fuzzy Set and Systems. 122, 237–243 (2001)
7. Dax, A.: The smallest correction of an inconsistent system of linear inequalities. Opti-
mization 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)
260 L. D. Popov
9. Skarin, V. D.: Regularized Lagrange function and correction methods for improper convex
programming problems. Proc. Steklov Inst. Math. Mathematical Programming. Regular-
ization 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 cor-
rection of improper linear programming problems of the first kind. Proceedings of the
Steklov Institute of Mathematics. 288:1, 173–179 (2015)
12. Fiacco, A.V.; McCormick, G.P.: Nonlinear Programming: Sequential Unconstrained Min-
imization 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: Math-
ematics 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 Con-
trol. 75:4, 622–631 (2014)