# Principal angles between subspaces

**URL:** <https://discourse.julialang.org/t/principal-angles-between-subspaces/101044>\
**Category:** Numerics\
**Tags:** linear-algebra\
**Created:** [July 1, 2023, 7:00am UTC](https://discourse.julialang.org/t/principal-angles-between-subspaces/101044 "2023-07-01T07:00:56Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![fph](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fph/32/17159_2.png) [@fph](https://discourse.julialang.org/u/fph)\
**Post date:** [July 1, 2023, 7:00am UTC](https://discourse.julialang.org/t/principal-angles-between-subspaces/101044/1 "2023-07-01T07:00:56Z")

</div>

Does Julia have a function to compute canonical / principal angles between subspaces, like [`scipy.linalg.subspace_angles`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.subspace_angles.html) or even just Matlab’s [`subspace`](https://it.mathworks.com/help/matlab/ref/subspace.html) 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.

---

<div class="post-metadata">

**Author:** ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)\
**Post date:** [July 1, 2023, 12:42pm UTC](https://discourse.julialang.org/t/principal-angles-between-subspaces/101044/2 "2023-07-01T12:42:27Z")

</div>

Untested, does it seem accurate?

```julia
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-repl
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.

---

<div class="post-metadata">

**Author:** ![fph](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fph/32/17159_2.png) [@fph](https://discourse.julialang.org/u/fph)\
**Post date:** [July 1, 2023, 12:53pm UTC](https://discourse.julialang.org/t/principal-angles-between-subspaces/101044/3 "2023-07-01T12:53:42Z")

</div>

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