中科院算法研究员带你学算法(1)——线性回归(1)

线性回归（1）

写在最前面

  这几年算法岗位的求职供远大于求，可谓是神仙打架，对于数学基础的考察越来越全面和深入，笔者在独角兽和体制内都工作过，参加过近百场面试，也辅导过数十位同学找到心仪的算法工作，想通过这一系列的教程体现机器学习算法中的数学原理，基本可以涵盖所有算法面试中的难题难点，提升面试时的成功率。
  那我们就首先从线性回归开始我们的学习。

Why We Use It

  线性模型是一种在pre computer age和当前都很有效的模型，形式简单又能保证一定的表达能力，同时提供了可解释性，是函数的一阶泰勒展开。从某种意义上说，简单的形式也保证了较好的抗过拟合能力。当输出连续的预测值时，此种线性模型被称为线性回归模型。

  除了上述优点，对线性回归模型输入进行transform，可以将模型的capacity大大扩展。

The Expression

  对于一个函数$f$，若$f({\sum }_i{\alpha }_i{x}_i)={\sum }_i{\alpha }_{i}f(x_i)$，则称此函数是一个线性函数，进一步可据此得到其表达形式

记输入为
$$X^T=\{X_1,X_2,...X_p\}$$
输出为$Y\in {\mathbf{R}}^N$，则有

$$Y=\sum _{k=1}^P{\beta }_k\times X_k+{\beta }_0$$

其中${\beta }_0$又被称为interception，截距项，可以通过的简单扩展，将原始的输入数据改写为$X^{\prime }=(1;X)$，则原式改写为
$$Y=\sum _{k=0}^P{\beta }_k\times X_k$$
这一步是一个简单的数学操作，进行形式上的统一。

  为了行文方便，下文中都是用$X$代替$X^{\prime }$，在很多场景下，得到的输入$X$都是包含了截距项的。

  给定一个数据集$X$，使用其对模型的参数进行估计，每一行是一个单独的sample。若记每一行的sample为$X_i$，记第$i$个样本的第$j$个特征取值为$X_{i,j}$，则$X$可表示为$X\in {\mathbf{R}}^{N\times (P+1)}$。
记

$$X=\left[\begin{array}{ccccc}{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}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ X_{N,0}& X_{N,1}& X_{N,2}& \cdots & X_{N,p}\end{array}\right]$$

且
$$Y_i=⟨X_i,\beta ⟩=\sum _{j=0}^P{\beta }_j{X}_{ij}$$

则
$$\begin{array}{rl}\hat{Y}& =(Y_1,Y_2...Y_N{)}^T\\ & =(⟨X_1,\beta ⟩,⟨X_2,\beta ⟩,...⟨X_N,\beta ⟩{)}^T\\ & =\left[\begin{array}{c}{\sum }_j{\beta }_j{X}_{1j}\\ {\sum }_j{\beta }_j{X}_{2j}\\ ⋮\\ {\sum }_j{\beta }_j{X}_{Nj}\end{array}\right]\end{array}$$

即预测值
$$\hat{Y}=X\beta \in {\mathbf{R}}^N$$

  所有data driven的任务都对数据有要求，此时假设数据集中的sample都是独立同分布的(i.i.d)，这是为了保证采样的数据集来自同一个分布且互不干扰。

  为了衡量其与真实值的差异大小，最常见的指标会选择RSS(residual sum of squares)，首先计算各个分量的差异，$(Y_i-{\hat{Y}}_i{)}^2$，将每个维度上的误差进行累加，使用此指标作为目标函数的优化算法被称为最小二乘法。

  此外，从方程求解的角度也可以对RSS进行理解和说明。

  若存在$m$个方程，即对于$\forall i\in [1,m],{\sum }_j{a}_{ij}{x}_j=b_i$。从几何的角度来看，这也是定义了$m$个超平面，同时满足约束的点即为交点的集合。

当矩阵
$$A=\left[\begin{array}{c}{a}_{11},a_{12},a_{13}...a_{1n}\\ a_{21},a_{22},a_{23}...a_{2n}\\ ⋮\\ a_{m1},a_{m2},a_{m3}...a_{mn}\end{array}\right]$$

