R/angle.R

#########################################
#' @title Standardized last principal angle
#' @description Standardised last principal angle between the subspaces generated by the columns of A and B.
#' @param A numerical matrix of size \emph{p} by \emph{k}
#' @param B numerical matrix of size \emph{q} by \emph{l}
#' @return Standardised last principal angle between A and B.
#' @author Tom Reynkens
#' @references Bjorck, A. and Golub, G. H. (1973), ``Numerical Methods for Computing Angles Between Linear Subspaces,'' \emph{Mathematics of Computation}, 27, 579--594.
#' @references Hubert, M., Rousseeuw, P. J., and Vanden Branden, K. (2005), ``ROBPCA: A New Approach to Robust Principal Component Analysis,'' \emph{Technometrics}, 47, 64--79.
#' @references Hubert, M., Reynkens, T., Schmitt, E. and Verdonck, T. (2016), ``Sparse PCA for High-Dimensional Data With Outliers,'' \emph{Technometrics}, 58, 424--434.

Angle <- function (A, B) {
  A <- pracma::orth(A)
  B <- pracma::orth(B)
  # Check rank and swap
  if (ncol(A) < ncol(B)) {
    tmp <- A
    A <- B
    B <- tmp
  }
  # Compute the projection according to Equation 13 in Bjorck and Golub (1973).
  B <- B - A%*%(t(A)%*%B)
  # Make sure its magnitude is less than 1.
  theta <- asin(min(1, max(svd(B)$d)))
  return(theta/(pi/2))
}

Try the ltsspca package in your browser

Any scripts or data that you put into this service are public.

ltsspca documentation built on Oct. 9, 2019, 5:06 p.m.