优化算法之最速梯度下降法、牛顿法、拟牛顿法

一、最速梯度下降法

我们知道常规的梯度下降法迭代法公式如下：
$${\theta }^{(k+1)}={\theta }^{(k)}-\eta \nabla f({\theta }^{(k)})$$
迭代公式中包含了 「方向」和「步长」两个参数。

1.1 何谓最速？

那么该如何理解最速梯度下降法中的 最速呢？

因为是要最快"下山"，找到越来越小的函数值，那么就希望每次迭代生成的$f({\theta }^{(k)}-\eta \nabla f({\theta }^{(k)}))$ 越小越好，如果找到使这个函数值最小的步长，则意味是最快速下降。

最速就是通过函数值最小，来求出对应的最优步长$\eta$ ，接着代入继续迭代。

最优步长的公式如下：

$${\eta }^{(k)}=argminf(f({\theta }^{(k)}-\eta \nabla f({\theta }^{(k)})))$$

1.2 最速梯度下降法的迭代过程

1.3 最速梯度下降法示例

$$\begin{gathered} 求解多元函数f(\theta )={\theta }_1^2+2{\theta }_2^2-2{\theta }_1{\theta }_2-4{\theta }_1的极小值点。 \\ 按照最速梯度下降法的迭代过程： \\ 第1次迭代过程如下： \\ 第一步，选取初始值k=0,{\theta }^{(0)}=(1,1) \\ 第二步，计算f({\theta }^{(k)})以及梯度\nabla f({\theta }^{(k)})。 \\ f({\theta }^{(k)})=f({\theta }_1=1,{\theta }_2=1)=-3 \\ \nabla f({\theta }^{(k)})=\left(\begin{array}{c}\frac{\partial f(\theta )}{\partial {\theta }_1},\frac{\partial f(\theta )}{\partial {\theta }_2}\end{array}\right)=\left(\begin{array}{c}2{\theta }_1-2{\theta }_2-4,4{\theta }_2-2{\theta }_1\end{array}\right)=(-4,2) \\ 第三步，计算最优步长\eta 。 \\ f(\eta )=f({\theta }^{(k)}-\eta \nabla f({\theta }^{(k)}))=f((1,1)-\eta (-4,2)) \\ =f((1+4\eta ,1-2\eta )) \\ =40{\eta }^2-20\eta -3 \\ 令f^{\prime }(\eta )=80\eta -20=0，得到\eta =\frac{1}{4} \\ 第四步、利用迭代公式进行参数更新。 \\ {\theta }^{(k+1)}={\theta }^{(k)}-{\eta }^{(k)}\nabla f({\theta }^{(k)})) \\ 此时k=0,那么{\theta }^{(1)}={\theta }^{(0)}-{\eta }^{(0)}\nabla f({\theta }^{(0)})) \\ =(1,1)-\frac{1}{4}\ast (-4,2) \\ =(2,\frac{1}{2}) \\ 第五步，计算||f({\theta }^{(k+1)})-f({\theta }^{(k)})||=||f({\theta }^{(1)})-f({\theta }^{(0)})||>精度 \\ 令k=k+1=1，转第三步继续迭代，更新参数，直到满足终止条件。 \end{gathered}$$

1.4 代码实现示例的迭代过程

我们用python中的sympy包来实现最速下降法求解上面的问题。

1.4.1 sympy包的几个函数

首先，我们来介绍我们将要用到的sympy包中的几个函数。结合jupyter notebook来给大家做个展示。

1、构建符号变量和符号函数

这个是用来构造两个符号变量$x_1,x_2$,就像代数中用字母代替变量一样。然后可以定义出我们的函数。

x_1 = symbols('x_1')
x_2 = symbols('x_2') 
              
fun = x_1**2 +  2 * x_2**2 -2 * x_1 * x_2 - 4 * x_1
fun  

可以看到jupyter notebook中直接就显示出了数学公式格式的形式，这是因为jupyter notebook中内嵌了LaTeX相关支持包的缘故。

2、对符号函数求导

下面计算$f(x_1,x_2)$的偏导数，代码很简单，用函数diff(函数, 变量)

grad_1 = diff(fun, x_1) 
grad_1

3、求函数值
有了符号函数，我们怎么知道自变量$x_1,x_2$取具体值的时候，符号函数的取值呢？

我们用函数函数.subs({变量1:变量1的取值, …}).evalf()，具体代码为

fun_value = fun.subs({x_1:1, x_2: 1}).evalf()
fun_value  

如果说只有一个变量已知，那就是下面的情况

fun_value = fun.subs({x_1:1}).evalf()
fun_value

