共轭梯度法详细推导分析
共轭梯度法是一种经典的优化算法。算法求解速度较快,虽然比梯度下降法复杂,但是比二阶方法简单。
一、引入
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)=21xTAx−bTx
其中 x x x为待优化变量, A A A为半正定矩阵(在线性代数中,正定矩阵为对称矩阵), b b b为已知变量。下标k表示优化步数,负梯度为
r k = − ( A x k − b ) r_k =-( Ax_k -b) rk=−(Axk−b)
假设最优变量为 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=x∗−xk
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.优化方向确定
假定第一次优化方向为初始负梯度方向
p 1 = r 1 = b − A x 1 p_1 = r_1 = b-Ax_1 p1=r1=b−Ax1
每一次优化方向与之前的优化方向正交,采用Gram-Schmidt方法进行向量正交化,每次优化方向根据当前步的梯度得出
p k = r k − ∑ i < k p i T A r k p i T A p i p i p_k = r_k-\sum_{i<k}\frac{p_i^TAr_k}{p_i^TAp_i}p_i pk=rk−i<k∑piTApipiTArkpi
便于后面证明,令 β 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}&=p_k^TA(x^*-x_{k+1})\\ &=p_k^TA(x^*-x_k+x_k-x_{k+1})\\ &=p_k^TA(e_k-\alpha_kp_k)\\ &=p_k^TAe_k-\alpha_kp_k^TAp_k=0 \end{aligned} pkTAek+1=pkTA(x∗−xk+1)=pkTA(x∗−xk+xk−xk+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&=\frac{p_k^TAe_k}{p_k^TAp_k}\\ &=\frac{p_k^TA(x^*-x_k)}{p_k^TAp_k}\\ &=\frac{p_k^T(Ax^*-Ax_k)}{p_k^TAp_k}\\ &=\frac{p_k^T(b-Ax_k)}{p_k^TAp_k}\\ &=\frac{p_k^Tr_k}{p_k^TAp_k}\\ \end{aligned} αk=pkTApkpkTAek=pkTApkpkTA(x∗−xk)=pkTApkpkT(Ax∗−Axk)=pkTApkpkT(b−Axk)=pkTApkpkTrk
上式在后文还会进一步化简。
三、三个推论
1.推论一
第k步计算的梯度 r k r_k rk和前k-1步的优化向量 { p i } i = 1 k − 1 \{p_i\}_{i=1}^{k-1} {pi}i=1k−1正交。
证明:
当 i < j i<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 &=p_i^T(Ax_j-b) \\ &=p_i^T(Ax_j-Ax^*)\\ &=p_i^TAe_j\\ &=p_i^TA(e_{i+1}-\sum_{k=1}^{j-1}\beta_kp_k)\\ &=0 \end{aligned} piTrj=piT(Axj−b)=piT(Axj−Ax∗)=piTAej=piTA(ei+1−k=1∑j−1β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=1k−1正交。
证明:
当 i < j i<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=1∑i−1β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=1k−2共轭正交。
证明:
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&=(b-Ax_{j+1})^Tr_i\\ &=(b-A(x_j+\alpha_jp_j))^Tr_i\\ &=(b-Ax_j-\alpha_j Ap_j)^Tr_i\\ &=(r_j-\alpha_jAp_j)^Tr_i\\ &=r_j^Tr_i-\alpha_jp_j^TAr_i \end{aligned} rj+1Tri=(b−Axj+1)Tri=(b−A(xj+αjpj))Tri=(b−Axj−αjApj)Tri=(rj−αjApj)Tri=rjTri−αjpjTAri
当 j = i − 1 j=i-1 j=i−1时, p j T A r i ≠ 0 p_j^TAr_i\ne 0 pjTAri̸=0。
当 j + 1 < i j+1<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}&=r_{k+1}-\frac{p_k^TAr_{k+1}}{p_k^TAk_k}p_k\\ &=r_{k+1}-\frac{(Ap_k)^Tr_{k+1}}{(Ap_k)^Tp_k}p_k\\ &=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\\ &=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\\ &=r_{t+1}+\frac{r_{k+1}^Tr_{k+1}}{r_k^Tr_k}p_k \end{aligned} pk+1=rk+1−pkTAkkpkTArk+1pk=rk+1−(Apk)Tpk(Apk)Trk+1pk=rk+1−(αrk−rk+1)Tpk(αrk−rk+1)Trk+1pk=rk+1−(αrk−rk+1)T(rk−βk−1pk−1)(αrk−rk+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 &= \frac{p_k^Tr_k}{p_k^tAp_k}\\ &=\frac{(r_k-\beta_{k-1}p_{k-1})^Tr_k}{p_k^tAp_k}\\ &=\frac{r_k^Tr_k}{p_k^TAp_k^T} \end{aligned} αk=pktApkpkTrk=pktApk(rk−βk−1pk−1)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}&=b-Ax_{k+1}\\ &=b-A(x_k+\alpha_kp_k)\\ &=b-Ax_k-\alpha_kAp_k\\ &=r_k-\alpha_kAp_k \end{aligned} rk+1=b−Axk+1=b−A(xk+αkpk)=b−Axk−αkApk=rk−αkApk
最终的推导结束。整理为如下的伪代码
五、伪代码
r 0 = b − A x 0 r_0=b-Ax_0 r0=b−Ax0
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
共轭梯度法详细推导分析相关推荐
- 关于罗德里格斯公式(Rodrigues‘sFormula)的详细推导过程
关于罗德里格斯公式[Rodrigues'sFormula]的详细推导过程 1 旋转向量 2 罗德里格斯公式 2.1 罗德里格斯公式定义 2.2 罗德里格斯公式推导 3 旋转矩阵到旋转向量的转换 1 旋 ...
- 机器学习强基计划6-2:详细推导马尔科夫随机场(MRF)及其应用(附例题)
目录 0 写在前面 1 无向概率图 2 马尔科夫随机场 3 马尔科夫独立性 4 例题分析 0 写在前面 机器学习强基计划聚焦深度和广度,加深对机器学习模型的理解与应用."深"在详细 ...
- 【转】流体动力学控制方程(详细推导)
[转]流体动力学控制方程(详细推导) diyhoos 原文链接:https://blog.csdn.net/qq_42020563/article/details/80940387 CFD建立在流体力 ...
- 【模拟电子技术Analog Electronics Technology 27】—— 非正弦波发生电路参数的详细计算分析(阈值电压和周期)
文章目录 1.矩形波发生电路 1.1产生矩形波的过程分析 1.2 参数计算(阈值电压,以及周期T[重要!]) 2. 占空比可调的矩形波发生电路 3.三角波发生电路 4.锯齿波发生电路 1.矩形波发生电 ...
- SMO算法详细推导(Sequential Minimal Optimization)
本文针对一般性的"软判断的核函数的对偶问题的SVM",形如下式: 上式问题所在:当采样点 xix_ixi 选取50000个点时,则基于核函数变量Θ(xi,xj)\bm{\Thet ...
- 【深度学习】感知器、线性神经网络案例应用、BP神经网络算法详细推导
感知器.线性神经网络.BP神经网络及手写数字识别 1. 单层感知器 1.1 感知器的介绍 1.2 感知器的学习规则 1.3 感知器单输入输出示例 1.4 学习率 η\etaη 1.5 模型训练收敛条件 ...
- 一文让你彻底搞懂最小二乘法(超详细推导)
要解决的问题 在工程应用中,我们经常会用一组观测数据去估计模型的参数,模型是我们根据先验知识定下的.比如我们有一组观测数据(xi,yi)(x_i,y_i)(xi,yi)(一维),通过一些数据分析我 ...
- 【转】卡尔曼滤波算法详细推导(相当值得一看)
转载自 卡尔曼滤波算法详细推导 这一篇对预备知识的介绍还是很好的,过程与原理讲解也很到位,应该是目前看到中文里最好的讲解了. 一.预备知识 1.协方差矩阵 是一个维列向量,是的期望 ...
- python 高斯烟羽模型_GPR(高斯过程回归)详细推导
GPR(高斯过程回归)详细推导 一.综述 GPR来源于线性模型,有两种方式可以推导出GPR,一种是weight space view,另外一种是function space view.两者考察方式假设 ...
最新文章
- Java8内存模型—永久代(PermGen)和元空间(Metaspace)
- 关于学习Python的一点学习总结(22->相关的迭代操作)
- 微信小程序实例开发教程之知乎新闻
- protobuf简单序列化反序列化示例
- 2007标注没有文字_应用技巧:CAD在机械工程制图中尺寸标注
- python字典有什么用_在Python中使用范围作为字典键,我有什么选...
- 一致性Hash(Consistent Hashing)原理剖析
- js常用reduce方法
- 贵州2021高考体考成绩查询,2021年贵州体育专业考试成绩查询网址:http://www.eaagz.org.cn/...
- Git 中文详细安装教程01(安装篇)
- 用JavaScript将数字转换为大写金额
- 【Flink】Flink ChildFirstClassLoader loadClassWithoutExceptionHandling 空指针
- 最新log4j2 远程代码执行漏洞(紧急扩散)
- 做游戏,学编程(C语言) 11 2048
- C/S模型之TCP协议
- 怎么提高计算机的桌面性能,电脑卡顿怎么办?4个小技巧帮你解决问题
- 笔记本装双系统!win10+Linux!所有的坑自己一个个爬过来,纪念一下。
- 苹果可以用android流量监控,iPhone怎么看流量统计?
- 勾号、叉号、圈号的收集
- AI人工智能可以做哪些课题的毕业设计毕设
热门文章
- BufferReader
- 迁移学习篇之如何迁移经典CNN网络-附迁移学习Alexnet,VGG,Googlenet,Resnet详细代码注释和方法-pytorch
- SE-Net:Squeeze-and-Excitation blocks
- autojs遍历当前页面所有控件_纯前端表格控件SpreadJS V14.0发布:组件化编辑器+数据透视表 - 葡萄城开发工具...
- 启动monitor白屏
- 安卓pdf阅读器_带笔的小尺寸BOOX Nova Pro电子书阅读器来了!
- taotao.sql文件(免积分分享)
- MySql8安装错误信息:The service already exists!
- 和Vue来一场美丽的邂逅
- 2012年9月ITbrand电信业4G技术品牌排行榜