multiple regression through coordinate descent

mrc.R

mrc <- function(X, y, max.iter = 1000, tol = .Machine$double.eps) {
  # Multiple regression using coordinate descent
  # Input: X matrix, y vector
  # Output: estimated beta vector
  
  # init ----
  bhat <- numeric(ncol(X))
  convergence <- FALSE
  i <- 0L
  
  # Coordinate descent ----
  # calculate the sum of squares of each x
  xss <- apply(X, 2, crossprod)
  
  while (!convergence && i < max.iter) {
    i <- i + 1
    bprev <- bhat # remember previous beta estimate
    
    for (j in 1:length(bhat)) {
      # calculate residuals w.r.t. remaining variables
      yres <- y - X[,-j] %*% bhat[-j]
      # calculate current bhat using univariate regression of y_res on x_j
      bhat[j] <- t(X[,j]) %*% yres / xss[j]
    }
    
    # calculate squared difference to determine convergence
    dif <- (bprev - bhat) %*% (bprev - bhat)
    if (dif < tol) convergence <- TRUE
  }
  
  if (!convergence) warning("Algorithm did not converge!")
  
  bhat
}