共轭梯度法详细推导分析

共轭梯度法是一种经典的优化算法。算法求解速度较快,虽然比梯度下降法复杂,但是比二阶方法简单。

一、引入

1. 优化模型建立

假定待优化的问题如下所示:
min ⁡ x f ( x ) = 1 2 x T A x − b T x \min_{x} f(x)=\frac{1}{2}x^TAx - b^Tx xminf(x)=21xTAxbTx
其中 x x x为待优化变量, A A A为半正定矩阵(在线性代数中,正定矩阵为对称矩阵), b b b为已知变量。下标k表示优化步数,负梯度为
r k = − ( A x k − b ) r_k =-( Ax_k -b) rk=(Axkb)
假设最优变量为 x ∗ x^* x,则优化问题可变为求方程 A x ∗ = b Ax^*=b Ax=b的解。梯度 r r r也可以称作每一步的残差。误差定义为 x x x与最优变量的差值
e k = x ∗ − x k e_k = x^* - x_k ek=xxk

2. 算法思想简述

虽然梯度下降法的每一步都是朝着局部最优的方向前进的,但是它在不同的迭代轮数中会选择非常近似的方向,说明这个方向的误差并没通过一次更新方向和步长更新完,在这个方向上还存在误差,因此参数更新的轨迹是锯齿状。共轭梯度法的思想是,选择一个优化方向后,本次选择的步长能够将这个方向的误差更新完,在以后的优化更新过程中不再需要朝这个方向更新了。由于每次将一个方向优化到了极小,后面的优化过程将不再影响之前优化方向上的极小值,所以理论上对N维问题求极小只用对N个方向都求出极小就行了。为了不影响之前优化方向上的更新量,需要每次优化方向共轭正交。假定每一步的优化方向用 p k p_k pk表示,可得共轭正交
p i A p j = 0 i ≠ j p_iAp_j = 0 \qquad i\ne j piApj=0i̸=j
由此可得,每一步优化后,当前的误差和刚才的优化方向共轭正交。
p k A e k + 1 = 0 p_kAe_{k+1}=0 pkAek+1=0
若为N维空间优化问题,则每次优化方向可以组成这个空间中的一组基底。 P = { p 1 , p 2 , … , p N } P=\{p_1,p_2,\dots,p_N\} P={p1,p2,,pN}

二、算法推导

算法只需要解决两个问题:

  1. 优化方向
  2. 优化步长

1.优化方向确定

假定第一次优化方向为初始负梯度方向
p 1 = r 1 = b − A x 1 p_1 = r_1 = b-Ax_1 p1=r1=bAx1
每一次优化方向与之前的优化方向正交,采用Gram-Schmidt方法进行向量正交化,每次优化方向根据当前步的梯度得出
p k = r k − ∑ i &lt; k p i T A r k p i T A p i p i p_k = r_k-\sum_{i&lt;k}\frac{p_i^TAr_k}{p_i^TAp_i}p_i pk=rki<kpiTApipiTArkpi
便于后面证明,令 β i = p i T A r k p i T A p i \beta_i=\frac{p_i^TAr_k}{p_i^TAp_i} βi=piTApipiTArk
上式在后面还会进一步优化,省去求和符号。

2.优化步长的选取

假定第k步的优化步长为 α k \alpha_k αk

方法一:

f ( x k + 1 ) = f ( x k + α k p k ) = g ( α k ) f(x_{k+1})=f(x_k+\alpha_kp_k)=g(\alpha_k) f(xk+1)=f(xk+αkpk)=g(αk),对 α k \alpha_k αk求导令导数为0可得 α k = p k T r k p k T A p k \alpha_k=\frac{p_k^Tr_k}{p_k^TAp_k} αk=pkTApkpkTrk

