【机器学习算法推导】线性回归，L1L2正则

线性回归，主要用于从数据中拟合出一条直线（或更高维的平面），这条直线能够很好地体现数据的特征，比如，它能够使得平面上的点都均匀地分布在这条直线上。

1.算法思想

  对于简单线性回归和多元线性回归，其算法过程是相同的，不同之处在于简单线性回归只有一个特征需要拟合，多元线性回归有n个。在线性回归中，我们的数据集拥有一定的参数和相应的输出，令$x^{(i)}$表示第i个数据的参数，$y^{(i)}$表示第i个数据的输出，令目标函数$H_{\theta }(x^{(i)})={\theta }_0+{\theta }_1{x}_1+{\theta }_2{x}_2+...+{\theta }_n{x}_n$，这是我们最终得到的函数，其中x1,x2,…,xn为我们拥有的一个个已知参数，比如我们的线性回归最终要预测房价，那么这里的x可能是房屋大小，离市中心的距离等等条件，函数的最终输出就是预测值。其中的θ都是未知的。不同的θ对应着不同的直线，为了得到最佳的拟合直线，我们的目标就是找到最佳的θ。因此线性回归的任务如下：
min ⁡ θ 1 , θ 2 , . . . θ n { 1 2 m ∑ i = 1 m ( h θ ( x ( i ) ) − y ( i ) ) 2 } \min \limits_{θ_1,θ_2,...θ_n}\{\frac{1}{2m}\sum_{i=1}^{m}(h_θ(x^{(i)})-y^{(i)})^2\} θ1​,θ2​,...θn​min​{2m1​i=1∑m​(hθ​(x(i))−y(i))2}
  其中$(h_{\theta }(x^{(i)})-y^{(i)}{)}^2$表示函数的预测值与实际值的平方误差(square error)。由于我们有若干的数据，因此对他们的平方误差计算后累加，由于相加后误差会变大，因此我们取一个均值除以m，最后再除以2的原因是，方便后面的计算（2可以被约掉），并且多除一个2并不会影响整个式子的意义。

2.损失函数

  定义$J({\theta }_0,{\theta }_1,...,{\theta }_n)=\frac{1}{2m}{\sum }_{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)}{)}^2$，其中的$(h_{\theta }(x^{(i)})-y^{(i)}{)}^2$作为损失函数，它定义了预测值与实际值之间的一种联系，当损失函数的值越小时，说明预测值与实际值越接近，从而我们的目标变为找到合适的θ使得损失函数的值越小越好。

3.梯度下降

  梯度下降(Gradient Descent)算法是机器学习中一种常用的算法，可以用来最小化我们的损失函数。损失函数并不止上面用到的这种，还包括其他各种各样的损失函数，都可以利用梯度下降来找到一个局部最优解。
  梯度下降算法的思想是这样的：对于一个由若干个θ形成的直线、曲线或者平面或其他形状，我们随机地给θ赋初始值，它必定位于图像上的某个位置，在那个位置上，寻找一个最陡的方向往下走，重复这个步骤，最终肯定能够到达某个点，并且没有比他更低的点，此时我们就找到了最优解。
  所谓的“最陡”怎么定义的呢？我们知道一条曲线，最陡的地方就是它的切线方向，因此我们可以通过求导的方式来得到它，对于多元函数，我们需要用到偏导来求。因此得到参数的更新公式如下：
$${\theta }_j:={\theta }_j-\alpha \frac{\partial }{\partial \theta }J({\theta }_0,{\theta }_1,...,{\theta }_n)$$

4.学习率

  其中α定义为学习率，直观上看就是每次更新时的步长。当α较小时，每次更新变化的数值都比较小，因此可能需要多次迭代才能到达最优解。当α较大时，每次更新变化的数值较大，可能会跳过最低点，甚至导致无法收敛的情况。通常a从0.001,0.01,0.1,1这些值里面尝试。

5.计算公式

  在简单线性回归的情况下，有$h_{\theta }(x^{(i)})={\theta }_0+{\theta }_{1}x1$。
  当j=0时，$$\frac{\partial }{\partial {\theta }_0}J({\theta }_0,{\theta }_1)=\frac{1}{m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)})$$
  当j=1时，$$\frac{\partial }{\partial {\theta }_1}J({\theta }_0,{\theta }_1)=\frac{1}{m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)})\cdot x^{(i)}$$
  此时${\theta }_0$和${\theta }_1$的更新公式如下：
$${\theta }_0={\theta }_0-\alpha \cdot \frac{1}{m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)})$$
$${\theta }_1={\theta }_1-\alpha \cdot \frac{1}{m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)})\cdot x^{(i)}$$
  在多元线性回归的情况下，计算方式也是一样的。

