  Updating a linear least squares solution can be critical for near real-time
signalprocessing applications. The Greville algorithm proposes a simple formula
for updating the pseudoinverse of a matrix A $\in$ R nxm with rank r. In this
paper, we explicitly derive a similar formula by maintaining a general rank
factorization, which we call rank-Greville. Based on this formula, we
implemented a recursive least squares algorithm exploiting the rank-deficiency
of A, achieving the update of the minimum-norm least-squares solution in O(mr)
operations and, therefore, solving the linear least-squares problem from
scratch in O(nmr) operations. We empirically confirmed that this algorithm
displays a better asymptotic time complexity than LAPACK solvers for
rank-deficient matrices. The numerical stability of rank-Greville was found to
be comparable to Cholesky-based solvers. Nonetheless, our implementation
supports exact numerical representations of rationals, due to its remarkable
algebraic simplicity.

Mathematical Software