共轭梯度法（算法分析+python代码解释）

背对人潮

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

阅读量5k

点赞数 13

分类专栏：最优化问题 # 无约束最优化方法文章标签：算法 python pycharm 线性代数机器学习

于 2023-11-03 11:45:25 首次发布

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

版权

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

24 篇文章 13 订阅

订阅专栏

无约束最优化方法

9 篇文章 0 订阅

订阅专栏

前言：

为了克服最速下降法在最优解附近收敛速度变慢和Newton法不具有整体收敛性以及需要计算Hesse矩阵的缺点，一种介于两者之间的方法被提出，这就是共轭方向法。这种方法不仅克服了上述缺点之外，还能保持上述两种方法各自的优点，即收敛性较好，并且可以在有限布内使二次函数达到最优，即具有二次终止性。

共轭方向法：

基本原理：假设 $f$ 是 n 元正定二次函数 $f(x)=\frac{1}{2}x^{T}Ax+b^{T}+c$ , A 正定，如果按最速下降法，则会发生锯齿现象。希望下次迭代就得到最优解，如果能够选定这样的搜索方向，那么对于二元二次函数，只需沿搜索方向 $d^{1},d^{2}$ 进行精确一维搜索就可以求到极小点 $x^{*}$ ,即 $x^{*}=x^{2}+\lambda _{2}d^{2}$ 。

$x^{*}$ 是 $f$ 的极小点，故 $x^{*}$ 是 $f$ 的驻点， $\bigtriangledown f(x^{*})=Ax^{*}+b=0$

$\Rightarrow 0=\bigtriangledown f(x^{*})=\bigtriangledown f(x^{2}+\lambda _{2}d^{2})$

$=A(x^{2}+\lambda _{2}d^{2})+b$

$=(Ax^{2}+b)+A\lambda _{2}d^{2}$

$=\bigtriangledown f(x^{2})+\lambda _{2}Ad^{2}$

等式两边同时左乘 $(d^{1})^{T}$ ,有：

$0=(d^{1})^{T}\bigtriangledown f(x^{2})+(d^{1})^{T}\lambda _{2}Ad^{2}$

由于 $(d^{1})^{T}\bigtriangledown f(x^{2})$ 相互正交，所以相乘为 0，

$\Rightarrow (d^{1})^{T}Ad^{2}=0$

$\Rightarrow d^{1},d^{2}$ 是 A 的共轭矩阵

共轭方向及其性质：

性质1： $R^{n}$ 中与 n 个线性无关的向量都正交的一定是零向量。

性质2： $R^{n}$ 中 A 共轭的非零向量组 $p^{1},p^{2},......p^{n}$ 是线性无关的。

性质3： $R^{n}$ 中互相共轭的非零向量的个数不超过 n 。

性质4：给定 n 元函数 $f(x)=\frac{1}{2}x^{T}Ax+b^{T}+c$ , $A=A^{T}$ 正定，设 n 维非零向量组 $p^{1},p^{2},......p^{n}$ 是 A 共轭向量组，从任意点 $x^{1}$ 出发，依次以 $p^{1},p^{2},......p^{n}$ 为搜索方向进行精确一维搜索，则

（1） $\bigtriangledown f(x^{k+1})$ 与 $p^{1},p^{2},......p^{n} (k=1,2,...,n)$ 正交。

（2）最多 n 次迭代必达到二次函数 $f$ 的极小点。

ok，了解了什么是共轭方向法，那接下来我们来看看什么是共轭梯度法！

共轭梯度法：

令 $f(x)=\frac{1}{2}x^{T}Ax+b^{T}+c$ , $A=A^{T}$ 正定，

（1）从任取初始点 $x^{1}$ 出发，沿负梯度方向进行精确一维搜索：

$p^{1}=-\bigtriangledown f(x^{1}),x^{2}=x^{1}+\lambda _{1}p^{1}$