则上述约束可写为
$$Ax=b$$

  当$A$满足$m=n$且满秩时，此矩阵可逆，则交点为$x=A^{-1}b$。

  若$m>n$，则这个情况被称为过约束，很可能不存在交点，则将求取的目标由交点改为到每个超平面的距离之和最小的点$x^{\prime }$，到任意超平面的距离为

  当超平面参数单位化，并使用平方误差代替绝对值之后，得到的优化目标等价于RSS。

  进一步考察其优化过程。将其改写为向量形式

记

$$obj=\sum _{i=1}^N⟨X_i\beta -Y_i⟩=(X\beta -y{)}^T(X\beta -y)$$

  这是一个关于$\beta$的二次函数，所以这是一个单调不减函数。对其进行一阶和二阶的求导，需要注意的是，使用链式求导法则时，每展开一层，进行一次左乘。

  此时二阶导大于等于0，和上述结论相匹配。

  当$X$列满秩时，$X^{T}X$是正定的。简单证明如下

  $X$列满秩，则$X\beta =0$无非零值，即记将$X$列分块，记为
$$X=[X_1^C{X}_2^C{X}_3^C...X_n^C]$$
  其中，$X_i^C$代表$X$的第$i$个列向量。

  又$\beta =[{\beta }_1,{\beta }_2...{\beta }_p{]}^T$，有
$$X\beta =\sum _{i=1}^p{X}_i^C\times {\beta }_i$$
  即此时的任何预测向量$X\beta$都是$X$列向量$X^C$的线性组合。

  $X$列满秩，即$X^C$之间线性无关，故不存在一个非零的$\beta$，使得$X\beta =0$。
  所以${\beta }^T{X}^{T}X\beta =(X\beta {)}^{T}X\beta =0$当且仅当$\beta =0$才成立，所以$X^{T}X$正定，进一步$X^{T}X$可逆。

  令一阶导数为0，有
$$\beta =(X^{T}X{)}^{-1}{X}^{T}Y$$

  则预测值$\hat{Y}=X(X^{T}X{)}^{-1}{X}^{T}Y$，其中$X(X^{T}X{)}^{-1}{X}^T$被称为hat matrix，因为这个矩阵给Y戴上了一个小帽子~

  由于使用的算法不同会导致得到的参数估计变化，所以对于$\beta$进行进一步明确，运用上标说明其所采用的算法，记为${\beta }^{ls}$。

又$X$的每一行是一个sample，则可将上式写成sample wise的形式。

将$X^T$列分块，将$X$和$Y$行分块，则有
$$\begin{array}{rl}\beta & =(X^{T}X{)}^{-1}{X}^{T}Y\\ & =\left(\right(X_1,X_2...X_N\left)\left(\begin{array}{c}{X}_1^T\\ X_2^T\\ ⋮\\ X_N^T\end{array}\right){)}^{-1}\right((X_1,X_2...X_N)\left(\begin{array}{c}{Y}_1\\ Y_2\\ ⋮\\ Y_N\end{array}\right))\\ & =\left(\sum _{i=1}^N{X}_i{X}_i^T{)}^{-1}\right(\sum _{i=1}^N{X}_i{Y}_i)\\ & =(N\mathbf{E}(X_i{X}_i^T){)}^{-1}(N\mathbf{E}(X_i{Y}_i))\\ & =(\mathbf{E}(X_i{X}_i^T){)}^{-1}(\mathbf{E}(X_i{Y}_i))\end{array}$$

对这样一个最优解，有另一种更有意思的理解角度。

有

上述的一阶导数为0，此时是一个0向量。

将导数值按照行进行分解得到如下形式

$$\left[\begin{array}{c}⟨X_1^C,X\beta -y⟩\\ ⟨X_2^C,X\beta -y⟩)\\ ⋮\\ ⟨X_N^C,X\beta -y⟩)\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ ⋮\\ 0\end{array}\right]$$

也就意味着$X^T$的每一个行向量，即$X$的每一个列向量点乘residual都为0，进一步有residual垂直于$X$的列空间，所以此时的最佳预测值$\hat{Y}$，实际是真实值$Y$在$X$列空间的投影。

  和上文结论相契合，预测值$\hat{Y}$实际是$X$列向量的线性组合，自然会落在$X$列向量，此处为$X_1$和$X_2$所张成的空间里。此时，真实值$Y$在其中的投影可以保证差值$Y-\hat{Y}$最小。

