利用Levenberg_Marquardt算法求解无约束的非线性最小二乘问题~

早就想写这篇文章，但是一直没抽出空，主要是画图比较麻烦，嘿嘿~ 现在介绍如何利用经典的Levenberg_Marquardt算法求解无约束的非线性最小二乘问题。Levenberg_Marquardt算法是以两位数学家命名的搜索算法，它比较于常见的最速下降（又被称作梯度下降），牛顿法等，具有较好的全局收敛性，所以得到了较多的重视与应用。

(转载请注明：http://blog.csdn.net/wsj998689aa/article/details/40826775, 作者：迷雾forest)

总所周知，在合理步长的前提下，最速下降法的搜索方向虽然始终为下降方向，但是收敛速度太慢（这个算法竟然叫做最速下降。。。）牛顿法虽然收敛速度得到了很大提升，但是搜索方向未必是函数下降方向。而Levenberg_Marquardt算法则吸取了二者的优点，通过引入一个简单的参数a，使其既能够快速的下降，又能够保证总体沿着下降方向进行搜索。并且，当陷入局部最优解时，可以调整搜索方向，跳出局部最优解。因此，该算法具有很大的概率收敛到全局最优解，即具有梦寐以求的全局收敛性。

再说说最小二乘问题，其数学形式如下：

(1)

其中，N就是观察数据个数，x当然就是所需求解的对象，如果f(x)为线性函数，如公式（2）所示，此时上述问题就是无约束的线性最小二乘问题，如果f(x)为非线性函数，例如指数函数等等，那么上述问题就被称作无约束的非线性最小二乘问题。

(2)

一般求解非线性最小二乘问题的思路都是一致的，那就是，通过泰勒展开对F(x)进行一阶泰勒近似(公式(3))，然后就顺利成章的把非线性最小二乘问题转换为线性最小二乘问题进行求解。注意，这里采用的是一阶泰勒近似，不是其它阶数，因为只有一阶的情况下，才可以得到线性函数形式，这点从公式（3）中可以明显看出。

(3)

之前做非线性最小二乘问题这块，百思不得其解，为什么线性最小二乘问题就可以直接通过求解析解的方式得到最优解，而非线性最小二乘问题就必须要迭代求取最优解，后来直接推导公式（1）才恍然大悟，原来非线性最小二乘问题，是无法构造出来一个矩阵线性方程组的，当然就无法直接求解，关于这一点，我觉得也值得好好的解释，以后有空写写。

由于关于Levenberg_Marquardt算法的算法流程文字描述比较拗口，所以我直接画了个流程图，会更直观一些，如下所示：

  

其中，Ak和fk分别代表雅克比矩阵和函数值向量，由于十分常见，这里我就不给出数学形式了。强烈提醒大家注意一点，在右侧分支那里，k并不增加，Ak和fk均为之前一步的，改变的只是参数a，我之前这里出现问题，导致我的LM算法无论如何都不收敛，让我郁闷非常，后来仔细回头看算法流程才发现自己的错误，所以我这次画了个图给大家，不要走我的弯路啊！

进一步解释一下LM算法为何牛逼，上面已经说过了，该算法吸取了最速下降法和牛顿法的优点，原因就在于那个参数a，当a趋于0时，算法实际上变成了牛顿法，当a很大时，算法实际上变成了最速下降法，剩下的，就介于二者之间，当发现陷入局部最优时，便会增大a值，使得参数变化加剧，达到跳出局部最优解的目的。

这里，我给出函数值变化的曲线图，通过这个曲线图我们明显可以发现，LM算法可以有效的避免陷入局部最优，并且，总体而言，保持下降方向。

最后，附上自己写的，Levenberg_Marquardt算法Matlab代码，有任何问题，请不吝指教，谢谢！！！不好意思，之前写的那个版本让我替换下去了，因为不够通用。

%% 函数功能:利用Levenberg_Marquardt算法求解无约束的非线性最小二乘问题
function [x_curr, iter] = Levenberg_Marquardt(Diff, f, var, x0, lamda, beita, eps)
    % 相关参数
    max_iter    = 100;
    flag        = 1;
    x_curr      = x0;

    for iter = 1 : max_iter
        if flag
            J   = double(subs(jacobian(Diff, var), var, x_curr));      % 雅可比矩阵
            H   = transpose(J) * J;                                    % 赫森矩阵    
            
            if iter == 1
                f_curr = double(subs(f, var, x_curr));
            end
        end

        % 计算新Hessen矩阵
        H_new = H + (lamda * eye(2, 2));

        % 计算步长d
        diff = double(subs(Diff, var, x_curr));
        d    = inv(H_new) * (J' * diff);
        
        % 更新参数
        x_new = x_curr - d;

        % 计算新函数值
        f_new = double(subs(f, var, x_new));
  
        if f_new < f_curr
            % 收敛条件（左侧分支）
            x_curr = x_new;
            if (norm(d) < eps)
                disp('收敛条件1')
                break;
            end
            
            f_curr = f_new;
            flag = 1;
            lamda = lamda / beita;
        else           
            % 收敛条件（右侧分支）
            if (norm(d) < eps)
                disp('收敛条件2')
                break;
            end
            
            flag = 0;
            lamda = lamda * beita;
        end
    end
end

补充：时隔一年，这一年中有不少朋友问过我关于这篇博客的一个疑问，大家总是不能正确使用Levenberg_Marquardt算法进行非线性最小二乘问题的求解，或者说得到的解总是不尽人意，其实很多人都犯了一个疏漏，就是公式（2），必须强调的是，（2）中的b指的不是函数的参数，而是观测数据，举个例子，如果线性函数是g(x) = a*x + b，那么f(x) = a*x + b - y，一些同学将f(x)错认为就是g(x)，这是错误的，一定要再减去观测数据项，希望大家看到这点补充之后不要再疏忽这一点。

一年了，时间过得真快~~~