4、求解方程的零点

为了寻找下降速度最快的方向，我们需要求解方程组。这里我们用到函数solve(函数，变量)

# 假设有如下函数
fun2 = 4 * x_1 ** 2 - 4
x_1_solution = solve(fun2, x_1)
x_1_solution

1.4.2 利用sympy包求解1.3案例

# pip install sympy -i http://mirrors.aliyun.com/pypi/simple/ --trusted-host mirrors.aliyun.com
from sympy import symbols, diff, solve
import numpy as np

x1 = symbols("x1")
x2 = symbols("x2")
λ =  symbols("λ")
f = x1 ** 2 + 2 * x2 ** 2 - 2 * x1 * x2 - 4 * x1


def gradient_descent(x1_init, x2_init, ε):
    # 1、求解x1及x1的偏导数
    grad_1 = diff(f, x1)
    grad_2 = diff(f, x2)

    x1_curr = x1_init
    x2_curr = x2_init

    count = -1

    while True:
        count = count + 1
        # 2、求解偏导数的零点值
        grad_1_value = grad_1.subs({x1: x1_curr, x2: x2_curr})
        grad_2_value = grad_2.subs({x1: x1_curr, x2: x2_curr})

        d1 = -grad_1_value
        d2 = -grad_2_value
        k = np.array([d1, d2], dtype=float)
        # 这里是求解梯度大小
        # ord，表示范数的种类，默认为2范数。2范数，表示向量中各个元素平方和 的 1/2 次方
        # axis, axis=0 表示按列向量来进行处理，求多个列向量的范数;
        #       axis =1 表示按行向量来进行处理，求多个行向量的范数
        norm_result = np.linalg.norm(k, ord=2, axis=0)
        # 3、计算f(λ)的表达式，f_new=f(λ)
        x1_new = x1_curr + λ * d1
        x2_new = x2_curr + λ * d2

        f_new = f.subs({x1: x1_new, x2: x2_new})
        # 4、对f(λ)求导，令导数为0，求取λ值
        grad_3 = diff(f_new, λ)
        λ_value = solve(grad_3, λ)[0]

        print("第{}次迭代".format(count), "λ的值为{}".format(λ_value))
        # 4、计算参数更新后的值
        x1_curr = x1_new.subs(λ, λ_value)
        x2_curr = x2_new.subs(λ, λ_value)
        print("当前x的值为[%f, %f]" % (float(x1_curr), float(x2_curr)))
        # 5、判断是否满足迭代结束条件
        if norm_result < ε:
            print("该精度下的最优解是[%f, %f]" % (float(x1_curr), float(x2_curr)))
            break

    return x1_curr, x2_curr


if __name__ == '__main__':
    result = gradient_descent(x1_init=1, x2_init=1, ε=0.01)

第0次迭代 λ的值为1/4
当前x的值为[2.000000, 0.500000]
第1次迭代 λ的值为1/2
当前x的值为[2.500000, 1.500000]
第2次迭代 λ的值为1/4
当前x的值为[3.000000, 1.250000]
第3次迭代 λ的值为1/2
当前x的值为[3.250000, 1.750000]
第4次迭代 λ的值为1/4
当前x的值为[3.500000, 1.625000]
第5次迭代 λ的值为1/2
当前x的值为[3.625000, 1.875000]
第6次迭代 λ的值为1/4
当前x的值为[3.750000, 1.812500]
第7次迭代 λ的值为1/2
当前x的值为[3.812500, 1.937500]
第8次迭代 λ的值为1/4
当前x的值为[3.875000, 1.906250]
第9次迭代 λ的值为1/2
当前x的值为[3.906250, 1.968750]
第10次迭代 λ的值为1/4
当前x的值为[3.937500, 1.953125]
第11次迭代 λ的值为1/2
当前x的值为[3.953125, 1.984375]
第12次迭代 λ的值为1/4
当前x的值为[3.968750, 1.976562]
第13次迭代 λ的值为1/2
当前x的值为[3.976562, 1.992188]
第14次迭代 λ的值为1/4
当前x的值为[3.984375, 1.988281]
第15次迭代 λ的值为1/2
当前x的值为[3.988281, 1.996094]
第16次迭代 λ的值为1/4
当前x的值为[3.992188, 1.994141]
第17次迭代 λ的值为1/2
当前x的值为[3.994141, 1.998047]
该精度下的最优解是[3.994141, 1.998047]

1.5 最速梯度下降的注意点

• 相邻两次迭代的搜索方向互相正交。

  • 最速下降法在两个相邻点之间的搜索方向是正交的
  • 最速下降法向极小点逼近是曲折前进的，这种现象被称之为锯齿现象
  • 除特殊的目标函数和极特殊的初始点外，锯齿现象都会发生

