在一般問(wèn)題的優(yōu)化中，最速下降法和共軛梯度法都是非常有用的經(jīng)典方法，但最速下降法往往以”之”字形下降，速度較慢，不能很快的達(dá)到最優(yōu)值，共軛梯度法則優(yōu)于最速下降法，在前面的某個(gè)文章中，我們給出了牛頓法和最速下降法的比較，牛頓法需要初值點(diǎn)在最優(yōu)點(diǎn)附近，條件較為苛刻。

算法來(lái)源：《數(shù)值最優(yōu)化方法》高立，P111

我們選用了64維的二次函數(shù)來(lái)作為驗(yàn)證函數(shù)，具體參見(jiàn)上書111頁(yè)。

采用的三種方法為：

共軛梯度方法（FR格式）、共軛梯度法（PRP格式）、最速下降法

# -*- coding: utf-8 -*-'''Created on Sat Oct 01 15:01:54 2016@author: zhangweiguo'''import sympy,numpyimport mathimport matplotlib.pyplot as plfrom mpl_toolkits.mplot3d import Axes3D as ax3import SD#這個(gè)文件里有最速下降法SD的方法，參見(jiàn)前面的博客#共軛梯度法FR、PRP兩種格式def CG_FR(x0,N,E,f,f_d): X=x0;Y=[];Y_d=[]; n = 1 ee = f_d(x0) e=(ee[0]**2+ee[1]**2)**0.5 d=-f_d(x0) Y.append(f(x0)[0,0]);Y_d.append(e) a=sympy.Symbol(’a’,real=True) print ’第%2s次迭代：e=%f’ % (n, e) while n<N and e>E: n=n+1 g1=f_d(x0) f1=f(x0+a*f_d(x0)) a0=sympy.solve(sympy.diff(f1[0,0],a,1)) x0=x0-d*a0 X=numpy.c_[X,x0];Y.append(f(x0)[0,0]) ee = f_d(x0) e = math.pow(math.pow(ee[0,0],2)+math.pow(ee[1,0],2),0.5) Y_d.append(e) g2=f_d(x0) beta=(numpy.dot(g2.T,g2))/numpy.dot(g1.T,g1) d=-f_d(x0)+beta*d print ’第%2s次迭代：e=%f’%(n,e) return X,Y,Y_ddef CG_PRP(x0,N,E,f,f_d): X=x0;Y=[];Y_d=[]; n = 1 ee = f_d(x0) e=(ee[0]**2+ee[1]**2)**0.5 d=-f_d(x0) Y.append(f(x0)[0,0]);Y_d.append(e) a=sympy.Symbol(’a’,real=True) print ’第%2s次迭代：e=%f’ % (n, e) while n<N and e>E: n=n+1 g1=f_d(x0) f1=f(x0+a*f_d(x0)) a0=sympy.solve(sympy.diff(f1[0,0],a,1)) x0=x0-d*a0 X=numpy.c_[X,x0];Y.append(f(x0)[0,0]) ee = f_d(x0) e = math.pow(math.pow(ee[0,0],2)+math.pow(ee[1,0],2),0.5) Y_d.append(e) g2=f_d(x0) beta=(numpy.dot(g2.T,g2-g1))/numpy.dot(g1.T,g1) d=-f_d(x0)+beta*d print ’第%2s次迭代：e=%f’%(n,e) return X,Y,Y_dif __name__==’__main__’: ’’’ G=numpy.array([[21.0,4.0],[4.0,15.0]]) #G=numpy.array([[21.0,4.0],[4.0,1.0]]) b=numpy.array([[2.0],[3.0]]) c=10.0 x0=numpy.array([[-10.0],[100.0]]) ’’’ m=4 T=6*numpy.eye(m) T[0,1]=-1;T[m-1,m-2]=-1 for i in xrange(1,m-1): T[i,i+1]=-1 T[i,i-1]=-1 W=numpy.zeros((m**2,m**2)) W[0:m,0:m]=T W[m**2-m:m**2,m**2-m:m**2]=T W[0:m,m:2*m]=-numpy.eye(m) W[m**2-m:m**2,m**2-2*m:m**2-m]=-numpy.eye(m) for i in xrange(1,m-1): W[i*m:(i+1)*m,i*m:(i+1)*m]=T W[i*m:(i+1)*m,i*m+m:(i+1)*m+m]=-numpy.eye(m) W[i*m:(i+1)*m,i*m-m:(i+1)*m-m]=-numpy.eye(m) mm=m**2 mmm=m**3 G=numpy.zeros((mmm,mmm)) G[0:mm,0:mm]=W;G[mmm-mm:mmm,mmm-mm:mmm]=W; G[0:mm,mm:2*mm]=-numpy.eye(mm) G[mmm-mm:mmm,mmm-2*mm:mmm-mm]=-numpy.eye(mm) for i in xrange(1,m-1): G[i*mm:(i+1)*mm,i*mm:(i+1)*mm]=W G[i*mm:(i+1)*mm,i*mm-mm:(i+1)*mm-mm]=-numpy.eye(mm) G[i*mm:(i+1)*mm,i*mm+mm:(i+1)*mm+mm]=-numpy.eye(mm) x_goal=numpy.ones((mmm,1)) b=-numpy.dot(G,x_goal) c=0 f = lambda x: 0.5 * (numpy.dot(numpy.dot(x.T, G), x)) + numpy.dot(b.T, x) + c f_d = lambda x: numpy.dot(G, x) + b x0=x_goal+numpy.random.rand(mmm,1)*100 N=100 E=10**(-6) print ’共軛梯度PR’ X1, Y1, Y_d1=CG_FR(x0,N,E,f,f_d) print ’共軛梯度PBR’ X2, Y2, Y_d2=CG_PRP(x0,N,E,f,f_d) figure1=pl.figure(’trend’) n1=len(Y1) n2=len(Y2) x1=numpy.arange(1,n1+1) x2=numpy.arange(1,n2+1) X3, Y3, Y_d3=SD.SD(x0,N,E,f,f_d) n3=len(Y3) x3=range(1,n3+1) pl.semilogy(x3,Y3,’g*’,markersize=10,label=’SD:’+str(n3)) pl.semilogy(x1,Y1,’r*’,markersize=10,label=’CG-FR:’+str(n1)) pl.semilogy(x2,Y2,’b*’,markersize=10,label=’CG-PRP:’+str(n2)) pl.legend() #圖像顯示了三種不同的方法各自迭代的次數(shù)與最優(yōu)值變化情況，共軛梯度方法是明顯優(yōu)于最速下降法的 pl.xlabel(’n’) pl.ylabel(’f(x)’) pl.show()

