机器学习算法理论与代码实现——线性回归

一 线性模型

给定由n个属性描述的列向量$\textbf{x}={(x^{(1)};x^{(2)};...;x^{(n)})}$，其中 $x^{(j)}$是$\textbf{x}$在第$j$个属性的取值。线性模型即为通过对属性进行线性组合的函数，即

$$f(\textbf{x})=w_0+w_1x^{(1)}+...+w_nx^{(n)}$$
写成向量形式如下：

$$f(\textbf{x})=\mathbf{w}^\mathrm{T}\mathbf{x}$$
其中 列向量 $\mathbf{w}=(w_0;w_1;...;w_n)$， 列向量 $\mathbf{x}=(1;x^{(1)};...;x^{(n)})$。 列向量 $\mathbf{w}$确定后，模型随之确定。
线性模型形式简单，易于建模；
直观表达了各属性在预测中的重要性，因此具有很好的可解释性；

二 线性回归模型求解

对于给定的数据集$\mathbf{D}=\left \{ (\mathbf{x_1},y_1),(\mathbf{x_2},y_2),...,(\mathbf{x_m},y_m)\right \}$,其中$\mathbf{x_i}=(x_i^{(1)};...;x_i^{(n)})$，$y_i$为第$i$个实例的实际值。“线性回归”试图学得一个线性模型以尽可能准确地预测实例的输出值，使之接近实际值。
关键问题是如何衡量两者之间的误差。这里采用均方误差作为性能度量，即利用最小二乘法来进行参数估计。

$$\underset{\mathbf{w}}{argmin}\mathbf{\mathit{J}}(\mathbf{w})=\frac{1}{2m}\sum_{i=1}^{m}(f(\mathbf{x}_i)-y_i)^2$$
其实这里实质是用极大似然函数来进行参数估计，在一定假设前提下推导出上述目标函数进行参数估计,推导过程如下。
根据中心极限定理，认为误差项 $\mathbf\xi$服从均值为零的高斯分布

