最小二乘法

作者：@houkai
本文为作者原创，转载请注明出处：https://www.cnblogs.com/houkai/p/6369870.html

最小二乘法用于在回归分析中，计算overdetemined system的近似解。overdetemined system是指等式数目多余未知量的问题。其思想是：最小化各个等式的误差的平方和。

最常见应用于数据拟合（如下），此时优化目标是最小化平方残差，残差即为观测值和模型计算的拟合值的差。

按照残差项是够是线性的，最小二乘法分为两类：线性最小二乘、非线性最小二乘。

线性最小二乘常见于回归分析，它有封闭的解（解具有等式形式，如 $a*x^{2}+b*x+c=0$ ）;非线性最小二乘通常是迭代寻优，每一次迭代中,将问题近似为线性。

问题描述

最小二乘问题的数学表示：

$\left\{ \begin{array}{ll} S=\sum_{i=1}^{n}r_i^{2}\\ r_{i}=y_{i}-f(x_{i},\beta)) \end{array} \right.$

其中 $f(x_{i},\beta)$ 为模型, 如直线拟合中: $f(x,\beta)=\beta_{0}+\beta_{1}*x$

$\beta$ 为模型参数向量,观察值为 $(x_{i},y_{i})\,i=1...n$

算法求解

$min(S)$ 为问题的目标函数，最小化平方和，即寻找梯度为0的解，设有m个参数，则有：

$\frac{\partial{S}}{\partial{\beta_{j}}}= 2\sum_{i=1}^{n}{r_{i}*\frac{\partial{r_{i}}}{\partial{\beta_{j}}}}=0$ ，其中 $j=1...m$ .

带入 $r_{i}=y_{i}-f(x_{i},\beta)$ , 得到:

$-2\sum_{i=1}^{n}{r_{i}* \frac{ \partial{f(x_{i},\beta)} }{ \partial{\beta_{j}} } }=0$ , $j=1...m$ （1）

线性最小二乘法

模型是参数的线性组合, $f(x,\beta)=\sum_{j=1}^{m}{\beta_{j}*\phi_{j}(x)}$ , 其中 $\phi_j$ 是x的函数（它是一个确定值）。所以：

$X_{ij}=\frac{\partial{f(x_{i},\beta)}}{\partial{\beta_{j}}}=\phi_{j}(x)$

将 $r_{i}=y_{i}-\sum_{j=1}^{m}{\beta_{j}*X_{ij}}$ 带入(1)式子，得到矩阵形式：

$X^{T}X\hat{\beta}=X^{T}y$ 其中 $X^{T}X$ 是正定矩阵,所有:

$\hat{\beta}=(X^{T}X)^{-1}X^{T}y$

非线性最小二乘法

并不存在上述类似的封闭解,而是采用数字算法来寻优 $\hat{\beta}$ , 为参数设置初值,然后进行迭代调优,直至收敛。

$\beta_{j}^{k+1}=\beta_{j}^{k}+\Delta\beta_{j}$ , k为迭代的次数

$\Delta\beta_{j}$ 称为 shift vector，位移向量。

每一次迭代，模型可以用关于 $\beta^{k}$ 的Taylor一阶展开进行线性近似：

$f(x_{i},\hat{\beta})=f(x_{i},\beta^{k})+\sum_{j=0}^{m}{ \frac{\partial{f(x_{i},\beta)}} {\partial{\beta_{j}}}*(\beta_{j}-\beta_{j}^{k}) }\\ =f(x_{i},\beta^{k})+\sum_{j=0}^{m}{J_{ij}*\Delta\beta_{j}}$

J是确定数值的Jacobian矩阵（独立于y和参数β）。

由此得：

$r_{i}=y_{i}-f(x_{i},\beta^{k})-\sum_{j=1}^{m}{J_{ij}\Delta\beta{j}}= \Delta{y_{i}}-\sum_{j=1}^{m}{J_{ij}\Delta\beta_{j}}$

最小化上式,梯度为0,得：

$-2\sum_{i=1}^{n}{r_{i}*J_{ij}}=0$ (2)

将 $r_{i}$ 带入(2)式，得到结果的矩阵形式：

$J^{T}J\Delta{\beta}=J^{T}\Delta{y}$

所以: $\Delta{\beta}=(J^{T}J)^{-1}J^{T}\Delta{y}$

这就是Gauss–Newton algorithm的等式。

The Gauss–Newton algorithm is used to solve non-linear least squares problems. It is a modification of Newton's method for finding a minimum of a function. Unlike Newton's method, the Gauss–Newton algorithm can only be used to minimize a sum of squared function values, but it has the advantage that second derivatives, which can be challenging to compute, are not required.

Gauss-Newton算法不同于Newton方法，它只用于最小化平方和问题，但是它具有类似二阶导的优势。

参考:

https://en.wikipedia.org/wiki/Least_squares

https://en.wikipedia.org/wiki/Gauss%E2%80%93Newton_algorithm

进一步扩展

最速下降法

在基本迭代公式中 $X_{k+1}=X_{k}+t_{k}P_{k}$ 中,每次迭代搜索方向 $P_{k}$ 取为目标函数f(x)的负梯度方向。

$P_{k}=-\nabla{f(x_{k})}$

而 $t_{k}$ 为最优步长，由此确定的算法称为最速下降法。

牛顿法

设最优问题为： $minf(X)$ ,其中f二阶可到，Hesse矩阵 $\nabla^{2}{f(X)}$ 正定。

k次迭代得到 $x_{k}$ , 将f(x)在 $x_{k}$ 处展开为二阶泰勒公式，得：

$f(X)≈Q(X)=f(X_{k})+\nabla{f(X_{k})^{T}}(X-X_{k})+\frac{1}{2}(X-X_{k})^{T}\nabla^{2}f(X_{k})(X-X_{k})$

显然Q(X)是正定二次函数，所以Q(X)是凸函数，且存在唯一局部极小值，得

$X-X_{k}=\frac{\nabla{f(X_{k})}} {\nabla^{2}f(X_{k})}$ $X-X_{k}=-\frac{\nabla{f(X_{k})}}{\nabla^{2}f(X_{k})}$

所以: $X=X_{k}- }[\nabla^{2}f(X_{k})]^{-1}{\nabla{f(X_{k})}$

相应有:

$P_{k}= -[\nabla^{2}f(X_{k})]^{-1}{\nabla{f(X_{k})}$

$X_{k+1}=X_{k}+t_{k}P_{k}$

对于二次函数, 一次迭代即可得到最优值。

Levenberg-Marquart算法

LM算法实现了在最速下降法和Inverse-Hessian算法之间进行平稳的变化。是求解最小二乘法最常用的方法。

将 $J^{T}J\Delta{\beta}=J^{T}\Delta{y}$

转化为： $(1+\lambda)*J^{T}J\Delta\beta=J^{T}\Delta{y}$

在L-M算法中，每一次迭代都是寻找一个合适的 $\lambda$ 值。算法开始时，通常取 $\lambda=10^{-3}$ ，若结算后的解 $\beta$ 导致误差减少，则接受的当前值，并在下一次迭代中以 $\lambda/10$ 代替 $\lambda$ 。若解导致误差的增大，则以 $10\lambda$ 代替 $\lambda$ 并重新求解增量方程。这个过程一直到求出一个使误差下降的 $\Delta{\beta}$ 为止，构成L-M算法的一次迭代。