【转】 Levenberg-Marquardt 算法快速入门教程

原帖地址：http://blog.csdn.net/a383201241/article/details/46299861

本文附的源程序是MATLAB代码，总共不到80行，实现了 求雅克比矩阵的解析解，演示了Levenberg-Marquardt最优化迭代过程，演示了如何求解拟合问题。本文用图文介绍了LM算法。转帖来自蜜蜂电脑

什么是最优化，可分为几大类？

答：Levenberg-Marquardt算法是最优化算法中的一种。最优化是寻找使得函数值最小的参数向量。它的应用领域非常广泛，如：经济学、管理优化、网络分析 、最优设计、机械或电子设计等等。

根据求导数的方法，可分为2大类。第一类，若f具有解析函数形式，知道x后求导数速度快。第二类，使用数值差分来求导数。根据使用模型不同，分为非约束最优化、约束最优化、最小二乘最优化。

什么是Levenberg-Marquardt算法？

它是使用最广泛的非线性最小二乘算法，中文为列文伯格-马夸尔特法。它是利用梯度求最大（小）值的算法，形象的说，属于“爬山”法的一种。它同时具有梯度法和牛顿法的优点。当λ很小时，步长等于牛顿法步长，当λ很大时，步长约等于梯度下降法的步长。在作者的科研项目中曾经使用过多次。图1显示了算法从起点，根据函数梯度信息，不断爬升直到最高点（最大值）的迭代过程。共进行了12步。（备注：图1中绿色线条为迭代过程，但是由于分辨率小，看得不太清楚，单击该图后可放大查看）。



图1 LM算法迭代过程形象描述

图1中，算法从山脚开始不断迭代。可以看到，它的寻优速度是比较快的，在山腰部分直接利用梯度大幅度提升（参见后文例子程序中lamda较小时），快到山顶时经过几次尝试（lamda较大时），最后达到顶峰（最大值点），算法终止。

如何快速学习LM算法？

学习该算法的主要困难是入门难。 要么国内中文教材太艰涩难懂，要么太抽象例子太少。目前，我看到的最好的英文入门教程是K. Madsen等人的《Methods for non-linear least squares problems》本来想把原文翻译一下，贴到这里。请让我偷个懒吧。能找到这里的读者，应该都是E文好手，我翻译得不清不楚，反而事倍功半了。

该文链接：http://www2.imm.dtu.dk/pubdb/public/publications.php? year=&pubtype=7&pubsubtype=§ion=1&cmd=full_view&lastndays=&order=author。或者直接下载pdf原文：http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/3215/pdf/imm3215.pdf

例子程序（MATLAB源程序）

本程序不到100行，实现了 求雅克比矩阵的解析解，Levenberg-Marquardt最优化迭代，演示了如何求解拟合问题。采用萧树铁主编的《数学试验》（第二版）（高等教育出版社）中p190例2（血药浓度）来演示。在MATLAB中可直接运行得到最优解。

<code class="language-matlab hljs  has-numbering" style="display: block; padding: 0px; background-color: transparent; color: inherit; box-sizing: border-box; font-family: 'Source Code Pro', monospace;font-size:undefined; white-space: pre; border-top-left-radius: 0px; border-top-right-radius: 0px; border-bottom-right-radius: 0px; border-bottom-left-radius: 0px; word-wrap: normal; background-position: initial initial; background-repeat: initial initial;">    <span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">% 计算函数f的雅克比矩阵，是解析式</span>
syms a b y x <span class="hljs-built_in" style="color: rgb(102, 0, 102); box-sizing: border-box;">real</span>;
f=a*<span class="hljs-built_in" style="color: rgb(102, 0, 102); box-sizing: border-box;">exp</span>(-b*x);
Jsym=jacobian(f,<span class="hljs-matrix" style="box-sizing: border-box;">[a b]</span>)

