一、最小二乘法的引出

我们考虑下面一个方程组的求解：
x1+x2=5x1−x2=3x1+x2=4x_1+x_2=5 \\ x_1-x_2=3 \\ x_1+x_2=4 \\ x1+x2=5x1−x2=3x1+x2=4
可以写成形如Ax=bAx=bAx=b的矩阵形式：
[111−111][x1x2]=[534]\begin{bmatrix} 1 & 1\\ 1 & -1\\ 1 & 1 \end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ \end{bmatrix} =\begin{bmatrix} 5\\ 3\\ 4 \end{bmatrix} ⎣⎡1111−11⎦⎤[x1x2]=⎣⎡534⎦⎤
由线性代数的知识可知，当方程个数m大于未知两个数n的时候，方程通常无解。但是对于Ax=bAx=bAx=b这个方程而言，我们仍然可以通过一些方法得到一个最接近解析解的一个向量，这个方法就是最小二乘法。
最小二乘法的抽象形式如下所示：
目标函数=∑(观测值−理论值)2目标函数=\sum (观测值-理论值)^2 目标函数=∑(观测值−理论值)2
观测值指的就是我们的多组样本，而理论值就是我们假设的拟合函数。目标函数就是我们常说的损失函数。我们的目标是得到使目标函最小化时的拟合函数模型。

我们举一个简单的例子，假设有m个测试样本：(xi,yi)(i=1,2,3,......m)(x_i,y_i)(i=1,2,3,......m)(xi,yi)(i=1,2,3,......m)，我们需要找到一个n维的变量x∗∈Rx^* \in \mathbb{R}x∗∈R,使得损失函数F(x)F(x)F(x)取得局部最小值：
F(x)=12∑i=1m(fi(x))2F(x)=\frac{1}{2} \sum_{i=1}^m (f_i(x))^2 F(x)=21i=1∑m(fi(x))2
其中fif_ifi是残差函数，即测量值和预测值之间的差，且有m>nm>nm>n。而局部最小值指任意的∣∣x−x∗∣∣<δ||x-x^*|| < \delta∣∣x−x∗∣∣<δ 有F(x∗)≤F(x)F(x^*) \leq F(x)F(x∗)≤F(x)。

二、线性最小二乘法

1.线性最小二乘的描述

对于mmm个方程求解nnn个未知数，有三种情况：

m=nm=nm=n且A为非奇异，则有唯一解，x=A−1bx=A^{-1}bx=A−1b
当m>nm>nm>n，约束的个数大于未知数的个数，称为超定方程（overdetermined）
m<nm<nm<n，负定/欠定问题(underdetermind)

通常我们遇到的问题都是超定问题，此时Ax=bAx=bAx=b的解是不存在的，从而转向解最小二乘问题：
J(x)=min∣∣Ax−b∣∣22J(x)=min||Ax-b||^2_2 J(x)=min∣∣Ax−b∣∣22
J(x)J(x)J(x)为凸函数，我们令一阶导数为0，可以得到：
x=(ATA)(−1)ATbx=(A^TA)^{(-1)}A^Tb x=(ATA)(−1)ATb

2.线性最小二乘特殊情况的求解

这时候我们考虑齐次方程Ax=0Ax=0Ax=0，当A行数大于列数的特殊情况：
min∣∣Ax∣∣min||Ax|| min∣∣Ax∣∣
此时的最小二乘解为ATAA^TAATA最小特征值对应的特征向量。我们可以通过以下方式进行求解；

3.线性最小二乘一般情况的求解

对于一个线性方程求解问题而言，最简单粗暴的方式是对ATAA^TAATA求逆，但是这种方式计算量极大，非常耗时，所以我们考虑用矩阵分解的方式加速线性最小二乘问题的求解。

矩阵分解是将矩阵拆解为数个矩阵的乘积，可分为三角分解、满秩分解、QR分解、Jordan分解、SVD（奇异值）分解等，常见的有三种：1)QR 分解法，2)Cholesky分解法，3)SVD分解法

