最优化方法Python计算：正定二次型基本共轭方向算法

设目标函数
$$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{b}\in {\text{R}}^n$为一常向量。即$f(\mathbf{x})$为一二次型函数，且$\nabla f(\mathbf{x})=\mathbf{Hx}-\mathbf{b}$，${\nabla }^{2}f(\mathbf{x})=\mathbf{H}$。由于$\mathbf{H}={\nabla }^{2}f(\mathbf{x})$正定，故$f(\mathbf{x})$有唯一的最小值点${\mathbf{x}}_0$（满足$\mathbf{o}=\nabla f({\mathbf{x}}_0)={\mathbf{Hx}}_0-\mathbf{b}$或${\mathbf{Hx}}_0=\mathbf{b}$）。
定理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{d}}_1,{\mathbf{d}}_2,\cdots ,{\mathbf{d}}_n$为$f(\mathbf{x})$的共轭方向向量（详见博文《最优化方法Python计算：构造正定矩阵的共轭向量》），任取${\mathbf{x}}_1\in {\text{R}}^n$，对$1\le k\le n$令${\alpha }_k=-\frac{\nabla f({\mathbf{x}}_k{)}^{\top }{\mathbf{d}}_k}{{\mathbf{d}}_k^{\top }{\mathbf{Hd}}_k}$，构造迭代式
$${\mathbf{x}}_{k+1}={\mathbf{x}}_k+{\alpha }_k{\mathbf{d}}_k$$
则存在$1\le m\le n+1$，使得${\mathbf{x}}_m={\mathbf{x}}_{m+1}=\cdots ={\mathbf{x}}_{n+1}={\mathbf{x}}_0$。其中${\mathbf{x}}_0$为$f(\mathbf{x})$的最优解点。
利用定理1，可建立计算二次型目标函数最优解的基本共轭方向算法。下列Python函数实现该算法。

import numpy as np                          #导入numpy
from scipy.optimize import OptimizeResult   #导入OptimizeResult
def conj(x1,H,b,c,gtol=1e-5):               #基本共轭方向算法实现函数
    d=conjugation(H)                        #构造H的共轭向量组
    xk=x1                                   #初始迭代点
    gk=(np.matmul(H,xk)-b)                  #计算梯度
    k=1
    while np.linalg.norm(gk)>=gtol:         #只要梯度不为0
        ak=-np.matmul(gk,d[:,k-1])/np.matmul(np.matmul(d[:,k-1],H),d[:,k-1])
        xk+=ak*d[:,k-1]
        gk=(np.matmul(H,xk)-b)
        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~15行定义实现基本共轭方向法的conj函数，参数x1表示初始迭代点${\mathbf{x}}_1$，H、b和c分别表示二次型函数
$$f(\mathbf{x})=\frac{1}{2}{\mathbf{x}}^{\top }\mathbf{Hx}-{\mathbf{x}}^{\top }\mathbf{b}+c$$
中的矩阵$\mathbf{H}$、向量$\mathbf{b}$和常数$c$，gtol表示容错误差$\epsilon$，缺省值为$10^{-5}$。
第4~7行执行初始化操作：第4行调用博文《最优化方法Python计算：构造正定矩阵的共轭向量》中定义的构造正定矩阵共轭向量组的函数conjugation，构造$\mathbf{H}$的共轭向量组
$${\mathbf{d}}_1,{\mathbf{d}}_2,\cdots ,{\mathbf{d}}_n$$
存储于d。第5行设置迭代点xk初始化为x1。第6行计算目标函数的当前梯度
$${\mathbf{g}}_1={\mathbf{Hx}}_1-\mathbf{b}$$
赋予gk。第7行将迭代次数k初始化为1。\par
第8~12行的while循环作迭代操作：第9行计算
$${\alpha }_k=-\frac{{\mathbf{g}}_k^{\top }{\mathbf{d}}_k}{{\mathbf{d}}_k^{\top }{\mathbf{Hd}}_k}$$
赋予ak。第10行用
$${\mathbf{x}}_k+{\alpha }_k{\mathbf{d}}_k$$
更新xk。第11行用xk更新梯度gk。第12行迭代次数k自增1。循环往复，直至条件
$$\Vert \mathbf{g}({\mathbf{x}}_k)\Vert <10^{-5}$$
（$\mathbf{g}({\mathbf{x}}_k)$近似为$\mathbf{o}$）成立。需要注意的是于Python数组的下标是从0开始编码的，故在第9、10行中k作为数组下标时表示为k-1。
第13~15行用$f({\mathbf{x}}_k)$，${\mathbf{x}}_k$及$k$构造OptimizeResult（第2行导入）对象，并返回。
例1 用函数conj计算正定二次型函数$f(\mathbf{x})=\frac{3}{2}{x}_1^2+2x_2^2+\frac{3}{2}{x}_3^2+x_1{x}_3+2x_2{x}_3-3x_1-x_3$的最优解。给定初始点${\mathbf{x}}_1=\left(\begin{array}{c}0\\ 0\\ 0\end{array}\right)$，容错误差$\epsilon =10^{-5}$。
解：$f(\mathbf{x})$的矩阵形式$f(\mathbf{x})=\frac{1}{2}{\mathbf{x}}^{\top }\mathbf{Hx}-{\mathbf{x}}^{\top }\mathbf{b}+c$中，$\mathbf{H}=\left(\begin{array}{ccc}3& 0& 1\\ 0& 4& 2\\ 1& 2& 3\end{array}\right)$，$\mathbf{b}=\left(\begin{array}{c}3\\ 0\\ 1\end{array}\right)$，$c=0$。下列代码完成计算。

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

借助代码内部的注释信息，不难理解本程序。需要注意的是，由于给定的容错误差为$\epsilon =10^{-5}$，故调用conj时使用参数gtol的缺省值。运行程序，输出

 fun: -1.5
 nit: 2
   x: array([1., 0., 0.])

这意味着在$\epsilon =10^{-5}$的容错误差下，迭代2次即算得$f(\mathbf{x})$的最优解${\mathbf{x}}_0=\left(\begin{array}{c}1\\ 0\\ 0\end{array}\right)$。