Best way to exactly solve a linear system (300x300) with integer coefficients

Question

I want to solve system of linear equations of size 300x300 (300 variables & 300 equations) exactly (not with floating point, aka dgesl, but with fractions in answer). All coefficients in this system are integer (e.g. 32 bit), there is only 1 solution to it. There are non-zero constant terms in right column (b).

A*x = b

• where A is matrix of coefficients, b is vector of constant terms.
• answer is x vector, given in rational numbers (fractions of pairs of very long integers).

The matrix A is dense (general case), but it can have up to 60 % of zero coefficients.

What is the best way (fastest algorithm for x86/x86_64) to solve such system?

Thanks!

PS typical answer of such systems have integers in fraction up to 50-80 digits long, so, please don't suggest anything based only on float/doubles. They have no needed precision.

Answer 1

The Dixon's algorithm is. It does a inversion of the matrix in modular arithmetic, based on some prime modulus P (this step is easy and exact). Then it solves the trivial system for each x to get a p-adic number representation of the solution. It can be converted to a rational number. Method was described in "Dixon J. D. Exact solution of linear equations using P-adic expansions // Numer. Math., 40, 1982, P. 137-141.", but this article is offline.

There is a article about Dixon's algorithm and its modifications online: http://www2.isye.gatech.edu/~dsteffy/papers/OSLifting.pdf

Author also have all described algorithms implemented and published at http://www2.isye.gatech.edu/~dsteffy/rational/

Another interesting approach is http://en.wikipedia.org/wiki/Montante's_method>Montante's method, which is a gauss variant for integer coefficient matrices. This variant don't create rational numbers in the matrix while converting it to triangular. Take a look at the last picture in link http://en.wikipedia.org/wiki/Montante's_method

Answer 2

Don't know if it's the best, but Scipy's linalg.solve works. You could look at the source code...

EDIT for the cases that I tried scipy.linalg.solve(A,b) * scipy.linalg.det(A) gives an exact result, e.g. gives a matrix corresponding to (1/det(A)) * x, so the actual fractional result is each element divided by that common denominator.

Answer 3

For the algorithm, straightforward Gaussian elimination should be fine. This part may actually be simpler than when dealing with floating-point numbers, since you don't have to worry about the numerical stability.

Depending on your matrix, you might be able to do a bit better that Gaussian elimination, e.g. if your matrix is symmetric, you can use a Cholesky factorization. But 60% sparseness isn't that much in the scheme of things. If it were tridiagonal, banded, etc. then you could try some specialized methods.

You will need a good rational number library to handle the actual arithmetic operations. GNU MP should fit the bill, and claims it can do rational numbers, but I've never used it.

Answer 4

This kind of thing can be done by symbolic computation packages like Maple, Mathematica, Maxima, etc.

If you need a subroutine to do this that you'll call from a larger program, then the answer will depend a lot on the programming language that you're using and what kinds of licensing you're willing to work with.

Answer 5

This was meant to be comment -- in Mathematica it takes me 0.5 seconds for 300x300 matrix with around 60% zeros and 40% positive integers uniformly sampled from 1-10.

To integrate Mathematica library with C, you need to use MathLink interface. I'm not sure if you can statically link it, it may require launching Mathematica's kernel along with your program in order to use the interface.

The relevant function is LinearSolve. Here's the complete code I used to run timing comparison

mat1 = RandomVariate[BernoulliDistribution[.4], {300, 300}];
mat2 = RandomInteger[{1, 10}, {300, 300}];
mat3 = mat1*mat2;
answer = RandomInteger[{1, 10}, 300];
result = LinearSolve[mat3, answer]

I don't know the algorithm they used. They say version 8 is about 1000 times faster that previous version, so it's probably something specialized and proprietary

Answer 6

Use Gaussian ellimination with only integers. Don't do the division. Instead to zero a paritucar cell in the matrix use euclid algorithm. You may need bignum library.

This is almost equivalent to using rational number gaussian ellimination but simpler.

Answer 7

The Gaussian elimination algorithm is simple to program.

Define a rational number object which can be a list of numerator factors and a list of denominator factors. Write functions for doing the four arithmetic operations using the rational number objects. In your Gaussian elimination program, substitute an array of the rational number objects for the array of integers, and use your functions for the calculations.

You may or may not want the number objects to automatically reduce to simplest terms after a computation.

I don't think it would be fast, but it's pretty straight forward, and I don't think the execution time would be long enough to be a bother, maybe a few seconds.