<span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">% 拟合用数据。参见《数学试验》，p190，例2</span>
data_1=<span class="hljs-matrix" style="box-sizing: border-box;">[<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">0.25</span> <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">0.5</span> <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1</span> <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1.5</span> <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">2</span> <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">3</span> <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">4</span> <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">6</span> <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">8</span>]</span>;
obs_1=<span class="hljs-matrix" style="box-sizing: border-box;">[<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">19.21</span> <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">18.15</span> <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">15.36</span> <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">14.10</span> <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">12.89</span> <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">9.32</span> <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">7.45</span> <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">5.24</span> <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">3.01</span>]</span>;

<span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">% 2. LM算法</span>
<span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">% 初始猜测s</span>
a0=<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">10</span>; b0=<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">0.5</span>;
y_init = a0*<span class="hljs-built_in" style="color: rgb(102, 0, 102); box-sizing: border-box;">exp</span>(-b0*data_1);
<span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">% 数据个数</span>
Ndata=<span class="hljs-built_in" style="color: rgb(102, 0, 102); box-sizing: border-box;">length</span>(obs_1);
<span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">% 参数维数</span>
Nparams=<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">2</span>;
<span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">% 迭代最大次数</span>
n_iters=<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">50</span>;
<span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">% LM算法的阻尼系数初值</span>
lamda=<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">0.01</span>;

<span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">% step1: 变量赋值</span>
updateJ=<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1</span>;
a_est=a0;
b_est=b0;

<span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">% step2: 迭代</span>
<span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">for</span> it=<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1</span>:n_iters
<span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">if</span> updateJ==<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1</span>
<span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">% 根据当前估计值，计算雅克比矩阵</span>
J=<span class="hljs-built_in" style="color: rgb(102, 0, 102); box-sizing: border-box;">zeros</span>(Ndata,Nparams);
<span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">for</span> <span class="hljs-built_in" style="color: rgb(102, 0, 102); box-sizing: border-box;">i</span>=<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1</span>:<span class="hljs-built_in" style="color: rgb(102, 0, 102); box-sizing: border-box;">length</span>(data_1)
J(<span class="hljs-built_in" style="color: rgb(102, 0, 102); box-sizing: border-box;">i</span>,:)=<span class="hljs-matrix" style="box-sizing: border-box;">[exp(-b_est*data_1(i)) -a_est*data_1(i)*exp(-b_est*data_1(i))]</span>;
<span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">end</span>
<span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">% 根据当前参数，得到函数值</span>
y_est = a_est*<span class="hljs-built_in" style="color: rgb(102, 0, 102); box-sizing: border-box;">exp</span>(-b_est*data_1);
<span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">% 计算误差</span>
d=obs_1-y_est;
<span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">% 计算（拟）海塞矩阵</span>
H=<span class="hljs-transposed_variable" style="box-sizing: border-box;">J'</span>*J;
<span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">% 若是第一次迭代，计算误差</span>
<span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">if</span> it==<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1</span>
e=<span class="hljs-built_in" style="color: rgb(102, 0, 102); box-sizing: border-box;">dot</span>(d,d);
<span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">end</span>
<span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">end</span>

<span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">% 根据阻尼系数lamda混合得到H矩阵</span>
H_lm=H+(lamda*<span class="hljs-built_in" style="color: rgb(102, 0, 102); box-sizing: border-box;">eye</span>(Nparams,Nparams));
<span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">% 计算步长dp，并根据步长计算新的可能的\参数估计值</span>
dp=inv(H_lm)*(<span class="hljs-transposed_variable" style="box-sizing: border-box;">J'</span>*d(:));
g = <span class="hljs-transposed_variable" style="box-sizing: border-box;">J'</span>*d(:);
a_lm=a_est+dp(<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1</span>);
b_lm=b_est+dp(<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">2</span>);
<span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">% 计算新的可能估计值对应的y和计算残差e</span>
y_est_lm = a_lm*<span class="hljs-built_in" style="color: rgb(102, 0, 102); box-sizing: border-box;">exp</span>(-b_lm*data_1);
d_lm=obs_1-y_est_lm;
e_lm=<span class="hljs-built_in" style="color: rgb(102, 0, 102); box-sizing: border-box;">dot</span>(d_lm,d_lm);
<span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">% 根据误差，决定如何更新参数和阻尼系数</span>
<span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">if</span> e_lm<e
lamda=lamda/<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">10</span>;
a_est=a_lm;
b_est=b_lm;
e=e_lm;
<span class="hljs-built_in" style="color: rgb(102, 0, 102); box-sizing: border-box;">disp</span>(e);
updateJ=<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1</span>;
<span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">else</span>
updateJ=<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">0</span>;
lamda=lamda*<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">10</span>;
<span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">end</span>
<span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">end</span>
<span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">%显示优化的结果</span>
a_est
b_est</code><ul class="pre-numbering" style="box-sizing: border-box; position: absolute; width: 50px; background-color: rgb(238, 238, 238); top: 0px; left: 0px; margin: 0px; padding: 6px 0px 40px; border-right-width: 1px; border-right-style: solid; border-right-color: rgb(221, 221, 221); list-style: none; text-align: right;"><li style="box-sizing: border-box; padding: 0px 5px;">1</li><li style="box-sizing: border-box; padding: 0px 5px;">2</li><li style="box-sizing: border-box; padding: 0px 5px;">3</li><li style="box-sizing: border-box; padding: 0px 5px;">4</li><li style="box-sizing: border-box; padding: 0px 5px;">5</li><li style="box-sizing: border-box; padding: 0px 5px;">6</li><li style="box-sizing: border-box; padding: 0px 5px;">7</li><li style="box-sizing: border-box; padding: 0px 5px;">8</li><li style="box-sizing: border-box; padding: 0px 5px;">9</li><li style="box-sizing: border-box; padding: 0px 5px;">10</li><li style="box-sizing: border-box; padding: 0px 5px;">11</li><li style="box-sizing: border-box; padding: 0px 5px;">12</li><li style="box-sizing: border-box; padding: 0px 5px;">13</li><li style="box-sizing: border-box; padding: 0px 5px;">14</li><li style="box-sizing: border-box; padding: 0px 5px;">15</li><li style="box-sizing: border-box; padding: 0px 5px;">16</li><li style="box-sizing: border-box; padding: 0px 5px;">17</li><li style="box-sizing: border-box; padding: 0px 5px;">18</li><li style="box-sizing: border-box; padding: 0px 5px;">19</li><li style="box-sizing: border-box; padding: 0px 5px;">20</li><li style="box-sizing: border-box; padding: 0px 5px;">21</li><li style="box-sizing: border-box; padding: 0px 5px;">22</li><li style="box-sizing: border-box; padding: 0px 5px;">23</li><li style="box-sizing: border-box; padding: 0px 5px;">24</li><li style="box-sizing: border-box; padding: 0px 5px;">25</li><li style="box-sizing: border-box; padding: 0px 5px;">26</li><li style="box-sizing: border-box; padding: 0px 5px;">27</li><li style="box-sizing: border-box; padding: 0px 5px;">28</li><li style="box-sizing: border-box; padding: 0px 5px;">29</li><li style="box-sizing: border-box; padding: 0px 5px;">30</li><li style="box-sizing: border-box; padding: 0px 5px;">31</li><li style="box-sizing: border-box; padding: 0px 5px;">32</li><li style="box-sizing: border-box; padding: 0px 5px;">33</li><li style="box-sizing: border-box; padding: 0px 5px;">34</li><li style="box-sizing: border-box; padding: 0px 5px;">35</li><li style="box-sizing: border-box; padding: 0px 5px;">36</li><li style="box-sizing: border-box; padding: 0px 5px;">37</li><li style="box-sizing: border-box; padding: 0px 5px;">38</li><li style="box-sizing: border-box; padding: 0px 5px;">39</li><li style="box-sizing: border-box; padding: 0px 5px;">40</li><li style="box-sizing: border-box; padding: 0px 5px;">41</li><li style="box-sizing: border-box; padding: 0px 5px;">42</li><li style="box-sizing: border-box; padding: 0px 5px;">43</li><li style="box-sizing: border-box; padding: 0px 5px;">44</li><li style="box-sizing: border-box; padding: 0px 5px;">45</li><li style="box-sizing: border-box; padding: 0px 5px;">46</li><li style="box-sizing: border-box; padding: 0px 5px;">47</li><li style="box-sizing: border-box; padding: 0px 5px;">48</li><li style="box-sizing: border-box; padding: 0px 5px;">49</li><li style="box-sizing: border-box; padding: 0px 5px;">50</li><li style="box-sizing: border-box; padding: 0px 5px;">51</li><li style="box-sizing: border-box; padding: 0px 5px;">52</li><li style="box-sizing: border-box; padding: 0px 5px;">53</li><li style="box-sizing: border-box; padding: 0px 5px;">54</li><li style="box-sizing: border-box; padding: 0px 5px;">55</li><li style="box-sizing: border-box; padding: 0px 5px;">56</li><li style="box-sizing: border-box; padding: 0px 5px;">57</li><li style="box-sizing: border-box; padding: 0px 5px;">58</li><li style="box-sizing: border-box; padding: 0px 5px;">59</li><li style="box-sizing: border-box; padding: 0px 5px;">60</li><li style="box-sizing: border-box; padding: 0px 5px;">61</li><li style="box-sizing: border-box; padding: 0px 5px;">62</li><li style="box-sizing: border-box; padding: 0px 5px;">63</li><li style="box-sizing: border-box; padding: 0px 5px;">64</li><li style="box-sizing: border-box; padding: 0px 5px;">65</li><li style="box-sizing: border-box; padding: 0px 5px;">66</li><li style="box-sizing: border-box; padding: 0px 5px;">67</li><li style="box-sizing: border-box; padding: 0px 5px;">68</li><li style="box-sizing: border-box; padding: 0px 5px;">69</li><li style="box-sizing: border-box; padding: 0px 5px;">70</li><li style="box-sizing: border-box; padding: 0px 5px;">71</li><li style="box-sizing: border-box; padding: 0px 5px;">72</li><li style="box-sizing: border-box; padding: 0px 5px;">73</li><li style="box-sizing: border-box; padding: 0px 5px;">74</li><li style="box-sizing: border-box; padding: 0px 5px;">75</li></ul>

本程序对应的C++实现，已经公开,点此进入。

演示程序求解的问题是《数学试验》（萧树铁,第二版,高等教育出版社）中p190例2。为了方便读者，提供该书籍的数据和目标函数照片（2012年4月1日）。