（一）线性回归算法

---

前言

我正在学习吴恩达的机器学习课程，想在学习的过程中写一写博客。一是相当于做了笔记，总结一下机器学习的知识；二是在总结的过程中也会发现很多细节的问题可以让自己收获更多；三是分享到网上，或许对别人有所帮助，也或许我的错误之处被别人发现，也可以促使我改正。

---

一、线性回归模型

提供一定的数据作为训练集，训练集中有m个训练数据，假设训练集中每个样本为$\left(x,y\right)$，其中$x=\left(x_1,x_2,\cdots ,x_n\right)$为特征向量，$y$为输出变量。我们将第$i$个样本记为$\left(x^{(i)},y^{(i)}\right)$，其中$x^{\left(i\right)}=\left(x_1^{\left(i\right)},x_2^{\left(i\right)},\cdots ,x_n^{\left(i\right)}\right)$。我们记：$$\begin{gathered} X=\left(\begin{array}{cccc}1& x_1^{\left(1\right)}& \cdots & x_n^{\left(1\right)}\\ 1& x_1^{\left(2\right)}& \cdots & x_n^{\left(2\right)}\\ ⋮& ⋮& & ⋮\\ 1& x_1^{\left(m\right)}& \cdots & x_n^{\left(m\right)}\end{array}\right)=\left(\begin{array}{cccc}{x}_0^{\left(1\right)}& x_1^{\left(1\right)}& \cdots & x_n^{\left(1\right)}\\ x_0^{\left(2\right)}& x_1^{\left(2\right)}& \cdots & x_n^{\left(2\right)}\\ ⋮& ⋮& & ⋮\\ x_0^{\left(m\right)}& x_1^{\left(m\right)}& \cdots & x_n^{\left(m\right)}\end{array}\right) \\ y=\left(\begin{array}{c}{y}^{\left(1\right)}\\ y^{\left(2\right)}\\ ⋮\\ y^{\left(m\right)}\end{array}\right) \end{gathered}$$
（注意上文中两个$y$的含义，虽然符号相同，但是含义不同）

为了下面的公式表示简单，我们重新定义$x=\left(1,x_1,x_2,\cdots ,x_m\right)$。

我们用这些数据来训练一个模型，用训练好的模型去预测新的数据。监督学习就是从训练集中“学习”，或称为“训练”，与人的学习类似，通过不断地学习，就可以知道这个特定问题的一些规律，然后就可以根据给定的特征向量$x$预测$y$。

假设特征向量和预测变量大致呈现出一种线性关系，即：$$\exists \theta ={\left({\theta }_0,{\theta }_1,\cdots ,{\theta }_n\right)}^T\in {\mathbf{R}}^{\mathbf{n}+1}s.t.{\theta }^{T}x\approx y$$
这里可以令$h_{\theta }\left(x\right)={\theta }^{T}x$，方便表示，同时也有利于和逻辑回归算法对比。

实际中很少有绝对的线性关系，所以上式使用约等于。但是当误差很小时，也能满足很多需求，我们的线性回归算法在学习时就是要想办法缩小以下误差：$$\sum _{i=1}^m{\left(h_{\theta }\left(x^{\left(i\right)}\right)-y^{\left(i\right)}\right)}^2$$
为了导函数比较简单，我们常常使用下式：$$J\left(\theta \right)=\frac{1}{2}\sum _{i=1}^m{\left(h_{\theta }\left(x^{\left(i\right)}\right)-y^{\left(i\right)}\right)}^2$$
我们把上式$J(\theta )$称为损失函数，使用矩阵的方式来表示就是：$$J\left(\theta \right)=\frac{1}{2}{\left(X\theta -y\right)}^T\left(X\theta -y\right)=\frac{1}{2}{\Vert X\theta -y\Vert }_2^2$$
我们的目标就是最小化这个损失函数，寻找使得损失函数最小的$\theta$。

---

二、求解方法

1.正规方程法

