ML02-逻辑回归（上部分）

2018-12-07 本文已影响57人杨强AT南京

本文主题-逻辑回归（上部分）：

逻辑回归的应用背景

逻辑回归的数学基础

逻辑回归的模型与推导

逻辑回归算法推导

梯度下降算法

逻辑回归梯度下降算法实现

一. 线性回归在分类上存在的问题

线性回归的条件

在上述例子中，我们看到线性回归可以做预测，也可以做分类，一切结果看起来都是很完美。
线性回归模型：

$y^{(i)}=x^{(i)}W + b^{(i)} + \epsilon^{(i)}$

其中 $\epsilon^{(i)}$ 的推导建立在一个假设之上：服从正态分布（Gaussian Distribution）

$p(x)=\dfrac{1}{\sqrt{2\pi} \sigma}exp(-\dfrac{(x- \mu)^2}{2 \sigma^2})$ 其中：μ表示期望（均数），σ表示标准差，σ2表示方差

记正态分布为 $N(μ，σ)$ ，标准正态分布是 $N( 0，1 )$

下面可视化正态分布：

% matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
#  绘图的坐标轴
figure=plt.figure('正态分布可视化', figsize=(6, 4))
ax = figure.add_axes([0.1, 0.1, 0.8, 0.8], label='正态分布函数曲线')
ax.set_xlabel('x值')
ax.set_ylabel('随机值：y')

# 设置绘制曲线范围
ax.set_xlim(left=-10,right=10)    # x取值范围
ax.set_ylim(bottom=0, top=0.5)    # 概率在[0, 1) 之间

# x取值范围
x=np.linspace( -10 ,10, 100,  dtype=np.float32 )

# 系数-方差
sigma=1
# 系数-均值
mu=0

# 正态分布常数系数
coefficient = 1.0 / ( np.sqrt( 2 * np.pi) * sigma )
# 正态分布指数
exponent = -(x-mu)**2/(2*sigma**2)

y=coefficient * np.exp( exponent )

ax.plot(x, y,color='r',label='正态分布曲线')
ax.legend()

figure.show(warn=False)

正态分布曲线

正态分布函数（概率密度函数）从正无穷到负无穷积分的概率为1。即频率的总和为100%
随机正态分布概率，就是其中一段区域的积分。

一个不适合使用线性回归分类的例子

使用随机方式构造一个数据集

% matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

cls1_data = np.random.uniform(low=10000,high=100000,size=50)
cls1_result = np.ones(shape=(50),dtype=np.int32)

cls2_data = np.random.uniform(low=0,high=30000,size=50)
cls2_result = np.zeros(shape=(50),dtype=np.int32)

# 绘制样本
figure=plt.figure('用于线性回归的分类样本', figsize=(6, 4))
ax = figure.add_axes([0.1, 0.1, 0.8, 0.8], label='2类样本')
ax.set_xlabel('取值区间')
ax.set_ylabel('分类：0与1')

# 设置绘制曲线范围
ax.set_xlim(left=0,right=100000)    # x取值范围
ax.set_ylim(bottom=-0.3, top=1.5)    # 概率在[0, 1) 之间

ax.scatter(cls1_data,cls1_result,color='r',marker='.')

ax.scatter(cls2_data,cls2_result,color='b',marker='.')

figure.show(warn=False)

随机数据集可视化

使用线性回归分类

# 采用sklearn模块来实现线性回归
from sklearn.linear_model import *
import numpy as np
x=np.hstack((cls1_data,cls2_data))
y=np.hstack((cls1_result,cls2_result))
x=x.reshape(-1, 1)

regression=LinearRegression()
# 训练
# 数据格式整理
# Reshape your data either using array.reshape(-1, 1) if your data has a single feature or array.reshape(1, -1) if it contains a single sample.
# 如果是一维特征，使用reshape处理单特征形状（-1, 1），如果只有一个样本，形状reshape(1,-1)
regression.fit(x, y)
# 测试
print('评估：',regression.score(x, y))
# 斜率
print('斜率：',regression.coef_)
# 截距
print('截距：',regression.intercept_ )

# 使用训练数据统计分类正确率
pre=regression.predict(x)
print(pre)
cls_a=pre[0:50]
cls_a=cls_a>0.5
a_num=cls_a.sum()

# 第二类：后50
cls_b=pre[50:]
cls_b=cls_b<0.5
b_num=cls_b.sum()
# 统计正确率
a_pct=cls_a.mean()
b_pct=cls_b.mean()
print("A类识别正确数：%d，识别正确率：%5.2f%%"%(a_num,a_pct*100))
print("B类识别正确数：%d，识别正确率：%5.2f%%"%(b_num,b_pct*100))
# 可视化分类情况
figure=plt.figure('分类情况', figsize=(6, 4))
ax = figure.add_axes([0.1, 0.1, 0.8, 0.8], label='2类样本')
ax.set_xlabel('取值区间')
ax.set_ylabel('分类：0与1')

