2019-04-10

共轭梯度

参考

主要参考这篇文章,非常推荐查看,文章最后还有个关于共轭梯度的图
Conjugate gradient method wiki
Gram-Schmidt方法 wiki
在介绍共轭梯度之前,我们先介绍一下Gram-Schmidt过程(革兰氏-施密特方法),看它如何对向量组进行正交归一化。

首先定义投影操作:{\displaystyle proj_u(v)=\frac{\langle u,v\rangle}{\langle u,u \rangle}u}
其中\langle u,v\rangle表示向量u,v的在R^n上的内积。该运算符将向量v正交投影到向量u所跨越的线上。如果u=0,我们定义proj_u(v)=0

Gram-Schmidt过程的工作原理如下:
u_1=v_1
u_2=v_2-proj_{u_1}(v_2)
u_3=v_3-proj_{u_1}(v_3)-proj_{u_2}(v_3)
\cdots
u_k=v_k-\sum\limits^{k-1}_{j=1}proj_{u_j}(v_k)

该方法进行如下运算:为了计算u_i,它首先将v_i正交投影到u_1,\cdots,u_{i-1}上,然后从向量u_i中减去为v_i与每个投影之间的差,从而保证与空间中的所有向量(u_1,\cdots,u_{i-1})正交。

图来源wiki
最终原理图如下所示:

optimization图23.gif

什么是共轭?
对矩阵A\in R^{n\times n}对称正定,而p_1,p_2,\cdots, p_mR^n中任意一组非零向量,若 p_i^TAp_j=0(i\ne j),则称向量组p_1,p_2,\cdots, p_m是相对于A共轭的。

由于矩阵是A是对称(A^T=A)且正定的,则关于其的共轭内积可定义为:
{\displaystyle \langle u,v\rangle_A:=\langle Au,v\rangle}=\langle u,A^Tv\rangle=\langle u,Av\rangle=u^TAv
(由于我们将要处理的是共轭而不是正交(A=I时为正交),所以在要运用Gram-Schmidt方法时需要运用共轭形式的内积。)

进入正题:

模型建立
假定待优化的函数:{\displaystyle \underset{x}{min}f(x)=\frac{1}{2}x^TAx-b^Tx}

(其中x为待优化变量,A为正定矩阵,b为已知变量。)

下标k表示优化步数,负梯度为:r_k=b-Ax_k

假设最优变量为x^*,则优化问题可变为求方程Ax^*=b的解。梯度r也可以称作为每一步的残差。
算法推导
优化方向确定:

假定第一次优化方向为初始负梯度方向:p_1=r_0

每一次优化方向与之前的优化方向正交,采用Gram-Schmidt方法进行向量正交化,每次优化方向根据当前步的梯度得出:
{\displaystyle p_k=r_k-\sum\limits^{k-1}_{i=1} \frac{p_i^TAr_k}{p^T_iAp_i}p_i}

再令{\displaystyle \beta_{ik}=\frac{p^T_iAr_k}{p^T_iAp_i}}(这个\beta计算出来是个实数),则有:{\displaystyle p_k=r_k-\sum\limits^{k-1}_{i=1}\beta_{ik} p_i}

优化步长确定:

假定第k步的优化步长为\alpha_k
f(x_{k+1})=f(x_k+\alpha_kp_k)求导令导数为零有:{\displaystyle \alpha_k= \frac{p_k^Tr_k}{p^T_kAp_k}}(行搜索)。

共轭梯度其实到这里已经证明结束了,但一般选择其化简式。所以接下来对等式进行化简:
误差e_k定义为x_k与最优变量x^*的差值:
{\displaystyle e_k=x^*-x_k=\sum\limits^n_{i=1}\alpha_ip_i-\sum\limits^k_{i=1}\alpha_ip_i=\sum\limits^n_{i=k+1}\alpha_ip_i}
对残差r_k有:
r_k=b-Ax_k=Ax^*-Ax_k=Ae_k

在得到最终的等式前,证明下面两个性质:

1.残差与之前的所有的方向正交

证明:
任取一个方向p_j,并且j\le k

{\displaystyle p_j^Tr_k=p_j^TAe_k=\sum\limits^n_{i=k+1}\alpha_i p_j^TAp_i}=0

