fastRUVIII.R

library(ruv)
library(rsvd) 

fastRUVIII = function(Y, M, ctl, k=NULL, eta=NULL, average=FALSE, fullalpha=NULL){
  # Assumes good input
  Y = RUV1(Y,eta,ctl)
  m = nrow(Y)
  Y0 = fast_residop(Y, M)
  fullalpha = diag(rsvd(Y0)$d) %*% t(rsvd(Y0)$v)
  alpha = fullalpha[1:k,,drop=FALSE]
  ac = alpha[,ctl,drop=FALSE]
  W = Y[,ctl]%*%t(ac)%*%solve(ac%*%t(ac))
  newY = Y - W%*%alpha
  return(list(newY = newY, fullalpha=fullalpha))
}

fast_residop <- function (A, B) {
  return(A - B %*% solve(t(B) %*% B) %*% (t(B) %*% A))
}