最小二乘法

前言

参考博客：半小时学习最小二乘法
最小二乘法在统计学的地位不必多言。本文的目的是全面地讲解最小二乘法，打好机器学习的基础。本文主要内容是最小二乘法的思想及在线性回归问题中的应用。后面的系列文章会继续讲解最小二乘的正则化。
至于非线性最小二乘和广义线性模型，如果以后有时间会进行整理。

核心思想

最小二乘法是勒让德( A. M. Legendre)于1805年在其著作《计算慧星轨道的新方法》中提出的。它的主要思想就是求解未知参数，使得理论值与观测值之差（即误差，或者说残差）的平方和达到最小：
$$E=\sum _{i=1}^n{e}_i^2=\sum _{i=1}^n{\left(y_i-\hat{y}\right)}^2$$
观测值 $y_i$ 就是我们的多组样本，理论值 $\hat{y}$ 就是我们的假设拟合函数。目标函数也就是在机器学习中 常说的损失函数 $E$ ，我们的目标是得到使目标函数最小化时候的参数。

所谓最小二乘，其实也可以叫做最小平方和，其目的就是通过最小化误差的平方和，使得拟合对象无限接近目标对象。换句话说，最小二乘法定义了一种函数的拟合标准，其目标是最小化误差的平方和。

直观理解

均方误差有非常好的几何意义，它对应了常用的欧几里德距离。在线性回归中，最小二乘法就是试图找到一条直线，使所有样本到直线的欧氏距离之和最小。

假设有一条直线 $y=ax+b$ ，要在这条直线上找到一点，距离 $\left(x_0,y_0\right)$ 这个点的距离最短。如果用绝 对值的方法寻找，也就是取 $\min \left(\left|y-y_0\right|+\left|x-x_0\right|\right)$ ，由于绝对值最小为0，所以最小的情况就是 $x=x_0$ 或者 $y=y_0$ 处，如下图1所示。

$$图1$$
如果用平方和的方法寻找，就是取min ${\left(y-y_0\right)}^2+{\left(x-x_0\right)}^2$ ，可以看出该式是两点间距离公式， 也就是距离的概念。那么最短的距离，就是点到直线的垂线，如下图2所示。

$$图2$$
事实上，线性回归中最小二乘法的解 $\theta ={\left(X^{T}X\right)}^{-1}{X}^{T}Y$ 也就是投影矩阵的公式: 将Y向量投影到X 构成的平面上。

对任意函数的通用解法

1. 列出损失函数$E={\sum }_{i=1}^n{e}_i^2={\sum }_{i=1}^n{\left(y_i-\hat{y}\right)}^2$
2. 根据损失函数对参数应用多元函数的求极值方法，直接求解函数最小值。而更常见的方法即是将损失函数$y_i$用$x_i$和参数表示，然后使用梯度下降算法。
3. 求得函数最小值的参数或待到梯度算法收敛，此时的参数即为所求这些个步骤说起来抽象，实际上这是在机器学习中应用最广泛的方法。但是对于后面的线性回归问题，有着更简洁的推导方法。

以算术平均值为例——为什么算术平均即是真值

可以说整部数理统计学的历史，就是对算术平均不断深入研究的历史。而最小二乘法可以解释为什么多次测量取算术平均的结果就是真值，比如估计身高可以测三次后取平均。

当我们对于某个末知量 $\theta$ 观测 $m$ 次，记每次的结果为 $x_i$
$$E=\sum _{i=1}^m{e}_i^2=\sum _{i=1}^m{\left(x_i-\theta \right)}^2$$
求导得
$$\sum _{i=1}^n-\left(x_i-\theta \right)=0$$
所以 $\theta =\bar{x}=\frac{{\sum }^{x_i}}{m}$

线性回归问题的定义

这部分我们简要回顾一下基本的线性代数知识。

我们知道一条直线，在三维空间中，可以用形如$ax+by+cz=d$的方程表示，而用矩阵的形式表达的话，即为$AX=B$，其中$B$为常数，而$A$和$X$都为n维向量。那么多个这样的方程联立在一起就被称为线性方程组，其中有两个基本问题：
1）方程组是否有解，即解的存在性问题？
2）如果有解，解的个数有多少个？

