最优化方法Python计算：正定二次型共轭梯度算法

用基本共轭方向法（详见博文《最优化方法Python计算：基本共轭方向算法》）计算正定二次型目标函数$f(\mathbf{x})=\frac{1}{2}{\mathbf{x}}^{\top }\mathbf{Hx}-{\mathbf{x}}^{\top }\mathbf{b}$，$\mathbf{x}\in {\text{R}}^n$的最优解点${\mathbf{x}}_0$，效率是很高的：至多迭代$n$次，并且初始点的选取是任意的。然而，共轭方向法需要事先计算正定阵$\mathbf{H}$的共轭向量组${\mathbf{d}}_1,{\mathbf{d}}_2,\cdots ,{\mathbf{d}}_n$。本文探讨一个在搜索过程中动态生成共轭向量${\mathbf{d}}_1,{\mathbf{d}}_2,\cdots ,{\mathbf{d}}_k$，$k=1,2,\cdots$的搜索方法。
定理1 二次型目标函数
$$f(\mathbf{x})=\frac{1}{2}{\mathbf{x}}^{\top }\mathbf{Hx}-{\mathbf{x}}^{\top }\mathbf{b},\mathbf{x}\in {\text{R}}^n$$
其中，$\mathbf{H}\in {\text{R}}^{n\times n}$为正定矩阵。任取${\mathbf{x}}_1\in {\text{R}}^n$，记$\nabla f({\mathbf{x}}_1)={\mathbf{g}}_1$，若${\mathbf{g}}_1\ne \mathbf{o}$，设${\mathbf{d}}_1=-{\mathbf{g}}_1,{\alpha }_1=-\frac{{\mathbf{g}}_1^{\top }{\mathbf{d}}_1}{{\mathbf{d}}_1^{\top }{\mathbf{Hd}}_1}$。对$1\le k\le n$，${\mathbf{x}}_{k+1}={\mathbf{x}}_k+{\alpha }_k{\mathbf{d}}_k$。假定对每个$k$，${\mathbf{g}}_{k+1}=\nabla f({\mathbf{x}}_k)\ne \mathbf{o}$，${\alpha }_k=-\frac{{\mathbf{g}}_k^{\top }{\mathbf{d}}_k}{{\mathbf{d}}_k^{\top }{\mathbf{Hd}}_k}$，${\beta }_k=\frac{{\mathbf{g}}_{k+1}^{\top }{\mathbf{Hd}}_k}{{\mathbf{d}}_k^{\top }{\mathbf{Hd}}_k}$，${\mathbf{d}}_{k+1}=-{\mathbf{g}}_k+{\beta }_k{\mathbf{d}}_k$，则${\mathbf{d}}_1,{\mathbf{d}}_2,\cdots ,{\mathbf{d}}_{k+1}$关于$\mathbf{H}$共轭，且则存在$1\le m\le n+1$，使得${\mathbf{x}}_m={\mathbf{x}}_0$。其中${\mathbf{x}}_0$为$f(\mathbf{x})$的最优解点。
利用定理1，我们将解正定二次型函数$f(\mathbf{x})=\frac{1}{2}{\mathbf{x}}^{\top }\mathbf{Hx}-{\mathbf{x}}^{\top }\mathbf{b}$最优解的基本共轭算法加以修改，由于正定矩阵$\mathbf{H}$每个共轭向量${\mathbf{d}}_{k+1}=-{\mathbf{g}}_k+{\beta }_k{\mathbf{d}}_k$均由梯度${\mathbf{g}}_k=\nabla f({\mathbf{x}}_k)$算得，故称为共轭梯度算法。以下的Python函数实现共轭梯度算法。

import numpy as np                                                      #导入numpy
from scipy.optimize import OptimizeResult                               #导入OptimizeResult
def conjG(x1,H,b,c,gtol=1e-5):                                          #实现共轭梯度算法的函数
    xk=x1                                                               #设置初始迭代点
    gk=(np.matmul(H,xk)-b)                                              #计算当前梯度
    dk=-gk                                                              #搜索方向
    k=1
    while np.linalg.norm(gk)>=gtol:                                     #只要梯度不为0
        ak=-np.matmul(gk,dk)/np.matmul(np.matmul(dk,H),dk)              #搜索步长
        xk+=ak*dk                                                       #更新迭代点
        gk=(np.matmul(H,xk)-b)                                          #计算当前梯度
        bk=np.matmul(np.matmul(gk,H),dk)/np.matmul(np.matmul(dk,H),dk)  #计算beta
        dk=-gk+bk*dk                                                    #计算搜索方向
        k+=1
    bestx=xk
    besty=(np.matmul(np.matmul(xk,H),xk)/2-np.matmul(b,xk))+c
    return OptimizeResult(fun=besty, x=bestx, nit=k)

