What is the most efficient way of obtaining the orthogonal column space of a matrix? (many more rows than columns)
Q, = svd( X; alg = LinearAlgebra.QRIteration() )
gives me what I want, but is there something more efficient?
Relatedly, why isn’t there a command similar to nullspace? (or is there?)
Maybe this?
julia> using LinearAlgebra
help?> nullspace
search: nullspace
nullspace(M; atol::Real=0, rtol::Real=atol>0 ? 0 : n*ϵ)
nullspace(M, rtol::Real) = nullspace(M; rtol=rtol) # to be deprecated in Julia 2.0
Computes a basis for the nullspace of M by including the singular vectors
of M whose singular values have magnitudes greater than max(atol,
rtol*σ₁), where σ₁ is M's largest singular value.
By default, the relative tolerance rtol is n*ϵ, where n is the size of
the smallest dimension of M, and ϵ is the eps of the element type of M.
Examples
≡≡≡≡≡≡≡≡≡≡
julia> M = [1 0 0; 0 1 0; 0 0 0]
3×3 Matrix{Int64}:
1 0 0
0 1 0
0 0 0
julia> nullspace(M)
3×1 Matrix{Float64}:
0.0
0.0
1.0
julia> nullspace(M, rtol=3)
3×3 Matrix{Float64}:
0.0 1.0 0.0
1.0 0.0 0.0
0.0 0.0 1.0
julia> nullspace(M, atol=0.95)
3×1 Matrix{Float64}:
0.0
0.0
1.0
Sorry for being unclear. I know about nullspace; I wanted something similar to that to get the column space.
Typically, the QR decomposition uses the Gram-Schmidt process (it computes an orthonormal basis of a set of vectors) to compute Q. So calling Gram-Schmidt would be enough (but computing the R matrix in the QR decomposition really is cheap).
To compute the nullspace, you can apply the QR decomposition to the transposed matrix: linear algebra - Null-space of a rectangular dense matrix - Computational Science Stack Exchange
I should also have understood that had I read you post properly 
qr(X).Q
Probably the reduced QR decomposition is what you want: A = QR, with m \times n matrix Q and n \times n matrix R for m \times n matrix A. The columns of Q span the columns of A. The Householder algorithm for computing the QR decomposition is better than Gram-Schmidt, in that it produces a Q that is unitary to machine precision and results in Ax=b solves that are backwards stable. Julia’s qr algorithm surely uses Householder.
Ah, qr(X).Q [:,1:size(X,2)] is an order of magnitude slower than the svd route, but using @view seems to do the trick, thanks.
It does. What I was missing there was that it also appeared to be computing the null space, but the reason qr(X).Q[:,1:size(X,2)] is so much slower is the subsequent indexing.
QR factors are stored in an implicit form that makes indexing slow, but they can be quickly multiplied by vectors. To get the first n columns, try multiplying Q by the first n columns of the identity matrix.
I think I have a PR that fixes this, but it stalled out. I really should revive it though because it’s a dumb performance cliff.
Probably another reason for the slowness of qr(X).Q is that Householder QR decomp usually does not compute Q explicitly, but rather stores a series of Householder reflector vectors. When you ask for Q, it has to be constructed from the reflector vectors.
EDIT: what @stevengj said!
In particular, try qr(X).Q * Matrix(I, size(X)...), and you should find that it is much faster than an SVD.
Yes, it’s pretty common to want to slice a Q matrix and we really should provide a fast algorithm for it.
Thanks, that’s very helpful and I agree.
Having a function like
colspace(X) = qr(X).Q * Matrix(I, size(X)...)
to complement nullspace would be welcome, also.
You need a more complicated implementation (and probably the SVD, with tolerance parameters similar to nullspace) if you want to handle the case of rank-deficient matrices.