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

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).

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.

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.

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.

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.