更多深度文章，请关注云计算频道：https://yq.aliyun.com/cloud
前言
在本篇博客中，我们要介绍的是牛顿法的原理，并且将之应用到实际的逻辑回归问题中。逻辑回归的主要知识点包括伯努利分布的对数似然和用来平滑的sigmoid函数。
我们还要介绍Hession，这是一个二阶偏导的方阵。看完了本片博客，您就知道如何使用Hession结合梯度来实现牛顿法。
和之前的博客一样，我们这篇也将从牛顿法的整体概述、数学推导以及编程实现几个方面展开。最终将理论和实践的结合，灵活运用牛顿法解决逻辑回归问题。
推荐学习的基础知识
- 导数和链式法则（微积分）
- 偏导数和梯度（多元微积分）
- 基本的向量操作（线性代数）
- 对Numpy有基本了解
- 独立概率的基本认识
数据
我们的本次实验采用的数据集是南波士顿的房产数据集。数据集中的每一行包含了一座房屋的价值，并且用一个布尔类型的数值表示是否有超过两个的浴室。
$$ \widehat{x} =HomeValue = <550000.00,600000.00,...578000.00>^{T} $$
$$ \widehat{y} = MoreThan2Bathrooms = <1,0,0,...1>^{T} $$
south_boston_bedroom_logreg_raw
数学模型
首先，我们需要掌握基本的线性回归模型。当我们给出房屋的价格的时候，该模型可以用二元分类器的方式为我们判断该房屋是否有两个以上的浴室。
接下来我们需要解决的是特征和权重的线性组合问题。我们需要将这个线性组合封装在一个函数中，而这个函数必须是平滑的，且值域是[0,1]。
按照惯例，我们选择满足以上要求的Logistic函数，它的公式如下所示：
$$ h(x)=1/(1+e^{-z}), $$
$$ z=\theta_{1} x+ \theta_{2} $$
Note:我们增加了\(theta_{2}\)参数来增加模型的灵活度，尽管我们的数据只有一维，但在这里我们使用了一个二维模型。
sigmoid_function
在线性回归中，我们需要定义一个平方和的目标函数，逻辑回归也是它相似。我们用我们的假设函数h(x)来生成似然函数，最终使用似然函数来求得合适的权重。以下是具体的数学推导。
数学原理：定义似然函数
首先我们要定义一个概率质量函数：
$$ P(y=1|x;\theta)=h_{\theta}(x) $$
$$ P(y=0|x;\theta)=1-h_{\theta}(x) $$
Note:上边的第一个公式的等式左边读作“在已知参数\(theta\)和特征向量x的条件下，y=1的概率”，我们使用右边的假设函数来计算这个概率。
这两个公式也可以合并成一个：
$$ P(y|x;\theta)=h_{\theta}(x)^{y}(1-h_{\theta}(x)^{1-y} $$
下表显示了我们假设函数的输出与预测概率的关系。
| h(x) | y=0 | y=1 |
|---|---|---|
| 0.25 | 0.75 | 0.25 |
| 0.5 | 0.5 | 0.5 |
| 0.75 | 0.25 | 0.75 |
| 0.9 | 0.1 | 0.9 |
很明显，我们需要最大化等式的右边，因此我们需要引申出似然函数的概念。
以我的理解，似然函数就是表示给定特征向量\(widehat{x}\)能正确预测出y的可能性。当然，从统计的角度讲，似然和概率是不一样的，这里给出了详细的解释。
现在我们可以详细的谈一谈似然函数的事情啦！在这个项目中，我们应用在训练集中的似然函数很简单，我们将每一种情况的可能性乘起来，这样就得到了我们模型正确预测的累积似然，公式如下。
$$ L(\theta) = \prod_{n}^{i=1}p(y_{i}|x_{i};\theta) $$
更通用的，可以用以下的公式表示：
$$ L(\theta) = h_{\theta}(x_{i})^{y_{i}}(1-h_{\theta}(x_{i}))^{1-y_{i}} $$
显然，我们将n个概率相乘了，而相乘的结果必然是处于[0,1]区间的。n表示样本的个数，所以最终的概率可能是个\(10^{-n}\)数量级的一个概率。这简直就是一个糟糕的消息，因为计算的过程中可能会损失精度，而且如果用python来计算的话，我们最终看到的结果很有可能被四舍五入成为0。
我们的解决办法是对似然函数求对数，这样一来，连乘的形式就变成了连加，公式如下：
$$ \ell (\theta)=log(L(\theta)) $$
$$ \ell (\theta)=\sum_{i=1}^{n}y_{i}log(h_{\theta}(x_{i}))+(1-y_{i})log(1-h_{\theta}) $$
Note:这里给出对数函数的基本运算法则
$$ log(xy)=log(x)+log(y) $$
$$ log(x)^{n}=nlog(x) $$
这个函数被称为我们假设函数的对数似然函数，因为对数函数具有单调性，因此不影响我们下面的计算。
记住，当数据非常小的时候，我们之前假设的函数会变得非常不精确（因为数值太小，连乘运算之后误差很大），而对数似然就显现出良好的精度，对数似然的图如下所示，而接下来我们要做的就是最大化对数似然函数。
log_likelihood_normalized
数学原理：单变量牛顿法
在我们求对数似然函数的最大值之前，我们先来介绍一下牛顿法
牛顿法使用迭代的思想，它是求多项式函数根的一种算法。在简单的单变量情况下，牛顿法的实现步骤如下：
- 在\(p(x_{n},y_{n})\)处求出f(x)的切线\($y={f}'(x_{n})(x-x_{n})+f(x_{n})$\)
- 计算出切线在x轴上的截距
$$ 0={f}'(x_{n})(x-x_{n})+f(x_{n}) $$
$$ -f(x_{n})={f}'(x_{n})(x-x_{n}) $$
$$ x_{n+1}=x_{n}-\frac{f(n)}{{f}'(n)} $$
- 求出该位置的y值
$$ y_{n+1}=f(x_{n+1}) $$
- 如果\(y_{n+1}approx y{n}\)恭喜你，我们的函数收敛了
- 否则继续更新\(p(x_{n},y_{n})\)
\(x=x_{n+1},y=y_{n+1}\),返回第一步迭代
下面我们使用维基百科上的gif图来更直观的理解这一过程：
NewtonIteration_Ani
如果你理解了以上的步骤，我们可以把这一过程归纳为：
$$ x_{n+1}=x_{n}-\frac{f(n)}{{f}'(n)} $$
（直到\(x_{n}-x_{n+1}approx 0\)）
任何高中学习过微积分的人都能理解上面的内容。但是我们如何推广到多元的，也就是“n维”的情况呢？
数学原理：n维牛顿法
在多维的情况下，我们使用偏导数的向量，也就是梯度来代替一维情况下的导数。
如果你对梯度的概念有疑问，那么请看这篇帖子复习一下。
相应的，我们的更新规则也发生了改变，具体的公式如下所示：
$$ x_{n+1}=x_{n}-\frac{f(n)}{\triangledown {f}'(n)} $$
Note从一维推广到n维，其实就是将导数换成了梯度，也就是将\({f}'(n)$变成triangledown {f}'(n)\)
数学原理：用牛顿法最大化对数似然
还记得我们之前提到的，我们的目的是最大化对数似然函数，也就是最大化\(ell(theta)\)函数，因为log函数是单调上涨的，因此，只要保证对数似然函数是最大的，那么我们假设的函数\(H_{theta}(x)\)也就达到了最大值。
通常，求最大值的方式是对函数求偏导数，并让其等于0，求解\(theta_{1}和theta_{2}\)来算出极值,极值也就是我们要找的最大值了。
Note:我们的log函数是递增的，因此我们只有一个极值点，而本方法就是求这个极值点的唯一方法。
这个听起来是不是非常的熟悉，接下来，我们就要求让偏导数为0的\(theta_{1},theta_{2}\)两个参数的值了，这正是牛顿法解决的问题，让我们一起再来回顾下：
$$ x_{n+1}=x_{n}-\frac{f(n)}{\triangledown {f}'(n)} $$
我们用梯度来替换$f(x_{n})$,$\triangledown \l(\theta)$，这样我们的公式就会变成
$$ x_{n+1}=x_{n}-\frac{\triangledown \ell(\theta)}{?} $$
?表示什么？直觉告诉我们，我们需要取梯度向量的导数，就像我们以前取f(x)的导数一样。
接下来我们就需要Hessian矩阵
数学原理：Hessian矩阵
根据以前学的多元微积分的知识，我们应该知道，为了找到函数的二阶偏导数，我们取每一阶偏导数的偏导数。如果我们有n个参数，然后我们就会得到n×n阶导数。因此，Hessian矩阵是n*n阶偏导数的平方矩阵。
在我们的案例中，我们只有两个参数\(theta_{1},theta_{2}\)，我们的Hessian矩阵就是这样:
$$ H_{\ell(\theta)}=\begin{bmatrix} \frac{\partial^{2} \ell }{\partial {\theta_{1}}^{2}}& \frac{\partial^{2} \ell }{\partial \theta_{1}\partial \theta_{2}} \\ \frac{\partial^{2} \ell }{\partial \theta_{2}\partial \theta_{1}} & \frac{\partial^{2} \ell }{\partial {\theta_{2}}^{2}} \end{bmatrix}\quad $$
数学原理：整合
在我们的通过将Hessian矩阵替换进牛顿法的更新步骤，那么我们的公式就变成：
$$ \theta_{n+1}=\theta_{n}+H_{\ell(\theta)}^{-1}\bigtriangledown \ell(\theta) $$
Note在这里我们取Hessian矩阵的逆，而不是取它的倒数，因为它是一个矩阵。
为了简洁起见，这篇文章省去了梯度和Hessian矩阵的数学推导。如果你想了解下面推导的资源可以在：
- 对数似然函数的梯度的推导，参考Andrew Ng的讲义17-18页。
- Hessian矩阵的推导是非常直白的，当你计算出梯度的时候，求Sigmoid函数导数也可以参考Andrew Ng的讲义
\(ell(theta)\)的梯度是
$$ \begin{bmatrix} \sum_{i=1}^{n}(y_{i} - h_{\theta}(x_{i}))x_{i} \\ \sum_{i=1}^{n}(y_{i} - h_{\theta}(x_{i})) \end{bmatrix} $$
\(ell(theta)\)的Hessian矩阵是
$$ H_{\ell(\theta)}=\begin{bmatrix} \sum_{i=1}^{n}h_{\theta}(x_{i})(1-h_{\theta}(x_{i}))\theta_{1}\theta_{1} & \sum_{i=1}^{n}h_{\theta}(x_{i})(1-h_{\theta}(x_{i}))\theta_{1} \ \sum_{i=1}^{n}h_{\theta}(x_{i})(1-h_{\theta}(x_{i}))\theta_{1} & \sum_{i=1}^{n}h_{\theta}(x_{i})(1-h_{\theta}(x_{i}) \end{bmatrix}\quad $$
其中\(h_{theta}(x) = frac{1}{1 + e^{-z}} z = theta_{1}x + theta_{2}\)
牛顿法实现
我们首先定义我们的假设函数为Sigmoid函数。
def sigmoid(x, Θ_1, Θ_2):
z = (Θ_1*x + Θ_2).astype("float_")
return 1.0 / (1.0 + np.exp(-z)) 然后我们定义\(ell(theta)\),也就是我们的对数似然函数
def log_likelihood(x, y, Θ_1, Θ_2):
sigmoid_probs = sigmoid(x, Θ_1, Θ_2)
return np.sum(y * np.log(sigmoid_probs)
+ (1 - y) * np.log(1 - sigmoid_probs))最后，我们实现了对数似然函数的梯度和Hessian矩阵。
def gradient(x, y, Θ_1, Θ_2):
sigmoid_probs = sigmoid(x, Θ_1, Θ_2)
return np.array([[np.sum((y - sigmoid_probs) * x),
np.sum((y - sigmoid_probs) * 1)]])
def hessian(x, y, Θ_1, Θ_2):
sigmoid_probs = sigmoid(x, Θ_1, Θ_2)
d1 = np.sum((sigmoid_probs * (1 - sigmoid_probs)) * x * x)
d2 = np.sum((sigmoid_probs * (1 - sigmoid_probs)) * x * 1)
d3 = np.sum((sigmoid_probs * (1 - sigmoid_probs)) * 1 * 1)
H = np.array([[d1, d2],[d2, d3]])
return H我们在这4个函数中实现了所有的数学公式，我们创建了一个While循环，并使用牛顿的方法迭代，直到收敛到最大值为止。
def newtons_method(x, y):
"""
:param x (np.array(float)): Vector of Boston House Values in dollars
:param y (np.array(boolean)): Vector of Bools indicting if house has > 2 bedrooms:
:returns: np.array of logreg's parameters after convergence, [Θ_1, Θ_2]
"""
# Initialize log_likelihood & parameters
Θ_1 = 15.1
Θ_2 = -.4 # The intercept term
Δl = np.Infinity
l = log_likelihood(x, y, Θ_1, Θ_2)
# Convergence Conditions
δ = .0000000001
max_iterations = 15
i = 0
while abs(Δl) > δ and i < max_iterations:
i += 1
g = gradient(x, y, Θ_1, Θ_2)
hess = hessian(x, y, Θ_1, Θ_2)
H_inv = np.linalg.inv(hess)
# @ is syntactic sugar for np.dot(H_inv, g.T)¹
Δ = H_inv @ g.T
ΔΘ_1 = Δ[0][0]
ΔΘ_2 = Δ[1][0]
# Perform our update step
Θ_1 += ΔΘ_1
Θ_2 += ΔΘ_2
# Update the log-likelihood at each iteration
l_new = log_likelihood(x, y, Θ_1, Θ_2)
Δl = l - l_new
l = l_new
return np.array([Θ_1, Θ_2]) 可视化：牛顿法
让我们用可视化的方式看一看，对数似然函数的表面上绘制牛顿法的每一次迭代会发生什么：
newtons_method_ll
Note:第一次迭代我们用红色来标识，第二次迭代我们用橘黄色来标识...最后一次迭代的颜色是紫色
newtons_method_max
这个图形可以很直观的看到我们的“紫色”点是最大值的点，这就说明，我们方法已经收敛了！
可视化：我们的方案
我们拥有的数据集其实是一个一维数据集，通常来说，一维数据的可视化是将这些点放置在线上，然后在线上做一些划分的边界。可问题是我们数据集里的点太多，如果画在一条线上，数据点会串在一起。
所以，我就将这些点沿着y轴给拉了出来，并用不同的颜色表示不同的类别。另外，我还在图中做了三条辅助线来按百分比划分这些点。
south_boston_bedroom_logreg
结论
本篇已经涵盖了许多新的主题，比如Hessian矩阵、对数似然函数、Sigmoid函数。结合这些方法，我们已经能够实现牛顿的方法来解决logistic回归。
尽管我们的方法非常的具体，但是需要担心的是我们的牛顿法也有发散的可能，当然，这个超出了本文的讨论范围，如果感兴趣，可以查阅一些相关的资料。
作者介绍: Sean Harrington 全栈工程师、数据科学家，熟悉Python, Javascript, C, C++, Bash等语言
个人主页：http://thelaziestprogrammer.com/author/sharrington/
以上为译文
本文由北邮@爱可可-爱生活 老师推荐，阿里云云栖社区组织翻译。
文章原标题《Solving Logistic Regression with Newton's Method》，作者：Sean Harrington，译者：爱小乖，审校：6816816151。
文章为简译，更为详细的内容，请查看原文。
本文由用户为个人学习及研究之目的自行翻译发表，如发现侵犯原作者的版权，请与社区联系处理yqgroup@service.aliyun.com
