Skip to contents

Computes the Grassmann distance between two k-dimensional subspaces in R^p, represented by orthonormal basis matrices U and V (both p x k). The distance is `sqrt(sum(theta_i^2))`, where `theta_i` are the principal angles between the subspaces. Principal angles `theta_i = acos(s_i)`, where `s_i` are the singular values of `U^T V`.

Usage

grassmann_distance(U, V, tol = 1e-06, sv_floor = 1e-07)

Arguments

U

A numeric matrix (p x k) with orthonormal columns.

V

A numeric matrix (p x k) with orthonormal columns.

tol

Tolerance for checking orthonormality. Default 1e-6.

sv_floor

Smallest value for singular values before acos to avoid `acos(>1)` due to numerical precision. Values are clamped to `[-1+sv_floor, 1-sv_floor]`. Default 1e-7.

Value

The Grassmann distance (a non-negative scalar).

Examples

# Example from P. Edelman, T. A. Arias, and A. S. Edelman. "Geometry of algorithms
# with orthogonality constraints." SIAM Journal on Matrix Analysis and Applications, 1998.
p <- 3; k <- 2
U <- qr.Q(qr(matrix(rnorm(p*k), p, k)))
V <- qr.Q(qr(matrix(rnorm(p*k), p, k)))
# grassmann_distance(U, V)

# Identical subspaces
# grassmann_distance(U, U) # Should be 0

# Orthogonal subspaces (if k <= p/2)
if (k <= p/2) {
  # V_ortho <- qr.Q(qr(matrix(rnorm(p*k), p, k)), complete=TRUE)[, (k+1):(2*k), drop=FALSE]
  # if (ncol(V_ortho) == k) grassmann_distance(U, V_ortho) # Should be sqrt(k * (pi/2)^2)
}