? ? ? ?我們知道，除了微分方程，不少物理規(guī)律都可以用極值的形式來(lái)描述，比如光傳播時(shí)的光程極值原理、運(yùn)動(dòng)體系的最小作用量原理、材料計(jì)算中的密度泛函原理、穩(wěn)定結(jié)構(gòu)的自由能極小等。而在這些問(wèn)題中，我們都需要面對(duì)多變量問(wèn)題的極值求解。

? ? ? ? 在數(shù)值計(jì)算中，優(yōu)化相關(guān)的內(nèi)容非常豐富，這里先介紹常見(jiàn)的局域優(yōu)化，想了解更多的朋友可以參考《最優(yōu)化導(dǎo)論》，作者: Edwin K. P. Chong / Stanislaw H. Zak，An Introduction to Optimization,Foulth Edition，譯者: 孫志強(qiáng) / 白圣建 / 鄭永斌 / 劉偉。

（1） 函數(shù)極值和解方程（組）的關(guān)聯(lián)

? ? ? ?對(duì)于多變量函數(shù)，函數(shù)極值和之前提到的解方程（組）是互相關(guān)聯(lián)的。因此，函數(shù)極值的計(jì)算（尤其是局域優(yōu)化的思路）通常和初始值的選擇有關(guān)，而程序通常在梯度達(dá)到閾值時(shí)停止。如果要得到更好的極值，通常需要全局搜索的算法（比如模擬退火、遺傳-進(jìn)化算法、微粒群算法等），這個(gè)在后續(xù)的內(nèi)容再介紹。

（2） 庫(kù)函數(shù)fminbnd和fminunc的使用

? ? ? ?考慮一元函數(shù)?cos(3x).exp(-x)-1)在[-2,4]區(qū)間，從圖像曲線看，有3個(gè)極小值點(diǎn)。在定義好目標(biāo)函數(shù)fun?=?@(x)?cos(3*x).*exp(-x)-1?后，如果用fminbnd函數(shù)，語(yǔ)法是fminbnd(fun,-2,4)，這里的-2和4就是指定的區(qū)間，返回結(jié)果是?3.0343，是其中一個(gè)極小值對(duì)應(yīng)的自變量位置。如果用fminunc函數(shù)，fminunc(fun,-2)?這里的-2是初始點(diǎn)，取-2，1，4等作為初始點(diǎn)，分別得到極值位置-1.1544，0.9399，3.0343，從圖像上也是很好理解的，就是初始點(diǎn)接近，而且中間梯度單調(diào)，自然落入附近的極值，這是局域優(yōu)化常見(jiàn)的現(xiàn)象。

? ? 此外，fminunc?可以用于多變量函數(shù)。比如函數(shù) fun?=?@(x)?x(1).^2+x(2).^2-x(1)-x(2);?用 fminunc?(fun,?[1?1]) 計(jì)算極值，這里[1 1]是x_1和x_2的初值，結(jié)果返回?：?0.5000???0.5000，和我們解析推導(dǎo)結(jié)果一致。如果處理帶約束的問(wèn)題，需要用fmincon函數(shù)，在matlab里是自帶的，在octave中需要通過(guò)?pkg?load?optim?提前加載。

（3） 牛頓法、擬牛頓法

? ? ? 假定在點(diǎn)x_k附近進(jìn)行二階泰勒展開(kāi)，得到q(x)是原函數(shù)f(x)的近似，進(jìn)一步對(duì)q(x)求導(dǎo)，可以得到x_{k+1}和x_k的迭代公式。

? ? ? ?注意，拓展到高維時(shí)，f?'?是梯度，f?''?對(duì)應(yīng)海森矩陣。和之前的解方程是類似的。由于該方法每次要計(jì)算海森矩陣，20世紀(jì)60年代，Davidon-Fletcher-Powell提出DFP方法來(lái)簡(jiǎn)化，70年代，Broyden-Fletcher-Goldfard-Shanno等人提出改進(jìn)方法，這類方法統(tǒng)稱擬牛頓法。

