Posterior Constraint Selection for Nonnegative Linear Programming

Posterior constraint optimal selection techniques (COSTs) are developed for nonnegative linear programming problems (NNLPs), and a geometric interpretation is provided. The posterior approach is used in both a dynamic and non-dynamic active-set framework. The computational performance of these methods is compared with the CPLEX standard linear programming algorithms, with two most-violated constraint approaches, and with previously developed COST algorithms for large-scale problems.

(2) 0, ≥ x (3) where c and x are n-dimensional column vectors of objective coefficients and variables respectively; A is an m n × matrix ij a     with 1 n × row vectors , 1, , ; i i m =  a b is an 1 m × column vector; and 0 is an 1 n × vector of zeros.
The non-polynomial simplex methods and the polynomial interior-point barrier-function algorithms are currently the principal two-solution approaches for solving problem , P but for either there are problem instances for which it performs poorly [1].Since the principle use of LP in industrial applications is in binary and integer programming algorithms, however, pivoting algorithms with efficient post-optimality analysis are frequently preferable to interior-point methods.On the other hand, simplex methods often cannot solve large-scale LPs at a speed required by many current applications.The purpose here is to develop an approach for solving a certain class of LPs faster than existing methods.
In this paper we consider the nonnegative linear programming problem (NNLP), which is the special case of P with 0 i ≥ a but 0, 1, , ; NNLPs model a large number of linear programming applications such as determining an optimal driving path for navigation systems using traffic data [2], updating flight status due to weather conditions [3], and detecting errors in DNA sequences [4].NNLPs have the following two important properties: 1) the origin 1, , .min : 0 , 1, , Thus NNLPs have a bounded feasible region and bounded objective function if and only if no column of A is a zero vector.It follows that the boundedness of an NNLP objective function is easily verifiable without computation.

Background
We propose here an active-set method to solve nonnegative linear programming problems faster than current approaches.Our method divides the constraints of problem P into operative and inoperative constraints at each active-set iteration.Operative constraints are those active in the current relaxed subproblem , 1, 2, , r P r =  of P at iteration , r while the inoperative ones are constraints of the problem P not active in .Active-set methods have been studied by Stone [5], Thompson et al. [6], Adler et al. [7], Zeleny [8], Myers and Shih [9], Curet [10], and Bixby et al. [11], among others.The term ''constraint selection technique'' was introduced in [9], while the approaches of [7] and [8] illustrate two distinct classes of active-set methods.When the constraint selection metric for choosing violated inoperative constraints to be added to r P does not depend on the solution * r x , the asso- ciated active-set method is called a prior method.On the other hand, if the constraint selection at r P does depend on * r x , it is called a posterior method.Adler et al. [7] developed a prior method in which a violated inoperative constraint was chosen randomly at each active-set iteration.Zeleny [8] proposed a posterior method in which the inoperative constraint most violated by * r x was added.This method is a classical cutting-plane generation technique and is called VIOL here.VIOL is also used as a pricing rule in delayed column generation [12], as an approach for adding multiple constraints in the interior point cutting plane method of [13], and as part of the sifting algorithm of [11] for column generation.
More recently, Corley et al. [14] developed a geometric prior active-set method for P called the cosine simplex method.At each active-set iteration , r a single violated constraint maximizing the cosine of the angle between i a and c is added to the operative set for r P .This cosine constraint selection crite- rion is equivalent to the "most-obtuse-angle" pivot rule for the modified simplex algorithm introduced by Pan [15], where it was applied to the dual problem for P. Junior and Lins [16] also utilized a cosine criterion to choose an initial basis for the simplex algorithm on P resulting in a fewer number of simplex iterations.
References [17] [18] [19] [20] are most directly related to the current work and involve the authors here.In [17], Corley and Rosenberger proposed the constraint selection metric maximizing ( ) for NNLPs.RAD is a geometric constraint selection criterion for determining the constraints most likely to be binding at optimality.In the associated active-set algorithm of [18], all constraints of (2) are initially ordered by decreasing value of RAD prior to solving an initial bounded problem 0 P by the primal simplex.The dual simplex is then used when violated inoperative constraints are added according to their RAD value.In computational experiments, RAD proved superior to existing linear programming methods for NNLPs.A similar constraint selection metric GRAD was developed in [19] to solve general linear programs (LPs).Finally, in [20] a dynamic active-set method was developed for adding a varying number of violated constraints at r P based on progress at 1 .r P − It was incorporated into both RAD and GRAD to improve the computational results of [18] and [19].

