计算数学

数值优化

一维搜索之连续抛物线插值

在黄金分割法中，除了比较函数值大小之外，函数值本身并未在最优求解中用到。 这说明，我们并未充分利用现有的信息来加速迭代收敛。这一节介绍连续抛物线插值法 求解一维搜索问题，这是一种较少浪费函数值的方法。
下面阐述连续抛物线插值的基本思想，先在最小值附件选3个点$r,s,t$。计算这3个点的目标函数值$f$，且可以 给出过这3个点的抛物线 \[P\left( x \right) = f\left( r \right) + {d_1}\left( {x - r} \right) + {d_3}\left( {x - r} \right)\left( {x - s} \right)\] 其中 \[{d_1} = \frac{{f\left( s \right) - f\left( r \right)}}{{s - r}}\] \[{d_2} = \frac{{f\left( t \right) - f\left( s \right)}}{{t - s}}\] \[{d_3} = \frac{{{d_2} - {d_1}}}{{t - r}}\] 令\(P\left( x \right)\)可以给出抛物线的极小值公式 \[x = \frac{{r + s}}{2} - \frac{{\left( {f\left( s \right) - f\left( r \right)} \right)\left( {t - r} \right)\left( {t - s} \right)}}{{2\left[ {\left( {s - t} \right)\left( {f\left( t \right) - f\left( s \right)} \right) - \left( {f\left( s \right) - f\left( r \right)} \right)\left( {t - s} \right)} \right]}}\] 此点即为极小值新的近似。有些教材给的算法到此为止，这样的话只能给出较粗的极小值。 在实际应用时，可以重复以上过程，用新的x代替r，s，t中最远的或次最远的点， 反复迭代，即可以给出需求精度下的极小值。
下面给出连续抛物线求一维搜索算法：

• [Input] 最大允许迭代次数、误差容限，初始猜想极小点r，s，t 。
• [Output] 极值点的自变量与函数值。
• [STEP1] 对\(i = 1,2,3, \cdots \)，作如下循环 \[x = \frac{{r + s}}{2} - \frac{{\left( {f\left( s \right) - f\left( r \right)} \right)\left( {t - r} \right)\left( {t - s} \right)}}{{2\left[ {\left( {s - t} \right)\left( {f\left( t \right) - f\left( s \right)} \right) - \left( {f\left( s \right) - f\left( r \right)} \right)\left( {t - s} \right)} \right]}}\] \[t = s\] \[s = r\] \[r = x\] 如果$i$超过与设定的最大允许迭代次数，停止循环，并给出告警信息。 如果已经满足预设的精度需求，则停止循环。
• [STEP2] 输出计算结果，并停止该算法。

一般而言，连续抛物线插值法要快于黄金分割方法，因为它更有效的利用了函数值的信息。 用连续抛物线插值法求解函数 \[f\left( x \right) = {x^6} - 11{x^3} + 17{x^2} - 7x + 1\] 在区间\(\left[ {0,1} \right]\)上的极小值。 go语言版本代码如下。

package main
import (
	"fmt"
	."math"
)
func main() {
/*------------------------------------------------------
  Author  : Song Yezhi                       
  verison : 2021-10-7 18:25
  go build -gcflags "-N -l"  
     
------------------------------------------------------*/
	r :=0.0
	s :=0.7
	t :=1.0
	_ , _ = paraInterp(r,s,t,myfx)

}

func myfx(x float64) float64 {
/*------------------------------------------------------
  Author  : Song Yezhi                       
  verison : 2021-10-7 18:26
  go build -gcflags "-N -l"  
------------------------------------------------------*/
	fx := Pow(x,6)-11.0*Pow(x,3)+17.0*x*x - 7.0*x +1.0
	return fx
}

func paraInterp(r,s,t float64,funcX func(float64) float64) (x, fx float64) {
/*------------------------------------------------------
  Author  : Song Yezhi                       
  verison : 2021-10-7 18:28
  go build -gcflags "-N -l"  
     连续抛物型法一维搜索
  -----------------------------------------------------  
  Input  Parameters :
     r,s,t  ------初值
     funcX -----  目标函数
  Output Parameters :
     x---极小点
     fx--极小点函数值
------------------------------------------------------*/
	var imax int = 20
	var tol float64 = 1e-8
	
	var fs,ft,fr float64
	var tmp1,tmp2 float64
	
	fs = funcX(s)
	fr = funcX(r)
	ft = funcX(t)

	for i := 0; i < imax; i++ {		
		tmp1 = (fs -fr) * (t-r) * (t-s)
		tmp2 = 2.0*((s-r)*(ft-fs) -(fs-fr)*(t-s))
		
		x = (r+s)/2.0 - tmp1/tmp2
		
		if Abs(fs-fr) < tol {
		    break
		}
		
	    t = s
		s = r
		r = x
		
		ft = fs
		fs = fr
		fr = funcX(r)
		// *+++single function evaluation +++*   
		fmt.Printf("i= %4d  x= %12.7f  f(x)=%12.7f \n", i, x, ft)
	}
	return x,fx
}    