? ? ?此外，牛頓法對(duì)體系的二階導(dǎo)數(shù)有要求，如果類似拋物線等情況，迭代很快指向體系的平衡位置，但有時(shí)會(huì)偏離平衡位置越來(lái)越遠(yuǎn)（有興趣的朋友可以測(cè)試一下李納-瓊斯勢(shì)的情況）。

? ? ? 為求函數(shù)（單變量、多變量）極值，在局域優(yōu)化方法中，關(guān)鍵的是搜索方向和移動(dòng)步長(zhǎng)。在牛頓法中，可以認(rèn)為方向和步長(zhǎng)都乘在一起了，該方法最大的問(wèn)題就是每次要求二階導(dǎo)數(shù)。最速下降法的搜索方向是負(fù)梯度方向，步長(zhǎng)通過(guò)一元函數(shù)極值來(lái)確定，假設(shè)前進(jìn)步長(zhǎng)是alpha，新的位置是當(dāng)前位置加上alpha乘以負(fù)梯度，求一元函數(shù)的極小值確定alpha，然后得到新的位置；之后不斷迭代，在達(dá)到局域極小附近，梯度向量對(duì)應(yīng)的模逐漸減小，可以通過(guò)設(shè)置閾值使程序停止迭代。

? ? ? 附錄程序是用最速下降法求函數(shù)的最小值。這里選擇初始x=2.5,計(jì)算梯度gx和函數(shù)值y，如果gx向量的模沒(méi)有達(dá)到閾值，進(jìn)入循環(huán)，這里每次步長(zhǎng)alpha從2開(kāi)始，或者更大的數(shù)，之后不斷減半，直到對(duì)應(yīng)的函數(shù)值小于當(dāng)前數(shù)值。這里當(dāng)前函數(shù)值是y，更新后是ny。在初始點(diǎn)附近的3個(gè)點(diǎn)做插值（這3個(gè)點(diǎn)滿足第二點(diǎn)比另外兩點(diǎn)數(shù)值?。?，找到使函數(shù)值最小的步長(zhǎng)，更新自變量x，然后判斷該點(diǎn)梯度是否達(dá)到閾值。

? ? ? ?在高維問(wèn)題中（上圖右邊是二維情況），最速下降法在接近極小值附近會(huì)出現(xiàn)鋸齒形的迭代路徑，影響計(jì)算效率。共軛梯度法可以很好解決該問(wèn)題，搜索方向在梯度方向上加了修正，保證每次搜到方向關(guān)于海森矩陣共軛；通過(guò)Hestenes-Stiefel公式消去迭代公式中的海森矩陣，可以使效率極大提高。此外，對(duì)于最速下降法，沿著梯度做一維搜索步長(zhǎng)也是耗時(shí)的，可以根據(jù)梯度分量的大小適當(dāng)變步長(zhǎng)即可，在多數(shù)問(wèn)題時(shí)，評(píng)估目標(biāo)函數(shù)是很慢的。

附錄：

A. 牛頓法求函數(shù)極值

clear?

x(1)=2;

for i=2:10

? f1=4*x(i-1)^3-3;

? f2=12*x(i-1)^2;

? x(i)=x(i-1)-f1/f2;

end

y=x.^4-3*x-10;

plot(x,y,'o-')

hold on?

clear

x=linspace(-2,2);

y=x.^4-3*x-10;

plot(x,y)

B. 最速下降法求函數(shù)極值

clear

x=0.8;y=1/x^12-2/x^6;i=1;gx=-12/x^13+12/x^7;?

? f(i)=1/x^12-2/x^6;

? d(i)=x;

while abs(gx)>0.001&&i<100 %梯度向量模的閾值

? alpha=0.01;

? tp=[];ys=[];

? tp(1)=0;ys=y;

? x=x-gx*alpha;? %計(jì)算更新后的x和y

? gx=-12/x^13+12/x^7;?

? i=i+1;

? f(i)=1/x^12-2/x^6;

? d(i)=x;

end

plot(d,f,'o-')

hold on?

clear

x=linspace(0.8,2);

y=1./x.^12-2./x.^6;

plot(x,y)