LU三角分解：三角分解法是将原正方 (square) 矩阵分解成一个上三角形矩阵和一个下三角形矩阵.
QR分解：QR分解法是将矩阵分解成一个正规正交矩阵与上三角形矩阵,所以称为QR分解法。
奇异值分解：[U,S,V]=svd(A)，其中U和V分别代表两个正交矩阵，而S代表一对角矩阵。SVD是最可靠的分解法，但是它比QR 分解法更加耗时。
LLT分解：Cholesky 分解是把一个对称正定的矩阵表示成一个下三角矩阵L和其转置的乘积的分解。它要求矩阵的所有特征值必须大于零，故分解的下三角的对角元也是大于零的（LU三角分解法的变形）。
LDLT分解法：其中L为下三角形单位矩阵，D为对角矩阵，L^T为L的转置矩阵。LDLT分解法实际上是Cholesky分解法的改进，因为Cholesky分解法虽然不需要选主元，但其运算过程中涉及到开方问题，而LDLT分解法则避免了这一问题。

以下是eigen库中提供的矩阵分解求解线性最小二乘问题的方法，我们可以根据实际需求做出选择。

三、非线性最小二乘法

1.非线性最小二乘的描述

考虑一个非线性最小二乘模型如下：
F(x)=12∑i=1m(fi(x))2F(x)=\frac{1}{2}\sum^m_{i=1} (f_i(x))^2 F(x)=21i=1∑m(fi(x))2
我们假设损失函数F(x)F(x)F(x)是平滑可导的，所以我们在二阶泰勒展开有：
F(x+Δx)=F(x)+JΔx+12ΔxTHΔx+O(∣∣Δ∣∣3)F(x+\Delta x)=F(x)+J\Delta x+\frac{1}{2}\Delta x^TH\Delta x+O(||\Delta||^3) F(x+Δx)=F(x)+JΔx+21ΔxTHΔx+O(∣∣Δ∣∣3)
其中JJJ和HHH分别为损失函数FFF对变量xxx的一阶导和二阶导矩阵。

2.非线性最小二乘的求解

我们求解的思路是寻找一个下降方向使损失函数随xxx的迭代逐渐减小，直到xxx收敛到x∗x^*x∗时有:
F(x∗)≤F(x)x∈∀F(x^*) \leq F(x) \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ x \in \forall F(x∗)≤F(x) x∈∀
所以我们需要寻找下降方向的单位向量ddd和下降步长α\alphaα。假设我们的α\alphaα足够小，我们对损失函数F(x)F(x)F(x)进行一阶泰勒展开：
F(x+αd)≈F(x)+αJdF(x+\alpha d) \approx F(x) +\alpha Jd F(x+αd)≈F(x)+αJd
我们的下降方向只要满足：
Jd<0Jd < 0 Jd<0
而下降的步长可以通过线性搜索的方式寻找：
α∗=argminα>0F(x+αd)\alpha^*=argmin_{\alpha>0}{F(x+\alpha d)} α∗=argminα>0F(x+αd)

a. 最速下降法

从下降的方向条件可知：Jd=∣∣J∣∣cosθJd=||J||cos\thetaJd=∣∣J∣∣cosθ,θ\thetaθ表示方向和梯度方向的夹角，当θ=π\theta=\piθ=π时：
d=−JT∣∣J∣∣d=\frac{-J^T}{||J||} d=∣∣J∣∣−JT
即梯度的负方向为最速下降方向，此时的缺点为最优值附件震荡，收敛慢。

b.牛顿法

在局部最优点x∗x^*x∗附近时，若x+Δxx+\Delta xx+Δx是最优解，则损失函数对Δx\Delta xΔx的导数等于0：
∂∂x(F(x)+JΔx+12ΔTHΔx)=JT+HΔx=0\frac{\partial }{\partial x}(F(x)+J\Delta x+\frac{1}{2}\Delta^TH\Delta x)=J^T+H\Delta x=0 ∂x∂(F(x)+JΔx+21ΔTHΔx)=JT+HΔx=0
可以得到 : Δx=−H−1JT\Delta x=-H^{-1}J^TΔx=−H−1JT，显然这种方法缺点非常明显，需要计算二阶导矩阵，计算比较复杂

c.高斯牛顿法

