R/constrained-orthonormal-least-squares.R

In bonStats/gcreg: General Constraint REGression

#' The COLS algorithm.
#' 
#' Provided for advanced users. Using cpm() is recommended as main interface.
#' 
#' @param par initial starting point of algorithm.
#' @param Y (scaled) observation matrix.
#' @param X orthonormalised design matrix.
#' @param oracle_fun function that determines if a given point is inside the constrained space.
#' @param control list of controls generated by cols_control().
#' @param aim_gamma The unconstrained solution to descend towards, on the orthonormal scale, computed as \code{crossprod(X, Y)} if not provided.
#' @param save_gamma Save and return a list of the gamma iterates.
#' @param ... additional parameters passed to control if not specified.
#' @return constrained minimum solution.
#' @keywords internal
#' @export

# base implementation
optim_cols <- function(par, Y, X, oracle_fun, control = cols_control(...), aim_gamma, save_gamma = F, ...) {
  # all in orthonormal basis
  # oracle function should have converter within
  # orcale function is set up by a helper function which tries different points? Or needs to be specified.
  
  p <- ncol(X)
  
  # aim: OLS beta orthonormal design matrix
  if(missing(aim_gamma)){
    aim_gamma <- crossprod(X, Y)
  }
  curr_gamma <- matrix(par, ncol = 1)
  
  # line search function in 1 dimension
  ls_find <- gen_line_search(oracle_fun)
  
  # convergence counter
  conv_counter <- 0
  
  # while loop index
  i <- 0
  
  # step size intiate
  step_size <- control$step_start
  
  # decide which method to use and set function
  if(control$method == "best-step"){
    
    find_gamma <- function(cur, aim, loop_i) { # loop_i currently not used
      
      pot_new_gammas <- lapply(1:p, FUN = ls_find, cur = cur, aim = aim)
      
      pot_moves <- diag(sapply(pot_new_gammas, function(g) {g - cur}))
      
      ix <- which.min(2 * as.vector(cur - aim) * pot_moves + pot_moves^2)
      
      list(gamma = pot_new_gammas[[ix]], coord = ix)
      
    }
    
  } else if(control$method == "avoid-boundary"){
    
    find_gamma <- function(cur, aim, loop_i) { # loop_i currently not used
      
      pot_new_gammas <- lapply(1:p, FUN = ls_find, cur = cur, aim = aim)
      
      pot_moves <- diag(sapply(pot_new_gammas, function(g) {g - cur}))
      
      ix <- which.max(abs(pot_moves/(aim-cur)))
      
      list(gamma = pot_new_gammas[[ix]], coord = ix)
      
    }
    
  } else if(control$method == "once-best-step"){
    
    find_gamma <- function(cur, aim, loop_i) { # loop_i currently not used
      
      if ((loop_i %% (p+1)) == 0){
        
        pot_new_gammas <- lapply(1:p, FUN = ls_find, cur = cur, aim = aim)
        
        pot_moves <- diag(sapply(pot_new_gammas, function(g) {g - cur}))
        
        ix <- which.min(2 * as.vector(cur - aim) * pot_moves + pot_moves^2)
        
        cat("best ix=",ix,"\n")
        
        return(list(gamma = pot_new_gammas[[ix]], coord = ix))
        
      } else {
        
        ix <- (loop_i %% p) + 1
        
        cat("gen ix=",ix,"\n")
        
        return(list(gamma = ls_find(index = ix, cur = cur, aim = aim), coord = ix))
        
      }
    }
    
  } else if(control$method == "up-walk") {
    
    find_gamma <- function(cur, aim, loop_i) {
      
      ix <- (loop_i %% p) + 1
      
      list(gamma = ls_find(index = ix, cur = cur, aim = aim), coord = ix)
      
    }
    
  } else {
    # method = "down-walk"
    find_gamma <- function(cur, aim, loop_i) {
      
      ix <- p - (loop_i %% p)
      
      list(gamma = ls_find(index = ix, cur = cur, aim = aim), coord = ix)
      
    }
  }
  
  if(save_gamma) gamms <- list()
  
  # search
  while(i <= control$maxit & conv_counter < 2 * p){
    
    # finds the full move possible based on the function defined by find_gamma above
    full_move <- find_gamma(cur = curr_gamma, aim = aim_gamma, loop_i = i)
    
    # reduces move to specified step size (aside: would work for change in one or several coordinates with changes)
    new_gamma <- replace(curr_gamma, full_move$coord, curr_gamma[full_move$coord] * (1 - step_size) + with(full_move, gamma[coord]) * step_size)
    
    if(save_gamma) gamms[[i+1]] <- new_gamma
    
    if(all(abs(new_gamma - curr_gamma) < control$tol)) {
      conv_counter <- conv_counter + 1
    } else {
      conv_counter <- 0
    }
    
    # update gamma
    curr_gamma <- new_gamma
    
    # increment
    i <- i + 1
    
    # change step size
    if((i %% (2 * p)) == 0){
      step_size <- min(step_size + control$step_increment, 1)
    }
    
  }
  
  attr(new_gamma,"convg") <- i != control$maxit
  
  pot_new_gammas <- lapply(1:p, FUN = ls_find, cur = new_gamma, aim = aim_gamma)
  
  attr(new_gamma,"pot_moves") <- diag(sapply(pot_new_gammas, function(g) {g - new_gamma}))
  
  if(save_gamma) attr(new_gamma,"gamma_iter") <- gamms
  
  return(new_gamma)

}