Nullspace algorithm for a sparse matrix

Question

I am dealing with large, sparse, rational matrices that I need to determine the nullspace of. Currently, I have one that is about 12000x12000 (but not square), where one in every 2000ish elements is nonzero. Subsequent matrices will grow both in size and sparsity.

Most nullspace algorithms rely on bringing the matrix in (some form of) a reduced row-Echelon form, and then reading the nullspace vectors off of that form. Since these algorithms require large amounts of elementary row operations, I would think that using Compressed Row Storage is the best option, since it allows easy acces to full rows, which is what is needed for the row operations.

However, as a part of these operations, it is required to insert elements in the sparse matrix occasionally. Now, as I understand it, inserting elements in a CRS matrix becomes increasingly expensive when the size of the matrix grows. Are there other options for the storage or the nullspace algorithm that avoids this problem?

Answer 1

tl;dr : Use list-of-lists format, reorder your matrix, use sympy.

You can use a different sparse matrix storage scheme like the list-of-lists format, which allows for fast insertion of new matrix entries but suffers from poor cache coherency. You can then convert the inefficient storage scheme back to CRS once you're done finding the matrix's echelon form.

Permuting the entries of the matrix can make a huge difference by reducing the amount of fill-in that you encounter. The Cuthill-McKee ordering is fairly common and works by doing a breadth-first search of the underlying sparse graph, reordering the matrix in the order that the unknowns were visited. This tends to cluster all the unknowns closer to the diagonal of the matrix. Alternatively, you can use the approximate minimum degree algorithm, which you can play around with in MATLAB using the function symamd. The documentation for the deal.II library has a really awesome comparison of a mess of different reordering schemes but I can't seem to find it.

Ascertaining the connectivity structure of the reduced echelon form of the matrix is basically a problem of graph theory for which you needn't know the actual matrix entries (see Davis's book on sparse direct methods). So you can first find which entries are non-zero in the echelon form of the matrix, and only then perform the computations, content in the knowledge that you won't have to add any new non-zero entries.

To my knowledge, the more common computational libraries don't have support for exact rational arithmetic, only floating point. If you like Python, the sympy package can do exact rational arithmetic and has some support for sparse matrices.