# Nullspace algorithm for a sparse matrix

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?

• I guess you are trying to achieve this without any open source package or etc right? Feb 6, 2014 at 11:33
• Well, I am, but if there is anything out there that can help, I am open to that. There's another requirement I didn't mention though: I need to do the computation exactly. The elements of the matrix are all integers, so if I were to use a package, it would have to support computation with integers and thus rational output. Feb 6, 2014 at 11:51
• You should definitely revise the question to to add the parts about integer arithmetic! This will drive the responses more than anything else about it. Also, you should quantify large and very sparse. Feb 6, 2014 at 15:26
• Okay, thanks. I edited it in. I will be dealing with rational matrices in the future, so I listed that instead of integers. There shouldn't really be any diffeernce though. Feb 6, 2014 at 16:36

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.

• Thank you and sorry for the late reply. I am very unfamilliar with many things you mention, but I will look into them as much as I can. Feb 10, 2014 at 13:06
• Ok, feel free to ask if you need me to elaborate more, I'm happy to help! Feb 10, 2014 at 17:15
• Before trying to reduce fill-in, I'm first trying to find an algorithm that can find integer nullsapces. It turns out that most algorithms fail because they encounter very large numbers very quickly. A solution exists, but I think it might require more storage than just storing my initial matrix regularly. I think this is a separate question though, so I will post a new one. Feb 12, 2014 at 14:49