线性回归

线性回归和Python代码

一、问题

1.1 线性模型基本形式

给定d个属性描述的示例 $\mathbf{x}=(x_1;x_2;...;x_d)$，其中$x_i$是$\mathbf{x}$在第$i$个属性上的取值，线性模型(linear model)试图学得一个通过属性的线性组合来进行预测的函数，即

$$f(\mathbf{x})={\omega }_1{x}_1+{\omega }_2{x}_2+...+{\omega }_d{x}_d+b\text{\hspace{1em}(1)}$$
一般向量形式写成

$$f(\mathbf{x})={\omega }^T\mathbf{x}+b\text{\hspace{1em}(2)}$$
其中$\omega =({\omega }_1;{\omega }_2;...;{\omega }_d)$

1.2 线性回归

1. 给定数据集$\mathbf{D}$ = {(${\mathbf{x}}_1$,$y_1$),(${\mathbf{x}}_2$,$y_2$),…,(${\mathbf{x}}_m$,$y_m$)}，其中 ${\mathbf{x}}_i$ = {$x_{i1}$;$x_{i2}$;…;$x_{id}$}，$y_i$ $\in$ $\mathbb{R}$
   线性回归问题(linear regression)试图学习模型
   $$f(\mathbf{x})={\omega }^T\mathbf{x}+b$$

对离散属性，若属性值间存在“序”(order)关系，可通过连续化将其转化为连续值，例如二值属性“身高”的取值“高” “矮”可转化为{1.0, 0.0}，三值属性高度的取值“高” “中” “低”可转化为{1.0, 0.5, 0.0}；若属性间不存在序关系，假定有k个属性值，则通常转化为k维向量，例如“瓜类”的取值“西瓜” “南瓜” “黄瓜” 可转化为(0, 0, 1),(0, 1, 0),(1, 0, 0)

2. 使得$f(x_i)\simeq y_i$

该问题也被称作多元线性回归(multivariate linear regression)

2. 对于每个$x_i$，其预测值$y_i={\omega }^T{\mathbf{x}}_i+b$,采用均方误差（平方损失函数），则在训练集$\mathbf{D}$上，模型的损失函数为:
   $$L(f)=\sum _{n=1}^N(f(x_i)-y_i{)}^2=\sum _{n=1}^N(y_i-{\omega }^T{\mathbf{x}}_i-b{)}^2\text{\hspace{1em}(3)}$$
   优化目标是损失函数最小化，即
   $$({\omega }^{\ast },b^{\ast })=\underset{(\omega ,b)}{\mathrm{argmin}}\sum _{n=1}^N(y_i-{\omega }^T{\mathbf{x}}_i-b{)}^2\text{\hspace{1em}(4)}$$

二、求解

• 可以用梯度下降法来求解上述最优化问题的数值解，但是实际上该最优化问题可以通过最小二乘法获得解析解。

• 为方便讨论，我们把$\omega$和$b$吸收入向量形式$\omega =(\omega ;\mathbf{b})$，相应的把数据集$D$表示为一个$m\times (d+1)$大小的矩阵$\mathbf{X}$，其中每行对应于一个示例，该行前d个元素对应于示例的d个属性值，最后一个元素恒置为1，即$$\mathbf{X}=\left(\begin{array}{cccccc}& x_{11}& x_{12}& \cdots & x_{1d}& 1\\ & x_{21}& x_{22}& \cdots & x_{2d}& 1\\ & ⋮& ⋮& \ddots & ⋮& ⋮\\ & x_{m1}& x_{m2}& \cdots & x_{md}& 1\end{array}\right)=\left(\begin{array}{ccc}& x_1^T& 1\\ & x_2^T& 1\\ & ⋮& ⋮\\ & x_m^T& 1\end{array}\right)$$

  再把标记也写成向量形式$\mathbf{y}=(y_1;y_2;\dots ;y_m)$，类似于式(4)，有

  $${\hat{\omega }}^{\ast }=\underset{\hat{\omega }}{\mathrm{argmin}}(\mathbf{y}-\mathbf{X}\hat{\omega }{)}^T(\mathbf{y}-\mathbf{X}\hat{\omega })\text{\hspace{1em}(5)}$$

  令$E_{\hat{\omega }}=(\mathbf{y}-\mathbf{X}\hat{\omega }{)}^T(\mathbf{y}-\mathbf{X}\hat{\omega })$，对$\hat{\omega }$求导得到

  $$\frac{\partial E_{\hat{\omega }}}{\partial \hat{\omega }}=2{\mathbf{X}}^T(\mathbf{X}\hat{\omega }-y)\text{\hspace{1em}(6)}$$

  令上式为0可得$\hat{\omega }$最优解，但由于涉及矩阵逆的计算