(2)若 $\bigtriangledown f(x^{2})=0$ ，停止，否则在 $-\bigtriangledown f(x^{2})$ 和 $p^{1}$ 张成的正交锥中找一个向量 $p^{2}$ ，即令 $p^{2}=-\bigtriangledown f(x^{2})+\alpha p^{1}$ ，使得 $(p^{2})^{T}Ap^{1}=0$

$(p^{2})^{T}Ap^{1}=0\Rightarrow (-\bigtriangledown f(x^{2})+\alpha p^{1})Ap^{1}=0$

$\Rightarrow \alpha _{1}=\frac{\bigtriangledown f(x^{2})^{T}Ap^{1}}{(p^{1})^{T}Ap^{1}}$

（3）在 $x^{2}$ 处沿 $p^{2}$ 方向进行精确一维搜索：

$x^{3}=x^{2}+\lambda _{2}p^{2}$

（4）以此类推；

（5）若 $\bigtriangledown f(x^{k+1})=0$ ，停止，否则在 $-\bigtriangledown f(x^{k+1})$ 和 $p^{k}$ 张成的正交锥中找一个向量 $p^{k+1}$ ，即令 $p^{k+1}=-\bigtriangledown f(x^{k+1})+\alpha_{k} p^{k}$ ，使得 $(p^{k+1})^{T}Ap^{k}=0$

$(p^{k+1})^{T}Ap^{k}=0\Rightarrow (-\bigtriangledown f(x^{k+1})+\alpha _{k}p^{k})Ap^{k}=0$

$\Rightarrow \alpha _{k}=\frac{\bigtriangledown f(x^{k+1})^{T}Ap^{k}}{(p^{k})^{T}Ap^{k}}$

这样就构造了一组向量组 $p^{1},p^{2},......p^{n}$ ，为A的共轭向量组。

由此，可引申出以下定理。

定理1：

设向量组 $p^{1},p^{2},......p^{n}$ 是由上述方法产生的向量组，向量组 $g_{1},g_{2},...,g_{n}$ 是由各点的梯度生成的向量组（ $g_{k}=\bigtriangledown f(x^{k})$ ） ,则

（1） $p^{1},p^{2},......p^{n}$ 是正交向量组；

（2） $g_{1},g_{2},...,g_{n}$ 是 A 的共轭向量组。

（注：为保证方向的共轭性，初始方向取负梯度方向）

定理2：

令 $f(x)=\frac{1}{2}x^{T}Ax+b^{T}+c$ , $A=A^{T}$ 正定，向量组 $p^{1},p^{2},......p^{n}$ 是由上述方法构造的 A 的共轭向量组， $g_{k}=\bigtriangledown f(x^{k})$ ，则以下几个计算公式等价：

（1） $\alpha _{k}=\frac{\bigtriangledown f(x^{k+1})^{T}Ap^{k}}{(p^{k})^{T}Ap^{k}}=\frac{(g_{k+1})^{T}Ap^{k}}{(p^{k})^{T}Ap^{k}}$ （Daniel共轭梯度法）

（2） $\alpha _{k}=\frac{(g_{k+1})^{T}(g_{k+1}-g_{k})}{(p^{k})^{T}(g_{k+1}-g_{k})}$ （SW共轭梯度法）

（3） $\alpha _{k}=\frac{(g_{k+1})^{T}g_{k+1}}{(p^{k})^{T}g_{k}}$ （DM共轭梯度法）

（4） $\alpha _{k}=\frac{\left | \left | g_{k+1} \right | \right |^{2}}{\left | \left | g_{k} \right | \right |^{2}}$ （FR共轭梯度法）

（5） $\alpha _{k}=\frac{(g_{k+1})^{T}(g_{k+1}-g_{k})}{(g_{k})^{T}g_{k}}$ （PPR共轭梯度法）

下面可以用一个反例来理解：

例题：用共轭方向法求下列问题的极值。

