L-BFGS-B局部極小化算法在Fortran中的使用簡例
L-BFGS-B局部極小化算法在Fortran中的使用簡例
文/Sobereva@北京科音 2020-Mar-3
本文簡要介紹一下L-BFGS-B局部極小化算法的常識，并通過實例介紹怎么在Fortran中使用這個方法解決實際問題。
1 相關知識
牛頓法是最重要的做多元函數局部極小化的算法之一，BFGS（Broyden–Fletcher–Goldfarb–Shanno algorithm）是準牛頓法的一種實現，通過近似方式估計Hessian矩陣的逆矩陣而避免了牛頓法那樣需要對其精確計算帶來的不便和耗時。L-BFGS（limited-memory BFGS）不需要像BFGS那樣需要構造和記錄完整的Hessian的逆矩陣，使得消耗內存量隨被優化的變量數只是線性增長，因此尤為適合用于變量數甚巨的場合。L-BFGS-B（B代表bound）則進一步使L-BFGS支持了極小化過程中對變量施加約束。相關知識看https://en.wikipedia.org/wiki/Limited-memory_BFGS。目前BFGS或L-BFGS已經被應用得極為普遍，諸如量子化學程序里做幾何優化主要就是基于BFGS的思想。
從耗時來說，牛頓法>BFGS>L-BFGS，而從優化效率來說（達到同樣精度所需步數），也是牛頓法>BFGS>L-BFGS（但差距不算特別大）。
L-BFGS-B的作者直接提供了實現L-BFGS-B算法的Fortran 77的代碼，見http://users.iems.northwestern.edu/~nocedal/software.html#lbfgs。但由于Fortran 77的限制，其子程序調用起來相當麻煩。為了使得其使用盡可能簡便，我寫了一個Fortran 90的界面。下面通過一個實際例子，講解怎么具體使用L-BFGS-B對一個簡單的非線性函數做極小化。
2 例子
L-BFGS-B源文件（.f后綴的），我寫的界面（LBFGSB_module.f90）以及本文的例子代碼（test_LBFGSB.f90）都可以在這里得到：http://www.shanxitv.org/attach/538/file.rar。將文件包里所有文件都一起編譯即可得到可執行文件。可見L-BFGS-B子程序利用了BLAS和LINPACK庫的部分子程序（后者已經被LAPACK庫取代）。
此例我們要對下面的二元函數做極小化，對兩個變量的解析導數公式也給出了：
雖然L-BFGS-B支持設置約束條件（即各個變量上下限），但是本例我們不考慮約束條件，相當于只是做L-BFGS。
本例的代碼如下，我已經寫得盡可能容易理解。
!A code to illustrate the use of L-BFGS-B local minimization routine
!Written by Tian Lu (sobereva@sina.com), 2020-Mar-3
program test_LBFGSB
use LBFGSB_module
implicit real*8(a-h,o-z)
character*60 task
real*8,allocatable :: X(:),g(:)
nvar=2
allocate(X(nvar),g(nvar))
maxcyc=50 !Maximum cycles
ftol=1D-7 !Tolerance of function variation
gtol=1D-5 !Tolerance of absolute maximum gradient
X=0 !Initial guess
task="START"
icyc=0
do while(task=="START".or.task(1:2)=="FG".or.task=="NEW_X")
call LBFGS(X,f,g,task,ncorr=3)
if (task(1:2)=="FG") then !Need evaluate function value and gradient
icyc=icyc+1
if (icyc>maxcyc) exit
call fval(nvar,X,f)
call fgrad(nvar,X,g)
write(*,"(' Iteration',i10)") icyc
if (icyc==1) then
write(*,"(' Function value:',f16.10)") f
else
fvar=f-fold
write(*,"(' Function value:',f16.10,' Variation:',f16.10)") f,fvar
end if
fold=f
gmax=maxval(g)
grms=dsqrt(sum(g**2))
write(*,"(' Maximum gradient:',f16.10,' RMS gradient:',f16.10)") gmax,grms
if (icyc>1) then
if (abs(fvar)<ftol.and.abs(gmax)<gtol) exit
end if
else if (task=="NEW_X") then !Variables have been updated
write(*,*) "New variables:"
write(*,"(6f12.8)") X(:)
end if
write(*,*)
end do
if (task=="ERROR") then
write(*,*) "Error was encountered during L-BFGS-B running!"
else if (icyc==maxcyc+1) then
write(*,"(a,i6,a)") " L-BFGS-B iteration was unconverged after",icyc," cycles!"
else
write(*,"(/,a,i6,a)") " Convergence criterions have been satisfied after",icyc," cycles!"
write(*,*) "Final variable values:"
write(*,"(6f12.8)") X(:)
write(*,"(a,f16.10)") " Final function value:",f
end if
read(*,*)
end program
subroutine fval(nvar,X,f)
integer nvar
real*8 :: X(nvar),f
f=(X(1)-2)**2+(X(2)-3.5D0)**2+X(1)*X(2)
end subroutine
subroutine fgrad(nvar,X,g)
integer nvar
real*8 :: X(nvar),g(nvar)
g(1)=2*(X(1)-2)+X(2)
g(2)=2*(X(2)-3.5D0)+X(1)
end subroutine
運行輸出如下
Iteration 1
Function value: 16.2500000000
Maximum gradient: -4.0000000000 RMS gradient: 8.0622577483
Iteration 2
Function value: 9.6185114825 Variation: -6.6314885175
Maximum gradient: -2.1394789812 RMS gradient: 5.2254408980
New variables:
0.49613894 0.86824314
Iteration 3
Function value: 4.3223858291 Variation: -5.2961256533
Maximum gradient: 0.8008108585 RMS gradient: 0.9223347926
New variables:
1.01974264 2.76132558
Iteration 4
Function value: 4.0640697287 Variation: -0.2583161005
Maximum gradient: 0.5716770197 RMS gradient: 0.5931409754
New variables:
0.76715766 3.03736169
Iteration 5
Function value: 3.9166672452 Variation: -0.1474024835
Maximum gradient: -0.0012552698 RMS gradient: 0.0018594564
New variables:
0.33283721 3.33295376
Iteration 6
Function value: 3.9166666667 Variation: -0.0000005785
Maximum gradient: -0.0000002118 RMS gradient: 0.0000003765
New variables:
0.33333330 3.33333320
Iteration 7
Function value: 3.9166666667 Variation: -0.0000000000
Maximum gradient: 0.0000000000 RMS gradient: 0.0000000000
Convergence criterions have been satisfied after 7 cycles!
Final variable values:
0.33333333 3.33333333
Final function value: 3.9166666667
可見收斂非常順利。函數變化程度、最大梯度、方均根梯度都隨迭代迅速下降。得到的極小點處的兩個變量值為0.33333333和3.33333333，與手算得到的解析解1/3和10/3精確吻合，故運行得非常成功。并且可見極小點處的函數值是3.9166666667。
下面解釋一下這個例子代碼
開頭必須寫use LBFGSB_module，這樣才能用我寫的叫做LBFGS的子程序。這個子程序是原版L-BFGS-B代碼的setulb子程序的一個wrapper，在迭代過程中不斷被調用，通過接收當前位置的函數值和梯度后給出新的位置。
fval是用戶寫的計算函數值的子程序。fgrad是用戶寫的計算函數梯度矢量的子程序。nvar是變量數目，X是傳入的變量矢量。如果被極小化的函數沒法算解析導數，用有限差分計算數值導數來定義fgrad也完全可以。
本例我們將變量的初值設為了(0,0)。由于此函數就一個極小點，所以初值非常隨意。
maxcyc變量設的是迭代次數上限，每次迭代都需要計算一次函數值和梯度。ftol和gtol變量分別設定函數值變化的絕對值和最大梯度的絕對值的收斂限，由代碼可見它倆同時滿足時迭代就宣告收斂。
task是運行狀態字符串變量，必須是60個字符。初始狀態必須是START。當LBFGS子程序返回的task是"FG"時，代表現在需要對當前的位置計算函數值和梯度并傳入LBFGS；當返回值是"NEW_X"，說明這次LBFGS更新了位置矢量X。實際上，整個do循環過程中，是"FG"和"NEW_X"狀態不斷交替，因此每兩輪循環計算一次函數值和梯度（算作一次迭代）。
LBFGS子程序的ncorr是可選參數，它控制校正數，數值通常建議在3~10之間，取值影響收斂效率，設多大最合適視具體問題而定，可以進行實測。
如果想在極小化過程中施加約束，需要向LBFGS子程序傳入記錄各個變量下限和上限的Xlower和Xupper矢量，都是nvar個元素。還要傳入boundstat矢量設置各變量的約束狀態，含義見LBFGSB_module.f90里的注釋。
對LBFGS子程序還可以設置info變量，用于控制輸出信息的詳細程度。此例沒有設它，相當于用了默認的info=-1，此時LBFGS子程序它自己在運行過程中不會輸出任何信息。可用選項見LBFGSB_module.f90里的注釋。
順帶一提，之前在《無需導數的局部極小化算法NEWUOA在FORTRAN中的使用簡例》（http://www.shanxitv.org/536）一文中筆者也用過這個例子。NEWUOA方法的好處是不需要導數，但達到收斂所需的迭代次數高于L-BFGS-B，穩健程度也不及它。但對于沒有解析導數的情況，用NEWUOA算法從實際耗時上可能更劃算。