# 设置绘制曲线范围
ax.set_xlim(left=0,right=100000)    # x取值范围
ax.set_ylim(bottom=-0.3, top=1.5)    # 概率在[0, 1) 之间

ax.scatter(cls1_data,cls1_result,color='r',marker='.')

ax.scatter(cls2_data,cls2_result,color='b',marker='.')

ax.plot(x, regression.predict(x),color='y', label='回归直线')

ax.legend()

figure.show(warn=False)

评估： 0.555891654577714
斜率： [1.34255212e-05]
截距： 0.01655814267519684
[0.84828292 0.58331511 1.11550652 1.05660159 0.5944528  1.15643254
 1.14064579 0.94602928 0.48950987 0.48060621 0.5319122  0.71212052
 0.50981994 0.7330047  0.62845884 1.21112727 0.54341076 0.21257363
 0.26175963 1.02209569 1.23536988 0.53998925 0.98959512 0.17659679
 0.4907314  0.64503974 0.18405636 0.76383753 0.97628677 1.35292452
 0.20352834 0.98268824 1.06256155 1.32459732 0.99089918 0.46460692
 0.98906349 0.91117129 0.2198635  0.59830641 0.8756055  1.2736329
 0.83592414 1.06118727 0.65489012 0.8660766  0.61240252 1.33904091
 0.47394198 1.02521    0.05723705 0.37106916 0.33210881 0.33587503
 0.19013076 0.24586235 0.36920796 0.14182335 0.37130829 0.25797109
 0.2270329  0.27208044 0.15928228 0.1021966  0.17905334 0.22717218
 0.0948184  0.13249571 0.38909531 0.05946123 0.35638356 0.1328463
 0.0309184  0.41697969 0.08458613 0.11541662 0.08907712 0.02159664
 0.273264   0.13543039 0.27976214 0.38864573 0.02586998 0.35168836
 0.10483192 0.32039087 0.09642115 0.30624559 0.35838634 0.35280365
 0.15571477 0.38641816 0.11186056 0.14146531 0.38562691 0.06455996
 0.23412796 0.22797953 0.38572513 0.25240354]
A类识别正确数：39，识别正确率：78.00%
B类识别正确数：50，识别正确率：100.00%

随机数据集的拟合可视化

明显的两类样本是可分的，使用线性回归得到的分类效果明显存在较大误差。这是因为我们的样本违背了线性回归的三个假设：
|-自变量与因变量是线性关系；
|-自变量与误差项相互独立；
|-误差项服从正态分布；

二、逻辑回归的数学基础

从分类开始谈起
某个样本属于A类还是B类，从结果上讲就是值为0，还是值为1。但影响这个分类的是由一些因素决定的。我们从数学角度可以使用向量表示这些因素（这些因素就影响某个样本属于A类还是B类）：

$x=(x_1,x_2,\dots,x_n)$
其中 $x$ 就是表示一个样本，样本 $x$ 具有n个影响分类的特征。如果考虑偏置项，则可以增加一个份量1。
$x=(x_1,x_2,\dots,x_n,1)$
建立分类的逻辑模型
我们假设有两套标准判定样本所属的分类，使用数学函数表示如下：
|- $y_A=f(x)$ 样本x属于A的可能性；
|- $y_B=g(x)$ 样本x属于B的可能性；

这样我们就可以建立一个分类模型：

$y=\begin{cases} 1,\quad y_A>y_B\\ 0,\quad y_A\leqslant y_B\\ \end{cases}$

当 $y=1$ ，则样本 $x$ 属于A类；当 $y=0$ ，则样本 $x$ 属于B类；

可以把上述模型表示为：

$y=\begin{cases} 1,\quad y_A-y_B>0\\ 0,\quad 其他\\ \end{cases}$

分类逻辑模型的概率分析基础

如果假设 $y_A,y_B$ 与 $x$ 是线性关系，同时考虑 $y_A,y_B$ 的误差都独立服从正态分布，可以把 $y_A,y_B$ 表示如下：

|- $y_A=xW_A+\epsilon_A$
|- $y_B=xW_B+\epsilon_B$

其中 $\epsilon_A,\epsilon_B$ 是服从独立分布的误差项，可以假设服从正态分布。
记 $z=y_A-y_B$ ，则：
|- $z=x(W_A-W_B)+(\epsilon_A-\epsilon_B)$

从而分类逻辑模型可以表示如下：
|- $y=\begin{cases} 1,\quad z>0\\ 0,\quad 其他\\ \end{cases}$
其中 $W=W_A-W_B,\epsilon=\epsilon_A-\epsilon_B$
$z=xW+\epsilon$

则样本X属于A的概率可以表示为：
|- $P(y=1)=P(z>0)=P(xW+\epsilon>0)=P(\epsilon>-xW)$
从正态分布可以继续推导：
|- $P(y=1)=1-P(\epsilon\leqslant-xW)=1-F_{\epsilon}(-xW)$
其中 $F_{\epsilon}$ 是变量 $\epsilon$ 的累积分布函数； $P(y=1)$ 表示样本属于A的概率
probit模型
上述推导的公式在学术上称为probit模型，建立的回归模型也称proit回归。

