Experiments with Active-Set LP Algorithms Allowing Basis Deﬁciency

: An interesting question for linear programming (LP) algorithms is how to deal with solutions in which the number of nonzero variables is less than the number of rows of the matrix in standard form. An approach is that of basis deﬁciency-allowing (BDA) simplex variations, which work with a subset of independent columns of the coefﬁcient matrix in standard form, wherein the basis is not necessarily represented by a square matrix. We describe one such algorithm with several variants. The research question deals with studying the computational behaviour by using small, extreme cases. For these instances, we must wonder which parameter setting or variants are more appropriate. We compare the setting of two nonsimplex active-set methods with Holmström’s TomLab L P S IMPLEX V 3.0 commercial sparse primal simplex commercial implementation. All of them update a sparse QR factorization in M ATLAB . The ﬁrst two implementations require fewer iterations and provide better solution quality and running time.


Introduction
Since the introduction of the simplex method for linear programming (LP) in 1947, there has been interest in viewing algorithms from a nonlinear optimisation point of view.Most well known and developed are the ideas of interior point methods, established at the time of the first barrier methods of Dikin in 1967.Although the questions on the simplex method originate from its beginning in 1947, there is still a lot of ongoing research; in the last three years one can find publications such as [1][2][3][4][5][6][7][8], just to cite a few.Various aspects of the method are investigated in these works, where Im and Wolkowicz [6] specifically also focus on degeneracy and related issues.However, these recent works neither work with a nonsquare basis, nor do they include a sparse linear algebra procedure to efficiently solve the sparse least squares subproblems that we need in our algorithmic approach.
A generalization of the well-known dual simplex method can be found in [9,10], which discuss relaxing operations toward nonsimplex steps.One can find more recent considerations in this line in [11,12] with the incorporation of steepest-edge pivoting rules (see also [13,14]).Gill and Murray [9] presented a nonsimplex active-set algorithm for solving a problem (P) starting with a feasible point x 0 .Descriptions and analyses can also be found in [14][15][16].We investigated and developed a nonsimplex active-set (NSA) method for a dual linear program in standard form (D), [17].A challenge for the development is to have to deal with high degeneracy, which for a simplex-like method provides basisdeficiency, i.e., the dual variable set has less than n variables.Therefore, we aimed for the development of a method, which is also called a basis deficiency-allowing (BDA) simplex variation (see [11,12,18]).
This paper presents a study of the computational behavior of a MATLAB implementation on a set of extreme cases and sparse inequality-constrained linear programs customized from the literature.We compare two nonsimplex active-set methods with the commercial solver of Holmström called TomLab LPSIMPLEX V3.0 sparse primal simplex implementation.Various numerical linear algebra procedures are used in MATLAB to keep a sparse QR factorization of the tall-and-skinny basis matrix updated.To elaborate upon the experiments, we introduce the notation and the main algorithm and its implementation in Section 2. We report on computational findings in Section 3. Section 4 summarises our findings.

Materials and Methods
Consider the nonsymmetric primal dual relation in linear programming, wherein we deviate from the usual notation here in exchanging b and c, x and y, n and m, and (P) and (D) in the notation of e.g., [14].Now consider where A ∈ R n×m with m ≥ n and rank(A) = n.Let F and G denote the feasible region of (P) and (D), respectively.Here, c and b are, respectively, the primal and dual gradient of the cost functions to be optimized, whereas a j and b j are, respectively, the coefficient vector and right-hand side value of the jth primal constraint a T j x − b j ≥ 0. The idea is the following.We start with a dual feasible point y 0 .From there, we separate the index set [1 : m] := {1, 2, . . ., m} into an ordered basic set B k with m k := |B k | elements, with m k ≤ n and its complement representing the columns in N k := [1 : m] \ B k , which can have more than m − n zero elements.Let A k ∈ R n×m k and N k ∈ R n×(m−m k ) be submatrices of A formed by the columns corresponding to B k and N k , respectively.Notice that rank(A k ) = m k .We will use the notation R(A) to represent the range of the columns in A, while null(A) denotes its null space, and |A| denotes the number of its nonzeros.Column j of A is denoted by a j .Similarly, a T i is row i of matrix A T .Moreover, b B k represents an m k vector with the elements of B k of b and b N k has the nonbasic elements.The index k will be left out if it is clear from the context, • 2 2 denotes the square of the Euclidean 2-norm, and symbol § will be used in citations as an abbreviation of chapter/section for brevity purposes.