方法二:
p k T A e k + 1 = p k T A ( x ∗ − x k + 1 ) = p k T A ( x ∗ − x k + x k − x k + 1 ) = p k T A ( e k − α k p k ) = p k T A e k − α k p k T A p k = 0 \begin{aligned} p_k^TAe_{k+1}&amp;=p_k^TA(x^*-x_{k+1})\\ &amp;=p_k^TA(x^*-x_k+x_k-x_{k+1})\\ &amp;=p_k^TA(e_k-\alpha_kp_k)\\ &amp;=p_k^TAe_k-\alpha_kp_k^TAp_k=0 \end{aligned} pkTAek+1=pkTA(xxk+1)=pkTA(xxk+xkxk+1)=pkTA(ekαkpk)=pkTAekαkpkTApk=0
可得
α k = p k T A e k p k T A p k = p k T A ( x ∗ − x k ) p k T A p k = p k T ( A x ∗ − A x k ) p k T A p k = p k T ( b − A x k ) p k T A p k = p k T r k p k T A p k \begin{aligned} \alpha_k&amp;=\frac{p_k^TAe_k}{p_k^TAp_k}\\ &amp;=\frac{p_k^TA(x^*-x_k)}{p_k^TAp_k}\\ &amp;=\frac{p_k^T(Ax^*-Ax_k)}{p_k^TAp_k}\\ &amp;=\frac{p_k^T(b-Ax_k)}{p_k^TAp_k}\\ &amp;=\frac{p_k^Tr_k}{p_k^TAp_k}\\ \end{aligned} αk=pkTApkpkTAek=pkTApkpkTA(xxk)=pkTApkpkT(AxAxk)=pkTApkpkT(bAxk)=pkTApkpkTrk
上式在后文还会进一步化简。

三、三个推论

1.推论一

第k步计算的梯度 r k r_k rk和前k-1步的优化向量 { p i } i = 1 k − 1 \{p_i\}_{i=1}^{k-1} {pi}i=1k1正交。

证明:
i &lt; j i&lt;j i<j
p i T r j = p i T ( A x j − b ) = p i T ( A x j − A x ∗ ) = p i T A e j = p i T A ( e i + 1 − ∑ k = 1 j − 1 β k p k ) = 0 \begin{aligned} p_i^Tr_j &amp;=p_i^T(Ax_j-b) \\ &amp;=p_i^T(Ax_j-Ax^*)\\ &amp;=p_i^TAe_j\\ &amp;=p_i^TA(e_{i+1}-\sum_{k=1}^{j-1}\beta_kp_k)\\ &amp;=0 \end{aligned} piTrj=piT(Axjb)=piT(AxjAx)=piTAej=piTA(ei+1k=1j1βkpk)=0

2.推论二

第k步计算的梯度 r k r_k rk和前k-1步的梯度 { r i } i = 1 k − 1 \{r_i\}_{i=1}^{k-1} {ri}i=1k1正交。

证明:
i &lt; j i&lt;j i<j
r i T r j = ( p i + ∑ k = 1 i − 1 β k p k ) r j = 0 \begin{aligned} r_i^Tr_j=(p_i+\sum_{k=1}^{i-1}\beta_kp_k)r_j=0 \end{aligned} riTrj=(pi+k=1i1βkpk)rj=0

3.推论三

第k步计算的梯度 r k r_k rk和前k-2步的优化向量 { p i } i = 1 k − 2 \{p_i\}_{i=1}^{k-2} {pi}i=1k2共轭正交。

证明:
r j + 1 T r i = ( b − A x j + 1 ) T r i = ( b − A ( x j + α j p j ) ) T r i = ( b − A x j − α j A p j ) T r i = ( r j − α j A p j ) T r i = r j T r i − α j p j T A r i \begin{aligned} r_{j+1}^Tr_i&amp;=(b-Ax_{j+1})^Tr_i\\ &amp;=(b-A(x_j+\alpha_jp_j))^Tr_i\\ &amp;=(b-Ax_j-\alpha_j Ap_j)^Tr_i\\ &amp;=(r_j-\alpha_jAp_j)^Tr_i\\ &amp;=r_j^Tr_i-\alpha_jp_j^TAr_i \end{aligned} rj+1Tri=(bAxj+1)Tri=(bA(xj+αjpj))Tri=(bAxjαjApj)Tri=(rjαjApj)Tri=rjTriαjpjTAri
j = i − 1 j=i-1 j=i1时, p j T A r i ≠ 0 p_j^TAr_i\ne 0 pjTAri̸=0

j + 1 &lt; i j+1&lt;i j+1<i时, p j T A r i = 0 p_j^TAr_i= 0 pjTAri=0

四、最终简化

算法在三中基本推导完毕,但是在工程应用中如果每次进行 p k p_k pk的正交化需要对之前所有的优化向量求解 β \beta β,现简化如下:

1. 优化方向简化

