# 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)
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

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

1 Like