$$P(\xi_i)=\frac{1}{\sqrt{2\pi }\sigma }exp(-\frac{\xi_i^2}{2\sigma ^2})$$
$$P(y_i|\mathbf{x_i};\mathbf{w})=\frac{1}{\sqrt{2\pi }\sigma }exp(-\frac{(y_i-\mathbf{x}_i)^2}{2\sigma ^2})$$
由上可得似然函数，
$$\mathbf{L}(\mathbf{w})=\prod_{i=1}^{n}P(y_i|\mathbf{x_i};\mathbf{w})$$
取对数，得对数似然函数
$$\mathbf{L}(\mathbf{w})=nlog\frac{1}{\sqrt{2\pi }\sigma }-\frac{1}{2\sigma ^2}\sum_{i=1}^{m}(f(\mathbf{x}_i)-y_i)^2$$
对上式取极大值等价于下式取极小值
$$\mathbf{\mathit{J}}(\mathbf{w})=\frac{1}{2m}\sum_{i=1}^{m}(f(\mathbf{x}_i)-y_i)^2$$
推导完毕。
模型求解方法：矩阵直接求解和梯度下降法。
1 矩阵求解法
对于数据集 $D$中的每个实例组成一个矩阵，矩阵形式如下： $$\mathbf{X}=\begin{pmatrix} 1 & x_1^{(1)} & ... & x_1^{(n)}\\ 1& x_2^{(1)} & ... & x_2^{(n)}\\ .& . & ... & .\\ 1& x_m^{(1)} & . & x_m^{(n)} \end{pmatrix}=\begin{pmatrix} 1 &\mathbf{x}_1^T \\ 1 &\mathbf{x}_2^T\\ .& .\\ 1& \mathbf{x}_m^T \end{pmatrix}$$
对应的实际值写成 列向量形式 $\mathbf{y}=(y_1;y_2;...;y_m)$,则有
$$\mathbf{\hat{w}}^*=\underset{\hat{\mathbf{w}}}{argmin}(\mathbf{y}-\mathbf{X\hat{w}})^T(\mathbf{y}-\mathbf{X\hat{w}})$$ 上式argmin后面部分对 $\hat{\mathbf{w}}$求导，另之等于零， $$\mathbf{\hat{w}}^*=(\mathbf{X^TX})^{-1}\mathbf{X^T}\mathbf{y}$$ 令 $\mathbf{\hat{x}}_i=(1;\mathbf{x}_i)$从而得到线性模型 $f(\mathbf{\hat{x}}_i)=\mathbf{\hat{x}}_i^T\mathbf{\hat{w}}^*$,或者 $f(\mathbf{\hat{x}}_i)=\mathbf{\hat{w}}^{*T}\mathbf{\hat{x}}_i$。
但是，现实情况中 $\mathbf{X^TX}往往不可逆$，通常原因有两种，一是高度共线性；二是数据特征过多而训练数据较少，此时可以通过正则化来解决。
2 梯度下降法
梯度下降是一种常用的一阶优化方法，是求解无约束优化问题的经典方法之一。对于连续可微函数上某一点，有各个方向导数，沿梯度方向的方向导数达到最大值，也就是说,梯度的方向是函数在这点增长最快的方向。
因此，我们可以得到如下结论：函数在某点的梯度是这样一个向量，它的方向与取得最大方向导数的方向一致，而它的模为方向导数的最大值。
所以我们可以沿反梯度方向不断一步一步迭代，得到局部极小点。当目标函数为凸函数时，局部极小点就是全局最小点，此时梯度下降法可确保收敛到全局最优解。
将损失函数对 列向量 $\mathbf{w}$求导，得到 $w_j$的偏导：
$$\frac{\partial \mathbf{J(w)}}{\partial w_j}=\frac{\partial }{\partial w_j}\frac{1}{2m}\sum_{i=1}^{m}(f(\mathbf{x}_i)-y_i)^2=\frac{1}{m}\sum_{i=1}^{m}(f(\mathbf{x}_i)-y_i)\mathbf{x}_i^{(j)},j=0,1,2,...,n$$ 然后对各个分量都以下面形式更新 $w_j$： $$w_j=w_j-\alpha\frac{1}{m} \sum_{i=1}^{m}(f(\mathbf{x}_i)-y_i)\mathbf{x}_i^{(j)}$$ 有公式可以看出对于每一个分量进行一次迭代时计算了所有训练样本数据，这种称为 批量梯度下降。因此在数据量很大的时候，每次迭代都要遍历训练集一遍，开销会很大。
为改善上述情况，可以在每次迭代仅选择一个训练样本去计算代价函数的梯度，然后更新参数。即使是大规模数据集，随机梯度下降法也会很快收敛。这种方法称为 随机梯度下降。此时有， $$w_j=w_j-\alpha (f(\mathbf{x}_i)-y_i)\mathbf{x}_i^{(j)}$$
比较：
批量梯度收敛速度慢，随机梯度收敛速度快。
批量梯度是在θ更新前对所有样例汇总误差，而随机梯度下降的权值是通过考查某个样本来更新的
批量梯度的开销大，随机梯度的开销小。

三 线性回归代码实现

