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.

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.