最小二乘法与梯度下降法、牛顿法、高斯-牛顿法的理解

最小二乘法与梯度下降法、牛顿法、高斯-牛顿法的理解

最小二乘法

我们知道，当我们解$n$元线性方程组时，如果有恰好有$n$个方程（假设线性无关），那么可以得出对应的唯一解；二档方程个数大于$n$时，如何确定未知数的值呢？

再举一个更具体的例子：

例如我们收集到了$N$组有关立定跳远成绩（$y$）和身高（$x$）的数据 { ( x i , y j ) ∣ i = 1 , 2 , ⋯   , N } \{ (x_i,y_j)| i=1,2,\cdots, N \} {(xi​,yj​)∣i=1,2,⋯,N}，并希望探寻二者之间的关系关系。假设通过画散点图我们看出二者之间呈线性关系，因而用$y=ax+b$的线性模型进行回归。这本质上也是一个解方程的问题，方程的个数，即数据量（一般）远大于变量个数（2）。我们显然不能直接取其中两点并连线作为拟合曲线，这样的话其余的数据的价值完全没有体现出来，且这样拟合的结果很可能与数据整体呈现的分布趋势大相径庭。而**最小二乘法（Least Square Method，LSM）**就可以实现这种数据分布“整体衡量”或称“全局最优”的拟合。

最小二乘法的核心思想是计算误差函数（Loss Function）
$$\epsilon =\sum _{i=1}^N{\left({\hat{y}}_i-y_i\right)}^2=\sum _{i=1}^N{\left(a{x}_i+b-y_i\right)}^2\text{\hspace{1em}(1.1)}$$
并令其最小，此时对应的参数$a,b$即为所求。