$P(y=1)=1-F_{\epsilon}(-xW)$

$F_{\epsilon}(-xW)$ 是正态分布函数的累积函数，上述累积分布函数，在服从正态分布的时候比较麻烦，因为正态分布的累积函数还没有解析表达式能够表达。从而其参数估计非常麻烦，如果需要应用，则需要简化（做近似处理）。
为了解决正态分布累积函数的问题，正态分布的累积函数的计算居然是通过查表的形式提供运算结果，大家如果想不起，可以查阅自己的高中数据或者大学数学书。
正态分布的近似处理
正态分布表示：
$p(x)=\dfrac{1}{\sqrt{2\pi} \sigma}exp(-\dfrac{(x-\mu)^2}{2 \sigma^2})$ 其中：μ表示期望（均数），σ表示标准差，σ2表示方差

记正态分布为 $N(μ，σ)$ ，标准正态分布是 $N( 0，1 )$

我们可以简化下正态分布累积函数的表示：
$F_{\epsilon}(x)=\phi(\dfrac{x-\mu}{\sigma})$
因为 $\dfrac{x-\mu}{\sigma}$ 是线性，所以只需要考虑标准正态分布。
在数学上存在一个逻辑分布，与正态分布非常相似，
逻辑分布与标准正态分布
下面使用可视化来认识下逻辑分布函数与正态分布函数的近似度。
|- 正态分布： $p(x)=\dfrac{1}{\sqrt{2\pi}\sigma}exp(-\dfrac{(x-\mu)^2}{2\sigma^2})$
|- 逻辑分布： $p(x)=\dfrac{e^{-x}}{(1+e^{-x})^2}$

逻辑分布函数与正态分布的区别就在于：逻辑分布有累积函数，其累积函数如下：
|- 逻辑分布累积函数： $S(x)=\dfrac{1}{1+e^{-x}}$

下面是逻辑分布函数，正态分布函数的比较：

% matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
#  绘图的坐标轴
figure=plt.figure('正态分布与逻辑分布', figsize=(10, 4))
ax = figure.add_axes([0.1, 0.1, 0.8, 0.8], label='逻辑分布与正态分布')
ax.set_xlabel('x值')
ax.set_ylabel('y值')

# 设置绘制曲线范围
ax.set_xlim(left=-10,right=10)    # x取值范围
ax.set_ylim(bottom=0, top=0.5)    # 概率在[0, 1) 之间

# 标准正态分布函数--------------------
# x取值范围
x_n=np.linspace( -10 ,10, 100,  dtype=np.float32 )
# 系数-方差
sigma=1
# 系数-均值
mu=0
# 正态分布常数系数
coefficient = 1.0 / ( np.sqrt( 2 * np.pi) * sigma )
# 正态分布指数
exponent = -(x_n-mu)**2/(2*sigma**2)
y_n=coefficient * np.exp( exponent )
ax.plot(x_n, y_n,color='r',label='正态分布')
# 逻辑分布函数--------------------
x_l=np.linspace( -10 ,10, 100,  dtype=np.float32 )
y_l=np.exp( -x_l ) / ( 1 + np.exp( -x_l ) )**2
ax.plot(x_l, y_l,color='b',label='逻辑分布')

ax.legend()

figure.show(warn=False)

正态分布与逻辑分布比较

逻辑分布与标准正态分布的累积函数
标准正态分布没有累积函数，所以我们采用scipy的积分函数来绘制。
使用积分表示正态分布累积函数的公式如下：
|- $F_{\epsilon}(x)=\int_{-\infty}^{x}\dfrac{1}{\sqrt{2\pi}\sigma}exp(-\dfrac{(t-\mu)^2}{2\sigma^2})\mathrm{d} t$ ，

% matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import scipy.integrate as integrate

#  定义正态函数，用于积分运算
def  normal(t,m,s):
    # 系数-方差
    sigma=s
    # 系数-均值
    mu=m
    # 正态分布常数系数
    coefficient = 1.0 / ( np.sqrt( 2 * np.pi) * sigma )
    # 正态分布指数
    exponent = -(t-mu)**2/(2*sigma**2)
    re=coefficient * np.exp( exponent )
    return re

def cumulative(x):
    re=integrate.quad(     
        normal,         # 积分函数
        -np.inf,         # 积分下限
        x,              # 积分上限
        args=(0.0, 1.0)   # 传递给函数的参数（除第一个参数外，按照顺序来）
    )
    return re[0]     #第一个是积分，第二个是误差上限

#  绘图的坐标轴
figure=plt.figure('正态分布与逻辑分布', figsize=(10, 4))
ax = figure.add_axes([0.1, 0.1, 0.8, 0.8], label='逻辑分布与正态分布')
ax.set_xlabel('x值')
ax.set_ylabel('y值')

