Sparse and Exact Linear Programming

Introduction and Motivation

Many of todays linear programming solvers rely on floating point arithmetic. State-of-the-art LP solvers --- such as CPLEX --- do not guarantee the correctness of their results. Almost always the answers are correct, but it might happen, that your solver finds a solution that is not optimal or not even feasible. Maybe it finds a solution, although the problem is infeasible or unbounded. Thus, how can you be sure that an answer is the correct solution to your problem?

Revised Simplex Algorithm

There are several techniques to solve a given LP; our approach bases on the Revised Simplex Algorithm. Let a linear programming problem

maximize c^Tx
subject to Ax = b
x >= 0

be given in standard form. We use standard pivoting techniques to solve that problem. At every step in the pivoting algorithm we only need to know the current basis inverse matrix A_B^-1 - and of course the given matrix A and the vectors b and c. Once we know the basis index set B, we can inverse the basis matrix and calculate all necessary data, such as current objective value, current basis solution or current dual solution, as well as the whole or a part of the simplex tableau. But of course we will not invert a matrix every round during the Simplex algorithm. Instead we decompose the basis matrix A_.B into a lower triangular matrix L and an upper triangular matrix U.

During one pivoting round we only exchange one basis variable with one non-basis variable. This means for our decomposition matrix U that only one column is changing. Exploiting this and the fact that we can order the basis and non-basis variables as we like, we can reorder the rows and columns of the matrices having as less work as possible for the updating procedure.

Sparsity and Exactness

If the constraint matrix A is sparse, the chances to do further required row- and column permutation in the matrices L and U arises to avoid divisions by zero, but additional computational efficiency can be obtained by making further row- and/or column-permutations too. The goal is to minimize the number of single calculations in the entire algorithm; especially multiplications with an exact number type (see below) are expensive. Therefore keep all matrices (and maybe vectors) as sparse as possible to minimize the cost for the computation of parts of the simplex tableau and of basis and non-basis values.

In general it is not easy to find the best decomposition strategy keeping the matrices sparse. There are heuristics to minimize the fill-in, but these methods operate only locally. Indeed, it has been shown, that the problem of finding the optimal strategy to decompose a matrix into sparse matrices, is NP-complete.

But there is one more problem: We want to solve LP-problems exactly. Therefore inexact floating point computations should be avoided. As usual in complexity theory, we may assume that all input data is integer or equivalently rational. A possible approach is to use rational numbers, which guarantee the correctness of the result; but this method has tremendous cost in time and space. A multiplication is not longer a constant time operation, but depends on the length of the involved numbers. To avoid exponential growth of intermediate numbers we can use the Euclidean Algorithm to normalize a quotient. But this results in an algorithm in which not even one pivot step can be calculated in polynomial time.

This problem can be solved by using either Q-pivoting or by applying Sylvester's Identity to develop an integral LU-decomposition. These techniques execute fraction free divisions at some point in the procedure; and you can prove, that your numbers grow only polynomially. With this new technique the intermediate matrices and vectors are still integral.

Decomposition Strategies

The goals of a decomposition strategy for the new method are slightly different than for the common LU-decomposition. Numerical stability is no longer an issue. Instead, keeping the size of the matrix elements as small as possible and trying to keep the matrices as sparse as possible is important.

We developed and adapted some decomposition strategies; some only concerning sparsity, some controlling the growth of the numbers. Most of the developed strategies are greedy; figuratively speaking they look only one step into the future during the decomposition procedure. We looked at several objectives, such as minimizing the fill-in or minimizing the maximum number in one decomposition step.

To get rid of the locality we invented a strategy, that chooses the decomposition procedure for which the next k-1 decomposition steps yield the sparsest possible matrix. But as results showed, only for very small k and for not too dense matrices, this k-Lookahead Search technique is applicable.

Last but not least we also implemented a mixed strategy, that uses the slower but more accurate k-Lookahead Search as long as the matrix is not too dense and then switches to a fast greedy heuristic. This approach was also motivated by the result that the beginning of the decomposition procedure is very important. If you do bad work here, it propagates to the end. So the little extra work investing in the beginning of the procedure is worthwhile compared with the result, since in the Simplex Algorithm the density of the matrix is really important.

A random strategy serves as a benchmark for the other strategies.

Results

Surprisingly the size of the maximum element is not an issue. All strategies - also the random strategy - yield more or less the same size of the maximum element. We showed, that the size of the maximum element increases by 8 bits per decomposition step.

All decomposition methods have been implemented and tested against several objectives: Sparsity (or density), maximum size of an element and cost for each decomposition step. As reference we set up the random strategy. In the following picture you can see the evolution of the work cost during the decomposition procedure for a 150x150 random generated sparse and regular matrix. We measured the number of operations weighted by the length of the numbers. The overall cost includes the strategy cost and the effective procedure cost. Not very surprisingly the Lookahead method is very expensive, because of its strategy cost. But as you can see the difference to the others is not very big in the first half. On the other hand the random strategy's procedure cost outweighs the strategy cost (which is negligible).

See also ...

... Term Project by Dominic Jost