用基本共轭方向法（详见博文《最优化方法Python计算：基本共轭方向算法》）计算正定二次型目标函数 f ( x ) = 1 2 x ⊤ H x − x ⊤ b f(\boldsymbol{x})=\frac{1}{2}\boldsymbol{x}^\top\boldsymbol{Hx}-\boldsymbol{x}^\top\boldsymbol{b} f(x)=21x⊤Hx−x⊤b， x ∈ R n \boldsymbol{x}\in\text{ℝ}^n x∈Rn的最优解点 x 0 \boldsymbol{x}_0 x0，效率是很高的：至多迭代 n n n次，并且初始点的选取是任意的。然而，共轭方向法需要事先计算正定阵 H \boldsymbol{H} H的共轭向量组 d 1 , d 2 , ⋯ , d n \boldsymbol{d}_1,\boldsymbol{d}_2,\cdots,\boldsymbol{d}_n d1,d2,⋯,dn。本文探讨一个在搜索过程中动态生成共轭向量 d 1 , d 2 , ⋯ , d k \boldsymbol{d}_1,\boldsymbol{d}_2,\cdots,\boldsymbol{d}_k d1,d2,⋯,dk， k = 1 , 2 , ⋯ k=1,2,\cdots k=1,2,⋯的搜索方法。
定理1 二次型目标函数
f ( x ) = 1 2 x ⊤ H x − x ⊤ b , x ∈ R n f(\boldsymbol{x})=\frac{1}{2}\boldsymbol{x}^\top\boldsymbol{Hx}-\boldsymbol{x}^\top\boldsymbol{b},\boldsymbol{x}\in\text{ℝ}^n f(x)=21x⊤Hx−x⊤b,x∈Rn
其中， H ∈ R n × n \boldsymbol{H}\in\text{ℝ}^{n\times n} H∈Rn×n为正定矩阵。任取 x 1 ∈ R n \boldsymbol{x}_1\in\text{ℝ}^n x1∈Rn，记 ∇ f ( x 1 ) = g 1 \nabla f(\boldsymbol{x}_1)=\boldsymbol{g}_1 ∇f(x1)=g1，若 g 1 ≠ o \boldsymbol{g}_1\not=\boldsymbol{o} g1=o，设 d 1 = − g 1 , α 1 = − g 1 ⊤ d 1 d 1 ⊤ H d 1 \boldsymbol{d}_1=-\boldsymbol{g}_1,\alpha_1=-\frac{\boldsymbol{g}_1^\top\boldsymbol{d}_1}{\boldsymbol{d}_1^\top\boldsymbol{Hd}_1} d1=−g1,α1=−d1⊤Hd1g1⊤d1。对 1 ≤ k ≤ n 1\leq k\leq n 1≤k≤n， x k + 1 = x k + α k d k \boldsymbol{x}_{k+1}=\boldsymbol{x}_k+\alpha_k\boldsymbol{d}_k xk+1=xk+αkdk。假定对每个 k k k， g k + 1 = ∇ f ( x k ) ≠ o \boldsymbol{g}_{k+1}=\nabla f(\boldsymbol{x}_k)\not=\boldsymbol{o} gk+1=∇f(xk)=o， α k = − g k ⊤ d k d k ⊤ H d k \alpha_k=-\frac{\boldsymbol{g}_k^\top\boldsymbol{d}_k}{\boldsymbol{d}_k^\top\boldsymbol{Hd}_k} αk=−dk⊤Hdkgk⊤dk， β k = g k + 1 ⊤ H d k d k ⊤ H d k \beta_k=\frac{\boldsymbol{g}_{k+1}^\top\boldsymbol{Hd}_{k}}{\boldsymbol{d}_k^\top\boldsymbol{Hd}_k} βk=dk⊤Hdkgk+1⊤Hdk， d k + 1 = − g k + β k d k \boldsymbol{d}_{k+1}=-\boldsymbol{g}_k+\beta_k\boldsymbol{d}_k dk+1=−gk+βkdk，则 d 1 , d 2 , ⋯ , d k + 1 \boldsymbol{d}_1,\boldsymbol{d}_2,\cdots,\boldsymbol{d}_{k+1} d1,d2,⋯,dk+1关于 H \boldsymbol{H} H共轭，且则存在 1 ≤ m ≤ n + 1 1\leq m\leq n+1 1≤m≤n+1，使得 x m = x 0 \boldsymbol{x}_{m}=\boldsymbol{x}_0 xm=x0。其中 x 0 \boldsymbol{x}_0 x0为 f ( x ) f(\boldsymbol{x}) f(x)的最优解点。
利用定理1，我们将解正定二次型函数 f ( x ) = 1 2 x ⊤ H x − x ⊤ b f(\boldsymbol{x})=\frac{1}{2}\boldsymbol{x}^\top\boldsymbol{Hx}-\boldsymbol{x}^\top\boldsymbol{b} f(x)=21x⊤Hx−x⊤b最优解的基本共轭算法加以修改，由于正定矩阵 H \boldsymbol{H} H每个共轭向量 d k + 1 = − g k + β k d k \boldsymbol{d}_{k+1}=-\boldsymbol{g}_k+\beta_k\boldsymbol{d}_k dk+1=−gk+βkdk均由梯度 g k = ∇ f ( x k ) \boldsymbol{g}_{k}=\nabla f(\boldsymbol{x}_k) gk=∇f(xk)算得，故称为共轭梯度算法。以下的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表示初始迭代点 x 1 \boldsymbol{x}_1 x1，H，b，c分别表示函数 f ( x ) f(\boldsymbol{x}) f(x)表达式中的矩阵 H \boldsymbol{H} H，向量 b \boldsymbol{b} b和常数 c c c，gtol表示容错误差 ε \varepsilon ε，缺省值为 1 0 − 5 10^{-5} 10−5。
第4~7行进行初始化操作：第4行用x1初始化表示迭代点 x k \boldsymbol{x}_k xk的xk。第5行计算二次型函数 f ( x ) f(\boldsymbol{x}) f(x)的梯度
g 1 = ∇ f ( x 1 ) = H x 1 − b \boldsymbol{g}_1=\nabla f(\boldsymbol{x}_1)=\boldsymbol{Hx}_1-\boldsymbol{b} g1=∇f(x1)=Hx1−b
赋予gk。第6行按定理1计算搜索方向
d 1 = − ∇ f ( x 1 ) = − g 1 \boldsymbol{d}_1=-\nabla f(\boldsymbol{x}_1)=-\boldsymbol{g}_1 d1=−∇f(x1)=−g1
赋予dk。第7行将迭代次数k初始化为1。\par
第8~14行的while循环执行迭代：第9行按定理1计算
α k = − g k ⊤ d k d k H d k \alpha_k=\frac{-\boldsymbol{g}_k^\top\boldsymbol{d}_k}{\boldsymbol{d}_k\boldsymbol{Hd}_k} αk=dkHdk−gk⊤dk
赋予ak。第10行计算迭代点
x k + 1 = x k + α k d k \boldsymbol{x}_{k+1}=\boldsymbol{x}_k+\alpha_k\boldsymbol{d}_k xk+1=xk+αkdk
更新xk。第11行计算
g k + 1 = H x k + 1 − b \boldsymbol{g}_{k+1}=\boldsymbol{Hx}_{k+1}-\boldsymbol{b} gk+1=Hxk+1−b
更新gk。第12行按定理1计算组合系数
β k = g k + 1 ⊤ H d k d k ⊤ H d k \beta_k=\frac{\boldsymbol{g}_{k+1}^\top\boldsymbol{Hd}_k}{\boldsymbol{d}_k^\top\boldsymbol{Hd}_k} βk=dk⊤Hdkgk+1⊤Hdk
赋予bk。第13行按定理1计算共额方向
d k + 1 = − g k + 1 + β k d k \boldsymbol{d}_{k+1}=-\boldsymbol{g}_{k+1}+\beta_k\boldsymbol{d}_k dk+1=−gk+1+βkdk
更新dk。第14行将迭代次数自增1。循环往复，直至条件
∥ g k ∥ < ε \lVert\boldsymbol{g}_k\rVert<\varepsilon ∥gk∥<ε
成立为止。
第15~17行用 f ( x k ) = 1 2 x k ⊤ H x k − b ⊤ x k + c f(\boldsymbol{x}_k)=\frac{1}{2}\boldsymbol{x}_k^\top\boldsymbol{Hx}_k-\boldsymbol{b}^\top\boldsymbol{x}_k+c f(xk)=21xk⊤Hxk−b⊤xk+c， x k \boldsymbol{x}_k xk和 k k k构造OptimizeResult对象（第2行导入）并返回。
例1 利用共轭梯度法计算正定二次型函数 f ( x ) = 5 2 x 1 2 + x 2 2 − 3 x 1 x 2 − x 2 − 7 f(\boldsymbol{x})=\frac{5}{2}x_1^2+x_2^2-3x_1x_2-x_2-7 f(x)=25x12+x22−3x1x2−x2−7的最优解点 x 0 ∈ R 2 \boldsymbol{x}_0\in\text{ℝ}^2 x0∈R2，给定初始点 x 1 = ( 0 0 ) \boldsymbol{x}_1=\begin{pmatrix}0\\0\end{pmatrix} x1=(00)。
解：目标函数的矩阵形式为
f ( x ) = 1 2 x ⊤ H x − x ⊤ b + c f(\boldsymbol{x})=\frac{1}{2}\boldsymbol{x}^\top\boldsymbol{Hx}-\boldsymbol{x}^\top\boldsymbol{b}+c f(x)=21x⊤Hx−x⊤b+c
其中， x = ( x 1 x 2 ) ∈ R 2 \boldsymbol{x}=\begin{pmatrix}x_1\\x_2\end{pmatrix}\in\text{ℝ}^2 x=(x1x2)∈R2， H = ( 5 − 3 − 3 2 ) \boldsymbol{H}=\begin{pmatrix}5&-3\\-3&2\end{pmatrix} H=(5−3−32)， b = ( 0 1 ) \boldsymbol{b}=\begin{pmatrix}0\\1\end{pmatrix} b=(01)， c = 7 c=7 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.])
这意味着在 ε = 1 0 − 5 \varepsilon=10^{-5} ε=10−5的容错误差下，共轭梯度算法迭代3次算得正定二次型函数 f ( x ) = 5 2 x 1 2 + x 2 2 − 3 x 1 x 2 − x 2 − 7 f(\boldsymbol{x})=\frac{5}{2}x_1^2+x_2^2-3x_1x_2-x_2-7 f(x)=25x12+x22−3x1x2−x2−7的最优解 x 0 = ( 3 5 ) \boldsymbol{x}_0=\begin{pmatrix}3\\5\end{pmatrix} x0=(35)。