# 设置绘制曲线范围
ax.set_xlim(left=-10,right=10)    # x取值范围
ax.set_ylim(bottom=-0.1, top=1.1)    # 概率在[0, 1) 之间

# 标准正态分布函数--------------------
# x取值范围
x_n=np.linspace( -10 ,10, 100,  dtype=np.float32 )
y_n=[cumulative(x) for  x in x_n ]
ax.plot(x_n, y_n,color='r',label='正态分布累积函数')
# 逻辑分布函数--------------------
x_l=np.linspace( -10 ,10, 100,  dtype=np.float32 )
y_l=1.0 / ( 1 + np.exp( -x_l ) )
ax.plot(x_l, y_l,color='b',label='逻辑分布累积函数')

ax.legend()
ax.grid(b=True)

figure.show(warn=False)

逻辑分布与正态分布的累积函数曲线

最佳的逻辑分布
如果逻辑分布参数做一些调整，可以最佳接近。下面是最佳逻辑分布累积函数表示：
|- $S(x)=\dfrac{1}{1+e^{-1.702x}}$
下面是可视化后的直观效果。

% matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import scipy.integrate as integrate

#  定义正态函数，用于积分运算
def  normal(t,m,s):
    # 系数-方差
    sigma=s
    # 系数-均值
    mu=m
    # 正态分布常数系数
    coefficient = 1.0 / ( np.sqrt( 2 * np.pi) * sigma )
    # 正态分布指数
    exponent = -(t-mu)**2/(2*sigma**2)
    re=coefficient * np.exp( exponent )
    return re

def cumulative(x):
    re=integrate.quad(     
        normal,         # 积分函数
        -np.inf,         # 积分下限
        x,              # 积分上限
        args=(0.0, 1.0)   # 传递给函数的参数（除第一个参数外，按照顺序来）
    )
    return re[0]     #第一个是积分，第二个是误差上限

#  绘图的坐标轴
figure=plt.figure('正态分布与逻辑分布', figsize=(10, 4))
ax = figure.add_axes([0.1, 0.1, 0.8, 0.8], label='逻辑分布与正态分布')
ax.set_xlabel('x值')
ax.set_ylabel('y值')

# 设置绘制曲线范围
ax.set_xlim(left=-10,right=10)    # x取值范围
ax.set_ylim(bottom=-0.1, top=1.1)    # 概率在[0, 1) 之间

# 标准正态分布函数--------------------
# x取值范围
x_n=np.linspace( -10 ,10, 100,  dtype=np.float32 )
y_n=[cumulative(x) for  x in x_n ]
ax.plot(x_n, y_n,color='r',label='正态分布累积函数 $F_{\epsilon}(x)=\int_{-\infty}^{x}\dfrac{1}{\sqrt{2\pi}\sigma}exp(-\dfrac{(t-\mu)^2}{2\sigma^2})\mathrm{d} t$')
# 逻辑分布函数--------------------
x_l=np.linspace( -10 ,10, 100,  dtype=np.float32 )
y_l=1.0 / ( 1 + np.exp( -1.702*x_l ) )
ax.plot(x_l, y_l,color='b',label='逻辑分布累积函数  $S(x)=\dfrac{1}{1+e^{-1.702x}}$ ')

ax.legend()
ax.grid(b=True)

figure.show(warn=False)

逻辑分布与正态分布累积函数的最佳逼近1.702

9.sigmoid函数
$S(x)=\dfrac{1}{1+e^{-1.702x}}$ 就是鼎鼎大名的sigmoid函数，该函数具备许多良好的性质。尤其是深度学习中的大杀器！经常会用到。
函数的特点：
|- 函数的导数可以用自身表示： $S^\prime(x)=S(x)(1-S(x))$
|- 连续处处可微
|- 值在[0,1)范围
|- 单调递增
|- 非线性
|- 平滑性
|- 原点附近近似identity（f(x)≈x）或者先线性性。
函数的缺点：
|- 运算量大（后面使用梯度下降算法可以提现出来，尤其误差传递的时候）
|- 容易出现梯度消失的情况，从而不利于样本训练，甚至完不成梯度训练。（大家可以从逻辑分布函数看得出来，sigmoid函数的导数就是逻辑分布函数，逻辑分布函数从0开始，就逐步趋向于0，容易产生梯度消失）

三、逻辑回归模型与模型推导

逻辑回归模型
从上述逻辑分布与正态分布的相似性，可以得到逻辑回归模型:

|- （1）proit回归模型： $p(y=1)=1-F_{\epsilon}(-xW)$
|- （2）正态分布简写模型： $F_{\epsilon}(x)=\phi(\dfrac{x-\mu}{\sigma}) \qquad\to\qquad F_{\epsilon}(\sigma x+\mu)=\phi(x) \approx S(1.702x)$

