最小二乘：梯度下降法、牛顿法、高斯-牛顿法

1.线性预测器最佳预测系数线性方程推导

2.最小二乘法

2.1 简介

最小二乘法（又称最小平方法）是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。

2.2 举例

比如说，我们知道夏天来了，天气一热我们就想吃冰棍，那么温度（x）和冰棍的销量（y）应该有着某种线性关系，我们假设这种线性关系为:
$f(x)=ax+b$f(x)=ax+b
通过最小二乘法的思想，总误差为：
$ε= ∑_{i=1}^{N}[f(x_i)-y_i]^2=∑_{i=1}^{N}[ax_i+b-y_i]^2$ε=i=1∑N​[f(xi​)−yi​]2=i=1∑N​[axi​+b−yi​]2
不同的a,b 会导致不同的 ε，根据多元微积分的知识，当：
$\begin{cases} \frac{\partial}{\partial a}\epsilon=2\sum (ax_i+b-y_i)x_i=0\\ \quad\\ \frac{\partial}{\partial b}\epsilon=2\sum (ax_i+b-y_i)=0\end{cases}$⎩⎪⎨⎪⎧​∂a∂​ϵ=2∑(axi​+b−yi​)xi​=0∂b∂​ϵ=2∑(axi​+b−yi​)=0​

这个时候 ε 取最小值。
对于a,b 而言，上述方程组为线性方程组，用之前的数据解出来：
$\begin{cases} a\approx m\\ \quad\\ b\approx n\end{cases}$⎩⎪⎨⎪⎧​a≈mb≈n​

从而求得了最佳匹配的a和b，也就是数m和n。

上面我们假设了 $f(x)=ax+b$f(x)=ax+b 进行求解，我们其实也可以选用 $f(x)=ax^2+bx+c$f(x)=ax2+bx+c 来进行拟合，结果与上式不同。其实，最小二乘的方法，还有多种，下面介绍其中三种的基本原理。

2.3 梯度下降法

2.3.1 梯度

在微积分里面，对多元函数的参数求偏导数，把求得的各个参数的偏导数以向量的形式写出来，就是梯度。比如函数f(x,y), 分别对x,y求偏导数，求得的梯度向量就是$(\frac{∂f}{∂x}, \frac{∂f}{∂y})$(∂x∂f​,∂y∂f​)T,简称grad f(x,y)或者▽f(x,y)。对于在点(x0,y0)的具体梯度向量就是$(\frac{∂f}{∂x_0}, \frac{∂f}{∂y_0})$(∂x0​∂f​,∂y0​∂f​)T
那么这个梯度向量求出来有什么意义呢？他的意义从几何意义上讲，就是函数变化增加最快的地方。具体来说，对于函数f(x,y),在点$(x_0,y_0)$(x0​,y0​)，沿着梯度向量的方向就是$(\frac{∂f}{∂x_0}, \frac{∂f}{∂y_0})$(∂x0​∂f​,∂y0​∂f​)T的方向是f(x,y)增加最快的地方。或者说，沿着梯度向量的方向，更加容易找到函数的最大值。反过来说，沿着梯度向量相反的方向，梯度减少最快，也就是更加容易找到函数的最小值。

2.3.2 算法详解

首先来看看梯度下降的一个直观的解释。比如我们在一座大山上的某处位置，由于我们不知道怎么下山，于是决定走一步算一步，也就是在每走到一个位置的时候，求解当前位置的梯度，沿着梯度的负方向，也就是当前最陡峭的位置向下走一步，然后继续求解当前位置梯度，向这一步所在位置沿着最陡峭最易下山的位置走一步。这样一步步的走下去，一直走到觉得我们已经到了山脚。当然这样走下去，有可能我们不能走到山脚，而是到了某一个局部的山峰低处。

从上面的解释可以看出，梯度下降不一定能够找到全局的最优解，有可能是一个局部最优解。当然，如果损失函数是凸函数，梯度下降法得到的解就一定是全局最优解。

梯度下降法的先决条件： 确认优化模型的假设函数和损失函数

假设我们的优化模型采用线性回归，假设函数表示为$h_\theta(x_1, x_2, ...x_n) = \theta_0 + \theta_{1}x_1 + ... + \theta_{n}x_{n}$hθ​(x1​,x2​,...xn​)=θ0​+θ1​x1​+...+θn​xn​

其中，θi为模型参数，$x_i$xi​为每个样本的n个特征值

对应于上面的假设函数，损失函数为：
$J(\theta_0, \theta_1..., \theta_n) = \frac{1}{2m}\sum\limits_{j=1}^{m}(h_\theta(x_0^{(j)}, x_1^{(j)}, ...x_n^{(j)}) - y_j)^2$J(θ0​,θ1​...,θn​)=2m1​j=1∑m​(hθ​(x0(j)​,x1(j)​,...xn(j)​)−yj​)2

