机器学习01:线性回归

线性回归的基本概念

线性模型的基本形式

给定由$n$个属性描述的示例$\mathbf{x}=(x_1;x_2;...;x_n)$,其中$x_i$是$\mathbf{x}$在第$i$个属性上的取值,线性模型试图学得一个通过属性的线性组合来进行预测的函数,即
$$h(\mathbf{x})=w_1{x}_1+w_2{x}_2+...+w_n{x}_n+b$$

定义权重向量$\mathbf{w}=(w_1;w_2;...;w_n)$,表示每个属性在预测结果中所占的权重,上式可以写成
$$h(\mathbf{x})={\mathbf{w}}^T\mathbf{x}+b$$

线性回归的最终目的,就是要求得权重向量$\mathbf{w}$和截距$b$.

线性模型的优点

线性回归(linear regreression)是机器学习里面最基本的模型,它主要的优点有如下两条:

1. 线性模型形式简单,易于建模,但却蕴含着机器学习中的一些重要的基本思想,许多功能更为强大的非线性模型可在线性模型的基础上引入层级结构或高维映射得到.
2. 线性模型具有很好的可解释性(comprehensibility),器权重向量可以直观表达各属性在预测中的重要性.

线性回归的损失函数

损失函数(又叫误差函数)反映了预测结果与实际结果之间的差别.对于线性回归模型,其损失函数为均方误差.
$$\begin{array}{rl}J(\mathbf{w})& =\sum _{i=1}^m(h({\mathbf{x}}_i)-{\mathbf{y}}_i{)}^2\\ & =\sum _{i=1}^m({\mathbf{y}}_i-{\mathbf{w}}^T{\mathbf{x}}_i-b{)}^2\end{array}$$

下面是该损失函数的推导过程:

定义${\mathbf{\epsilon }}_i=f({\mathbf{x}}_i)-{\mathbf{y}}_i$,表示预测值与实际值之间的误差.根据线性回归模型的基本假设,$\epsilon$是独立同分布的,根据中心极限定理,$\epsilon$应该服从均值为$0$,方差为${\sigma }^2$的高斯分布.

将误差的定义式带入高斯分布的概率密度函数
$$p({\mathbf{\epsilon }}_i)=\frac{1}{\sqrt{2\pi }\sigma }\exp \left(-\frac{{\epsilon }_i^2}{2{\sigma }^2}\right)$$

即
$$p({\mathbf{y}}_i|{\mathbf{x}}_i;\mathbf{w})=\frac{1}{\sqrt{2\pi }\sigma }\exp \left(-\frac{({\mathbf{y}}_i-h({\mathbf{x}}_i){)}^2}{2{\sigma }^2}\right)$$

由于各个样本是独立的,它们的联合概率密度就是各自的概率密度的乘积,即似然函数
$$\begin{array}{rl}L(\mathbf{w})& =\prod _{i=1}^{m}p({\mathbf{y}}_i|{\mathbf{x}}_i;\mathbf{w})\\ & =\prod _{i=1}^m\frac{1}{\sqrt{2\pi }\sigma }\exp \left(-\frac{{\left({\mathbf{y}}_i-{\mathbf{w}}^T{\mathbf{x}}_i\right)}^2}{2{\sigma }^2}\right)\end{array}$$

对上式取对数,得到
$$\begin{array}{rl}\ell (\mathbf{w})& =\log L(\mathbf{w})\\ & =\log \prod _{i=1}^m\frac{1}{\sqrt{2\pi }\sigma }\exp \left(-\frac{({\mathbf{y}}_i-{\mathbf{w}}^T{\mathbf{x}}_i{)}^2}{2{\sigma }^2}\right)\\ & =\sum _{i=1}^m\log \frac{1}{\sqrt{2\pi }\sigma }\exp \left(-\frac{({\mathbf{y}}_i-{\mathbf{w}}^T{\mathbf{x}}_i{)}^2}{2{\sigma }^2}\right)\\ & =m\log \frac{1}{\sqrt{2\pi }\sigma }-\frac{1}{{\sigma }^2}\cdot \frac{1}{2}\sum _{i=1}^m{\left({\mathbf{y}}_i-{\mathbf{w}}^T{\mathbf{x}}_i\right)}^2\end{array}$$