在porit回归模型基础上，简单推导，可以定义如下模型：
|- $p(y=1) =\dfrac{1}{1+e^{-xW}}$
上上述公式就是逻辑回归模型（省略了具体的推导还有一部分假设）。
发生比
我们可以定义 $\dfrac{p(y=1)}{1-p(y=1)}$ 为发生比。然后简单取自然对数，得到如下公式：
$ln \dfrac{p(y=1)}{1-p(y=1)}=xW$

逻辑回归是建立在一个假设之上：事件发生比实际上是一个线性模型。
似然函数
逻辑回归的算法还是要从模型参数的估计开始。
从统计来说，数据出现的概率越大越好，这就是最大似然估计法，由于数据的取值是离散的0，1，从而数据取值的概率分布函数可以表示为(这里用到了二项分布概率累积函数的计算方法)：
|- $p(y)={p(y=1)}^{y}{p(y=0)}^{1-y}$
取自然对数：
|- $ln\ p(y)=y\ ln\ p(y=1) + (1-y)\ ln\ (1- p(y=1))$

考虑所有样本，而且每个样本独立，所以可以对所有样本形成下面公式：
|- $p(Y)=\prod\limits_ip(y_i)$ 其中 $i$ 表示遍历所有样本， $y_i$ 表示其中第 $i$ 个样本。
把所有样本的概率分布函数展开：
|- $p(Y)=\prod\limits_i {p(y_i=1)}^{y_i}{p(y_i=0)}^{1-y_i}$
替换为概率分布函数：sigmoid函数，可以得到如下公式：
|- $p(Y)=\prod\limits_i {h(X_i)}^{y_i}{(1-h(X_i))}^{1-y_i}$
其中 $h(X_i)$ 表示sigmoid函数 $h(X_i) =\dfrac{1}{1+e^{-X_iW}}$ ， $X_i$ 表示第 $i$ 个样本。
上述公式就是W的似然函数，可以更加明确的按照条件概率方式记为：
|- $L(W)=p(Y|W)=\prod\limits_i {h(X_i)}^{y_i}{(1-h(X_i))}^{1-y_i}$
参数估计公式推导
使用最大似然估计法，可以得到W的参数估计公式（求对数）：
|- $l(W)=ln(L(W))=\sum\limits_{i}y_i\ ln(h(X_i)) +(1-y_i)\ln(1-h(X_i))$
机器学习的损失函数定义
由于使似然函数值最大，相反考虑，定义损失函数如下：
$Loss= - ln(L(W))=-l(W)$
这样就使得问题编程求最小值，我们可以定义成逻辑回归模型的损失函数。使用算法找到 $Loss$ 的最小值即可。

四、逻辑回归模型定义

决策函数
$y_i=S(X_iW)$

其中 $S(x)= \dfrac{1}{1+e^{-x}}$ 是sigmoid函数
其中 $X_i=(x_1,x_2,\dots,x_n,1)$ 是第 $i$ 个样本数据。
其中 $W$ 是样本的线性系数，也可以理解成样本特征的重要性加权。
机器学习中的误差定义
$Loss=J(W)=-l(W)=-ln(L(W))=-\sum\limits_{i}y_i\ ln(h(X_i)) +(1-y_i)\ln(1-h(X_i))$
上述公式被称为交叉熵（Cross Entropy）
其中 $h(X_i)=S(X_iW)= \dfrac{1}{1+e^{-X_iW}}$

五、逻辑回归的算法推导

逻辑回归最终还是找到一个W，使得分类误差最小，也就是使得损失函数最小，所以逻辑回归属于于广义线性回归的一种。
逻辑回归算法也是找到一个W，使得对所有样本 $J(W)$ 函数的值最小，由于考虑问题的思路不同，这里的 $J(W)$ 不再像线性回归函数那样：通过推导直接得到 $W$ 的求解公式。
下面我们尝试使用求导的方式，来推导能否求解&J(W)&的最小值，从而得到W的解。

推导前的一个已知知识：sigmoid函数（也称S曲线函数）有一个显著特性就是其导数可以使用自身表示。
|- $S^\prime(x)=S(x)(1-S(x))$

求偏导数

$\nabla_W J(W)=\nabla_W(-\sum\limits_{i}y_i\ ln(h(X_i)) +(1-y_i)\ln(1-h(X_i)))$

$\nabla_W J(W)=\dfrac{\partial J(W)}{\partial W}=-\sum\limits_i (\dfrac{y_i}{h(X_i)} + \dfrac{1-y_i}{1-h(Xi)} ) \dfrac{\partial h(X_i)}{\partial W}$

$\nabla_W J(W)=-\sum\limits_i (\dfrac{y_i}{h(X_i)} + \dfrac{1-y_i}{1-h(Xi)} ) (h(X_i)(1-h(X_i))) \dfrac{\partial X_iW)}{\partial W}$

上述推导基于一个数学知识：链式求偏导公式。下面继续推导：