求解时，可以分别对$a,b$求偏导再求零点，从而求得最佳匹配的$a,b$：
$$\{\begin{array}{l}\frac{\partial \epsilon }{\partial a}=2{\sum }_{i=1}^N\left(a{x}_i+b-y_i\right)x_i=0\\ \frac{\partial \epsilon }{\partial b}=2{\sum }_{i=1}^N\left(a{x}_i+b-y_i\right)=0\end{array}\text{\hspace{1em}(1.2)}$$
对于线性情况，我们可以将对应关系写为矩阵形式$\mathbf{A}\vec{x}=\vec{b}$，并得出闭合的解析解
$${\left({\mathbf{A}}^{\mathrm{T}}\mathbf{A}\right)}^{-1}{\mathbf{A}}^{\mathrm{T}}\vec{b}\text{\hspace{1em}(1.3)}$$
对应的MATLAB命令为x = A\b。

实际上，除此之外，对于最小化式$(1.1)$的方法，还有多种，下面尽量以通俗简明的方式，介绍其中三种的基本原理。

梯度下降法

上面的例子中，只有两个未知数，尚且可以使用偏导数置零的方法求解；但当未知数增加时，计算量将会快速增加。此外，LSM只对线性情况有解析解，非线性情况无法求解。在这时我们可以使用梯度下降法（Gradient Descent Method）。其本质是一种迭代求解目标函数最小值的方法，应用在LSM的问题中就可以看作是一种更简单的进行最后一步求解的方法。

基本原理

对于任意标量函数$u(x,y,z)$（这里以三元函数为例），其梯度定义为：
KaTeX parse error: No such environment: align* at position 8: \begin{̲a̲l̲i̲g̲n̲*̲}̲ \operatorname{…
我们知道，梯度代表函数在给定点增加最快的方向，那么，如果我们以一个较小的步长，沿着梯度的负方向逐步前进，就可以逐渐逼近最小值点，这就是梯度下降法的基本思想。

这可以做如下类比：一个人在山坡上，想要以最快的速度下山，那么他就可以寻找山体最陡峭的方向，一小步一小步（防止摔倒）走到山谷处。

这一过程可以用数学工具描述为：
$${\Theta }^{i+1}={\Theta }^i-\alpha \nabla J(\Theta )\text{\hspace{1em}(2.1.2)}$$
其中：

• $\Theta =({\theta }_0,{\theta }_1,\cdots ,{\theta }_n)$；
• $\alpha$为一个较小步长或学习率（较小为了保证不错过最小值点）；
• ${\Theta }^i,{\Theta }^{i+1}$分别代表当前位置和进行计算后的位置。

下面我们来以几个例子具体说明。

单变量的情形

我们先从单变量的简单情形开始讨论，并先将问题与最小二乘法的背景脱离。假设有一个单变量的函数
$$J(\theta )={\theta }^2+2\text{\hspace{1em}(2.2.1)}$$
这是一个简单的抛物线，现在我们要求其最小值。

我们先求其导数（一维情形下，梯度退化为导数）
$$J^{\prime }(\theta )=2\theta \text{\hspace{1em}(2.2.2)}$$
假设迭代起点为${\theta }^0=0$，步长设为$\alpha =0.4$，那么我们可以进行如下的迭代计算：
$$\begin{array}{rl}{\theta }^0& =1\\ {\theta }^1& ={\theta }^0-\alpha J^{\prime }(\theta ){|}_{\theta ={\theta }^0}=1-0.4\times 2\times 1=0.2\\ {\theta }^2& ={\theta }^1-\alpha J^{\prime }(\theta ){|}_{\theta ={\theta }^1}=0.2-0.4\times 2\times 0.2=0.04\\ {\theta }^3& =0.08\\ {\theta }^4& =0.0016\\ \cdots & \end{array}\text{\hspace{1em}(2.1.3)}$$
经过四次迭代，我们基本已经达到了最小值点。

多变量的情形

接下来我们来讨论第一节LSM问题的求解方法：

我们将线性模型的参数$a,b$记为${\theta }_0,{\theta }_1$。同样地，计算误差函数（这里称之为代价函数$J$）：
$$J({\theta }_0,{\theta }_1)=\frac{1}{2N}\sum _{i=1}^N{\left({\hat{y}}_i-y_i\right)}^2$$
其中$\frac{1}{2}$是为了方便之后，与求导形成的“2”相抵消。

然后进行梯度下降法的求解（多变量情形与单变量完全类似，只不过需要将导数换为梯度，且标量计算变为向量计算）：

任意选取一个迭代起点${\Theta }^0=({\theta }_0^0,{\theta }_1^0)$，并选取一个合适大小的步长$\alpha$：
$$\begin{array}{rl}({\theta }_0^1,{\theta }_1^1)& =({\theta }_0^0,{\theta }_1^0)-\alpha J^{\prime }(\Theta ){|}_{\Theta =({\theta }_0^0,{\theta }_1^0)}\\ ({\theta }_0^2,{\theta }_1^2)& =({\theta }_0^1,{\theta }_1^1)-\alpha J^{\prime }(\Theta ){|}_{\Theta =({\theta }_0^1,{\theta }_1^1)}\\ \cdots & \end{array}\text{\hspace{1em}(2.1.4)}$$
即可逐步逼近最小值点，最终对应的${\theta }_0,{\theta }_1$即为拟合参数。

牛顿法

在前面我们提到，为了令$(1.1)$式最小，我们通过求偏导数，转化为求$(1.2)$式的零点。我们可以将这一问题抽象：

现希望求解某一方程$f=0$的根，但并不是所有方程都有求根公式（或求根公式非常复杂），那么我们就可以使用**牛顿法（Newton’s Method）**进行迭代求解。

二维情形

若要求解方程$f(x)=0$的根，由泰勒公式
$$f(x)=\frac{f\left(x_0\right)}{0!}+\frac{f^{\prime }\left(x_0\right)}{1!}\left(x-x_0\right)+\frac{f^{\prime \prime }\left(x_0\right)}{2!}{\left(x-x_0\right)}^2+\dots +\frac{f^{(n)}\left(x_0\right)}{n!}{\left(x-x_0\right)}^n+R_n(x)\text{\hspace{1em}(3.1.1)}$$
将$f(x)$在一个比较接近零点的$x_0$处进行一阶泰勒展开：
$$f(x)\approx f(x_0)+f^{\prime }(x_0)(x-x_0)\text{\hspace{1em}(3.1.2)}$$
并求解$f(x)=0$，得
$$x=x_1=x_0-\frac{f(x_0)}{f^{\prime }(x_0)}\text{\hspace{1em}(3.1.3)}$$
需要注意，这只是一个近似解，但通常可以让$x_1$比$x_0$更加接近真实的零点。继续迭代，我们便可以得出
$$x_{n+1}=x_n-\frac{f(x_n)}{f^{\prime }(x_n)},n\ge 0\text{\hspace{1em}(3.1.4)}$$
当迭代次数足够多时，就可以无限逼近真实零点。需要注意的是，使用牛顿法在求极值时，如果初始点选取不好，则可能不收敛于极小点。

而在最优化问题中，我们一般要求函数$f$的最大、最小值，这可以转化为求解$f’=0$的问题，这与前面完全类似，只不过泰勒展开到二阶即可。对应的迭代公式为：
$$x_{n+1}=x_n-\frac{f^{\prime }(x_n)}{f^{\prime \prime }(x_n)},n\ge 0\text{\hspace{1em}(3.1.5)}$$
这恰恰与前面LSM的思想不谋而合，因而可以用于好解决尤其是非线性情况下的最小二乘问题。

高维情形

与二维情形原理类似，这里直接给出高维情形的迭代公式（梯度取代了一阶导数，而Hessian矩阵取代了二阶导数）：
$$X_{n+1}=X_n-{\mathbf{H}}^{-1}\left(f(X_n)\right)\nabla f(X_n),n\ge 0\text{\hspace{1em}(3.2.1)}$$
其中$\mathbf{H}$为Hessian矩阵：
$$\mathbf{H}(f)=\left[\begin{array}{cccc}\frac{{\partial }^{2}f}{\partial x_1^2}& \frac{{\partial }^{2}f}{\partial x_1\partial x_2}& \cdots & \frac{{\partial }^{2}f}{\partial x_1\partial x_n}\\ \frac{{\partial }^{2}f}{\partial x_2\partial x_1}& \frac{{\partial }^{2}f}{\partial x_2^2}& \cdots & \frac{{\partial }^{2}f}{\partial x_2\partial x_n}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{{\partial }^{2}f}{\partial x_n\partial x_1}& \frac{{\partial }^{2}f}{\partial x_n\partial x_2}& \cdots & \frac{{\partial }^{2}f}{\partial x_n^2}\end{array}\right]\text{\hspace{1em}(3.2.2)}$$

与梯度下降法的比较

• 本质上，牛顿法是二阶收敛，梯度下降是一阶收敛，所以牛顿法收敛更快。例如，我们想找一条最短的路径走到山谷，梯度下降法每次只从当前所处位置选一个坡度最大的方向走一步，而牛顿法在选择方向时，不仅会考虑坡度是否够大，还会考虑你走了一步之后，坡度是否会变得更大。可以说牛顿法比梯度下降法看得更远一点，能更快地走到最底部（牛顿法目光更加长远，少走弯路；相对而言，梯度下降法只考虑了局部的最优，没有全局思想）；

  
  
  红色曲线为利用牛顿法求解，绿色曲线为利用梯度下降法求解

• 但是求解每一步都需要求解目标函数的Hessian矩阵的逆矩阵，比较复杂。

高斯-牛顿法

高斯-牛顿法（Gauss-Newton Method）实际上是牛顿法的在求解非线性最小二乘问题时的一个特例。我们的目的是，在给定$N$组数据点 { ( x i , y j ) ∣ i = 1 , 2 , ⋯   , N } \{ (x_i,y_j)| i=1,2,\cdots, N \} {(xi​,yj​)∣i=1,2,⋯,N}时，拟合一条具有$M$个参数$\mathbf{a}=[a_1,\cdots ,a_M{]}^{\mathrm{T}}$非线性曲线
$$y=f(x,a_1,\cdots ,a_M)\triangleq f(x,\mathbf{a})\text{\hspace{1em}(4.1)}$$
我们定义$f_i(\mathbf{a})=f(x_i,\mathbf{a})$，残差$r_i=y_i-f(x_i,\mathbf{a})=y_i-f(\mathbf{a})(i=1,\cdots ,M)$，并可以将之写为矩阵形式：
$$\mathbf{r}=\left[\begin{array}{c}{r}_1\\ ⋮\\ r_N\end{array}\right]=\left[\begin{array}{c}{y}_1-f_1(\mathbf{a})\\ ⋮\\ y_N-f_N(\mathbf{a})\end{array}\right]=\left[\begin{array}{c}{y}_1\\ ⋮\\ y_N\end{array}\right]-\left[\begin{array}{c}{f}_1(\mathbf{a})\\ ⋮\\ f_N(\mathbf{a})\end{array}\right]=\mathbf{y}-\mathbf{f}(\mathbf{a})\text{\hspace{1em}(4.2)}$$
那么，误差函数就可以写为
$$\begin{array}{r}\epsilon (\mathbf{a})=\sum _{i=1}^N{r}_i^2={\mathbf{r}}^{\mathrm{T}}\mathbf{r}=\Vert \mathbf{r}{\Vert }^2=\Vert \mathbf{y}-\mathbf{f}(\mathbf{a}){\Vert }^2\end{array}\text{\hspace{1em}(4.3)}$$
使得$\epsilon (\mathbf{a})$最小的最佳参数$\mathbf{a}$应满足梯度向量为0：
$$\begin{array}{rl}\frac{\partial }{\partial a_j}\epsilon (\mathbf{a})& =\frac{\partial }{\partial a_j}\sum _{i=1}^N{\left[y_i-f_i(\mathbf{a})\right]}^2=-2\sum _{i=1}^N\left[y_i-f_i(\mathbf{a})\right]\frac{\partial f_i(\mathbf{a})}{\partial a_j}\\ & =-2\sum _{i=1}\left[y_i-f_i(\mathbf{a})\right]J_{ij}=0\end{array}(j=1,\cdots ,M)\text{\hspace{1em}(4.4-1)}$$
对应的向量形式：
$$\mathbf{g}(\epsilon (\mathbf{a}))=\frac{d}{d\mathbf{a}}\epsilon (\mathbf{a})=\frac{d}{d\mathbf{a}}\Vert \mathbf{y}-\mathbf{f}(\mathbf{a}){\Vert }^2=-2{\mathbf{J}}^T(\mathbf{y}-\mathbf{f}(\mathbf{a}))=0\text{\hspace{1em}(4.4-2)}$$
其中$\mathbf{J}$是Jacobian矩阵
$$\mathbf{J}=\left[\begin{array}{cccc}\frac{\partial f_1(\mathbf{a})}{\partial a_1}& \frac{\partial f_1(\mathbf{a})}{\partial a_2}& \cdots & \frac{\partial f_1(\mathbf{a})}{\partial a_M}\\ \frac{\partial f_2(\mathbf{a})}{\partial a_1}& \frac{\partial f_2(\mathbf{a})}{\partial a_2}& \cdots & \frac{\partial f_2(\mathbf{a})}{\partial a_M}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{\partial f_N(\mathbf{a})}{\partial a_1}& \frac{\partial f_N(\mathbf{a})}{\partial a_2}& \cdots & \frac{\partial f_N(\mathbf{a})}{\partial a_M}\end{array}\right]\text{\hspace{1em}(4.5)}$$
但求导置零的方法不一定能得到闭合形式的解，因而我们采用迭代的方法：
$${\mathbf{a}}_{n+1}={\mathbf{a}}_n+\Delta \mathbf{a}\text{\hspace{1em}(4.6)}$$
为了找到$\Delta \mathbf{a}={\mathbf{a}}_{n+1}-{\mathbf{a}}_n={\left[\Delta a_1,\cdots ,\Delta a_M\right]}^{\mathrm{T}}$，我们考虑$f_i({\mathbf{a}}_{n+1})$在${\mathbf{a}}_n$处的一阶泰勒展开式：
$$f_i\left({\mathbf{a}}_{n+1}\right)\approx f_i\left({\mathbf{a}}_n\right)+\sum _{k=1}^M\frac{\partial f_i\left({\mathbf{a}}_n\right)}{\partial a_k}\Delta a_k=f_i\left({\mathbf{a}}_n\right)+\sum _{k=1}^M{J}_{ik}\Delta a_k,(i=1,\cdots ,N)\text{\hspace{1em}(4.7-1)}$$
也就是：
$$\mathbf{f}\left({\mathbf{a}}_{n+1}\right)\approx \mathbf{f}\left({\mathbf{a}}_n\right)+\mathbf{J}\Delta \mathbf{a}\text{\hspace{1em}(4.7-2)}$$
将这一关系代入$(4.4-1)$中：
$$\begin{array}{rl}\sum _{i=1}^N\left[y_i-f_i\left({\mathbf{a}}_{n+1}\right)\right]J_{ij}& \approx \sum _{i=1}^N\left[y_i-f_i\left({\mathbf{a}}_n\right)-\sum _{k=1}^M{J}_{ik}\Delta a_k\right]J_{ij}=0\\ \sum _{i=1}^N{J}_{ij}\sum _{k=1}^M{J}_{ik}\Delta a_k& =\sum _{i=1}^N{J}_{ij}\left[y_i-f_i\left({\mathbf{a}}_n\right)\right]\end{array}(j=1,\cdots ,M)\text{\hspace{1em}(4.8-1)}$$
写为矩阵形式：
$$\left({\mathbf{J}}^{\mathrm{T}}\mathbf{J}\right)△\mathbf{a}={\mathbf{J}}^{\mathrm{T}}\left(\mathbf{y}-\mathbf{f}\left({\mathbf{a}}_n\right)\right)\text{\hspace{1em}(4.8-2)}$$
这样的好处是可以方便地对$\mathbf{a}$进行求解：
$$\Delta \mathbf{a}={\mathbf{a}}_{n+1}-{\mathbf{a}}_n={\left({\mathbf{J}}^T\mathbf{J}\right)}^{-1}{\mathbf{J}}^T\left(\mathbf{y}-\mathbf{f}\left({\mathbf{a}}_n\right)\right)={\mathbf{J}}^{-}\left(\mathbf{y}-\mathbf{f}\left({\mathbf{a}}_n\right)\right)\text{\hspace{1em}(4.9)}$$
其中${\mathbf{J}}^{-}={\left({\mathbf{J}}^{\mathrm{T}}\mathbf{J}\right)}^{-1}{\mathbf{J}}^{\mathrm{T}}$为$\mathbf{J}$的伪逆矩阵。这样我们就得到了迭代表达式：
$${\mathbf{a}}_{n+1}={\mathbf{a}}_n+\Delta \mathbf{a}={\mathbf{a}}_n+{\mathbf{J}}^{-}\left(\mathbf{y}-\mathbf{f}\left({\mathbf{a}}_n\right)\right)={\mathbf{a}}_n-{\mathbf{J}}^{-}\left(\mathbf{f}\left({\mathbf{a}}_n\right)-\mathbf{y}\right)\text{\hspace{1em}(4.10)}$$
有了迭代公式，我们自然就可以更简便地求解导数置零的问题，应用于最小二乘问题中也就可以大大简化计算并解决LSM无法解决的非线性问题。

参考文献

• 最小二乘法和梯度下降法有哪些区别？ - 任妍Carol的回答 - 知乎

• 梯度下降算法原理讲解——机器学习

• Newton法（牛顿法 Newton Method）

• 【math】梯度下降法(梯度下降法，牛顿法，高斯牛顿法，Levenberg-Marquardt算法)

• Gauss-Newton algorithm for nonlinear models