算法相关参数初始化：主要是初始化$\theta_0, \theta_1..., θ_n$θ0​,θ1​...,θn​,算法终止距离????以及步长????。在没有任何先验知识的时候，我们可以先将所有的????初始化为0，将步长初始化为1。

算法过程：

1）确定当前位置的损失函数的梯度，对于????????,其梯度表达式如下:
$\frac{\partial}{\partial\theta_i}J(\theta_0, \theta_1..., θ_n)$∂θi​∂​J(θ0​,θ1​...,θn​)

2）用步长乘以损失函数的梯度，得到当前位置下降的距离 $\alpha\frac{\partial}{\partial\theta_i}J(\theta_0, \theta_1..., θ_n)$α∂θi​∂​J(θ0​,θ1​...,θn​)

3）确定是否所有的 ???????? ,梯度下降的距离都小于????，如果小于????则算法终止，当前所有的 ???????? (i=0,1,…n) 即为最终结果。否则进入步骤4.

4）更新所有的????，对于 ????????，其更新表达式如下。更新完毕后继续转入步骤1.
$\theta_i = \theta_i - \alpha\frac{\partial}{\partial\theta_i}J(\theta_0, \theta_1..., θ_n)$θi​=θi​−α∂θi​∂​J(θ0​,θ1​...,θn​)

2.3.3 算法调优

在使用梯度下降时，一些地方是可以进行调优的，以使梯度下降的

1. 算法的步长选择。在前面的算法描述中，取步长为1，但是实际上取值取决于数据样本，可以多取一些值，从大到小，分别运行算法，看看迭代效果，如果损失函数在变小，说明取值有效，否则要增大步长。步长太大，会导致迭代过快，甚至有可能错过最优解。步长太小，迭代速度太慢，很长时间算法都不能结束。所以算法的步长需要多次运行后才能得到一个较为优的值。

2. 算法参数的初始值选择。 初始值不同，获得的最小值也有可能不同，因此梯度下降求得的只是局部最小值；由于有局部最优解的风险，需要多次用不同初始值运行算法，关键损失函数的最小值，选择损失函数最小化的初值。

3. 归一化。由于样本不同特征的取值范围不一样，可能导致迭代很慢，为了减少特征取值的影响，可以对特征数据归一化

2.4 牛顿法

在前面的冰棍例子中我们提到，我们可以通过求偏导数来求得a,b的最优解，但这存在局限性，对于非线性优化问题，牛顿法提供了一种求解的办法。

假设任务是优化一个目标函数f，求函数f=0的问题，并不是所有的方程都有求根公式，或者求根公式很复杂，导致求解困难。利用牛顿法，可以迭代求解。

利用泰勒公式，在$x_0$x0​处展开，且展开到一阶:
$f(x) = f(x_0)+(x－x_0)f'(x_0)$f(x)=f(x0​)+(x－x0​)f′(x0​)