• 在一般情况下，当用最速下降法寻找极小点时，其搜索路径呈直角锯齿状（如下图）

  • 在开头几步，目标函数下降较快；
  • 但在接近极小点时，收敛速度长久不理想了。特别适当目标函数的等值 线为比较扁平的椭圆时，收敛就更慢了。
  • 因此，在实用中常用最速下降法和其他方法联合应用，在前期使用最速下降法，而在接近极小值点时，可改用收敛较快的其他方法。

二、牛顿法

2.1 利用牛顿法求解方程

我们先利用牛顿法求解一个方程。
$$g(x)=3x^5-4x^4+6x^3+4x-4$$
我们把曲线画出来，零点貌似就在区间[0.6, 0.8]之间。

• 我们取一个比较大的值，「令 $x_0=1.4$」, 求得$g(x_0)=18.8$。这个值可是相当的大，距离零点挺遥远的，这就是我们最初所在的初始值了。
• 接着往下找。我们现在沿着$x_0$这个点做这个曲线的切线，也就是下图中红色的切线，这个点处的切线与 轴有了一个新的交点，那么可以将新的交点记作下一个迭代值$x_1$，可以求得$x_1=1.045$。
• 那么$g(x_1)=6$，可以发现$x_1$距离$g(x)=0$近了。
• 继续采用同样的方法寻找下一个点，依旧是沿着曲线做切线，也就是我们所得到蓝色的切线，然后就又找到了一个与x轴的交点$x_2=0.787$ 。
• 我们看到$x_2$所对应的函数值$g(x_2)$更小了，比我们刚才那个$g(x_1)=6$还要小，更接近于 0。
• 不停地重复下去，最终就能够找到我们想要的那个零点。

假设我们设定迭代的阈值为0.0001，下表就展示了整个迭代过程，可以看到最终求到的零点为0.6618。

迭代次数	x	|g(x)|
0	1.4000	18.8323
1	1.0447	5.9879
2	0.7873	1.4482
3	0.6769	0.1549
4	0.6620	0.0023
5	0.6618	0.0000

回顾之前整个过程，每次都是在相应的点处求切线，然后找到切线与$x$轴的交点，所以关键点有两个。

• 第一步找到迭代点处的切线；
• 第二步找到切线与$x$轴的交点。

2.2 牛顿法的迭代原理

第1步、用迭代点和方程导函数确定切线。

• 第$k$次的迭代点为$(x^{(k)},g(x^{(k)}))$，过该点做曲线的切线。

• 切线的斜率$g^{\prime }(x^{(k)})$由该点的导函数确定。

• 切线肯定过迭代点，那么切线方程为：$y-g(x^{(k)})=g^{\prime }(x^{(k)})(x-x^{(k)})$

• 整理一下，得到$y=g(x^{(k)})+g^{\prime }(x^{(k)})(x-x^{(k)})$

第2步、用切线与x轴的交点计算函数值。

• 接着我们寻找切线方程与$x$轴的交点：$0=g(x^{(k)})+g^{\prime }(x^{(k)})(x-x^{(k)})$

• 所得交点为新的迭代点:

$$x^{(k+1)}=x^{(k)}-\frac{g(x^{(k)})}{g^{\prime }(x^{(k)})}$$

• 代入方程中，就可以得到第$k+1$次的迭代点为$(x^{(k+1)},g(x^{(k+1)}))$
• 重复以上的迭代步骤

第3步、设置停止迭代的条件

可以从两个方面来设置停止条件

一是设置迭代次数，比如可以设置迭代次数不超过500次。

二是设置阈值精度，比如可以设置精度只要是在0.0001以内，就可以输出零点。

代码示例：

我们利用迭代原理，用代码来模拟一下之前方程的求解过程。

import math


def fun(x0, func, dfunc):
    eps = 1e-10
    max_iter = 500
    step = 0
    x = [0] * (max_iter + 1)
    x[0] = x0
    # 打印初始值
    print('epoch = 0, x = %f, f(x) = %f' % (x0, func(x0)))
    for n in range(max_iter):
        # 利用迭代原理中的公式求取新的迭代点坐标x
        x[n + 1] = x[n] - (func(x[n]) / dfunc(x[n]))
        # 将迭代点x坐标，代入函数得到y
        y = func(x[n + 1])
        step = n + 1
        print('epoch = %d, x = %f, f(x) = %f' % (n + 1, x[n + 1], y) )
        if abs(y) <= eps:
            break
    # 返回零点值
    return x[step]