损失函数是一个凸函数，直接对函数关于$\theta$求导可得：$$\frac{\partial J\left(\theta \right)}{\partial \theta }=X^T\left(X\theta -y\right)$$
令上述导数等于向量$0$（$n+1$维列向量）：$$\frac{\partial J\left(\theta \right)}{\partial \theta }=X^T\left(X\theta -y\right)=X^{T}X\theta -X^{T}y=0$$得出：$$\theta ={\left(X^{T}X\right)}^{-1}{X}^{T}y$$
这样就可以得出全局最优解了，这叫做正规方程解法。或许在使用这个公式时会遇到$X^{T}X$不可逆的情况，这个时候通常是选择的特征中有线性相关性比较强的特征，这时可以尝试去寻找线性相关性强的特征，删掉一些特征大概率会解决这个问题。当然也可以使用$X^{T}X$的伪逆来代替的逆，这样得到的答案也是可靠的。

这里顺便再计算一下偏导数：$$\frac{\partial J\left(\theta \right)}{\partial {\theta }_k}=\sum _{i=1}^m\left(h_{\theta }\left(x^{\left(i\right)}\right)-y^{\left(i\right)}\right)x_k^{\left(i\right)}$$

2.梯度下降算法

对应于正规方程的解法，还可以使用迭代解法去求解全局最优解。使用梯度下降法也可以寻找全局最优解。梯度的方向是函数在该点上升最快的方向，只要自变量向梯度的反方向移动，函数就会下降。该问题的损失函数$J\left(\theta \right)$的梯度如下：$${\nabla }_{\theta }J\left(\theta \right)=\frac{\partial J\left(\theta \right)}{\partial \theta }=X^T\left(X\theta -y\right)$$
该问题的梯度下降法步骤如下：

1. 初始化参数$\theta$ ;
2. 不断执行下列操作：$$\begin{gathered} \theta ⟵\theta -\alpha X^T\left(X\theta -y\right)\text{或} \\ {\theta }_k⟵{\theta }_k-\alpha \sum _{i=1}^m\left(h_{\theta }\left(x^{\left(i\right)}\right)-y^{\left(i\right)}\right)x_k^{\left(i\right)}\left(k=0,1,\cdots ,n\right) \end{gathered}$$直到$\theta$的变化很小时停止。

这就是梯度下降法的大致步骤，步骤中$\alpha$（$\alpha >0$）称为学习率，控制$\theta$沿着梯度的反方向移动的步长。学习率的取值要恰当，当学习率$\alpha$太大，移动的步长就太大，最终可能导致$\theta$在全局最优解附近反复横跳，甚至远离全局最优解；当学习率$\alpha$太小，迭代次数通常也会增大，可能会导致收敛速度太慢。

梯度下降法看起来比正规方程法要复杂一些，因为正规方程只需要计算一下表达式就可以直接求出最优的$\theta$，用代码实现正规方程算法也只需要一条计算的语句就可以求出最优的$\theta$，而梯度下降算法需要考虑到学习率的问题，还要不断迭代，用代码实现起来也比正规方程稍微复杂。不过有时候梯度下降算法却有着优势。当样本的特征数目很大时，正规方程的求逆运算就会耗费较多的时间，这时候可能梯度下降训练速度比正规方程更快。通常特征数目大于10000时，正规方程的运算代价就会比较高，这时选择梯度下降较好。

上面的梯度下降算法在每次迭代时使用训练集中所有的数据进行计算，我们把它称为批量梯度下降法，但是当训练集样本数很多时，批量梯度下降法耗费的时间也会很多。这时我们可以对训练样本单独做梯度下降，这样会减少训练时间，我们把这样的方法称为随机梯度下降法。

一个样本$\left(x,y\right)$的损失函数为：$$j\left(\theta \right)=\frac{1}{2}{\left(h_{\theta }\left(x\right)-y\right)}^2$$梯度为：$${\nabla }_{\theta }j\left(\theta \right)=\left(h_{\theta }\left(x\right)-y\right)x$$
因此该问题的随机梯度下降法步骤如下：

1. 将训练集中的样本顺序打乱；
2. 初始化参数$\theta$ ;
3. 对于$i=1,2,\cdots ,m$ ，依次进行如下计算：$$\theta ⟵\theta -\alpha \left(h_{\theta }\left(x^{\left(i\right)}\right)-y^{\left(i\right)}\right)x^{\left(i\right)}$$
4. 重复上一步骤几次（通常 $1\sim 10$ 次）。

