Finding the Kernel of a Matrix

I'm working on an optimization problem (coding in Java) in which, should various celestial bodies align the wrong way, I may need to compute the rank of a real matrix $A$ and, if it's less than full rank, a basis for its kernel. (Actually, I could get by with just one nonzero vector in the kernel, but I'm greedy.)

So I spent a couple of days doing Google searches to see which open-source linear algebra libraries do what. My ideal library would be easy to install (no compiling from source), would support sparse matrices (which mine will be), would make it easy to find the kernel (which turned out not to be a given), and would run fast (without using my computer's GPU, which some of the libraries do). My search led me to install and test several libraries, only to discover that some did not support sparse matrices and some did certain things in rather peculiar ways, designed to make it hard if not impossible to drill down to the basis of the kernel. (One particularly annoying library threw an exception because the matrix I had just constructed, in the immediate preceding line using a sparse matrix constructor, was not considered sparse.) I found something that works, and I thought I'd document it here. If I find something better in the future, I'll post that as well.

Before proceeding, let me point out the Java Matrix Benchmark, which provides useful benchmarking information (and links to) a number of linear algebra packages.

What works for me, at least for now, involves the Apache Commons Mathematics library. This is one of the more commonly used (no pun intended) mathematics libraries in Java-land. Commons Math supports sparse matrices, at least for storage. (I'm not sure if computational operations, other than add/subtract/multiply, exploit sparsity.) It also does both QR and SVD decompositions. A number of responses on Q&A sites suggested using either QR (faster) or SVD (more numerically stable) decomposition of the matrix $A$ to get to its kernel. I opted for a QR decomposition. As I found out the hard way, though, not all QR decompositions are created equal.

Cutting to the chase scene, the key is to do a "rank-revealing" QR decomposition, which means using the RRQRDecomposition class, not the QRDecomposition class. What you decompose is actually $A^T$, the transpose of $A$. So if $A$ is an $m \times n$ matrix, the decomposition looks like$$A^T P = Q R,$$where

• $P$ is an $m \times m$ pivot matrix,
• $Q$ is an $n \times n$ orthonormal matrix (i.e., $Q^T Q = I$), and
• $R$ is an $n \times m$ upper triangular matrix.

If $A$ has rank $r$, the last $n - r$ columns of $Q$ provide a basis for the kernel of $A$. (If $r = n$, $A$ is full column rank and the kernel is just $\{0\}$.)

I wrote a little test program (one short Java file) to make sure I was doing things correctly. It generates a random matrix, decomposes it, and confirms that the last however many columns of $Q$ really belong to the kernel of the matrix. If you want to see things in action, you can get the code from the blog's GitLab repository. You'll need to have a recent version of the Commons Math library (I used 3.6.1) on your class path. There are various parameters you can play with: a random seed; the dimensions of $A$; how dense $A$ should be; a rounding tolerance (how close to 0 counts as 0); and a flag which, if set true, tells the matrix generator to replace one column of $A$ with a random linear combination of the others (just to ensure that $A$ does not have full column rank).