Principal angles between subspaces

Does Julia have a function to compute canonical / principal angles between subspaces, like scipy.linalg.subspace_angles or even just Matlab’s subspace which computes the maximum angle only?

The computation is essentially acos.(svd(qr(A).Q' * qr(B).Q)), but some additional tricks are required so that it is accurate when the angles are small.

Untested, does it seem accurate?

module AngleBetweenSubspaces

import LinearAlgebra

svd_u(m::AbstractMatrix) = LinearAlgebra.svd(m).U

The `rank` argument should be the rank of `m`.
range_orthonormal_basis(m::AbstractMatrix, rank::Int) = svd_u(m)[:, begin:(begin + rank - 1)]

sv_min(m::AbstractMatrix) = last(LinearAlgebra.svdvals(m))
sv_max(m::AbstractMatrix) = first(LinearAlgebra.svdvals(m))

function largest_principal_angle_unchecked(
  a::AbstractMatrix, b::AbstractMatrix,
  rank_a::Int, rank_b::Int,
  # Inspired by the Octave `subspace` code and Andrew Knyazev's Matlab `subspace` code.

  A = range_orthonormal_basis(a, rank_a)
  B = range_orthonormal_basis(b, rank_b)

  c = A' * B

  cos = clamp(sv_min(c), false, true)

  if true < 2*cos^2
    d = (
      (size(A, 2) < size(B, 2)) ?
      A - B * c' :
      B - A * c
    sin = clamp(sv_max(d), false, true)

The largest principal angle between the ranges of `a` and `b`.

The matrix `a` has rank `rank_a`, similarly with `b`.
function largest_principal_angle(
  a::AbstractMatrix, b::AbstractMatrix,
  rank_a::Int = LinearAlgebra.rank(a),
  rank_b::Int = LinearAlgebra.rank(b),
  (size(a, 1) == size(b, 1)) || error("the row counts must match for the two matrices")
  largest_principal_angle_unchecked(a, b, rank_a, rank_b)


Usage example:

julia> include("/home/nsajko/AngleBetweenSubspaces.jl")

julia> const ang = AngleBetweenSubspaces.largest_principal_angle
largest_principal_angle (generic function with 3 methods)

julia> ang(rand(7, 2), rand(7, 7))

Regarding my implementation, if one doesn’t provide the ranks of the two matrices as parameters, there are two LinearAlgebra.rank calls that slow things down considerably. LinearAlgebra.rank in fact calls LinearAlgebra.svd with the same parameters as in the main part of the code, so this additional work is actually duplicated. The reason I did it like that was to avoid thinking about what’s a good tolerance for determining the rank of a matrix.


Thanks, great work! If it matches Knyazev’s implementation then that should be the state-of-the art algorithm.

1 Like