R/nostra.R

In AliceGiampino/gradient: Tidy package

#set.seed(8675309)
n = 1000
x1 = rnorm(n)
x2 = rnorm(n)

#input
y = 1 + .5*x1 + .2*x2 + rnorm(n)
x <- cbind(1, x1, x2)
stepsize <- 1e-5
tolerance <- 1e-10
maxit <- 10000
b <- c(0.1, 0.1, 0.1)
verbose <- TRUE


#2
i <- 0

while(i < maxit){

  if(abs(max(b)) > 1e300) stop("Beta trooooppo grande: diminuisci la step size!")

  grad <- 2*t(x)%*%(x%*%b-y)
  bp <- b
  b <- bp - grad*stepsize

  if(max(abs(b - bp)) < tolerance){
    if(verbose) message("Tolerance reached")
    break()
  }

  i <- i + 1
}





# Steepest Descent --------------------------------------------------------


#input
y = 1 + .5*x1 + .2*x2 + rnorm(n)
x <- cbind(1, x1, x2)
stepsize <- 1e-5
tolerance <- 1e-10
maxit <- 10000
b <- c(0.1, 0.1, 0.1)
verbose <- TRUE


#2
i <- 0

while(i < maxit){

  if(abs(max(b)) > 1e300) stop("Beta trooooppo grande: diminuisci la step size!")

  hessian <- 4*t(x)%*%x
  grad <- 2*t(x)%*%(x%*%b-y)
  bp <- b
  step <- sqrt(sum(grad^2))/(t(grad)%*%hessian%*%grad)
  b <- bp - grad%*%step

  if(max(abs(b - bp)) < tolerance){
    if(verbose) message("Tolerance reached")
    break()
  }

  i <- i + 1
}
b



# Functions ---------------------------------------------------------------

# this could be your function signature
linear_gd_optim <- function(par,             # beta(0)
                            x,               # data predictors
                            y,               # response variable
                            tolerance=1e-3,  # tolerance
                            maxit=1000,      # max iteration, not to run forever
                            stepsize=1e-3,   # stepsize parameter
                            verbose=T        # should the function write messages during
) {

  i <- 0

  while(i < maxit){

    if(abs(max(b)) > 1e300) stop("Beta too large: decrease the step size!")

    grad <- 2*t(x)%*%(x%*%b-y)
    bp <- b
    b <- bp - grad*stepsize

    if(max(abs(b - bp)) < tolerance){
      if(verbose) message("Tolerance reached")
      break()
    }

    i <- i + 1
  }

  return(b)
}


linear_sd_optim <- function(par,             # beta(0)
                            x,               # data predictors
                            y,               # response variable
                            tolerance=1e-3,  # tolerance
                            maxit=1000,      # max iteration, not to run forever
                            stepsize=1e-3,   # stepsize parameter
                            verbose=T        # should the function write messages during
) {

  i <- 0

  while(i < maxit){

    if(abs(max(b)) > 1e300) stop("Beta too large: decrease the step size!")

    hessian <- 4*t(x)%*%x
    grad <- 2*t(x)%*%(x%*%b-y)
    bp <- b
    step <- sqrt(sum(grad^2))/(t(grad)%*%hessian%*%grad)
    b <- bp - grad%*%step

    if(max(abs(b - bp)) < tolerance){
      if(verbose) message("Tolerance reached")
      break()
    }

    i <- i + 1
  }

  return(b)
}