非线性最小二乘

转载来自z​​​​​​​非线性最小二乘 - 知乎 (zhihu.com)，如有问题告知删除

当我们用一个模型 ϕ(t) 来描述现实中的一系列数据时，模型的预测结果与实际的测量结果总会存在一定偏差，这一偏差就称为残差。非线性最小二乘的目的就是，调整模型的参数，使得总的残差最小。

 

对于图上y是实际值， ϕ(t) 是计算得到的值，首先，对于每个测量结果， yj−ϕ(x;tj) 表示测量值和预测值之差，也就是残差。残差当然是越接近0越好。但是，为什么我们最小化残差的平方和呢？为什么不是绝对值之和？或者四次方之和？

最直接的原因是平方操作不显著增加多项式次数，且处处可导。相比之下，绝对值操作在原点不可导，四次方阶数太高，不易求解。

除此之外，平方和还有更深层次上的含义。

最大似然估计

通常，我们可以认为多个测量结果之间是独立同分布的，也就是某次测量不影响其它测量的结果。另一方面，假设我们设计的模型的确能够很好地描述真实的物理事件，那么，任意一次测量的误差就应该是一个随机误差，不妨认为该误差满足如下的高斯分布

 我们都知道，高斯分布的对称轴位置概率最大。恰好， gσ(ϵ) 的对称轴在原点，而我们期望的优化方向正是朝向原点。于是，最小化残差就可以等价于最大化所有测量的联合概率。由于测量是独立同分布的，所以可以直接把每个测量的 gσ(ϵj) 乘起来，得到

最大化上面式子得到：

 

由于 m 和 σ 都是定值，因此最大化 p(y;x,σ) 等价于最小化 ，正好是式(10.7)。

所以，无论从理论上还是实际操作上，残差的平方和都是最佳的目标函数。

接下来介绍几种求解非线性最小二乘的方法

1.高斯-牛顿法

其中， rj 是第 j 个残差。在该定义下，目标函数的一阶导和二阶导可以用残差的雅克比矩阵和海森矩阵来表示

 

其中， J(x) 是残差的雅克比矩阵， ∇2r(x) 是残差的海森矩阵。

牛顿法的每次迭代需要解线性方程 

将式(10.4)和(10.5)代入该线性方程，同时，舍弃式(10.5)中的第二项，可以得到

这一方程比标准的牛顿法要简便许多，因为不需要事先计算海森矩阵，只需要残差和雅克比矩阵即可。而且，这样做是有道理的，舍弃的二阶项通常很小，因为残差本身就很小，且越接近极小值点残差越小，因此不会影响整个算法的最终精度。

2.列文伯格-马夸尔特法 

这个算法还没好好学过，先跳过

这个方法还涉及信赖域法，下次一起看

求解线性方程

 

求解线性方程的关键在于避免对系数矩阵求逆，因为求逆的计算量太大。

最常用的线性方程解法是Cholesky分解，将 JTJ 分解为上三角矩阵与其转置的乘积，即

此时，线性方程就可以逐行求解，计算量大幅降低。但是，Cholesky分解的稳定性受条件数影响， JTJ 的条件数是 J 条件数的平方，这在某些不稳定的情况下会使得分解结果不可用。

另一种方法是QR分解，它的稳定性只受 J 条件数的影响，通常比Cholesky分解的精度更高。

第三种方法，也是精度最高的方法是SVD分解。即使雅克比矩阵 J 是病态的，只需要将奇异值为0对应的部分忽略，仍然可以求得稳定的解。

日后再好好看看这三种分解方法

正交距离

回顾非线性最小二乘问题的定义，可以发现残差仅仅指观测值与模型预测值之差。但实际上，在时间维度也会存在残差。我们可以通过作图来明显看出这二者的差别。对于同一个模型和同样的观测数据，我们定义的残差可以用下图中虚线线段表示

显然，这些线段都与y轴平行，表示它们只处于观测维度上。如果同时考虑观测和时间两个维度，我们可以像如下的方式定义残差 

这些残差既不平行于y轴，也不平行于t轴，而是与模型曲线垂直，表示了模型上的点到观测点的最短距离。

由于这些残差与模型垂直，因此这种方法称为正交距离回归（Orthogonal Distance Regression）。此时，模型不再表示为 ϕ(x;t) ，而是表示为 ϕ(x;t+δ) 。其中， δ 为时间维度上的残差。目标函数修改为同时最小化观测残差和时间残差

 其中， ϵj 仍然是观测残差 yj−ϕ(x;tj+δj) 。后续的求解方法和普通的非线性最小二乘是一样的，只不过优化变量增加了 m 个，计算量略微增大。