Overview
In this paper a posterior constraint selection metric NVRAD is developed for NNLPs.NVRAD may be considered as a posterior version of RAD.The posterior NVRAD is then implemented in the dynamic framework of [20].It should be noted that a constraint selection metric and the associated active-set method are identified by the same name -in this case NVRAD.For the active-set method NVRAD, we provide extensive computational extensive computational experiments to show that it solves NNLPs faster than other computational methods, including RAD and various versions of the existing posterior active-set method VIOL described above.More specifically, in Section 2 we state the posterior constraint selection metric NVRAD and provide a geometric interpretation.A dynamic version of NVRAD for NNLPs is then developed.In Section 2 we extend NVRAD to a hybrid approach HYBR, where RAD and NVRAD are alternated.In Section 3, computational results are presented.NVRAD is shown to be significantly faster for NNLPs than all CPLEX solvers, as well as faster than VIOL and RAD.HYBR appears slightly faster than NVRAD.In Section 4, we present conclusions.Throughout the paper, both a constraint selection metric and the associated active-set algorithm are identified by the same name-RAD or NVRAD, for example.The use the term should be clear from context.The active-set algorithm itself is called a COST, i.e., a "Constraint Optimal Selection Technique".

Definition and Interpretation
Let * r x be the current optimal solution for some r P with a perpendicular dis- tance x a a to a violated hyperplane .P For this reason, in [18] the inoperative constraint maximizing ( ) c is deemed the best constraint to add to r P based on prior information.We combine this prior information (4) with the posterior information on the right side of ( 5) by multiplying them to give Equation ( 6) thus incorporates global information from RAD with local information at * r x , and our posterior active-set method adds to r P an inopera- tive constraint * i for which ( ) We mention that the term 2 i b in the denominator of (7) works better than simply .
i b This fact was established in computational results not reported here but obtained to support the above derivation.

The Dynamic Active-Set Algorithm
A dynamic version of RAD was developed by the authors in [20].A similar ap-proach is now used for NVRAD.Let * r x be the optimal extreme point for r P , with r θ the angle between * r x and .c Then is nonnegative since r P is also an NNLP.We would like to decrease r θ at each active-set iteration so that * r x points more in the same direction as the gradient c of the objective function in (1).We adapt our dynamic heuristic of [20] that adds a varying number of violated inoperative constraints to r P according to the progress made As our ideal goal, let 0 When 0 r θ = , it follows from ( 9) that ( ) as a measure of the performance of our active-set method at iteration .r The value of ( ) * r r δ x decreases as r θ decreases.Such a decrease usually occurs as * r x approaches an optimal extreme point of P itself.
The dynamic COST NVRAD for solving NNLPs is described as follows.Constraints are initially ordered by the RAD constraint selection metric (4).To construct 0 P we choose constraints from (2) in descending order of RAD (since there is no * r x ) until the 0 A matrix of has no 0 column, i.e., until each variable j x has an 0. ij a > These selected constraints become the constraints of 0 P , and we say that the variables are covered by the inequality constraints of the initial problem.0 P is then solved by the primal simplex to achieve an initial solu- tion * 0 , x and With [ ] ⋅ denoting the greatest integer function, let ( ) ( ) where 1 200.

ϕ =
The value of r ϕ is an upper bound on the possible number of violated constraints that can be added at active-set iteration .r The actual number added is The active-set function r ϕ increases at every iteration since the optimal value of the objective function for r P is usually less affected by a constraint with a small value of ( 6) than one with a large value.
Hence, more violated constraints should be added as r increases.Equation The pseudocode for dynamic NVRAD algorithm is as follows.

A Hybrid Approach
A reasonable conjecture is that that combining the global information of RAD and the local information of NVRAD might be advantageous.Therefore, we will also consider an approach that alternates the dynamic RAD and NVRAD metrics in a single algorithm at even and odd iterations, respectively, to yield a hybrid COST designated here as HYBR.The results obtained for HYBR demonstrate that combining posterior and prior COSTs may be superior to either a prior or posterior approach by itself.

Computational Experiments
Dynamic NVRAD is compared in this section with the CPLEX primal simplex, dual simplex, and barrier methods.It is also compared with the prior active-set method RAD and the standard posterior active-set method VIOL, as well as to a normalized version of VIOL called NVIOL that was superior to VIOL in computational results not reported here.Both dynamic and multi-bound, multi-cut versions of NVRAD were compared to dynamic and multi-bound, multi-cut versions of the other active-set methods for insight into the individual merits of the dynamic and posterior approaches.

