计算数学
变尺度之DFP方法
在牛顿法求函数极值时，需要计算二阶偏导数矩阵，而目标函数的偏导数矩阵可能是非正定的，这样的话容易造成计算失败。 为了克服以上缺点，拟牛顿法被提出用于处理多元函数的极值问题。 拟牛顿法的基本思想是用不包含二阶导数的矩阵近似牛顿法中Hesse矩阵的逆矩阵。由于构造逆矩阵的方法不同，出现了种类繁多的拟 牛顿算法系列。 其中影响最大的当属本节介绍的DFP方法与下一节介绍的BFGS方法。DFP方法矩阵校正格式为 \[{{\bf{H}}_{k + 1}}{\bf{ = }}{{\bf{H}}_k}{\bf{ + }}\frac{{{{\bf{p}}_k}{\bf{p}}_k^T}}{{{\bf{p}}_k^T{{\bf{q}}_k}}}{\bf{ - }}\frac{{{{\bf{H}}_k}{{\bf{q}}_k}{\bf{q}}_k^T{{\bf{H}}_k}}}{{{\bf{q}}_k^T{{\bf{H}}_k}{{\bf{q}}_k}}}\] 这里略去原理推导，直接给出算法。[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, ...做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{ = }}{{\bf{H}}_k}{\bf{ + }}\frac{{{{\bf{p}}_k}{\bf{p}}_k^T}}{{{\bf{p}}_k^T{{\bf{q}}_k}}}{\bf{ - }}\frac{{{{\bf{H}}_k}{{\bf{q}}_k}{\bf{q}}_k^T{{\bf{H}}_k}}}{{{\bf{q}}_k^T{{\bf{H}}_k}{{\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 cDFP
/*--------------------------------------class comment
Version : V1.0
Coded by : syz
Date : 2011-07-11 15:24:20 *星期一*
----------------------------------------------------
Desciption : 变尺度之DFP方法
*
parameters :
* fout-------输出文件对象
*
Methods :
* DFP-------DFP方法函数
--------------------------------------------------*/
{
public StreamWriter fout=new StreamWriter("fout.txt");
public void DFP(out VEC x,out double fx,VEC x0)
/*------------------------------------------ comment
Author : syz
Date : 2011-07-11 15:27:17
---------------------------------------------------
Desciption : DFP方法
*
Post Script :
* 这里的矩阵与向量运算较为复杂，我们把复杂的矩阵向量运算
* 通过算符重载的方法写入矩阵和向量类
* 如此，这里便可以像MATLAB那样较为方便的使用各种矩阵运算方法
* 注意在DFP中，没有矩阵求逆或解线性方程这一步，正是变尺度法的优势
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);
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;
//这里运算符多次重载
//意义也有很多重
//这里之所以可以如此简洁
//是因为通过算符重载把复杂的矩阵运算
//放在了矩阵与向量类中
H = H + (p ^ p) / (p | q)
- (((H | q) ^ q) | H )/ (q |( H | q));
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;
}
}