What is the most efficient way of obtaining the orthogonal column space of a matrix?

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

2 Likes

I should also have understood that had I read you post properly :sweat_smile:

1 Like

qr(X).Q

2 Likes

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.

1 Like

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.

2 Likes

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.

2 Likes

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!

1 Like

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.

6 Likes

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.

1 Like

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.

2 Likes