R/l2e_regression_convex_MM.R

In jocelynchi/L2E-package-demo: Robust Structured Regression via the L2 Criterion

#' L2E convex regression - MM
#'
#' \code{l2e_regression_convex_MM} performs L2E convex regression via block coordinate descent
#' with MM for updating beta and modified Newton for updating tau.
#'
#' @param y response
#' @param beta initial vector of regression coefficients
#' @param tau initial precision estimate
#' @param max_iter maximum number of iterations
#' @param tol relative tolerance
#' @param Show.Time Report the computing time
#' @return Returns a list object containing the estimates for beta (vector) and tau (scalar),
#' the number of outer block descent iterations until convergence (scalar),
#' and the number of inner iterations per outer iteration for updating beta and eta (vectors)
#' @export
#' @examples
#' set.seed(12345)
#' n <- 200
#' tau <- 1
#' x <- seq(-2, 2, length.out=n)
#' f <- x^4 + x
#' y <- f + (1/tau)*rnorm(n)
#'
#' ## Clean
#' plot(x, y, pch=16, cex.lab=1.5, cex.axis=1.5, cex.sub=1.5, col='gray')
#' lines(x, f, lwd=3)
#'
#' tau <- 1
#' b <- y
#' sol <- l2e_regression_convex_MM(y, b, tau)
#'
#' plot(x, y, pch=16, cex.lab=1.5, cex.axis=1.5, cex.sub=1.5, col='gray')
#' lines(x, f, lwd=3)
#' cvx <- fitted(cobs::conreg(y, convex=TRUE))
#' lines(x, cvx, col='blue', lwd=3)
#' lines(x, sol$beta, col='red', lwd=3)
#'
#' ## Contaminated
#' ix <- 0:9
#' y[45 + ix] <- 14 + rnorm(10)
#'
#' plot(x, y, pch=16, cex.lab=1.5, cex.axis=1.5, cex.sub=1.5, col='gray')
#' lines(x, f, lwd=3)
#'
#' tau <- 1
#' b <- y
#' sol <- l2e_regression_convex_MM(y, b, tau)
#'
#' plot(x, y, pch=16, cex.lab=1.5, cex.axis=1.5, cex.sub=1.5, col='gray')
#' lines(x, f, lwd=3)
#' cvx <- fitted(cobs::conreg(y, convex=TRUE))
#' lines(x, cvx, col='blue', lwd=3)
#' lines(x, sol$beta, col='red', lwd=3)
#'
l2e_regression_convex_MM <- function(y,beta,tau,max_iter=1e2,tol=1e-4,Show.Time=TRUE) {

  if (tau <= 0) stop("Entered non-positive initial tau")

  time <- proc.time()

  # save the inner loop iters
  iter_beta <- rep(0, max_iter)
  iter_eta <- rep(0, max_iter)

  for (i in 1:max_iter) {
    # update beta
    beta_last <- beta
    res_beta <- update_beta_MM_convex(y,beta_last,tau,max_iter=1e2,tol=1e-4)
    beta <- res_beta$beta
    iter_beta[i] <- res_beta$iter

    # update tau/eta
    r <- y - beta
    eta_last <- log(tau)  # get eta as in line 9
    res_eta <- update_eta_bktk(r,eta_last, tol=tol) # update eta as in line 10-12
    eta <- res_eta$eta
    tau <- exp(eta) # update tau as in line 13
    iter_eta[i] <- res_eta$iter


    A <- norm(as.matrix(beta_last-beta),'f') < tol*(1 + norm(as.matrix(beta_last),'f'))
    B <- abs(eta_last-eta) < tol*(1 + abs(eta_last))  # tau_last could be negative, so need abs() here
    if (A & B) break
  }
  if(Show.Time) print(proc.time() - time)

  return(list(beta=beta,tau=tau, iter=i,
              iter_beta = iter_beta[1:i], iter_eta = iter_eta[1:i]))
}