Principal angles between subspaces

Specific Domains Numerics
1

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.

2

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)
    asin(sin)
  else
    acos(cos)
  end 
end

"""
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)
end

end

Usage example:

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

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

julia> ang(rand(7, 2), rand(7, 7))
8.407572180002955e-16

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.

2 Likes
3

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

1 Like