Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I'm looking for a library that performs matrix operations on large sparse matrices w/o sacrificing numerical stability. Matrices will be 1000+ by 1000+ and values of the matrix will be between 0 and 1000. I will be performing the index calculus algorithm so I will be generating (sparse) row vectors of the matrix serially. As I develop each row, I will need to test for linear independence. Once I fill my matrix with the desired number of linearly independent vectors, I will then need to transform the matrix into reduced row echelon form.

The problem now is that my implementation uses Gaussian elimination to determine linear independence (ensuring row echelon form once all my row vectors have been found). However, given the density and size of the matrix, this means the entries in each new row become exponentially larger over time, as the lcm of the leading entries must be found in order to perform cancellation. Finding the reduced form of the matrix further exacerbates the problem.

So my question is, is there an algorithm, or better yet an implementation, that can test linear independence and solve the reduced row echelon form while keeping the entries as small as possible? An efficient test for linear independence is especially important since in the index calculus algorithm it is performed by far the most.

share|improve this question
up vote 4 down vote accepted

You can work modulo a number of large primes to get the results modulo these primes, then check whether there are rationals with few enough digits satisfying these congruences. If yes, you can check by a matrix-vector multiply whether the approximation found is exact. This can be turned into an exact decision algorithm.

However, if the determiant of the matrix has a size of the order of $10^{1000}$ (quite possible in your scenario), you'll generically not find solutions whose components need less that a few thousand digits.

related links:

share|improve this answer
Thanks, the first paper looks promising! – jgonagle Oct 21 '12 at 21:05

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.