|
梯度下降和牛頓法都是用來求解最優(yōu)化問題的，在機器學習中應用甚廣，尤其是梯度下降，它在各種分類和回歸模型的求解中都會用到。那什么是梯度下降，什么是牛頓法呢，本文就帶你一探究竟。 梯度下降假定最優(yōu)化的目標是求解使得 取最小值時的待估參數(shù) ，我們來看梯度下降是如何來求解待估參數(shù)的： 1 初始化待估參數(shù)初始化，可以認為設定一個初值，當然也可以讓計算機隨機生成。 2 循環(huán)迭代采用如下公式對待估參數(shù)就行循環(huán)迭代： 
其中， 為一階偏導數(shù)即所謂的梯度，可以看到每次迭代，會向梯度相反的方向移動 ，這里  為步長，如果步長太小，迭代可能導致太慢，如果步長太大，可能可能跳過局部最小值，不能保證收斂，所以計算時需要選取合適的進行迭代運算。 3 計算完成當 不再減小或者達到預設的循環(huán)上界時，計算終止。 實際上，上述梯度下降算法為批量梯度下降，本文以僅此為例來講解，因為當你理解之后你會發(fā)現(xiàn)，其他類型的梯度下降算法均為此算法的變種。
牛頓法待優(yōu)化問題同上述梯度下降算法。 1 初始化待估參數(shù)同上。 2 循環(huán)迭代循環(huán)迭代的公式如下： 
其中， 為一階偏導數(shù)，為二階導，這個公式實際上是在 處的二階泰勒展開式的變形，由于它考慮了目標函數(shù)的二階導，因此它的收斂速度也更快，但是它對目標函數(shù)的要求也更嚴格，需要存在二階導。而且從上式中還可以看到它與梯度下降的另外一個區(qū)別是不需要設置步長參數(shù)了。 3 計算完成當不再減小或者達到預設的循環(huán)上界時，計算終止。 R語言實現(xiàn)了解了兩種算法的基本概念后，我們嘗試封裝R函數(shù)，來實現(xiàn)對一元線性回歸模型系數(shù)的估計。1 樣例數(shù)據(jù)采用蒙特卡洛模擬生成一組變量數(shù)據(jù)x和y： # generate random data
# y ~ a + b*x
set.seed(2017)
n <- 100
x <- runif(n,1,10)
y <- x + rnorm(n)
plot(x,y,pch=20,main = "Monte Carlo Simulation for Simple Linear Regression")
可以看到這兩組數(shù)據(jù)之間存在明顯的線性相關。于是我們可以采用一元線性回歸模型去擬合，形如：
求解上述模型就是要估計出上式中截距項和斜率的大小。接下來我們就分別采用梯度下降和牛頓法來求解上述模型。 2 編寫函數(shù)構(gòu)造普通線性回歸的損失函數(shù)：
梯度下降和牛頓算法的目的就是要極小化上述殘差平方和。上代碼： # gradient descent algorithm
gd <- function(x,y,a0,a1,alpha=0.01,tol=1e-5,M=5000){
# Args:
# x -- independent variable
# y -- dependent variable
# a0 -- initial value for intercept
# a1 -- initial value for slope
# alpha -- step size
# tol -- tolerance
# M -- Maximum Iterations
# Returns:
# Iterated value sequence, is a dataframe.
i <- 1
res <- data.frame(a0=a0,a1=a1)
repeat{
J0 <- 1/2*sum((a0+a1*x-y)^2)
a0_hat <- a0 - alpha*mean(a0+a1*x-y)
a1_hat <- a1 - alpha*mean((a0+a1*x-y)*x)
a0 <- a0_hat
a1 <- a1_hat
J1 <- 1/2*sum((a0_hat+a1_hat*x-y)^2)
res <- rbind(res,data.frame(a0=a0,a1=a1))
if( abs(J1-J0) < tol | i >= M )
break
i <- i + 1
}
return(res)
}
# Newton's method
nt <- function(x,y,a0,a1,tol=1e-5,M=500){
# Args:
# x -- independent variable
# y -- dependent variable
# a0 -- initial value for intercept
# a1 -- initial value for slope
# tol -- tolerance
# M -- Maximum Iterations
# Returns:
# Iterated value sequence, is a dataframe.
i <- 1
res <- data.frame(a0=a0,a1=a1)
repeat{
J0 <- 1/2*sum((a0+a1*x-y)^2)
a0_hat <- a0 - mean(a0+a1*x-y)
a1_hat <- a1 - mean((a0+a1*x-y)*x)/mean(x^2)
a0 <- a0_hat
a1 <- a1_hat
J1 <- 1/2*sum((a0_hat+a1_hat*x-y)^2)
res <- rbind(res,data.frame(a0=a0,a1=a1))
if( abs(J1-J0) < tol | i >= M )
break
i <- i + 1
}
return(res)
}
3 對比算法效果查看梯度下降和牛頓算法的效果，并將其與系統(tǒng)自帶函數(shù)的估計結(jié)果做對比： # compare two algorithms
res.gd <- gd(x = x,y = y,a0 = 0,a1 = 0,tol=1e-8)
tail(res.gd,1)
## a0 a1
## 2964 0.2515242 0.9429152
res.nt <- nt(x = x,y = y,a0 = 0,a1 = 0,tol=1e-8)
tail(res.nt,1)
## a0 a1
## 123 0.2520766 0.9428287
# use lm function
fit <- lm(y~x)
(a <- coef(fit))
## (Intercept) x
## 0.2520778 0.9428331
可以看到，梯度下降經(jīng)過2964步后收斂，牛頓算法經(jīng)過僅僅123步就收斂。這是因為牛頓算法考慮了二階導，即梯度的變化，所以他在迭代次數(shù)上顯得比梯度下降更有優(yōu)勢。 3 迭代路徑可視化# plot for gradient descent algorithm
nr <- nrow(res.gd)
ind <- c(1:10,1:floor(nr/100)*100,nr)
a0 <- res.gd$a0[ind]
a1 <- res.gd$a1[ind]
plot(a0,a1,cex=0.5,xlim = c(0,0.3),ylim = c(0,1),pch=20,main = "Gradient Descent Algorithm")
arrows(head(a0,-1),head(a1,-1),tail(a0,-1),tail(a1,-1),length=0.05)
points(a[1],a[2],col=2,cex=1.5,pch=20)
# plot for Newton's method
a0 <- res.nt$a0
a1 <- res.nt$a1
plot(a0,a1,xlim = c(0,6),ylim = c(0,1),pch=20,main = "Newton's method")
arrows(head(a0,-1),head(a1,-1),tail(a0,-1),tail(a1,-1),length=0.08)
points(a[1],a[2],col=2,cex=1.5,pch=20)
總結(jié)梯度下降算法和牛頓法各有優(yōu)劣： 梯度下降法考慮了目標函數(shù)的一階偏導、以負梯度方向作為搜索方向去尋找最小值，需要認定合適的步長，迭代次數(shù)多，容易陷入局部最小值。 牛頓法同時考慮了目標函數(shù)的一、二階偏導數(shù)，相當于考慮了梯度的梯度，所以能確定合適的搜索方向加快收斂，但牛頓法要求也更嚴格，目標函數(shù)必須存在連續(xù)的一、二階偏導數(shù)，且海森矩陣必須正定，迭代次數(shù)雖少，但是當待估參數(shù)較多時，單次的運算量會遠大于梯度下降算法。 所以在具體的數(shù)據(jù)分析和挖掘中，需要權衡利弊后再選擇更為合適的算法。
|