• 当${\mathbf{X}}^T\mathbf{X}$为满秩矩阵(full-rank matrix)或正定矩阵(positive definite matrix)时，令式(6)为0可得

  $${\hat{\omega }}^{\ast }=({\mathbf{X}}^T\mathbf{X}{)}^{-1}{\mathbf{X}}^T\mathbf{y}\text{\hspace{1em}(7)}$$

  其中$({\mathbf{X}}^T\mathbf{X}{)}^{-1}$是矩阵${\mathbf{X}}^T\mathbf{X}$的逆矩阵，最终学得的多元线性回归模型为

  $$f(x_i)=x_i^T（{\mathbf{X}}^T\mathbf{X}{)}^{-1}{\mathbf{X}}^T\mathbf{y}\text{\hspace{1em}(8)}$$

• 当${\mathbf{X}}^T\mathbf{X}$不是满秩矩阵，此时存在多个解析解，他们都能使得均方误差最小化。究竟选择哪个解作为输出，由算法的偏好决定。

  比如$N<n$（样本数量小于特征种类的数量），根据$\mathbf{X}$的秩小于等于$N,n$中的最小值，即小于等于$N$（矩阵的秩一定小于等于矩阵的行数和列数）； 而矩阵${\mathbf{X}}^T\mathbf{X}$是$n\times n$大小的，它的秩一定小于等于$N$，因此不是满秩矩阵。

  常见的做法是引入正则化项

  • $L_1$正则化：此时称作Lasso Regression

    $${\hat{\omega }}^{\ast }=\underset{\hat{\omega }}{\mathrm{argmin}}\left[(\mathbf{y}-\mathbf{X}\hat{\omega }{)}^T(\mathbf{y}-\mathbf{X}\hat{\omega })+\lambda {\Vert \hat{\omega }\Vert }_1\right]\text{\hspace{1em}(9)}$$

  • $L_2$正则化：此时称作Ridge Regression

    $${\hat{\omega }}^{\ast }=\underset{\hat{\omega }}{\mathrm{argmin}}\left[(\mathbf{y}-\mathbf{X}\hat{\omega }{)}^T(\mathbf{y}-\mathbf{X}\hat{\omega })+\lambda {\Vert \hat{\omega }\Vert }_2^2\right]\text{\hspace{1em}(10)}$$

  • 同时包含$L_1$，$L_2$正则化，此时称作Elastic Net

    $${\hat{\omega }}^{\ast }=\underset{\hat{\omega }}{\mathrm{argmin}}\left[(\mathbf{y}-\mathbf{X}\hat{\omega }{)}^T(\mathbf{y}-\mathbf{X}\hat{\omega })+\lambda \rho {\Vert \hat{\omega }\Vert }_1+\frac{\lambda (1-\rho )}{2}{\Vert \hat{\omega }\Vert }_2^2\right]\text{\hspace{1em}(11)}$$

    其中

    • $\lambda >0$为正则化系数，调整正则化项与训练误差的比例
    • $0⩽\rho ⩽1$ 为比例系数，调整$L_1$正则化与$L_2$正则化的比例

三、广义线性模型

为便于观察，我们把线性模型简写为

$$y={\omega }^T\mathbf{x}+b\text{\hspace{1em}(12)}$$

如果将输出标记的对数作为线性模型逼近的目标，即

$$lny={\omega }^T\mathbf{x}+b\text{\hspace{1em}(13)}$$

这就是对数线性回归(log-linear regression)

更一般的，考虑单调可微函数$g(\cdot )$，令