要使得最大似然函数$L(\mathbf{w})$最大,就要使$J(\mathbf{w})=\sum _{i=1}^m({\mathbf{y}}_i-{\mathbf{w}}^T{\mathbf{x}}_i{)}^2$最小

线性回归的拟合方法

线性回归最常见的拟合方法有两种:最小二乘法(least square method)和梯度下降法(gradient descent).

最小二乘法

根据数学推导,权重矩阵$\mathbf{w}$的最优解
$$\mathbf{w}=({\mathbf{X}}^T\mathbf{X}{)}^{-1}{\mathbf{X}}^T\mathbf{y}$$

最小二乘法的理论推导

为便于使用最小二乘法对$\mathbf{w}$和$b$进行估计,我们把$\mathbf{w}$和$b$纳入向量形式$\hat{\mathbf{w}}=(\mathbf{w};b)$,相应的,把数据集$D$表示为一个$d\times (n+1)$大小的矩阵$\mathbf{X}$,其中每行对应一个示例前$n$个元素对应于示例的$n$个属性值,最后一个元素恒置为1,即
$$\mathbf{X}=\left(\begin{array}{ccccc}{x}_{11}& x_{12}& \cdots & x_{1n}& 1\\ x_{21}& x_{22}& \cdots & x_{2n}& 1\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ x_{m1}& x_{m2}& \cdots & x_{mn}& 1\end{array}\right)=\left(\begin{array}{cc}{\mathbf{x}}_1^T& 1\\ {\mathbf{x}}_2^T& 1\\ ⋮& ⋮\\ {\mathbf{x}}_m^T& 1\end{array}\right)$$

再把标记也写成向量形式$\mathbf{y}=(y_1;y_2;...;y_m)$,将$J(\hat{\mathbf{w}})$写成向量的形式,有
$$\begin{array}{rl}J(\hat{\mathbf{w}})& =\sum _{i=1}^m(h({\mathbf{x}}_i)-{\mathbf{y}}_i{)}^2\\ & =(\mathbf{X}\hat{\mathbf{w}}-\mathbf{y}{)}^T(\mathbf{X}\hat{\mathbf{w}}-\mathbf{y})\end{array}$$

将$J(\hat{\mathbf{w}})$对$\hat{\mathbf{w}}$求导
$$\begin{array}{rl}\frac{\partial J(\hat{\mathbf{w}})}{\partial \hat{\mathbf{w}}}& =2{\mathbf{X}}^T(\mathbf{X}\hat{\mathbf{w}}-\mathbf{y})\end{array}$$

因为$J(\hat{\mathbf{w}})$是一个凸函数,其最小值取值在拐点处.故令上式结果等于$0$,得到
$$\hat{\mathbf{w}}=({\mathbf{X}}^T\mathbf{X}{)}^{-1}{\mathbf{X}}^T\mathbf{y}$$

最小二乘法的代码实现

使用Python代码实现最小二乘法:

import numpy as np

# 创建X1,X2,Y为数据集,其中的点大致满足 Y=3*X1+2*X2+4 的规律
X1 = np.random.rand(100, 1)
X2 = np.random.rand(100, 1)
Y = 3 * X1 + 2 * X2 + 4 + np.random.randn(100, 1)
X_b = np.c_[X1, X2, np.ones((100, 1))]  # 构建增广X矩阵

# 使用解析法求解omega
omega_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(Y)
print(omega_best)

计算得到$\hat{\mathbf{w}}=[3.71349706,2.27034494,3.48867339]$.

梯度下降法

梯度下降法的原理