def func(x):
    return 3 * math.pow(x, 5) - 4 * math.pow(x, 4) + 6 * math.pow(x, 3) + 4 * x - 4


def dfunc(x):
    return 15 * math.pow(x, 4) - 16 * math.pow(x, 3) + 18 * math.pow(x, 2) + 4


if __name__ == '__main__':
    x0 = 1.4
    fun(x0, func, dfunc)

epoch = 0, x = 1.400000, f(x) = 18.832320
epoch = 1, x = 1.044673, f(x) = 5.987861
epoch = 2, x = 0.787330, f(x) = 1.448242
epoch = 3, x = 0.676887, f(x) = 0.154938
epoch = 4, x = 0.662038, f(x) = 0.002285
epoch = 5, x = 0.661812, f(x) = 0.000001
epoch = 6, x = 0.661812, f(x) = 0.000000

2.3 牛顿法的注意事项

• 牛顿法不是总收敛的
• 牛顿法对初始值非常敏感，不同的初始值是否收敛就存在着不同
• 在牛顿法中不同的初始值，收敛值可能不同

我们可以对于第3点举个例子：
$$\begin{gathered} g(x)=x^3-2x^2-11x+12 \\ {g}^{\prime }(x)=3x^2-4x-11 \end{gathered}$$
还是应用上述代码，进行求解

import math


def fun(x0, func, dfunc):
    eps = 1e-10
    max_iter = 500
    step = 0
    x = [0] * (max_iter + 1)
    x[0] = x0
    # 打印初始值
    print('epoch = 0, x = %f, f(x) = %f' % (x0, func(x0)))
    for n in range(max_iter):
        # 利用公式求取新的迭代点坐标x
        x[n + 1] = x[n] - (func(x[n]) / dfunc(x[n]))
        # 将迭代点x坐标，代入函数得到y
        y = func(x[n + 1])
        step = n + 1
        print('epoch = %d, x = %f, f(x) = %f' % (n + 1, x[n + 1], y) )
        if abs(y) <= eps:
            break
    # 返回零点值
    return x[step]

def func(x):
    return math.pow(x, 3) - 2 * math.pow(x, 2) - 11 * x + 12


def dfunc(x):
    return  3 * math.pow(x, 2) - 4 * x - 11


if __name__ == '__main__':
    x0 = 2.35284172
    # x0 = 2.35287525
    fun(x0, func, dfunc)

• 当初始值为2.35284172时候，得到epoch = 32, x = -3.000000, f(x) = 0.000000
• 当初始值为2.35287525时候，得到epoch = 29, x = 4.000000, f(x) = 0.000000
• 我们可以画出函数图像，发现g(x)的零点有多个。

2.4 利用牛顿法求极值

如果函数可微，我们想找到局部极值点，只要去找导函数为零的位置就可以了。

由此，求解极值的问题就可以转化为寻找导函数零点的问题。

2.4.1 利用牛顿法求一元极值点

$$\begin{gathered} 假如我们要求f(x)=\frac{1}{4}{x}^4-\frac{1}{8}x的极值点，已经知道存在极小值，假如现在想求极小值所对应的位置求出来。 \\ 很简单，求导函数就行。g(x)=f^{\prime }(x)=x^3-\frac{1}{8} \\ 这样就回到了「用牛顿法求零点」的方法中。思路很简单，分成两步走。 \\ 第一步、用迭代点和方程导函数确定切线 \\ 第二步、用切线与x轴的交点计算函数值 \\ 迭代公式如下： \\ x=x^{(k)}-\frac{g(x(k))}{g^{\prime }(x(k))} \\ 用原函数f(x)表示如下： \\ x=x^{(k)}-\frac{f^{\prime }(x(k))}{f^{\prime \prime }(x(k))} \\ 利用上面公式就可以很容易的求出f(x)的极小值。 \end{gathered}$$

用牛顿法求极值的算法：

例题求解：

按照算法步骤我们求得：
$$\begin{gathered} g(x)=f^{\prime }(x)=x^3-\frac{1}{8} \\ {g}^{\prime }(x)=f^{\prime \prime }(x)=3x^2 \end{gathered}$$
我们利用代码求解

import math