6.局部最优与全局最优

  由于上面的梯度下降算法是先累计m次的结果后再进行参数的更新，因此也称为批量梯度下降(Batch Gradient Descent)。它有一个明显的缺点，就是容易陷入局部最优解的情况，当我们的参数构成的图像并不完全是凸的时候，不同的初始值可能会陷入到不同的局部最优解中。当梯度下降算法接近一个最优解时，导数值会越来越小，最终收敛，但是这个局部最优解很可能不是一个全局最优解，虽然对于简单线性回归来说并没有这个问题，因为它的图像就是一个简单的凸平面，局部最优就是全局最优。这里涉及到凸优化理论，不做深入讨论，待之后再进行学习吧。

7.特征缩放（feature scaling）

  我们的参数，通常有各种各样的范围，比如身高的范围在150到200。体重的范围在80到200，年龄到范围在0到150。当这些参数作用到损失函数上时，会导致我们的特征等高线图变得臃肿（实在不知道怎么形容）…在这样的图形上，局部下降需要迭代更多的次数，才能够到达最优解所在的位置。
  为了防止这种情况，使用特征缩放将参数归一化到一个固定的范围，如-1到1之间。当所有特征都位于这个范围内时，特征等高线图会变得平滑，也更容易收敛。

8.特征归一化

  通过梯度下降法求解的模型通常是需要归一化的，包括线性回归，逻辑回归，SVM，神经网络等模型，它使得特征之间的更新速度更为一致，容易更快地通过梯度下降找到最优解。

8.1 线性函数归一化

  一种常见的特征缩放技巧，将特征$x_i$缩放为$$x_{norm}=\frac{x_i-x_{min}}{x_{max}-x_{min}}$$，从而使$0\le x_{norm}\le 1$

8.2 零均值归一化

  另一种归一化的方式，它会将原始数据映射到均值为0，标准差为1的分布上。假设原始特征的均值为$\mu$，标准差为$\sigma$，那么归一化公式为将$$z=\frac{x-\mu }{\sigma }$$

9.标准方程法（正规方程法）

  使用梯度下降来求解线性回归，需要较多次的迭代，还需要选择合适的a，才能够收敛到一个合适的值。使用标准方程法，可以一次性求解θ的最优值，下面看看标准方程法的具体过程。

10.标准方程法推导

10.1 从损失函数角度

  令线性回归的解析式如下
$$f(x)=W^{T}x$$
  其中W，x为列向量。
  线性回归的损失函数如下
$$\begin{array}{rl}L(w)& =\sum _{i=1}^n||w^T{x}_i-y_i|{|}^2\\ & =(W^T{X}^T-Y^T)(XW-Y)\\ & =W^T{X}^{T}XW-W^T{X}^{T}Y-Y^{T}XW+Y^{T}Y\\ & =W^T{X}^{T}XW-2W^T{X}^{T}Y+Y^{T}Y\end{array}$$
  则W的最优解$\hat{W}$为
$$\hat{W}=\arg \min L(W)$$
  通过求导解其最优解
$$\begin{array}{r}\frac{\partial L(W)}{\partial W}=2X^{T}XW-2X^{T}Y\end{array}$$
  令偏导为0，则
$$2X^{T}XW=2X^{T}Y⟺W=(X^{T}X{)}^{-1}{X}^{T}Y$$
  得到最终的解析式。

10.2 从线性代数角度

  对上面的解析式稍作变形，使其变为以下形式
$$y=x^{T}W$$
  变形后，矩阵X中，每一列代表原始样本中的一条数据，转置之后就变成了每一行代表原先的一条数据。假如原先每个数据为p*1的向量，总共有n个数据，则此时每一列为所有样本的特征x1，x2，x3…,xp,如果将每一列看做是样本点，那么我们就有了p个n维向量，矩阵Y也是已知的，我们需要求解$X^{T}W=Y$，找到合适的W使得$X^{T}W=Y$成立。根据线性代数的知识我们知道$X^{T}W=Y$只有当$Y$中每个值是对应$X_i$每个值的线性组合，此方程才有解。因为W矩阵反映了他们之间是怎样的一种线性组合。换句话说，Y向量必须存在于X矩阵列向量形成的列空间中，方程才能有解。但是这个不一定成立，所以我们想办法，将向量Y映射到X的列空间中。
  从几何的角度上来看，实际上就是我们要在列空间中寻找一个向量，使Y向量离这条线是最近的，或者说向量Y离这个空间最近，那我们要找的也就是向量Y在这个空间内的投影P了。这个投影P实际上就是由列向量的线性组合得到，即$$P=X^{T}W$$，由于Y-P可以得到法向量，并且法向量与平面内的向量垂直，则有
