最小二乘法

最小二乘法是用来做函数拟合或者求函数极值的方法。

最小二乘法的原理与要解决的问题

最小二乘法是由勒让德在 19 世纪发现的，原理的一般形式很简单，当然发现的过程是非常艰难的。形式如下式：

$$目标函数=\sum (观测值−理论值)^2$$

观测值就是我们的多组样本，理论值就是我们的假设拟合函数, 目标函数也就是在机器学习中常说的损失函数，我们的目标是得到使目标函数最小化时候的拟合函数的模型。举一个最简单的线性回归的简单例子，比如我们有 m 个只有一个特征的样本：

$$(x^{(1)},y^{(1)}), (x^{(2)},y^{(2)}),...(x^{(m)},y^{(m)})$$

样本采用下面的拟合函数：

$$h_\theta(x) = \theta_0 + \theta_1 x$$

这样我们的样本有一个特征 x，对应的拟合函数有两个参数$\theta_0$和$\theta_1$需要求出。

我们的目标函数为：

$$J(\theta_0, \theta_1) = \sum\limits_{i=1}^{m}(y^{(i)} - h_\theta(x^{(i)})^2 = \sum\limits_{i=1}^{m}(y^{(i)} - \theta_0 - \theta_1 x^{(i)})^2$$

为了使$J(\theta_0, \theta_1)$的值最小，使用最小二乘法

最小二乘法的代数法解法

假设$J(a, b) = \sum_{i=1}^m(y^{(i)} - ax^{(i)} - b)^2$, 要获得该函数的极值，即对该函数的各个未知分量求导，使其导数为 0。即：

$\frac {\partial J(a, b)}{\partial a} = 0$， $\frac {\partial J(a, b)}{\partial b} = 0$

对 b 求导：

$$\begin{align} \frac {\partial J(a, b)}{\partial b} = & \sum_{i=1}^m2(y^{(i)} - ax^{(i)} - b)(-1) = 0 \\ & \sum_{i=1}^m(y^{(i)} - ax^{(i)} - b) = 0 \\ & \sum_{i=1}^my^{(i)} - a\sum_{i=1}^mx^{(i)} - \sum_{i=1}^mb = 0\\ & \sum_{i=1}^my^{(i)} - a\sum_{i=1}^mx^{(i)} - mb = 0\\ & \sum_{i=1}^my^{(i)} - a\sum_{i=1}^mx^{(i)} = mb \Rightarrow b = \bar y - a\bar x (\bar x和\bar y分别为x和y的均值)\\ \end{align}$$

对 a 求导：

$$\begin{align} \frac {\partial J(a, b)}{\partial a} = & \sum_{i=1}^m2(y^{(i)} - ax^{(i)} - b)(-x^{(i)}) = 0 \\ & \sum_{i=1}^m(y^{(i)} - ax^{(i)} - b)x^{(i)} = 0 \\ & \sum_{i=1}^m(y^{(i)} - ax^{(i)} - \bar y + a\bar x)x^{(i)} = 0\\ & \sum_{i=1}^m(x^{(i)}y^{(i)} - a(x^{(i)})^2 - x^{(i)}\bar y + a\bar xx^{(i)}) = 0\\ & \sum_{i=1}^m(x^{(i)}y^{(i)} - x^{(i)}\bar y - a(x^{(i)})^2 + a\bar xx^{(i)}) = 0\\ & \sum_{i=1}^m(x^{(i)}y^{(i)} - x^{(i)}\bar y) - \sum_{i=1}^m(a(x^{(i)})^2 - a\bar xx^{(i)}) = 0\\ & \sum_{i=1}^m(x^{(i)}y^{(i)} - x^{(i)}\bar y) - a\sum_{i=1}^m((x^{(i)})^2 - \bar xx^{(i)}) = 0\\ & a = \frac{\sum_{i=1}^m(x^{(i)}y^{(i)} - x^{(i)}\bar y)}{\sum_{i=1}^m((x^{(i)})^2 - \bar xx^{(i)})} \end{align}$$

对于公式$\sum_{i=1}^mx^{(i)}\bar y$ 可以进一步推导

$$\begin{align} \sum_{i=1}^mx^{(i)}\bar y & = \bar y\sum_{i=1}^mx^{(i)} \\ & = m\bar y\bar x = \sum_{i=1}^m\bar x\bar y\\ & = \bar x \sum_{i=1}^my^{(i)}\\ & = \sum_{i=1}^m\bar xy^{(i)} \\ \end{align}$$

带入 a 中可得到

$$\begin{align} a = \frac{\sum_{i=1}^m(x^{(i)}y^{(i)} - x^{(i)}\bar y)- \bar xy^{(i)} + \bar x \bar y)}{\sum_{i=1}^m((x^{(i)})^2 - \bar xx^{(i)} - \bar xx^{(i)} + \bar x^2)} & = \frac{\sum_{i=1}^m(x^{(i)} - \bar x)(y^{(i)} - \bar y)}{\sum_{i=1}^m(x^{(i)} - \bar x)^2} \end{align}$$

总结: 通过上面的推导可以得到 a 和 b 的值

$$a = \frac{\sum_{i=1}^m(x^{(i)} - \bar x)(y^{(i)} - \bar y)}{\sum_{i=1}^m(x^{(i)} - \bar x)^2}, b = \bar y - a\bar x$$

