牛顿法解机器学习中的Logistic回归

这仍然是近期系列文章中的一篇。在这一个系列中，我打算把机器学习中的Logistic回归从原理到应用详细串起来。最初我们介绍了在Python中利用Scikit-Learn来建立Logistic回归分类器的方法

Python机器学习之Logistic回归

此后，我们对上述文章进行了更深一层的探讨，介绍了利用Logistic回归在自然语言处理中的应用（对微博进行Sentiment Analysis）

自然语言处理实战之微博情感偏向分析

从应用角度介绍了Logistic回归之后，我们又从源头介绍了Logistic回归的数学原理

机器学习之详解Logistic回归

在这篇文章的最后，我们得到了一个似然函数为

然后我们的目标是求出使这一似然函数值最大的参数估计，于是对函数取对数，再求一阶偏导数得到 

从这个偏导数出发，在实际中有很多方法可以解决前面的似然函数最大化参数估计问题，而我们这里要介绍的就是其中非常重要的一种，即牛顿法和拟牛顿法。

牛顿迭代法解方程的简单回顾

现代计算中涉及大量的工程计算问题，这些计算问题往往很少采用我们通常在求解计算题甚至是考试时所采用的方法，因为计算机最擅长的无非就是“重复执行大量的简单任务”，所以数值计算方面的迭代法在计算机时代便有了很大的作用。例如我们之前介绍过的利用“Jacobi迭代法与Gauss-Seidel迭代法”求解方程组的方法。另外一个例子则是利用牛顿迭代法近似求解方程的方法。牛顿迭代又称为牛顿-拉夫逊（拉弗森）方法（Newton-Raphson method）。如果读者对这部分内容感兴趣，可以详细参阅

牛顿迭代法与一道经典编程问题

我们在此做简要回顾。有时候某些方程的求根公式可能很复杂（甚至有些方程可能没有求根公式），导致求解困难。这时便可利用牛顿法进行迭代求解。

假设我们要求解方程f(x)=0的根，首先随便找一个初始值x0，如果x0不是解，做一个经过(x0,f(x0)) 这个点的切线，与x轴的交点为x1。同样的道理，如果x1不是解，做一个经过(x1,f(x1))这个点的切线，与x轴的交点为x2。 以此类推。以这样的方式得到的xi会无限趋近于 f(x)=0 的解。

判断xi是否是f(x)=0的解有两种方法：一是直接计算f(xi)的值判断是否为0，二是判断前后两个解xi和xi−1是否无限接近。经过(xi,f(xi))这个点的切线方程为（注意这也是一元函数的一阶泰勒展式） 

其中，f′(x)为f(x)的导数。令切线方程等于 0，即可求出

于是乎我们就得到了一个迭代公式，而且它必然在 f(x∗)=0处收敛，其中x∗就是方程的根，由此便可对方程进行迭代求根。

牛顿法在最优化问题中的应用

假设当前任务是优化一个目标函数 f，也就是求该函数的极大值或极小值问题，可以转化为求解函数 f 的导数 f′=0 的问题，这样求可以把优化问题看成方程 f′=0 求解问题。剩下的问题就和前面提到的牛顿迭代法求解很相似了。

这次为了求解方程 f′=0 的根，把原函数 f(x) 的做泰勒展开，展开到二阶形式（注意之前是一阶）：

当且仅当 Δx 无线趋近于0时，（可以舍得后面的无穷小项）使得等式成立。此时上式等价与： 

注意因为Δx 无线趋近于0，前面的的常数 1/2 将不再起作用，可以将其一并忽略，即 

求解 

得出迭代公式 

在之前的文章中我们也提到最优化问题除了用牛顿法来解之外，还可以用梯度下降法来解。但是通常来说，牛顿法可以利用到曲线本身的信息，比梯度下降法更容易收敛，即迭代更少次数。

再次联系到我们之前给出的一篇文章“Hessian矩阵与多元函数极值”，对于一个多维向量 X, 以及在点 X0 的邻域内有连续二阶偏导数的多元函数 f(X) ，可以写出该函数在点 X0 处的（二阶）泰勒展开式

其中，o(∥X−X0∥2) 是高阶无穷小表示的皮亚诺余项。而 Hf(X0) 是一个Hessian矩阵。依据之前的思路，忽略掉无穷小项，写出迭代公式即为 

由此，高维情况依然可以用牛顿迭代求解，但是问题是Hessian矩阵引入的复杂性，使得牛顿迭代求解的难度大大增加。所以人们又提出了所谓的拟牛顿法（Quasi-Newton method），不再直接计算Hessian矩阵，而是每一步的时候使用梯度向量更新Hessian矩阵的近似。这一点我们后续还会再讨论。

牛顿迭代法求解Logistic回归

现在回归到最开始的那个问题上。我们已经求出了Logistic回归的似然函数的一阶偏导数

由于lnL(w)是一个多元函数，变量是 w=w0,w1,⋯,wn ，所以根据多元函数求极值问题的规则，易知极值点处的导数一定均为零，所以一共需要列出n+1个方程，联立解出所有的参数。下面列出方程组如下 

当然，在具体解方程组之前需要用Hessian矩阵来判断极值的存在性。求Hessian矩阵就得先求二阶偏导，即

显然可以用Hessian矩阵来表示以上多元函数的二阶偏导数，于是有 

所以得到Hessian矩阵H=XTAX，可以看出矩阵A是负定的。线性代数的知识告诉我们，如果A是负定的，那么Hessian矩阵H也是负定的。也就是说多元函数存在局部极大值，这刚好符合我们的要求的最大似然估计相吻合。于是我们确信可以用牛顿迭代法来继续求解最优化问题。

对于多元函数求解零点，同样可以用牛顿迭代法，对于当前讨论的Logistic回归，可以得到如下迭代式 

其中H是Hessian矩阵，U的表达式如下 

由于Hessian矩阵H是对称负定的，将矩阵A提取一个负号出来，得到 

然后Hessian矩阵H变为H′=XTA′X，这样H′就是对称正定的了。那么现在牛顿迭代公式变为 

现在我们需要考虑如何快速地算得(H′)−1U，即解方程组(H′)−1X=U，通常的做法是直接用高斯消元法求解，但是这种方法的效率一遍比较低。而当前我们可以利用的一个有利条件是H′是对称正定的，所以可以用Cholesky矩阵分解法来解。Cholesky分解原理可以参考线性代数方面的书籍，此处不再赘述。

到这里，牛顿迭代法求解Logistic回归的原理已经介绍完了，但正如前面所提过的，在这个过程中因为要对Hessian矩阵求逆，计算量还是很大。而在后续的文章里我们会再来详细探讨拟牛顿法原理及应用，它是针对牛顿法的弱点进行了改进，更具实践应用价值。