julia> using LinearAlgebra
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.
julia> M = [1 0 0; 0 1 0; 0 0 0]
1 0 0
0 1 0
0 0 0
julia> nullspace(M, rtol=3)
0.0 1.0 0.0
1.0 0.0 0.0
0.0 0.0 1.0
julia> nullspace(M, atol=0.95)
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.
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.