def fun(x0, func, dfunc, ddfunc):
    eps = 0.0001
    max_iter = 500
    step = 0
    x = [0] * (max_iter + 1)
    x[0] = x0
    # 打印初始值
    print('epoch = 0, x = %f, |g(x)| = %f, f(x) = %f' % (x0, abs(dfunc(x[0])), func(x0)))
    for n in range(max_iter):
        # 利用公式求取新的迭代点坐标x
        x[n + 1] = x[n] - (dfunc(x[n]) / ddfunc(x[n]))
        # 将迭代点x坐标，代入函数得到y
        y = func(x[n + 1])
        step = n + 1
        print('epoch = %d, x = %f, |g(x)| = %f,  f(x) = %f' % (n + 1, x[n + 1], abs(dfunc(x[n + 1])), y) )
        if abs(dfunc(x[n + 1])) <= eps:
            break
    # 返回零点值
    return x[step]

def func(x):
    return (1 / 4) * math.pow(x, 4) - (1 / 8) * x


def dfunc(x):
    return  math.pow(x, 3) - (1 / 8)


def ddfunc(x):
    return  3 * math.pow(x, 2)

if __name__ == '__main__':
    # 设定初始值为2
    x0 = 2.0
    fun(x0, func, dfunc, ddfunc)

epoch = 0, x = 2.000000, |g(x)| = 7.875000, f(x) = 3.750000
epoch = 1, x = 1.343750, |g(x)| = 2.301361,  f(x) = 0.647137
epoch = 2, x = 0.918909, |g(x)| = 0.650921,  f(x) = 0.063386
epoch = 3, x = 0.661951, |g(x)| = 0.165053,  f(x) = -0.034744
epoch = 4, x = 0.536391, |g(x)| = 0.029328,  f(x) = -0.046354
epoch = 5, x = 0.502413, |g(x)| = 0.001819,  f(x) = -0.046873
epoch = 6, x = 0.500012, |g(x)| = 0.000009,  f(x) = -0.046875

迭代停止，求出了对应的极小值点为0.500012时对应的函数值。

2.4.2 多元函数求极值点

那么如果是要求多元的，也就是涉及到多个参数该怎么办呢？
$$\begin{gathered} 已知f(x)=f(x_1,x_2,...,x_n)，我们求解f(x)对向量\overrightarrow{x}=\left(\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right)的导数 \\ 按照分子布局，可以求得\frac{\partial f(x)}{\partial \overrightarrow{x}}=\left(\begin{array}{c}\frac{\partial f(x)}{\partial x_1},\frac{\partial f(x)}{\partial x_2},\cdots ,\frac{\partial f(x)}{\partial x_n}\end{array}\right) \\ 假如我们对此导数，再次求导，即二阶导数： \\ 就得到了海森矩阵H_f=\left(\begin{array}{cccc}\frac{{\partial }^{2}f(x)}{\partial x_1\partial x_1}& \frac{{\partial }^{2}f(x)}{\partial x_1\partial x_2}& \ddots & \frac{{\partial }^{2}f(x)}{\partial x_1\partial x_n}\\ \frac{{\partial }^{2}f(x)}{\partial x_2\partial x_1}& \frac{{\partial }^{2}f(x)}{\partial x_2\partial x_2}& \cdots & \frac{{\partial }^{2}f(x)}{\partial x_2\partial x_n}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{{\partial }^{2}f(x)}{\partial x_n\partial x_1}& \frac{{\partial }^{2}f(x)}{\partial x_n\partial x_2}& \cdots & \frac{{\partial }^{2}f(x)}{\partial x_n\partial x_n}\end{array}\right) \end{gathered}$$
用牛顿法求解多元极值问题：

例题求解：
$$\begin{gathered} f(x_1,x_2)=5x_1^4+4x_1^2{x}_2-x_1{x}_2^3+4x_2^4-x_1 \\ f(x_1,x_2)梯度=\left(\begin{array}{c}\frac{\partial f(x)}{\partial x_1},\frac{\partial f(x)}{\partial x_2}\end{array}\right)=\left(\begin{array}{c}20x_1^3+8x_1{x}_2-x_2^3-1,4x_1^2-3x_1{x}_2^2+16x_2^3\end{array}\right) \\ 求解海森矩阵H_f=\left(\begin{array}{cc}\frac{{\partial }^{2}f(x)}{\partial x_1\partial x_1}& \frac{{\partial }^{2}f(x)}{\partial x_1\partial x_2}\\ \frac{{\partial }^{2}f(x)}{\partial x_2\partial x_1}& \frac{{\partial }^{2}f(x)}{\partial x_2\partial x_2}\end{array}\right)=\left(\begin{array}{cc}60x_1^2+8x_2& 8x_1-3x_2^2\\ 8x_1-3x_2^3& -6x_1{x}_2+48x_2^2\end{array}\right) \\ 多元矩阵的迭代公式如下： \\ {x}^{(k+1)}=x^{(k)}-\nabla f(x_1^k,x_2^k)H_f^{-1}(x_1^k,x_2^k) \end{gathered}$$
下面我们就可以采用之前所说的步骤进行迭代计算。