为了回答这两个问题，我们令线性方程组中的方程都为线性无关方程（可以通过简单地初等行变换），方程个数为m，特征数（参数数量）为n。在这样的假设下，如果m=n，则一行方程对应一个参数的解，此时方程组有唯一解。如果m>n，则方程无解，因为存在互相矛盾的两个方程。如果m<n，则方程和参数不能一一对应，存在无穷多解。而最小二乘法是在m>n时可以使用的，其通过让残差平方和最小，找到互相矛盾解之间的近似解。

使用最小二乘求解线性回归问题

对于线性回归问题，当然可以使用求导的代数方法来找到损失函数的最小值。但矩阵法比代数法要简洁，所以现在很多书和机器学习库都是用的矩阵法来做最小二乘法，本文这里介绍一下如何使用矩阵法求解线性回归问题。

对于函数 $h_{\theta }\left(x_1,x_2,\dots x_n\right)={\theta }_0+{\theta }_1{x}_1+\dots +{\theta }_n{x}_n$ ，我们将其的矩阵形式记作:
$$X\theta =Y$$

$$x^{(i)}=\left(\begin{array}{c}{x}_0^{(i)}\\ x_1^{(i)}\\ ⋮\\ x_n^{(i)}\end{array}\right)\theta =\left(\begin{array}{c}{\theta }_0\\ {\theta }_1\\ ⋮\\ {\theta }_n\end{array}\right)X=\left(\begin{array}{c}{\left(x^{(1)}\right)}^T\\ {\left(x^{(2)}\right)}^T\\ ⋮\\ {\left(x^{(m)}\right)}^T\end{array}\right)Y=\left(\begin{array}{c}{y}^{(1)}\\ y^{(2)}\\ ⋮\\ y^{(m)}\end{array}\right)$$
故损失函数根据定义将$Y$用$X$和$\theta$代替：（系数1/2是为了简化计算添加的，求迹前和求迹后值不变）
$$J(\theta )=\frac{1}{2}\sum _{i=1}^m{\left(h_{\theta }\left(x^{(i)}\right)-y^{(i)}\right)}^2=\frac{1}{2}\mathrm{tr}\left[(X\theta -Y{)}^T(X\theta -Y)\right]$$
应用矩阵迹的计算公式：
$$\begin{array}{rl}\frac{\partial J(\theta )}{\partial \theta }& =\frac{1}{2}\cdot \frac{\partial \mathrm{tr}\left({\theta }^T{X}^{T}X\theta -{\theta }^T{X}^{T}Y-Y^{T}X\theta +Y^{T}Y\right)}{\partial \theta }\\ & =\frac{1}{2}\cdot \left[\frac{\partial \mathrm{tr}\left(\theta I{\theta }^T{X}^{T}X\right)}{\partial \theta }-\frac{\partial \mathrm{tr}\left({\theta }^T{X}^{T}Y\right)}{\partial \theta }-\frac{\partial \mathrm{tr}\left(\theta Y^{T}X\right)}{\partial \theta }\right]\\ & =\frac{1}{2}\cdot \left[X^{T}X\theta I+{\left(X^{T}X\right)}^T\theta I^T-X^{T}Y-{\left(Y^{T}X\right)}^T\right]\\ & =X^{T}X\theta -X^{T}Y\end{array}$$
令上式为 0 ，解得 $\theta ={\left(X^{T}X\right)}^{-1}{X}^{T}Y$

Note：矩阵求导坑多，使用迹来计算比较方便

线性回归的t检验
记n为回归方程的特征个数，m为样本数
$$\begin{array}{rl}& S_{\text{回 }}=\sum _{i=1}^m(\hat{y}-\bar{y})\\ & S_{\text{剩 }}=\sum _{i=1}^m\left(y_i-\hat{y}\right)\end{array}$$
总平方和（SST）可分解为回归平方和（SSR）与残差平方和（SSE）两部
$$\begin{array}{c}MSR=SSR/k\\ MSE=SSE/(n-k-1)\\ F=\frac{MSR}{MSE}=\frac{S_{\text{回 }}/k}{S_{\text{剩 }}/(n-k-1)}\end{array}$$
若用样本计算的F>F0.05(k,n−k−1)，则拒绝H0，则回归方程在显著性水平α＝0.05下是显著的

