机器学习-线性回归

线性回归

基本介绍 {#sec:basic}

假定数据集合中有$m$个样本点，对于每一个样本点$x_i$，具有值$y_i$，且 x = { x i ∣ i ∈ { 0 , 1 , … , m − 1 } } x = \{x_i | i \in \{0, 1, \ldots, m-1\} \} x={xi​∣i∈{0,1,…,m−1}}与结果 y = { y i ∣ i ∈ { 0 , 1 , … , m − 1 } } y = \{y_i |i \in \{0, 1, \ldots, m-1\} \} y={yi​∣i∈{0,1,…,m−1}}呈现线性关系。
我们的目标是根据现有的$m$个样本，拟合出一条直线，即
$$y=w\ast x+b\text{\hspace{1em}(1)}$$
使得根据这条直线得到的结果与真实的结果误差最小。那么如何来衡量这个误差呢？我们引入平方损失函数，即对于样本$x_i$，我们有误差$J_i$
$$J_i=(w\ast x_i+b-y_i{)}^2$$ 因此，全部$m$个样本的平均误差为：
$$J=\frac{1}{m}\sum _{i=0}^{m-1}{J}_i$$
从数学上表述为，我们要求得$w$和$b$ ，使得$j$最小，即
$$\arg \min \frac{1}{m}\sum _{i=0}^{m-1}{J}_i$$

针对于最小化特定方程的，我们区分为几个section来介绍。

普通求导 {#par:simple}

根据多元函数求极值得方法，直接分别对$w$和$b$求偏导，且使得偏导为0，即：
$$\frac{\partial J}{\partial w}=0,\frac{\partial J}{\partial b}=0$$
只有两个未知数，两个等式，因此可以解出$w,b$。

但是，这里面假定了样本点$x_i$只有一列属性，也就是说每一个样本都可以在一个二维图像中描绘出来。当一个样本具有$n$列属性的时候，上述就会有$n+1$个等式，即：
$$\frac{\partial J}{\partial w_0}=0,\dots ,\frac{\partial J}{\partial w_{n-1}}=0,\frac{\partial J}{\partial b}=0$$
这样联立线性方程组求解非常复杂.

向量求导

当未特殊说明，向量一律为列向量。
对于具有多个属性的$x_i$，我们拟合的公式$(1)$可以转化为： $$y=w^{T}x+b$$
为了统一化，我们将$w,b$合并为向量$\theta$，并给$x_i$增加额外的一个属性，值为1.

这样我们有： $$y={\theta }^{T}x=x^T\theta$$

此时，平均误差公式可以进一步简化为： $$\begin{array}{rl}J(\theta )& =\frac{1}{m}\sum _{i=0}^{m-1}(x_i^T\theta -y_i{)}^2\\ & =\frac{1}{m}(X\theta -Y{)}^T(X\theta -Y)\end{array}$$
其中$X=\left[\begin{array}{c}{x}_0^T\\ x_1^T\\ ⋮\\ x_{m-1}^T\end{array}\right]$，且 $Y=\left[\begin{array}{c}{y}_0\\ y_1\\ ⋮\\ y_{m-1}\end{array}\right]$

我们的目标是求得最小化$J(\theta )$时的 $\theta$值，这可通过向量求导：
$$\frac{\partial J(\theta )}{\partial \theta }=0\text{\hspace{1em}(2)}$$

为了对公式$(2)$进行求解，我们首先引入几个基本向量求导公式：
$$\begin{array}{rl}\frac{\partial x^{T}a}{\partial x}& =\frac{\partial a{x}^T}{\partial x}=a\\ \frac{\partial x^{T}Bx}{\partial x}& =(B+B^T)x\end{array}$$
对于上述的基本求导公式，它的证明只有一个思想：数对向量求导，相当于此数对向量中的各个元素逐个求导。上述的两个公式都是对列向量$x$求导，因此结果仍然是一个列向量。

上述公式的简单记忆方法： 前导不变，后导转置，公式表达为：
$$\frac{\partial x^{T}a}{\partial x}=a,\frac{\partial bx}{\partial x}=b^T$$
注意： 这里的$b$是行向量。

令$Z(\theta )=(X\theta -Y{)}^T(X\theta -Y)$，
则$J(\theta )=\frac{1}{m}Z(\theta )$，且 $$\begin{array}{rl}Z(\theta )& =({\theta }^T{X}^T-Y^T)(X\theta -Y)\\ & ={\theta }^T{X}^{T}X\theta -{\theta }^T{X}^{T}Y-Y^{T}X\theta +Y^{T}Y\end{array}$$

因此我们有： $$\begin{array}{rl}\frac{\partial (Z(\theta ))}{\partial \theta }& =(X^{T}X+X^{T}X)\theta -X^{T}Y-X^{T}Y+0\\ & =2X^{T}X\theta -2X^{T}Y\end{array}$$

代入$Z(\theta )$到公式$(2)$，我们有： $$\begin{array}{rl}\frac{\partial (\frac{1}{m}Z(\theta ))}{\partial \theta }& =0\\ \text{即：}\frac{1}{m}(2X^{T}X\theta -2X^{T}Y)& =0\\ \text{即：}{X}^{T}X\theta -X^{T}Y& =0\end{array}$$

如果$X^{T}X$可逆，那么我们有： $$\theta =(X^{T}X{)}^{-1}{X}^{T}Y$$

上述的优化方法需要计算矩阵的逆，当数据量很小而且可逆的时候，速度快，效果好，但是并不适用于大量高维度的数据。是否有一些数学中的迭代的方式能不断的逼近最小值呢？这个答案有很多，下面将尽力逐个介绍。

梯度下降

梯度下降，顾名思义就是顺着梯度的方向的反方向下滑，直到滑动到某一个最低点为止。
我们可以把它想象成一个小山，如下图所示（
很抱歉是用latex写的，画图采用的tikz，第一次用不熟练，很丑），假定我们的起始点为start，然后顺着下山的方向，一步一步的就会滑落到它左边的最低点。

这可能会有一个问题： 如果它的右边有一个更低的最低点呢？
它有可能落不到全局最低，但是可以达到局部最低。不过，如果我们的函数是凸函数，也就是说只有一个极值点的时候，它的局部最优也就是全局最优了。

梯度下降的``小山''图，从start点开始下落

原理篇

一维直观篇

为了与高等数学教材的保持一致，采用书里面的记法，在$△x\to 0$时，我们有：
$$f(x+△x)=f(x)+△x\cdot f^{\prime }(x)+\text{o}(△x)$$
此后我们简单的忽略$\text{o}(△x)$这一项。

令上图中的start点的横坐标为$x_0$，代入$x=x_0$，我们有:
$$f(x_0+△x)=f(x_0)+△x\cdot f^{\prime }(x_0)$$

为了使它向着小山下方（左面）滑动，我们需要有$f(x_0+△x)<f(x_0)$，即
$$△x\cdot f^{\prime }(x_0)<0$$

那么$△x$需要满足什么条件呢？我们令$g=f^{\prime }(x),\sigma =△x$，我们有：
$$△x\cdot f^{\prime }(x_0)=g\cdot \sigma$$
很显然，我们只需要$\sigma$取得$g$的相反的方向，例如$\sigma =-0.01g$，即可满足$g\cdot \sigma <0$。请原谅我，虽然进行$g$和$\sigma$的定义看起来多此一举，但请相信我，当将其映射到高维的时候，它能更好的帮助理解。

负梯度：这里面的$g$是斜率，其实也是梯度。$\sigma$与$g$变化方向相反，因此我们称之为$x$沿着梯度下降的方向，也就是负梯度方向。

学习率：
上面我们说$\sigma$只需要与梯度的方向相反，那么$\sigma =-0.1g,\sigma =-0.01g$都可以满足要求，因此我们使用变量$\alpha$表达学习率的概念，即$\sigma =-\alpha \cdot g$。

你应该已经知道学习率为何不能太大的原因了吧，上面的描述都是基于$\sigma =△x\to 0$的情况下，进行的推导，在实际中我们通常不会选用太大的$\alpha$，不过太小的$\alpha$意味着学习速度变慢，也就是说需要迭代更多的步数才可以到达最小值。

假定学习率$\alpha =0.01$，从而我们有迭代公式：
$$x_1=x_0-0.01g,\dots ,x_{n+1}=x_n-0.01g$$

利用如上的递推公式，我们就可以不断的使得$f(x_{n+1})<f(x_n)$，即不断下滑到最低点。

二维梯度篇

我先从二维逐步扩展到高维的形式，二维的微分形式:
$$f(x+△x,y+△y)=f(x,y)+△x{f}^{\prime }(x)+△y{f}^{\prime }(y)+\text{o}(\sqrt{{△}^{2}x+{△}^{2}y})$$

我们的目标是使得在对于自变量$x,y$进行相应的$△x,△y$变化后，新的$f(x+△x,y+△y)$能够更小，也就是说使得：
$$△x{f}^{\prime }(x)+△y{f}^{\prime }(y)+\text{o}(\sqrt{{△}^{2}x+{△}^{2}y})<0$$

那$△x,△y$需要满足什么条件呢？
我们令$g^T=(f^{\prime }(x),f^{\prime }(y)),{\sigma }^T=(△x,△y)$，我们有：
$$△x{f}^{\prime }(x)+△y{f}^{\prime }(y)=g^T\cdot \sigma =||g||\cdot ||\sigma ||\cdot cos(\theta )$$
很显然，我们应该使得$cos(\theta )=-1$，也就是说$\sigma$的方向与梯度$g$的方向相反。

如果加上学习率，我们可以令$\sigma =-\alpha \cdot g$.
同样，与一维的情况一样，也是沿着负梯度的方向。

高维梯度篇

我们终于来到了高维情况，直接利用前面叙述的参数，梯度$g$，学习率$\alpha$，自变量变化$\sigma$，注意$\sigma$是一个向量，它代表各个自变量参数的变化情况，可简单把它想象成 { △ x , △ y , ⋯   } \{\triangle x, \triangle y, \cdots\} {△x,△y,⋯}。
我们有：
$$f(X+\sigma )=f(X)+g^T\cdot \sigma +\text{o}(||\sigma ||)$$

同样，我们应该有$\sigma =-\alpha \cdot g$。

梯度下降简单实现demo篇

看到这里，有没有想实现一个梯度下降的冲动，反正我是有的，因此，我用python写了一个非常非常基本的二维的梯度下降，来实现线性回归。

请稍微回到我们在section 1.1中介绍的平面图上的m个点$(x_i,y_i)$，且呈现线性关系，我们的平均误差方程为：
$$J(w,b)=\frac{1}{m}\sum _{i=0}^{m-1}{J}_i$$
其中$J_i(w,b)=(w\ast x_i+b-y_i{)}^2$。我们的目标是解出$w,b$以使得$J$最小。

梯度下降的4个要素：

1. 初始值 $$w_0=0,b_0=0$$

2. 学习率 $$\alpha =0.01\to 0.0001$$

3. 梯度 $$\begin{array}{rl}\frac{\partial J(w,b)}{\partial w}& =\frac{1}{m}\sum _{i=0}^{m-1}2(w{x}_i+b-y_i)x_i\\ \frac{\partial J(w,b)}{\partial b}& =\frac{1}{m}\sum _{i=0}^{m-1}2(w{x}_i+b-y_i)\end{array}$$

4. 梯度更新 $$\begin{array}{rl}{w}_{n+1}& =w_n-\alpha \cdot \frac{\partial J(w,b)}{\partial w}\\ b_{n+1}& =b_n-\alpha \cdot \frac{\partial J(w,b)}{\partial b}\end{array}$$

我们使用的数据点如下面的散点图所示，它是使用$w=2,b=3$的公式在 x ∈ { 0 , 1 , … , 99 } x \in \{0, 1, \ldots, 99 \} x∈{0,1,…,99}中加入一定的随机数得到的。

线性回归之梯度下降数据分布图

生出散点图的python代码如下

import numpy as np
import random
import matplotlib.pyplot as plt
import pdb

random.seed(10)
x = np.arange(100)
w, b = 2, 3
y = np.array([w * i + b + random.randint(-5,5) for i in x])
one = np.ones(len(x))
plt.plot(x, y, '.')
plt.show()

利用梯度下降的4个要素编写的代码如下：

#alpha and g
alpha = 0.0003
w0, b0 = 0, 0

def gradient_w(w, b):
    return np.average([2*(w*xi + b - yi)* xi for (xi, yi) in list(zip(x,y))])
def gradient_b(w, b):
    return np.average([2*(w*xi + b - yi) for (xi, yi) in list(zip(x,y))])
    
for i in range(100000):
    w1 = w0 - alpha * gradient_w(w0, b0)
    b1 = b0 - alpha * gradient_b(w0, b0)
    w0, b0 = w1, b1
plt.plot(x, y, 'k+')
plt.plot(x, w0 * x + b0, 'r')
plt.show()

最终迭代出的结果为： $w1=1.99700571226b1=3.19821704612$

画出的曲线如下图所示

线性回归之梯度下降曲线结果

代码其实几分钟就写完了，但是一直都拟合不出来，当时我设置的学习率为$alpha=0.001$，还没有迭代多少步，就出现了值无穷大(nan)的情况，我一直以为是代码有问题，检查了一遍又一遍，还是感觉代码没问题，后来直到调整了学习率之后，才能够正确拟合出结果。

学习率的影响：
学习率过小，导致学习速度过慢，如果设置了最大迭代次数，那么很有可能还没有学到最终结果的时候，就已经终止了。不过我们可以看到它的误差，即$J(w,b)$是呈现一个不断下降的趋势。
学习率过大，虽然学习速度上去了，但是很有可能跳出了最优解，这可能会导致算法最终并不会收敛，误差通常会溢出。

我从网上盗取了两张图，来分别展示学习率大小对结果的影响。

线性回归之梯度下降不同学习率小（左）和大（右）的情况

对于学习率引起的问题，我给上述的梯度下降代码加入了误差计算并存储，并画出误差图。

线性回归之梯度下降不同学习率小（左）和大（右）的误差情况

在使用学习率$\alpha =0.00001$，迭代2000次后，我们得到$w1=2.04427897398,b1=0.0626663170812$，可以看出这与实际的$w=2,b=3$还是有不小的差距的，其误差画出曲线如误差图中的左图所示。
在使用学习率$\alpha =0.001$，迭代2000次后，我们得到$w1=nan,b1=nan,error=nan$，它的误差如误差图的右图所示。

根据误差图，就可以很容易指导我们是学习率不够，迭代次数不够，还是学习率过大。当误差依然处于下降的趋势，我们的迭代次数通常是不够的，如果时间不允许，那么可以稍微调整学习率，使用更大的学习率来加快收敛，但是不能过大，以免出现误差不收敛。

牛顿法

牛顿法与梯度下降法具有很大的相似性，区别在于梯度下降是采用的一阶导数，而牛顿法采用的二阶导数。

一维直观篇

请拿上我们的高等数学中的泰勒展开，采用书里面的记法，在$△x\to 0$时，我们有：
$$f(x+△x)=f(x)+△x\cdot f^{\prime }(x)+\frac{1}{2}{△}^{2}x\cdot f^{\prime \prime }(x)+\text{o}({△}^{2}x)$$

同样的，我们需要求得这样的$△x$，以使得$f(x+△x)$尽可能的小。回顾在梯度下降中，我们只是将泰勒展开到前面的两项，然后得到$△x=-\alpha \cdot f^{\prime }(x)$可以使得$f(x+△x)$呈现下降趋势，即沿着负梯度的方向。

在牛顿法中，我们的泰勒展开到了$△x$的平方项，这使得我们可以换个思路以最小化$f(x+△x)$：将$△x$当成未知量，公式$f(x)+△x\cdot f^{\prime }(x)+\frac{1}{2}{△}^{2}x\cdot f^{\prime \prime }(x)$是一个一元二次方程，曲线的开口向上，因此我们有$△x=-\frac{f^{\prime }(x)}{f^{\prime \prime }(x)}$使得$f(x+△x)$取得最小值。
它与梯度下降很大的不同在于梯度下降的$△x$只要求方向与梯度方向相反即可，而牛顿法却可以精确的求出$△x$的值。

高维原理篇

假定有n个自变量参数 { x 1 , x 2 , … , x n } \{x_1, x_2, \ldots, x_n\} {x1​,x2​,…,xn​}，，自变量变化$\sigma$，注意$\sigma$是一个向量，它代表各个自变量参数的变化情况，可简单把它想象成 { △ x 1 , △ x 2 , ⋯   } \{\triangle x_1, \triangle x_2, \cdots\} {△x1​,△x2​,⋯}，这里我们还需要额外的一个$G$，它表示对自变量的二阶导数。

我们有学习率$\alpha$(这是一个数)，梯度向量$g$: $$g=\left[\begin{array}{c}\frac{\partial f}{\partial x_1}\\ \frac{\partial f}{\partial x_2}\\ ⋮\\ \frac{\partial f}{\partial x_n}\end{array}\right]$$ 自变量的变化$\sigma$: $$\sigma =\left[\begin{array}{c}△x_1\\ △x_2\\ ⋮\\ △x_n\end{array}\right]$$ 二阶导数矩阵hessian矩阵$G$： $$G=\left[\begin{array}{c}\frac{{\partial }^{2}f}{\partial x_1^2},\frac{{\partial }^{2}f}{\partial x_1\partial x_2},\dots ,\frac{{\partial }^{2}f}{\partial x_1\partial x_n}\\ \frac{{\partial }^{2}f}{\partial x_2\partial x_1},\frac{{\partial }^{2}f}{\partial x_2^2},\dots ,\frac{{\partial }^{2}f}{\partial x_2\partial x_n}\\ ⋮\\ \frac{{\partial }^{2}f}{\partial x_n\partial x_1},\frac{{\partial }^{2}f}{\partial x_n\partial x_2},\dots ,\frac{{\partial }^{2}f}{\partial x_n^2}\end{array}\right]$$ 可以看出$G$是对称矩阵，即$G=G^T$。

我们有：
$$f(X+\sigma )=f(X)+g^T\cdot \sigma +\frac{1}{2}{\sigma }^T\cdot G\cdot \sigma +\text{o}(||\sigma |{|}^2)$$

利用前面学到的向量求导公式，对$\sigma$求导，使得导数为0，即
$$\begin{array}{rl}\frac{\partial f(X+\sigma )}{\partial \sigma }& =g+\frac{1}{2}(G+G^T)\cdot \sigma \\ & =g+G\cdot \sigma \\ & =0\end{array}$$ 因此，我们有: $$\sigma =-G^{-1}\cdot g$$
然后我们就可以利用得到的$\sigma$，得到新的$X_{new}=X+\sigma$，进行下一步迭代。

牛顿法代码Demo篇

计算$g$和$G$的代码：

def gradient_w(w, b):
    return np.average([2*(w*xi + b - yi)* xi for (xi, yi) in list(zip(x,y))])
def gradient_b(w, b):
    return np.average([2*(w*xi + b - yi) for (xi, yi) in list(zip(x,y))])

def gradient_w_w(w, b):
    return np.average([2 * xi * xi for (xi, yi) in list(zip(x, y))])
def gradient_w_b(w, b):
    return np.average([2 * xi for (xi, yi) in list(zip(x, y))])
def gradient_b_w(w, b):
    return np.average([2 * xi for (xi, yi) in list(zip(x,y))])
def gradient_b_b(w, b):
    return np.average([2 for (xi, yi) in list(zip(x,y))])
def error(w, b):
    return np.average([np.square(w*xi + b - yi) for (xi, yi) in list(zip(x,y))])

牛顿迭代的代码：

erros = [[], []]
for i in range(2000):
    g = np.mat([gradient_w(w0, b0), gradient_b(w0, b0)]).T
    G = np.mat([[gradient_w_w(w0, b0), gradient_w_b(w0, b0)], [gradient_b_w(w0, b0), gradient_b_b(w0, b0)]])
    
    sigma = -G.I * g
    w1 = w0 + sigma[0, 0]
    b1 = b0 + sigma[1, 0]
    
    erros[0].append(i)  
    erros[1].append(error(w1, b1))
    w0, b0 = w1, b1
    
plt.plot(x, y, 'k+')
plt.plot(x, w0 * x + b0, 'r')
plt.show()

从误差结果看，它一次迭代就可以得到最小值，后经过@景宽提醒，牛顿法本身计算的就是二阶泰勒展开的值。我们上述的直线方程只有两个系数，它的三阶泰勒值就是0，因此二阶泰勒展开就是它的最小值了。

最小二乘法的最大似然解释

前面在Section基本介绍中，我们直接用平方损失函数来最小化来得到最优直线。那么，你是否有疑问，为什么平方损失函数最小就可以得到最优的直线，而不是绝对值损失，或者四次方损失呢？

这其实涉及到概率论中的最大似然估计.
上述的问题，可以转化为：对于$m$个点形成的集合 D = { d 1 , d 2 , … , d m } D = \{d1, d2, \ldots, d_m\} D={d1,d2,…,dm​}，我们需要拟合一条直线$h$，以使得拟合的直线最切合这个数据集合$D$，也就是说我们要最大化概率$P(h|D)$。

由贝叶斯公式 $$\begin{array}{rl}P(h|D)& =\frac{P(D|h)\cdot P(h)}{P(D)}\propto P(D|h)\end{array}$$
如果各个点相互独立，则
$$P(h|D)\propto P(d_1|h)\cdot P(d_2|h)\cdot \dots \cdot P(d_n|h)$$

对于$P(d_i|h)$表示的是在给定直线$h$的情况下，具有点$d_i$的概率。我们假设各个点的误差$△y_i=\overline{y_i}-y_i$服从均值为0，方差为$\sigma$的正态分布，也就是说，在给定直线$h$的情况下，点$d_i$出现的概率为服从：
$$P(d_i|h)\propto e^{-\frac{{△}^2{y}_i}{2{\sigma }^2}}$$

于是，我们有:
$$P(h|D)\propto e^{{\sum }_{i=1}^m(-\frac{{△}^2{y}_i}{2{\sigma }^2})}$$

从而最大化$P(h|D)$，等价于最小化${\sum }_{i=1}^m({△}^2{y}_i)$，也就是前面说的最小化平方误差。

另一种说法：
对于样本点$d_i$，对样本点的预测取决于参数$\theta$的选取，且我们有
$$△y_i=\overline{y_i}-y_i$$
其中$\overline{y_i}$是$y_i$的预测值，$y_i$为真实值，$△y_i$为误差。

一般我们认为误差服从正态分布，即 $$\begin{array}{r}△y_i\sim \text{N}(0,\sigma )\end{array}$$
现在我们需要选取$\theta$，以使得取得$y_i$的概率最大，即
$$L(\theta )=\prod _{i=1}^{m}P(y_i|\theta )=\prod _{i=1}^m\frac{1}{\sqrt{2\pi }\sigma }\exp (-\frac{{△}^2{y}_i}{2\sigma })$$
两边取$log$，乘积变求和，因此同样转换为最小化${\sum }_{i=1}^m({△}^2{y}_i)$，也就是前面说的最小化平方误差。

参考文献：

https://www.cnblogs.com/21207-iHome/p/5222993.html

http://blog.youkuaiyun.com/ying_xu/article/details/51240291

https://www.cnblogs.com/happylion/p/4172632.html

http://blog.youkuaiyun.com/lydyangliu/article/details/9208635

http://www.fuzihao.org/blog/2014/06/13/为什么最小二乘法对误差的估计要用平方/

https://www.zhihu.com/question/20447622