多维非线性最优化牛顿下山法

对于多元函数极值问题描述如下 \begin{equation} \min f\left( {\bf{x}} \right) \end{equation} 其中 \begin{equation} {\bf{x}} = \left( \begin{array}{l} {x_1}\\ {x_2}\\ \vdots \\ {x_n} \end{array} \right) \end{equation} 多元函数极值问题有多种数值方法求解，本节先介绍一个基础且重要的方法-牛顿法。 牛顿法的基本思想是把函数在局部用泰勒级数的二次展开取近似原函数， 而二次函数的极值问题是容易求得的，故而逐步迭代，最终得到满足精度需求的方法。 对目标函数进行泰勒展开，保留至二次项 \[f\left( {\bf{x}} \right) \approx \phi \left( {\bf{x}} \right) = f\left( {{{\bf{x}}_k}} \right) + {\left[ {\nabla f\left( {{{\bf{x}}_k}} \right)} \right]^T}\left( {{\bf{x}} - {{\bf{x}}_k}} \right) + \frac{1}{2}{\left( {{\bf{x}} - {{\bf{x}}_k}} \right)^T}{\nabla ^2}f\left( {{{\bf{x}}_k}} \right)\left( {{\bf{x}} - {{\bf{x}}_k}} \right)\] 求其驻值 \[\nabla \phi \left( {\bf{x}} \right) = {\bf{0}}\] 即 \[\nabla f\left( {{{\bf{x}}_k}} \right) + {\nabla ^2}f\left( {{{\bf{x}}_k}} \right)\left( {{\bf{x}} - {{\bf{x}}_k}} \right) = {\bf{0}}\] 从而有下面的迭代格式 \[\left\{ \begin{array}{l} {\bf{H}}\left( {{{\bf{x}}_k}} \right){\bf{v = - }}\nabla {\bf{f}}\left( {{{\bf{x}}_k}} \right)\\ {{\bf{x}}_{k + 1}} = {{\bf{x}}_k} + {\bf{v}} \end{array} \right.\] 其中$$\nabla {\bf{f}}\left( {{{\bf{x}}_k}} \right)$$为梯度向量 \[\nabla {\bf{f}}\left( {{{\bf{x}}_k}} \right) = \left[ \begin{array}{l} \frac{{\partial f}}{{\partial {x_1}}}\left( {{x_1}, \cdots ,{x_n}} \right)\\ \quad \quad \vdots \\ \frac{{\partial f}}{{\partial {x_n}}}\left( {{x_1}, \cdots ,{x_n}} \right) \end{array} \right]\] $${\bf{H}}\left( {{{\bf{x}}_k}} \right)$$ 为Hessian矩阵 \[{\bf{H}} = \left[ {\begin{array}{*{20}{c}} {\frac{{{\partial ^2}f}}{{\partial {x_1}\partial {x_1}}}}& \cdots &{\frac{{{\partial ^2}f}}{{\partial {x_1}\partial {x_n}}}}\\ \vdots & \ddots & \vdots \\ {\frac{{{\partial ^2}f}}{{\partial {x_n}\partial {x_1}}}}& \cdots &{\frac{{{\partial ^2}f}}{{\partial {x_n}\partial {x_n}}}} \end{array}} \right]\] c#版本代码如下。

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.IO;
class cnewton
/*--------------------------------------class comment
 Version   :  V1.0
 Coded by  :  syz
 Date      :  2011-07-11 08:46:04          *星期一*
----------------------------------------------------
Desciption :  牛顿下山法计算最优化问题
Methods    :
 *    newton--------牛顿下山方法函数
 *    pivot---------选主元消去法解线性方程
 *    uptri---------上三角方程回代方法
 *    func----------目标函数
 *    grad----------梯度函数
 *    hess----------hessian矩阵
--------------------------------------------------*/
{
    public StreamWriter fout = new StreamWriter("fout.txt");
    //输出文件
    public void newton(out VEC x, out double fx,VEC x0)
    /*------------------------------------------ comment
    Author      : syz 
    Date        : 2011-07-11 08:52:58         
    ---------------------------------------------------
    Desciption  :  牛顿下山法计算非线性最优化
     *   注意：牛顿法形式上需要计算Hessian矩阵的逆，
     *   实际计算时直接解算线性方程
    paramtegers :
     * x------极小值点
     * fx-----极值点的函数值
     * x0-----初值向量
     * 
     * ~~~~以下量在程序开头部分控制，也可以把控制量放在参数中，由外部控制
     * imax------最大迭代次数
     * tol-------误差容限
    -------------------------------------------------*/
    {
        int imax = 50;
        //最大允许迭代次数
        double tol = 1e-7;
        //误差容限
        VEC dx, x1, df;
        MAT H;
        //这里皆不指定维数大小
        for (int k = 0; k < imax; k++)
        {
            grad(out df, x0);
            //求的梯度向量
            hess(out H, x0);
            //求Hessian矩阵
            //实际方程为  H dx= -df
            //这里转一次负号，重载算符
            //df = 0 - df;
            pivot(H, df, out dx);
            x1 = x0 - dx;
            //------------------------
            //向量的模，算符重载
            double err = ~dx;
            //满足精度时退出循环
            if (err < tol)
                break;
            fout.WriteLine("{0,5:D3}{1,12:F6}{2,12:F6}",
                k, x1[0], x1[1]);
            x0 = x1;
        }
        x = x0;
        func(out fx, x);
    }
    public void pivot(MAT A, VEC b, out VEC x)
    {
        /*------------------------------------------ comment
        Author      : syz 
        Date        : 2011-06-20 10:59:59         
        ---------------------------------------------------
        Desciption  :
         *      消去法解线性方程组
        Post Script :
         *      要用到上三角方程的回代方法
        paramtegers :
         *     A---方程组系数矩阵
         *     b---右向量
         *     x---计算结果
        -------------------------------------------------*/
        int n, i, j, k;
        n = b.dim;
        //方程维数
        double tmp;
        //增广矩阵
        MAT Ab = new MAT(n, n + 1);
        //赋值给增广矩阵
        for (i = 0; i < n; i++)
        {
            for (j = 0; j < n; j++)
                Ab.ele[i, j] = A.ele[i, j];
            Ab.ele[i, n] = b.ele[i];
        }
        double elmax;
        int id_max;
        double tmp1, tmp2;
        //这段是高斯消去法的核心
        for (k = 0; k <= n - 2; k++)
        {
            //主元素初值
            elmax = Math.Abs(Ab.ele[k, k]);
            //主元素标号
            id_max = k;
            //下面一段为找出主元素所在的行
            //即获得主元素标号的值
            for (i = k + 1; i <= n - 1; i++)
            {
                if (Math.Abs(Ab.ele[i, k]) > elmax)
                {
                    elmax = Math.Abs(Ab.ele[i, k]);
                    //注意这里把该元素的绝对值复制给最大元素
                    id_max = i;
                }
            }
            //n+1个元素，该循环为交换两行元素
            for (i = 0; i < n + 1; i++)
            {
                tmp1 = Ab.ele[k, i];
                tmp2 = Ab.ele[id_max, i];

                Ab.ele[k, i] = tmp2;
                Ab.ele[id_max, i] = tmp1;
            }
            //交换后实施消去
            for (i = k + 1; i <= n - 1; i++)
            {
                tmp = Ab.ele[i, k] / Ab.ele[k, k];
                for (j = 0; j < n + 1; j++)
                    Ab.ele[i, j] = Ab.ele[i, j] - tmp * Ab.ele[k, j];
            }
        }
        //在C#中使用传值方法，不会改变原始变量的值
        for (i = 0; i < n; i++)
        {
            for (j = 0; j < n; j++)
                A.ele[i, j] = Ab.ele[i, j];
            b.ele[i] = Ab.ele[i, n];
        }
        uptri(A, b, out x);
    }

    public void uptri(MAT A, VEC b, out VEC x)
    {
        /*------------------------------------------ comment
        Author      : syz 
        Date        : 2011-06-17 10:04:27         
        ---------------------------------------------------
        Desciption  :
         *         上三角矩阵方程解法
        Post Script :
         *         要用到MAT 与 VEC 类
        paramtegers :
         *      A-----系数矩阵
         *      b-----右向量
         *      x-----方程的解
        -------------------------------------------------*/
        int N;
        N = A.dim1;
        //获取方程维数
        x = new VEC(N);
        x.ele[N - 1] = b.ele[N - 1] / A.ele[N - 1, N - 1];
        for (int i = N - 2; i >= 0; i--)
        {
            x.ele[i] = b.ele[i];
            for (int j = i + 1; j <= N - 1; j++)
                x.ele[i] = x.ele[i] - A.ele[i, j] * x.ele[j];
            x.ele[i] = x.ele[i] / A.ele[i, i];
        }
    }
}

变尺度之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;
    }
}
   

拟牛顿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] 结束。

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