我之前没有想通随机梯度下降法和批量梯度下降法的区别，我现在想通了，是一个微小的问题被我察觉了。
我觉得下面这样的算法也是可以提高效率的，下面的算法和上面的算法主要区别在于：上面的算法一定可以用到所有样本点的信息，而下面的算法并不一定会用到所有样本点的信息：

• 初始化参数 $\theta$ ；

• 随机选择一个样本点 $\left(x^{\left(k\right)},y^{\left(k\right)}\right)$ ，执行下列操作：
  $$\theta ⟵\theta -\alpha \left(h_{\theta }\left(x^{\left(k\right)}\right)-y^{\left(k\right)}\right)x^{\left(k\right)}$$

• 重复上一步骤多次。

批量梯度下降法一定可以达到全局最优解，但是随机梯度下降法通常不能找到全局最优解，不过通常与全局最优解也比较接近，可能也会取得不错的效果。

---

三、概率解释

将线性回归算法的训练过程转化为了优化目标函数的过程，这样做还可以运用概率统计方面的知识去解释。

我们定义符号：$\left(x,y\right)$是一个随机变量，而训练集中的每个样本$\left(x^{\left(i\right)},y^{\left(i\right)}\right)$是随机变量$\left(x,y\right)$的取值。

我们的目的是寻找参数使得对于每组的样本$\left(x^{\left(i\right)},y^{\left(i\right)}\right)$，都有${\theta }^T{x}^{\left(i\right)}$与$y^{\left(i\right)}$的差别不大，模型中我们使用了所有样本平方误差的总和来衡量差别不大，为什么不用绝对误差总和来衡量呢？为什么不用别的方式来衡量呢？下面就是答案。

对于每组样本，存在一个误差项${\epsilon }^{\left(i\right)}$，使得：
$$y^{\left(i\right)}={\theta }^T{x}^{\left(i\right)}+{\epsilon }^{\left(i\right)}$$
根据中心极限定理，可以假设${\epsilon }^{\left(i\right)}$是一个随机变量，并且${\epsilon }^{\left(i\right)}\sim N\left(0,{\sigma }^2\right)$，那么对于给定的$x^{\left(i\right)}$，$y\sim N\left({\theta }^T{x}^{(i)},{\sigma }^2\right)$，那么对于每组样本$\left(x^{\left(i\right)},y^{\left(i\right)}\right)$，$y$的概率密度函数取值为：$$p\left(y=y^{\left(i\right)}|x=x^{\left(i\right)};\theta \right)=\frac{1}{\sqrt{2\pi }\sigma }\exp \left(-\frac{{\left(y^{\left(i\right)}-{\theta }^T{x}^{\left(i\right)}\right)}^2}{2{\sigma }^2}\right)$$
这样就可以运用极大似然估计了，似然函数：$$\begin{gathered} L\left(\theta ,x^{\left(1\right)},\cdots ,x^{\left(m\right)},y^{\left(1\right)},\cdots ,y^{\left(m\right)}\right)=\prod _{i=1}^{m}p\left(y=y^{\left(i\right)}|x=x^{\left(i\right)};\theta \right) \\ =\prod _{i=1}^m\frac{1}{\sqrt{2\pi }\sigma }\exp \left(-\frac{{\left(y^{\left(i\right)}-{\theta }^T{x}^{\left(i\right)}\right)}^2}{2{\sigma }^2}\right) \end{gathered}$$
对数似然函数为：$$\ln L\left(\theta \right)=m\ln \frac{1}{\sqrt{2\pi }\sigma }-\frac{1}{2{\sigma }^2}\sum _{i=1}^m{\left(y^{\left(i\right)}-{\theta }^T{x}^{\left(i\right)}\right)}^2$$
最大化该对数似然函数，就等价于最小化下面的函数：$$J\left(\theta \right)=\frac{1}{2}\sum _{i=1}^m{\left({\theta }^T{x}^{\left(i\right)}-y^{\left(i\right)}\right)}^2$$
这就是该问题的概率解释。

