机器学习（回归一）——线性回归-最小二乘

从这篇博客开始将介绍机器学习中的比较重要的内容——回归，期间会穿插着涉及到的知识点，如正则化。

什么是回归（及线性回归）

概念

说到回归，一般都是指线性回归（linear regression），但从广义上来说线性回归是回归家族中的一种，也是最简单的一种。简单来说， 假设现在有一些数据点，我们用一条直线对这些点进行拟合（该线称为最佳拟合直线），这个拟合过程就称作回归。

回归算法是一种比较常用的机器学习算法，属于有监督算法，用来建立“解释”变量(自变量X)和观测值(因变量Y)之间的关系；从机器学习的角度来讲，用于构建一个算法模型(函数)来做属性(X)与标签(Y)之间的映射关系，在算法的学习过程中，试图寻找一个函数 $h:R^2\to R$ 使得参数之间的关系拟合性最好。

回归算法中算法(函数)的最终结果是一个连续的数据值，输入值(属性值)是一个d维度的属性/数值向量。

“回归”一词的来历 :
今天我们所知道的回归是由达尔文（Charles Darwin）的表兄弟Francis Galton发明的。Galton于1877年完成了第一次回归预测，目的是根据上一代豌豆种子（双亲）的尺寸来预测下一代豌豆种子（孩子）的尺寸。Galton在大量对象上应用了回归分析，甚至包括人的身高。他注意到，如果双亲的高度比平均高度高，他们的子女也倾向于比平均高度高，但尚不及双亲。孩子的高度向着平均高度回退（回归）。Galton在多项研究上都注意到这个现象，所以尽管这个英文单词跟数值预测没有任何关系，但这种研究方法仍被称作回归 。

举例

举个房价预测的例子，现在有房屋面积及其对应的租赁价格如下

房屋面积($m^2$)	租赁价格(1000￥)
10	0.8
15	1
20	1.8
30	2
50	3.2
60	3
60	3.1
70	3.5

请问，如果现在有一个房屋面积为55平，最终的租赁价格是多少比较合适?
比如近似满足y = ax + b。求得a b，代入55，就可得到大概价格。大致的函数图像如下：

当不仅仅考虑面积，还要考虑房间数数量。可在三维坐标系中表示。即：$h(x)={\theta }_0+{\theta }_1{x}_1+{\theta }_2{x}_2$

房屋面积($m^2$)	房间数量	租赁价格(1000￥)
10	1	0.8
20	1	1.8
30	1	2.2
30	2	2.5
70	3	5.5
70	2	5.2
…	…	…

如果考虑的因素更多，也就是有更多的维度呢？此时的表达式：
$$\begin{array}{rl}{h}_{\theta }(x)& ={\theta }_0+{\theta }_1{x}_1+\dots +{\theta }_n{x}_n\\ & ={\theta }_{0}1+{\theta }_1{x}_1+\dots +{\theta }_n{x}_n\\ & ={\theta }_0{x}_0+{\theta }_1{x}_1+\dots +{\theta }_n{x}_n\\ & =\sum _{i=0}^n{\theta }_i{x}_i\\ & ={\mathbf{\theta }}^{T}x\end{array}$$

最终要求是计算出$\theta$ 的值，并选择最优的$\theta$值构成算法公式。其中$x_0=1$，${\theta }_0$为偏置项，也就是截距。
ps：（线性回归的本质）

1. 线性：θ最高次项为一次（$\mathbf{\theta }$是我们最终要求得的值）
2. 回归：Y（或者 $h_{\theta }(x)$）连续

最小二乘法公式推导

误差

上面例子中找到的最终表达式就一定是正确的吗？不一定，往往存在误差，即每个样本的误差。就像上图中，那个平面是我们找到表达式，每个红点是真实的值，他们之间就可能会存在误差，记为${\epsilon }_n$。