$$X^T(Y-P)=X^T(Y-X^{T}W)=0$$
  解得
$$W=(X^{T}X{)}^{-1}{X}^{T}Y$$
  当$X^{T}X$可逆时，便可以得到$W=(X^{T}X{)}^{-1}{X}^{T}Y$。这就是标准方程法最终的形式。当$X^{T}X$不可逆时，可以使用伪逆来计算，虽然这种情况很少出现。如果出现了这种情况，说明维度高，样本少，会出现过拟合的现象。

11.L2正则

  正则项的出现，是为了抑制过拟合的现象。L2正则表示为$$||w|{|}_2^2=w_1^2+w_2^2+...+w_n^2=w^{T}w$$
  在模型后面加上L2正则作为惩罚项，可以使得总体的w偏小，不至于过拟合。加入正则项后的损失函数如下
$$\begin{array}{rl}J(W)& =\sum _{i=1}^n||w^T{x}_i-y_i|{|}^2+\lambda W^{T}W\\ & =W^T{X}^{T}XW-2W^T{X}^{T}Y+Y^{T}Y+\lambda W^{T}W\\ & =W^T(X^{T}X-\lambda I)W-2W^T{X}^{T}Y+Y^{T}Y\end{array}$$
  则有
$$\hat{w}=\arg \underset{W}{\min }J(W)$$
  求导并令导数为0，可得
$$\begin{gathered} \frac{\partial J(W)}{\partial W}=2(X^{T}X+\lambda I)W-2X^{T}Y=0 \\ \hat{W}=(X^{T}X+\lambda I{)}^{-1}{X}^{T}Y \end{gathered}$$
  由于$X^{T}X$为半正定矩阵（由定义可证），加上$\lambda$倍对角矩阵之后，一定是正定矩阵，而正定矩阵一定有逆矩阵，所以它起到了抑制过拟合的效果。

12.概率视角看线性回归

  我们将样本x看做是独立同分布的随机变量的不同取值，y是输出，定义一个高斯噪声$ϵ\sim N(0,{\sigma }^2)$，则
$$y=f(x)+ϵ=w^{T}x+ϵ$$
  则有
$$y|x,w\sim N(w^{T}x,{\sigma }^2)$$
  则其似然函数为
$$\begin{gathered} \mathcal{L}(w)=\log P(Y|X,w)=\log \prod _{i=1}^{n}P(y_i|x_i;w)=\sum _{i=1}^n\log P(y_i|x_i;w) \\ =\sum _{i=1}^n\left(\log \frac{1}{\sqrt{2\pi }\sigma }+\log \exp (-\frac{(y_i-w^T{x}_i{)}^2}{2{\sigma }^2})\right) \\ =\sum _{i=1}^n\left(\log \frac{1}{\sqrt{2\pi }\sigma }-\frac{(y_i-w^T{x}_i{)}^2}{2{\sigma }^2}\right) \end{gathered}$$
  则
$$\begin{array}{rl}\hat{w}& =\arg \underset{w}{\max }\mathcal{L}\\ & =\arg \underset{w}{\max }\left(-\frac{(y_i-w^T{x}_i{)}^2}{2{\sigma }^2}\right)\\ & =\arg \underset{w}{\min }{\left(y_i-w^T{x}_i\right)}^2\end{array}$$
  可以发现上式实际上等价于最小二乘估计的定义一模一样。所以最小二乘估计实际上隐含了一个它的噪声服从正态分布这样的一个假设！所以在线性回归中，有$$LSE=MLE+高斯噪声$$
  在添加正则化项之后，我们给W也来一个先验分布$W\sim N(0,{\sigma }_0^2)$，由贝叶斯公式
$$p(w|y)=\frac{p(y|w)p(w)}{p(y)}$$
  有
$$\begin{array}{rl}\hat{W}& =\arg \underset{W}{\max }\log P(y|w)P(w)\\ & =\arg \underset{W}{\max }\log (\frac{1}{\sqrt{2\pi }\sigma }\frac{1}{\sqrt{2\pi }{\sigma }_0})+(-\frac{(y-w^{T}x{)}^2}{2{\sigma }^2}-\frac{||w|{|}^2}{2{\sigma }_0^2})\\ & =\arg \underset{W}{\min }(\frac{(y-w^{T}x{)}^2}{2{\sigma }^2}+\frac{||w|{|}^2}{2{\sigma }_0^2})\\ & =\arg \underset{W}{\min }\left((y-w^{T}x{)}^2+\frac{{\sigma }^2}{{\sigma }_0^2}||w|{|}^2\right)\end{array}$$
  把$\frac{{\sigma }^2}{{\sigma }_0^2}$看做是$\lambda$的话，我们发现这个式子跟加了L2正则的最小二乘估计的模型是一样的。所以我们又得到一个结论
