计算数学

content

拟牛顿BFGS方法

上一节已经提到DFP与BFGS都是著名的拟牛顿算法系列，BFGS矩阵校正格式为 $H_{k + 1} = (I - \frac{p_{k} q_{k}^{T}}{p_{k}^{T} q_{k}}) H_{k} (I - \frac{p_{k} q_{k}^{T}}{p_{k}^{T} q_{k}}) + \frac{p_{k} p_{k}^{T}}{p_{k}^{T} q_{k}}$ BFGS算法与DFP类似，这里给出其算法。

[STEP1] 给定迭代初值 $x^{(1)} \in ℜ^{n}$ ，误差容限TOL，最大允许迭代次数IMAX。
[STEP2] 设置 $H_{1} = I_{n}$ 为单位矩阵，并计算出初值处的梯度 $g_{1} = \nabla f (x^{(1)})$ 若 $k < I M A X$ ，则对 $k = 1, 2, \dots$ 做STEP3-STEP8。
[STEP3] $d^{(k)} = - H_{k} g_{k}$
[STEP4] 从 $x^{(k)}$ 出发，沿方向 $d^{(k)}$ 搜索，求步长 $λ_{k}$ ，使得 $f (x^{(k)} + λ_{k} d^{(k)}) = min_{λ > 0} (x^{(k)} + λ d^{(k)})$ 这一步由一维搜索实现。求出 $λ_{k}$ 之后，令 $x^{(k + 1)} = x^{(k)} + λ d^{(k)}$
[STEP5] 若 $‖ \nabla f (x^{(k + 1)}) ‖ ⩽ T O L$ 则停止退出循环。最终优化结果取 $x^{(k + 1)}$ 。
[STEP6] 计算新的梯度 $g_{k + 1} = \nabla f (x^{(k + 1)})$
[STEP7] 令 $p^{(k)} = x^{(k + 1)} - x^{(k)}$ $q^{(k)} = g_{k + 1} - g_{k}$ 计算 $H_{k + 1} = (I - \frac{p_{k} q_{k}^{T}}{p_{k}^{T} q_{k}}) H_{k} (I - \frac{p_{k} q_{k}^{T}}{p_{k}^{T} q_{k}}) + \frac{p_{k} p_{k}^{T}}{p_{k}^{T} q_{k}}$
[STEP8] 更新迭代变量，令 $x^{(k)} = x^{(k + 1)}$ $g_{k} = g_{k + 1}$
[STEP9] 结束。

许多数值实验证实，BFGS比DFP有更好的稳定性。

其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矩阵，并且具有二次终止性，对于一般情形有超线性收敛速度。