面积为10平和面积为70平之间有关系吗？没有。所以产生的误差${\epsilon }_1$和${\epsilon }_2$有关系吗？没有，也就是独立的。（独立同分布）$\epsilon$是在各个特征上的误差的叠加。$\mathbf{\epsilon }=({\epsilon }_1,{\epsilon }_2,\dots ,{\epsilon }_n$）根据中心极限定理，服从高斯正态分布，而且是均值为0。
$$y^{(i)}={\mathbf{\theta }}^T{x}^{(i)}+{\epsilon }^{(i)}\text{\hspace{1em}(1)}$$

对于每个样本 $x^{(i)}$ 代入 h(x) 得到的估计值 ${\hat{y}}^{(i)}$ 与真实值 $y^{(i)}$ 存在的差值 ${\epsilon }^{(i)}$

总结一下，就是误差 ${\epsilon }^{(i)}(1\le i\le n)$ 是独立同分布的，服从均值为0，方差为某定值 ${\sigma }^2$ 的高斯正态分布。（原因：中心极限定理）。实际问题中，很多随机现象可以看做众多因素的独立影响的综合反应，往往服从正态分布.

观察上面的（1）式，$y^{(i)}$ 就是真实值，并且 ${\epsilon }^{(i)}$ 满足均值为0的高斯正态分布，所以：
$$p({\epsilon }^{(i)})=\frac{1}{\sigma \sqrt{2\pi }}{e}^{-\frac{{\left({\epsilon }^{(i)}\right)}^2}{2{\sigma }^2}}\text{\hspace{1em}(2)}$$
把（1）式代入（2）式，得：
$$p(y^{(i)}|x^{(i)};\mathbf{\theta })=\frac{1}{\sigma \sqrt{2\pi }}exp\left(-\frac{{\left(y^{(i)}-{\mathbf{\theta }}^T{x}^{(i)}\right)}^2}{2{\sigma }^2}\right)\text{\hspace{1em}(3)}$$
注：$e^x$ 与 $exp(x)$ 是一样的。

似然函数

构建似然函数。为什么要这样构建，可以参考《【最透彻】13_最大似然估计【小元老师】考研数学，概率论与数理统计》，小元老师讲得确实很透彻。
$$\begin{array}{rl}L(\theta )& =\prod _{i=1}^{m}p(y^{(i)}|x^{(i)};\mathbf{\theta })\\ & =\prod _{i=1}^m\frac{1}{\sigma \sqrt{2\pi }}exp\left(-\frac{{\left(y^{(i)}-{\mathbf{\theta }}^T{x}^{(i)}\right)}^2}{2{\sigma }^2}\right)\end{array}\text{\hspace{1em}(4)}$$

对数似然

初始目的是让似然函数 $L(\theta )$ 最大化，但不方便计算。现在把（4）式取对数（因为取对数不改变单调性，即$L(\theta )$ 和 $l(\theta )$单调性一致，并且能把上面的“累积”变成“累和”，从而简化运算）
$$\begin{array}{rl}l(\theta )& =\log L(\theta )\\ & =\log \prod _{i=1}^m\frac{1}{\sigma \sqrt{2\pi }}exp\left(-\frac{{\left(y^{(i)}-{\mathbf{\theta }}^T{x}^{(i)}\right)}^2}{2{\sigma }^2}\right)\\ & =m\log \frac{1}{\sigma \sqrt{2\pi }}-\frac{1}{{\sigma }^2}\bullet \frac{1}{2}\sum _{i=1}^m{\left(y^{(i)}-{\mathbf{\theta }}^T{x}^{(i)}\right)}^2\end{array}\text{\hspace{1em}(5)}$$

目标函数

通过上面取对数，把似然函数 $L(\theta )$ 最大化问题转换成$l(\theta )$最大化问题。

因为（5）式中 $m\log \frac{1}{\sigma \sqrt{2\pi }}$ 与 $\frac{1}{{\sigma }^2}$ 是定值，所以 $l(\theta )$的大小取决于 $\frac{1}{2}{\sum }_{i=1}^m{\left(y^{(i)}-{\mathbf{\theta }}^T{x}^{(i)}\right)}^2$。当 $\frac{1}{2}{\sum }_{i=1}^m{\left(y^{(i)}-{\mathbf{\theta }}^T{x}^{(i)}\right)}^2$ 最小时 $l(\theta )$ 取得最大值。所以，似然函数 $L(\theta )$ 最大化等价于 $\frac{1}{2}{\sum }_{i=1}^m{\left(y^{(i)}-{\mathbf{\theta }}^T{x}^{(i)}\right)}^2$ 最小化。即得到最终的目标函数:
$$loss(y_j,{\hat{y}}_j)=J(\theta )=\frac{1}{2}\sum _{i=1}^m{\left(h_{\theta }\left(x^{(i)}\right)-y^{(i)}\right)}^2\text{\hspace{1em}(6)}$$
我们的最终目的，就是让目标函数$J(\theta )$最小化。（这里给出$\frac{1}{2}$的目的就是为了方便求导时化简用）

求解 $\theta$

我们的基本思路就是，因为 $J(\theta )$ 是凸函数，在导数为0时，取得极值，所以重点就是求导。
$$J(\theta )=\frac{1}{2}\sum _{i=1}^m{\left(h_{\theta }\left(x^{(i)}\right)-y^{(i)}\right)}^2=\frac{1}{2}(X\mathbf{\theta }-\mathbf{y}{)}^T(X\mathbf{\theta }-\mathbf{y})\text{\hspace{1em}(7)}$$
对（7）式（（7）式的详细推导在文章最后有介绍） $J(\theta )$ 求导（用到了矩阵求导，详细内容文章最后有介绍）：
$$\begin{array}{rl}{\nabla }_{\theta }J(\theta )& ={\nabla }_{\theta }\left(\frac{1}{2}(X\mathbf{\theta }-\mathbf{y}{)}^T(X\mathbf{\theta }-\mathbf{y})\right)\\ & ={\nabla }_{\theta }\left(\frac{1}{2}({\mathbf{\theta }}^T{X}^T-{\mathbf{y}}^T)(X\mathbf{\theta }-\mathbf{y})\right)\\ & ={\nabla }_{\theta }\left(\frac{1}{2}({\mathbf{\theta }}^T{X}^{T}X\mathbf{\theta }-{\mathbf{\theta }}^T{X}^T\mathbf{y}-{\mathbf{y}}^{T}X\mathbf{\theta }+{\mathbf{y}}^T\mathbf{y})\right)\\ & =\frac{1}{2}\left(2X^{T}X\mathbf{\theta }-X^T\mathbf{y}-X^T\mathbf{y}\right)\\ & =X^{T}X\mathbf{\theta }-X^T\mathbf{y}\end{array}\text{\hspace{1em}(8)}$$
导求为0，即 $X^{T}X\mathbf{\theta }-X^T\mathbf{y}=0$，求得最优解
$$\begin{gathered} X^{T}X\mathbf{\theta }-X^T\mathbf{y}=0 \\ \Downarrow \\ {X}^{T}X\mathbf{\theta }=X^T\mathbf{y} \\ \Downarrow \\ {\left(X^{T}X\right)}^{-1}{X}^{T}X\mathbf{\theta }={\left(X^{T}X\right)}^{-1}{X}^T\mathbf{y} \\ \Downarrow \\ \mathbf{\theta }={\left(X^{T}X\right)}^{-1}{X}^T\mathbf{y} \end{gathered}$$
所以，最小二乘法（最小二乘法的定义文章最后有介绍）的参数最优解为：
$$\underset{\mathbf{\theta }}{\min }J(\theta )={\left(X^{T}X\right)}^{-1}{X}^T\mathbf{y}\text{\hspace{1em}(9)}$$

参数解析式优化

最小二乘法的使用要求矩阵 $X^{T}X$ 是可逆的；为了防止不可逆或者过拟合的问题存在，可以增加额外数据影响，使得最终的矩阵是可逆的：
$$\mathbf{\theta }={\left(X^{T}X+\lambda I\right)}^{-1}{X}^T\mathbf{y}\text{\hspace{1em}(10)}$$
（10）式（文章最后有推导，来帮助理解）其实是引入L2正则化的岭回归（后续博客会有专门介绍）的解。
由于矩阵的逆的求解是一个难处，尤其对数据量特别大的时候。所以实际中很少用最小二乘法。

总结

上面讲的线性回归，直观来讲就是找出自变量X和因变量Y的线性关系，也就我们要训练出来的就是这样一个模型：
$$h_{\theta }(x)={\theta }_0+{\theta }_1{x}_1+\dots +{\theta }_n{x}_n$$
实际上就是训练出 $\mathbf{\theta }$ 值。通过模型，就能得到预测值，我们最终要做的就是预测值和真实值之间的差最小，可通过误差来反映，得到一个损失函数：
$$J(\theta )=\frac{1}{m}\sum _{i=1}^m\left(y^{(i)}-{\hat{y}}^{(i)}\right)$$
上文中是通过最大似然估计得到的目标函数：
$$J(\theta )=\frac{1}{2}\sum _{i=1}^m{\left(h_{\theta }\left(x^{(i)}\right)-y^{(i)}\right)}^2$$
效果是一样的，这也是最小二乘法。
由于目标函数是一个凸函数的形式，那么在极值处取得最小值，也就是导数为0的时候。最终得到：
$$\underset{\mathbf{\theta }}{\min }J(\theta )={\left(X^{T}X\right)}^{-1}{X}^T\mathbf{y}$$
但有个要求，就是 $X^{T}X$ 必须可逆，此时保证不了，但正定矩阵一定可逆。而 $X^{T}X$ 是一个半正定矩阵，可以给他的对角线加一个偏移量，得到一个正定矩阵，即：
$$\mathbf{\theta }={\left(X^{T}X+\lambda I\right)}^{-1}{X}^T\mathbf{y}$$
这个式子，也是加上L2正则，即岭回归求解得到的式子，因此这个式子也能解决过拟合的问题。

---

---

补充

式子（7）的推导

对于列向量 $\mathbf{\alpha }=(a_1,a_2,a_3,\dots ,a_m{)}^T$ :
$$\begin{array}{r}{\mathbf{\alpha }}^T\mathbf{\alpha }=\mathbf{\alpha }\bullet \mathbf{\alpha }=a_1^2+a_2^2+a_3^2+\dots +a_m^2=\sum _{i=1}^m{a}_i^2\end{array}$$
$\mathbf{\alpha }\bullet \mathbf{\alpha }$为向量点乘，也称为数量积或内积。

注意，我们在机器学习公式推导中用到的数据表示：

$x^{(i)}$：表示第 i 个样本；
$x_i$：x 向量的第 i 维度的值。

对于真实的值构成的列向量 $\mathbf{y}=(y^{(1)},y^{(2)},y^{(3)},\dots ,y^{(m)}{)}^T$
以及对应的预测值构成的列向量 $\hat{\mathbf{y}}=({\hat{y}}^{(1)},{\hat{y}}^{(2)},{\hat{3}}^{(1)},\dots ,{\hat{y}}^{(m)}{)}^T$
求得的参数解列向量 $\mathbf{\theta }=({\theta }_0,{\theta }_1,{\theta }_2,\dots ,{\theta }_n{)}^T$
而样本X（m条样本，样本的维度为n）为：
$$X=\left(\begin{array}{ccccc}{x}_0^1& x_1^1& x_2^1& \cdots & x_n^1\\ x_0^2& x_1^2& x_2^2& \cdots & x_n^2\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ x_0^m& x_1^m& x_2^m& \cdots & x_n^m\end{array}\right)$$
其中$\theta$是未知数，可以写出关于$\theta$的方程组:
$$\{\begin{array}{llc}{x}_0^{(1)}{\theta }_0+x_1^{(1)}{\theta }_1+x_2^{(1)}{\theta }_2+& \cdots & +x_n^{(1)}{\theta }_n={\hat{y}}^{(1)}\\ x_0^{(2)}{\theta }_0+x_1^{(2)}{\theta }_1+x_2^{(2)}{\theta }_2+& \cdots & +x_n^{(2)}{\theta }_n={\hat{y}}^{(2)}\\ & \cdots & \\ x_0^{(m)}{\theta }_0+x_1^{(m)}{\theta }_1+x_2^{(m)}{\theta }_2+& \cdots & +x_n^{(m)}{\theta }_n={\hat{y}}^{(m)}\end{array}$$
也就是：$$X\mathbf{\theta }=\hat{\mathbf{y}}$$
并且：
$$h_{\theta }\left(x^{(i)}\right)=x_0^{(i)}{\theta }_0+x_1^{(i)}{\theta }_1+x_2^{(i)}{\theta }_2+\cdots +x_n^{(i)}{\theta }_n=x^{(i)}\theta ={\hat{y}}^{(i)}$$
所以：
$$\begin{array}{rl}J(\theta )& =\frac{1}{2}\sum _{i=1}^m{\left(h_{\theta }\left(x^{(i)}\right)-y^{(i)}\right)}^2\\ & =\frac{1}{2}\sum _{i=1}^m{\left({\hat{y}}^{(i)}-y^{(i)}\right)}^2\\ & =\frac{1}{2}\left({\left({\hat{y}}^{(1)}-y^{(1)}\right)}^2+{\left({\hat{y}}^{(1)}-y^{(1)}\right)}^2+\cdots +{\left({\hat{y}}^{(n)}-y^{(n)}\right)}^2\right)\\ & =\frac{1}{2}{\left(\hat{\mathbf{y}}-\mathbf{y}\right)}^2\\ & =\frac{1}{2}\left(\hat{\mathbf{y}}-\mathbf{y}\right)\bullet \left(\hat{\mathbf{y}}-\mathbf{y}\right)\\ & =\frac{1}{2}{\left(\hat{\mathbf{y}}-\mathbf{y}\right)}^T\left(\hat{\mathbf{y}}-\mathbf{y}\right)\\ & =\frac{1}{2}(X\mathbf{\theta }-\mathbf{y}{)}^T(X\mathbf{\theta }-\mathbf{y})\end{array}$$

矩阵求导

式子（8）中用到了矩阵求导，常见的几个求导公式：

详细内容可参考《矩阵求导法则与性质》

对于上面的第2个式子，$f(x)={\mathbf{x}}^{T}A\mathbf{x}$，当 $A$ 为对称矩阵，即 $A=A^T$时：
$$\frac{\partial f(x)}{\partial \mathbf{x}}=\frac{\partial \left({\mathbf{x}}^{T}A\mathbf{x}\right)}{\partial \mathbf{x}}=A\mathbf{x}+A^T\mathbf{x}=2A\mathbf{x}$$
在式子（8）中，$X^{T}X$ 为对称矩阵，所以：
$$\frac{\partial ({\mathbf{\theta }}^T{X}^{T}X\mathbf{\theta })}{\partial \mathbf{\theta }}=2X^{T}X\mathbf{\theta }$$

对于式于（8）中的 ${\mathbf{\theta }}^T{X}^T\mathbf{y}$ ，可以参考上面第4个式子中的第一个，即：
$$\frac{\partial ({\mathbf{\theta }}^T{X}^T\mathbf{y})}{\partial \mathbf{\theta }}=X^T\mathbf{y}$$
对于式子（8）中的 ${\mathbf{y}}^{T}X\mathbf{\theta }$，可以参考上面的第3个式子，即：
$$\frac{\partial ({\mathbf{y}}^{T}X\mathbf{\theta })}{\partial \mathbf{\theta }}=\frac{\partial \left((X^T\mathbf{y}{)}^T\mathbf{\theta }\right)}{\partial \mathbf{\theta }}=X^T\mathbf{y}$$
而式子（8）的 ${\mathbf{y}}^T\mathbf{y}$ 为常数，求导为0，即：
$$\frac{\partial ({\mathbf{y}}^T\mathbf{y})}{\partial \mathbf{\theta }}=0$$
所以有式子（8）中的：
$${\nabla }_{\theta }\left(\frac{1}{2}({\mathbf{\theta }}^T{X}^{T}X\mathbf{\theta }-{\mathbf{\theta }}^T{X}^T\mathbf{y}-{\mathbf{y}}^{T}X\mathbf{\theta }+{\mathbf{y}}^T\mathbf{y})\right)=\frac{1}{2}\left(2X^{T}X\mathbf{\theta }-X^T\mathbf{y}-X^T\mathbf{y}\right)$$

式子（10）式的推导

上面的式子（9）：$$\underset{\mathbf{\theta }}{\min }J(\theta )={\left(X^{T}X\right)}^{-1}{X}^T\mathbf{y}$$
对于 ${\left(X^{T}X\right)}^{-1}$ ，如果存在，则必须行列式 $|X^{T}X|$ > 0
而 $|X^{T}X|$ 无法保证大于 0，可进行一下变换，引入非零向量 $\mathbf{\mu }$
$${\mathbf{\mu }}^T{X}^{T}X\mathbf{\mu }=(X\mathbf{\mu }{)}^T(X\mathbf{\mu })=(X\mathbf{\mu })\bullet (X\mathbf{\mu })\ge 0$$
可得 $X^{T}X$ 为半正定矩阵。而对于非零向量 $\mathbf{\mu }$，有 ${\mathbf{\mu }}^T\mathbf{\mu }$ > 0
所以
$$\begin{gathered} {\mathbf{\mu }}^T{X}^{T}X\mathbf{\mu }+{\mathbf{\mu }}^T\mathbf{\mu }>0 \\ {\mathbf{\mu }}^T{X}^{T}X\mathbf{\mu }+{\mathbf{\mu }}^{T}I\mathbf{\mu }>0 \\ {\mathbf{\mu }}^T{X}^{T}X\mathbf{\mu }+{\mathbf{\mu }}^T\lambda I\mathbf{\mu }>0 \\ {\mathbf{\mu }}^T(X^{T}X+\lambda I)\mathbf{\mu }>0 \end{gathered}$$
所以 $X^{T}X+\lambda I$ 为正定矩阵，即可逆。（因：正定矩阵可逆）
所以式子（9）转换为了式子（10）
$$\mathbf{\theta }={\left(X^{T}X+\lambda I\right)}^{-1}{X}^T\mathbf{y}$$
这样就一定能得最优解。同时，引入了 “$\lambda I$” 能有效防止过拟合。

对于式子（9），还可以从以下角度来理解（仅仅是为了帮助记忆）：
原式子 $X\mathbf{\theta }=\mathbf{y}$， 要解出 $\mathbf{\theta }$ ，必须消掉 X，即：
$$X^{-1}X\mathbf{\theta }=X^{-1}\mathbf{y}$$
如果 X 可逆，首先得满足 X 为方阵，而 $X^{T}X$ 为方阵， 所以，原式子 $X\mathbf{\theta }=\mathbf{y}$ 两边同乘上 $X^T$：
$$X^{T}X\mathbf{\theta }=X^T\mathbf{y}$$
现在出现了 $X^{T}X$ ，两边再同乘 $(X^{T}X{)}^{-1}$ :
$$(X^{T}X{)}^{-1}(X^{T}X)\mathbf{\theta }=(X^{T}X{)}^{-1}{X}^T\mathbf{y}$$
得到最终的解：
$$\mathbf{\theta }=(X^{T}X{)}^{-1}{X}^T\mathbf{y}$$

最小二乘法

我们都学过求解关于 $\mathbf{x}$ 的方程 $A\mathbf{x}=\mathbf{y}$ ，线性回归中是关于 $\theta$ 的方程 $X\mathbf{\theta }=\mathbf{y}$ ，都是一样的。方程有时是没有解或多解的情况，这时，我们采取如下方针：

• 无解的情况下，至少给出一个答案 $\mathbf{x}$ ，使得 $A\mathbf{x}$ 尽量“接近 $\mathbf{y}$”
• 有很多解的情况下，从中选出一个最“合理”的解

要注意，这里我们说的“接近”“合理”都是没有定义的。根据要解决的问题的不同，对这些概念进行度量的方法也有所不同，这就不单单是数学能解决的问题了。

实际中我们经常采用以下标准：

• $A\mathbf{x}-\mathbf{y}$ 的长度越小，我们认为越“接近”
• $\mathbf{x}$ 的长度越小，我们认为越“合理”

（上面提到的长度，就是模长的意思）在这样的标准下，基于以上方针求解问题答案的方法，称为最小二乘法。最小二乘法作为常用的道具，在奇异值分解、广义逆矩阵等的计算中发挥着重要作用。
很显然，线性回归属于上面的第一种情况，是无解的情况。我们的目标就是解出 $\mathbf{\theta }$ 使得 $X\mathbf{\theta }=\hat{\mathbf{y}}$ 尽量接近 $\mathbf{y}$，也就是使误差 $\mathbf{\epsilon }$ 尽量小。