How to find the linearly independent columns (rows) of a matrix

General Usage
14

Yes, I think you can use something like:

using LinearAlgebra
function independent_cols(A; rtol=size(A,2)*eps(eltype(A)))
    QR = qr(A, ColumnNorm())
    atol = rtol * abs(QR.R[1,1])
    # just use QR as alternative to rank(A) (which computes the SVD):
    rank = something(findfirst(i -> abs(QR.R[i,i]) <= atol, 1:size(A,2)), size(A,2)+1) - 1
    return QR.p[1:rank]
end

to get a list of the independent columns up to some tolerance.

PS. For the algorithm in the linked paper, note that the expression R_1 = C_1^T (C_1 C_1^T)^{-1} is a badly conditioned way to compute the pseudo-inverse of C_1 — you would be generally better off using the QR factorization (of C^T or C_1^T, and you already have the former from computing the independent rows of C so you might as well re-use it) or the SVD for this. In general, the formulations in that paper appear oriented towards exact arithmetic and seem like they would benefit from some re-formulation using SVD or pivoted QR.

PPS. Seems like rank(QR) should be implemented. I filed an issue: rank(::QRPivoted) method? · Issue #53214 · JuliaLang/julia · GitHub

3 Likes