我们这里使用代码实现：

import math
import numpy as np

def fun(x10, x20, func, dfunc, ddfunc):
    eps = 0.0001
    max_iter = 500
    step = 0
    # 初始化
    x = np.zeros((max_iter + 1, 2))
    x[0][0] = x10
    x[0][1] = x20
    # 保留历次计算的f(x)值
    y_arr = [0] * (max_iter + 1)
    y_arr[0] = func(x10, x20)
    # 打印初始值
    print('epoch = 0, x1 = %.4f, x1 = %.4f, f(x) = %.4f' % (x10, x20 , func(x10, x20) ))
    for n in range(max_iter):
        # 利用公式求取新的迭代点坐标x
        g_f = dfunc(x[n][0], x[n][1])    # 梯度
        h_f_i = ddfunc(x[n][0], x[n][1]) # 海森矩阵的逆
        eta = np.dot(g_f, h_f_i)
        x[n + 1] = x[n] - eta
        # 将迭代点x坐标，代入函数得到y
        x1 = x[n + 1][0]
        x2 = x[n + 1][1]
        y = func(x1, x2)
        y_arr[n + 1] = y
        step = n + 1
        diff = abs(y - y_arr[n])
        print('epoch = %d, x1 = %.4f, x1 = %.4f, f(x) = %.4f, |f(x1, x2)_k+1 - f(x1, x2)_k| = %.4f' % ( n + 1, x1, x2, func(x1, x2), diff))
        if diff <= eps:
            break
    # 返回零点值
    return x[step]


def func(x1, x2):
    # 原函数
    return 5 * math.pow(x1, 4) + 4 * math.pow(x1, 2) * x2 - x1 * math.pow(x2, 3) + 4 * math.pow(x2, 4) - x1


def dfunc(x1, x2):
    # 求解梯度
    return  np.array([
        20 * math.pow(x1, 3) + 8 * x1 * x2 - math.pow(x2, 3) - 1,
        4 * math.pow(x1, 2) - 3 * x1 * math.pow(x2, 2) + 16 * math.pow(x2, 3)
        ])


def ddfunc(x1, x2):
    H_f = np.mat(
        [
            [60*x1*x1 + 8 * x2, 8*x1 - 3 * x2 * x2],
            [8*x1 - 3 * x2 * x2, -6 * x1 * x2 + 48 * x2 * x2]
        ]
    )
    return H_f.I  # 求解海森矩阵的逆

if __name__ == '__main__':
    x10 = 1.0
    x20 = 1.0
    fun(x10, x20 , func, dfunc, ddfunc)

epoch = 0, x1 = 1.0000, x1 = 1.0000, f(x) = 11.0000
epoch = 1, x1 = 0.6443, x1 = 0.6376, f(x) = 1.7700, |f(x1, x2)_k+1 - f(x1, x2)_k| = 9.2300
epoch = 2, x1 = 0.4306, x1 = 0.3923, f(x) = 0.1011, |f(x1, x2)_k+1 - f(x1, x2)_k| = 1.6689
epoch = 3, x1 = 0.3388, x1 = 0.1986, f(x) = -0.1782, |f(x1, x2)_k+1 - f(x1, x2)_k| = 0.2793
epoch = 4, x1 = 0.5001, x1 = -0.4477, f(x) = -0.4296, |f(x1, x2)_k+1 - f(x1, x2)_k| = 0.2515
epoch = 5, x1 = 0.4974, x1 = -0.3797, f(x) = -0.4567, |f(x1, x2)_k+1 - f(x1, x2)_k| = 0.0271
epoch = 6, x1 = 0.4926, x1 = -0.3650, f(x) = -0.4575, |f(x1, x2)_k+1 - f(x1, x2)_k| = 0.0008
epoch = 7, x1 = 0.4923, x1 = -0.3643, f(x) = -0.4575, |f(x1, x2)_k+1 - f(x1, x2)_k| = 0.0000

可以看出到第7次的时候，差值结果小于精度阈值，意味着找到了极小值点。

注意：

当然这里的难点就在于要求出逆的海森矩阵，这里的例子因为比较简单可以求出，但是如果数值计算量过大，就不容易啦，那就要使用新的方法-拟牛顿法了。

三、拟牛顿法

• 用牛顿法来求多元极值点时，需要用到海森矩阵。

• 难点就在于要求出海森矩阵的逆，如果数值计算量过大，由于是二阶偏导，那在计算过程中就会复杂很多，同时，如果海森矩阵的行列式为零 ，则计算不出逆矩阵，因此这时就需要拟牛顿法

