Exact Computation of Basic Solutions for Linear Programming
Exact linear programming has many applications such as finding provably
solutions in standard application areas, computer assisted proofs and subroutines
for exact integer programming. The recent proof of the Kepler Conjecture by
Thomas Hales relies on solving thousands of linear programs exactly.
The fastest known technique for solving linear programs exactly is performing
simplex algorithm with floating point arithmetic , then computing, verifying and
correcting the final solution with symbolic computation. This leaves solving sparse
rational systems of equations as the major bottleneck. We investigate methods to
solve sparse rational systems of linear equations exactly, and study their behavior
computationally on a large test set arising from linear programming problems.
Exact Linear Programming
Linear programming is a technique for
optimizing a linear objective function
over a set of linear inequalities. Geometrically,
the inequalities represent a convex
Standard form LP
The most commonly used technique for linear programming is the simplex
A simplified geometric description is:
•Start at a vertex of the polyhedron
•Follow edges of the polyhedron to better vertices
•Terminate at locally optimal vertex
At each vertex of the polyhedron the corresponding solution vector x can be
by solving a system of linear equations. This is called a basic solution, and
the subsystem of A defining this solution is called a basis.
If the simplex method is performed using floating point computation, the
the final solution provides a structural description of the solution as the solution to
a system of equations. The coefficients of these equations come from the original
problem description, so the final solution can be computed exactly.
Linear programming duality theory provides a tool to certify optimality of an
solution. Every linear programming problem has a related dual problem and if
either is feasible they share the same optimal objective value.
Linear Program (Primal)
If x and y are feasible, then they are optimal if and only if cTx = bTy. An
basis of the primal problem (vertex of the polyhedron) structurally describes both
the primal optimal solution, and the dual optimal solution y. Given the optimal
basis, calculating and certifying the exact optimal solution involves solving two
systems of linear equations.
QSopt ex: A Rational LP Solver
The QSopt ex software  (by Applegate, et al.) does the following.
•Apply simplex in floating point arithmetic
•Compute basic solution exactly from final basis
•If solution is not optimal, increase precision
•Repeat procedure in higher precision, starting at current basis
This strategy proves very effective. Much of the simplex algorithm is done in
point. Even when a suboptimal basis is returned due to floating point errors,
it is often close to the optimal basis, reducing the extended precision computation .
The final output includes the optimal solution and an exact certificate of optimality.
Solving Sparse Rational Systems
Several methods are known for solving sparse systems of rational equations Ax
exactly, beyond applying standard techniques using rational arithmetic [3, 4]. If an
accurate enough approximate solution is calculated, the exact rational solution can
be reconstructed using diophantine approximation. If the problem is scaled to have
integer coefficients and solved modulo a large enough integer, the exact rational solution
can be calculated using rational number reconstruction . Diophantine approximation
and rational number reconstruction are performed by calculating continued
fraction convergents using the Extended Euclidean Algorithm.
Exact LU Decomposition
•Solve system directly using an LU decomposition over Q
Sparse decomposition algorithms focus on permuting rows and columns of the
to avoid fill-in and arithmetic operations .
•Compute floating point LU decomposition of A
•Iteratively refine approximate solution
•Use rational number reconstruction to find exact solution
Using mixed precision iterative refinement an approximate solution is
suitable for rational reconstruction/diophantine approximation to be applied.
•Use Wiedemann’s method to solve Ax = b over finite field
•Use p-adic lifting to find solution modulo large number
•Use rational number reconstruction to find exact solution
Wiedemann’s method  for solving systems of equations over a finite field
black box algorithm, it only accesses the matrix for matrix-vector multiplication.
Compute LU decomposition of A over Zp
Use p-adic lifting to find solution modulo large number
Use rational number reconstruction to find exact solution
The finite field solves are done using a sparse LU
factorization, designed to minimize
We implemented the four methods in C. We rely on the GNU
Library to store rational numbers and large integers. The codes for LU factorization
over the rational numbers, finite fields and double precision floating point numbers
all share a common structure and were adapted from the QSopt software, which is
designed for very sparse matrices. The rational reconstruction codes employ output
sensitive lifting, Lehmer’s algorithm and other techniques to increase speed.
The linear-programming research community is fortunate to
have several publicly-available
libraries of test instances. We collected instances from several test libraries,
including NETLIB, MIPLIB 3.0 and many others. These collections are
comprised of instances gathered from business and industrial applications, and from
academic studies. The instances were were given to QSopt ex and the optimal basis
matrix, or the last basis encountered after 24 hours was recorded. After some
preprocessing and reduction the final set contains 276 instances, with dimensions
ranging from 100 to over 50,000. For each instance we have both the square basis
matrix and the corresponding right-hand-side vector.
We have given a performance profile comparing the
relationships. A performance
profile plots the number of instances solved within a factor x of the fastest solve
time among the methods. The vertical axis represents the number of instances.
The horizontal axis gives the solve-time ratios. Ratios are compared instead of
averaging solution times , which ranged from hours to seconds.
This table gives the solve time ratios on subsets of the
problems partitioned by dimension,
solution size and density (nonzero entries per row). The solution bitsize
is the largest number of bits required to represent a component of the solution vector.
We consider the ratio with respect to the Dixon method, and give the geometric
average of the ratios. This table allows us to observe the effect of various problem
characteristics on the relative performance of the methods.
We evaluated several methods for symbolically solving
large sparse rational systems
of linear equations arising from a wide range of linear programming problems.
Dixon’s method and iterative refinement had similar performance, beating
the other methods. Iterative refinement was slightly faster in our experience, but
Dixon’s method avoids numerical difficulties by using finite field arithmetic. We
found Wiedemann’s algorithm to be much slower than the other methods tested on
this problem set. This set of real world problems also provides a useful test set for
future studies. Further computational results and analysis are contained in .