Algorithm
The algorithm works like the simplex method in determining in each iteration a nonbasic variable that enters the basis and a basic variable(s) leaving the basis.The big difference is that it can handle basis deficiency.This means that m k ≤ n and m k may reduce and grow.A pseudocode is given in Algorithm 1.

•
Notice that in line 4, if we have a complementary pair of primal dual feasible points, then x k is optimal for (P).• Line 15 determines step size and basis-leaving variable(s), unless in line 12 we found that we have an unbounded dual objective function.
The algorithm does not require a primal feasible starting point.In that sense, it is an exterior method.In earlier studies, we have shown [19] that the algorithm is equivalent to the primal BDA simplex variation given by Pan in [11].However, we put emphasis on a description, which is easy to understand due to a geometrical interpretation and its independence of implementation details.The algorithm also fits in the primal-feasibility search loop of the Sagitta method described by Santos-Palomo [20], which is related to a loop used in dual active-set methods for quadratic programming.Moreover, Li [21] extended the algorithm to deal with upper and lower bounds on y.
Practically, one obtains the primal-feasibility phase I of Dax's 1978 [22] idea for Rosen's 1960 paper [10] when one replaces the min-ratio test by a most-obtuse-angle row rule to determine the leaving variable and removing an objective function consideration.In this context, we focused earlier [23] on the classical example of Powell [24], which illustrates the cycling behaviour of a most-obtuse-angle simplex pivoting rule described in [12,25].
1: Set k ← 0, y 0 feasible point of (D) with corresponding B 0 and A 0 2:  Step-size τ ← min q∈B k { The algorithm ends with an optimal solution for (P) and (D).However, the solution does not necessarily correspond to a square basis.Notice that the algorithm maintains dual feasibility and complementary slackness.In fact, two algorithmic strategies now emerge.One of them, Reliable Sagitta, relies on starting Algorithm 1 from a dual feasible point.The other one, Sparse Sagitta-abbreviated as Sagitta-starts Algorithm 1 as soon as a solution of Ay = c is found with a suitable, linearly independent subset of columns of A where y is possibly violating nonnegativity.These two alternatives, along with the chosen touchstone for a fair comparison, are outlined in Section 2.2.

Implementation
To run the algorithm, we first find a feasible dual vertex, i.e., a vector y ≥ 0 such that Ay = c with columns of A corresponding to positive components of y being linearly independent.This is called a phase-I procedure.One option is to use the procedure of the nonnegative least squares (NNLS) approach described in [19,26].One can apply the NNLS algorithm [27] or Dax algorithm [28] to the positive semidefinite quadratic problem min which implies solving a sequence of unconstrained least squares problems of the form min A k z − c 2 (NNLS) or of the form min A k w − s k 2 (Dax).Such an approach does not require artificial variables and obtains a feasible direction of the feasible region F (i.e., a so-called direction of F ) or a dual basic feasible solution with m k ≤ n.It appears that the computation can be described in terms of a primal null-space descent direction [29].Instead of (2), one can focus on its Wolfe dual, which is a least-distance, strictly convex quadratic problem where v := A T c.This also implies solving a sequence of least squares problems.To solve (3), one can use any of the primal, dual, and primal dual active-set methods introduced in [30].Note that a primal direction d := u − c is easy to generate.
Once a dual feasible vertex y * is obtained, a phase-II procedure is performed to reach optimality.A traditional way to continue is to artificially enlarge B from m k to n columns and then apply the primal simplex method to (D) starting from a degenerate dual vertex of G.However, the interesting approach we follow is not to artificially enlarge B, but to use y * as the starting point for a nonsimplex active-set method for linear programs in standard form as sketched in Algorithm 1.We shall refer to this first algorithm for solving linear programs, such as Reliable Sagitta (RELSAG V1.1), for which termination is guaranteed under a dual nondegeneracy assumption (see [17]).
The second variant of phase I to solve (2) adapts the inner loop, called Loop B in [27] and step S2 in [29].This loop sequentially removes indices from B of those columns of A k corresponding to negative multiplier least-squares estimates ȳB i of actual Lagrange multipliers y B i for (P).We substitute this loop by the following steps to determine the search direction: Compute temporal variables ȳB solving min Then we set y ← [ ȳB ; 0], and take d ← A k ȳB − c and go on to prepare the next iteration.Notice that ȳB is not constrained to be nonnegative.This NNLS modification finishes with a dual solution y * verifying Ay = c, but nonnegativity does not hold in general.Hence the primal-feasibility search loop is not ensured to start from a feasible solution y ∈ G, and termination is not guaranteed.In the event that this loop ends, we can check whether y ≥ 0 stops with the optimal solution, or deletes a negative element y B i to restart the whole process.Note that we have delayed the elimination of negative multipliers until they have become actual multipliers.We shall refer to this second algorithm for solving linear programs such as Sparse Sagitta (SPASAG V1.1).
The comparison of the two nonsimplex active-set methods described above was done with the commercial code of Holmström TomLab LPSIMPLEX V3.0 [31] sparse primal simplex implementation, applied to solve (D).Details of this MATLAB solver can be found in http://tomlab.biz/products/base/solvers/lpSimplex.php(accessed on 15 November 2022).This code was previously known as LPSOLVE, but should not be confused with the freely available compiled code lp_solve from http://lpsolve.sourceforge.net/5.5/(accessed on 15 November 2022).To make results comparable, the input of the problems in all the solvers are taken the same way.However, LPSIMPLEX deals with the columns of A, whereas RELSAG and SPASAG deal with the rows of A T ; hence all of them are under the same conditions when solving (P) and (D).
We also considered comparison with other available software, which appeared less appropriate.

1.
LINPROG V2.1 [32] included in MATLAB Optimization Toolbox v2.1.Its active-set method does not allow sparse arithmetic as described above.Other available MATLAB codes like those of Morgan [33] have not been considered because he ultimately rejected the sparse dynamic LU updating that he implemented in favour of a systematic LU restart made from scratch on the basis matrix after each column exchange.

2.
We did not compare our code to commercial or freely available compiled codes like MOSEK [34] or MINOS [35], because the running time measurements are not comparable with those of our interpreted MATLAB code.Although solution time is not the only metric that we have used here (number of iterations, quality of solutions), we consider that our solvers must be at least comparable in running time with implementations developed under the same conditions.

3.
It is not our aim to use the MATLAB environment to exclude alternative implementations like CPLEX, XPRESS-MP, COIN-OR, or lp_solve, to name just a few.We plan to compare against them with larger problems when we develop our own implementation.Moreover, although it is part of the folklore that sparse QR-based methods can suffer more fill-in than similar LU based ones, we expect that easier updating and downdating techniques as those of CHOLMOD described recently by Davis and Hager in [36] can be used to show the scalability of the developed implementations to larger models.

Results
This section first describes details of the used settings in Section 3.1 and then focuses on the question of how to make use of specific sparse matrix computations.We discuss results for small, extreme cases from the literature in Section 3.2 to find out which settings are more appropriate for which type of instances.As we will show, this might provide some advantage compared to a general commercial solver like Tomlab.In Section 3.3, we deal with features to handle sparse instances.In earlier work, we also reported on computational tests with 23 out of the smallest 31 NETLIB LP problems, which were originally put forward in standard form (D) (i.e., those with no simple bounds on our dual variables), which are currently being analyzed by Im and Wolkowicz [6].The details of these tests can be found in [37], and a summary of these tests has been included in [17].

Settings
The used implementation (called Sagitta) is intended to find a minimum feasible point of (P), without using an initial feasible point.Default values for the parameters are  Our implementation does not make use of simple bounds.Furthermore, it uses a null-space steepest-descent direction d as the search direction.A constraint a T j x ≥ b j is said to be opposite to direction d if a T j d < 0. In the presence of descent direction d (see Table 1), one of the available addition criteria is to add the constraint a T j x ≥ b j most opposite with Euclidean normalization to the descent direction if o(1)='t' by default.The implementation chooses the opposite −c of the cost vector if o(1)='b', or chooses the difference between the descent direction and the cost vector if o(1)='p'.The selection is made among those constraints opposite to d in the last two cases.The implementation can also select not to normalize using the options o(1)='s', o(1)='a', and o(1)='o', respectively.It is worth noting that all criteria given in Table 1 have a strong geometrical interpretation, trying to take advantage of the fact that both d and −c are descent directions.Setting o(5) controls the choice of the addition criterion in the absence of a descent direction (see Table 2).It chooses the most violated constraint with 2 norm if o(5)='t' by default or without normalization if o(5)='s'.A constraint deletion (Table 3) can be performed both in the presence or in the absence of a descent direction.The former (only available for RELSAG) chooses the constraint with the most negative multiplier estimate with 2 norm if o(5)='t' by default or without normalization if o(5)='s'.The latter (only available for SPASAG) works in a similar way by using actual multipliers instead of multiplier estimates when there is no violated constraint.Moreover, in a constraint exchange this eventuality may be needed in the case of primal feasibility; the constraint to be deleted is chosen in both algorithms by using a min-ratio test similar to that used in the primal simplex method [20,37].
The criteria shown in Table 2 can be thought of as being different ways by which to weigh the residues of the violated constraints.When all of them have unit weight, the classical Dantzig rule applies o(5)='s' and when the weight for constraint i is a i −1 , rule o(5)='t' is selected.In this way, steepest-edge rules work with weight Note that, in the simplex case, a i ∈ R(A k ) for all nonbasic i because of the regularity of A k , but the compatibility of A k δ i B = a i is not guaranteed in the nonsimplex case, and hence we have to rely on the pseudoinverse A † k := (A T k A k ) −1 A T k (see e.g., [39] (p.17, Equation (1.2.27))).Thus, One can approximate the weights of the steepest edges by 1 Notice that the factor A † k −1 is shared for all i.This reasoning implies that rule o(5)='t' can be thought of as a rough approximate steepest-edge rule.An easier rule to approximate 1 δ i in a simplex context has been described by Świe ¸tanowski [40].Both nonsimplex active-set methods have been implemented by using range-space techniques by maintaining the sparse QR factorization of A k .Specifically, the purpose of the modified sparse least squares (MSLS) toolbox v1.3 (that we have developed and described in [29] ( §3)) is the sparse updating and downdating of the Cholesky factor R k of A T k A k .The orthonormal factor of the QR factorization of the full column rank matrix k .The triangular factor R k recurs in sparse form, and is recomputed from scratch by using the sparse MATLAB QR factorization primitive qr when a refactorization is triggered due to an accumulation of rounding errors.Our toolbox consists of the routines SqrIni (initialization of data structures), SqrIns (addition of a constraint to the working set), SqrDel (deletion of a constraint from the working set), SqrRei (refactorization of A k ), Givens (computation of a Givens rotation) and EjemMSLS (application example).
By using MATLAB notation, we store the working set B in a vector ACT.A call to SqrIni(mode,scal) is used to initialize the use of NaNs as a lower triangular static data structure L, such that L(ACT,ACT) is an upper bound for the transposed Cholesky factor of A(:,ACT)'*A(:,ACT) independent of what ACT is.When mode='dmperm', both the rows and columns of A (and consequently b and c) are reordered for A to be in lower block triangular form (LBTF). Hence, the static structure is tightly set up in L(ACT,ACT).With mode='colmmd', only the columns of A (and consequently b) are reordered for the static structure set up in L(ACT,ACT) to be an upper bound for that Cholesky factor.Scaling techniques can be a priori applied (scal='noscale', by default), like that in which all columns of A and c have a Euclidean norm unit (scal='txtbook'), or like that used in CPLEX, where the columns of A are normalized in ∞ norm, and then the same is done with the rows of A, when scal='cplex92' (see Bixby [41] (p.271)).
Because we are not using the orthogonal factor due to sparsity considerations, the least squares subproblems are solved by using corrected seminormal equations (CSNE) (see e.g., [39] (p.126)).Thus, the descent direction d is computed by solving (4) with CSNE.This requires us to perform an iterative refinement step after having found the solution ȳB of the seminormal equations Then, d = A k ȳB − c is determined as the orthogonal projection of −c onto null(A T k ).This means that d is a null-space steepest-descent direction.The same technique is used to solve the compatible system A k y B = c when c ∈ R(A k ).To test whether v ∈ R n is in R(A k ), we check whether the residual of min A k δ B − v is zero.As a solution x of the undetermined system A T k x = b B we take that of the minimum norm by using corrected seminormal equations of the second kind (e.g., [39] (pp.7, 70, 126)).This requires an iterative refinement step after having found a solution z of the seminormal equations of the second kind, The commercial program TomLab LPSIMPLEX V3.0 [31] enhances a refined MATLAB implementation of the primal simplex method, in which sparse matrices can be handled in the solution process.A sparse QR factorization of the basis matrix is maintained with an explicit orthogonal factor.The updating and downdating is done with customized slight modifications of the dense MATLAB primitives qrinsert and qrdelete, to facilitate dealing with sparse operands.Notice that the authors of the software claim it handles sparsity exploitation better.Default settings were chosen for LPSIMPLEX.

Extreme Cases
The question we have is what options of Sagitta are necessary to solve extreme cases from the literature.This relates to the general question in the investigation of optimization techniques of which type of algorithms are more appropriate to solve which type (defined by their characteristics) of optimization problems.The investigation does not require big instances, but extremely designed cases like the well-known Klee-Minty box.We report on our finding here.
Gass mentioned in [42] that for any new LP code, it is good to find out which instances may work badly.He specifically designed textbook examples.We gathered 27 from a wealth of different sources in a MATLAB routine called ExampSag, which can be found in the Appendix A.
We solved a first set of problems involving primal degenerate (problems 9 and 10) and dual degenerate (problems 14 to 16) problems, corresponding, respectively, to [14] (p.61), [38] (p.31), and [43] (p.136).Sagitta solved all of them without any problem and avoided cycling in the last three cases by using minimum-norm solutions instead of basic ones when solving the underdetermined compatible systems.It is interesting to note that in the latter two cases, cycling appears if we use basic solutions.This illustrates perfectly that the possibility of a cycle cannot be ruled out for some pivoting rules.We looked for cases unfavourable for the EXPAND technique of solving the degeneration of [44].We found them in instances of [45] ( §2, §4), where the Sagitta code appeared effective in obtaining correct results (problems 19 to 22) for those with a nonzero cost vector.
The literature provides instances that are problematic for interior point methods (Mészáros [46] (pp.8, 9) and Stojković and Stanimirović [47] (pp.436, 437)).These are problems 23 to 26 in the Appendix A. Basically, they illustrate that approximate arithmetic may lead to challenges for accepting the optimal solution.
The results of Sagitta are given in Table 4.For the first problem with label 0, the setting 'annns' provides active set {1, 3, 10, 12, 14, 16, 17, 20} with small values for y 10 , y 14 , and y 20 .The setting 'ansns' results in active set {1, 2, 3, 10, 12, 16, 17, 20} with small values for y 2 , y 10 , and y 20 .For the second problem with label 1, Sagitta cannot find the optimal solution with settings 'annns' and 'bnsnt'.Setting 'ansns' gives active set {1, 3, 6, 12, 14, 16, 17, 20, 21}.We focus now on problems designed by Powell [48] to provide challenges for interior point methods with m ≥ 3. We focus on the problem for m ∈ {3, 4, . . ., 500}.We run Sagitta with tolerance √ ε and settings o(1)='b' and o(1)='t' by using o(5)='t' in both cases.For the cases where m is a multiple of 4, the algorithm requires only one iteration.For the cases where m is not a multiple of 4, the setting o(1)='b' requires two iterations.This good behaviour was also predicted by Sherali et al. [49] and Künzi and Tzschach [50], Pan [51], as it fits the feasible set.However, the setting o(1)='t', requires approximately log 2 (mA + B) iterations with B ≈ 10.2 and A ≈ 0.9.Notice that this increases to nine iterations for m = 499.Klee and Minty [52] started a worst-case example for the simplex method, where depending on the pivot rule chosen, the number of iterations is exponential in the dimension n of the problem.We focus on the description given by Paparrizos et al. [53] with three variants.For 0 < ε ≤ 1  3 , the first formulation of the problem is The second formulation follows from reordering the columns of matrix A: The third formulation adds a parameter µ := 1 ε ≥ 3: We introduce a fourth variant (also with parameter µ) in order to obtain the second formulation counterpart with the modified right-hand side.This seems to be a natural generalization, but was not explicitly given by Paparrizos et al. [53]: 1) , i ∈ {1, . . ., n} y ≥ 0.
In the procedure KleMinSp in the Appendix A, we extend to the possibility to convert the problem into standard form by adding dual variables in a negative or positive sense; e.g., KleMinSp(3,10,-3) is the dual of the problem appearing in [43] (p.270) and in [14] (p.61).We varied parameter µ ∈ N and activation criteria.The results of the experiments are given in Table 5.
It is clear that without normalisation of the residuals (Dantzig), several instances show an exponential number of iterations.Normalization provides for all combinations a linear number of iterations of an order between n and 2n.The influence of the position taken by adding the dual slack variables is also interesting.With the setting o='tnnns' for µ ≥ 5, with mode = −3, n + 2 n − 1 iterations are required.Using mode = 3 only requires n.Such differences may disappear with more sophisticated tie-breaking activation criteria.We use for tie-breaking the minimum index criterion (Bland).
Goldfarb [54] designed another deformation of the cube such that steepest-edge rules also require an exponential number of iterations.The parametric problem of Goldfarb [54], with n ≥ 3, β > 2 and δ > 2β, is given by The generation of the problem is provided by the procedure GoldfaSp in the Appendix A. The results obtained by Sagitta without scaling and with 'dmperm' as ordering method, option setting 'tnnnt', and tolerance √ ε do not provide an exponential growth in the dimension of the number of iterations (see Table 6).The computational time is given in centiseconds (CntSecs), which represent hundredths of a second.This illustrates the requirement of careful numerical procedure and the recoveries.As a final extreme case, we consider the work of Clausen [55] to design an exponential LP, with α = 4/5, β = 5/4, and γ = 5.The idea of the designed problem is to have an exponential number of iterations when using the slack basis for both the simplex method for LPs in standard form and the simplex method for LPs in inequality form applied to the dual of the given problem.The procedure ClauseSp in the Appendix A generates the parameters.Let c, A, and b be the data of the original problem max{b T y : Ay ≤ c, y ≥ 0} and let c, Ā, and b be the data of the problem max{ bT ȳ : Ā ȳ = c, ȳ ≥ 0} at the output of the routine; then we can select the following variants: The variants facilitate applying Sagitta to problem min{ cT x : ĀT x ≥ b} in unfavourable circumstances for the simplex method for LP in standard form (if mode=1) and for the simplex method for LP problems in the inequality form (if mode=0).The other modes play the same role with the original problem of minimization instead of maximization, as noticed by Murty [56].
The results obtained by neither scaling nor preordering and with tolerances to √ ε are shown in Table 7 for varying values of o(5) and o(1).We use for the number of iterations The number of iterations for modes 0 and 2 on the one hand, and 1 and 3 on the other hand are very similar, separated by a comma in the last column.The last two experiment rows of the table reflect exponential behaviour.For the first three rows, the behaviour of Sagitta with o(3)=n was surprisingly polynomial, requiring at most n + 2 iterations by using option setting 'tnnnt' in each mode.Preliminary experiments were carried out on classical problems with sparse matrices to compare the Sagitta sparse treatment with the one used by TomLab LpSolve v3.0.Initial problems of interest at the time of the appearance of the simplex method [57] (first published in) in 1947 were diet problems, such as the favorite Stigler 9-variable, 77-constraint [58] (p.3), [59] (pp.551-567) problem.Instead, we focus on the problems by Andrus and Schäferkotter [60] Procedure AndSchSp in the Appendix A generates the problem data given the chosen input parameters.The number of iterations and the computation time (in centiseconds) for different values of n can be found in Table 8.We chose 'dmperm' as the ordering option as it allows A to be banded, and we did not scale 'noscale' the problem, using default options for Sagitta with all tolerances to √ ε.LpSolve needs to perform n + 1 iterations and was unable to solve problems with n = 1000 and n = 2500 after 9 h of execution due to time limits and memory problems.Sagitta uses typically n iterations and took 20 s compared to almost 6 min that LpSolve needed to solve the problem with n = 500.The problems with n = 1000 and n = 2500 were solved in almost a one and a half minutes and just over 13 min, respectively.Note that Sagitta needed n iterations, whereas TomLab needed n + 1.To test it, they generated sparse random testing problems where the constraints are tangential to the unit ball.Consider a matrix A ∈ R n×m with Gaussian-distributed elements generated in MATLAB according to sprandn(n,m,dens,1).Consequently, A is a nonsingular, wellconditioned matrix with random normally distributed elements.Moreover, let index number g be uniformly selected from {1, 2, . . ., m}.In the Appendix A, procedure ChPaSaSp generates a sparse problems of the form max y g , y ∈ R m s. t. ∑ m j=1 A ij y j = 1 , i ∈ {1, . . ., n} y ≥ 0.
As in Chen et al. [62] ( §8), for each of the combinations of interest of number of rows n, number of columns, m and density d, we generate blocks of 20 problems and average the number of required iterations.Moreover, we also average the elapsed time (in hundredths of a second).In Sagitta, we set 'colmmd' as the ordering option, and we do not scale 'noscale' the problem by using default options.
As combinations of interest we first chose the same combinations as Chen et al. [62] §8, and obtained the results shown in Table 10.One can observe that, although Sagitta had a slight advantage in terms of number of iterations, TomLab is clearly faster.This is due to the fact that the problems are not large (n ≤ 150) and above all due to the fact that the maximum skeleton L is not sparse at all, i.e. |L| is between 75% and 85%.To know the behavior for sparser combinations, i.e., |L| does not exceed 10%, we reran the experiments for such instances.The results are given in Table 11.The results of these sparse experiments show that that Sagitta has a slight advantage in terms of number of iterations, but is also faster.The dimension n = 250 is the largest of those solved with TomLab.It took Sagitta 4 s to solve the problem compared to the 83 s required by TomLab.The two largest problems were solved by Sagitta in n iterations requiring 14 and 57 s, respectively.

