计算数学
拟牛顿BFGS方法
上一节已经提到DFP与BFGS都是著名的拟牛顿算法系列，BFGS矩阵校正格式为 \[{{\bf{H}}_{k + 1}}{\bf{ = }}\left( {{\bf{I - }}\frac{{{{\bf{p}}_k}{\bf{q}}_k^T}}{{{\bf{p}}_k^T{{\bf{q}}_k}}}} \right){{\bf{H}}_k}\left( {{\bf{I - }}\frac{{{{\bf{p}}_k}{\bf{q}}_k^T}}{{{\bf{p}}_k^T{{\bf{q}}_k}}}} \right){\bf{ + }}\frac{{{{\bf{p}}_k}{\bf{p}}_k^T}}{{{\bf{p}}_k^T{{\bf{q}}_k}}}\] BFGS算法与DFP类似，这里给出其算法。
- [STEP1] 给定迭代初值\({{\bf{x}}^{\left( 1 \right)}} \in {\Re^n}\)，误差容限TOL，最大允许迭代次数IMAX。
- [STEP2] 设置\({{\bf{H}}_1} = {{\bf{I}}_n}\)为单位矩阵，并计算出初值处的梯度 \[{{\bf{g}}_1} = \nabla f\left( {{{\bf{x}}^{\left( 1 \right)}}} \right)\] 若\(k < IMAX\)，则对\(k = 1,2, \cdots \)做STEP3-STEP8。
- [STEP3] $${{\bf{d}}^{\left( k \right)}}{\bf{ = - }}{{\bf{H}}_k}{{\bf{g}}_k}$$
- [STEP4] 从\({{\bf{x}}^{\left( k \right)}}\)出发，沿方向\({{\bf{d}}^{\left( k \right)}}\)搜索， 求步长\({\lambda _k}\)，使得 \[f\left( {{{\bf{x}}^{\left( k \right)}} + {\lambda _k}{{\bf{d}}^{\left( k \right)}}} \right) = \mathop {\min }\limits_{\lambda > 0} \left( {{{\bf{x}}^{\left( k \right)}} + \lambda {{\bf{d}}^{\left( k \right)}}} \right)\] 这一步由一维搜索实现。 求出\({{\lambda _k}}\)之后，令 \[{{\bf{x}}^{\left( {k + 1} \right)}} = {{\bf{x}}^{\left( k \right)}} + \lambda {{\bf{d}}^{\left( k \right)}}\]
- [STEP5] 若 \[\left\| {\nabla f\left( {{{\mathbf{x}}^{\left( {k + 1} \right)}}} \right)} \right\| \leqslant TOL\] 则停止退出循环。最终优化结果取\({{\bf{x}}^{\left( {k + 1} \right)}}\)。
- [STEP6] 计算新的梯度 \[{{\bf{g}}_{k + 1}} = \nabla f\left( {{{\bf{x}}^{\left( {k + 1} \right)}}} \right)\]
- [STEP7] 令 \[{{\bf{p}}^{\left( k \right)}} = {{\bf{x}}^{\left( {k + 1} \right)}} - {{\bf{x}}^{\left( k \right)}}\] \[{{\bf{q}}^{\left( k \right)}} = {{\bf{g}}_{k + 1}} - {{\bf{g}}_k}\] 计算 \[{{\bf{H}}_{k + 1}}{\bf{ = }}\left( {{\bf{I - }}\frac{{{{\bf{p}}_k}{\bf{q}}_k^T}}{{{\bf{p}}_k^T{{\bf{q}}_k}}}} \right){{\bf{H}}_k}\left( {{\bf{I - }}\frac{{{{\bf{p}}_k}{\bf{q}}_k^T}}{{{\bf{p}}_k^T{{\bf{q}}_k}}}} \right){\bf{ + }}\frac{{{{\bf{p}}_k}{\bf{p}}_k^T}}{{{\bf{p}}_k^T{{\bf{q}}_k}}}\]
- [STEP8] 更新迭代变量，令 \[{{\bf{x}}^{\left( k \right)}} = {{\bf{x}}^{\left( {k + 1} \right)}}\] \[{{\bf{g}}_k} = {{\bf{g}}_{k + 1}}\]
- [STEP9] 结束。
其c#代码为如下
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.IO;
class cBFGS
/*--------------------------------------class comment
Version : V1.0
Coded by : syz
Date : 2011-07-11 15:24:20 *星期一*
----------------------------------------------------
Desciption : 变尺度之BFGS方法
*
parameters :
* BFGS各种矩阵运算方法通过矩阵向量类实现
* 请参考作者编写的矩阵向量类
*
Methods :
* BFGS-------BFGS方法
* goldsec----一维搜索方法
* func-------目标函数
* grad-------梯度函数
* drive------驱动函数
--------------------------------------------------*/
{
public StreamWriter fout = new StreamWriter("fout.txt");
public void BFGS(out VEC x, out double fx, VEC x0)
/*------------------------------------------ comment
Author : syz
Date : 2011-07-11 15:27:17
---------------------------------------------------
Desciption : 拟牛顿之BFGS方法
*
Post Script :
*
paramtegers :
* x--------极值点自变量
* fx-------函数极小值
* x0-------初值
*~~~~~~~~~~~~~以下变量放在函数内控制
* imax-------最大允许迭代次数
* tol--------误差容限
-------------------------------------------------*/
{
int imax = 200;
//最大允许迭代次数
double tol = 1e-8;
//误差容限
int N = x0.dim;
//自变量维数
MAT H = new MAT(N, N);
for (int i = 0; i < N; i++)
H[i, i] = 1.0;
//设置H矩阵对角元素
VEC g0, g1, p, q, d;
x = new VEC(N);
grad(out g0, x0);
MAT E;
E = H; //E为单位矩阵
MAT TMP1 = new MAT(N, N);
for (int k = 0; k < imax; k++)
{
//这里为 H*g 矩阵和向量的乘积
//采用运算符重载，重载函数见作者编写的
//矩阵向量类
d = -(H | g0);
double namda;
goldsec(out namda, 0, 1, x0, d);
x = x0 + namda * d;
//调用梯度函数
grad(out g1, x);
p = x - x0;
q = g1 - g0;
//这里运算符多次重载
//意义也有很多重
//这里之所以可以如此简洁
//是因为通过算符重载把复杂的矩阵运算
//放在了矩阵与向量类中
TMP1 = (p ^ q) / (p | q);
TMP1 = E - TMP1;
H = TMP1 | H | (~TMP1)
+ (p ^ q) / (p | q);
//-------以上为BFGS核心部分
fout.WriteLine("{0,5:D3}{1,12:F6}{2,12:F6}",
k, x[0], x[1]);
x0 = x;
g0 = g1;
//算符重载，计算向量的2范数
double err = ~p;
if (err < tol)
break;
}
func(out fx, x);
}
public void goldsec(out double namda, double a, double b, VEC XX, VEC d)
/*------------------------------------------ comment
Author : syz
Date : 2011-07-05 12:11:39
---------------------------------------------------
Desciption : 黄金分割法一维搜索
*
Post Script :
* 控制变量----
* 1.最大迭代次数
* 2.误差容限
* 这两个量放在程序内控制
paramtegers :
* fx----函数极小值
* x-----极小值时候的自变量
* [a,b]--搜索区间
-------------------------------------------------*/
{
const int imax = 200;
//最大允许迭代次数
const double tol = 1e-7;
//误差容限
//这黄金分割比在外面算好，不带入循环
double g = (Math.Sqrt(5) - 1) / 2;
double x1 = a + (1 - g) * (b - a);
double x2 = a + g * (b - a);
double f1, f2;
func(out f1, XX + x1 * d);
func(out f2, XX + x2 * d);
for (int k = 0; k < imax; k++)
{
if (f1 < f2)
{
b = x2;
x2 = x1;
x1 = a + (1 - g) * (b - a);
f2 = f1;
//注意这里的形式，真正的变量是x1,XX和d是参数
func(out f1, XX + x1 * d);
}
else
{
a = x1;
x1 = x2;
x2 = a + g * (b - a);
f1 = f2;
//真正的变量是x2,XX和d是参数
func(out f2, XX + x2 * d);
}
double err = Math.Abs(b - a);
if (err < tol) break;
}
namda = (a + b) / 2;
}
}
拟牛顿法是非线性无约束最优化方法中最为有效的一类方法。该系列方法有许多优点， 如该系列算法不需要计算Hesse矩阵，并且具有二次终止性，对于一般情形有超线性收敛速度。