• 拟牛顿法的核心就是找到替代品来取代求逆的海森矩阵

• 拟牛顿法的思路是不计算目标函数的Hessian矩阵然后求逆矩阵，而是通过其他手段得到一个近似Hessian矩阵逆的矩阵。具体做法是构造一个近似Hessian矩阵或其逆矩阵的正定对称矩阵，用该矩阵进行牛顿法的迭代。

3.1 泰勒展开公式

• 我们知道梯度下降法是通过一阶泰勒展开得到的。

• 其实，牛顿法的实质其实是对它的目标函数$f(x)$进行二阶泰勒展开，这里选择的点就是$x^k$。

泰勒展开公式如下：

$$\begin{gathered} 我们对f(x)在x^{(k)}进行二阶泰勒展开(这里写「约等号」原因是没有写上后面的余项)： \\ f(x)\approx f(x^{(k)})+g_k^T(x-x^{(k)})+\frac{1}{2}(x-x^{(k)}{)}^{T}H(x^{(k)})(x-x^{(k)}) \\ 接着我们对目标函数求导数： \\ \frac{\partial f(x)}{\partial x}\approx 0+g_k(x^{(k)})+\frac{1}{2}\ast 2H(x^{(k)})(x-x^{(k)}) \\ =g_k(x^{(k)})+H(x^{(k)})(x-x^{(k)}) \end{gathered}$$
在式子的左边恰好就是梯度函数$g(x)$ ，它代表了向量，而不是数值。

「第一项」 $f(x^{(k)})$其实是指第k次迭代值后计算出的**「具体的数值」**，求导为零。

「第二项」 包含了x ，求导后就是在第k个迭代点处的梯度向量。

「第三项」 相当于是一个向量平方的形式。向量平方$z^{T}Az$求导得到$(A+A^T)z$，如果矩阵A是对称的，那么就与数值平方求导类似，为$2Az$。(不熟悉矩阵求导的可以参考：矩阵的导数运算(理解分子布局、分母布局))

3.2 逆牛顿法的条件

$$\begin{gathered} g(x)=\frac{\partial f(x)}{\partial x}=g(x^{(k)})+H(x^{(k)})(x-x^{(k)}) \\ 我们令x=x^{(k+1)}，那么 \\ g(x^{(k+1)})=g(x^{(k)})+H(x^{(k)})(x^{(k+1)}-x^{(k)}) \\ 可以写成：g(x^{(k+1)})-g(x^{(k)})=H(x^{(k)})(x^{(k+1)}-x^{(k)}) \end{gathered}$$

• 「左式」是两个导函数的差值
• 「右式」是两个相邻迭代点之间的差值。

$$\begin{gathered} 令左式g_{k+1}-g_k=y_k,右式：x^{(k+1)}-x^{(k)}={\delta }_k \\ 上面式子可以写成： \\ {y}_k=H_k{\delta }_k或者H_k^{-1}{y}_k={\delta }_k \end{gathered}$$

• 这两个式子就是拟牛顿法的条件。
• 这两个条件这就是拟牛顿法的精髓。

3.3 拟牛顿法DFP算法

• DFP算法是三个人名的首字母，分别对应的是 Davidon、Fletcher、Powell，这个算法就是由他们三人得到的

• DFP算法的精髓是找到求解海森矩阵的逆的替代品

• DFP算法其实就是条件$H_k^{-1}{y}_k={\delta }_k$推导而来的

3.3.1 DFP算法的推导过程

讲解一下第二步时，$P_k$矩阵的构造：
$$\begin{gathered} 在求让P_k{y}_k={\delta }_k成立的P_k时候，不能直接做除法，即P_k=\frac{{\delta }_k}{y_k} \\ 这是因为{\delta }_k和y_k都是向量，向量不能直接相除。 \\ 假设：此时对于x有3个维度，当k=1时，那么意味着计算出的{\delta }_1便是一个(3,1)的矩阵, \\ 此时y_k也是一个(3,1)的矩阵，代入到构造的矩阵P_k=\frac{{\delta }_k{\delta }_k^T}{{\delta }_k^T{y}_k}中 \\ 此时分子就是(3,3)的矩阵，分母就变成了一个常数，因此P_k是一个(3,3)的矩阵 \end{gathered}$$

3.3.2 DFP算法的迭代过程

3.3.3 DFP算法代码实现

我们利用DFP算法来求解$f(x)=100(x_1^2-x_2^2{)}^2+(x_1-1{)}^2$的极小值点。

"""
Newton法求二元函数
"""
import numpy as np
from sympy import symbols

