深入浅出最优化(6) 最小二乘问题的特殊方法

1 基本数学表达

在前面3节中，我们使用了不同下降方法来求解同一个非线性最小二乘问题，但其实非线性最小二乘问题只是这些下降方法能够求解的问题当中的一个特例。接下来要介绍的方法，将是专门针对非线性最小二乘问题设计的，具有非常优良的相性。由于非线性最小二乘问题的应用广泛，这个方法的介绍也是尤为重要的。

下面来看一些基本的数学表达：

梯度： $g=\nabla f(x)$
黑森矩阵： $G=\nabla^2 f(x)$
雅可比矩阵： $J(x)=\left[\begin{matrix}\frac{\partial y_1}{\partial x_1}&...&\frac{\partial y_1}{\partial x_n}\\\vdots&\ddots&\vdots\\\frac{\partial y_m}{\partial x_1}&...&\frac{\partial y_m}{\partial x_n}\end{matrix}\right]$

雅可比矩阵作用：如果 $P$ 是 $R^n$ 中的一点， $F$ 在 $P$ 点可微分，那么在这一点的导数由 $J_F(P)$ 给出，在此情况下，由 $F(P)$ 描述的线性算子即接近点 $P$ 的 $F$ 的最优线性逼近： $F(x)\approx F(P)+J_F(P)(X-P)$
残差：表示实际观测值与估计值（拟合值）之间的差
多维无约束牛顿法： $x_{k+1}=x_k-G_k^{-1}g_k$

2 高斯-牛顿法

对于非线性最小二乘问题，假设有 $n$ 个待定系数（ $x\in R^n$ ）和 $m$ 组样本（ $m\geq n$ ，这是确定待定系数的必要条件），则残差函数： $r_i(x)=y_i-\tilde{f}(x;t_i),i=1,...,m$ 。

设 $r(x)=\left[\begin{matrix}r_1(x)\\\vdots\\r_m(x)\end{matrix}\right]$

则非线性最小二乘问题的目标函数为 $min f(x)=min\{\frac{1}{2}\displaystyle\sum_{i=1}^m[r_i(x)]^2\}=min\{\frac{1}{2}r(x)^Tr(x)\}$

同时我们也可以得到 $r(x)$ 的雅可比矩阵：

$J_{r}(x)=\left[\begin{matrix}\frac{\partial r_1}{\partial x_1}&...&\frac{\partial r_1}{\partial x_n}\\\vdots&\ddots&\vdots\\\frac{\partial r_m}{\partial x_1}&...&\frac{\partial r_m}{\partial x_n}\end{matrix}\right]=\left[\begin{matrix}\nabla r_1^T\\\vdots\\\nabla r_m^T\end{matrix}\right]$

用牛顿法来优化该函数，对 $f$ 求梯度得 $g=\nabla f=2\displaystyle\sum_{i=1}^mr_i\frac{\partial r_i}{\partial x_j}$ ，再对所有元分别求导得 $G_{jk}=2\displaystyle\sum_{i=1}^m(\frac{\partial r_i}{\partial x_j}\frac{\partial r_i}{\partial x_k}+r_i\frac{\partial^2r_i}{\partial x_j\partial x_k})$ 。由于二阶导数非常小，一般可以忽略，则 $g(x)=\displaystyle\sum_{i=1}^mr_i(x)\nabla r_i(x)=J_r(x)^Tr(x),G(x)=J_r(x)^TJ_r(x)$ ，相应的迭代式被改写为 $x_{k+1}=x_k+d_k,d_k=-(J_r^TJ_r)^{-1}J_r^Tr$ 。

假如采用阻尼牛顿法来搜索，则步骤如下所示：

给定解的初始估计 $x_0$ ， $\epsilon>0,k=0$
如果 $x_k$ 满足终止准则，则停止迭代
解方程组 $J_k^TJ_kd_k=-J_k^Tr_k$
计算步长 $\alpha_k$
$x_k=x_k+\alpha_kd_k$ ， $k=k+1$ ，转步2

高斯-牛顿法具有和阻尼牛顿法相似的收敛特性，且通过将黑森矩阵用 $J_r(x)^TJ_r(x)$ 近似，在黑森矩阵非奇异时，可以保证产生方向是下降方向
对于线性最小二次问题一步收敛到最优解

3 LMF方法

高斯-牛顿法可以说是牛顿法的一个改进方案，但并非是万能的。若雅可比矩阵奇异，则d无解，这说明单纯使用雅可比矩阵不能很好地近似原函数的梯度与黑森矩阵了。不解方程组 $J_k^TJ_kd_k=-J_k^Tr_k$ ，而是考虑最小二乘问题 $(J_k^TJ_k+v_kI)d=-J_k^Tr_k$ ，即使用一个因子 $v_kI$ 来去奇异。这个因子通过信赖域方法（LM法）来确定，根据雅可比矩阵近似的函数变化率与真实函数变化率的差别来增大或减小。

根据我们之前的近似 $g(x)=\displaystyle\sum_{i=1}^mr_i(x)\nabla r_i(x)=J_r(x)^Tr(x),G(x)=J_r(x)^TJ_r(x)$ 。将 $minf(x_k+d)$ 近似为 $minq_k(d)=\frac{1}{2}||J_kd+r_k||^2$ ，其中 $q_k(d)=\frac{1}{2}(J_kd+r_k)^T(J_kd+r_k)=\frac{1}{2}d^TJ_k^TJ_kd+d^T(J_k^Tr_k)+\frac{1}{2}r_k^Tr_k$ 。

