Newton及阻尼Newton法（算法分析+python代码解释）

背对人潮

已于 2023-11-12 15:06:41 修改

阅读量1.6k

点赞数 1

分类专栏：最优化问题 # 无约束最优化方法文章标签：算法机器学习 python pycharm 1024程序员节

于 2023-10-31 14:20:03 首次发布

本文链接：https://blog.csdn.net/m0_61209712/article/details/134137865

版权

最优化问题同时被 2 个专栏收录

24 篇文章 13 订阅

订阅专栏

无约束最优化方法

9 篇文章 0 订阅

订阅专栏

前言：

由于最速下降法的搜索路径在靠近最优解附近呈锯齿状，所以搜索很慢，究其原因，应该是与搜索方向由目标函数的线性近似函数来决定有关。事实上，由一阶泰勒公式，有：

$f(x^{k}+\lambda _{k}p^{k})=f(x^{k})+\lambda _{k}\bigtriangledown f(x^{k})^{T}p^{k}+o(\left | \left |\lambda _{k}p^{k} \right | \right |)\approx f(x^{k})+\lambda _{k}\bigtriangledown f(x^{k})^{T}p^{k}$

由于 $p^{k}=-\bigtriangledown f(x^{k})$ ,所以 $f(x^{k}+\lambda _{k}p^{k})-f(x^{k})\approx -\lambda _{k}\left | \left | \bigtriangledown f(x^{k}) \right | \right |^{2}< 0$ ，这就是说在最速下降法中目标函数下降的幅度其实是由线性近似多项式函数的下降幅度来决定的，由此我们考虑，如果能利用函数的二次近似多项式来做目标函数的近似表达式，或许算法的收敛速度可以提高。因此，对于有较好的解析性质的目标函数，如具有二阶连续偏导数的目标函数，可以利用其二次近似函数的性质来界定新的搜索方向，或者以目标函数的二次近似函数的极小点作为新目标函数本身的近似解，如果不满足计算精度，则再求该点处目标函数的二次近似函数，进一步求新的近似函数的极小点，如此反复，直至近似点满足事先给定的精度为止，这就是牛顿法的基本思想。

Newton法：

算法步骤：

步骤1：选定初始点 $x^{1}$ ，计算 $f_{1}=f(x_{1})$ ,k=1；

步骤2：如果 $\left | \left | \bigtriangledown f(x^{k}) \right | \right |\leq \epsilon$ ，算法停止， $x^{*}=x^{k}$ ,否则，转步骤3；

步骤3：计算搜索方向 $d^{k}=-(\bigtriangledown ^{2}f(x^{k}))^{-1}\bigtriangledown f(x^{k})$

步骤4：令 $x^{k+1}=x^{k}+d^{k}$ ，k=k+1,转步骤2

（说明：搜索方向 $d^{k}$ 可通过 $\bigtriangledown ^{2}f(x^{k})d^{k}+\bigtriangledown f(x^{k})=0$ 求出，步骤4中定的步长为1）

例题分析：

下面我们通过一道例题来加深一下对算法的认识。

例题：用Newton法求 $f(x_{1},x_{2})=x_{1}^{2}+25x_{2}^{2}$ 的极小点。

解：取初始点 $x^{1}=(2,2)^{T}$

计算：

$\bigtriangledown f(x^{1})=\binom{2x_{1}}{50x_{2}}|_{x^{1}}=\binom{4}{100}$

$\bigtriangledown ^{2}f(x^{1})=\begin{pmatrix} 2 & 0\\ 0 & 50 \end{pmatrix}$

代入Newton迭代公式：

$x^{2}=x^{1}-(\bigtriangledown ^{2}f(x^{1}))^{-1}\bigtriangledown f(x^{1})=\binom{2}{2}-\begin{pmatrix} \frac{1}{2} &0 \\0 &\frac{1}{50} \end{pmatrix}\binom{4}{100}=\binom{0}{0}$

得到问题的极小点 $x^{*}=x^{2}$ .

结果分析：

由以上例题可以看出过程中一步迭代就达到了最优解，经过推论，我们发现，当f为正定二次函数时，用Newton法从任意初始点可一步达到极小点。

分析：

Newton法的优点在于当初始点离极小点很近时，算法收敛速度快；算法简单，不需要进行一维搜索，对正定二次函数迭代一次就可得到极小点。但是缺点也较为明显：当初始点远离最优解时，算法可能不收敛；Hesse矩阵非奇异时，算法不可行；计算量和存储量大。