最優(yōu)值變化趨勢(shì)：

從圖中可以看出，最速下降法SD的迭代次數(shù)是最多的，在與共軛梯度（FR與PRP兩種方法）的比較中，明顯較差。

補(bǔ)充知識(shí)：python實(shí)現(xiàn)牛頓迭代法和二分法求平方根，精確到小數(shù)點(diǎn)后無(wú)限多位-4

首先來(lái)看一下牛頓迭代法求平方根的過(guò)程：計(jì)算3的平方根

如圖，是求根號(hào)3的牛頓迭代法過(guò)程。這里使用的初始迭代值（也就是猜測(cè)值）為1，其實(shí)可以為任何值最終都能得到結(jié)果。每次開(kāi)始，先檢測(cè)猜測(cè)值是否合理，不合理時(shí)，用上面的平均值來(lái)?yè)Q掉猜測(cè)值，依次繼續(xù)迭代，直到猜測(cè)值合理。

原理：現(xiàn)在取一個(gè)猜測(cè)值 a, 如果猜測(cè)值合理的話，那么就有a^2=x，即x/a=a ,x為被開(kāi)方數(shù)。不合理的話呢，就用表中的猜測(cè)值和商的平均值來(lái)?yè)Q掉猜測(cè)值。當(dāng)不合理時(shí)，比如 a>真實(shí)值，那么x/a<真實(shí)值，這時(shí)候取a 與 x/a 的平均值來(lái)代替a的話，那么新的a就會(huì)比原來(lái)的a要更接近真實(shí)值。同理有 a<真實(shí)值 的情況。于是，這樣不斷迭代下去最終是一個(gè)a不斷收斂到真實(shí)值的一個(gè)過(guò)程。于是不斷迭代就能得到真實(shí)值，證明了迭代法是正確的。

