一、首先把两个作者 JORGE NOCEDAL 给的文件看明白sdrive.f和lbfgs.f
sdrive.f - is a simple driver
lbfgs.f - contains all the supporting routines
1、lbfgs.f文件包括LBFGS算法和运行函数等
主体函数为：
subroutine lbfgs(n, m, x, f, g, diagco, diag, iprint, eps, xtol, w, iflag)
int n, m, iprint[2], iflag;
double x[n], g[n], diag[n], w[n*(2*m+1)+2*m];
double eps, xtol;
bool diagco;
使用LBFGS方法，能够有效地解决大规模变量的问题。算法开始定义一个对角矩阵Hk0，之后应用多次迭代，每次迭代使用前M个BFGS的结果更新这个矩阵，直到得到Hk矩阵，这个矩阵近似等于海塞矩阵的逆矩阵。这里的M是用户指定的，会影响到运行中所需的存储空间大小。用户同样可以提供一个不满足默认选择的初始对角矩阵。算法的详细描述可以阅读"On the limited memory BFGS method C for large scale optimization" by D. Liu and J. Nocedal。
用户需要自己计算给定的函数F的值和梯度G，为了便于用户控制这些计算可以使用逆向交流（这个还没看到，看到再写？），通过iflag实现重复调用函数运算。
每次迭代都会定义步长，定义使用CSRCH的微变换MCVSRCH。
调用语句为call lbfgs(n, m, x, f, g, diagco, diag, iprint, eps, xtol, w, iflag);
下面详细解释参数：
n:整型变量，变量个数 n>0
m:整型变量，BFGS更新时使用的修正个数 3<=m<=7
x:双精度浮点数组变量，长度为n，初始化为可能的解得向量 X=(x1, x2,..., xn)，程序结束时存储最优解结果向量
f:双精度浮点变量，初始化和iflag=1时为X带入的函数值
g:双精度浮点数组变量，长度为n，初始化呵iflag=1时为带入X后的函数的梯度
diagco:逻辑变量，如果用户希望每次迭代前提供初始化的对角矩阵设定为true，否则使用默认矩阵设为false。如果设为true， 每次迭代返回iflag=2且初始对角矩阵需要通过数组diag设置。
diag:双精度浮点数组，长度为n，如果diagco=true，初始化和iflag=2时，用户需要输入初始矩阵的数值，且每个值必须非负
iprint:整型数组变量，长度为2，用户设定
iprint[0]表示输出的频率--<0不生成输出，=0只在第一次和最后一次迭代后输出，>0每次迭代输出
iprint[1]表示输出类型----=0输出迭代次数，函数评估，函数值，归一后的梯度值和步长
=1除了以上之外，还输出最初点的变量向量和梯度向量
=2除了以上之外，还输出结果变量向量
=3除了以上之外，还输出结果梯度向量
eps:正的双精度浮点变量，用户设定，定义找到解终止迭代的误差范围 ||g||<epsmax(1,||x||);
xtol:正的双精度浮点变量，用户设定，机器精度，迭代在这个精度以下时自动终止
w:双精度浮点数组变量，长度为n(2m+1)+2m，作为LBFGS的运行空间
iflag:整型变量，如果返回结果为负表述出现错误，返回0表示无错终止，返回1用户必须给函数f和梯度g赋值，返回2用户需要提供初始对角矩阵
下面的错误信息就懒得翻译了，比较容易理解反正
The following negative values of IFLAG, detecting an error,
C are possible:
C
C IFLAG=-1 The line search routine MCSRCH failed. The
C parameter INFO provides more detailed information
C (see also the documentation of MCSRCH):
C
C INFO = 0 IMPROPER INPUT PARAMETERS.
C
C INFO = 2 RELATIVE WIDTH OF THE INTERVAL OF
C UNCERTAINTY IS AT MOST XTOL.
C
C INFO = 3 MORE THAN 20 FUNCTION EVALUATIONS WERE
C REQUIRED AT THE PRESENT ITERATION.
C
C INFO = 4 THE STEP IS TOO SMALL.
C
C INFO = 5 THE STEP IS TOO LARGE.
C
C INFO = 6 ROUNDING ERRORS PREVENT FURTHER PROGRESS.
C THERE MAY NOT BE A STEP WHICH SATISFIES
C THE SUFFICIENT DECREASE AND CURVATURE
C CONDITIONS. TOLERANCES MAY BE TOO SMALL.
C
C
C IFLAG=-2 The i-th diagonal element of the diagonal inverse
C Hessian approximation, given in DIAG, is not
C positive.
C
C IFLAG=-3 Improper input parameters for LBFGS (N or M are
C not positive).
on the driver:
程序要包括一个全局数据块lb2
common /lb3/mp, lp, gtol, stpmin, stpmax;
mp:默认为6的整型变量
lp:默认为6的整型变量
gtol:双精度浮点变量默认为0.9, 可以设为较低的数=0.1，不能小于1.D-0.4
stpmin stpmax:非负双精度浮点变量，定义行搜索的步长的上下界，默认为1.D-20~1.D+20
machine dependencies:xtol, stpmin, stpmax
general information:
直接调用:daxpy, ddot, lb1, mcsrch
输入输出:无输入，输出MP的特征信息和LP的错误信息
/****************************************
w 说明：
0~n存储梯度和其他临时信息
n+1~n+m存储scalars rho?
n+m+1~n+2m存储alpha，在计算h*g的时候使用
n+2m+1~n+2m+nm存储最后m次搜索的步骤结果
n+2m+nm+1~n+2m+2nm存储最后的梯度偏差
****************************************/
后面主要是代码可以看懂就不写了，直接找的有人已经翻译了-----------摘自http://qxred.ycool.com/post.1634524.html
DOUBLE PRECISION GTOL,ONE,ZERO,GNORM,DDOT,STP1,FTOL,STPMIN,
. STPMAX,STP,YS,YY,SQ,YR,BETA,XNORM
INTEGER MP,LP,ITER,NFUN,POINT,ISPT,IYPT,MAXFEV,INFO,
. BOUND,NPT,CP,I,NFEV,INMC,IYCN,ISCN
LOGICAL FINISH
C
SAVE
DATA ONE,ZERO/1.0D+0,0.0D+0/
C
C INITIALIZE
C ----------
C
IF(IFLAG.EQ.0) GO TO 10
GO TO (172,100) IFLAG
C IFLAG一开始为0,如果IFLAG=1表示上次MCSEARCH(172)需要计算g(a),g'(a)
C 判断是否满足wolfe conditions,IFLAG=2表示需要用户提供H,在CRF中,不需要用
C 户提供H
10 ITER= 0
IF(N.LE.0.OR.M.LE.0) GO TO 196
IF(GTOL.LE.1.D-04) THEN
C GTOL不能太小,默认0.9
IF(LP.GT.0) WRITE(LP,245)
GTOL=9.D-01
ENDIF
NFUN= 1
C 迭代次数
POINT= 0
FINISH= .FALSE.
IF(DIAGCO) THEN
DO 30 I=1,N
30 IF (DIAG(I).LE.ZERO) GO TO 195
ELSE
C 自动设置H
DO 40 I=1,N
40 DIAG(I)= 1.0D0
ENDIF
C
C THE WORK VECTOR W IS DIVIDED AS FOLLOWS:
C ---------------------------------------
C THE FIRST N LOCATIONS ARE USED TO STORE THE GRADIENT AND
C OTHER TEMPORARY INFORMATION.
C LOCATIONS (N+1)...(N+M) STORE THE SCALARS RHO.
C LOCATIONS (N+M+1)...(N+2M) STORE THE NUMBERS ALPHA USED
C IN THE FORMULA THAT COMPUTES H*G.
C LOCATIONS (N+2M+1)...(N+2M+NM) STORE THE LAST M SEARCH
C STEPS.
C LOCATIONS (N+2M+NM+1)...(N+2M+2NM) STORE THE LAST M
C GRADIENT DIFFERENCES.
C
C THE SEARCH STEPS AND GRADIENT DIFFERENCES ARE STORED IN A
C CIRCULAR ORDER CONTROLLED BY THE PARAMETER POINT.
C
C W=[f'_k,RHO,ALPHA,S,Y]
C 其中保证第一次计算得到的值放在第一个，后面的依次循环存放
C 比如在RHO_1在第二次迭代中计算得到，那存放RHO_i的地址为(k-2)%m+1
C 而y_1在第二次迭代中计算得到，那么存放Y_i的地址为IYPT+(k-2)%m+1
ISPT= N+2*M
IYPT= ISPT+N*M
DO 50 I=1,N
50 W(ISPT+I)= -G(I)*DIAG(I)
C q=-H*g(0)
GNORM= DSQRT(DDOT(N,G,1,G,1))
C |g(0)|
STP1= ONE/GNORM
C STP1=1/|g(0)|初始化第一次试探步长
C PARAMETERS FOR LINE SEARCH ROUTINE
C
FTOL= 1.0D-4
MAXFEV= 20
C
IF(IPRINT(1).GE.0) CALL LB1(IPRINT,ITER,NFUN,
* GNORM,N,M,X,F,G,STP,FINISH)
C
C --------------------
C MAIN ITERATION LOOP
C --------------------
C
80 ITER= ITER+1
C k=ITER
INFO=0
BOUND=ITER-1
IF(ITER.EQ.1) GO TO 165
C 第一次试探，用取a_l=0,a_t=STP1进入MCSEARCH搜索
IF (ITER .GT. M)BOUND=M
C BOUND=min(M,k-1)表示回退的次数，比如k=2时，只需要回退1次
YS= DDOT(N,W(IYPT+NPT+1),1,W(ISPT+NPT+1),1)
C NPT=(k-1),y_{k-1}*s_{k-1}
IF(.NOT.DIAGCO) THEN
C 无需用户给出H,lbfgs自动估计H,CRF中DIAGCO=FALSE
YY= DDOT(N,W(IYPT+NPT+1),1,W(IYPT+NPT+1),1)
DO 90 I=1,N
90 DIAG(I)= YS/YY
ELSE
IFLAG=2
RETURN
ENDIF
100 CONTINUE
IF(DIAGCO) THEN
DO 110 I=1,N
110 IF (DIAG(I).LE.ZERO) GO TO 195
ENDIF
C
C COMPUTE -H*G USING THE FORMULA GIVEN IN: Nocedal, J. 1980,
C "Updating quasi-Newton matrices with limited storage",
C Mathematics of Computation, Vol.24, No.151, pp. 773-782.
C ---------------------------------------------------------
C
CP= POINT
C CP=(k-2)%m+1,这里k>=2,CP从1开始，代表了ALPHA_{k-1},RHO_{k-1}
C 循环存放的序号,比如k=2时，第一次计算RHO,RHO_1的序号为CP=1
C 而(k-1)%m代表了s_{k},y_{k}的存放序号，比如s_1的偏移为0
IF (POINT.EQ.0) CP=M
C 如果k=m+1,那么CP=m
W(N+CP)= ONE/YS
C RHO_{k-1}=1/[y_{k-1}*s_{k-1}]，计算并保存RHO_{k-1}
DO 112 I=1,N
112 W(I)= -G(I)
C q=-f'_k
CP= POINT
DO 125 I= 1,BOUND
CP=CP-1
IF (CP.EQ. -1)CP=M-1
C CP=(k-1-I)%m
C CP代表了s_{k-I},y_{k-I}的存放序号
C s_{k-I}从W(ISPT+CP*N+1)存放到W(ISPT+CP*N+N)
SQ= DDOT(N,W(ISPT+CP*N+1),1,W,1)
C s_{k-I}*y_{k-I}
INMC=N+M+CP+1
IYCN=IYPT+CP*N
W(INMC)= W(N+CP+1)*SQ
C ALPHA_{k-I}=RHO_{k-I}*s_{k-I}*y_{k-I}
CALL DAXPY(N,-W(INMC),W(IYCN+1),1,W,1)
C q=q-ALPHA_{k-I}*y_{k-I}
125 CONTINUE
C CP=POINT-BOUND=max{k-m,0}
DO 130 I=1,N
130 W(I)=DIAG(I)*W(I)
C r=H0*q
DO 145 I=1,BOUND
C FOR i=k-m,...,k-1
YR= DDOT(N,W(IYPT+CP*N+1),1,W,1)
BETA= W(N+CP+1)*YR
C BETA=RHO_{k-M-1+I}*y_{k-M-1+I}*r
INMC=N+M+CP+1
BETA= W(INMC)-BETA
ISCN=ISPT+CP*N
CALL DAXPY(N,BETA,W(ISCN+1),1,W,1)
C r=r+s_{k-M-1+I}*BETA
CP=CP+1
C CP+1是ALPHA,RHO的偏移，CP是y,s的偏移
IF (CP.EQ.M)CP=0
145 CONTINUE
C
C STORE THE NEW SEARCH DIRECTION
C ------------------------------
C r=-Hk*f'_k
DO 160 I=1,N
160 W(ISPT+POINT*N+I)= W(I)
C y_k=r
C OBTAIN THE ONE-DIMENSIONAL MINIMIZER OF THE FUNCTION
C BY USING THE LINE SEARCH ROUTINE MCSRCH
C ----------------------------------------------------
165 NFEV=0
STP=ONE
IF (ITER.EQ.1) STP=STP1
DO 170 I=1,N
170 W(I)=G(I)
172 CONTINUE
CALL MCSRCH(N,X,F,G,W(ISPT+POINT*N+1),STP,FTOL,
* XTOL,MAXFEV,INFO,NFEV,DIAG)
IF (INFO .EQ. -1) THEN
C 需要g'(a),g(a)
IFLAG=1
RETURN
ENDIF
IF (INFO .NE. 1) GO TO 190
NFUN= NFUN + NFEV
C
C COMPUTE THE NEW STEP AND GRADIENT CHANGE
C -----------------------------------------
C
NPT=POINT*N
DO 175 I=1,N
W(ISPT+NPT+I)= STP*W(ISPT+NPT+I)
C s_k
175 W(IYPT+NPT+I)= G(I)-W(I)
C y_k
POINT=POINT+1
IF (POINT.EQ.M)POINT=0
C
C TERMINATION TEST
C ----------------
C
GNORM= DSQRT(DDOT(N,G,1,G,1))
XNORM= DSQRT(DDOT(N,X,1,X,1))
XNORM= DMAX1(1.0D0,XNORM)
IF (GNORM/XNORM .LE. EPS) FINISH=.TRUE.
C
IF(IPRINT(1).GE.0) CALL LB1(IPRINT,ITER,NFUN,
* GNORM,N,M,X,F,G,STP,FINISH)
IF (FINISH) THEN
IFLAG=0
RETURN
ENDIF
GO TO 80
C
C ------------------------------------------------------------
C END OF MAIN ITERATION LOOP. ERROR EXITS.
C ------------------------------------------------------------
C
190 IFLAG=-1
IF(LP.GT.0) WRITE(LP,200) INFO
RETURN
195 IFLAG=-2
IF(LP.GT.0) WRITE(LP,235) I
RETURN
196 IFLAG= -3
IF(LP.GT.0) WRITE(LP,240)
C
C FORMATS
C -------
C
200 FORMAT(/' IFLAG= -1 ',/' LINE SEARCH FAILED. SEE'
. ' DOCUMENTATION OF ROUTINE MCSRCH',/' ERROR RETURN'
. ' OF LINE SEARCH: INFO= ',I2,/
. ' POSSIBLE CAUSES: FUNCTION OR GRADIENT ARE INCORRECT',/,
. ' OR INCORRECT TOLERANCES')
235 FORMAT(/' IFLAG= -2',/' THE',I5,'-TH DIAGONAL ELEMENT OF THE',/,
. ' INVERSE HESSIAN APPROXIMATION IS NOT POSITIVE')
240 FORMAT(/' IFLAG= -3',/' IMPROPER INPUT PARAMETERS (N OR M',
. ' ARE NOT POSITIVE)')
245 FORMAT(/' GTOL IS LESS THAN OR EQUAL TO 1.D-04',
. / ' IT HAS BEEN RESET TO 9.D-01')
RETURN
END
C
C LAST LINE OF SUBROUTINE LBFGS
C
C
2、sdrive.f包含一个简单的测试问题，函数最优解为X=(1,1...1)，f=0
下面一边写出源代码一边说明
program sdrive //这里ndim即为X的维数
parameter(ndim = 2000, msave = 7, nwork = ndim*(2*msave+1)+2*msave)
double x[ndim], g[ndim], diag[ndim], w[nwork];
double f, eps, xtol, gtol, t1, t2, stpmin, stpmax;
int iprint[2], iflag, icall, n, m, mp, lp, j;
bool diagco;
external lb2
common /lb3/mp, lp, gtol, stpmin, stpmax;
n = 100;
m = 5;
iprint[0] = 1;
iprint[1] = 0;
diagco = false;
eps = 1.0D-5;
xtol = 1.0D-16;
icall = 0;
iflage = 0;
for(j=1; j<=n; j+=2){ x[j] = -1.2D0; x[j+1] = 1.D0;}
20:continue;
f = 0.D0;
for(j=1; j<=n; j+=2){
t1 = 1.D0-x[j];
t2 = 1.D1*(x[j+1]-x[j]*2);
g[j+1] = 2.D1*t2;
g[j] = -2.D0*(x[j]*g[j+1]+t1);
f=f+t1*2+t2*2;
}
call lbfgs(n, m, x, f, g, diag, iprint, eps, xtol, w, iflag);
if(flag<=0) go to 50
icall = icall+1;
if(icall>2000) go to 50;
go to 20;
50:continue;
end
http://www.alglib.net/optimization/lbfgs.php
在这个网址能够下到alglib库，里面包含很多方法，有多个版本的，其中有LBFGS的方法
http://www.alglib.net/translator/man/manual.cpp.html#unit_minlbfgs
这个页面有对minlbfgs subpackage的详细介绍以及给了的两个例子
minlbfgs subpackage包括两个类：minlbfgsreport minlbfgsstate
函数有：minlbfgscreate minlbfgsoptimize minlbfgsrestartfrom minlbfgsresults
minlbfgsresultsbuf minlbfgssetcholeskypreconditioner minlbfgssetcond
minlbfgssetdefaultpreconditioner minlbfgssetstpmax minlbfgssetxrep
函数说明：
minlbgfscreate:
用户使用：
1、初始化算法的状态调用函数MinLBFGSCreate()
2、调整参数调用函数MinLBFGSSetCond() MinLBFGSSetStpMax()及其他参数
3、调用函数MinLBFGSOptimize()计算F/G
4、调用MinLBFGSResults()得到结果
5、选择性调用：MinLBFGSRestartFrom()解决另一个问题，问题同样为N/M，但是以另一个指针和功能函数开始
输入参数：
N-问题规模
M-BFGS算法中要求的海塞矩阵更新数，3<=M<=7 M<=N
x-初始化结果估计值 N维
输出参数：
state-存储算法状态
注意：
1、通过MinLBFGSSetCond()函数调整停止条件
2、如果目标函数包括exp()或者其他快增长的功能函数，优化算法可能导致溢出，使用MinLBFGSSetStpMax()函数限制算法的步长，然而，L-BFGS算法很少需要这样的调整
void alglib::minlbfgscreate()(
ae_int_t m,
real_1d_array x,
minlbfgsstate &state);
void alglib::minlbfgscreate()(
ae_int_t n,
ae_int_t m,
real_1d_array x,
minlbfgsstate &state);
minlbfgsoptimize function:
state-算法状态
grad-调用计算函数计算在点x的函数和梯度值
rep-选择性的调用，每次迭代后调回，可以设置为NULL
ptr-选择性使用的指针，可以为NULL
void minlbfgsoptimize(
minlbfgstate &state,
void (*grad)(const real_1d_array &x, double &func, real_1d_array &grad, void *ptr),
void *ptr = NULL);
minlbfgsrestartfrom function:
state-同上
x-新的开始指针
void alglib::minlbfgsrestartfrom(minlbfgsstate state, real_1d_array x);
minlbfgsresults function:
input:
state-算法状态
output:
x-结果，0~(n-1)
rep-优化反馈--- 结束类型:
*-2有最优解 *-1制定了不正确的参数 *1相关函数的改进只是epsF
*2相关步长只是EpsX *4梯度规范化只是EpsG *5 使用MaxIts步长
*7停止约束太过严格，还有提升的空间
--- 迭代次数
---函数值
void alglib::minlbfgsresults(
minlbfgsstate state,
real_1d_array &x,
minlbfgsreport &rep);
minlbfgsresultsbuf function:
用于为X[]提前申请缓冲区存储，如果buffer的大小过小能够重新设置大小。函数主要是为了用于解决循环算法内部的因为循环导致向量的重定位太过庞大
void alglib::minlbfgsresultsbuf(
minlbfgsstate state,
real_1d_array &x,
minlbfgsreport &rep);
minlbfgssetcholeskypreconditioner function:(设置乔姆斯基运行条件)
条件的修正：海塞矩阵的估计需要用到乔姆斯基的因式分解
input:
state-同上
p-三角的调节器，海塞矩阵的估计的乔姆斯基因式分解向量[0..N-1, 0..N-1]
IsUpper-高级或者低级P是否给定
void alglib::minlbfgssetcholeskypreconditioner(
minlbfgsstate state,
real_2d_array p,
bool isupper);
minlbfgssetcond function:
设置L-BFGS优化算法的停止条件
input:
state-同上
EpsG->=0 如果满足条件||G||<EpsG 则子程序终止返回 G为梯度
EpsF->=0 如果满足条件|F(k+1)-F(k)|<=EpsF*max{|F(k)|,|F(k+1)|,1}
EpsX->=0 如果满足条件在k+1次迭代后，|X(k+1)-X(k)|<=EpsX则子程序终止返回
MaxIts-最大迭代次数，如果为0，则不限制迭代次数
如果EpsG=0, EpsF=0, EpsX=0 且 MaxIts=0 则使用默认设置
void alglib::minlbfgssetcond(
minlbfgsstate state,
double epsg,
double epsf,
double epsx,
ae_int_t maxits);
minlbfgssetdefaultpreconditioner function:
条件调整函数，使用默认条件设置
input:
state-同上
void alglib::minlbfgssetdefaultpreconditioner(minlbfgsstate state);
minlbfgssetstpmax function:
设置最大步长
input:
state-同上
StpMax-最大步长，>=0，设置0.0即为默认不限制步长
当函数包括exp等快速增长的函数时需要设置防止溢出
void alglib::minlbfgssetstpmax(minlbfgsstate state, double stpmax);
minlbfgssetxrep function:
函数调整输出report
input:
state-同上
NeedXRep-report是否需要，true则调用rep()
void alglib::minlbfgssetxrep(minlbfgsstate state, bool needxrep);
相关链接：http://blog.csdn.net/settingsun1225/article/details/6142739
http://blog.csdn.net/settingsun1225/article/details/6142853
http://qxred.ycool.com/post.1634524.html