$$LSE(regularized)⟺MAP(高斯噪声+参数w服从高斯分布)$$

13.L2正则分析

  进一步分析，令$w^{\ast }$为未正则化的目标函数取得最小训练误差时的权重向量，即$w^{\ast }=\arg {\min }_{w}J(w)$，在$w^{\ast }$领域附近对目标函数做泰勒展开，有
$$\hat{J}(\theta )=J(w^{\ast })+\frac{1}{2}(w-w^{\ast }{)}^{T}H(w-w^{\ast })$$
  因为$w^{\ast }$为最优，因此一阶项为0。对上式求导，有
$${\nabla }_w\hat{J}(w)=H(w-w^{\ast })=0$$
  用$\tilde{w}$表示添加正则项后的最优解，则上式变为
$$a\tilde{w}+H(\tilde{w}-w^{\ast })=0⟺\tilde{w}=(H+aI{)}^{-1}H{w}^{\ast }$$
  由于H为实对称矩阵，对其进行分解有
$$H=Q\Lambda Q^T$$
  则有
$$\begin{array}{rl}\tilde{w}& =(Q\Lambda Q^T+aI{)}^{-1}Q\Lambda Q^T{w}^{\ast }\\ & =(Q(\Lambda +aI)Q^T{)}^{-1}Q\Lambda Q^T{w}^{\ast }\\ & =Q(\Lambda +aI{)}^{-1}\Lambda Q^T{w}^{\ast }\end{array}$$
  可以发现权重衰减的效果其实是沿着由H的特征向量所定义的轴缩放$w^{\ast }$，缩放因子为$\frac{{\lambda }_i}{{\lambda }_i+a}$。对于远大于a的特征值，对其应用缩放因子之后，影响并不会很大，而对于远小于a的特征值，将会收缩到几乎为零。只有在显著减小目标函数方向上的参数会保留得相对完好。在无助于目标函数减小的方向（对应 Hessian 矩阵较小的特征值）上改变参数不会显著增加梯度。这种不重要方向对应的分量会在训练过程中因正则化而衰减掉。

14.L1正则分析

  对模型的L1正则定义如下
$$||w|{|}_1=\sum _i|w_i|$$
  正则化的目标函数如下
$$\tilde{J}(w)=J(w)+a||w|{|}_1$$
  对应的梯度为
$${\nabla }_w\tilde{J}(w)=asign(w)+{\nabla }_{m}J(w)$$
  在完全一般化的海塞矩阵的情况下，无法直接得到清晰的代数表达式，因此我们进一步简化Hessian矩阵为对角矩阵（去除输入特征之间的相关性，那么这一假设成立）。那么我们可以将L1正则化目标函数的二次泰勒展开近似分解成关于参数的求和:
$$\tilde{J}(w)=J(w^{\ast })+\sum _i\left[\frac{1}{2}{H}_{ii}(w_i-w_i^{\ast }{)}^2+a|w_i|\right]$$
  下面的解析解可以最小化这个近似代价函数
w i = s i g n ( w i ∗ ) max ⁡ { ∣ w i ∗ ∣ − a H i i } w_i=sign(w_i^*)\max\{|w_i^*|-\frac{a}{H_{ii}}\} wi​=sign(wi∗​)max{∣wi∗​∣−Hii​a​}
  当$w_i^{\ast }$大于0时，将最优值移动$\frac{a}{H_{ii}}$往0逼近，若是移动后变号了就直接赋值为0，小于0时也同样。
  相比L2正则化，L1正则化会产生更稀疏的解。由L1正则化导出的稀疏性质已经被广泛地用于特征选择。

15.标准方程法vs梯度下降法

  标准方程法虽然可以一次性求出最优解，也不需要选择合适的学习率，但是在求$(X^{T}X{)}^{-1}$时需要耗费大量时间，当我们的特征数较多时，往往不推荐使用。通常当特征数小于10000时较为合适。
  梯度下降法虽然需要迭代，需要选取合适的学习率，但是它的适用范围较广，在其他众多机器学习算法中，梯度下降算法往往是适用的，虽然迭代次数上可能有些差异，但是总体上表现良好。而标准方程法只适用于求线性回归的情况。