Sample Nullspace Using Colt

Question

I'm writing Java and using colt as my matrix library and would like to find a (any) vector in the kernel of a matrix. I can do this in python using sympy as follows:

def kernel(A, n):
    if A.rows == 0:
        return Matrix([1]*n)

    R, pivots = A.rref()
    Ap = A.extract(range(A.rows), pivots)
    bp = Matrix([0]*Ap.rows)

    free = list(set(range(n)) - set(pivots))
    for i in free:
        bp -= A[:, i]

    xp = Ap.LUsolve(bp)
    x = [1]*n

    for i in range(len(pivots)):
        x[pivots[i]] = xp[i]

    return Matrix(x)

Using sympy I can call nullspace to get the entire nullspace or use rref to get the pivots used when reducing to row-echelon form and from that find a single vector in the nullspace myself. I cant find a function in Colt to calculate the nullspace and trapezoidalLower doesn't return the pivots.

Am I left to write my own rref or does someone know a higher level way of achieving this with Colt?

Answer 1

The answer is WHATEVER YOU DO DONT USE RREF in java. Converting to reduced echelon form turns out to have lots of comparisons to 0. If the value is 0 we do one thing. If the value is very close to 0, but not quite 0, we do something completely different (like divide by the value). This means one unstable algorithm.

Instead we can use QR Decomposition, which happens to be implemented in colt.