|
|
|
|
|
|
Numerous applications in optimization require solving large scale linear programs (LPs), for which the Simplex Algorithm is one of the most widely used solution methods. Each step of this algorithm requires solving (at least) two large, very sparse and nonsymmetric linear equation systems, one with the current basis matrix and one with its transpose. In most existing computer codes these solves are performed directly, based on an LU factorization of the basis matrix.
We evaluated the applicability of modern Krylov subspace methods for the linear systems arising in the Simplex algorithm. We concluded that in this LP context these methods do not offer a competitive alternative to the established direct methods. To get an impression why, consider figure (F1) below: it depicts the spectrum of an arbitrarily chosen basis matrix of the problem momentum1 from the MIPLIB. The eigenvalues form a cluster at zero and enclose the origin circularly, so that one can expect rather unfavorable convergence behavior of Krylov subspace methods. In our report we argue that even if we had a very good preconditioner at hand, direct methods would still be the method of choice. Although this particular basis matrix was picked arbitrarily, the situation is typical for most LPs we have investigated.
|
|
| (F1) | (F2) |
We then focussed on the analysis of the existing linear algebra kernel based on direct methods. For our experiments, we used the SoPlex LP solver from ZIB, which implements a commonly used algorithmic framework for the kernel. More precisely, the linear algebra kernel of such simplex based solvers has three components: routines for the LU factorization, routines for the solution of sparse triangular systems and an LU update mechanism, that allows for reuse of an existing LU factorization over many simplex iterations. The overall time spent in the linear algebra kernel typically accounts for the major portion of actual run time. In particular, in order to minimize the time spent for the solution of the triangular systems stemming from the factorization and to maximize the effectiveness of the LU update procedure, it is crucial to minimize the number of fill-in elements created during the factorization process. To reduce the fill-in, most existing LP codes use the Markowitz pivoting heuristic, which origins back to the 1950's. We demonstrate numerically that for LP basis matrices, despite the immense improvements of more recently developed factorization techniques, Markowitz pivoting produces LU factors close to optimal in terms of fill-in and consistently outperforms modern LU codes.
Figure (F2) above gives an impression of the effectiveness of the Markowitz heuristic for LP basis matrices: Every LP was solved to optimality (or until a certain time limit was hit). Whenever a factorization of the basis matrix was computed, the ratio of number of nonzero elements in the factorization and the number of nonzero elements in the basis matrix was recorded. These ratios were averaged over all factorizations performed within the run, and sorted into "bins". So (F2) shows an histogram of the average relative fill in for all LP's of our test set. It is immediately clear that the vast majority of the LP's suffer very few fill-in (for most LP's, the ratio is less than 1.25). As we explain in our report, the reason for this is twofold: First, many LP basis matrices contain a large triangular part which doesn't need factorization at all. Secondly, the remaining part is typically quite small, completely unstructured and very sparse. In contrast to modern LU techniques, Markowitz performs pivoting for sparsity as much as possible, thereby deemphasizing aspects like locality of data access, potential parallelization and numerical stability. Our experiments thus confirm the common perception of the superiority of Markowitz pivoting over other techniques in the LP context.