$\nabla_W J(W)=-\sum\limits_i (y_i(1-h(X_i)) - (1-y_i)h(X_i) ) X_i$

上述公式注意最后的推导，我们直接采用了向量求导，直接得到 $X_i$ 。
下面对前面几项做多项式拆分，得到如下公式：

$\nabla_W J(W)=-\sum\limits_i (y_i- y_i\ h(X_i) + y_i\ h(X_i) -h(X_i) ) X_i$
$\nabla_W J(W)=-\sum\limits_i (y_i -h(X_i) ) X_i$
$\nabla_W J(W)=\sum\limits_i (h(X_i)-y_i ) X_i$

求解W
按照原来的套路，利用极值求解W，由于其中 $h(X_i)$ 是一个指数函数，加上求和等线性因素，上述公式按照导数=0的方式求解非常麻烦，应该说无法使用公式推导求的W的值。下面我们来回顾下，数学中求极值算法。
从最大似然函数的推导结果看问题
最大似然函数
|-均方差(损失函数)
|-线性（Closed-Form）
|-最小二乘法
|-非线性
|-梯度下降法
|-牛顿迭代法
|-拟牛顿法
|-坐标下降法
|-交叉熵(损失函数)
|-非线性
|-梯度下降法
|-牛顿迭代法
|-拟牛顿法
|-坐标下降法

常见的损失误差：
|-1. 铰链损失（Hinge Loss）：主要用于支持向量机（SVM）中；
|-2. 互熵损失（Cross Entropy Loss，Softmax Loss ）：用于Logistic 回归与Softmax 分类中；
|-3. 平方损失（Square Loss）：主要是最小二乘法（OLS）中；
|-4. 指数损失（Exponential Loss）：主要用于Adaboost 集成学习算法中；
|-5. 其他损失（如0-1损失，绝对值损失）

六、梯度下降算法

梯度下降算法的前置条件
|-梯度下降就是让一个损失函数最小，所以梯度下降首先需要一个损失函数。
|-梯度下降使得损失函数最小，必须是因为某个量，在我们的机器学习中，这个量都是特征（ $X_i$ ）权重的影响量；在线性回归与逻辑回归中就是线性系数 $W$ 。

（1）线性回归中，损失函数是：
|- $J(W)=\dfrac{1}{2}\sum\limits_{i=1}^{m}(x^{(i)}W-y^{(i)})^2$
（2）逻辑回归中，损失函数是：
|- $J(W)=-l(W)=-ln(L(W))=-\sum\limits_{i}y_i\ ln(h(X_i)) +(1-y_i)\ln(1-h(X_i))$
|-上述公式被称为交叉熵（Cross Entropy）
|-其中 $h(X_i)=S(X_iW)= \dfrac{1}{1+e^{-X_iW}}$

梯度下降的数学极值与凸优化理论
求最小值，大家第一个想到的就是极值定理。
|-极大值与极小值就是函数在其定义域的某些局部区域所达到的相对最大值或相对最小值。当函数在其定义域的某一点的值大于该点周围任何点的值时，称函数在该点有极大值; 当函数在其定义域的某一点的值小于该点周围任何点的值时，称函数在该点有极小值。这里的极大和极小只具有局部意义。
|-极值点：极值点只能在函数不可导的点或导数为零的点上取得。

极值定理

对于一元可微函数𝒇(𝑥)，它在某点x0有极值的充分必要条件是𝒇(𝑥)在𝑥0的某邻域上一阶可导，在𝑥0处二阶可导，且𝒇'(𝑥0)=0，𝒇' '(𝑥0)≠0，那么：
       1）若𝒇' '(𝑥0)<0，则f在x0取得极大值；

       2）若𝒇' '(𝑥0)>0，则f在x0取得极小值。
求解极值，直接解𝒇 ' (𝑥0)=0方程，可以得到极值点。

最值与极值
因为极值是局部意义上的，最值是全局意义上的；在算法过程中，容易导致求出的极值在局部最大，在全局不是最大。为了保证极值就是最值，需要函数满足凸性。

凸集

对于集合 $X$ 满足：任给的 $x,y∈X$ ,总有 $λx+(1−λ)y∈X$ ,对于任意的 $λ∈(0,1)$ ，我们称X是凸集。

凸函数

对于一个函数 $f(x)$ ，满足：任给的 $x,y∈X$ （ $X$ 是函数定义域）。总有 $f(λx+(1−λ)y))≤λf(x)+(1−λ)f(y)$ ,对于任意的 $λ∈(0,1)$ ，我们称 $f(x)$ 为凸函数。
换成几何描述就是，函数的几何图形中任意两个点的连线上的点都在函数的几何图形中。

均方差函数的可视化
均方差函数本身是凸的，如果从权重考虑，则发现该函数在逻辑回归中是非凸的。
从输出的图形来看是凸的。

