1. 牛顿算法原理
与梯度下降算法一样，牛顿拉斐逊算法也是优化算法的一种，可用于对多维参数的优化问题，同时也可用于求解方程的根。牛顿算法在优化算法中有着重要的地位，很多拟牛顿算法都是从该算法沿伸出来的，例如，BFGS算法、L-BFGS, DFP, Broyden算法。在梯度下降法中，该方法主要利用目标函数局部可微可导的性质，而牛顿法则是利用局部的一阶和二阶偏导信息，求得出近似函数的全局最小值。相比最速下降法，牛顿法带有一定对全局的预测性，收敛速度也更快。
对于多维参数非线性优化问题，假设目标函数为均方误差函数为\(f=J(\theta_0,\theta_1,\ldots,\theta_n)=\sum_{i=1}^{m}(\theta^T \cdot x^{(i)} - y^{(i)})^2\)，优化任务是求函数的极值，则可以转化为求解目标函数\(f\)的导数\(f'=0\)的问题，即求解方程或方程组\(f'= \frac{\partial f}{\partial \theta} = 0\)的可行解。牛顿法的具体推导过程如下：
对于求解方程\(f'=0\)的根，我们可以把\(f(x)\)按泰勒级数展开，展开成二阶形式：
\[f(\theta+\vartriangle \theta)=f(\theta)+f'(\theta) \vartriangle \theta + \frac{1}{2}f''(\theta) \vartriangle \theta^2\]当且仅当\(\vartriangle x \rightarrow 0\)时，上式成立。然后，目标求解方程\(f'=0\)的根等价于\(f'=\frac{\partial f}{\partial \theta}=\frac{f(\theta + \vartriangle \theta) - f(\theta)}{\vartriangle \theta}=f'(\theta) + \frac{1}{2} f''(\theta) \vartriangle \theta = 0\)。
因此，对于均方误差函数的优化问题，我们可以获得\(\vartriangle \theta = - \frac{2 f'(\theta)}{f''(\theta)}\)，进而每次迭代更新的\(\theta\)值为：
\[\theta_{n+1} := \theta_n - \alpha /(2m) * \frac{2 f'(\theta)}{f''(\theta)}=\theta_n - \alpha*\frac{\sum^m_{i=1}\,(\theta^T \cdot x^{(i)} - y^{(i)})x^{(i)}}{\sum^m_{i=1}\, x^{(i)^2}}\]上面求得的二阶导数当且仅当参数为一维时才成立，当参数为多维时，牛顿算法的迭代公式可以写成：
\(\theta_{n+1} := \theta_n - [Hf(\theta)]^{-1} \nabla f(\theta)\)，其中，\(H\)是hessian矩阵，定义为：
\[H(f)= \left( \begin{matrix} \frac{\partial^2 f}{\partial \theta^2_1} & \frac{\partial^2 f}{\partial \theta_1 \partial \theta_2} \cdots & \frac{\partial^2 f}{\partial \theta_1 \partial \theta_n} \\\ \frac{\partial^2 f}{\partial \theta_2 \partial \theta_1} & \frac{\partial^2 f}{\partial \theta^2_2} \cdots & \frac{\partial^2 f}{\partial \theta_2 \partial \theta_n} \\\ \vdots & \ddots & \vdots \\\ \frac{\partial^2 f}{\partial \theta_n \partial \theta_1} & \frac{\partial^2 f}{\partial \theta_n \partial \theta_2} \cdots & \frac{\partial^2 f}{\partial \theta^2_n} \\\ \end{matrix}\right)\]而梯度\(\nabla f(\theta)\)则为Jacobi矩阵的特例，梯度的Jacobi矩阵就是hessian矩阵。
综上所述，与梯度下降算法相比，牛顿算法需要多计算目标函数的二阶导数和二阶偏导数，即hessian矩阵，还要计算hessian矩阵的逆，所以，虽然牛顿法可以用于计算高维参数的近似值，但是求解hessian矩阵及其逆运算带来的复杂性使得牛顿算法求解难度和复杂度大大增加，并且当hessian矩阵不可逆时，牛顿算法则会无解。同时，如果初始值离局部极小值太远，泰勒展开公式并不能对原函数进行良好的逼近。
鉴于以上问题，我们一般采用拟牛顿法（Quasi-Newton Method）进行目标函数的优化，此时，不再需要计算hessian矩阵，而是每一步迭代使用梯度向量更新hessian矩阵，获得hessian矩阵的近似值。具体拟牛顿算法的细节将会在以后的博文中分享。
2. 牛顿算法实例
与梯度下降算法一样，这里仍然采用iris数据中花萼的长度（Sepal.Length）来预测花瓣的长度（Petal.Length），其中\(x\)为自变量，\(y\)为因变量。实验分别采用梯度下降算法、和牛顿优化算法对比估计的结果：
#读取数据x和y
x <- matrix(iris[which(iris$Species != "setosa"), ]$Sepal.Length, ncol=1)
y <- matrix(iris[which(iris$Species != "setosa"), ]$Petal.Length, ncol=1)
#绘制二维数据的散点图
plot(x, y,
main="versicolor和virginica花的散点图",
sub="花萼的长度和花瓣的长度",
col="blue",
family='STKaiti' #mac下需要增加该命令，否则不显示中文字符
)
#梯度下降和牛顿拉斐逊算法
##1. 计算代价函数
computeCost <- function(independent, dependent, variable){
m = length(dependent)
J = 0
J = 1/(2*m)*sum((independent%*%variable-y)^2)
return (J)
}
##2. 牛顿拉斐逊算法
newton <- function(independent, dependent, variable, num_iters){
m = length(dependent)
n = length(variable[,1])
J_history = rep(0, num_iters)
T_history = matrix(0, nrow = num_iters, ncol = n)
hessian = matrix(0, nrow=n, ncol=n)
## 计算hessian矩阵
hessian[1,1] = sum(independent[,1]^2)
hessian[1,2] = sum(independent[,1]*independent[,2])
hessian[2,1] = sum(independent[,2]*independent[,1])
hessian[2,2] = sum(independent[,2]^2)
## 计算jacobi矩阵，即梯度向量
jacobi = t(independent) %*% (independent %*% variable - dependent)
for (i in 1:num_iters){
J_history[i] = computeCost(independent, dependent, variable)
##更新theta
variable = variable - (solve(hessian) %*% jacobi) / m
T_history[i,] = variable
## 设定停止迭代条件
if (i > 1){
thres = J_history[i-1] - J_history[i] ##注意是上一次减去本次计算的代价函数值
}
else{
thres = J_history[i] - 0
}
if (thres < 1e-05 | thres <= 0) ##如果代价函数前后变化值小于1e-05或者发生增长
break;
}
result <- list(variable, J_history, T_history)
return (result)
}
##3. 梯度下降算法
gradient <- function(independent, dependent, variable, alpha, num_iters){
m = length(dependent)
n = length(variable[,1])
J_history = rep(0, num_iters)
T_history = matrix(0, nrow = num_iters, ncol = n)
for (i in 1:num_iters){
J_history[i] = computeCost(independent, dependent, variable)
variable = variable - alpha / m * (t(x) %*% (independent %*% variable - y))
T_history[i,] = variable
if (i > 1){
thres = J_history[i-1] - J_history[i]
}
else{
thres = J_history[i] - 0
}
if (thres < 1e-05 | thres <= 0)
break;
}
result <- list(variable, J_history, T_history)
return (result)
}
#训练LR模型
xLength <- length(x[,1])
intercept <- matrix(rep(1, xLength), ncol=1)
x <- cbind(intercept, x)
theta <- matrix(rep(0, 2), ncol=1)
alpha <- 0.001
num_iters <- 1000
newton.result <- newton(x, y, theta, num_iters)
abline(newton.result[[1]][1], newton.result[[1]][2])
gradient.result <- gradient(x, y, theta, alpha, num_iters)
abline(gradient.result[[1]][1], gradient.result[[1]][2], col="red")
legend("topleft", legend=c("newton", "gradient-descent"),
bty='n', col=c("black","red"),
lty = c(1, 1), cex=0.8, y.intersp=0.8)
#学习曲线绘制
newton.cost <- newton.result[[2]][which(newton.result[[2]] > 0)]
plot(1:length(newton.cost), newton.cost,
main="牛顿拉斐逊算法学习曲线",
sub="迭代次数为1:100",
family='STKaiti')
gradient.cost <- gradient.result[[2]][which(gradient.result[[2]] > 0)]
plot(1:length(gradient.cost), gradient.cost,
main="梯度下降算法学习曲线",
sub="迭代次数为1:100",
family='STKaiti')
#输出结果
print(sprintf("牛顿拉斐逊算法求的theta值和代价函数值为：(%g,%g)和%g",
newton.result[[1]][1],
newton.result[[1]][2],
min(newton.cost)))
print(sprintf("梯度下降算法求的theta值和代价函数值为：(%g,%g)和%g",
gradient.result[[1]][1],
gradient.result[[1]][2],
min(gradient.cost)))
从输出结果可以看出，牛顿算法迭代的最终代价函数值为0.1058，要小于梯度下降算法的最终代价函数值0.12，且牛顿计算得到的\(\theta\)值与前一篇博文(梯度下降算法)[http://www.hanlongfei.com/%E5%87%B8%E4%BC%98%E5%8C%96/2015/07/29/gradient/]中采用单纯形和模拟退火算法迭代得到的估计值更为接近，优化效果要好于梯度下降算法。
从迭代次数和学习曲线可以看出，梯度下降算法在迭代效率上要低于牛顿拉斐逊算法，牛顿算法在100次左右即可逼近最优值，而梯度下降算法在学习率为0.001情况下需要140-150次左右才能达到次优值（代价函数要大一些）。这里需要指出的是，适当调整学习率也可能会更快的逼近最优值，迭代次数小于牛顿法，但是，其效率还是要小于拟牛顿算法（L-BFGS-B）。