对残差函数f(x)f(x)f(x)一阶泰勒展开：
f(x+Δx)≈ζ(Δx)=f(x)+JΔxf(x+\Delta x) \approx \zeta(\Delta x)=f(x)+J\Delta x f(x+Δx)≈ζ(Δx)=f(x)+JΔx
需要注意的是这里的JJJ代表的是残差函数f(x)f(x)f(x)的雅可比矩阵，我们带入损失函数有：
F(x+Δx)≈L(Δx)=12ζ(Δx)Tζ(Δx)=F(x)+ΔxTJTf+12ΔxTJTJΔxF(x+\Delta x) \approx L(\Delta x) = \frac{1}{2}\zeta(\Delta x)^T\zeta(\Delta x) \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ =F(x)+\Delta x^TJ^Tf+\frac{1}{2}\Delta x^TJ^TJ\Delta x F(x+Δx)≈L(Δx)=21ζ(Δx)Tζ(Δx) =F(x)+ΔxTJTf+21ΔxTJTJΔx
令上述公式的一阶导为0，我们可以得到：
(JTJ)Δxgn=−JTf(J^TJ) \Delta x_{gn}=-J^Tf (JTJ)Δxgn=−JTf
常表述为：HΔxgn=bH\Delta x_{gn}=bHΔxgn=b

d.LM法

LM法是高斯牛顿法和最速下降法的一种结合：
(JTJ+uI)Δxlm=−JTfwithu≥0(J^TJ+uI)\Delta x_{lm} =-J^Tf \ \ \ \ \ \ \ \ \ with \ \ \ u\geq0 (JTJ+uI)Δxlm=−JTf with u≥0
相当于在求解过程中引入阻尼因子，其作用为：

u>0u>0u>0保证(JTJ+uI)(J^TJ+uI)(JTJ+uI)正定，迭代朝着下降的方向进行。
uuu比较大时，Δxlm=−1uJTf=−1uF′(x)T\Delta x_{lm}=-\frac{1}{u}J^Tf=-\frac{1}{u}F'(x)^TΔxlm=−u1JTf=−u1F′(x)T,接近最速下降法。
uuu比较小时，Δxlm≈Δxgn\Delta x_{lm} \approx \Delta x_{gn}Δxlm≈Δxgn,接近高斯牛顿法。

一文详解线性最小二乘与非线性最小二乘相关推荐

详解线性结构和非线性结构
一.线性结构 1.线性结构作为最常用的数据结构,其特点是数据元素之间存在一对一的线性关系. 2.线性结构有两种不同的存储结构,即顺序存储结构和链式存储结构.顺序存储的线性表称为顺序表,顺序表中的存储元 ...
一文详解编程中的随机数
一文详解编程中的随机数随机数的类型真随机数生成器 TRNG - True Random Number Generator 伪随机数生成器 PRNG - Pseudo Random Number G ...
一文详解基于测距的空间定位算法
一文详解基于测距的空间定位算法文章目录一文详解基于测距的空间定位算法 0 定位算法分类 0.1 基于测距与非基于测距的定位算法 0.2 集中式与分布式定位算法 0.3 绝对与相对定位算法 0.4 ...
一文详解JavaBean 看这篇就够了
一文详解JavaBean 看这篇就够了 JavaBean的历史渊源 JavaBean的定义(通俗版) JavaBean应用 < jsp:useBean > < jsp:getProp ...
【卷积神经网络结构专题】一文详解AlexNet(附代码实现)
关注上方"深度学习技术前沿",选择"星标公众号", 资源干货,第一时间送达! [导读]本文是卷积神经网络结构系列专题第二篇文章,前面我们已经介绍了第一个真正意义 ...
一文详解 YOLO 2 与 YOLO 9000 目标检测系统
一文详解 YOLO 2 与 YOLO 9000 目标检测系统 from 雷锋网雷锋网 AI 科技评论按:YOLO 是 Joseph Redmon 和 Ali Farhadi 等人于 2015 年提出 ...
一文详解决策树算法模型
AI有道一个有情怀的公众号上文我们主要介绍了Adaptive Boosting.AdaBoost演算法通过调整每笔资料的权重,得到不同的hypotheses,然后将不同的hypothesis乘以不 ...
「软件项目管理」一文详解软件配置管理计划
一文详解软件配置管理计划前言一.配置管理概述 1. 配置管理(SCM)定义 2. 软件配置项目(SCI) 3. 基线 4. 软件配置控制委员会(SCCB) 二.软件配置管理过程 1. 管理过程 2 ...
「软件项目管理」一文详解软件项目质量计划
一文详解软件项目质量计划

一文详解线性最小二乘与非线性最小二乘