$f(x)=x_{1}^{2}+\frac{1}{2}x_{2}^{2}+\frac{1}{2}x_{3}^{2}$ ,已知初始点 $x^{1}=\begin{pmatrix} 1\\1 \\1 \end{pmatrix}$ ， $p^{1}=\begin{pmatrix} -1\\-2 \\0 \end{pmatrix}$ .

解：

由于初始方向为下降方向，但不是负梯度方向，所得出的 $p^{1},p^{2},p^{3}$ 并不是 A 的共轭向量组，导致 n 元二次函数最多 n 次迭代即可达到极小点的结论不成立。

收敛性分析：

设 $f$ ： $R^{n}\rightarrow R$ 具有一阶连续偏导数， $x_{0}\in R$ ，记 $\alpha = f(x_{0})$ ，假设水平集 $L(f,\alpha )$ 有界，令 { ${x_{k}}$ } 是由最速下降法产生的点列，则：

（1）当 { ${x_{k}}$ } 是有穷点列时，其中最后一个点是 $f$ 的驻点；

（2）当 { ${x_{k}}$ } 是无穷点列时，它必有极限点，并且任一极限点都是 $f$ 的驻点。

FR共轭梯度法：

FR共轭梯度法是共轭梯度法中最为常见，也是用的最多的一个，那么，下面我们就来展开谈谈FR共轭梯度法的具体做法。

算法步骤：

步骤1：选定初始点 $x^{1}$

步骤2：如果 $\left | \left | g_{1} \right | \right |\leq \varepsilon$ ，算法停止， $x^{*}=x^{1}$ ，否则转步骤3.

步骤3：取 $p^{1}=-g_{1},k=1$

步骤4：精确一维搜索找最佳步长 $\lambda _{k}$ ，令 $x^{k+1}=x^{k}+\lambda _{k}p^{k}$ .

步骤5：如果 $\left | \left | g_{k+1} \right | \right |\leq \varepsilon$ ，算法停止， $x^{*}=x^{k+1}$ ，否则转步骤6.

步骤6：如果 k=n ，令 $x^{1}=x^{K+1},p^{1}=-g_{k+1},k=1$ ，转步骤4；否则转步骤7.

步骤7：计算 $\alpha _{k}=\frac{\left | \left | g_{k+1} \right | \right |^{2}}{\left | \left | g_{k} \right | \right |^{2}}$ ， $p^{k+1}=-g_{k+1}+ \alpha _{k}p^{k}$ ，k=k+1，转步骤4.

（注：计算过程中存在一定的误差，会使 n 步迭代得不到正定二次函数的极小点， $R^{n}$ 中共轭方向最多有 n 个，n 步后构造的搜索方向不再是共轭的，会大大降低收敛速度，故需要重新开始，将 $x^{n+1}$ 作为新的 $x^{1}$ 。）

例题分析：

同样的，我们来用一道例题加深对该算法的理解

例题：用 FR 共轭梯度法求 $f(x)=x_{1}^{2}+2x_{2}^{2}-4x_{1}-2x_{1}x_{2}$ 的极小点，已知初始点 $(1,1)^{T}$ ,迭代精度 $\varepsilon$ =0.001.

解：

1）第一次迭代：沿负梯度方向搜索

$g_{1}=\bigtriangledown f(x^{1})=\binom{2x_{1}-4-2x_{2}}{4x_{2}-2x_{1}}=\binom{-4}{2}$ ， $p^{1}=-g_{1}=\binom{4}{-2}$

$x^{2}=x^{1}+\lambda p^{1}=\binom{1+4\lambda }{1-2\lambda }$

$\phi (\lambda )=(1+4\lambda )^{2}+2(1-2\lambda )^{2}-4(1+4\lambda )-2(1+4\lambda )(1-2\lambda )$

令 $\phi '(\lambda )=(40\lambda ^{2}-20\lambda -3)'=80\lambda -20=0$