2.残差与之前的所有的残差正交
如果残差与搜索方向形成向量基正交,那么残差也会形成向量基。这意味着任何残差都与它先前的残差正交。
证明如下:
k+1\le j时:
{\displaystyle p_{k+1}=r_k-\sum\limits_{i\le k}\beta_{ik}p_i}
右乘r_j有
{\displaystyle p_{k+1}^Tr_j=r_k^Tr_j-\sum\limits_{i\le k}\beta_{ik}p_i^Tr_j}
0=r_k^Tr_j+0

基于残差的梯度步长:(求简化的\alpha)
由之前可得:r_{k+1}^Tr_k=0
又由于:r_{k+1}=b-Ax_{k+1}则:
(b-Ax_{k+1})^Tr_k=0
(b-A(x_k+\alpha_{k+1}p_{k+1}))^Tr_k=0
(b-Ax_k)^Tr_k-\alpha_{k+1}(Ap_{k+1}^T)r_k=0
r_k^Tr_k-\alpha_{k+1}(Ap_{k+1})^Tr_k=0
{\displaystyle \alpha_{k+1}=\frac{r_k^Tr_k}{p_{k+1}^TAr_k}} (A=A^T)
又因为:{\displaystyle p_{k+1}=r_k-\sum\limits_{i\le k}\beta_i p_i}
{\displaystyle \alpha_{k+1}=\frac{r_k^Tr_k}{p_{k+1}^TAp_{k+1}}}

正交化步骤仅取决于先前的搜索方向::(求简化的\beta)
x_{k+1}=x_k+\alpha_{k+1}p_{k+1}(更新x_k的值)
b-Ax_{k+1}=b-Ax_k-\alpha_{k+1}Ap_{k+1}
(\alpha_{k+1}的值为实数)
r_{k+1}=r_k-\alpha_{k+1}Ap_{k+1}
前乘r_j^T值,求新的\beta值:
r_j^Tr_{k+1}=r_j^Tr_k-\alpha_{k+1}r_j^TAp_{k+1}
\alpha_{k+1}r_j^TAp_{k+1}=-r_j^Tr_{k+1}+r_j^Tr_k
又因为:{\displaystyle \beta_{ik}=\frac{r_k^TAp_i}{p_i^TAp_i}}
\alpha_{k+1}\beta_{k+1,j}p_{k+1}^TAp_{k+1}=-r_j^Tr_{k+1}+r_j^Tr_k
\beta_{k+1,j}r_k^Tr_k=-r_j^Tr_{k+1}+r_j^Tr_k
{\displaystyle \alpha_{k+1}=\frac{r_k^Tr_k}{p^T_{k+1}Ap_{k+1}}}代替\alpha_{k+1}p_{k+1}^TAp_{k+1}有:

得到{\displaystyle \beta_{k+1,j}=\frac{-r_j^Tr_{k+1}+r_j^Tr_k}{r_k^Tr_k}}
变化一下得:{\displaystyle \beta_{k,j}=\frac{-r_j^Tr_{k}+r_j^Tr_{k-1}}{r_{k-1}^Tr_{k-1}}}
又因为:{\displaystyle p_{k+1}=r_k-\sum\limits_{i\le k}\beta_{ik}p_i=r_k-\sum\limits_{j\le k}\beta_{k,j}p_i}
可简化为:{\displaystyle p_{k+1}=r_k-\beta_{k,k}p_k}
其中:{\displaystyle \beta_{k,k}=\frac{-r_k^Tr_k}{r_{k-1}^Tr_{k-1}}}

伪代码:
在简化后:
r_0=b-Ax_0
p_1=r_0
k=1
repeat
{\displaystyle \alpha_k=\frac{r_{k-1}^Tr_{k-1}}{p_k^TAp_k}}
x_k=x_{k-1}+\alpha_kp_k
r_k=r_{k-1}-\alpha_kAp_k
{\displaystyle \beta_k=\frac{r_k^Tr_k}{r_{k-1}^Tr_{k-1}}}
p_{k+1}=r_k+\beta_kp_k
k=k+1
until r_{k-1}<\epsilon
(残差小于某个值便可以跳出来了。)
return x_{k-1}

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 217,185评论 6 503
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 92,652评论 3 393
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 163,524评论 0 353
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,339评论 1 293
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,387评论 6 391
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,287评论 1 301
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,130评论 3 418
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,985评论 0 275
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,420评论 1 313
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,617评论 3 334
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,779评论 1 348
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,477评论 5 345
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,088评论 3 328
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,716评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,857评论 1 269
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,876评论 2 370
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,700评论 2 354

推荐阅读更多精彩内容