R/daarem_lasso_gaussian_ngtp2.R

In nchenderson/nidaarem: Nesterov Initialized Damped Anderson Acceleration with Epsilon Monotonicity

daarem_lasso_gaussian_ngtp2 <- function(par, X, y, lambda, stplngth, nlag, a1, kappa, maxiter, 
                                        tol, mtol, cycl.mon.tol) {
  ## DAAREM for Lasso problem where there is no monitoring of the objective function.
  num.params <- ncol(X)
  lasso.pen <- lambda
  Fdiff <- Xdiff <- matrix(0.0, nrow=num.params, ncol=nlag)
  resid_vals <- objfn_track <- rep(NA, maxiter + 2)
  kr <- .01
  if(length(mtol)==2) {
    rho <- .95
  } else {
    rho <- mtol[1]
  }
  
  yty <- sum(y*y)
  Xty <- crossprod(X, y)
  
  xold <- par
  xnew <- SoftThresh(xold + stplngth*(Xty - crossprod(X, X%*%xold)), lambda=lambda*stplngth)
  fold <- xnew - xold
  objfn_track[1] <-  LassoObjFn_ngtp(xold, X, Xty, yty, lasso.pen)
  resid_vals[1] <- sqrt(crossprod(fold))
  fnew <- SoftThresh(xnew + stplngth*(Xty - crossprod(X,X%*%xnew)), lambda=lambda*stplngth) - xnew
  resid_vals[2] <- sqrt(crossprod(fnew))
  ss.resids <- resid_vals[2]
  objfn_track[2] <-  LassoObjFn_ngtp(xnew, X, Xty, yty, lasso.pen)
  krg <- ss.resids    
  
  k <- 1
  count <- 0
  shrink.count <- 0
  shrink.target <- 1/(1 + a1^kappa)
  
  lambda.ridge <- 100000
  r.penalty <- 0
  conv <- TRUE
  num.em <- 0  ## number of EM fallbacks
  n.aa <- 0
  DD <- 10
  n.viol <- 0
  while(k < maxiter) {
    count <- count + 1
    
    Fdiff[,count] <- fnew - fold
    Xdiff[,count] <- xnew - xold
    
    np <- count
    Ftmp <- matrix(Fdiff[,1:np], nrow=length(fnew), ncol=np)
    Xtmp <- matrix(Xdiff[,1:np], nrow=length(fnew), ncol=np)
    
    tmp <- La.svd(Ftmp)
    dvec <- tmp$d
    uy <- crossprod(tmp$u, fnew)
    uy.sq <- uy*uy
    
    max.d <- max(tmp$d)
    min.d <- min(tmp$d)
    cond.number <- ifelse(max.d==min.d, 1, max.d/min.d)  ## to take care of cases where max.d=min.d=0
    if(cond.number > 1e14) {
      shrink.count <- shrink.count - 2
    }
    
    ### Still need to compute Ftf
    Ftf <- sqrt(sum(as.vector(crossprod(Ftmp, fnew))^2))
    tmp_lam <- DampingFind(uy.sq, dvec, a1, kappa, shrink.count, Ftf, lambda.start=lambda.ridge, r.start=r.penalty)
    lambda.ridge <- tmp_lam$lambda
    r.penalty <- tmp_lam$rr
    
    dd <- (dvec*uy)/(dvec^2 + lambda.ridge)
    gamma_vec <- crossprod(tmp$vt, dd)
    
    xbar <- xnew - drop(Xtmp%*%gamma_vec)
    fbar <- fnew - drop(Ftmp%*%gamma_vec)
    
    x.propose <- xbar + fbar
    
    ftmp <- SoftThresh(x.propose + stplngth*(Xty - crossprod(X,X%*%x.propose)), lambda=lambda*stplngth) - x.propose
    ss.tmp <- sqrt(crossprod(ftmp))
    
    theor.val <- (DD*krg)/((n.aa + 1)^(1 + .05))
    if(ss.tmp <= min(ss.resids + rho^k, theor.val)) {
      ## Increase delta
      fold <- fnew
      xold <- xnew
      
      xnew <- x.propose
      fnew <- ftmp
      ss.resids <- ss.tmp
      
      shrink.count <- shrink.count + 1
      n.aa <- n.aa + 1
    } else {
      ## Keep delta the same
      fold <- fnew
      xold <- xnew
      
      xnew <- SoftThresh(xold + stplngth*(Xty - crossprod(X,X%*%xold)), lambda=lambda*stplngth)
      fnew <- SoftThresh(xnew + stplngth*(Xty - crossprod(X,X%*%xnew)), lambda=lambda*stplngth) - xnew
      ss.resids <- sqrt(crossprod(fnew))
    }
    if(ss.resids < tol & count==nlag) break
    
    resid_vals[k + 2] <- ss.resids
    objfn_track[k + 2] <-  LassoObjFn_ngtp(xnew, X, Xty, yty, lasso.pen)
    if(count==nlag) {
      count <- 0
      if(length(mtol)==2) {
        rho = ifelse(rho==mtol[1], mtol[2], mtol[1])
      } 
      ## restart count
    }
    shrink.target <-  1/(1 + a1^(kappa - shrink.count))
    k <- k+1
  }
  value.obj <- LassoObjFn_ngtp(xnew, X, Xty, yty, lasso.pen)
  resid_vals <- resid_vals[!is.na(resid_vals)]
  objfn_track <- objfn_track[!is.na(objfn_track)]
  
  if(k >= maxiter) {
    conv <- FALSE
    warning("Algorithm did not converge")
  }
  return(list(par=c(xnew), fpevals = k, value.objfn=value.obj, objfn.track=objfn_track,
              convergence=conv, resid.track=resid_vals, stplngth=stplngth, num.aa=n.aa, n.viol=n.viol))
}