Алгоритмы линейного регрессионного градиентного спуска в R дают разные результаты

Question

Я пытаюсь реализовать линейную регрессию в R с нуля без использования каких-либо пакетов или библиотек, используя следующие данные:

Хранилище данных машинного обучения UCI, Bike-Sharing-Dataset

Линейная регрессия была достаточно простой, вот код:

data <- read.csv("Bike-Sharing-Dataset/hour.csv")

# Select the useable features
data1 <- data[, c("season", "mnth", "hr", "holiday", "weekday", "workingday", "weathersit", "temp", "atemp", "hum", "windspeed", "cnt")]

# Split the data
trainingObs<-sample(nrow(data1),0.70*nrow(data1),replace=FALSE)

# Create the training dataset
trainingDS<-data1[trainingObs,]

# Create the test dataset
testDS<-data1[-trainingObs,]

x0 <- rep(1, nrow(trainingDS)) # column of 1's
x1 <- trainingDS[, c("season", "mnth", "hr", "holiday", "weekday", "workingday", "weathersit", "temp", "atemp", "hum", "windspeed")]

# create the x- matrix of explanatory variables
x <- as.matrix(cbind(x0,x1))

# create the y-matrix of dependent variables

y <- as.matrix(trainingDS$cnt)
m <- nrow(y)

solve(t(x)%*%x)%*%t(x)%*%y 

Следующим шагом является реализация пакетного обновления с градиентным спуском, и здесь я сталкиваюсь с проблемами. Я не знаю, откуда приходят ошибки или как их исправить, но код работает. Проблема в том, что производимые значения радикально отличаются от результатов регрессии, и я не знаю, почему.

Я реализовал две версии градиентного спуска пакетного обновления: результаты обоих алгоритмов отличаются друг от друга и от результатов регрессии:

# Gradient descent 1
gradientDesc <- function(x, y, learn_rate, conv_threshold, n, max_iter) {
  plot(x, y, col = "blue", pch = 20)
  m <- runif(1, 0, 1)
  c <- runif(1, 0, 1)
  yhat <- m * x + c
  MSE <- sum((y - yhat) ^ 2) / n
  converged = F
  iterations = 0
  while(converged == F) {
    ## Implement the gradient descent algorithm
    m_new <- m - learn_rate * ((1 / n) * (sum((yhat - y) * x)))
    c_new <- c - learn_rate * ((1 / n) * (sum(yhat - y)))
    m <- m_new
    c <- c_new
    yhat <- m * x + c
    MSE_new <- sum((y - yhat) ^ 2) / n
    if(MSE - MSE_new <= conv_threshold) {
      abline(c, m) 
      converged = T
      return(paste("Optimal intercept:", c, "Optimal slope:", m))
    }
    iterations = iterations + 1
    if(iterations > max_iter) { 
      abline(c, m) 
      converged = T
      return(paste("Optimal intercept:", c, "Optimal slope:", m))
    }
  }
  return(paste("MSE=", MSE))
}

А ТАКЖЕ:

grad <- function(x, y, theta) { # note that for readability, I redefined theta as a column vector
  gradient <-  1/m* t(x) %*% (x %*% theta - y) 
  return(gradient)
}
grad.descent <- function(x, maxit, alpha){
  theta <- matrix(rep(0, length=ncol(x)), ncol = 1)
  for (i in 1:maxit) {
    theta <- theta - alpha  * grad(x, y, theta)   
  }
  return(theta)
}

Если бы кто-то мог объяснить, почему эти две функции дают разные результаты, я был бы очень признателен. Я также хочу убедиться, что я на самом деле правильно выполняю градиентный спуск.

Наконец, как я могу отобразить результаты спуска с различными скоростями обучения и наложить эти данные на результаты самой регрессии?

РЕДАКТИРОВАТЬ Вот результаты запуска двух алгоритмов с альфа = .005 и 10000 итераций:

1)

> gradientDesc(trainingDS, y, 0.005, 0.001, 32, 10000)
TEXT_SHOW_BACKTRACE environmental variable.
[1] "Optimal intercept: 2183458.95872599 Optimal slope: 62417773.0184353"

2)

> print(grad.descent(x, 10000, .005))
                   [,1]
x0            8.3681113
season       19.8399837
mnth         -0.3515479
hr            8.0269388
holiday     -16.2429750
weekday       1.9615369
workingday    7.6063719
weathersit  -12.0611266
temp        157.5315413
atemp       138.8019732
hum        -162.7948299
windspeed    31.5442471

Answer 1

Чтобы дать вам пример того, как написать такие функции немного лучше, рассмотрите следующее:

gradientDesc <- function(x, y, learn_rate, conv_threshold, max_iter) {
  n <- nrow(x) 
  m <- runif(ncol(x), 0, 1) # m is a vector of dimension ncol(x), 1
  yhat <- x %*% m # since x already contains a constant, no need to add another one

  MSE <- sum((y - yhat) ^ 2) / n

  converged = F
  iterations = 0

  while(converged == F) {
    m <- m - learn_rate * ( 1/n * t(x) %*% (yhat - y))
    yhat <- x %*% m
    MSE_new <- sum((y - yhat) ^ 2) / n

    if( abs(MSE - MSE_new) <= conv_threshold) {
      converged = T
    }
    iterations = iterations + 1
    MSE <- MSE_new

    if(iterations >= max_iter) break
  }
  return(list(converged = converged, 
              num_iterations = iterations, 
              MSE = MSE_new, 
              coefs = m) )
}

Для сравнения:

ols <- solve(t(x)%*%x)%*%t(x)%*%y 

Сейчас,

out <- gradientDesc(x,y, 0.005, 1e-7, 200000)

data.frame(ols, out$coefs)
                    ols    out.coefs
x0           33.0663095   35.2995589
season       18.5603565   18.5779534
mnth         -0.1441603   -0.1458521
hr            7.4374031    7.4420685
holiday     -21.0608520  -21.3284449
weekday       1.5115838    1.4813259
workingday    5.9953383    5.9643950
weathersit   -0.2990723   -0.4073493
temp        100.0719903  147.1157262
atemp       226.9828394  174.0260534
hum        -225.7411524 -225.2686640
windspeed    12.3671942    9.5792498

Вот, x относится к вашему x как определено в вашем первом фрагменте кода. Обратите внимание на сходство коэффициентов. Тем не менее, также обратите внимание, что

out$converged
[1] FALSE

так что вы можете увеличить точность, увеличив количество итераций или поиграв с размером шага. Это также может помочь в первую очередь масштабировать ваши переменные.