%  matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure('3D可视化',figsize=(8,6))
# 创建3D图形
# 方式一：
# ax = Axes3D(fig)
# 方式二：
ax = fig.add_subplot(111, projection='3d')
# J(W)对多个样本绘制
# 年龄
x = np.loadtxt('ex2x.dat')
# 身高
y = np.loadtxt('ex2y.dat')


w = np.arange(-2, 2, 0.02)
b = np.arange(-2, 2, 0.02)
w, b = np.meshgrid(w, b)    # x-y 平面的网格

z=0
for i in range(len(x)):
    z += ( ( x[i]*w + b - y[i])**2 ) 
z/=len(x)

# rstride:行之间的跨度  cstride:列之间的跨度
ax.plot_surface(w, b, z, rstride=1, cstride=1, cmap=plt.get_cmap('rainbow'),label='三维抛物面')
# 设置图像z轴的显示范围，x、y轴设置方式相同
#ax.set_xlim(-2,2)
#ax.set_ylim(-2,2)
#ax.set_zlim(0,4)

plt.show()

均方差函数的3D可视化

均方差采用S曲线函数输出的可视化
从输出的图形来看，非凸，所以在逻辑回归的损失函数，采用均方差公式训练会存在错误收敛的问题。

%  matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

fig = plt.figure('3D可视化',figsize=(8,6))
# 创建3D图形
# 方式一：
# ax = Axes3D(fig)
# 方式二：
ax = fig.add_subplot(111, projection='3d')
# J(W)对多个样本绘制
cls1_data = np.random.uniform(low=0,high=100,size=50)
cls1_result = np.ones(shape=(50),dtype=np.int32)

cls2_data = np.random.uniform(low=0,high=3,size=50)
cls2_result = np.zeros(shape=(50),dtype=np.int32)

x=np.hstack((cls1_data,cls2_data))
y=np.hstack((cls1_result,cls2_result))

w = np.arange(-1, 1, 0.04)
b = np.arange(-1, 1, 0.04)
w, b = np.meshgrid(w, b)    # x-y 平面的网格


s=lambda  p: 1.0/(1+np.exp(-p))

z=0
for i in range(len(x)):
    z += ( ( s(x[i]*w + b) - y[i] )**2 ) 
z/=len(x)
    
# rstride:行之间的跨度  cstride:列之间的跨度
ax.plot_surface(w, b, z, rstride=1, cstride=1, cmap=plt.get_cmap('rainbow'),label='三维抛物面')

#等高线绘制（直观想象3D图形）
cset = ax.contour(w, b, z, zdir='z', offset=-0.2, cmap=cm.coolwarm)
cset = ax.contour(w, b, z, zdir='x', offset=-1, cmap=cm.coolwarm)
cset = ax.contour(w, b, z, zdir='y', offset=1, cmap=cm.coolwarm)

ax.set_xlim(-1.1,1.1)
ax.set_xlim(-1.1,1.1)
ax.set_zlim(-0.2,0.6)
plt.show()

采用S曲线的均方差3D可视化

交叉熵函数的可视化
$J(W)=-ln(L(W))=-\sum\limits_{i}y_i\ ln(h(X_i)) +(1-y_i)\ln(1-h(X_i))$
|-其中 $h(X_i)=S(X_iW)= \dfrac{1}{1+e^{-X_iW}}$

在下面例子中，为了避免计算精度产生的问题，我们对样本数据同比缩小，限制在0.1之间（一种正则化思维），这样可以避免数据计算溢出的问题（无穷大与无穷小的问题）
从下面可视化效果上，看的出损失函数是凸的。实际上从数学推理上也可以证明，上面的公式非凸，这个公式是凸的。

%  matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

fig = plt.figure('3D可视化',figsize=(8,6))
# 创建3D图形
# 方式一：
# ax = Axes3D(fig)
# 方式二：
ax = fig.add_subplot(111, projection='3d')
# J(W)对多个样本绘制
cls1_data = np.random.uniform(low=0.0,high=1.0,size=50)
cls1_result = np.ones(shape=(50),dtype=np.float32)

cls2_data = np.random.uniform(low=0.01,high=0.04,size=50)
cls2_result = np.zeros(shape=(50),dtype=np.float32)

x=np.hstack((cls1_data,cls2_data))
y=np.hstack((cls1_result,cls2_result))

w = np.arange(-1, 1, 0.04,dtype=np.float64)
b = np.arange(-1, 1, 0.04,dtype=np.float64)
w, b = np.meshgrid(w, b)    # x-y 平面的网格

s=lambda  p: 1.0/(1.0+np.exp(-p))

for i in range(len(x)):
    z+=y[i]*np.log(s(x[i]*w+b))+(1.0-y[i])*(np.log(1.0-s(x[i]*w+b)))
z=-z
z/=len(x)

# rstride:行之间的跨度  cstride:列之间的跨度
ax.plot_surface(w, b, z, rstride=1, cstride=1, cmap=plt.get_cmap('rainbow'),label='三维抛物面')