最小二乘法的适用场景
当样本量m很少，小于特征数n的时候，这时拟合方程是欠定的，需要使用LASSO。当m=n时，用方程组求解。当m>n时，拟合方程是超定的，我们可以使用最小二乘法。

局限性

• 首先，最小二乘法需要计算(XTX)−1逆矩阵，有可能逆矩阵不存在，这样就没有办法直接用最小二乘法。
• 第二，当样本特征n非常的大的时候，计算逆矩阵是一个非常耗时的工作，甚至不可行。建议不超过10000个特征。
• 第三，如果拟合函数不是线性的，这时无法使用最小二乘法，需要通过一些技巧转化为线性才能使用。

最小二乘法和M估计
在统计数据时，难免会遇到异常值，即人为误差。而这种误差对结果的影响远比系统误差大，比如将1记录成10。所以我们使用稳健性来评价一个方法对异常值的敏感程度。

最小二乘法是一种稳健性较差的方法，原因在于其目标函数是误差的平方，是一个增长很快的函数。
所以不难想到，对于$E={\Sigma }^f(x_i)E={\Sigma }^f(x_i)$，我们可以取$f(x)=|x|$来减小函数的增长速度。

统计学家休伯将这一想法用于对单个末知量 $\theta$ 参数估计的情况，误差满足正态分布 $x_i=\theta +e_i$ ，就给 定 $\rho$ 函数:
取定函数 $\rho$ ，找出使函数 $M(\theta )={\sum }_{i=1}^m\rho \left(x_i-\theta \right)$ 达到最小的 $\hat{\theta }$ ，将其作为 $\theta$ 的估计值。 $\hat{\theta }$ 称为 $\theta$ 的 $M$ 估计。

M估计是一类估计，主要包括 $\rho (u)=u^2$ 的最小二乘法和 $\rho (u)=|x|$ 的最小一乘法。M估计也可以和 最小二乘法一样，推广到多元线性回归，称为稳健回归，但是因为难于计算等局限，应用并不广泛。

Note: 最小一乘法对末知参数 $\theta$ 的估计值 $\hat{\theta }=x_i$ 的中位数

最小二乘法和正则化

当 ${\left(X^{T}X\right)}^{-1}$ 不存在，即 $X_{T}X$ 不满秩时， $\theta$ 无唯一解。

故考虑在原先的A的最小二乘估计中加一个小扰动 $\lambda I$ ，使原先无法求广义逆的情况变成可以求出其广 义逆，使得问题稳定并得以求解。有:
$$\hat{\theta }={\left(X^{T}X+\lambda I\right)}^{-1}{X}^{T}Y$$
而此时对应的损失函数为:
$$J(\theta )=\sum _{i=1}^m{\left(y_i-{\theta }^T{x}_i\right)}^2+\lambda \Vert \theta {\Vert }_2^2$$

上式称为岭回归（ridge regression），通过引入L2范数正则化。当然也可以将L2范数替换为L1范数。对应有
$$J(\theta )=\sum _{i=1}^m{\left(y_i-{\theta }^T{x}_i\right)}^2+\lambda \Vert \theta {\Vert }_1$$

上式称为LASSO。对于L2范数，本质上其实是对误差的高斯先验。而L1范数则对应于误差的Laplace先验。

最小二乘法的理论证明

拉普拉斯指出随机误差应满足正态分布，而高斯创造性地发明并使用极大似然法证明了最小二乘法。

故测量误差服从高斯分布的情况下， 最小二乘法等价于极大似然估计。

对于任何拟合方程都有: $y=\hat{y}+e$
因为 $e\sim N\left(0,{\sigma }^2\right)$, 故 $y\sim N\left(\hat{y},{\sigma }^2\right)$
由极大似然估计， $L(\theta )={\prod }_{i}f\left(y_i\right)$
$$L=\frac{1}{(\sqrt{2\pi }\sigma {)}^n}\exp \left\{-\frac{1}{2{\sigma }^2}\sum _{i=1}^n{\left(y_i-\hat{y}\right)}^2\right\}=\frac{1}{(\sqrt{2\pi }\sigma {)}^n}\exp -\frac{1}{2{\sigma }^2}\sum _{i=1}^n{e}_i^2$$

Note：数学的发展史很多时候是不符合逻辑顺序的。事实上，高斯当时是循环论证最小二乘法的，推理有缺陷。而后拉普拉斯才弥补了这一点。