Problem Instances
Five sets of NNLPs from [18] are used to evaluate the performance of the dynamic posterior COST NVRAD.Each of Sets 1 -4 contains 105 randomly generated NNLPs at 21 different density levels ranging from 0.005 to 1, and four ratios of ( m constraints)/( n variables) ranging from 200 to 1.The ratios for Sets 1 -4 are 200, 20, 2, and 1, respectively.For each of Sets 1 -4, there are five problem instances per combination of density level and ratio.In these problem sets, randomly generated real numbers between 1 and 5, 1 and 10, and 1 and 10 were assigned to the elements of , , A b and , c respectively.To prevent any constraint of P from having the form of an upper bound on some variable, each constraint is required to have at least two nonzero ij a .Next, problem Set 5 of NNLPs is a set of large-scale problems with 5000 variables and 1,000,000 constraints.In this set, real numbers between 1 and 100 are assigned to the elements of b and c with densities p ranging from 0.0004 to 0.06.Again, each con-straint is required to have at least two nonzero ij a .

CPLEX Preprocessing
Two CPLEX parameters for solving linear programming are discussed here.The preprocessing pre-solve indicator (PREIND) and the preprocessing dual setting (PREDUAL) are the two parameters that CPLEX uses for solving linear programming.Preprocessing pre-solver is enabled with the parameter setting PREIND = 1 (ON), which reduces both the number of variables and the constraints before any type of algorithm is used.The pre-solver routine in CPLEX is disabled by setting PREIND = 0 (OFF).The second preprocessing parameter in CPLEX affecting the computational speed is PREDUAL.By setting parameter PREDUAL = 0 (ON) or PREDUAL = −1 (OFF), CPLEX automatically selects whether to solve the dual of the original LP or not, respectively.Both PREIND and PREDUAL were turned off for CPLEX when CPLEX was used as part of NVRAD or HYBR.However, all computational results reported here for any individual CPLEX solver had PREIND and PREDUAL turned on.In other words, our NVRAD was compared to CPLEX at its fastest setting.CPLEX would choose automatically whether to solve either the primal or dual, whichever seemed best.Moreover, preprocessing would substantially reduce the size of any problem P by removing appropriate rows or columns of the constraint matrix A before applying the primal simplex, dual simplex, or interior-point barrier method.In fact, much of the speed of the CPLEX solvers is due to its proprietary preprocessing routines.

Computational Results
The experiments were performed on an Intel ® Core TM 2 Duo X9650 3.00 GHz processor with a Linux 64-bit operating system and 8 GB of RAM.The COST NVRAD uses the IBM CPLEX 12.5 callable library to solve 0 P by the primal simplex and then , 1, 2, r P r =  by the dual simplex when selected constraints are added to The results of Table 1 for Set 1 compare NVRAD to VIOL, as well as to both a dynamic and non-dynamic version of NVIOL.In addition, the dynamic NVRAD described in Section 2.2 was compared to a non-dynamic NVRAD that applies the multi-cut and multi-bound technique of [18].The dynamic version was significantly faster.The efficacy of the dynamic approach was further demonstrated by the fact that in higher density problems a dynamic version of NVIOL was up to 21 times faster than the multi-cut, multi-bound NVIOL.Overall, dynamic NVRAD was faster than VIOL and NVIOL on every problem instance.
In Table 2, the CPU times of the test problems solved by dynamic NVRAD are compared with those for RAD.In problem Set 1, RAD is slightly faster than NVRAD over all densities and averages 3.98 compared to 4.55 seconds.However, in problem Set 2, the average computation times for RAD and dynamic Table 1.CPU times for multi-cut, multi-bound and dynamic active-set approaches on problem Set 1 for random NNLPs with 1000 variables, 200,000 constraints, and 1 5,  3 presents the CPU times for problem Sets 1 -4 solved by dynamic versions of both RAD and HYBR.In Table 3 HYBR is superior to RAD.Moreover, a comparison of Table 3 with Table 2 shows that HYBR is also slightly better than dynamic NVRAD on these problem sets.Such observations suggest that combining the global information of RAD and the local information of NVRAD gives a superior performance than either RAD or NVRAD by itself.We note further that HYBR can probably be improved.However, it is not our goal to seek the optimal combination of RAD and NVRAD in HYBR since an optimal combination would likely differ depending on various factors such as density and the ratio m/n.Table 4, taken from [18], provides a comparison of the posterior COST NVRAD with the standard CPLEX solvers.Comparing the results of Table 4 for the CPLEX solvers with the results for NVRAD in Table 2 shows that NVRAD was significantly faster across virtually all ratios m n and all densities.For example, the primal simplex was the most robust CPLEX solver, but on the average across all densities the primal simplex took approximately 3 to 14 times more CPU time for the different rations m n than NVRAD.For the dual simplex, the av- erage CPU across all densities was approximately 15 to 50 times greater than NVRAD over the different ratios.However, the CPLEX barrier method was slightly faster than NVRAD in problem instances with 20 m n = and with densities less than 0.02.On the other hand, when the density reached 0.08 for 20 m n = , NVRAD was already more than ten times faster than the barrier solver.Furthermore, note that average CPU times in Table 4 greater than 3000 seconds (50 minutes) at any density were not reported.This situation occurred for the CPLEX barrier solver for the ratios 1, 2, 20, and 200 with densities at least 0.3, 0.4, 0.5, and 0.75, respectively.
Finally, for large-scale, low-density test problems with 5000 n = and 1, 000, 000.m = Table 5 compares dynamic NVRAD to multi-cut and multi-bound RAD, VIOL, NVIOL, and NVRAD, as well as to the CPLEX primal simplex, dual simplex, and barrier solvers.Only the prior COST RAD was competitive.NVRAD averaged 63.45 seconds overall as compared to 71.79 for RAD.It should be noted that the highest density used in problem Set 5 was 0.0600 since the CPLEX solvers could not solve denser problems of such magnitude in a reasonable amount of time.Average CPU times greater than 2400 seconds (40 minutes) at any density were not reported in Table 5.This situation occurred beginning at some individual threshold density level for each CPLEX solver.

