四参数拟合算法之高斯牛顿法

介绍

 前面提到了牛顿法，那其实相当于求根的算法。跟一般最小二乘法的区别是，它并没有显示的最小二乘目标式子。
$$F(A,B,C,D,x_i)=0$$

 下面提到的高斯牛顿法，则要正式引入最小二乘法的目标式子。首先对牛顿法做一次更深入的展开。

牛顿法

 牛顿法在用于一元方程求根的时候，只需要做一阶泰勒展开，这个时候，用到的是迭代点的导数信息找到下一个迭代点。在多维的情况下面，则用到了梯度信息。在求解的问题的时候，假如初始点选得不好，可能会遇到梯度消失的问题，最后会导致收敛失败。但这好像是优化算法的一个通病，并非牛顿法的专利。言归正传，下面解释一下牛顿法在四参数数据拟合中应用。
对于给定的m组数据 $(x_i,y_i)$,
$$\begin{gathered} w=[A,B,C,D{]}^T \\ f(w;x_i)=D+\frac{A-D}{1+(x_i/C{)}^B} \\ {r}_i(w)=y_i-f(w;x_i) \end{gathered}$$
$r_i(w)\in \mathbb{R}$称为剩余值，或者叫残差，偏差，是一个标量，$r(w)\in {\mathbb{R}}^{m\times 1}$表示剩余函数，是一个向量。$f$为待拟合的公式
求解的目标是使得残差平方最小，如下
$$minF(w)=\frac{1}{2}r(w{)}^{T}r(w)=\frac{1}{2}\sum _{i=1}^m||y_i-f(w;x_i)|{|}^2$$
牛顿法的第一步是求得$F{}^{\prime }(w)$，令$F{}^{\prime }(w)=0$,再对其做泰勒展开
$$\begin{gathered} F{}^{\prime }(w+\Delta w)\approx F{}^{\prime }(w)+F{}^{\prime \prime }(w)\Delta w=0 \\ \Rightarrow \Delta w=-(F{}^{\prime \prime }(w){)}^{-1}F{}^{\prime }(w) \end{gathered}$$
总共需要对F做二次求导，这里需要引入jacob矩阵
w 关于r的雅可比矩阵如下
$$J=\left[\begin{array}{ccc}\frac{\partial r_1}{\partial A}& \cdots & \frac{\partial r_1}{\partial D}\\ ⋮& \ddots & ⋮\\ \frac{\partial r_m}{\partial A}& \cdots & \frac{\partial r_m}{\partial D}\end{array}\right]=\left[\begin{array}{c}{r}_1^{{}^{\prime }}(w)\\ ⋮\\ r_m^{{}^{\prime }}(w)\end{array}\right]$$

F的梯度信息g(w)
$$g(w)=F^{{}^{\prime }}(w)=\sum _{i=1}^m{r}_i(w)r_i^{{}^{\prime }}(w)=J^{T}r(w)$$
$r_i(w)\in \mathbb{R}$为标量，$r_i(w{)}^{{}^{\prime }}\in {\mathbb{R}}^{m\times 1}$为偏导，或者说梯度，是一个向量，最后得到，$g(w)\in {\mathbb{R}}^{m\times 1}$
F的海赛阵为
$$G(w)=\sum _{i=1}^m(r_i^{{}^{\prime }}(w)r_i^{{}^{\prime }}(w{)}^T+r_i(w)r_i^{{}^{\prime \prime }}(w))=J^{T}J+S(w)$$
$$\Delta w=-(F{}^{\prime \prime }(w){)}^{-1}F{}^{\prime }(w)=-(J^{T}J+S(w){)}^{-1}{J}^{T}r(w)$$
$S(w)\in {\mathbb{R}}^{n\times n}$通常运算量巨大，这部分常用两种处理，一种是忽略，即小剩余算法，另外一种是近似，大剩余算法。
高斯牛顿法忽略了其高阶项S(x),相当于目标式子的一阶泰勒展开
$$\Delta w=-(F{}^{\prime \prime }(w){)}^{-1}F{}^{\prime }(w)=-(J^{T}J{)}^{-1}{J}^{T}r(w)$$
这个做法其实等同于之前的牛顿法。当残量$r(w^{\ast })$比较小时，这种方法比较有效，反之，算法就比较难收敛。另外，从上公式里面可以看到，$\Delta w$的计算要求$J$为列满秩，否则求得的结果可能不稳定，导致算法难以收敛。 为了克服这个问题，其中一个手段是改善矩阵$J^{T}J$的变态程度。常用的处理就是在方程中加入扰动，也就是岭回归做的事情。从另外一个角度来说，也是加入$\Delta w$的2范式惩罚。为什么这么做呢，因为对于泰勒展开而言，无疑$\Delta w$越大，展开式越不准确。由此，便得到了Levenberg-Marquardt算法。

Matlab Code

>> [A,B,C,D] =-10 36.9 1974.7 431.0

下面的代码和之前代码基本一致，只不过去掉了matlab 的symbolic运算，速度快的飞起

clear;
load census;
x1 = cdate ;
y1 = pop ;
lamda = 0.3;
%init a,b,c,d
d = max(y1)+1;
a = min(y1)-0.1;
y2 = log((y1-d)./(a-y1));
x2 = log(x1);
[curve2,gof2] = fit(x2,y2, 'poly1');
b = -curve2.p1;
c = exp(curve2.p2/b);
%newton
for k = 1:1:1000
    fy = [];
    JacobiMatrix = ones(length(x1),4);
    w=[a,b,c,d];
    [res,R,fit] = evaluateFit(y1,x1,w);
    if R>0.997
        return;
    end
    for i = 1:1:length(x1)
        fy(i,:) =  d+(a-d)/(1+(x1(i)/c)^b);
        JacobiMatrix(i,1) = calc_pA(x1(i),a,b,c,d);
        JacobiMatrix(i,2) = calc_pB(x1(i),a,b,c,d);
        JacobiMatrix(i,3) = calc_pC(x1(i),a,b,c,d);
        JacobiMatrix(i,4) = calc_pD(x1(i),a,b,c,d);
    end
    delta_w = inv((JacobiMatrix)'*(JacobiMatrix))*JacobiMatrix'*(y1-fy);
    %update a b c d
    a = a +lamda*delta_w(1);
    b = b+lamda*delta_w(2);
    c = c +lamda*delta_w(3);
    d= d +lamda*delta_w(4);
end

function [res,R2,fit] = evaluateFit(y,x,w)
fit = getFittingValue(x,w);
res =  norm(y-fit)/sqrt(length(fit));
yu =  mean(y);
R2 = 1 - norm(y-fit)^2/norm(y - yu)^2;
end
function fit = getFittingValue(x,w)
len = length(x);
fit = ones(len,1);
for i = 1:1:len
    fit(i)  = hypothesis(x(i),w);
end
end
function val  = hypothesis(x,w)
a = w(1);b= w(2);c= w(3);d= w(4);
val = d+(a-d)/(1+(x/c)^b);
end
function val = calc_pA(x,A,B,C,D)
val = 1/((x/C)^B + 1);
end
function val = calc_pB(x,A,B,C,D)
val = -(log(x/C)*(A - D)*(x/C)^B)/((x/C)^B + 1)^2;
end
function val = calc_pC(x,A,B,C,D)
val = (B*x*(A - D)*(x/C)^(B - 1))/(C^2*((x/C)^B + 1)^2);
end
function val = calc_pD(x,A,B,C,D)
val = 1 - 1/((x/C)^B + 1);
end