那么我们下一步介绍改进后的Newton法——阻尼Newton法。

收敛性分析：

假设目标函数具有二阶连续偏导数，并且 Hesse 矩阵在最优值点 $\bar{x}$ 的一个领域 $N_{\delta }(\bar{x})$ 内是利普希茨连续的，即存在常熟 L> 0 ，使得 $\left | \left | \bigtriangledown ^{2}f(x)-\bigtriangledown ^{2}f(y) \right | \right |\leq L\left | \left | x-y \right | \right |$

$\forall x,y\in N_{\delta }(\bar{x})$ ，如果 f(x) 在点 $\bar{x}$ 处满足 $\bigtriangledown f(\bar{x})=0$ ， $\bigtriangledown^{2} f(\bar{x})$ 正定，那么对于牛顿法的迭代式具有一下结论：

（1）如果初始点 $x_{0}$ 与最优点 $\bar{x}$ 足够近，则牛顿法产生的点列 { $x_{k}$ } 收敛于 $\bar{x}$ 。

（2）点列 { $x_{k}$ } 收敛到 $\bar{x}$ 的速度是二阶的。

阻尼Newton法：

在Newton的基础上，将迭代公式 $x^{k+1}=x^{k}-(\bigtriangledown ^{2}f(x^{k}))^{-1}\bigtriangledown f(x^{k})$ ，加入精确一维搜索： $minf(x^{k}+\lambda d^{k})$ ，求得最佳步长 $\lambda _{k}$ ，得到下一个迭代点 $x^{k+1}=x^{k}+\lambda _{k}d^{k}$ ，由此就可以改进Newton法的收敛性。

例题分析：

下面我们通过一道例题来加深一下对算法的认识。

例题：用阻尼Newton法求下列无约束优化的极小点，已知 $f(x)=x_{1}^{2}+2x_{2}^{2}-2x_{1}x_{2}-4x_{1}$ ， $x^{1}=(1,1)^{T},\varepsilon =0.1$ .

解：

$\bigtriangledown f(x^{1})=\binom{2x_{1}-2x_{2}-4}{-2x_{1}+4x_{2}}|_{x^{1}}=\binom{-4}{2}$

$\bigtriangledown ^{2}f(x^{1})=\begin{pmatrix} 2 & -2\\ -2 & 4 \end{pmatrix}$

故 $d^{1}=-(\bigtriangledown ^{2}f(x^{1}))^{-1}\bigtriangledown f(x^{1})=-\begin{pmatrix} 2 &-2 \\-2 &4 \end{pmatrix}^{-1}\binom{-4}{2}=\binom{3}{1}$

$x^{1}=\binom{1}{1}$ ， $d^{1}=\binom{3}{1}$ ， $x^{1}+\lambda d^{1}=\binom{1+3\lambda }{1+\lambda }$

$\phi (\lambda )=f(x^{1}+\lambda d^{1})=f(1+3\lambda ,1+\lambda )=(1+3\lambda )^{2}+2(1+\lambda )^{2}-2(1+3\lambda )(1+\lambda )-4(1+3\lambda )$

令 $\phi '(\lambda )=0\Rightarrow \lambda _{1}=1$

$x^{2}=x^{1}+\lambda _{1}d^{1}=\binom{4 }{2 }$

$\bigtriangledown ^{2}f(x^{1})=\binom{0}{0}$

得到问题的极小点 $x^{*}=x^{2}$ = $\binom{4}{2}$

算法代码：

'''
阻尼Newton算法
2023.10.19
'''
import numpy as np
from sympy import *

x1 = symbols("x1")
x2 = symbols("x2")
λ = symbols("λ")  # 阻尼Newton算法，步长 λ
# f = x1**2 + 2*x2**2 - 2*x1*x2 - 4*x1
# f = x1**2 + 8*x2**2 - 2*x1*x2 + 4*x1 + 5
f = (x1 - x2) ** 2 + x2 ** 2 - 4 * x1


