本文介绍了线性回归的用途、表达形式及其在机器学习中的应用。线性回归模型简单、可解释性强，适合预测连续输出。通过特征工程扩展输入，可以提升模型能力。在多输出情况下，线性回归的参数矩阵各列可独立优化。
写在最前面
这几年算法岗位的求职供远大于求，可谓是神仙打架，对于数学基础的考察越来越全面和深入，笔者在独角兽和体制内都工作过，参加过近百场面试，也辅导过数十位同学找到心仪的算法工作，想通过这一系列的教程体现机器学习算法中的数学原理，基本可以涵盖所有算法面试中的难题难点，提升面试时的成功率。
那我们就首先从线性回归开始我们的学习。
Why We Use It
线性模型是一种在pre computer age和当前都很有效的模型，形式简单又能保证一定的表达能力，同时提供了可解释性，是函数的一阶泰勒展开。从某种意义上说，简单的形式也保证了较好的抗过拟合能力。当输出连续的预测值时，此种线性模型被称为线性回归模型。
除了上述优点，对线性回归模型输入进行transform，可以将模型的capacity大大扩展。
The Expression
对于一个函数 f f f，若 f ( ∑ i α i x i ) = ∑ i α i f ( x i ) f(\sum_i \alpha_ix_i) = \sum_i\alpha_if(x_i) f(∑iαixi)=∑iαif(xi)，则称此函数是一个线性函数，进一步可据此得到其表达形式
记输入为
X
T
=
{
X
1
,
X
2
,
.
.
.
X
p
}
(1)
X^T = \{X_1, X_2, ...X_p\} \tag{1}
XT={X1,X2,...Xp}(1)
输出为
Y
∈
R
N
Y \in R^N
Y∈RN，则有
Y = ∑ k = 1 P β k × X k + β 0 (2) Y=\sum_{k=1}^{P}\beta_k\times X_k + \beta_0 \tag{2} Y=k=1∑Pβk×Xk+β0(2)
其中
β
0
\beta_0
β0又被称为interception，截距项，可以通过的简单扩展，将原始的输入数据改写为
X
′
=
(
1
;
X
)
X'=(1;X)
X′=(1;X)，则原式改写为
Y
=
∑
k
=
0
P
β
k
×
X
k
(3)
Y=\sum_{k=0}^{P}\beta_k\times X_k \tag{3}
Y=k=0∑Pβk×Xk(3)
这一步是一个简单的数学操作，进行形式上的统一。
为了行文方便，下文中都是用 X X X代替 X ′ X' X′，在很多场景下，得到的输入 X X X都是包含了截距项的。
给定一个数据集 X X X，使用其对模型的参数进行估计，每一行是一个单独的sample。若记每一行的sample为 X i X_i Xi，记第 i i i个样本的第 j j j个特征取值为 X i , j X_{i,j} Xi,j，则 X X X可表示为 X ∈ R N × ( P + 1 ) X \in R^{N\times(P+1)} X∈RN×(P+1)。记
X = [ X 1 , 0 X 1 , 1 X 1 , 2 ⋯ X 1 , p X 2 , 0 X 2 , 1 X 2 , 2 ⋯ X 2 , p ⋮ ⋮ ⋮ ⋱ ⋮ X N , 0 X N , 1 X N , 2 ⋯ X N , p ] (4) X = \left[ \begin{matrix} X_{1,0}&X_{1,1}&X_{1,2}&\cdots &X_{1,p}\\ X_{2,0}&X_{2,1}&X_{2,2}&\cdots &X_{2,p}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ X_{N,0}&X_{N,1}&X_{N,2}&\cdots &X_{N,p} \end{matrix} \right] \tag{4} X=⎣⎢⎢⎢⎡X1,0X2,0⋮XN,0X1,1X2,1⋮XN,1X1,2X2,2⋮XN,2⋯⋯⋱⋯X1,pX2,p⋮XN,p⎦⎥⎥⎥⎤(4)
且
Y
i
=
⟨
X
i
,
β
⟩
=
∑
j
=
0
P
β
j
X
i
j
(5)
Y_i= \langle X_i, \beta\rangle=\sum_{j=0}^{P} \beta_jX_{ij} \tag{5}
Yi=⟨Xi,β⟩=j=0∑PβjXij(5)
则
Y
^
=
(
Y
1
,
Y
2
.
.
.
Y
N
)
T
=
(
⟨
X
1
,
β
⟩
,
⟨
X
2
,
β
⟩
,
.
.
.
⟨
X
N
,
β
⟩
)
T
=
[
∑
j
β
j
X
1
j
∑
j
β
j
X
2
j
⋮
∑
j
β
j
X
N
j
]
(6)
\begin{aligned} \hat{Y} &= (Y_1, Y_2...Y_N)^T \\ &= (\langle X_1, \beta\rangle, \langle X_2, \beta\rangle,...\langle X_N, \beta\rangle)^T \\ &= \left[ \begin{matrix} \sum_j \beta_jX_{1j}\\ \sum_j \beta_jX_{2j}\\ \vdots\\ \sum_j \beta_jX_{Nj} \end{matrix} \right] \end{aligned} \tag{6}
Y^=(Y1,Y2...YN)T=(⟨X1,β⟩,⟨X2,β⟩,...⟨XN,β⟩)T=⎣⎢⎢⎢⎡∑jβjX1j∑jβjX2j⋮∑jβjXNj⎦⎥⎥⎥⎤(6)
即预测值
Y
^
=
X
β
∈
R
N
(7)
\hat{Y} =X\beta \in R^{N} \tag{7}
Y^=Xβ∈RN(7)
所有data driven的任务都对数据有要求，此时假设数据集中的sample都是独立同分布的(i.i.d)，这是为了保证采样的数据集来自同一个分布且互不干扰。
为了衡量其与真实值的差异大小，最常见的指标会选择RSS(residual sum of squares)，首先计算各个分量的差异，
(
Y
i
−
Y
^
i
)
2
(Y_i - \hat{Y}_i)^2
(Yi−Y^i)2，将每个维度上的误差进行累加，使用此指标作为目标函数的优化算法被称为最小二乘法。
此外，从方程求解的角度也可以对RSS进行理解和说明。
若存在 m m m个方程，即对于 ∀ i ∈ [ 1 , m ] , ∑ j a i j x j = b i \forall{i} \in [1,m], \sum_j a_{ij}x_j = b_i ∀i∈[1,m],∑jaijxj=bi。从几何的角度来看，这也是定义了 m m m个超平面，同时满足约束的点即为交点的集合。
当矩阵
A
=
[
a
11
,
a
12
,
a
13
.
.
.
a
1
n
a
21
,
a
22
,
a
23
.
.
.
a
2
n
⋮
a
m
1
,
a
m
2
,
a
m
3
.
.
.
a
m
n
]
(8)
A = \left[ \begin{matrix} a_{11}, a_{12}, a_{13}...a_{1n} \\ a_{21}, a_{22}, a_{23}...a_{2n} \\ \vdots \\ a_{m1}, a_{m2}, a_{m3}...a_{mn} \end{matrix} \right] \tag{8}
A=⎣⎢⎢⎢⎡a11,a12,a13...a1na21,a22,a23...a2n⋮am1,am2,am3...amn⎦⎥⎥⎥⎤(8)
则上述约束可写为
A
x
=
b
(9)
Ax = b \tag{9}
Ax=b(9)
当 A A A满足 m = n m = n m=n且满秩时，此矩阵可逆，则交点为 x = A − 1 b x = A^{-1}b x=A−1b。
若
m
>
n
m > n
m>n，则这个情况被称为过约束，很可能不存在交点，则将求取的目标由交点改为到每个超平面的距离之和最小的点
x
′
x'
x′，到任意超平面的距离为
∣
∑
j
a
i
j
x
j
′
−
b
i
∣
∑
j
a
i
j
2
(10)
\frac{|\sum_j a_{ij}x_j' - b_i|}{\sqrt{\sum_j a_{ij}^2}} \tag{10}
∑jaij2∣∑jaijxj′−bi∣(10)
当超平面参数单位化，并使用平方误差代替绝对值之后，得到的优化目标等价于RSS。
进一步考察其优化过程。将其改写为向量形式
记
o b j = ∑ i = 1 N ⟨ X i β − Y i ⟩ = ( X β − y ) T ( X β − y ) (11) obj = \sum_{i=1}^{N} \langle X_i\beta - Y_i \rangle = (X\beta - y)^T(X\beta - y) \tag{11} obj=i=1∑N⟨Xiβ−Yi⟩=(Xβ−y)T(Xβ−y)(11)
这是一个关于 β \beta β的二次函数，所以这是一个单调不减函数。对其进行一阶和二阶的求导，需要注意的是，使用链式求导法则时，每展开一层，进行一次左乘。
∂ o b j ∂ β = 2 X T ( X β − y ) (12) \frac{\partial obj}{\partial \beta} = 2 X ^{T}(X\beta - y) \tag{12} ∂β∂obj=2XT(Xβ−y)(12)
∂ 2 o b j ∂ β 2 = 2 X T X (13) \frac{\partial^2 obj}{\partial \beta^2} = 2X^TX \tag{13} ∂β2∂2obj=2XTX(13)
此时二阶导大于等于0，和上述结论相匹配。
当 X X X列满秩时， X T X X^TX XTX是正定的。简单证明如下
X
X
X列满秩，则
X
β
=
0
X\beta=0
Xβ=0无非零值，即记将
X
X
X列分块，记为
X
=
[
X
1
C
X
2
C
X
3
C
.
.
.
X
n
C
]
(14)
X = [X_1^C \ X_2^C\ X_3^C ...\ X_n^C] \tag{14}
X=[X1C X2C X3C... XnC](14)
其中，
X
i
C
X_i^C
XiC代表
X
X
X的第
i
i
i个列向量。
又
β
=
[
β
1
,
β
2
.
.
.
β
p
]
T
\beta = [\beta_1, \beta_2...\beta_p]^T
β=[β1,β2...βp]T，有
X
β
=
∑
i
=
1
p
X
i
C
×
β
i
(15)
X \beta = \sum_{i=1}^p X_i^C \times \beta_i \tag{15}
Xβ=i=1∑pXiC×βi(15)
即此时的任何预测向量
X
β
X \beta
Xβ都是
X
X
X列向量
X
C
X^C
XC的线性组合。
X
X
X列满秩，即
X
C
X^C
XC之间线性无关，故不存在一个非零的
β
\beta
β，使得
X
β
=
0
X\beta=0
Xβ=0。
所以
β
T
X
T
X
β
=
(
X
β
)
T
X
β
=
0
\beta^TX^TX\beta=(X \beta)^TX \beta = 0
βTXTXβ=(Xβ)TXβ=0当且仅当
β
=
0
\beta = 0
β=0才成立，所以
X
T
X
X^TX
XTX正定，进一步
X
T
X
X^TX
XTX可逆。
令一阶导数为0，有
β
=
(
X
T
X
)
−
1
X
T
Y
(16)
\beta = (X^TX)^{-1}X^TY \tag{16}
β=(XTX)−1XTY(16)
则预测值 Y ^ = X ( X T X ) − 1 X T Y \hat{Y} = X(X^TX)^{-1}X^TY Y^=X(XTX)−1XTY，其中 X ( X T X ) − 1 X T X(X^TX)^{-1}X^T X(XTX)−1XT被称为hat matrix，因为这个矩阵给Y戴上了一个小帽子~
由于使用的算法不同会导致得到的参数估计变化，所以对于 β \beta β进行进一步明确，运用上标说明其所采用的算法，记为 β l s \beta^{ls} βls。
又
X
X
X的每一行是一个sample，则可将上式写成sample wise的形式。
将
X
T
X^T
XT列分块，将
X
X
X和
Y
Y
Y行分块，则有
β
=
(
X
T
X
)
−
1
X
T
Y
=
(
(
X
1
,
X
2
.
.
.
X
N
)
(
X
1
T
X
2
T
⋮
X
N
T
)
)
−
1
(
(
X
1
,
X
2
.
.
.
X
N
)
(
Y
1
Y
2
⋮
Y
N
)
)
=
(
∑
i
=
1
N
X
i
X
i
T
)
−
1
(
∑
i
=
1
N
X
i
Y
i
)
=
(
N
E
(
X
i
X
i
T
)
)
−
1
(
N
E
(
X
i
Y
i
)
)
=
(
E
(
X
i
X
i
T
)
)
−
1
(
E
(
X
i
Y
i
)
)
(17)
\begin{aligned} \beta &= (X^TX)^{-1}X^TY \\ &= \bigg(\big(X_1, X_2...X_N\big) \left( \begin{matrix} X_1^T \\ X_2^T \\ \vdots \\ X_N^T \end{matrix} \right)\bigg)^{-1} \bigg(\big(X_1, X_2...X_N\big) \left( \begin{matrix} Y_1 \\ Y_2 \\ \vdots \\ Y_N \end{matrix} \right)\bigg) \\ &= \bigg(\sum_{i=1}^N X_iX_i^T\bigg)^{-1} \bigg(\sum_{i=1}^N X_iY_i\bigg)\\ &= \bigg(N\bold{E}(X_iX_i^T)\bigg)^{-1}\bigg(N\bold{E}(X_iY_i)\bigg) \\ & = \bigg(\bold{E}(X_iX_i^T)\bigg)^{-1}\bigg(\bold{E}(X_iY_i)\bigg) \end{aligned} \tag{17}
β=(XTX)−1XTY=((X1,X2...XN)⎝⎜⎜⎜⎛X1TX2T⋮XNT⎠⎟⎟⎟⎞)−1((X1,X2...XN)⎝⎜⎜⎜⎛Y1Y2⋮YN⎠⎟⎟⎟⎞)=(i=1∑NXiXiT)−1(i=1∑NXiYi)=(NE(XiXiT))−1(NE(XiYi))=(E(XiXiT))−1(E(XiYi))(17)
对这样一个最优解，有另一种更有意思的理解角度。
有
∂
o
b
j
∂
β
=
2
X
T
(
X
β
−
y
)
(18)
\frac{\partial obj}{\partial \beta} = 2 X ^{T}(X\beta - y) \tag{18}
∂β∂obj=2XT(Xβ−y)(18)
上述的一阶导数为0，此时是一个0向量。
将导数值按照行进行分解得到如下形式
[ ⟨ X 1 C , X β − y ⟩ ⟨ X 2 C , X β − y ⟩ ) ⋮ ⟨ X N C , X β − y ⟩ ) ] = [ 0 0 ⋮ 0 ] (19) \left[ \begin{matrix} \langle X_1^C, X\beta - y\rangle\\ \langle X_2^C, X\beta - y\rangle)\\ \vdots\\ \langle X_N^C, X\beta - y\rangle) \end{matrix} \right] =\left[ \begin{matrix} 0\\ 0\\ \vdots\\ 0 \end{matrix} \right] \tag{19} ⎣⎢⎢⎢⎡⟨X1C,Xβ−y⟩⟨X2C,Xβ−y⟩)⋮⟨XNC,Xβ−y⟩)⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡00⋮0⎦⎥⎥⎥⎤(19)
也就意味着 X T X^T XT的每一个行向量，即 X X X的每一个列向量点乘residual都为0，进一步有residual垂直于 X X X的列空间，所以此时的最佳预测值 Y ^ \hat{Y} Y^，实际是真实值 Y Y Y在 X X X列空间的投影。
和上文结论相契合，预测值 Y ^ \hat{Y} Y^实际是 X X X列向量的线性组合，自然会落在 X X X列向量，此处为 X 1 X_1 X1和 X 2 X_2 X2所张成的空间里。此时，真实值 Y Y Y在其中的投影可以保证差值 Y − Y ^ Y-\hat{Y} Y−Y^最小。
The Formations of Linear Regression
上文中给出了线性回归的具体形式，但有一个点没有进行讨论，什么是输入？
按照上文的说法，对于一个单独的样本，它的输入是一个向量，记为 X T = { X 1 , X 2 , . . . X p } X^T = \{X_1, X_2, ...X_p\} XT={X1,X2,...Xp}，这里的每一个feature代表的是什么，亦或者说，线性回归中的线性，指的是针对谁的线性关系
一般来说，在获得输入和将数据feed给模型之前，会有一个特征工程(feature engineering) 存在，将处理之后的输入特征feed到模型进行拟合工作，所以实际上，使用的所谓线性模型，是在特征工程之后得到的结果中保持了线性关系，而线性关系简单从形式上来说，就是加权求和。
所以数值输入元素的
- 交叉项
- 指数项
- 数学操作的结果
类别输入元素的
one-hot编码
也可以作为线性回归的输入项。
亦或者说，这种线性关系保持在参数
β
s
\beta s
βs之间。
这就意味着
Y ( X ∣ β ) = β 1 exp ( x 1 ) + β 2 sin ( x 2 ) + β 3 x 1 x 3 3 (20) Y(X|\beta) = \beta_1 \exp(x_1)+\beta_2 \sin(x_2)+\beta_3x_1x_3^3 \tag{20} Y(X∣β)=β1exp(x1)+β2sin(x2)+β3x1x33(20)
也是一个线性回归模型。
每一种transform都被称作一种basis function。
若记特征工程之后得到的feature全体为
H
H
H，则线性回归模型从
X
X
X的角度可写为如下形式：
Y
=
∑
h
j
∈
H
w
j
h
j
(
X
)
(21)
Y = \sum_{h_j \in H} w_jh_j(X) \tag{21}
Y=hj∈H∑wjhj(X)(21)
可以看出，进行线性组合的基底由原先的
X
X
X列向量变为了所有transform的输出，这是一种基底的扩展，使用basis expansion之间的线性代替
X
X
X列向量之间的线性，故上述特征工程也被称为basis expansion in X。
ESL中对于这个操作有着如下的描述
The core idea … is to augment/replace the vector of inputs X with additional variables, which are transformations of X, and then use linear models in this new space of derived input features.
在完成了basis expansion之后，剩余操作和origin linear model没有差异。
若
∃
i
,
h
i
(
X
)
=
X
1
,
X
2
.
.
.
X
p
\exists i, h_i(X)=X_1, X_2...X_p
∃i,hi(X)=X1,X2...Xp
则此时包含了原始的线性回归形式。
若
∃
i
,
h
i
(
X
)
=
X
1
2
,
X
2
2
,
.
.
.
X
1
X
2
,
.
.
.
X
1
X
p
,
.
.
.
X
p
2
\exists i, h_i(X)=X_1^2, X_2^2,...X_1X_2,...X_1X_p,...X_p^2
∃i,hi(X)=X12,X22,...X1X2,...X1Xp,...Xp2
则对模型进行了degree为2的多项式处理。
因此，通过扩充basis expansion的输出
H
H
H，可以增加模型的capacity，加大获得最佳线性模型的概率。但是从模型复杂度的角度来看，需要对模型进行一定的限制。
可以通过
- 先验的特征预选
- subset selection
- 对 w j w_j wj进行正则控制参与拟合的特征数
进行复杂度的控制。
通过上述操作，极大地扩充了模型的capacity，但是这种操作是global的，所以对于整个domain都会进行相同的处理。对于一些domain specific或是分段的函数形式，这种处理方式是不合理的，故进一步提出piecewise polynomial的概念，即将整个domain进行分块，在每一个sub space都分别拟合一个local的polynomial的函数。
Y = ∑ i f i ( x ) I ( x ∈ s u b s p a c e i ) Y = \sum_i f_i(x)I(x \in subspace_i) Y=i∑fi(x)I(x∈subspacei)
通过限制函数取值和指定阶导数在整个区间连续，可以对函数的取值进行一定的限制。
以一个order-K的piecewise polynomial(即
x
x
x的幂次最大为
K
−
1
K-1
K−1，因为存在截距项幂次为0)模型为例，若其从0阶到
K
−
2
K-2
K−2阶的导数都连续，则此piecewise polynomial也被称作order-K的spline。
记以 M i , M 2 , . . . M n M_i,M_2,...M_n Mi,M2,...Mn分割为 n + 1 n+1 n+1个子区间，则对此模型拟合可依赖如下basis：
h
j
(
X
)
=
X
(
j
−
1
)
,
j
=
1...
K
h_j(X)=X^{(j-1)}, j = 1...K
hj(X)=X(j−1),j=1...K
h
K
+
i
(
X
)
(
X
−
M
i
)
+
K
−
1
,
i
=
1
,
2...
n
h_{K+i}(X)(X - M_i)_+^{K-1}, i = 1,2...n
hK+i(X)(X−Mi)+K−1,i=1,2...n
Multiple Outputs
上述的所有的论述都是针对单一输出的回归模型，将模型所对应的空间进行扩展，当输出 Y ∈ R k Y \in R^k Y∈Rk时进行如下讨论
对于包含了 N N N个已进行中心化处理的样本的训练集，此时输入按照行进行组织可得到同样尺度的 X ∈ R N × p X \in R^{N \times p} X∈RN×p，但是输出 Y = [ Y 1 , Y 2 . . . Y k ] Y=[Y_1, Y_2...Y_k] Y=[Y1,Y2...Yk]的形状变化为 R N × k R^{N\times k} RN×k，参数尺度变化为 B ∈ R p × k B \in R^{p \times k} B∈Rp×k。
依然使用mse进行优化，需要对每一个输出维度的损失进行加和，则目标函数可写为
o b j = t r [ ( X B − Y ) ( X B − Y ) T ] (22) obj = tr[(XB - Y)(XB - Y)^T] \tag{22} obj=tr[(XB−Y)(XB−Y)T](22)
引入几个重要的迹相关的求导性质
d
(
t
r
(
X
)
)
=
t
r
(
d
X
)
(23)
d(tr(X)) = tr(dX) \tag{23}
d(tr(X))=tr(dX)(23)
d
(
t
r
(
A
X
T
B
)
)
d
X
=
B
A
(24)
\frac{d(tr(AX^TB))}{dX} = BA \tag{24}
dXd(tr(AXTB))=BA(24)
d
(
t
r
(
A
X
B
)
)
d
X
=
A
T
B
T
(25)
\frac{d(tr(AXB))}{dX} = A^TB^T \tag{25}
dXd(tr(AXB))=ATBT(25)
d
(
t
r
(
A
X
M
X
T
C
)
)
d
X
=
C
A
X
M
+
A
T
C
T
X
M
T
(26)
\frac{d(tr(AXMX^TC))}{dX} = CAXM + A^TC^TXM^T \tag{26}
dXd(tr(AXMXTC))=CAXM+ATCTXMT(26)
故，结合公式26，将
M
=
I
,
X
=
B
,
A
=
X
,
C
=
X
T
M=I,X=B,A=X,C=X^T
M=I,X=B,A=X,C=XT代入，有
∂
[
t
r
(
X
B
B
T
X
T
)
]
∂
B
=
X
T
X
B
+
X
T
X
B
(27)
\frac{\partial[tr(XBB^TX^T)]}{\partial B} = X^TXB+X^TXB \tag{27}
∂B∂[tr(XBBTXT)]=XTXB+XTXB(27)
令导数为0有
∂ o b j ∂ B = ∂ [ t r ( X B B T X T ) − t r ( X B Y T ) − t r ( Y B T X T ) ] ∂ B = 2 X T X B − 2 X T Y = 0 (28) \begin{aligned} \frac{\partial obj}{\partial B}&=\frac{\partial[tr(XBB^TX^T) - tr(XBY^T) -tr(YB^TX^T)]}{\partial B} \\ &=2 X^TXB - 2X^TY = 0 \end{aligned} \tag{28} ∂B∂obj=∂B∂[tr(XBBTXT)−tr(XBYT)−tr(YBTXT)]=2XTXB−2XTY=0(28)
故此时 B = ( X T X ) − 1 X T Y B = (X^TX)^{-1}X^TY B=(XTX)−1XTY，是单输出回归模型的拓展， B B B的每一个列向量互不干扰，可独立优化，记 B = [ B 1 , B 2 . . . B k ] B = [B_1, B_2...B_k] B=[B1,B2...Bk]。有
B i = ( X T X ) − 1 X T Y i (29) B_i = (X^TX)^{-1}X^TY_i \tag{29} Bi=(XTX)−1XTYi(29)
扫一扫
3860
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
