|
梯度下降和牛顿法都是用来求解最优化问题的，在机器学习中应用甚广，尤其是梯度下降，它在各种分类和回归模型的求解中都会用到。那什么是梯度下降，什么是牛顿法呢，本文就带你一探究竟。 梯度下降假定最优化的目标是求解使得 取最小值时的待估参数 ，我们来看梯度下降是如何来求解待估参数的： 1 初始化待估参数初始化，可以认为设定一个初值，当然也可以让计算机随机生成。 2 循环迭代采用如下公式对待估参数就行循环迭代： 
其中， 为一阶偏导数即所谓的梯度，可以看到每次迭代，会向梯度相反的方向移动 ，这里  为步长，如果步长太小，迭代可能导致太慢，如果步长太大，可能可能跳过局部最小值，不能保证收敛，所以计算时需要选取合适的进行迭代运算。 3 计算完成当 不再减小或者达到预设的循环上界时，计算终止。 实际上，上述梯度下降算法为批量梯度下降，本文以仅此为例来讲解，因为当你理解之后你会发现，其他类型的梯度下降算法均为此算法的变种。
牛顿法待优化问题同上述梯度下降算法。 1 初始化待估参数同上。 2 循环迭代循环迭代的公式如下： 
其中， 为一阶偏导数，为二阶导，这个公式实际上是在 处的二阶泰勒展开式的变形，由于它考虑了目标函数的二阶导，因此它的收敛速度也更快，但是它对目标函数的要求也更严格，需要存在二阶导。而且从上式中还可以看到它与梯度下降的另外一个区别是不需要设置步长参数了。 3 计算完成当不再减小或者达到预设的循环上界时，计算终止。 R语言实现了解了两种算法的基本概念后，我们尝试封装R函数，来实现对一元线性回归模型系数的估计。1 样例数据采用蒙特卡洛模拟生成一组变量数据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")
可以看到这两组数据之间存在明显的线性相关。于是我们可以采用一元线性回归模型去拟合，形如：
求解上述模型就是要估计出上式中截距项和斜率的大小。接下来我们就分别采用梯度下降和牛顿法来求解上述模型。 2 编写函数构造普通线性回归的损失函数：
梯度下降和牛顿算法的目的就是要极小化上述残差平方和。上代码： # 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 对比算法效果查看梯度下降和牛顿算法的效果，并将其与系统自带函数的估计结果做对比： # 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
可以看到，梯度下降经过2964步后收敛，牛顿算法经过仅仅123步就收敛。这是因为牛顿算法考虑了二阶导，即梯度的变化，所以他在迭代次数上显得比梯度下降更有优势。 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)
总结梯度下降算法和牛顿法各有优劣： 梯度下降法考虑了目标函数的一阶偏导、以负梯度方向作为搜索方向去寻找最小值，需要认定合适的步长，迭代次数多，容易陷入局部最小值。 牛顿法同时考虑了目标函数的一、二阶偏导数，相当于考虑了梯度的梯度，所以能确定合适的搜索方向加快收敛，但牛顿法要求也更严格，目标函数必须存在连续的一、二阶偏导数，且海森矩阵必须正定，迭代次数虽少，但是当待估参数较多时，单次的运算量会远大于梯度下降算法。 所以在具体的数据分析和挖掘中，需要权衡利弊后再选择更为合适的算法。
|