def gradient_descent(x1_init, x2_init, ε):
    # 一阶求导
    first_grad_1 = diff(f, x1)
    first_grad_2 = diff(f, x2)
    # 二阶求导
    second_grad_1 = diff(first_grad_1, x1)
    second_grad_2 = diff(first_grad_2, x2)
    second_grad_12 = diff(first_grad_1, x2)

    x1_curr = x1_init
    x2_curr = x2_init

    count = 0

    while True:
        count = count + 1
        first_grad_1_value = first_grad_1.subs({x1: x1_curr, x2: x2_curr})
        first_grad_2_value = first_grad_2.subs({x1: x1_curr, x2: x2_curr})
        second_grad_1_value = second_grad_1.subs({x1: x1_curr, x2: x2_curr})
        second_grad_2_value = second_grad_2.subs({x1: x1_curr, x2: x2_curr})
        second_grad_12_value = second_grad_12.subs({x1: x1_curr, x2: x2_curr})

        print("第{}次迭代".format(count - 1))
        print("此时的解为[%f, %f]" % (float(x1_curr), float(x2_curr)))

        arr = [[second_grad_1_value, second_grad_12_value], [second_grad_12_value, second_grad_2_value]]
        arr = np.array(arr, dtype=np.float64)
        arr = np.linalg.inv(arr)

        d1 = -1 * (arr[0, 0] * first_grad_1_value + arr[0, 1] * first_grad_2_value)
        d2 = -1 * (arr[1, 0] * first_grad_1_value + arr[1, 1] * first_grad_2_value)

        k = np.array([first_grad_1_value, first_grad_2_value], dtype=float)
        norm_result = np.linalg.norm(k, ord=2, axis=0)

        # x1_curr = x1_curr + d1
        # x2_curr = x2_curr + d2

        x1_curr = x1_curr + λ * d1
        x2_curr = x2_curr + λ * d2

        f_new = f.subs({x1: x1_curr, x2: x2_curr})

        grad_3 = diff(f_new, λ)
        if grad_3 == 0:
            print("该精度下的最优解是[%f, %f]" % (float(x1_curr), float(x2_curr)))
            break
        else:
            λ_value = solve(grad_3, λ)[0]
            print("λ的值为{}".format(λ_value))

            x1_curr = x1_curr.subs(λ, λ_value)
            x2_curr = x2_curr.subs(λ, λ_value)

            if norm_result < ε:
                print("该精度下的最优解是[%f, %f]" % (float(x1_curr), float(x2_curr)))
                break

    return x1_curr, x2_curr


result = gradient_descent(x1_init=1, x2_init=1, ε=0.1)

用该代码所运算的结果为：

第0次迭代
此时的解为[1.000000, 1.000000]
λ的值为1.00000000000000
第1次迭代
此时的解为[4.000000, 2.000000]
该精度下的最优解是[4.000000, 2.000000]

总结：

1）针对Hesse矩阵计算量的改进：

为减小工作量，取m（正整数），使每m次迭代使用同一个Hesse阵，迭代公式变成：

$x^{k_{m}+j+1}=x^{k_{m}+j}-(\bigtriangledown ^{2}f(x^{k_{m}}))^{-1}\bigtriangledown f(x^{k_{m}+j})$

2)针对Hesse矩阵非正定和奇异的情况

取: $d^{k}=\left\{\begin{matrix} -(\bigtriangledown ^{2}f(x^{k}))^{-1}\bigtriangledown f(x^{k}) \\-\bigtriangledown f(x^{k})& \end{matrix}\right.$ ,并采用下列非精确一维搜索（Armijo-Goldstein准则）求 $\lambda _{k}$ $f(x^{k}+\lambda _{k}d^{k})\leq f(x^{k})+\delta \lambda _{k}\bigtriangledown f(x^{k})^{T}d^{k}$

$f(x^{k}+\lambda _{k}d^{k})\geq f(x^{k})+(1-\delta )\lambda _{k}\bigtriangledown f(x^{k})^{T}d^{k}$ ,其中 $\delta \in (0,\frac{1}{2})$

（行文中若有纰漏，希望大家指正）。

背对人潮

关注

1
点赞
踩
15

收藏

觉得还不错? 一键收藏
打赏
1
评论
Newton及阻尼Newton法（算法分析+python代码解释）

因此，对于有较好的解析性质的目标函数，如具有二阶连续偏导数的目标函数，可以利用其二次近似函数的性质来界定新的搜索方向，或者以目标函数的二次近似函数的极小点作为新目标函数本身的近似解，如果不满足计算精度，则再求该点处目标函数的二次近似函数，进一步求新的近似函数的极小点，如此反复，直至近似点满足事先给定的精度为止，这就是牛顿法的基本思想。
复制链接

扫一扫