Conclusion
An efficient posterior COST called NVRAD was developed here for NNLPs to utilize both prior global information and posterior local information.The associated constraint selection metric NVRAD is a heuristic, so a geometric interpretation was presented to offer insight into its performance.NVRAD's inherent active-set efficiency was enhanced by a dynamic approach varying the number of constraints added at each iteration.In addition to NVRAD, adynamic active-set approach HYBR was also proposed.HYBR alternates between the posterior method NVRAD and the prior method RAD.To check their performance, both NVRAD and HYBR were used to solve five sets of large-scale NNLPs.Dynamic NVRAD outperformed the previously developed COST RAD, as well as the standard posterior cutting-plane method VIOL.Dynamic NVRAD significantly outperformed the CPLEX primal simplex, dual simplex, and barrier solvers.On

rP 1 
In our active-set method we iteratively solve of P after adding one or more violated inoperative con- straints from (2) to 1 r P − until the solution * r x to r P is a solution to P .
side of (5) can be interpreted from the left side of (5) as selecting a violated constraint giving the deepest cut based on information derived from * r x .But from[18], the expression i i b a c on the right side of (4) is the distance from the origin to the hyperplane , c i.e., the direction of steepest ascent for the objective function (1) of the NNLP .

the number of constraints of problem P violated by * r x . Then at iteration 1 r
calculated; and the percentage of im- provement made in reducing the angle between vectors * r x and c at iteration r is measured by

( 13 )
represents one approach for doing so.If > In other words, a much larger number and perhaps all of the remaining violated constraints could be added.NVRAD stops when 0, r γ = i.e., when there are no more violated constraints.

Step 3 -
Optimized false ←9: end whileStep 2-Using the primal simplex method, obtain an optimal * 0 x for the ini- tial problem.Perform the following iterations until an answer to problem P is found.

1 rP
− .The CPU times shown in the tables below represent the average computation time of five problem instances at each density level.
Solve the following r P by the dual simplex method to obtain * .
over all densities are 19.07 and 16.86 seconds, respectively.For Set 3, dynamic NVRAD is superior to RAD averaging 38.91 seconds compared to 41.87 seconds.Similarly, for Set 4 the averages are 41.87 for NVRAD as compared to 46.98 for RAD.Thus the results of Table2affirm NVRAD's ability to add appropriate constraints at each iteration.The results for Set 1 simply reflect how well the prior COST RAD performs when m is very much larger than .n ++ Average of 5 instances at each density.−− Used CPLEX presolve = OFF and predual = OFF.NVRAD

Table 2 .
CPU times for multi-cut, multi-bound and dynamic active-set approaches on problem Sets 1 -4 for random NNLPs with 1 5, ++ Average of 5 instances of LP at each density.−− Used CPLEX presolve = OFF and predual = OFF.

Table 5 .
CPU times for NVRAD versus RAD, VIOL, NVIOL, and the CPLEX solvers on problem Set 5 for random NNLPs with 5000 variables, 1,000,000 constraints, and 1 5, ++ Average of 5 instances at each density.b Runs with CPU times > 2400 s are not reported.−− Used CPLEX presolve = OFF and predual = OFF.+ Used CPLEX presolve = ON and predual = ON.