这个方法很容易推广到多个样本特征的线性拟合。拟合函数表示为$h_\theta(x_1, x_2, ...x_n) = \theta_0 + \theta_{1}x_1 + ... + \theta_{n}x_{n}$, 其中$\theta_i$ (i = 0,1,2... n)为模型参数，$x_i$ (i = 0,1,2... n)为每个样本的 n 个特征值。这个表示可以简化，我们增加一个特征$x_0 = 1$ ，这样拟合函数表示为：

$$J(\theta_0, \theta_1..., \theta_n) = \sum\limits_{j=1}^{m}(h_\theta(x_0^{(j)}), x_1^{(j)}, ...x_n^{(j)})) - y^{(j)}))^2 = \sum\limits_{j=1}^{m}(\sum\limits_{i=0}^{n}\theta_{i}x_{i}^{(j)}- y^{(j)})^2$$

利用损失函数分别对$\theta_i$(i=0,1,...n)求导,并令导数为 0 可得：

$$\sum\limits_{j=0}^{m}(\sum\limits_{i=0}^{n}(\theta_{i}x_{i}^{(j)} - y^{(j)})x_i^{(j)}= 0 (i=0,1,...n)$$

这样我们得到一个N+1元一次方程组，这个方程组有N+1个方程，求解这个方程，就可以得到所有的N+1个未知的$\theta$

最小二乘法的矩阵法解法

矩阵法比代数法要简洁，且矩阵运算可以取代循环，所以现在很多书和机器学习库都是用的矩阵法来做最小二乘法。这里用上面的多元线性回归例子来描述矩阵法解法。

假设函数$h_\theta(x_1, x_2, ...x_n) = \theta_0 + \theta_{1}x_1 + ... + \theta_{n-1}x_{n-1}$的矩阵表达方式为：

$$h_\mathbf{\theta}(\mathbf{x}) = \mathbf{X\theta}$$

其中， 假设函数$h_\mathbf{\theta}(\mathbf{X})$为 mx1 的向量,$\mathbf{\theta}$为 nx1 的向量，里面有 n 个代数法的模型参数。$\mathbf{X}$为 mxn 维的矩阵。m 代表样本的个数，n 代表样本的特征数。

损失函数定义为:

$$J(\mathbf\theta) = \frac{1}{2}(\mathbf{X\theta} - \mathbf{Y})^T(\mathbf{X\theta} - \mathbf{Y})$$

其中$\mathbf{Y}$是样本的输出向量，维度为 mx1. $\frac{1}{2}$在这主要是为了求导后系数为 1，方便计算。

根据最小二乘法的原理，我们要对这个损失函数对$\theta$向量求导取 0。结果如下式：

$$\frac{\partial}{\partial\mathbf\theta}J(\mathbf\theta) = \mathbf{X}^T(\mathbf{X\theta} - \mathbf{Y}) = 0$$

这里面用到了矩阵求导链式法则，和两个矩阵求导的公式。

公式 1：$\frac{\partial}{\partial\mathbf{x}}(\mathbf{x^Tx}) =2\mathbf{x}\;\;x为向量$ 公式 2: $\nabla_Xf(AX+B) = A^T\nabla_Yf,\;\; Y=AX+B,\;\;f(Y)为标量$

对上述求导等式整理后可得： $\mathbf{X^{T}X\theta} = \mathbf{X^{T}Y}$

两边同时左乘$(\mathbf{X^{T}X})^{-1}$可得：

$\mathbf{\theta} = (\mathbf{X^{T}X})^{-1}\mathbf{X^{T}Y}$

这样就求出了$\theta$向量表达式的公式，免去了代数法一个个去求导的麻烦。只要给了数据,我们就可以用 $\mathbf{\theta} = (\mathbf{X^{T}X})^{-1}\mathbf{X^{T}Y}$ 算出$\theta$

最小二乘法的局限性和适用场景

最小二乘法适用简洁高效，比梯度下降这样的迭代法似乎方便很多。但是这里我们就聊聊最小二乘法的局限性。

首先，最小二乘法需要计算$\mathbf{X^{T}X}$的逆矩阵，有可能它的逆矩阵不存在，这样就没有办法直接用最小二乘法了，此时梯度下降法仍然可以使用。当然，我们可以通过对样本数据进行整理，去掉冗余特征。让$\mathbf{X^{T}X}$的行列式不为 0，然后继续使用最小二乘法。

第二，当样本特征 n 非常的大的时候，计算$\mathbf{X^{T}X}$的逆矩阵是一个非常耗时的工作（nxn 的矩阵求逆），甚至不可行。此时以梯度下降为代表的迭代法仍然可以使用。那这个 n 到底多大就不适合最小二乘法呢？如果你没有很多的分布式大数据计算资源，建议超过 10000 个特征就用迭代法吧。或者通过主成分分析降低特征的维度后再用最小二乘法。

第三，如果拟合函数不是线性的，这时无法使用最小二乘法，需要通过一些技巧转化为线性才能使用，此时梯度下降仍然可以用。

第四，讲一些特殊情况。当样本量 m 很少，小于特征数 n 的时候，这时拟合方程是欠定的，常用的优化方法都无法去拟合数据。当样本量 m 等于特征数 n 的时候，用方程组求解就可以了。当 m 大于 n 时，拟合方程是超定的，也就是我们常用与最小二乘法的场景了。