Julia - Find basis of column space of matrix

Viewed 310

I have a (n, m) matrix A which has rank r < minimum([n, m]). I would like to find a basis of r vectors spanning the column/row space. How can I do that?

Here's a how one could generate the data. Since the rank is 3 I would expect to be able to find a basis of 3 vectors.

using LinearAlgebra: svd, rank, Diagonal

function generate_rank_r_matrix(r; n=20, m=10)
    U, S, V = svd(randn(n, m))
    return U[:, 1:r] * Diagonal(S[1:r]) * V[:, 1:r]'
end

A = generate_rank_r_matrix(3)
rank(A) 

Attempt

I know that the Q matrix of the qr decomposition essentially contains some basis for the column space but I am unsure about how to find the r vectors only.

Perhaps one could use the SVD?

2 Answers

RowEchelon is right what you need! All you need is drop zero rows after calculating row echelon form

You could definitely use the SVD. But the QR decomposition is generally cheaper.

For your example, the result of the QR decomposition has two attributes Q and R and the diagonal elements of R tell you how many elements of Q are interesting as your desired column basis. You may need to account for permutations if the decomposition used fancy pivoting. Sparse matrices will almost certainly make use of that to avoid in-fill.

julia> x = qr(A);

julia> x.R
10×10 Matrix{Float64}:
 -3.74988  -1.84617  -0.690789   0.572068     …  -2.25879      -1.66442
  0.0      -1.35248  -0.951006  -0.424635        -1.81135       1.73008
  0.0       0.0      -1.01584    1.49595         -3.20437       0.880599
  0.0       0.0       0.0       -8.03415e-16      1.38778e-15  -1.11022e-16
  0.0       0.0       0.0        0.0             -3.88578e-16   2.63678e-16
  0.0       0.0       0.0        0.0          …  -2.91378e-16   4.50689e-16
  0.0       0.0       0.0        0.0             -4.09934e-17   8.76055e-17
  0.0       0.0       0.0        0.0             -2.66033e-16   1.40809e-17
  0.0       0.0       0.0        0.0             -6.19269e-16   1.2146e-16
  0.0       0.0       0.0        0.0              0.0           3.00838e-16
julia> diag(x.R)
10-element Vector{Float64}:
 -3.7498805345096264
 -1.3524838792645477
 -1.0158363404176274
 -8.034148497889382e-16
 -4.903462616556098e-16
 -9.260239239557565e-16
  1.0342057190512138e-15
 -2.5223684462590304e-16
 -6.192691559225543e-16
  3.008378301264802e-16

Please note that all of this information is available in a very readable form in the documentation. To find that, type '?' in the REPL and then qr.

Related