Coding Gradient boosted machines in 100 lines of code
Motivation
有很多機器學習演算法。你不可能學習所有演算法的機制，但是，許多演算法都是從最成熟的演算法中發展出來的，例如：普通最小二乘法，梯度增強法，支援向量機，基於樹的演算法和神經網路。在STATWORX，我們每天討論演算法，以評估它們對特定專案或問題的有用性。無論如何，理解這些核心演算法是文獻中大多數機器學習演算法的關鍵。
雖然我喜歡閱讀機器學習研究論文，但數學有時難以理解。這就是為什麼我喜歡自己在R中實現演算法的原因。
在我隨後的兩篇部落格文章中，我將在不到150行的R程式碼中介紹兩種機器學習演算法。演算法將涵蓋所有核心機制，同時非常通用。你可以在我的GitHub上找到所有程式碼。
image.png
這篇博文將向您介紹梯度提升機器。 當然，有很多偉大的文章在理論上解釋了梯度增強，並附有一個動手例項。 這不是本博文的目的。 如果ni對演算法的理論和演變感興趣，我強烈建議您閱讀Bühlmann和Hothorn（2007）的論文。（1）本博文的目的是通過編寫簡單的R程式碼來建立演算法理論。 你不需要任何先前的演算法知識。
Gathering all puzzle pieces
# this is our loss function # y - the target # yhat - the fitted values loss <- function(y, yhat) return(mean(1/2*(y-yhat)^2))
我們的目標是找到這個損失函式的最小值，這是對均方誤差的適應。 該演算法利用了這樣的事實，即我們的損失函式中至少沒有斜率。 函式的梯度，即相對於yhat的負偏導數，描述了這個斜率。 因此，我們實際上可以重新設計我們的目標，即搜尋零或接近零的梯度，以最終實現我們減少損失的目標。
# the derivative of our loss function is the gradient
# y - the target
# yhat - the fitted values
gradient <- function(y, yhat) {return(y - yhat)}
現在是演算法的有趣部分。 在我們的例子中，梯度與殘差u = y - yhat一致。 請記住，我們希望漸變為零或接近零。 換句話說，我們希望狀態y = yhat。 該演算法通過簡單地迭代擬合（增強）殘差（梯度）而不是y本身來利用這一事實。 這意味著我們在每次迭代中更新漸變以改善我們的擬合。 這一步描述了我們的運動到最小的損失函式。 一旦我們看到理論在行動，我們將稍微闡述一下這個觀察。
The algorithm
# grad_boost
#' Fit a boosted linear model
#'
#' @param formula an object of class formula
#' @param data a data.frame or matrix
#' @param nu a numeric within the range of [0,1], the learning rate
#' @param stop a numeric value determining the total boosting iterations
#' @param loss.fun a convex loss function which is continuously differentiable
#' @param grad.fun a function which computes the gradient of the loss function
#' @param yhat.init a numeric value determining the starting point
#'
#' @return \itemize{
#' \item theta - this is our estimator
#' \item u - this is the last gradient value
#' \item fit - our fitted values, i.e. X %*% theta
#' \item formula - the underlying formula
#' \item data - the underlying data
#' }
#' @export
#'
#' @examples # Complete runthrough see: www.github.com/andrebleier/cheapml
grad_boost <- function(formula, data, nu = 0.01, stop,
grad.fun, loss.fun, yhat.init = 0) {
# coerce to data.frame
data <- as.data.frame(data)
# handle formula
formula <- terms.formula(formula)
# get the design matrix
X <- model.matrix(formula, data)
# extract target
y <- data[, as.character(formula)[2]]
# initial fit
fit <- yhat.init
# initialize the gradient with yhat.init
u <- grad.fun(y = y, yhat = fit)
# initial working parameter (our theta)
# this is just an empty body
theta <- rep(0, ncol(X))
# track loss
loss <- c()
# boost from 1 to stop
for (i in 1:stop) {
# estimate working parameter (via OLS)
# theta_i = (X'X)^-1 X'y
# This is the (base procedure) where you can literally place
# any optimization algorithm
base_prod <- lm.fit(x = X, y = u)
theta_i <- coef(base_prod)
# update our theta
theta <- theta + nu*as.vector(theta_i)
# new fitted values
fit <- fit + nu * fitted(base_prod)
# update gradient
u <- grad.fun(y = y, yhat = fit)
# update loss
loss <- append(loss, loss.fun(y = y, yhat = fit))
}
# naming the estimator
names(theta) <- colnames(X)
# define return
return(list(theta = theta, u = u, fit = fit, loss = loss,
formula = formula, data = data))
}
for迴圈使得編碼少於30行程式碼。 其他一切只是輸入/輸出處理等等。
# coerce to data.frame data <- as.data.frame(data) # handle formula formula <- terms.formula(formula) # get the design matrix X <- model.matrix(formula, data) # extract target y <- data[, as.character(formula)[2]] # initial fit fit <- f.init # initialize the gradient with yhat.init u <- grad.fun(y = y, yhat = fit) # initial working parameter (our theta) # this is just an empty body theta <- rep(0, ncol(X)) # track loss loss <- c()
他的第一部分只是輸入處理，但我們可以看到我們必須初始化漸變，從最佳猜測擬合yhat.init開始，例如， 0 初始化之後，我們需要一些容器來儲存我們的損失和我們的估算器theta。
下一部分是關於我們的優化，我們迭代地接近我們預先指定的損失函式的最小值。
# boost from 1 to stop
for (i in 1:stop) {
# estimate working parameter (via OLS)
# theta_i = (X'X)^-1 X'y
# This is the (base procedure) where you can literally place
# any optimization algorithm
base_prod <- lm.fit(x = X, y = u)
theta_i <- coef(base_prod)
# update our theta
theta <- theta + nu*as.vector(theta_i)
# new fitted values
fit <- fit + nu * fitted(base_prod)
# update gradient
u <- grad.fun(y = y, yhat = fit)
# update loss
loss <- append(loss, loss.fun(y = y, yhat = fit))
}
Unveiling the magic
最初我已經說過梯度增強是機器學習演算法的載體，而不是機器學習演算法本身。從本質上講，這是因為我們可以在這裡使用任何其他優化函式而不是lm.fit。大多數梯度增強機器使用基於樹的演算法，例如， xgboost。這使得梯度增強機器成為非常獨特的機器學習演算法。
References
Bühlmann, P., & Hothorn, T. (2007). Boosting algorithms: Regularization, prediction and model fitting. Statistical Science, 22(4), 477-505
Friedman, J. H. (2001). Greedy function approximation: a gradient boosting machine. Annals of statistics, 1189-1232.