$\Rightarrow \lambda =\frac{1}{4}$

$\therefore x^{2}=\binom{2}{\frac{1}{2}},g_{2}=\bigtriangledown f(x^{2})=\binom{-1}{-2}$

不满足精度，继续迭代。

2）第二次迭代：沿负梯度方向搜索

$\alpha _{1}=\frac{\left | \left | g_{2} \right | \right |^{2}}{\left | \left | g_{1} \right | \right |^{2}}=\frac{1}{4}$ ， $p^{2}=-g_{2}+\alpha _{1}p^{1}=\binom{2}{\frac{3}{2}}$

$x^{3}=x^{2}+\lambda p^{2}=\binom{2+2\lambda }{\frac{1}{2}+\frac{3}{2}\lambda }$

$\phi (\lambda )=(2+2\lambda )^{2}+2({\frac{1}{2}+\frac{3}{2}\lambda })^{2}-4(2+2\lambda )-2(2+2\lambda )({\frac{1}{2}+\frac{3}{2}\lambda })$

令 $\phi '(\lambda )=0$

$\Rightarrow \lambda =1$

$\therefore x^{3}=\binom{4}{2},g_{3}=\bigtriangledown f(x^{3})=\binom{0}{0}$

$\left | \left | g_{3} \right | \right |\leq \varepsilon$ ,得到最优解 $x^{*}=x^{3}=\binom{4}{2}$

算法代码：

'''
FR共轭方向算法(二元二次函数)
2023.10.20
'''
import numpy as np
from sympy import symbols, diff, solve

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


def fletcher_reeves(x1_init, x2_init, ε):
    # 一阶求导
    grad_1 = diff(f, x1)
    grad_2 = diff(f, x2)

    x1_curr = x1_init
    x2_curr = x2_init

    first_grad_1_value = grad_1.subs({x1: x1_curr, x2: x2_curr})
    first_grad_2_value = grad_2.subs({x1: x1_curr, x2: x2_curr})

    g1 = np.array([first_grad_1_value, first_grad_2_value], dtype=float)
    norm_result_1 = np.linalg.norm(g1, ord=2, axis=0)
    if norm_result_1 <= ε:
        print("算法停止，该精度下的最优解是[%f,%f]" % (float(x1_curr), float(x2_curr)))
    else:
        x1_new = x1_curr - λ * first_grad_1_value
        x2_new = x2_curr - λ * first_grad_2_value

        f_new = f.subs({x1: x1_new, x2: x2_new})
        grad_3 = diff(f_new, λ)
        λ_value = solve(grad_3, λ)[0]

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

        print("第{}次迭代".format(1), "当前最优解为[%f,%f]" % (float(x1_curr), float(x2_curr)))

        second_grad_1_value = grad_1.subs({x1: x1_curr, x2: x2_curr})
        second_grad_2_value = grad_2.subs({x1: x1_curr, x2: x2_curr})

        g2 = np.array([second_grad_1_value, second_grad_2_value], dtype=float)
        norm_result_2 = np.linalg.norm(g2, ord=2, axis=0)
        if norm_result_2 <= ε:
            print("算法停止，该精度下的最优解是[%f,%f]" % (float(x1_curr), float(x2_curr)))
        else:
            α1 = norm_result_2 ** 2 / norm_result_1 ** 2
            p2_1 = -1 * second_grad_1_value - α1 * first_grad_1_value
            p2_2 = -1 * second_grad_2_value - α1 * first_grad_2_value

            x1_new = x1_curr + λ * p2_1
            x2_new = x2_curr + λ * p2_2

            f_new = f.subs({x1: x1_new, x2: x2_new})
            grad_4 = diff(f_new, λ)
            λ_value = solve(grad_4, λ)[0]

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

            print("第{}次迭代".format(2), "当前最优解为[%f,%f]" % (float(x1_curr), float(x2_curr)))

            grad_1_value = grad_1.subs({x1: x1_curr, x2: x2_curr})
            grad_2_value = grad_2.subs({x1: x1_curr, x2: x2_curr})

            g3 = np.array([grad_1_value, grad_2_value], dtype=float)
            norm_result_3 = np.linalg.norm(g3, ord=2, axis=0)
            if norm_result_3 <= ε:
                print("算法停止，该精度下的最优解是[%f,%f]" % (float(x1_curr), float(x2_curr)))

    return x1_curr, x2_curr