The Formations of Linear Regression

  上文中给出了线性回归的具体形式，但有一个点没有进行讨论，什么是输入？

  按照上文的说法，对于一个单独的样本，它的输入是一个向量，记为$X^T=\{X_1,X_2,...X_p\}$，这里的每一个feature代表的是什么，亦或者说，线性回归中的线性，指的是针对谁的线性关系

  一般来说，在获得输入和将数据feed给模型之前，会有一个特征工程(feature engineering) 存在，将处理之后的输入特征feed到模型进行拟合工作，所以实际上，使用的所谓线性模型，是在特征工程之后得到的结果中保持了线性关系，而线性关系简单从形式上来说，就是加权求和。

所以数值输入元素的

• 交叉项
• 指数项
• 数学操作的结果

类别输入元素的

• one-hot编码
  也可以作为线性回归的输入项。

亦或者说，这种线性关系保持在参数$\beta s$之间。
这就意味着

$$Y(X|\beta )={\beta }_1\exp (x_1)+{\beta }_2\sin (x_2)+{\beta }_3{x}_1{x}_3^3$$

也是一个线性回归模型。

每一种transform都被称作一种basis function。

若记特征工程之后得到的feature全体为$H$，则线性回归模型从$X$的角度可写为如下形式：
$$Y=\sum _{h_j\in H}{w}_j{h}_j(X)$$

可以看出，进行线性组合的基底由原先的$X$列向量变为了所有transform的输出，这是一种基底的扩展，使用basis expansion之间的线性代替$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没有差异。

若
$$\exists i,h_i(X)=X_1,X_2...X_p$$
则此时包含了原始的线性回归形式。

若
$$\exists i,h_i(X)=X_1^2,X_2^2,...X_1{X}_2,...X_1{X}_p,...X_p^2$$
则对模型进行了degree为2的多项式处理。

因此，通过扩充basis expansion的输出$H$，可以增加模型的capacity，加大获得最佳线性模型的概率。但是从模型复杂度的角度来看，需要对模型进行一定的限制。

可以通过

• 先验的特征预选
• subset selection
• 对$w_j$进行正则控制参与拟合的特征数

进行复杂度的控制。

  通过上述操作，极大地扩充了模型的capacity，但是这种操作是global的，所以对于整个domain都会进行相同的处理。对于一些domain specific或是分段的函数形式，这种处理方式是不合理的，故进一步提出piecewise polynomial的概念，即将整个domain进行分块，在每一个sub space都分别拟合一个local的polynomial的函数。

$$Y=\sum _i{f}_i(x)I(x\in subspac{e}_i)$$

  通过限制函数取值和指定阶导数在整个区间连续，可以对函数的取值进行一定的限制。

  以一个order-K的piecewise polynomial(即$x$的幂次最大为$K-1$，因为存在截距项幂次为0)模型为例，若其从0阶到$K-2$阶的导数都连续，则此piecewise polynomial也被称作order-K的spline。

记以$M_1,M_2,...M_n$分割为$n+1$个子区间，则对此模型拟合可依赖如下basis function：

$$h_j(X)=X^{(j-1)},j=1...K$$
$$h_{K+i}(X)(X-M_i{)}_{+}^{K-1},i=1,2...n$$

Multiple Outputs

  上述的所有的论述都是针对单一输出的回归模型，将模型所对应的空间进行扩展，当输出$Y\in {\mathbf{R}}^k$时进行如下讨论

  对于包含了$N$个已进行中心化处理的样本的训练集，此时输入按照行进行组织可得到同样尺度的$X\in {\mathbf{R}}^{N\times p}$，但是输出$Y=[Y_1,Y_2...Y_k]$的形状变化为${\mathbf{R}}^{N\times k}$，参数尺度变化为$B\in {\mathbf{R}}^{p\times k}$。

  依然使用mse进行优化，需要对每一个输出维度的损失进行加和，则目标函数可写为

$$obj=tr[(XB-Y)(XB-Y{)}^T]$$

引入几个重要的迹相关的求导性质

故，结合上式，将$M=I,X=B,A=X,C=X^T$代入，有

令导数为0有

故此时$B=(X^{T}X{)}^{-1}{X}^{T}Y$，是单输出回归模型的拓展，$B$的每一个列向量互不干扰，可独立优化，记$B=[B_1,B_2...B_k]$。有

$$B_i=(X^{T}X{)}^{-1}{X}^T{Y}_i$$

保持了和单输出类似的形式。