import numpy as np
class MyLinearRegression():

    def __init__(self, n_iterations=10000, learning_rate=0.0001, regularization=None, gradient=True):
        '''初始化。是否正则化及L1L2的选择；选用梯度下降法还是正规方程法。梯度下降学习率以及迭代次数'''
        self.n_iterations = n_iterations
        self.learning_rate = learning_rate
        self.gradient = gradient
        if regularization == None:
            self.regularization = lambda x: 0
            self.regularization.grad = lambda x: 0
        else:
            self.regularization = regularization

    def initialize_weights(self, n_features):
        '''初始化权重.初始化模型参数,加入w0'''
        limit = np.sqrt(1 / n_features)
        w = np.random.uniform(-limit, limit, (n_features, 1))              #二维数组，n行一列。
        b = 0
        self.w = np.insert(w, 0, b, axis=0)                                

    def fit(self,X,y,):
        '''进行拟合'''
        m_samples, n_features = X.shape                                      # !!!
        self.initialize_weights(n_features)
        X = np.insert(X, 0, 1, axis=1)                                      #二维数组，每行前面加上元素1
        y = np.reshape(y, (m_samples, 1))                                    #二维数组，m 行一列
        self.training_errors = []
        if self.gradient == True:                                            #批量梯度下降
            for i in range(self.n_iterations):
                y_pred = X.dot(self.w)
                loss = np.mean(0.5 * (y_pred - y) ** 2)/m_samples + self.regularization(self.w)  # 矩阵运算
                '''mean()函数功能：求取均值
                经常操作的参数为axis，以m * n矩阵举例：
                axis 不设置值，对 m*n 个数求均值，返回一个实数
                axis = 0：压缩行，对各列求均值，返回 1* n 矩阵
                axis =1 ：压缩列，对各行求均值，返回 m *1 矩阵
                np.mean(X,axis=0或者1,keepdims=True)
                '''
                self.training_errors.append(loss)
                w_grad = X.T.dot(y_pred - y) + self.regularization.grad(self.w)  # (y_pred - y).T.dot(X)，计算梯度
                self.w = self.w - self.learning_rate * w_grad  # 更新权值w
        else:
            # 正规方程
            X = np.matrix(X)
            y = np.matrix(y)
            X_T_X = X.T.dot(X)
            X_T_X_I_X_T = X_T_X.I.dot(X.T)
            X_T_X_I_X_T_X_T_y = X_T_X_I_X_T.dot(y)
            self.w = X_T_X_I_X_T_X_T_y

    def predict(self, X):
        X = np.insert(X, 0, 1, axis=1)
        y_pred = X.dot(self.w)
        return y_pred
'''以二元为例，进行拟合'''
lr = MyLinearRegression()
X = np.array([[1,2],[2,4],[50,3],[23,59],[10,45],[10,61]])
y = np.array([3,6,53,82,55,71])
lr.fit(X,y)
y_test = lr.predict(np.array([[1,40],[2,6]]))
print(y_test)

这里用的是批量梯度下降法，其实也可以用随机梯度下降与小批量梯度下降，这里不再赘述。并且代码里省略了L1与L2正则化代码。
L1正则化下的损失函数

$$\mathbf{\mathit{J}}(\mathbf{w})=\frac{1}{2m}\sum_{i=1}^{m}(f(\mathbf{x}_i)-y_i)^2+\lambda \left \| W \right \|_{1}$$
其中 $\left \| W \right \|_{1}=\sum w_j$
L1正则化下的损失函数 $$\mathbf{\mathit{J}}(\mathbf{w})=\frac{1}{2m}\sum_{i=1}^{m}(f(\mathbf{x}_i)-y_i)^2+\lambda \left \| W \right \|_{2}^2$$
其中 $\left \| W \right \|_{2}=\sqrt{\sum w_j^2}$
L1正则化、L2正则化也称为Lasso正则化、Ridge正则化，其中 λλ 为模型的超参数。
L2 regularizer ：使得模型的解偏向于 norm 较小的 W，通过限制 W 的 norm 的大小实现了对模型空间的限制，从而在一定程度上避免了 overfitting 。不过 ridge regression 并不具有产生稀疏解的能力，得到的系数 仍然需要数据中的所有特征才能计算预测结果，从计算量上来说并没有得到改观。
L1 regularizer ： 它的优良性质是能产生稀疏性，导致 W 中许多项变成零。 稀疏的解除了计算量上的好处之外，更重要的是更具有“可解释性”。
关于正则化参考[ https://www.jianshu.com/p/a47c46153326]