梯度即函数下降速度最快的方向,梯度在各方向上的分量值为函数对该分量的导数.设二元函数$z=f(x,y)$在平面区域$D$上具有一阶连续偏导数,则向量$\left\{\frac{\partial f}{\partial x},\frac{\partial f}{\partial y}\right\}=f_x(x,y)\cdot \mathbf{i}+f_y(x,y)\cdot \mathbf{j}$为$f(x,y)$的梯度,记为$gradf(x,y)$或$\nabla f(x,y)$.即
$$gradf(x,y)=\nabla f(x,y)=\left\{\frac{\partial f}{\partial x},\frac{\partial f}{\partial y}\right\}=f_x(x,y)\cdot \mathbf{i}+f_y(x,y)\cdot \mathbf{j}$$

梯度下降法即通过对线性回归的损失函数$J(\mathbf{w})$求梯度,每次迭代都使$\hat{\mathbf{w}}$沿着梯度方向下降,最终能找到函数的局部最小值;因为该损失函数是一个凸函数,我们找到的也就是损失函数的全局最小值.

梯度下降法的步骤如下:

1. 初始化${\hat{\mathbf{w}}}_0=(w_1;w_2;...;w_n)$
2. 计算当前${\hat{\mathbf{w}}}_t$处的梯度$\nabla J({\hat{\mathbf{w}}}_t)=(\frac{\partial J}{\partial w_1},\frac{\partial J}{\partial w_2},...,\frac{\partial J}{\partial w_n})$
3. 调整${\hat{\mathbf{w}}}_{t+1}={\hat{\mathbf{w}}}_t-\alpha \cdot \nabla J({\hat{\mathbf{w}}}_t$,其中$\alpha$为学习率,是一个超参数,它的选取可以影响迭代所需次数:若$\alpha$过大,则会发生振荡,若$\alpha$过小,则算法收敛过慢.
4. 重复第2,3步直到梯度$\nabla J({\hat{\mathbf{w}}}_t)$小于某阈值${\nabla }_m$,该阈值${\nabla }_m$也是一个超参数,同样会影响算法的收敛.

梯度下降法的代码实现

梯度下降法有三种常见实现: 分别为批量梯度下降算法,随机梯度下降算法和小批量梯度下降算法.它们的区别别在于迭代过程中计算梯度时所使用的样本条数不同.

批量梯度下降算法

批量梯度下降算法(batch gradient descent)在计算梯度时,先对每条样本均求取梯度,最后选取其平均值代入公式.即
$${\hat{\mathbf{w}}}_{t+1}={\hat{\mathbf{w}}}_t-\frac{1}{m}\alpha {\mathbf{X}}^T(\mathbf{X}\widehat{{\mathbf{w}}_t}-\mathbf{y})$$

其优点为计算准确,更可能收敛到最优解,缺点为计算量较大.

import numpy as np

# 创建X1,X2,Y为数据集,其中的点大致满足 Y=3*X1+2*X2+4 的规律
m = 100    # 样本数
X1 = np.random.rand(m, 1)
X2 = np.random.rand(m, 1)
Y = 3 * X1 + 2 * X2 + 4
X_b = np.c_[X1, X2, np.ones((m, 1))]  # 构建增广X矩阵

# 使用批量梯度下降算法求解omega
learning_rate = 0.1     # 定义学习率
iteration_num = 1000    # 定义迭代轮数

omega = np.random.randn(3, 1)   # 初始化权重向量
for i in range(iteration_num):  # 迭代
    # 计算所有样本梯度的均值
    gradient = 1 / m * X_b.T.dot(X_b.dot(omega) - Y)
    omega = omega - learning_rate * gradient

print(omega)

可以看到,随着迭代轮数的增加,准确率升高.

随机梯度下降算法

随机梯度下降算法(stochastic gradient descent)在计算梯度时,随机取一个样本,对其计算梯度并代入公式.即
$${\hat{\mathbf{w}}}_{t+1}={\hat{\mathbf{w}}}_t-\frac{1}{m}\alpha {\mathbf{x}}_i^T({\mathbf{x}}_i\widehat{{\mathbf{w}}_t}-y_i)$$

其优点为计算量小,且有时可以跳出局部最小值.

import numpy as np

# 创建X1,X2,Y为数据集,其中的点大致满足 Y=3*X1+2*X2+4 的规律
m = 100     # 样本数
X1 = np.random.rand(m, 1)
X2 = np.random.rand(m, 1)
Y = 3 * X1 + 2 * X2 + 4
X_b = np.c_[X1, X2, np.ones((m, 1))]  # 构建增广X矩阵

# 使用随机梯度下降算法求解omega
learning_rate = 0.1  # 定义学习率
iteration_num = 1000  # 定义迭代轮数

omega = np.random.randn(3, 1)  # 初始化权重向量
for i in range(iteration_num):  # 迭代
    # 选取随机一个样本计算梯度
    random_index = np.random.randint(m)
    x_i = X_b[random_index:random_index + 1]
    y_i = Y[random_index:random_index + 1]
    gradient = x_i.T * (x_i.dot(omega) - y_i)
    omega = omega - learning_rate * gradient

print(omega)

下图直观展示了批量梯度下降算法和随机梯度下降算法的区别:(左为批量梯度下降算法,右为随机梯度下降算法)

小批量梯度下降算法

小批量梯度下降算法(mini-batch gradient descent)计算题梯时的策略介于批量梯度下降算法和随机梯度下降算法之间,选取样本中的一部分进行计算求取梯度.

过拟合

过拟合(overfitting)是指训练出的模型过于复杂,导致你和函数完美的符合训练集,但泛化能力较差,不能推广到新的数据集.

解决过拟合的办法:正则化

${\ell }_1$正则和${\ell }_2$正则

在线性模型中,为了防止过拟合,我们希望所有属性的权重的绝对值$||w||$越小越好,因此,我们使用权重构造惩罚项,加在损失函数中,可以有效防止过拟合现象.我们可以构造出两种惩罚项,分别为${\ell }_1$正则(${\ell }_1$-norm)和${\ell }_2$正则(${\ell }_2$-norm).

• ${\ell }_1$正则: 所有维度上的权重$w_i$的绝对值之和.
  $$||w|{|}_1=\frac{1}{m}\sum _{i=1}^n{w}_i$$

• ${\ell }_2$正则: 所有维度上的权重$w_i$的平方之和.
  $$||w|{|}_2=\frac{1}{m}\sum _{i=1}^n{w}_i^2$$

${\ell }_1$正则和${\ell }_2$正则的作用

• 使用${\ell }_1$正则,可以使得到权重向量$\mathbf{w}$结果稀疏化,即有的$w_i$接近于0,而有的$w_i$接近于1.这可以用来降维.
• 使用${\ell }_2$正则,可以使得权重向量$\mathbf{w}$的结果整体偏小.

下面以两个维度的情况为例,说明两种正则为何起到不同的作用:

对于带有${\ell }_1$和${\ell }_2$正则的损失函数,等价于下面两种形式:
$$\begin{gathered} \underset{w}{\min }\frac{1}{n}||y-Xw|{|}^2,st||w|{|}_1<c \\ \underset{w}{\min }\frac{1}{n}||y-Xw|{|}^2,st||w|{|}_2<c \end{gathered}$$

我们将($w_1$,$w_2$)的解空间画出,其中彩色的线是最小二乘表达式的等值线,黑色的线是正则约束的等值线,彩色线与黑色线相间的地方即为最优解.

因为${\ell }_1$正则约束性项的等值线带有棱角,所以最优解很可能落在棱角上;而${\ell }_2$正则正相反.

带有正则项的回归

通过在损失函数上加以不同的正则回归项,我们可以构造出三种回归模型: 岭回归(Ridge Regression,使用${\ell }_2$正则项),Lasso回归(Lasso Regression,使用${\ell }_1$正则项),Elastic-Net回归(Elastic-Net Regression,使用${\ell }_1$和${\ell }_2$正则项).

岭回归

岭回归的损失函数如下:
$$\underset{w}{\min }||Xw-y|{|}_2^2+\alpha ||w|{|}_2^2$$

使用sklearn求解岭回归的代码如下:

from sklearn import linear_model

reg = linear_model.Ridge(alpha=.5)	# 设置超参数alpha:l2正则的权重
reg.fit([[0, 0], [0, 0], [1, 1]], [0, .1, 1]) 
# Ridge(alpha=0.5, copy_X=True, fit_intercept=True, max_iter=None,
#      normalize=False, random_state=None, solver='auto', tol=0.001)

reg.coef_
# array([0.34545455, 0.34545455])

reg.intercept_ 
# 0.13636...

Lasso回归

Lasso回归的损失函数如下:
$$\underset{w}{\min }\frac{1}{2n_{samples}}||Xw-y|{|}_2^2+\alpha ||w|{|}_1$$

使用sklearn求解Lasso回归的代码如下:

from sklearn import linear_model
reg = linear_model.Lasso(alpha=0.1)	# 设置超参数alpha:l1正则的权重
reg.fit([[0, 0], [1, 1]], [0, 1])  
# Lasso(alpha=0.1, copy_X=True, fit_intercept=True, max_iter=1000,
#   normalize=False, positive=False, precompute=False, random_state=None,
#   selection='cyclic', tol=0.0001, warm_start=False)
reg.predict([[1, 1]])
# array([0.8])

Elastic-Net回归

Elastic-Net回归的损失函数如下:
$$\underset{w}{\min }\frac{1}{2n_{samples}}||Xw-y|{|}_2^2+\alpha \rho ||w|{|}_1+\frac{\alpha (1-\rho )}{2}||w|{|}_2^2$$

多项式回归

可以通过升维,使用线性回归来求解多项式回归.

对于下述多项式:
$$\hat{y}(w,x)=w_0+w_1{x}_1+w_2{x}_2+w_3{x}_1{x}_2+w_4{x}_1^2+w_5{x}_2^2$$

我们可以将$x$向量升维以构造$z$向量
$$z=[x_1;x_2;x_1{x}_2;x_1^2;x_2^2]$$

则$y$与$z$成线性关系,可以使用线性回归求解
$$\hat{y}(w,x)=w_0+w_1{z}_1+w_2{z}_2+w_3{z}_3+w_4{z}_4+w_5{z}_5$$

我们可以使用sklearn.preprocessing.PolynomialFeatures类实现原始数据的升维,示例如下:

from sklearn.preprocessing import PolynomialFeatures
import numpy as np
X = np.arange(6).reshape(3, 2)
# array([[0, 1],
#        [2, 3],
#        [4, 5]])

poly = PolynomialFeatures(degree=2)
poly.fit_transform(X)
# array([[ 1.,  0.,  1.,  0.,  0.,  1.],
#        [ 1.,  2.,  3.,  4.,  6.,  9.],
#        [ 1.,  4.,  5., 16., 20., 25.]])

使用升维进行多项式回归的完整代码如下:

import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression

# 构造数据集
X = 6 * np.random.rand(100, 1) - 3
y = 2 + X + 0.5 * X ** 2 + np.random.randn(100, 1)
plt.plot(X, y, 'b.')

degree_color = {1: 'g-', 2: 'r+', 10: 'y*'}
# 分别以degree=1,2,10进行多项式你和
for degree in degree_color.keys():
    # 将样本数据升维
    poly_features = PolynomialFeatures(degree=degree, include_bias=False)  # 指定生成的高维数据不包含0次项
    X_poly = poly_features.fit_transform(X)

    # 对升维后的数据使用线性回归
    lin_reg = LinearRegression(fit_intercept=True)  # 指定拟合斜率
    lin_reg.fit(X_poly, y)
    print(lin_reg.intercept_, lin_reg.coef_)

    # 输出并绘图
    y_predict = lin_reg.predict(X_poly)
    plt.plot(X, y_predict, degree_color[degree])

plt.legend(["degree=%d" % degree for degree in degree_color.keys()])
plt.show()

得到拟合结果如下:

degree=1:
intercept_=[3.35004007],coef_=[[0.84419104]]

degree=2:
intercept_=[1.98123375],coef_=[[1.00078942 0.5193602 ]]

degree=10:
intercept_=[1.83489589],coef_=[[ 1.36019108 0.26399949 -0.67077483 0.60043164 0.25872215 -0.25575867 -0.03607122 0.03928163 0.00170283 -0.00200932]]