线性回归和Python代码
一、问题
1.1 线性模型基本形式
给定d个属性描述的示例 x = ( x 1 ; x 2 ; . . . ; x d ) \mathbf{x} = (x_{1};x_{2};...;x_{d}) x=(x1;x2;...;xd),其中 x i x_{i} xi是 x \mathbf{x} x在第 i i i个属性上的取值,线性模型(linear model)试图学得一个通过属性的线性组合来进行预测的函数,即
f
(
x
)
=
ω
1
x
1
+
ω
2
x
2
+
.
.
.
+
ω
d
x
d
+
b
(1)
f(\mathbf{x}) = \omega_{1}x_{1} + \omega_{2}x_{2}\ +\ ...\ + \ \omega_{d}x_{d} + b \tag{1}
f(x)=ω1x1+ω2x2 + ... + ωdxd+b(1)
一般向量形式写成
f
(
x
)
=
ω
T
x
+
b
(2)
f(\mathbf{x}) = \mathbf{\omega}^{T}\mathbf{x} + b \tag{2}
f(x)=ωTx+b(2)
其中
ω
=
(
ω
1
;
ω
2
;
.
.
.
;
ω
d
)
\mathbf{\omega} = (\omega_{1};\omega_{2};...;\omega_{d})
ω=(ω1;ω2;...;ωd)
1.2 线性回归
- 给定数据集
D
\mathbf{D}
D = {(
x
1
\mathbf{x}_{1}
x1,
y
1
y_{1}
y1),(
x
2
\mathbf{x}_{2}
x2,
y
2
y_{2}
y2),…,(
x
m
\mathbf{x}_{m}
xm,
y
m
y_{m}
ym)},其中
x
i
\mathbf{x}_{i}
xi = {
x
i
1
x_{i1}
xi1;
x
i
2
x_{i2}
xi2;…;
x
i
d
x_{id}
xid},
y
i
y_{i}
yi
∈
\in
∈
R
\mathbb{R}
R
线性回归问题(linear regression)试图学习模型
f ( x ) = ω T x + b f(\mathbf{x}) = \mathbf{\omega}^{T}\mathbf{x} + b f(x)=ωTx+b
对离散属性,若属性值间存在“序”(order)关系,可通过连续化将其转化为连续值,例如二值属性“身高”的取值“高” “矮”可转化为{1.0, 0.0},三值属性高度的取值“高” “中” “低”可转化为{1.0, 0.5, 0.0};若属性间不存在序关系,假定有k个属性值,则通常转化为k维向量,例如“瓜类”的取值“西瓜” “南瓜” “黄瓜” 可转化为(0, 0, 1),(0, 1, 0),(1, 0, 0)
- 使得 f ( x i ) ≃ y i \ f(x_{i}) \simeq y_{i} f(xi)≃yi
该问题也被称作多元线性回归(multivariate linear regression)
- 对于每个
x
i
x_{i}
xi,其预测值
y
i
=
ω
T
x
i
+
b
y_{i} = \mathbf{\omega}^{T}\mathbf{x}_{i} + b
yi=ωTxi+b,采用均方误差(平方损失函数),则在训练集
D
\mathbf{D}
D上,模型的损失函数为:
L ( f ) = ∑ n = 1 N ( f ( x i ) − y i ) 2 = ∑ n = 1 N ( y i − ω T x i − b ) 2 (3) L(f) = \sum_{n=1}^N(f(x_{i}) - y_{i})^{2} =\sum_{n=1}^N(y_{i} - \mathbf{\omega}^{T}\mathbf{x}_{i} - b)^{2} \tag{3} L(f)=n=1∑N(f(xi)−yi)2=n=1∑N(yi−ωTxi−b)2(3)
优化目标是损失函数最小化,即
( ω ∗ , b ∗ ) = a r g min ( ω , b ) ∑ n = 1 N ( y i − ω T x i − b ) 2 (4) (\mathbf{\omega}^{*}, b^{*}) = \underset{(\mathbf{\omega},b)}{arg\min}\sum_{n=1}^N(y_{i} - \mathbf{\omega}^{T}\mathbf{x}_i - b)^{2} \tag{4} (ω∗,b∗)=(ω,b)argminn=1∑N(yi−ωTxi−b)2(4)
二、求解
-
可以用梯度下降法来求解上述最优化问题的数值解,但是实际上该最优化问题可以通过最小二乘法获得解析解。
-
为方便讨论,我们把 ω \mathbf{\omega} ω和 b b b吸收入向量形式 ω = ( ω ; b ) \mathbf{\omega} = (\mathbf{\omega};\mathbf{b}) ω=(ω;b),相应的把数据集 D D D表示为一个 m × ( d + 1 ) m\times(d+1) m×(d+1)大小的矩阵 X \mathbf{X} X,其中每行对应于一个示例,该行前d个元素对应于示例的d个属性值,最后一个元素恒置为1,即 X = ( x 11 x 12 ⋯ x 1 d 1 x 21 x 22 ⋯ x 2 d 1 ⋮ ⋮ ⋱ ⋮ ⋮ x m 1 x m 2 ⋯ x m d 1 ) = ( x 1 T 1 x 2 T 1 ⋮ ⋮ x m T 1 ) \mathbf{X} = \begin{pmatrix}&x_{11}&x_{12}&\cdots&x_{1d}&1 \\ &x_{21}&x_{22}&\cdots &x_{2d}&1 \\ &\vdots&\vdots&\ddots&\vdots &\vdots \\&x_{m1}&x_{m2}&\cdots &x_{md}&1 \end{pmatrix} = \begin{pmatrix}&x^{T}_{1}&1 \\ &x^{T}_{2}&1 \\&\vdots&\vdots \\ &x^{T}_{m}&1\end{pmatrix} X=⎝⎜⎜⎜⎛x11x21⋮xm1x12x22⋮xm2⋯⋯⋱⋯x1dx2d⋮xmd11⋮1⎠⎟⎟⎟⎞=⎝⎜⎜⎜⎛x1Tx2T⋮xmT11⋮1⎠⎟⎟⎟⎞
再把标记也写成向量形式 y = ( y 1 ; y 2 ; … ; y m ) \mathbf{y} = (y_{1};y_{2};…;y_{m}) y=(y1;y2;…;ym),类似于式(4),有
ω ^ ∗ = a r g min ω ^ ( y − X ω ^ ) T ( y − X ω ^ ) (5) \hat{\mathbf{\omega}}^{*} = \underset{\hat{\mathbf{\omega}}}{arg\min}(\mathbf{y} - \mathbf{X}\mathbf{\hat{ \mathbf{\omega}}})^T(\mathbf{y} - \mathbf{X}\mathbf{\hat{\mathbf{\omega}}}) \tag{5} ω^∗=ω^argmin(y−Xω^)T(y−Xω^)(5)
令 E ω ^ = ( y − X ω ^ ) T ( y − X ω ^ ) E_{\hat{\mathbf{\omega}}} = (\mathbf{y} - \mathbf{X}\mathbf{\hat{\mathbf{\omega}}})^T(\mathbf{y} - \mathbf{X}\mathbf{\hat{\mathbf{\omega}}}) Eω^=(y−Xω^)T(y−Xω^),对 ω ^ \hat{\omega} ω^求导得到
∂ E ω ^ ∂ ω ^ = 2 X T ( X ω ^ − y ) (6) \frac{\partial E_{\hat{\mathbf{\omega}}}}{\partial \hat{\mathbf{\omega}}} = 2\mathbf{X}^T(\mathbf{X\hat{\omega}}-y) \tag{6} ∂ω^∂Eω^=2XT(Xω^−y)(6)
令上式为0可得 ω ^ \hat{\mathbf{\omega}} ω^最优解,但由于涉及矩阵逆的计算
-
当 X T X \mathbf{X}^T\mathbf{X} XTX为满秩矩阵(full-rank matrix)或正定矩阵(positive definite matrix)时,令式(6)为0可得
ω ^ ∗ = ( X T X ) − 1 X T y (7) \hat{\mathbf{\omega}}^{*} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y} \tag{7} ω^∗=(XTX)−1XTy(7)
其中 ( X T X ) − 1 (\mathbf{X}^T\mathbf{X})^{-1} (XTX)−1是矩阵 X T X \mathbf{X}^T\mathbf{X} XTX的逆矩阵,最终学得的多元线性回归模型为
f ( x i ) = x i T ( X T X ) − 1 X T y (8) f(x_{i}) = x_{i}^T(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y} \tag{8} f(xi)=xiT(XTX)−1XTy(8)
-
当 X T X \mathbf{X}^T\mathbf{X} XTX不是满秩矩阵,此时存在多个解析解,他们都能使得均方误差最小化。究竟选择哪个解作为输出,由算法的偏好决定。
比如 N < n N < n N<n(样本数量小于特征种类的数量),根据 X \mathbf{X} X的秩小于等于 N , n N,n N,n中的最小值,即小于等于 N N N(矩阵的秩一定小于等于矩阵的行数和列数); 而矩阵 X T X \mathbf{X}^T\mathbf{X} XTX是 n × n n\times n n×n大小的,它的秩一定小于等于 N N N,因此不是满秩矩阵。
常见的做法是引入正则化项
-
L 1 L_{1} L1正则化:此时称作Lasso Regression
ω ^ ∗ = a r g min ω ^ [ ( y − X ω ^ ) T ( y − X ω ^ ) + λ ∥ ω ^ ∥ 1 ] (9) \hat{\mathbf{\omega}}^{*} = \underset{\hat{\mathbf{\omega}}}{arg\min}\left [ (\mathbf{y} - \mathbf{X}\mathbf{\hat{\omega}})^T(\mathbf{y} - \mathbf{X}\mathbf{\hat{\omega}})+ \lambda\left \| \mathbf{\hat{\omega}}\right \|_{1}\right ] \tag{9} ω^∗=ω^argmin[(y−Xω^)T(y−Xω^)+λ∥ω^∥1](9)
-
L 2 L_{2} L2正则化:此时称作Ridge Regression
ω ^ ∗ = a r g min ω ^ [ ( y − X ω ^ ) T ( y − X ω ^ ) + λ ∥ ω ^ ∥ 2 2 ] (10) \hat{\mathbf{\omega}}^{*} = \underset{\hat{\mathbf{\omega}}}{arg\min}\left [ (\mathbf{y} - \mathbf{X}\mathbf{\hat{\omega}})^T(\mathbf{y} - \mathbf{X}\mathbf{\hat{\omega}})+ \lambda\left \| \mathbf{\hat{\omega}}\right \|_{2}^{2}\right ] \tag{10} ω^∗=ω^argmin[(y−Xω^)T(y−Xω^)+λ∥ω^∥22](10)
-
同时包含 L 1 L_{1} L1, L 2 L_{2} L2正则化,此时称作Elastic Net
ω ^ ∗ = a r g min ω ^ [ ( y − X ω ^ ) T ( y − X ω ^ ) + λ ρ ∥ ω ^ ∥ 1 + λ ( 1 − ρ ) 2 ∥ ω ^ ∥ 2 2 ] (11) \hat{\mathbf{\omega}}^{*} = \underset{\hat{\mathbf{\omega}}}{arg\min}\left [ (\mathbf{y} - \mathbf{X}\mathbf{\hat {\omega}})^T(\mathbf{y} - \mathbf{X}\mathbf{\hat{\omega}})+ \lambda\rho\left \| \mathbf{\hat{\omega}}\right \|_{1} + \frac{\lambda(1-\rho)}{2}\left \| \mathbf{\hat{\omega}}\right \|_{2}^{2}\right ] \tag{11} ω^∗=ω^argmin[(y−Xω^)T(y−Xω^)+λρ∥ω^∥1+2λ(1−ρ)∥ω^∥22](11)
其中
- λ > 0 \lambda > 0 λ>0为正则化系数,调整正则化项与训练误差的比例
- 0 ⩽ ρ ⩽ 1 0\leqslant\rho\leqslant1 0⩽ρ⩽1 为比例系数,调整 L 1 L_{1} L1正则化与 L 2 L_{2} L2正则化的比例
-
三、广义线性模型
为便于观察,我们把线性模型简写为
y = ω T x + b (12) y = \mathbf{\omega}^{T}\mathbf{x} + b \tag{12} y=ωTx+b(12)
如果将输出标记的对数作为线性模型逼近的目标,即
l n y = ω T x + b (13) lny = \mathbf{\omega}^{T}\mathbf{x} + b \tag{13} lny=ωTx+b(13)
这就是对数线性回归(log-linear regression)
更一般的,考虑单调可微函数 g ( ⋅ ) g(·) g(⋅),令
y = g − 1 ( ω T x + b ) (14) y = g_{-1}(\mathbf{\omega}^{T}\mathbf{x} + b) \tag{14} y=g−1(ωTx+b)(14)
这样得到的模型称为广义线性模型(generalized-linear model),其中函数 g ( ⋅ ) g(·) g(⋅)称为联系函数(link function)。对数线性回归是广义线性模型在 g ( ⋅ ) = l n ( ⋅ ) g(·) = ln(·) g(⋅)=ln(⋅)时的特例。
四、公式推导
对于(6)推导
∂ E ω ^ ∂ ω ^ = 2 X T ( X ω ^ − y ) \frac{\partial E_{\hat{\mathbf{\omega}}}}{\partial \hat{\mathbf{\omega}}} = 2\mathbf{X}^T(\mathbf{X\hat{\omega}}-y) ∂ω^∂Eω^=2XT(Xω^−y)
将 E ω ^ = ( y − X ω ^ ) T ( y − X ω ^ ) E_{\hat{\mathbf{\omega}}} = (\mathbf{y} - \mathbf{X}\mathbf{\hat{\omega}})^T(\mathbf{y} - \mathbf{X}\mathbf{\hat{\omega}}) Eω^=(y−Xω^)T(y−Xω^)展开可得
E ω ^ = y T y − y T X ω − ω T X T y + ω T X T X ω E_{\hat{\mathbf{\omega}}} = \mathbf{y}^T\mathbf{y} - \mathbf{y}^T\mathbf{X}\mathbf{\omega} - \mathbf{\omega}^T\mathbf{X}^T\mathbf{y} + \mathbf{\omega}^T\mathbf{X}^T\mathbf{X}\mathbf{\omega} Eω^=yTy−yTXω−ωTXTy+ωTXTXω
对 ω ^ \hat{\omega} ω^求导可得
∂ E ω ^ ∂ ω ^ = ∂ y T y ∂ ω ^ − ∂ y T X ω ∂ ω ^ − ∂ ω T X T y ∂ ω ^ + ∂ ω T X T X ω ∂ ω ^ \frac{\partial E_{\hat{\mathbf{\omega}}}}{\partial \hat{\mathbf{\omega}}} = \frac{\partial \mathbf{y}^T\mathbf{y}}{\partial \hat{\mathbf{\omega}}} - \frac{\partial \mathbf{y}^T\mathbf{X}\mathbf{\omega}}{\partial \hat{\mathbf{\omega}}} - \frac{\partial \mathbf{\omega}^T\mathbf{X}^T\mathbf{y}}{\partial \hat{\mathbf{\omega}}} + \frac{\partial \mathbf{\omega}^T\mathbf{X}^T\mathbf{X}\mathbf{\omega}}{\partial \hat{\mathbf{\omega}}} ∂ω^∂Eω^=∂ω^∂yTy−∂ω^∂yTXω−∂ω^∂ωTXTy+∂ω^∂ωTXTXω
由向量的求导公式可得
∂ E ω ^ ∂ ω ^ = 0 − X T y − X T y + ( X T X + X X T ) ω ^ \frac{\partial E_{\hat{\mathbf{\omega}}}}{\partial \hat{\mathbf{\omega}}} = \mathbf{0} - \mathbf{X}^T\mathbf{y} - \mathbf{X}^T\mathbf{y} + (\mathbf{X}^T\mathbf{X} + \mathbf{X}\mathbf{X}^T)\hat{\mathbf{\omega}} ∂ω^∂Eω^=0−XTy−XTy+(XTX+XXT)ω^
∂ E ω ^ ∂ ω ^ = 2 X T ( X ω ^ − y ) \frac{\partial E_{\hat{\mathbf{\omega}}}}{\partial \hat{\mathbf{\omega}}} = 2\mathbf{X}^T(\mathbf{X\hat{\omega}}-y) ∂ω^∂Eω^=2XT(Xω^−y)
∂ X T A X ∂ X = ( A + A T ) X \frac{\partial \mathbf{X}^T\mathbf{A}\mathbf{X}}{\partial \mathbf{X}} = (\mathbf{A} + \mathbf{A}^T)\mathbf{X} ∂X∂XTAX=(A+AT)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
参考