#等高线绘制（直观想象3D图形）
cset = ax.contour(w, b, z, zdir='z', offset=0.5, cmap=cm.coolwarm)
cset = ax.contour(w, b, z, zdir='x', offset=-1.1, cmap=cm.coolwarm)
cset = ax.contour(w, b, z, zdir='y', offset=1.1, cmap=cm.coolwarm)

ax.set_xlim(-1.1,1.1)
ax.set_xlim(-1.1,1.1)
ax.set_zlim(0.5,1)

plt.show()

交叉熵3D可视化

交叉熵函数凸性证明
为了证明交叉熵函数是凸的，只需要求二阶导数就可以（这里用到凸优定理）：
上面已经求出一阶导数： $\nabla_W J(W)=\sum\limits_i (h(X_i)-y_i ) X_i$
|- $\dfrac{\partial^2J(w)}{\partial{W}\partial{W^T}}=\dfrac{\partial{\sum\limits_i (h(X_i)-y_i ) X_i}}{\partial{W^T}}$
|- $\dfrac{\partial^2J(w)}{\partial{W}\partial{W^T}}=\sum\limits_i {\dfrac{\partial{(h(X_i)-y_i)}}{\partial{W^T}}X_i }$
|- $\dfrac{\partial^2J(w)}{\partial{W}\partial{W^T}}=\sum\limits_i {\dfrac{\partial{(h(X_i)-y_i)}}{\partial{W^T}}X_i }$
|- $\dfrac{\partial^2J(w)}{\partial{W}\partial{W^T}}=\sum\limits_i {\dfrac{\partial{h(X_i)}}{\partial{W^T}}X_i }$
|- $\dfrac{\partial^2J(w)}{\partial{W}\partial{W^T}}=\sum\limits_i {\dfrac{e^{X_iW}{X_i}^T}{{(1+e^{X_iW})}^2}X_i}$
|- $\dfrac{\partial^2J(w)}{\partial{W}\partial{W^T}}=\sum\limits_i {\dfrac{e^{X_iW}}{{(1+e^{X_iW})}^2}{X_i}^T X_i}$
由于 $\dfrac{\partial^2J(w)}{\partial{W}\partial{W^T}}$ 是半负定矩阵，所以 $L(W)$ 是凸函数，有最小值。

最小值梯度下降求解
由于 $J(W)$ 函数我们从可视化直观观察，还是通过数学证明，都知道是凸函数，并具有最小值（注意对似然函数来讲是最大值，我们添加了一个负号- ）。由于难以使用公式符号推导求解，所以这里采用梯度下降求解。
梯度下降求解的原理
|-（1）取一个初始随机的 $W$ （这里包含偏置项）；
|-（2）根据逻辑回归模型的决策函数 $y_i=S(X_iW)$ 计算出预测值；
|-（3）如果预测值与实际值有差异，则调整 $W$ 的值；
|-（4）调整值一定要往最小值点方向调整；
|-（5）n次以后，我们可以认为 $W$ 是使得损失函数最小的稳定值。

梯度下降的方向与速度

梯度下降的方向
怎么保证调整的值是往最小值方向调整的呢？用下图的抛物线（凸函数）来充分说明梯度方向的选择（最小值情况）。
|- 为了趋向最小值点，在最小值左边应该是增加 $W$ ，在最小值右边应该是减少 $W$ （W是x轴）
|- 如果去 $J(W)$ 的导数作为 $W$ 的调整值，就可以很方便地决定方向，因为在最小值左边导数为负，载最小值右边导数为正。
|- 如果取导数作为调整值，则最小值左边减去导数值（ $W$ 增加），最小值右边加上导数值（ $W$ 减少）。

梯度下降的速度
从下图可以看见，如果 $W$ 的值调整过大，在最值点会摆动很大，导致无法准确得到最好的 $W$ ，为了控制在最小值点摆动产生的过拟合，可以对调整值加上一个参数 $\eta$ 来控制。
|-参数 $\eta$ 越大，训练速度快，取得的训练值容易过拟合。
|-参数 $\eta$ 越小，训练速度慢，容易取得最好训练值。

%  matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

fig = plt.figure('梯度下降说明',figsize=(8,6))
ax = fig.add_subplot(111,label='梯度下降说明',xlabel='X-轴（W值）',ylabel='Y-轴（损失值：$J(W)$）')
x = np.linspace(-2.0, 2.0, 20)
y = x**2
ax.plot(x, y,label='二次凸曲线',color=(1,0,0,1),marker='.',markersize=10,markerfacecolor=(0,0,1,1),markeredgecolor=(0,0,1,1))
for i in range(9):
    ax.arrow(x=x[i],y=y[i],dx=x[i+1]-x[i],dy=y[i+1]-y[i],width=0.02,color=(0,0,0,1) ,head_width=0.1,head_length=0.05)

plt.show()

梯度下降示意图

实现部分，限于篇幅，请阅读下部分