result = fletcher_reeves(1, 1, 0.001)

该代码所运算的结果为：

第1次迭代 当前最优解为[2.000000,0.500000]
第2次迭代 当前最优解为[4.000000,2.000000]
算法停止，该精度下的最优解是[4.000000,2.000000]

简便代码：

'''
FR共轭方向算法   (二元二次函数)
2023.11.8
'''
from sympy import symbols, diff, solve

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

def fletcher_reeves(x1_init, x2_init, ε):
    # 一阶求导
    x1_curr = x1_init
    x2_curr = x2_init
    x = [x1_init, x2_init]
    g = [diff(f, x1), diff(f, x2)]
    p1 = [-g[0].subs({x1: x1_curr, x2: x2_curr}), -g[1].subs({x1: x1_curr, x2: x2_curr})]
    count = 0
    while True:

        g_value = [g[0].subs({x1: x[0], x2: x[1]}), g[1].subs({x1: x[0], x2: x[1]})]
        norm_result = (g_value[0] ** 2 + g_value[1] ** 2) ** 0.5

        if norm_result <= ε:
            print("算法停止，该精度下的最优解是[%f,%f]" % (x[0], x[1]))
            return x
        x_next = [x[0] + λ * p1[0], x[1] + λ * p1[1]]
        f_new = f.subs({x1: x_next[0], x2: x_next[1]})
        grad = diff(f_new, λ)
        λ_value = solve(grad, λ)[0]

        x = [x_next[0].subs(λ, λ_value), x_next[1].subs(λ, λ_value)]
        g_value = [g[0].subs({x1: x[0], x2: x[1]}), g[1].subs({x1: x[0], x2: x[1]})]
        α = (g_value[0] ** 2 + g_value[1] ** 2) / (g[0].subs({x1: x1_curr, x2: x2_curr}) ** 2 + g[1].subs({x1: x1_curr, x2: x2_curr}) ** 2)
        p1 = [-g_value[0] + α * p1[0], -g_value[1] + α * p1[1]]
        count = count + 1
        print("第{}次迭代".format(count))

result = fletcher_reeves(1, 1, 0.001)

该代码所运算的结果为：

第1次迭代
第2次迭代
算法停止，该精度下的最优解是[4.000000,2.000000]

总结：

FR 共轭梯度法的特点

（1)对凸函数全局收敛（下降算法）；

（2）计算公式简单，不用求Hesse矩阵或者逆矩阵，计算量小，存储量小，每步迭代只需存储若干向量，适用于大规模问题；

（3）具有二次收敛性；

（4）Crowder 和 Wolfe 证明，共轭梯度法的收敛速度不坏于最速下降法。如果初始方向不用负梯度方向，则其收敛速度可能像线性收敛速率那样慢。

（5）整体来看，PPR的收敛结果是最好的。

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

背对人潮

关注

13
点赞
踩
51

收藏

觉得还不错? 一键收藏
打赏
6
评论
共轭梯度法（算法分析+python代码解释）

为了克服最速下降法在最优解附近收敛速度变慢和Newton法不具有整体收敛性以及需要计算Hesse矩阵的缺点，一种介于两者之间的方法被提出，这就是共轭方向法。这种方法不仅克服了上述缺点之外，还能保持上述两种方法各自的优点，即收敛性较好，并且可以在有限布内使二次函数达到最优，即具有二次终止性。
复制链接

扫一扫