最速下降法
定理 1
考虑线性方程组：
Ax=b,A∈Rn×n,A对称正定
和二次泛函
φ(x)=xTAx−2bTx
有定理 1：设 A 对称正定，求方程组 Ax=b 的解等价与求二次泛函 φ(x) 的极小值点。
即求 x∗∈Rn 使得 φ(x∗)=x∈Rminφ(x)
证明：
必要性
设 x∗ 是 φ(x∗) 在Rn 上的极小值点，则 grad(φ(x∗))=∇φ(x)=0
grad(φ(x∗))=⎝⎜⎜⎜⎜⎜⎛∂x1∂φ∂x2∂φ⋮∂xn∂φ⎠⎟⎟⎟⎟⎟⎞
下面开始求 φ(x) :
φ(x)=xTAx−2bTx=[x1,x2,⋯,xn]⎣⎢⎢⎢⎢⎡a11a21⋮an1a12a22⋮an2⋯⋯⋯a1na2n⋮ann⎦⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎡x1x2⋮xn⎦⎥⎥⎥⎥⎤−2[b1,b2,⋯,bn]⎣⎢⎢⎢⎢⎡x1x2⋮xn⎦⎥⎥⎥⎥⎤=[i=1∑nai1xi,i=2∑nai2xi,⋯,i=n∑nainxi]⎣⎢⎢⎢⎢⎡x1x2⋮xn⎦⎥⎥⎥⎥⎤−2i=1∑nbixi=i=1∑nai1xix1+i=1∑nai2xix2+⋯+i=1∑nainxixn−2i=1∑nbixi=j=1∑ni=1∑naijxixj−2i=1∑nbixi
求偏导：
∂xk∂φ(x)=ak1x1+ak2x2+ak,k−1xk−1+xk∂φ(xk)+ak,k+1xk+1+⋯+ak,nxn−2bk
上式中，把 xk∂φ(xk) 单独拎出来计算：
xk∂φ(xk)=xk∂(i=1∑naikxixk)=i=1,i=k∑naikxi+2akkxk
又由于 A 矩阵是对称的，即 A=AT,aki=aik
整理一下，得到：
∂xK∂φ(x)=2i=1∑nakixi−2bk
那么，
grad(φ(x∗))=⎝⎜⎜⎜⎜⎜⎛∂x1∂φ∂x2∂φ⋮∂xn∂φ⎠⎟⎟⎟⎟⎟⎞=2⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛i=1∑na1ixii=1∑na2ixi⋮i=1∑nanixi⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞−2⎝⎜⎜⎜⎜⎛b1b2⋮bn⎠⎟⎟⎟⎟⎞=2(Ax−b)=−2r,(r=b−Ax)
所以，如果 Ax∗=b ，那么，x∗ 就是 Ax=b 的解
充分性
若 Ax∗=b ，则 ∀y∈Rn ，有
φ(x∗+y)=(x∗+y)TA(x∗+y)−2bT(x∗+y)=x∗TAx∗+x∗TAy+yTAx∗+yTAy−2bTx∗−2bTy=x∗TAx∗+bTy+yTb+yTAy−2bTx∗−2bTy(AT=A,Ax∗=b，所以x∗TAy=x∗TATy=Ax∗T=bTy)=x∗TAx∗+yTAy−2bTx∗(bTy和yTb是相等的，都是数字，所以可以跟后面的−2bTy消掉)=φ(x∗)一定为非负数+yTAy≥φ(x∗)
∴x∗ 使得 φ(x) 达到最小
证毕
计算思路
通过前面的定理，那么用最速下降法求方程组就相当于求二次泛函的极小值点，求这个极小值点的方式类似于盲人下山
（1）对给定的初始向量 x0 ，确定一个下山方向 P0 ，沿着直线 x=x0+αP0 寻找点 x1=x0+α0P0 使得 φ(x0+α0P0)≤φ(x0+αP0),α∈R
（2）类似地，找点 x2=x1+α1P1，使得 φ(x2)=φ(x1+α1P1)≤φ(x1+αP1),α∈R
…
（k）找点 xk=xk−1+αk−1Pk−1 ，使得 φ(xk)=φ(xk−1+αk−1Pk−1)≤φ(xk−1+αPk−1),α∈R
…
直到 ∣∣rk∣∣=∣∣b−Axk∣∣≤tol ，tol 是能够容许的误差，通常为一个极小的值
在计算过程中，需要确定步长 αk 和搜索方向 Pk
首先确定搜索方向，应该是 φ(x) 减少速度最快的方向，也就是负梯度方向：Pk=rk=b−Axk
然后确定步长 αk ，使得 f(αk)=αminf(α)=αminφ(xk+αPk)=φ(xk+αkPk)
所以 αk 应为满足 f′(α)=0 的解
下面确定 αk,Pk
f(α)=φ(xk+αPk)=(xk+αPk)TA(xk+αPk)−2bT(xk+αPk)=xkTAxk这两项其实是一样的，可以用转置的定义推导得到+αxkTAPk+αPkTAxk+α2PkTAPk−2bTxk−2αbTPk=xkTAxk+2αxkTAPk+α2PkTAPk−2bTxk−2αbTPk=xkTAxk−2αxkTAPk+α2PkTAPk−2bTxk(rk=b−Axk,2αxkTAPk−2αbTPk=2α((Axk)T−bT)Pk=−2αrkTpk)=α2PkTAPk−2αrkTpk+φ(xk)
解方程得又由于∴:f′(α)=2αPkTAPk−2rkTPk=0:αk=PkTAPkrkTPk:Pk=rk=b−Axkαk=PkTAPkrkTrk
验证不等式 φ(xk+1)=φ(xk+αkPk)<φ(xk) 是否满足
φ(xk+1−φ(xk))=φ(xk+αkPk)−φ(xk)=(xk+αkPk)TA(xk+αkPk)−2bT(xk+αkPk)−(xkTAxk−2bTxk)=xkTAxk+αkxkTAPk+αkPkTAxk+αk2PkTAPk−2bTxk−2αkbTPk−xkTAxk+2bTxk=αkxkTAPk+αkPkTAxk+αk2PkTAPk−2αkbTPk=2αkxkTAPk+αk2PkTAPk−2αkbTPk=αk2PkTAPk−2αkrkTPk=(PkTAPk)2(rkTrk)2PkTAPk−2PkTAPkrkTrkrkTPk=PkTAPk(rkTrk)2−2PkTAPk(rkTrk)2=−PkTAPk(rkTrk)2<0(只要rkTPk=0)
若 rk=0,则:φ(xk+1=φxk)，即解为xk
一些向量的性质
PkTPk+1=0
其中，Pk+1=b−Axk+1=b−A(xk+αkPk)=b−Axk−αkAPk=rk−αkAPk
那么:PkTPk+1=PkTrk−PkTAPkrkTArkrkTrk=rkTrk−rkTrk=0
这一性质的几何含义：如果在二维上，表现为两次下降的方向是垂直的（正交的）
具体计算步骤
步0：任选初始点 x0 ，置 k=0
步1：计算 rk=b−Axk （下降方向为负梯度方向）
步2：若 rk=0 ，则停止
步3：计算 αk=rkTArkrkTrk （步长因子）
步4：计算 xk+1=xk+αkPk （新的迭代点）
步5：k=k+1 ，转步1
补充
定理2：设 A∈Rn×n 对称正定，A 的特征值为：
λn≥λn−1≥⋯≥λ1>0
设 {xk} 为由求解对称正定方程组 Ax=b 的最速下降法得到的向量序列，则有：
∣∣xk−x∗∣∣A≤(λn+λ1λn−λ1)k∣∣x0−x∗∣∣A
其中：Ax∗=b
定义：设 A∈Rn×n 对称正定，则称 ∣∣x∣∣A=(xTAx)21 为 Rn 上的 A−范数
注：可以用范数定义证明 ∣∣∗∣∣A 是范数
从定理的结论上看，当 λ1λn≫1 时，迭代序列收敛速度很慢
缺点：从局部上看，最速下降法的每一步迭代都是最优的，但是从整体上看却未必，为获得全局最优的迭代方法，需要采用共轭梯度法
共轭梯度法
推导
（1）∀x0∈Rn,r0=b−Ax0 ，仍取 P0=r0 。x0 沿 P0 方向作直线搜索
使得 φ(x0+α0P0)=αminφ(x0+αP0),α∈R 证明方法与最速下降法类似
解之，得 α0=r0TAr0r0Tr0>0 【 通过求f′(α)=0 中的 α 得到 】
∴x1=x0+α0P0,r1=b−Ax1
注：P0Tr1=P0T(b−Ax0−α0AP0)=P0Tr0−α0P0TAP0r0=P00
r0Tr1=P0Tr1=0(∗)
（2）从第二步起，下山方向就不再选取负梯度方向，而是在二维平面 π2={x=x1+ξr1+ηP0∣ξ,η∈R} 内找出使函数下降最快的方向作为新的下山方向 P1 和新的步长 α1 。
φ(ξ,η)=φ(x1+ξr1+ηP0)=(x1+ξr1+ηP0)TA(x1+ξr1+ηP0)−2bT(x1+ξr1+ηP0)=x1TAξr1+ξr1TAx1+ξ2r1TAr1+ξr1TAηP0+ηP0TAξr1−2bTξr1
对其求偏导，得到
ξφ(ξ,η)∂ξ∂φ(ξ,η)∂ξ∂φ(ξ,η)=2(ξr1TAr1+ηr1TAP0−r1Tr1)=2(ξr1TAP0+ηP0TAP0)
令（19）（20）为 0 ，即得 φ 在 π2 内的唯一极小值解
x=x1+ξ0r1+η0P0
其中：ξ0,η0 满足（19）=（20）= 0
即求方程组：{ξr1TAr1+ηr1TAP0=r1Tr1ξr1TAP0+ηP0TAP0 【 如果这里直接求解 ξ0,η0 比较困难，可以先求出 β0=ξ0η0=−P0TAP0r1TAP0 ，得出方向，再求方向的系数 】
（21）中， 若 r1=0则ξ0=0 ，故可取 x2=x1+α1P1 【 P1=r1+β0P0 作为新的方向】
因 α1 也满足 φ(x1+α1P1)=αminφ(x1+αP1),α∈R
故，得到： α1=P1TAP1r1TP1
由 P1=r1+β0P0 ，得 α1=P1TAP1r1T+β0P0=P1TAP1r1Tr1 【 由 (∗)，r1 和 P0 垂直，即 r1TP0=(P0Tr1)T=0 】
综上：x2=x1+α1P1,P1=r1+β0P0,r2=b−Ax2
一些向量的性质
P1TAP0=(r1+β0P0)TAP0=r1TAP0+β0P0TAP0=r1TAP0−P0TAP0r1TAP0P0TAP0=0
P0Tr2=P0T(r1b−Ax1−α1AP1)=P0Tr1−α1P0TAP1=0
P1Tr2=P1T(b−Ax1−α1AP1)=P1Tr1−α1P1TAP1=P1Tr1−P1TAP1r1Tr1P1TAP1=P1Tr1−r1Tr1=(r1+β0P0)Tr1−r1Tr1=β0P0Tr1=β0⋅0=0
r0Tr2=P0Tr2=0
r1Tr2=r1T(b−Ax1−α1AP1)=r1T(r1−α1AP1)=r1Tr1−r1Tα1AP1=r1Tr1−α1r1TAP1=r1Tr1−P1TAP1r1TP1r1TAP1=r1Tr1−(r1+β0P0)TAP1r1TP1r1TAP1=r1Tr1−r1TAP1+β0P0TAP1r1TP1r1TAP1=r1Tr1−r1TAP1r1Tr1r1TAP1=0
（3）若 r2=0 ，则重复第二步，总结为：
xk+1Pkβk−1αkrk+1=xk+αkPk=rk+βk−1Pk−1=−Pk−1TAPk−1rkTAPk−1=rk−1Trk−1rkTrk−1=PkTAPkrkTrk=b−Axk+1=rk−αkAPk
如此，便得到共轭梯度法的计算步骤
计算步骤
推导出来的计算步骤如下：
x0=初值
r0=b−Ax0;k=0
whilerk=0【实际上，应该用范数是否小于设定值来判断
k=k+1
ifk=1P0=r0
else
βk−2=rk−2Trk−2rk−1Trk−1
Pk−1=rk−1+βk−2Pk−2
end
αk−1=Pk−1TAPk−1rk−1Trk−1
xk=xk−1+αk−1Pk−1
rk=rk−1−αk−1APk−1
end
x=xk
但是在编写程序过程中，应当按照如下方式编写：
x0=初值
r0=b−Ax0;k=0
while∣∣rk∣∣≥tol(tol为一设定的很小的值)
k=k+1
ifk=1P0=r0
else
βk−2=rk−2Trk−2rk−1Trk−1
Pk−1=rk−1+βk−2Pk−2
end
temp=APk−1
αk−1=Pk−1T⋅temprk−1Trk−1
xk=xk−1+αk−1Pk−1
rk=rk−1−αk−1⋅temp
end
x=xk
补充
定理3：由共轭梯度法得到的向量组 {ri} 和 {Pi} 具有下面的性质：
（1）PITrj=0,0≤i<j≤k
（2）riTrj=0,i=j,0≤i,j≤k
（3）PiTAPj=0,i=j,0≤i,j≤k
（4）span{r0,r1,⋯,rk}=span{P0,P1,⋯,Pk}=K(A,r0,k+1)
K(A,r0,k+1)=span{r0,Ar0,⋯,Akr0}通常称为krylov子空间
A−共轭 ：设矩阵 A∈Rn×n 对称正定，若非零向量 P0,P1,⋯,Pk−1 满足 PiTAPj=0,i=j,i,j=0,⋯,k−1 ，则称向量 P0,P1,⋯,Pk−1 是 A−共轭 的，或称 P0,P1,⋯,Pk−1 是 A−共轭 方向
评论区