求解f(x)=0，得：
$x = x_1=x_0－\frac{f(x_0)}{f'(x_0)}$x=x1​=x0​－f′(x0​)f(x0​)​

因为这是利用泰勒公式的一阶展开，$f(x) = f(x_0)+(x－x_0)f'(x_0)$f(x)=f(x0​)+(x－x0​)f′(x0​)并不是完全相等，而是近似相等，这里求得的$x_1$x1​并不能让f(x)=0，只能说$f(x_1)$f(x1​)的值比$f(x_0)$f(x0​)更接近f(x)=0，于是，迭代求解的想法就很自然了，可以进而推出 $x(n+1)=x(n)－f(x(n))/f'(x(n))$x(n+1)=x(n)－f(x(n))/f′(x(n))，通过迭代，这个式子必然在f(x*)=0的时候收敛。整个过程如下图：

在上面讨论的是2维情况，高维情况的牛顿迭代公式是：

其中，H为Hessian矩阵：

一般认为牛顿法可以利用到曲线本身的信息，比梯度下降法更容易收敛（迭代更少次数），如下图是一个最小化一个目标方程的例子，红色曲线是利用牛顿法迭代求解，绿色曲线是利用梯度下降法求解。可以看出，牛顿法“少走了弯路”，但是问题在于Hessian矩阵引入的复杂性，使得牛顿迭代求解的难度大大增加，但是已经有了解决这个问题的办法就是Quasi-Newton methond（拟牛顿法）,不再直接计算Hessian矩阵，而是每一步的时候使用梯度向量更新Hessian矩阵的近似。

拟牛顿法
因为我们要选择一个矩阵来代替海森矩阵的逆，那么我们首先要研究一下海森矩阵需要具有什么样的特征才能保证牛顿法成功的应用。通过上面的描述我们知道
上式我们称之为拟牛顿条件。
因此，对于我们所选择的替代矩阵$G_k$Gk​，需要满足两个条件：
1、拟牛顿条件，即$G_k(f′(x_k+1)−f′(x_k))=x_{k+1}−x_k；$Gk​(f′(xk​+1)−f′(xk​))=xk+1​−xk​；
2、要保证$G_k$Gk​为正定矩阵，这是因为只有正定才能保证牛顿法的搜索方向是向下搜索的
假设$y_k=f′(x_{k+1})−f′(x_k)$yk​=f′(xk+1​)−f′(xk​)，$δ_k=x_{k+1}−x_k$δk​=xk+1​−xk​，因为每次迭代我们都需要更新替代矩阵$G_k$Gk​，下面介绍一种常用的迭代算法DFP(Davidon-Fletcher-Powell)

2.5 高斯-牛顿法

高斯–牛顿迭代法实际上是牛顿法的在求解非线性最小二乘问题时的一个特例子，它的基本思想是使用泰勒级数展开式去近似地代替非线性回归模型，然后通过多次迭代，多次修正回归系数，使回归系数不断逼近非线性回归模型的最佳回归系数，最后使原模型的残差平方和达到最小。

当目标函数为：
$S(\boldsymbol \beta)= \sum_{i=1}^m r_i^2(\boldsymbol β)$S(β)=i=1∑m​ri2​(β)

该函数趋于0，直接使用高维的牛顿法来求解，迭代公式为：

$\boldsymbol\beta^{(s+1)} = \boldsymbol\beta^{(s)} - \mathbf H^{-1} \mathbf g \,$β(s+1)=β(s)−H−1g

和牛顿法中的最优化问题高维迭代公式一样

目标函数可以简写：

$S = \sum_{i=1}^m r_i^2$S=i=1∑m​ri2​

梯度向量在方向上的分量：

$g_j=2\sum_{i=1}^m r_i\frac{\partial r_i}{\partial \beta_j}. （1）$gj​=2i=1∑m​ri​∂βj​∂ri​​.（1）

Hessian 矩阵的元素则直接在梯度向量的基础上求导：

$H_{jk}=2\sum_{i=1}^m \left(\frac{\partial r_i}{\partial \beta_j}\frac{\partial r_i}{\partial \beta_k}+r_i\frac{\partial^2 r_i}{\partial \beta_j \partial \beta_k} \right)$Hjk​=2i=1∑m​(∂βj​∂ri​​∂βk​∂ri​​+ri​∂βj​∂βk​∂2ri​​)

高斯牛顿法的一个小技巧是，将二次偏导省略，于是：

$H_{jk}\approx 2\sum_{i=1}^m J_{ij}J_{ik}（2）$Hjk​≈2i=1∑m​Jij​Jik​（2）

将(1)(2)改写成 矩阵相乘形式：

$\mathbf g=2\mathbf{J_r}^\top \mathbf{r}, \quad \mathbf{H} \approx 2 \mathbf{J_r}^\top \mathbf{J_r}.\,$g=2Jr​⊤r,H≈2Jr​⊤Jr​.

其中 Jr 是雅克比矩阵,r是残差$r_i$ri​ 组成的列向量。

代入牛顿法高维迭代方程的基本形式，得到高斯牛顿法迭代方程：

$\boldsymbol{\beta}^{(s+1)} = \boldsymbol\beta^{(s)}+\Delta;\quad \Delta = -\left( \mathbf{J_r}^\top \mathbf{J_r} \right)^{-1} \mathbf{J_r}^\top \mathbf{r}.$β(s+1)=β(s)+Δ;Δ=−(Jr​⊤Jr​)−1Jr​⊤r.

即

$\boldsymbol \beta^{(s+1)} = \boldsymbol \beta^{(s)} - \left(\mathbf{J_r}^\top \mathbf{J_r} \right)^{-1} \mathbf{ J_r} ^\top \mathbf{r}(\boldsymbol \beta^{(s)}) ， \mathbf{J_r} = \frac{\partial r_i (\boldsymbol \beta^{(s)})}{\partial \beta_j}$β(s+1)=β(s)−(Jr​⊤Jr​)−1Jr​⊤r(β(s))，Jr​=∂βj​∂ri​(β(s))​（wiki上的公式，我觉得分母写错了）

m ≥ n 是必要的，因为不然的话，$Jr^TJr$JrTJr肯定是不可逆的（利用rank(A’A) = rank(A)）

高斯牛顿法只能用于最小化平方和问题，但是优点是，不需要计算二阶导数。