$$y=g_{-1}({\omega }^T\mathbf{x}+b)\text{\hspace{1em}(14)}$$

这样得到的模型称为广义线性模型(generalized-linear model)，其中函数$g(\cdot )$称为联系函数(link function)。对数线性回归是广义线性模型在$g(\cdot )=ln(\cdot )$时的特例。

四、公式推导

对于(6)推导

$$\frac{\partial E_{\hat{\omega }}}{\partial \hat{\omega }}=2{\mathbf{X}}^T(\mathbf{X}\hat{\omega }-y)$$

将$E_{\hat{\omega }}=(\mathbf{y}-\mathbf{X}\hat{\omega }{)}^T(\mathbf{y}-\mathbf{X}\hat{\omega })$展开可得

$$E_{\hat{\omega }}={\mathbf{y}}^T\mathbf{y}-{\mathbf{y}}^T\mathbf{X}\omega -{\omega }^T{\mathbf{X}}^T\mathbf{y}+{\omega }^T{\mathbf{X}}^T\mathbf{X}\omega$$

对$\hat{\omega }$求导可得

$$\frac{\partial E_{\hat{\omega }}}{\partial \hat{\omega }}=\frac{\partial {\mathbf{y}}^T\mathbf{y}}{\partial \hat{\omega }}-\frac{\partial {\mathbf{y}}^T\mathbf{X}\omega }{\partial \hat{\omega }}-\frac{\partial {\omega }^T{\mathbf{X}}^T\mathbf{y}}{\partial \hat{\omega }}+\frac{\partial {\omega }^T{\mathbf{X}}^T\mathbf{X}\omega }{\partial \hat{\omega }}$$

由向量的求导公式可得

$$\frac{\partial E_{\hat{\omega }}}{\partial \hat{\omega }}=0-{\mathbf{X}}^T\mathbf{y}-{\mathbf{X}}^T\mathbf{y}+({\mathbf{X}}^T\mathbf{X}+\mathbf{X}{\mathbf{X}}^T)\hat{\omega }$$

$$\frac{\partial E_{\hat{\omega }}}{\partial \hat{\omega }}=2{\mathbf{X}}^T(\mathbf{X}\hat{\omega }-y)$$

$\frac{\partial {\mathbf{X}}^T\mathbf{A}\mathbf{X}}{\partial \mathbf{X}}=(\mathbf{A}+{\mathbf{A}}^T)\mathbf{X}$

五、Python实现代码

import numpy as np
import pandas as pd
from sklearn import preprocessing


class linearRegression:
    def __init__(self):
        self.w = None

    def fit(self, X: np.array, y: np.array) -> np.array:
        '''
        训练函数
        :param X: 特征数据
        :param y: 目标
        :return: 参数w
        '''
        m = len(y)
        y = y.reshape(-1, 1)
        X = np.hstack((X, np.ones((m, 1))))
        X, y = np.mat(X), np.mat(y)
        # 判断是否为满秩矩阵
        if np.linalg.det(X.T * X) == 0:
            print('XTX不是满秩矩阵')
            return
        self.w = (X.T * X).I * X.T * y

    def predict(self, X: np.array) -> np.array:
        '''
        预测函数
        :param X: 特征数据
        :return: 预测值
        '''
        m = X.shape[0]
        X = np.hstack((X, np.ones((m, 1))))
        res = np.zeros((m, 1))
        X = np.mat(X)
        for i in range(m):
            res[i] = X[i] * self.w
        return res


if __name__ == '__main__':
    f_path = './data.csv'
    df = pd.read_csv(f_path)

    array = np.array(df)
    X = array[:, :-1]
    y = array[:, -1]
    X = preprocessing.scale(X)

    X_train, y_train = X[:-10], y[:-10]
    X_test, y_test = X[-10:], y[-10:]

    lr = linearRegression()
    lr.fit(X_train, y_train)
    res = lr.predict(X_test)

    print(res)

代码及文件链接：https://github.com/yourSprite/machinelearninginaction/tree/master/machineLearning/linearRegression

参考

• 《机器学习》
• AI算法工程师手册
• 南瓜书PumpkinBook
• Matrix calculus