附上我的python代碼：

利用python整數(shù)運(yùn)算，python整數(shù)可以無(wú)限大，可以實(shí)現(xiàn)小數(shù)點(diǎn)后無(wú)限多位

#二分法求x的平方根小數(shù)點(diǎn)下任意K位數(shù)的精準(zhǔn)值，利用整數(shù)運(yùn)算 #思想：利用二分法，每次乘以10，取中間值，比較大小，從而定位精確值的范圍，將根擴(kuò)大10倍，則被開(kāi)方數(shù)擴(kuò)大100倍。 #quotient（商）牛頓迭代法：先猜測(cè)一個(gè)值，再求商，然后用猜測(cè)值和商的中間值代替猜測(cè)值，擴(kuò)大倍數(shù)，繼續(xù)進(jìn)行。

import mathfrom math import sqrt def check_precision(l,h,p,len1):#檢查是否達(dá)到了精確位 l=str(l);h=str(h) if len(l)<=len1+p or len(h)<=len1+p: return False for i in range(len1,p+len1):#檢查小數(shù)點(diǎn)后面的p個(gè)數(shù)是否相等 if l[i]!=h[i]: #當(dāng)l和h某一位不相等時(shí)，說(shuō)明沒(méi)有達(dá)到精確位 return False return True def print_result(x,len1,p): x=str(x) if len(x)-len1<p:#沒(méi)有達(dá)到要求的精度就已經(jīng)找出根 s=x[:len1]+'.'+x[len1:]+’0’*(p-len(x)+len1) else:s=x[:len1]+'.'+x[len1:len1+p] print(s) def binary_sqrt(x,p): x0=int(sqrt(x)) if x0*x0==x: #完全平方數(shù)直接開(kāi)方，不用繼續(xù)進(jìn)行 print_result(x0,len(str(x0)),p) return len1=len(str(x0))#找出整數(shù)部分的長(zhǎng)度 l=0;h=x while(not check_precision(l,h,p,len1)):#沒(méi)有達(dá)到精確位，繼續(xù)循環(huán) if not l==0:#第一次l=0，h=x時(shí)不用乘以10，直接取中間值 h=h*10 #l,h每次擴(kuò)大10倍 l=l*10 x=x*100 #x每次要擴(kuò)大100倍，因?yàn)槠椒? m=(l+h)//2 if m*m==x: return print_result(m,len1,p) elif m*m>x: h=m else: l=m return print_result(l,len1,p)#當(dāng)達(dá)到了要求的精度，直接返回l #牛頓迭代法求平方根def newton_sqrt(x,p): x0=int(sqrt(x)) if x0*x0==x: #完全平方數(shù)直接開(kāi)方，不用繼續(xù)進(jìn)行 print_result(x0,len(str(x0)),p) return len1=len(str(x0))#找出整數(shù)部分的長(zhǎng)度 g=1;q=x//g;g=(g+q)//2 while(not check_precision(g,q,p,len1)): x=x*100 g=g*10 q=x//g #求商 g=(g+q)//2 #更新猜測(cè)值為猜測(cè)值和商的中間值 return print_result(g,len1,p) while True: x=int(input('請(qǐng)輸入待開(kāi)方數(shù)：')) p=int(input('請(qǐng)輸入精度：')) print('binary_sqrt:',end='') binary_sqrt(x,p) print('newton_sqrt:',end='') newton_sqrt(x,p)

以上這篇基于Python共軛梯度法與最速下降法之間的對(duì)比就是小編分享給大家的全部?jī)?nèi)容了，希望能給大家一個(gè)參考，也希望大家多多支持好吧啦網(wǎng)。