---

四、python代码

仿照sklearn库，写出线性回归类。

import numpy as np

class LR(object):
    theta=0
    def __Normal_Equation(self, X, y): #正规方程解法
        self.theta=np.dot(np.dot(np.linalg.pinv(np.dot(X.T,X)),X.T),y)
    def __Gradient_Descent(self, X, y): #梯度下降解法（批量）
        alpha = 0.01  #学习率
        iternum = 0  #记录迭代次数
        while True:
            tmp = self.theta-alpha*np.dot(X.T,np.dot(X,self.theta)-y)
            if np.allclose(self.theta,tmp):
                break
            self.theta = tmp
            iternum += 1
        print("迭代次数：",iternum)
    def __SGD(self, X ,y):
        alpha = 0.1  #学习率
        idxs = np.arange(X.shape[0])
        np.random.shuffle(idxs)
        X, y = X[idxs], y[idxs]  # 打乱数据
        for i in range(10): # 外层循环迭代10次看看效果
            for k in range(X.shape[0]):
                self.theta = self.theta - alpha*(np.dot(self.theta,X[k,:])-y[k])*X[k,:]

            
    def fit(self, X_train, y_train, method = None):
        X, y=np.insert(X_train, 0, 1, axis=1), y_train  #这样写起来比较方便
        self.theta = np.zeros(X.shape[1]) #初始化
        if method == 'GD': #梯度下降解法
            self.__Gradient_Descent(X, y)
        elif method == 'NE':  #正规方程解法
            self.__Normal_Equation(X, y)
        elif method == 'SGD':  #随机梯度下降法
            self.__SGD(X, y)
        elif y.shape[0] <= 10000:  #这个时候用正规方程
            self.__Normal_Equation(X, y)
        else:                      #数据量大，用梯度下降
            self.__Gradient_Descent(X, y)
            
    def predict(self, X):  #预测
        return np.dot(np.insert(X, 0, 1, axis=1), self.theta)

自己创造数据利用上述类来画出回归曲线。

import matplotlib.pyplot as plt

def Create():
    return np.array([[1,2],[2,3.4],[1.5,3.2],[3.4,8],[2.5,5],[2.3,5.2]])

data=Create()  #导入数据
plt.scatter(data[:,0],data[:,1])  #可视化数据
X=data[:,0:data.shape[1]-1]
y=data[:,-1]  #数据中最后一列
np.random.seed(0)

reg = LR()
#正规方程解法
print("正规方程法：")
reg.fit(X, y)  #数据量小，默认正规方程法
theta = reg.theta
print("参数theta：",theta)
x=np.linspace(np.min(data[:,0])-1,np.max(data[:,0])+1,2)
plt.plot(x,theta[1]*x+theta[0],color='green')

#梯度下降解法（批量梯度下降法）
print("批量梯度下降法：")
reg.fit(X, y, method = 'GD')  #梯度下降法
theta = reg.theta
print("参数theta：",theta)
plt.plot(x,theta[1]*x+theta[0],color='red')

#随机梯度下降法
print("随机梯度下降法：")
reg.fit(X, y, method = 'SGD')
theta = reg.theta
print("参数theta：",theta)
plt.plot(x,theta[1]*x+theta[0],color='yellow')

plt.legend(['Normal Equation','Gradient Descent','SGD','data'])
plt.show()

结果如下：

正规方程法：
参数theta： [-0.74108602 2.4603556 ]
批量梯度下降法：
迭代次数： 1214
参数theta： [-0.73981544 2.45981353]
随机梯度下降法：
参数theta： [-0.25680749 2.20781238]

上图看不见正规方程的回归曲线是因为梯度下降的回归曲线将其挡住了。可以看到随机梯度下降法会较快得接近全局最优解，不过它往往不能收敛到全局最优。

总结

线性回归算法就是利用训练集中的样本学习出合适的参数$\theta$，这样对于新的特征向量，我们可以输出${\theta }^{T}x$作为预测值。对于如何求解参数$\theta$，我们介绍了两种方法，一种是直接求解的正规方程法，一种是迭代求解的梯度下降法，根据实际情况灵活选择求解方法。