$f(x)$ 的实际减少量为 $\Delta f_k=f(x_k)-f(x_k+d_k)$ ， $q(d)$ 的减少量为 $\Delta q_k=q(0)-q(d_k)$ 。

定义 $\gamma_k=\frac{\Delta f_k}{\Delta q_k}$ ，给出如下准则：

若 $\gamma_k<0.25$ ，说明不能很好地近似原函数，将 $v_k$ 放大为原来4倍后迭代下一步
若 $\gamma_k>0.75$ ，说明可以较好近似原函数，将 $v_k$ 缩小为原来 $\frac{1}{2}$ 后迭代下一步
若 $\gamma_k<0$ ，说明 $J_k^TJ_k+v_kI$ 可能负定，会导致产生的方向不是下降方向。此时不迭代下一步，而仅调整 $v_k$

LMF方法步骤如下：

给定 $x_0\in R^n$ ， $v_0>0,\epsilon>0,k=0$
若满足终止条件，停止迭代
求解 $(J_k^TJ_k+v_kI)d=-J_k^Tr_k$ ，得到 $d_k$
计算 $\gamma_k$
若 $\gamma_k<0.25$ ，则 $v_{k+1}=4v_k$ ，若 $\gamma>0.75$ ，则 $v_{k+1}=v_k/2$ ，否则 $v_{k+1}=v_k$
若 $\gamma\leq0$ ，则 $x_{k+1}=x_k$ ，否则 $x_{k+1}=x_k+d_k$ ，令 $k=k+1$ ，转步2

4 实战测试

对于本文集第一篇文章深入浅出最优化(1) 最优化问题概念与基本知识中提出的最小二乘问题， $x_1,x_2,x_3,x_4$ 的初值均在 $[-2,2]$ 的范围内随机生成，总共生成100组起点。统计迭代成功（在1000步内得到最优解且单次步长搜索迭代次数不超过1000次）的样本的平均迭代步数、平均迭代时间和得到的最优解及残差平方和最小值。

平均迭代步数	平均迭代时间	最优解	残差平方和最小值
7.2	0.06s	$x_1=0.1925~x_2=0.1959~x_3=0.1234~ x_4=0.1383$	$1.5378\times10^{-4}$

5 无约束最优化方法比较

在这里插入图片描述

7种方法分别是最速梯度下降法、牛顿-梯度下降混合法、BFGS拟牛顿法、DFP拟牛顿法、共轭梯度PRP+法、共轭梯度FR法和LMF法，对实战测试的迭代步数、计算时间和残差平方和最小值进行均值归一化后进行对比。

可以发现BFGS法、PRP+法和LMF法均有不俗的表现。其中LMF法是专门针对最小二乘问题的方法，在该问题上表现最为突出。

PRP+法作为共轭梯度法，在储存容量上相较BFGS法优势更大。因此虽然对于待定系数只有4个的测试问题，两者的差别不是非常明显，但遇到大规模最优化问题时，更推荐采用PRP+法。

代码实现

本博客所有代码在https://github.com/HarmoniaLeo/optimization-in-a-nutshell开源，如果帮助到你，请点个star，谢谢这对我真的很重要！

你可以在上面的GitHub链接或本文集的第一篇文章深入浅出最优化(1) 最优化问题概念与基本知识中找到Function.py和lagb.py

from Function import Function   #定义法求导工具
import numpy as np
from scipy import linalg
from lagb import *  #线性代数工具库

n=4 #待定系数数
y=np.array([0.1957,0.1947,0.1735,0.1600,0.0844,0.0627,0.0456,0.0342,0.0323,0.0235,0.0246])
t=np.array([4,2,1,0.5,0.25,0.1670,0.1250,0.1000,0.0833,0.0714,0.0625])
#样本列表


def myFunc(x):  #目标函数（残差平方和的0.5倍）
    return 1/2*np.sum(np.square(y-x[0]*(t**2+x[1]*t)/(t**2+x[2]*t+x[3])))

def r(x):   #残差函数
    yi = y[int(x[n])]
    ti = t[int(x[n])]
    fxt = x[0]*(ti**2+x[1]*ti)/(ti**2+x[2]*ti+x[3]) #原函数
    return yi - fxt

def Jaccobi(x):
    tar=Function(r)
    mat=np.empty((y.shape[0],n))
    for j in range(0,y.shape[0]):
        for i in range(0,n):
            mat[j][i]=tar.part(i,np.append(x,j))
    return mat

def q(J,R,d):
    return 0.5*muldot(turn(d),turn(J),J,d)+dot(turn(d),dot(turn(J),R))+0.5*dot(turn(R),R)

e=0.001
k=0
v=np.eye(n)
tar=Function(myFunc)
x=np.zeros(n)   #初值点
while tar.norm(np.concatenate((x,t)))>e:
    J=Jaccobi(x)
    A=dot(turn(J),J)+v
    R=np.empty(0)
    for i in range(0,t.shape[0]):
        R=np.append(R,r(np.append(x,i)))
    b=-dot(turn(J),R)
    d=linalg.solve(A,b)
    gamma=(tar.value(x)-tar.value(x+d))/(q(J,R,np.zeros(n))-q(J,R,d))
    if gamma<0.25:
        v*=4
    elif gamma>0.75:
        v/=2
    if gamma>0:
        x+=d
    k+=1
    print(k)
print(x)