# 首先定义二维向量内的元素
x1 = symbols("x1")
x2 = symbols("x2")


def jacobian(f, x):
    """
    求函数一阶导
    :param f: 原函数
    :param x: 初始值
    :return: 函数一阶导的值
    """
    grandient = np.array([400 * x[0] * (x[0] ** 2 - x[1]) + 2 * (x[0] - 1),
                          -200 * (x[0] ** 2 - x[1])], dtype=float)
    return grandient

def dfp_newton(f, x, iters):
    """
    实现DFP拟牛顿算法
    :param f: 原函数
    :param x: 初始值
    :param iters: 遍历的最大epoch
    :return: 最终更新完毕的x值
    """
    # 步长
    learning_rate = 1
    # 初始化正定矩阵
    G = np.eye(2)
    x_len = x.shape[0]
    # 一阶导g的第二范式的最小值（阈值）
    epsilon = 1e-5
    for i in range(1, iters):
        g = jacobian(f, x)
        if np.linalg.norm(g) < epsilon:
            break
        p = np.dot(G, g)
        # 更新x值
        x_new = x - p * learning_rate
        print("第" + str(i) + "次迭代后的结果为:", x_new)
        g_new = jacobian(f, x_new)
        y = g_new - g
        k = x_new - x
        Gy = np.dot(G, y)
        y_t_G = np.dot(y, G)
        yGy = np.dot(np.dot(y, G), y)
        # 更新G正定矩阵
        G = G + k.reshape([x_len, 1]) * k / np.dot(k, y) - Gy.reshape([x_len, 1]) * y_t_G / yGy
        x = x_new
    return x


if __name__ == "__main__":
    x = np.array([1, 9], dtype=float)
    f = 100 * (x1 ** 2 - x2 ** 2) ** 2 + (x1 - 1) ** 2
    print(dfp_newton(f, x, 1000))

3.4 拟牛顿法BFGS算法

• 和DFP一样，BFGS算法由是 Broyden、Fletcher、Goldfarb和Shanno创造。
• BFGS用到的是拟牛顿法中的条件$y_k=H_k{\delta }_k$推导而来的。
• BFGS算法不再是找海森矩阵的逆的替代品，而是直接是海森矩阵的替代品，即$H_k=B_k$

3.4.1 BFGS算法的推导过程

3.4.2 BFGS算法的迭代过程

3.4.3 BFGS算法代码实现

我们利用BFGS算法来求解$f(x)=100(x_1^2-x_2^2{)}^2+(x_1-1{)}^2$的极小值点。

"""
Newton法求二元函数
"""
import numpy as np
from sympy import symbols

# 首先定义二维向量内的元素
x1 = symbols("x1")
x2 = symbols("x2")


def jacobian(f, x):
    """
    求函数一阶导
    :param f: 原函数
    :param x: 初始值
    :return: 函数一阶导的值
    """
    grandient = np.array([400 * x[0] * (x[0] ** 2 - x[1]) + 2 * (x[0] - 1),
                          -200 * (x[0] ** 2 - x[1])], dtype=float)
    return grandient


def bfgs_newton(f, x, iters):
    """
    实现BFGS拟牛顿法
    :param f: 原函数
    :param x: 初始值
    :param iters: 遍历的最大epoch
    :return: 最终更新完毕的x值
    """
    # 步长。设为1才能收敛，小于1不能收敛
    learning_rate = 1
    # 初始化B正定矩阵
    B = np.eye(2)
    x_len = x.shape[0]
    # 一阶导g的第二范式的最小值（阈值）
    epsilon = 1e-5
    for i in range(1, iters):
        g = jacobian(f, x)
        if np.linalg.norm(g) < epsilon:
            break
        p = np.linalg.solve(B, g)
        # 更新x值
        x_new = x - p*learning_rate
        print("第" + str(i) + "次迭代后的结果为:", x_new)
        g_new = jacobian(f, x_new)
        y = g_new - g
        k = x_new - x
        y_t = y.reshape([x_len, 1])
        Bk = np.dot(B, k)
        k_t_B = np.dot(k, B)
        kBk = np.dot(np.dot(k, B), k)
        # 更新B正定矩阵。完全按照公式来计算
        B = B + y_t*y/np.dot(y, k) - Bk.reshape([x_len, 1]) * k_t_B / kBk
        x = x_new
    return x



if __name__ == "__main__":
    x = np.array([1, 9], dtype=float)
    f = 100 * (x1 ** 2 - x2 ** 2) ** 2 + (x1 - 1) ** 2
    print(bfgs_newton(f, x, 2000))