Discussion
A code has been presented for linear programming to take care of sparsity and basis deficiency.Extreme cases have been taken from the literature, and the performance has been investigated by using a commercial MATLAB solver Tomlab as benchmark.We structured this article as shown in Figure 1 and found that • the code is numerically robust (cf.Section 3.2), despite the fact that it dispenses with the orthogonal factor of the QR factorisation and in no case uses more than one iterative refinement step; • the computational effort in choosing suitable pivoting rules is quite encouraging; and • the code reveals advantages for problems that are really sparse (cf.Section 3.3), both in terms of iterations and computational time, compared to a commercial software with comparable features.
Questions on the simplex method originate from its beginning in 1947.However, there is still a lot of ongoing research; in the last three years one can find many publications, including [1][2][3][4][5][6][7][8], just to cite a few.Although they mainly focus on the simplex method (the second volume of forthcoming renewed edition of Pan's book [8] being the exception), some of the ideas they pointed out could be of interest to be incorporated into a deficient basis framework, thus reducing even more the necessary number of iterations.

else 14 :
Direction d k is initially unit vector e p and elements of B k follow from −δ 15:
§2 Algorithm and benchmark software description §1 Introduction with recent publications and main bibliographical sources §3.1 Numerical and algorithmical settings, along with sparsity-related options for preordering and scaling §3.2 Extreme Cases §3.3 Sparsity Exploitation Small LPs with primal or dual degeneracy and/or problematic for EXPAND, cf.Gill et al. (1989) and Hall & McKinnon (2004) Andrus & Schäferkotter (1996) diet problems (n ≤ 2500), cf.

Table 1 .
Constraint selection criteria if descent direction d exists.
o(1) Meaning 's' Most opposite to d no normalization 't' Most opposite to d with normalization 'a' Among opposite to d, most opposite to (−c) no normalization 'b' Among opposite to d, most opposite to (−c) with normalization 'o' Among opposite to d, most opposite to d + (−c) no normalization 'p' Among opposite to d, most opposite to d + (−c) with normalization

Table 2 .
Activation criteria in absence of descent direction.

Table 8 .
Andrus and Schäferkotter [60]diet problems.Qy ≤ 10 4 e , Q ∈ N n×n y ≥ 0, with e the all-ones vector and q ij ∈ N + randomly chosen in {1, 2, ..., 10 3 }.This causes matrix Q to be dense.The procedure KunQuaSp in the Appendix A generates the data of the problem from a nonsingular sparse Q of density 10 −5 and condition number 20 by means of round(1000*sprand(n,n,1e-5,1/20)), such that |Q| ≈ n.Running the Sagitta routine, we choose 'colmmd' as the ordering option and not to scale 'noscale' the problem and used default option settings.For each value of n, 10 problems were solved; the average number of iterations and the average computation time (in hundredths of a second) for different values of n are reported in Table9.Note that as before, Sagitta requires n iterations whereas TomLab needs n + 1.