由推论三可得
p k + 1 = r k + 1 − p k T A r k + 1 p k T A k k p k = r k + 1 − ( A p k ) T r k + 1 ( A p k ) T p k p k = r k + 1 − ( r k − r k + 1 α ) T r k + 1 ( r k − r k + 1 α ) T p k p k = r k + 1 − ( r k − r k + 1 α ) T r k + 1 ( r k − r k + 1 α ) T ( r k − β k − 1 p k − 1 ) p k = r t + 1 + r k + 1 T r k + 1 r k T r k p k \begin{aligned} p_{k+1}&amp;=r_{k+1}-\frac{p_k^TAr_{k+1}}{p_k^TAk_k}p_k\\ &amp;=r_{k+1}-\frac{(Ap_k)^Tr_{k+1}}{(Ap_k)^Tp_k}p_k\\ &amp;=r_{k+1}-\frac{(\frac{r_{k}-r_{k+1}}{\alpha})^Tr_{k+1}}{(\frac{r_{k}-r_{k+1}}{\alpha})^Tp_k}p_k\\ &amp;=r_{k+1}-\frac{(\frac{r_{k}-r_{k+1}}{\alpha})^Tr_{k+1}}{(\frac{r_{k}-r_{k+1}}{\alpha})^T(r_k-\beta_{k-1}p_{k-1})}p_k\\ &amp;=r_{t+1}+\frac{r_{k+1}^Tr_{k+1}}{r_k^Tr_k}p_k \end{aligned} pk+1=rk+1pkTAkkpkTArk+1pk=rk+1(Apk)Tpk(Apk)Trk+1pk=rk+1(αrkrk+1)Tpk(αrkrk+1)Trk+1pk=rk+1(αrkrk+1)T(rkβk1pk1)(αrkrk+1)Trk+1pk=rt+1+rkTrkrk+1Trk+1pk

2. 步长简化

第三个等式引用推论一
α k = p k T r k p k t A p k = ( r k − β k − 1 p k − 1 ) T r k p k t A p k = r k T r k p k T A p k T \begin{aligned} \alpha_k &amp;= \frac{p_k^Tr_k}{p_k^tAp_k}\\ &amp;=\frac{(r_k-\beta_{k-1}p_{k-1})^Tr_k}{p_k^tAp_k}\\ &amp;=\frac{r_k^Tr_k}{p_k^TAp_k^T} \end{aligned} αk=pktApkpkTrk=pktApk(rkβk1pk1)Trk=pkTApkTrkTrk

3. 梯度计算简化

r k + 1 = b − A x k + 1 = b − A ( x k + α k p k ) = b − A x k − α k A p k = r k − α k A p k \begin{aligned} r_{k+1}&amp;=b-Ax_{k+1}\\ &amp;=b-A(x_k+\alpha_kp_k)\\ &amp;=b-Ax_k-\alpha_kAp_k\\ &amp;=r_k-\alpha_kAp_k \end{aligned} rk+1=bAxk+1=bA(xk+αkpk)=bAxkαkApk=rkαkApk

最终的推导结束。整理为如下的伪代码

五、伪代码

r 0 = b − A x 0 r_0=b-Ax_0 r0=bAx0
p 0 = r 0 p_0=r_0 p0=r0
k=0
while
  α k = r k T r k p k T A p k \alpha_k=\frac{r_k^Tr_k}{p_k^TAp_k} αk=pkTApkrkTrk
  x k + 1 = x k + α k p k x_{k+1}=x_k+\alpha_kp_k xk+1=xk+αkpk
  r k + 1 = r k − α k A p k r_{k+1}=r_k-\alpha_kAp_k rk+1=rkαkApk
 if r k + 1 r_{k+1} rk+1 < ϵ \epsilon ϵ: break
  β k + 1 = r k + 1 T r k + 1 r k T r k \beta_{k+1}= \frac{r_{k+1}^Tr_{k+1}}{r_k^Tr_k} βk+1=rkTrkrk+1Trk+1
  p k + 1 = r k + 1 + β k p k p_{k+1}=r_{k+1}+\beta_kp_k pk+1=rk+1+βkpk
  k = k + 1 k=k+1 k=k+1
return x k + 1 x_{k+1} xk+1

评论 20 您还未登录,请先 登录 后发表或查看评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
©️2022 CSDN 皮肤主题:编程工作室 设计师:CSDN官方博客 返回首页

打赏作者

贰锤

你的鼓励将是我创作的最大动力

¥2 ¥4 ¥6 ¥10 ¥20
输入1-500的整数
余额支付 (余额:-- )
扫码支付
扫码支付:¥2
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值