最速下降法
定理 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,