關于最優化求解，吳軍有篇blog講的很不錯，證明和解釋都很清楚，具體可以參考http://www.cnblogs.com/joneswood/archive/2012/03/11/2390529.html。
這裡根據那篇blog的内容，主要講解運用最廣泛的LBFGS的算法思想和LBFGS源碼的求解實際的最優化問題。
理論部分
一般優化算法中，比較簡單的是梯度下降法，其主要思想為：
給定目标函數f(x),給定初始值x,沿着目标函數梯度方向下降，尋找最優解。
通過将目标函數，在初始值位置按照梯度下降方向，進行泰勒展開：
：代表第k個點的自變量（向量）；
:機關方向（向量），即|d|=1；
：步長（實數）；
：目标函數在Xk這一點的梯度（導數向量）；
：α的高階無窮小。
式[1]中的高階無窮小可以忽略，是以，要使[1]式取得最小值，應使
取到最小，由此可得，
取
時，目标函數下降得最快，這就是負梯度方向作為“最速下降”方向的由來。
由于梯度下降法隻用到了目标函數的一階導數資訊，由此想到将二階導數資訊運用到優化求解中，這就是牛頓法的思想。在
點處對目标函數進行泰勒展開，并隻取二階導數及其之前的幾項（忽略更高階的導數項），可得：
：目标函數在
這一點的梯度；
：目标函數在
這一點的Hessian矩陣（二階導數矩陣）。
由于極小值點必然是駐點，而駐點是一階導數為0的點，是以，對 r(X) 這個函數來說，要取到極小值，應該分析其一階導數。對X求一階導數，并令其等于0：
是以，下一點的計算方法：
在方向d上按步長1（1×d = d）移動到點X，方向d的計算方法：
牛頓法的基本步驟：每一步疊代過程中，通過解線性方程組得到搜尋方向，然後将自變量移動到下一個點，然後再計算是否符合收斂條件，不符合的話就一直按這個政策（解方程組→得到搜尋方向→移動點→檢驗收斂條件）繼續下去。由于牛頓法利用了目标函數的二階導數資訊，其收斂速度為二階收斂，是以它比最速下降法要快，但是其對一般問題不是整體收斂的，隻有當初始點充分接近極小點時，才有很好的收斂性。
牛頓法中，在每一次要得到新的搜尋方向的時候，都需要計算Hessian矩陣（二階導數矩陣）。在自變量維數非常大的時候，這個計算工作是非常耗時的，是以，拟牛頓法的誕生的意義就是：它采用了一定的方法來構造與Hessian矩陣相似的正定矩陣，而這個構造方法計算量比牛頓法小。
a. DFP算法
假設目标函數可以用二次函數進行近似（實際上很多函數可以用二次函數很好地近似）：
忽略高階無窮小部分，隻看前面3項，其中A為目标函數的Hessian矩陣。此式等号兩邊對X求導，可得：
于是，當 X=Xi 時，将[2]式兩邊均左乘Ai+1-1，有：
上式左右兩邊近似相等，但如果我們把它換成等号，并且用另一個矩陣H來代替上式中的A-1，則得到：
方程[4]就是拟牛頓方程，其中的矩陣H，就是Hessian矩陣的逆矩陣的一個近似矩陣.在疊代過程中生成的矩陣序列H0，H1，H2，……中，每一個矩陣Hi+1，都是由前一個矩陣Hi修正得到的。DFP算法的修正方法如下，設：
再設：
其中，m和n均為實數，v和w均為N維向量。将[6]代入[5]式，再将[5]式代入[4]式，可得：
得到的m，n，v，w值如下：
将[8]~[11]代入[6]式，然後再将[6]代入[5]式，就得到了Hessian矩陣的逆矩陣的近似陣H的計算方法：
通過上述描述，DFP算法的流程大緻如下：
已知初始正定矩陣H0，從一個初始點開始（疊代），用式子
來計算出下一個搜尋方向，并在該方向上求出可使目标函數極小化的步長α，然後用這個步長，将目前點挪到下一個點上，并檢測是否達到了程式終止的條件，如果沒有達到，則用上面所說的[13]式的方法計算出下一個修正矩陣H，并計算下一個搜尋方向……周而複始，直到達到程式終止條件。
值得注意的是，矩陣H正定是使目标函數值下降的條件，是以，它保持正定性很重要。可以證明，矩陣H保持正定的充分必要條件是：
在疊代過程中，這個條件也是容易滿足的。
b. BFGS算法
BFGS算法和DFP算法有些相似，但是形式上更加複雜一些。BFGS算法目前仍然被認為是最好的拟牛頓算法。BFGS算法中矩陣H的計算公式如下所示：
在[14]式中的最後一項（藍色部分）就是BFGS比DFP多出來的部分，其中w是一個n×１的向量。
在目标函數為二次型時，無論是DFP還是BFGS——也就是說，無論[14]式中有沒有最後一項——它們均可以使矩陣H在n步之内收斂于A-1。
代碼實驗：
網絡上面有很多BFGS的實作算法，用的比較多的是c++實作的源碼庫libbfgs,http://www.chokkan.org/software/liblbfgs/.
網上編譯即可以使用。
本文以求解 f(x1,x2) = 100(x2-x1^2)^2 + (1-x1)^2為例，對比matlab和libbfgs的結果。
在matlab中求解該優化問題，可得如下結果：
運用libbfgs庫，程式設計實作優化問題：
#include <stdio.h>
#include <lbfgs.h>
//優化回調函數，每次疊代運作，在這裡設定目标函數和梯度方程
static lbfgsfloatval_t evaluate(
void *instance,
const lbfgsfloatval_t *x,
lbfgsfloatval_t *g,
const int n,
const lbfgsfloatval_t step
)
{
lbfgsfloatval_t fx = 0.0;
//計算目标函數
fx = 100 * (x[1] - x[0] * x[0]) * (x[1] - x[0] * x[0]) + (1 - x[0]) * (1 - x[0]);
//計算梯度方程
g[0] = -(400 * x[0] * (x[1] - x[0] * x[0]) + 2 * (1 - x[0]));
g[1] = 200 * (x[1] - x[0] * x[0]);
return fx;
}
//疊代回調函數，每次疊代可以用來輸出疊代求解獲得的值
static int progress(
void *instance,
const lbfgsfloatval_t *x,
const lbfgsfloatval_t *g,
const lbfgsfloatval_t fx,
const lbfgsfloatval_t xnorm,
const lbfgsfloatval_t gnorm,
const lbfgsfloatval_t step,
int n,
int k,
int ls
)
{
//輸出目前疊代次數，目标值，x1,x2的值以及目标函數和梯度的誤差二範數等
printf("Iteration %d:\n", k);
printf(" fx = %f, x[0] = %f, x[1] = %f\n", fx, x[0], x[1]);
printf(" xnorm = %f, gnorm = %f, step = %f\n", xnorm, gnorm, step);
printf("\n");
return 0;
}
//定義目前參數的個數
#define N 2
int main(int argc, char *argv[])
{
lbfgsfloatval_t fx; //目标函數
lbfgsfloatval_t *x = lbfgs_malloc(N); //配置設定變量空間
lbfgs_parameter_t param; //參數值
if (x == NULL) {
printf("ERROR: Failed to allocate a memory block for variables.\n");
return 1;
}
/* 初始化參數. */
x[0] = -1.2;
x[1] = 2;
/* 初始化 L-BFGS 優化算法的設定參數. */
lbfgs_parameter_init(¶m);
/*
優化求解
*/
ret = lbfgs(N, x, &fx, evaluate, progress, NULL, ¶m);
/* 輸出最終的優化結果. */
printf("L-BFGS optimization terminated with status code = %d\n", ret);
printf(" fx = %f, x[0] = %f, x[1] = %f\n", fx, x[0], x[1]);
lbfgs_free(x); //釋放變量空間
return 0;
}
運作程式：
對比可得，L-BFGS程式經過36次疊代即可求得最優解，且跟matlab的求解結果一緻。