共轭梯度法
共軛梯度法
最速下降法以及牛頓法都具有其自身的局限性。本文將要介紹的共軛梯度法是介于最速下降法與牛頓法之間的一種無約束優化算法,它具有超線性的收斂速度,而且算法結構簡單,容易編程實現。此外,根最速下降法相類似,共軛梯度法只用到了目標函數及其梯度值,避免了二階導數的計算,從而降低了計算量和存儲量,因此它是求解無約束優化問題的一種比較有效而使用的算法。
一、共軛方向發
共軛方向法的基本思想是在求解nnn維正定二次目標函數極小點時產生一組共軛方向作為搜索方向,在線搜索條件下算法之多迭代nnn步即能求得極小點。京故宮適當修正后共軛方向法可以推廣到求解一般非二次目標函數情形。下面先介紹共軛方向的概念。
定義1:設GGG是nnn階對稱正定矩陣,若nnn維向量組d1,d2,?,dm(m≤n)滿足d_1,d_2,\cdots,d_m(m\le n)滿足d1?,d2?,?,dm?(m≤n)滿足diTGdj=0,i≠jd_i^TGd_j=0,i\neq jdiT?Gdj?=0,i?=j,則乘d1,d2,?,dmd_1,d_2,\cdots,d_md1?,d2?,?,dm?是GGG共軛的。
顯然,向量組的共軛是正交的推廣,即當G=IG=IG=I(單位陣)時,上述定義變成了向量組正交的定義。此外,不難證明,對稱正定矩陣GGG的共軛向量組必然是線性無關的。
下面我們考慮求解正定二次目標函數極小點的共軛方向法。設
min?f(x)=12xTGx+bTx+c(1)\min f(x)=\frac{1}{2}x^TGx+b^Tx+c\tag{1} minf(x)=21?xTGx+bTx+c(1)
其中GGG是nnn階對稱正定陣,bbb為nnn為常向量,ccc為常數。我們有下面的算法:
算法1:共軛方向法
0. 給定迭代精度0≤?≤10\le\epsilon\le 10≤?≤1和初始點x0x_0x0?。計算g0=?f(x0)g_0=\nabla f(x_0)g0?=?f(x0?)。選取初始方向d0d_0d0?使得d0Tg0<0d_0^Tg_0\lt 0d0T?g0?<0。令k=0k=0k=0.
該算法的收斂性證明略過,如果感興趣,可以去查找相應的數值優化專著。這里直接給出結論,在精確線搜索下,算法1求解正定二次目標函數極小化問題,之多在nnn步內即可求得其唯一極小點。這種能在有限步內求得二次函數極小點的性質通常稱為二次終止性。
二、共軛梯度法
共軛梯度法是在每一步迭代利用當前點處的最速下降方向來生成關于凸二次函數fff的海塞陣GGG的共軛方向,并建立求fff在Rn\mathbb{R^n}Rn上的極小點的方法。這一方法最早是由Hesteness和Stiefel于1952年為求解正定線性方程組而提出來的,后經Fletcher等人研究并應用于無約束優化問題取得了豐富的成果,共軛梯度法也因此成為當前求解無約束優化問題的重要算法類。
設函數如(1)式所定義,則fff的梯度和海塞矩陣為
g(x)=?f(x)=Gx+b,G(x)=?2f(x)=G(2)g(x)=\nabla f(x)=Gx+b,\quad G(x)=\nabla^2 f(x)=G\tag{2} g(x)=?f(x)=Gx+b,G(x)=?2f(x)=G(2)
下面我們討論算法(1)中共軛方向的構造。我們取初始方向d0d_0d0?為初始點x0x_0x0?處的負梯度方向,即
d0=??f(x0)=?g0(3)d_0=-\nabla f(x_0)=-g_0 \tag{3} d0?=??f(x0?)=?g0?(3)
從x0x_0x0?出發沿d0d_0d0?方向進行線搜索得到步長α0\alpha_0α0?,令
x1=x0+α0d0x_1=x_0+\alpha_0d_0 x1?=x0?+α0?d0?
其中α0\alpha_0α0?滿足條件
?f(x1)Td0=g1Td0(4)\nabla f(x_1)^Td_0=g_1^Td_0 \tag{4} ?f(x1?)Td0?=g1T?d0?(4)
在x1x_1x1?處,用fff在x1x_1x1?的負梯度方向?g1-g_1?g1?與d0d_0d0?的組合來生成d1d_1d1?,即
d1=?g1+β0d0(5)d_1=-g_1+\beta_0d_0 \tag{5} d1?=?g1?+β0?d0?(5)
然后選取系數β0\beta_0β0?使d1d_1d1?與d0d_0d0?關于G共軛,即令
d1TGd0=0(6)d_1^TGd_0 = 0 \tag{6} d1T?Gd0?=0(6)
來確定β0\beta_0β0?,將(6)代入(5)得
β0=g1TGd0d0TGd0(7)\beta_0 = \frac{g_1^TGd_0}{d_0^TGd_0} \tag{7} β0?=d0T?Gd0?g1T?Gd0??(7)
由(2)得
g1?g0=G(x1?x0)=α0Gd0(8)g_1-g_0=G(x_1-x_0)=\alpha_0Gd_0 \tag{8} g1??g0?=G(x1??x0?)=α0?Gd0?(8)
故由(3)~(5)可得
g2Tg0=0,g2Tg1=0,d0Tg0=?g0Tg0,d1Tg1=?g1Tg1g_2^Tg_0=0,g_2^Tg_1=0,d_0^Tg_0=-g_0^Tg_0,d_1^Tg_1=-g_1^Tg_1 g2T?g0?=0,g2T?g1?=0,d0T?g0?=?g0T?g0?,d1T?g1?=?g1T?g1?
現假設已得到互相共軛得搜索方向d1,d2,?,dk?1d_1,d_2,\cdots,d_{k-1}d1?,d2?,?,dk?1?,精確線搜索得到得步長為α0,α1,?,αk?1\alpha_0,\alpha_1,\cdots,\alpha_{k-1}α0?,α1?,?,αk?1?,且滿足
{dk?1TGdi=0,i=0,1,?,k?2,diTgi=?giTgi,i=0,1,?,k?1,gkTgi=0,gkTdi=0,i=0,1,?,k?1.(9)\left\{ \begin{array}{rcl} d_{k-1}^TGd_i=0, &i=0,1,\cdots,k-2,\\ d_i^Tg_i=-g_i^Tg_i,&i=0,1,\cdots,k-1,\\ g_k^Tg_i=0,g_k^Td_i=0,&i=0,1,\cdots,k-1. \end{array} \right. \tag{9} ????dk?1T?Gdi?=0,diT?gi?=?giT?gi?,gkT?gi?=0,gkT?di?=0,?i=0,1,?,k?2,i=0,1,?,k?1,i=0,1,?,k?1.?(9)
現令
dk=?gk+βk?1dk?1+∑i=0k?1βk(i)di(10)d_k=-g_k+\beta_{k-1}d_{k-1}+\sum_{i=0}^{k-1}\beta_{k}^{(i)}d_i\tag{10} dk?=?gk?+βk?1?dk?1?+i=0∑k?1?βk(i)?di?(10)
其中βk?1,βk(i)(i=0,1,?,k?2)\beta_{k-1},\beta_k^{(i)}(i=0,1,\cdots,k-2)βk?1?,βk(i)?(i=0,1,?,k?2)得選擇要滿足
dkTGdi=0,i=0,1,?,k?1(11)d_k^TGd_i=0,i=0,1,\cdots,k-1\tag{11} dkT?Gdi?=0,i=0,1,?,k?1(11)
用diTG(i=0,1,?,k?1)d_i^TG(i=0,1,\cdots,k-1)diT?G(i=0,1,?,k?1)左乘(10)得
βk?1=gkTGdk?1dk?1TGdk?1,βk(1)=gkTGdidiTGdi,i=0,1,?,k?2(12)\beta_{k-1}=\frac{g_k^TGd_{k-1}}{d_{k-1}^TGd_{k-1}},\beta_k^{(1)}=\frac{g_k^TGd_i}{d_i^TGd_i},i=0,1,\cdots,k-2\tag{12} βk?1?=dk?1T?Gdk?1?gkT?Gdk?1??,βk(1)?=diT?Gdi?gkT?Gdi??,i=0,1,?,k?2(12)
類似于(8),我們有
gi+1?gi=G(xi+1?xi)=αiGdi,i=0,1,?,k?1g_{i+1}-g_i=G(x_{i+1}-x_i)=\alpha_iGd_i,i=0,1,\cdots,k-1 gi+1??gi?=G(xi+1??xi?)=αi?Gdi?,i=0,1,?,k?1
及
αiGdi=gi+1?gi,i=0,1,?,k?1(13)\alpha_iGd_i=g_{i+1}-g_i,i=0,1,\cdots,k-1\tag{13} αi?Gdi?=gi+1??gi?,i=0,1,?,k?1(13)
于是由歸納法假設(9)可得
βk(i)=gkTGdidiTGdi=gkT(gi+1?gi)diT(gi+1?gi)=0,i=0,1,?,k?2.\beta_k^{(i)}=\frac{g_k^TGd_i}{d_i^TGd_i}=\frac{g_k^T(g_{i+1}-g_i)}{d_i^T(g_{i+1}-g_i)}=0,i=0,1,\cdots,k-2. βk(i)?=diT?Gdi?gkT?Gdi??=diT?(gi+1??gi?)gkT?(gi+1??gi?)?=0,i=0,1,?,k?2.
于是,第k步得搜索方向為
dk=?gk+βk?1dk?1,(14)d_k=-g_k+\beta_{k-1}d_{k-1},\tag{14} dk?=?gk?+βk?1?dk?1?,(14)
其中βk?1\beta_{k-1}βk?1?由(12)確定,即
βk?1=gkTGdk?1dk?1TGdk?1(15)\beta_{k-1}=\frac{g_k^TGd_{k-1}}{d_{k-1}^TGd_{k-1}}\tag{15} βk?1?=dk?1T?Gdk?1?gkT?Gdk?1??(15)
同時有dkTgk=?gkTgkd_k^Tg_k=-g_k^Tg_kdkT?gk?=?gkT?gk?。這樣確定了一組由負梯度方向形成得共軛方向,而把沿著這組方向進行迭代得方向稱為共軛梯度法。其證明過程這里略過。
下面我們給出共軛梯度法求解無約束優化問題(1)極小點得算法步驟
算法2:共軛梯度法
0. 給定迭代精度0≤?<10\le\epsilon\lt 10≤?<1和初始點x0x_0x0?。計算g0=?f(x0)g_0=\nabla f(x_0)g0?=?f(x0?)。令k=0k=0k=0
KaTeX parse error: Unknown column alignment: s at position 29: …\begin{array}{ls?c} -g_k,&k=0,\\…
其中當k≥1k\ge 1k≥1時,βk?1\beta_{k-1}βk?1?由(15)確定
計算公式(15)是由Fletcher和Reeves給出得,故稱之為FR公式,算法2也稱之為FR共軛梯度法。除FR工詩外,尚有下列著名公式:
βk=gk+1Tgk+1?dkTgk,(Dixon公式)βk=gk+1Tgk+1dkT(gk+1?gk),(Dai?Yuan公式)βk=gk+1T(gk+1?gk)dkT(gk+1?gk),(Crowder?Wolfe公式)βk=gk+1T(gk+1?gk)gkTgk,(Polak,Ribiere,Ployak,PRP公式)\begin{aligned} \beta_k&=\frac{g_{k+1}^Tg_{k+1}}{-d_k^Tg_k}, (Dixon公式) \\ \beta_k&=\frac{g_{k+1}^Tg_{k+1}}{d_k^T(g_{k+1}-g_k)}, (Dai-Yuan公式) \\ \beta_k&=\frac{g_{k+1}^T(g_{k+1}-g_k)}{d_k^T(g_{k+1}-g_k)}, (Crowder-Wolfe公式) \\ \beta_k&=\frac{g_{k+1}^T(g_{k+1}-g_k)}{g_k^Tg_k}, (Polak,Ribiere,Ployak,PRP公式) \\ \end{aligned} βk?βk?βk?βk??=?dkT?gk?gk+1T?gk+1??,(Dixon公式)=dkT?(gk+1??gk?)gk+1T?gk+1??,(Dai?Yuan公式)=dkT?(gk+1??gk?)gk+1T?(gk+1??gk?)?,(Crowder?Wolfe公式)=gkT?gk?gk+1T?(gk+1??gk?)?,(Polak,Ribiere,Ployak,PRP公式)?
三、共軛梯度法得matlab實現
在共軛梯度法得實際使用中,通常在迭代n步或n+1步之后,重新選取負梯度方向作為搜索方向,我們稱之為再開始共軛梯度法。這是因為對于一般非二次函數而言,n步迭代后共軛梯度法產生得搜索方向往往不再具有共軛性。而對于大規模問題,常常每m(m<n或m?n)m (m<n或m\ll n)m(m<n或m?n)步就進行再開始。此外,當搜索方向不是下降方向時,也插入負梯度方向所作為搜索方向。
這里給出基于Armijo-rule非精確線搜索得再開始FR共軛梯度法得matlab程序。
利用程序求解無約束優化問題
min?x∈R2f(x)=100(x12?x2)2+(x1?1)2\min_{x\in\mathbb{R^2}}\quad f(x)=100(x_1^2-x_2)^2+(x_1-1)^2 x∈R2min?f(x)=100(x12??x2?)2+(x1??1)2
該問題由精確解x?=(1,1)T,f(x?)=0x^*=(1,1)^T,f(x^*)=0x?=(1,1)T,f(x?)=0
求解main函數
x0 = [0, 0]'; epsilon = 1e-4;[fmin, xmin] = frcg('func', 'gfunc', x0, epsilon); fprintf('frcg: fmin = %f, xmin = (%f, %f)\n', fmin, xmin(1), xmin(2));[x, f] = fminsearch('func', x0); fprintf('build-in search: fmin = %f, xmin = (%f, %f)\n', f, x(1), x(2));函數定義以及梯度求解
function f = func(x)f = 100 * (x(1)^2 - x(2))^2 + (x(1) - 1)^2;end function grad = gfunc(x)grad = [400 * x(1) * (x(1)^2-x(2))+2*(x(1)-1); ...-200 * (x(1)^2-x(2))];end求解結果為:
frcg: fmin = 0.000000, xmin = (0.999921, 0.999841) build-in search: fmin = 0.000000, xmin = (1.000004, 1.000011)總結
- 上一篇: 字符设备驱动高级篇1——新接口介绍
- 下一篇: 扫地机器人滤网顺序_1分钟小课堂:扫地机