程序的第3~17行定义的conjG函数实现共轭梯度算法，参数x1表示初始迭代点${\mathbf{x}}_1$，H，b，c分别表示函数$f(\mathbf{x})$表达式中的矩阵$\mathbf{H}$，向量$\mathbf{b}$和常数$c$，gtol表示容错误差$\epsilon$，缺省值为$10^{-5}$。
第4~7行进行初始化操作：第4行用x1初始化表示迭代点${\mathbf{x}}_k$的xk。第5行计算二次型函数$f(\mathbf{x})$的梯度
$${\mathbf{g}}_1=\nabla f({\mathbf{x}}_1)={\mathbf{Hx}}_1-\mathbf{b}$$
赋予gk。第6行按定理1计算搜索方向
$${\mathbf{d}}_1=-\nabla f({\mathbf{x}}_1)=-{\mathbf{g}}_1$$
赋予dk。第7行将迭代次数k初始化为1。\par
第8~14行的while循环执行迭代：第9行按定理1计算
$${\alpha }_k=\frac{-{\mathbf{g}}_k^{\top }{\mathbf{d}}_k}{{\mathbf{d}}_k{\mathbf{Hd}}_k}$$
赋予ak。第10行计算迭代点
$${\mathbf{x}}_{k+1}={\mathbf{x}}_k+{\alpha }_k{\mathbf{d}}_k$$
更新xk。第11行计算
$${\mathbf{g}}_{k+1}={\mathbf{Hx}}_{k+1}-\mathbf{b}$$
更新gk。第12行按定理1计算组合系数
$${\beta }_k=\frac{{\mathbf{g}}_{k+1}^{\top }{\mathbf{Hd}}_k}{{\mathbf{d}}_k^{\top }{\mathbf{Hd}}_k}$$
赋予bk。第13行按定理1计算共额方向
$${\mathbf{d}}_{k+1}=-{\mathbf{g}}_{k+1}+{\beta }_k{\mathbf{d}}_k$$
更新dk。第14行将迭代次数自增1。循环往复，直至条件
$$\Vert {\mathbf{g}}_k\Vert <\epsilon$$
成立为止。
第15~17行用$f({\mathbf{x}}_k)=\frac{1}{2}{\mathbf{x}}_k^{\top }{\mathbf{Hx}}_k-{\mathbf{b}}^{\top }{\mathbf{x}}_k+c$，${\mathbf{x}}_k$和$k$构造OptimizeResult对象（第2行导入）并返回。
例1 利用共轭梯度法计算正定二次型函数$f(\mathbf{x})=\frac{5}{2}{x}_1^2+x_2^2-3x_1{x}_2-x_2-7$的最优解点${\mathbf{x}}_0\in {\text{R}}^2$，给定初始点${\mathbf{x}}_1=\left(\begin{array}{c}0\\ 0\end{array}\right)$。
解：目标函数的矩阵形式为
$$f(\mathbf{x})=\frac{1}{2}{\mathbf{x}}^{\top }\mathbf{Hx}-{\mathbf{x}}^{\top }\mathbf{b}+c$$
其中，$\mathbf{x}=\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)\in {\text{R}}^2$，$\mathbf{H}=\left(\begin{array}{cc}5& -3\\ -3& 2\end{array}\right)$，$\mathbf{b}=\left(\begin{array}{c}0\\ 1\end{array}\right)$，$c=7$。下列代码完成计算。

import numpy as np                          #导入numpy
H=np.array([[5, -3],                        #设置Hesse阵
            [-3, 2]],dtype='float')
b=np.array([0,1],dtype='float')             #设置向量
c=7                                         #设置常数
x=np.array([0,0],dtype='float')             #设置初始迭代点
print(conjG(x,H,b,c))                       #计算并输出最优解

利用代码内的注释信息容易理解本程序代码。运行程序，输出

 fun: 4.5
 nit: 3
   x: array([3., 5.])

这意味着在$\epsilon =10^{-5}$的容错误差下，共轭梯度算法迭代3次算得正定二次型函数$f(\mathbf{x})=\frac{5}{2}{x}_1^2+x_2^2-3x_1{x}_2-x_2-7$的最优解${\mathbf{x}}_0=\left(\begin{array}{c}3\\ 5\end{array}\right)$。
写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！