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

在上一篇博客中，自己介绍了Levenberg_Marquardt的算法流程，特点以及在非线性最小二乘问题上的应用，信赖域算法也是自己曾经研究过的算法，并且在姿态估计上进行了应用，比较下来，得到的精度和Levenberg_Marquardt算法不相上下，但是收敛速度却不如LM。本文介绍一下自己对信赖域算法的理解，与童鞋们分享一下。

相信做过机器学习的童鞋，一定使用过搜索算法，例如梯度下降法，牛顿法，拟牛顿法，BFGS算法，LM算法等，这些算法有一个共同的特点，它们的基本策略，也就是所谓的算法思想，都是“定方向，定步长”，不知道我这么概括准确不准确，比如说经典到骨子里的梯度下降，先选择梯度相反方向作为搜索方向，然后要么一维搜索，要么索性定个步长公式，沿着既定方向根据步长确定变量更新值，无一例外。然而，半路杀出了个程咬金，信赖域算法却不是这样干的，人家采用的方法是“定区域，定极值”的方法，来达到同样的目的。牛逼吧？我替你回答， 是的。这个毁三观的行为，造就了有趣的信赖域算法。

信赖域算法是在一个球形区域内进行搜索，具体说来，先定好一个球形区域，然后在这个球形区域的中心点进行二阶泰勒近似（和LM是不是很像？LM是一阶泰勒近似），如公式（1）所示，这个近似美其名曰“子问题”，由于二阶泰勒近似其实就是一个二次规划问题，可以直接求取极值，得到了子问题的最优解之后，下一步就要考察这个最优解（也就是参数变化量d）是否能够使得“原问题”不断的下降，这里，引入了一个比较巧妙的度量方式，用于度量函数值实际下降量和预测下降量之间的比值，如公式（2）所示，如果这个比值大于一定的数值，说明参数变化量d是正确的，进一步，调整球形的半径r，否则，就要摒弃这个d，减小半径r，重新搜索。精髓一句话，每一步迭代都会转化为子问题来处理，不断地更新参数，最终得到收敛解。

（1）

（2）

正如上一篇博客一样，担心大家听不懂我在说什么，我还是画个算法流程图，如下所示：

同样也附上算法代码，这里我还是用信赖域算法，去解决任何一个非线性最小二乘问题，输入变量的含义均已进行注释。

下一篇博客，我会比较LM与信赖域算法之间对于同一问题的处理速度和算法精度，我觉得这两个算法有必要好好的比较一下。

如果有bug或error，还请多多指教啊！小弟嗷嗷感激！！

[plain] view plain copy

1. %% 信赖域算法  
2. function[x, iter] = Trust_Region(f, var, x0, r0, miu, yita, eps)  
3. % 目标函数:                       f  
4. % 自变量向量：                    var  
5. % 初始点：                        x0  
6. % 初始信赖域半径：                r0  
7. % 初始参数：                      miu  
8. % 初始参数：                      yita  
9. % 精度：                          eps  
10. % 目标函数取最小值的自变量值：     x  
11. % 迭代次数：                      iter  
12. tol = 1;  
13. r   = r0;  
14. x   = x0;  
15.   
16. % 迭代次数  
17. iter = 1;  
18.   
19. while tol > eps     
20.     jacf    = jacobian(f, var);             % 雅可比矩阵  
21.     hessen  = jacobian(jacf, var);          % 赫森矩阵      
22.           
23.     % 计算目标函数值  
24.     fx      = double(subs(f, var, x));  
25.       
26.     % 计算雅克比矩阵值  
27.     v       = subs(jacf, var, x);  
28.       
29.     % 计算赫森矩阵值  
30.     pv      = subs(hessen, var, x);  
31.       
32.     tol     = double(norm(v));  
33.     M_2     = double(pv);                   % 二次规划 中的二次项矩阵      
34.     M_1     = double(transpose(v));         % 二次规划 中的一次项矩阵  
35.     lb      = -r * ones(length(var), 1);    % r为信赖域半径，自变量下界约束  
36.     ub      = r * ones(length(var), 1);     % 自变量上界约束  
37.       
38.     % 求解二次规划  
39.     [y, fy, EXITFLAG] = quadprog(M_2, M_1, [], [], [], [], lb, ub);  
40.           
41.     % 重新计算目标函数值  
42.     fx_n = double(subs(f, var, x + y));  
43.       
44.     % 计算目标函数实际下降与预测下降之比  
45.     p = (fx - fx_n) / (-fy);  
46.           
47.     if p <= miu % 目标函数实际下降不明显  
48.         r = 0.5 * r;      
49.     else % 目标函数值下降明显     
50.         x = x + y;    % 更新参数, 扩大信赖域半径  
51.           
52.         if p >= yita    % 如果 p > yita  
53.             r = 2 * r;  
54.         end  
55.     end  
56.           
57.     iter = iter + 1;  
58. end  
59. end  

感慨一下，matlab具有强大的符号计算能力，实在是太方便了，总有人说matlab不能算是编程语言，我是超级不同意的，matlab是一门强大的编程语言好不好。尤其是在算法方面。