trainLinearReg <- function(X, y, lambda) {
  
  #   Trains regularized linear regression given a dataset (X, y) 
  #   and a regularization parameter lambda
  #
  #   theta <- trainLinearReg(X, y, lambda) trains linear regression using
  #   the dataset (X, y) and regularization parameter lambda. Returns the
  #   trained parameters theta.
  
  if (class(X) == "numeric") 
    X <- t(X)
  
  theta <- rep(0, ncol(X))
  
  opt <-  optim(par = theta, 
                fn = LinearRCostFunctionReg,
                gr = LinearRGradientReg,
                method = "BFGS",
                X, y, lambda,
                control = list(maxit = 400))

  theta <- opt$par
  
  return(theta)
}