運籌優化（五）--線性規劃之內點法

Eason.wxd發表於2019-01-12
優化

近年來的內點演算法主要有三大類：

（1）投影尺度法,它是Karmarkar演算法的原型。這個方法要求問題具有特殊的單純形結構和最優目標值為零,在實際計算過程中, 需要用複雜的變換將實際問題轉換為這種標準形式。因此, 投影尺度法在實際中應用較少。

（2）仿射尺度法，這是一種已經比較成熟和廣泛的演算法。目前原仿射尺度法和對偶仿射尺度法雖然應用較多,但這兩種方法的多項式時間複雜性還不能從理論上得到證明。

（3）路徑跟蹤法,又稱為跟蹤中心軌跡法。這種方法將對數壁壘函式與牛頓法結合起來應用到線性規劃問題, 並且已經從理論上證明具有多項式時間複雜性。路徑跟蹤法收斂迅速, 魯棒性強, 對初值的選擇不敏感,現在已經被推廣應用到二次規劃領域,正被進一步發展為從複雜性角度研究一般非線性規劃的內點演算法,是目前最有發展潛力的一類內點演算法。

內點法是求解不等式約束最優化問題的一種十分有效方法，但不能處理等式約束。因為構造的內點懲罰函式是定義在可行域內的函式，而等式約束優化問題不存在可行域空間，因此，內點法不能用來求解等式約束優化問題。

投影尺度法

核心原理：

針對線性規劃的標準型：

min(max)\, \, c \cdot x

s.t.\: \: \: \: \: \: Ax = b \, \, \, \, \, \, \, x \geq 0

改進最快的移動方向是目標函式向量\Delta x = c,方向d在條件A\Delta x = 0上的投影，使得線性約束Ax = b成立，得到：

\Delta x = Pd

其中P是投影矩陣：P = (I - A^{^{T}}(AA^{T})^{-1}A)

仿射尺度變換

仿射尺度變換採用最簡單的策略遠離邊界，通過對模型進行重複尺度變換使當前解的變形與所有不等式約束距離相等。

仿射尺度變換之後，採用上面投影矩陣，計算改進方向。移動步長是向量d的長度倒數。

演算法不斷重複尺度變換，在投影方向上移動改進，計算新的可行解，直到目標函式值變化差異穩定。

障礙函式法（Barrier Method）

參考：https://zhuanlan.zhihu.com/p/32685234

障礙函式思想的內點法主要思想是：在可行域的邊界築起一道很高的“圍牆”，當迭代點靠近邊界時，目標函式徒然增大，以示懲罰，阻止迭代點穿越邊界，這樣就可以將最優解“檔”在可行域之內了。

約束優化演算法的基本思想是：通過引入效用函式的方法將約束優化問題轉換成無約束問題，再利用優化迭代過程不斷地更新效用函式，以使得演算法收斂。 這句話，可以用下面的推導來說明：

對於下面的不等式約束的優化問題： 
\min f(x), x \in R^n

s.t \quad g_{i}(x) \leq0, i=1,2,...,m

利用內點法進行求解時，構造懲罰函式的一般表示式為 
\varphi (X, r)=f(X)-r\sum_{i=1}^{m}\frac{1}{g_{i}(X)}

或者 
\varphi (X, r)=f(X)-r\sum_{i=1}^{m}{\ln[-g_{i}(X)]}

也就是通過這一步，將約束（\quad g_{i}(x)）優化問題轉化為無約束（\varphi (X, r)）優化問題，然後就是一個迭代逼近的，求極值的優化問題了。

更一般的表述：

線性規劃問題的一般形式為：

這個線性規劃問題可以重新表述為計算：

這裡，我們使用了一個indicator函式，定義為：

引入這個函式的意義在於可以將約束條件直接寫入到目標函式裡面，這樣我們直接求新的函式的極小值就可以了，而不必藉助於未知乘子。 但是這裡有一個問題，那就是indicator函式存在不可求導的點，因此在求函式極小值的時候我們沒法通過普通的微分法來確定函式的極小值。為了規避這個問題，我們可以用一個光滑的函式來近似這個indicator函式。一個不錯的選擇是用 ：

來代替indicator函式。 上式只有在u < 0的時候有定義，我們規定當u > 0的時候上式為正無窮，而且引數t > 0越大，函式Itu就越接近於Iu. 所以我們可以通過調節t的值來調節這個函式的近似程度。使用這個近似的indicator函式，我們新的的目標函式可以寫作 :

 因為線性函式是凸函式，並且indicator也是凸函式，所以fx是凸函式，因此我們可以很容易用凸優化的經典方法得到該函式的極小值。
演算法步驟，
取初始懲罰因子r^{(0)}>0，允許誤差\epsilon>0；
在可行域D內選取初始點X(0)，令k=1；
構造懲罰函式φ(X,r(k))，從X(k−1)點出發用無約束優化方法求懲罰函式φ(X,r(k))的極值點(X∗,r(k))；
檢查迭代終止準則：如果滿足
\|X^{*} r^{(k)}-X^{*} r^{(k-1)}\|\leq\epsilon_{1}=10^{-5}-10^{-7}
或者
\|\frac{\varphi (X^{*} ,r^{(k)})-\varphi (X^{*}, r^{(k-1)})}{\varphi (X^{*}, r^{(k-1)})}\|\leq\epsilon_{2}=10^{-3}-10^{-4}
則停止迭代計算，並以(X∗,r(k))作為原目標函式f(X)的約束最優解，否則轉入下一步；
取r(k+1)=cr(k)，X(0)=X∗r(k)，k=k+1，轉向步驟3。遞減係數c=0.1−0.5，通常取0.1。

內點懲罰函式法特點及其應用


懲罰函式定義於可行域內，序列迭代點在可行域內不斷趨於約束邊界上的最優點（這就是稱為內點法的原因，只適合求解具有不等式約束的優化問題
內點法求解案例
用內點法求下面約束優化問題的最優解，取迭代初始X0=[0,0]T，懲罰因子的初始值r0=1，收斂終止條件\|X^k - X^{k-1}\| \leq \varepsilon
\varepsilon = 0.01
\min f(X) = x_1^2 + x_1^2 - x_1x_2 - 10x_1 - 4x_2 + 60
構造內懲罰函式：\varphi(X, r) = x_1^2 + x_1^2 - x_1x_2 - 10x_1 - 4x_2 + 60 -r\ln(x_1 + x_2 -8)
用解析法求內懲罰函式的極小值
\nabla\varphi(X, r) = [2x_1 - x_2 - 10 - \frac{r}{x_1 + x_2 - 8} \quad 2x_2 - x_1 - 4 - \frac{r}{x_1 + x_2 - 8}]
令∇φ(X,r)=0得：\begin{align}2x_1 - x_2 - 10 - \frac{r}{x_1 + x_2 - 8} = 0 \\ 2x_2 - x_1 - 4 - \frac{r}{x_1 + x_2 - 8} = 0\end{align}
解得：

X^*_1(r) = \begin{bmatrix}\frac{13 + \sqrt{9+2r}}{2} & \frac{9 + \sqrt{9+2r}}{2}\end{bmatrix}^{\mathrm{T}}
∵g(X∗1(r))>0
∴ 捨去X∗1(r)
∵φ(X,r)為凸函式

∴無約束優化問題的最優解為X^*(r) = X^*_2(r) = \begin{bmatrix}\frac{13 - \sqrt{9+2r}}{2} & \frac{9 - \sqrt{9+2r}}{2}\end{bmatrix}^{\mathrm{T}}
求最優解
當r0=1時，X^*(r^0) = \begin{pmatrix}4.8417 & 2.8417\end{pmatrix}^{\mathrm{T}},\|X^*(r^0) - X^0\| = 5.6140 > \varepsilon
當r1=0.1時，X^*(r^1) = \begin{pmatrix}4.9834 & 2.9834\end{pmatrix}^{\mathrm{T}},\|X^*(r^1) - X^*(r^0)\| = 0.2004 > \varepsilon
當r2=0.01時，X^*(r^2) = \begin{pmatrix}4.9983 & 2.9983\end{pmatrix}^{\mathrm{T}},\|X^*(r^2) - X^*(r^1)\| = 0.0211 > \varepsilon
當r3=0.01時，X^*(r^3) = \begin{pmatrix}4.9998 & 2.9998\end{pmatrix}^{\mathrm{T}}，\|X^*(r^3) - X^*(r^2)\| = 0.0021 < \varepsilon
即X∗(r3)為最優解

原始對偶內點法（Primal Dual Interior Point Method）

下面一堆公式，不理解的化只要記住：

原始對偶內點法始終保持原始問題和對偶問題的解在每次迭代中嚴格可行，並在搜尋過程中系統性的減少互補鬆弛性的違背程度。其難點在於最優條件中的互補鬆弛約束，除非當前解最優，否則會存在對偶間隙，迭代的目的在於不斷減少當前的互補鬆弛結果與目標差異。
在 Primal Dual 的方法中，我們直接考慮一個標準的 LP 命題。
\min c^Tx, \quad \text{subject to } Ax=b, x\ge0 \tag9
它的對偶命題為：
\max b^T\lambda, \quad \text{subject to } A^T\lambda +s=c,s\ge0 \tag{10}
上面兩個命題的KKT條件為：
\begin{align} A^T\lambda+s&=c, \tag{11a} \\ Ax&=b, \tag{11b} \\ x_i s_i &= 0, i=1,...,n\tag{11c} \\(x,s) &\ge0.\tag{11d} \end{align}
為了後面的推導方便，將上面的KKT條件化為矩陣形式如下：
\begin{align} F(x,\lambda,s)=\begin{bmatrix} A^T\lambda+s-c\\ Ax-b\\ XSe\end{bmatrix} &= 0. \tag{12a} \\ (x,s)&\ge 0, \tag{12b}\end{align}
其中，
X=\text{diag} (x_1,x_2,...,x_n),\quad S=\text{diag} (s_1,s_2,...,s_n), \quad e = (1,1,...1)^T
如同一般的優化演算法，這裡需要定義一個量來檢驗當前的迭代點與最優點的差距。在Barrier Method中，使用 duality gap 的上界 \frac mt來檢驗的，在 Primal Dual Method 中，我們定義一個新的 duality measure 來進行某種衡量：
\mu=\frac 1n\sum_{i=1}^nx_is_i=\frac{x^Ts}n
注意：這裡的 μ與 Barrier Method 中的 μ 意義不同，為了與各自書中的表達方式對應，分別選用了各自書中的變數記法。
要求解原始的優化命題，需要做的就是去求解 (12) 這樣一個方程組，由於 F 陣第三行中 XS 一項的存在，因此這是一個非線性系統。要求解這樣一個非線性系統，一種實用的方法就是牛頓法。（注意，這裡說的牛頓法是一種求解非線性方程組的方法，與求解優化命題的牛頓法並不完全一樣，但核心思路是一致的，都是在迭代點處進行二階展開。）在當前點處求解一個前進方向，並通過迭代逼近非線性系統的解。
J(x,\lambda,s) \begin{bmatrix} \Delta x\\ \Delta\lambda \\ \Delta s\end{bmatrix}=-F(x,\lambda,s)
其中 JJ 是FF 陣的雅各比矩陣。我們將 FF 陣中的前兩行分別記為：
r_b=Ax-b,\quad r_c=A^T\lambda+s-c \tag{13}
那麼在每次迭代中需要求解的線性系統為：
\begin{bmatrix} 0 &A^T &I \\ A & 0 & 0 \\ S & 0 &X \end{bmatrix}\begin{bmatrix} \Delta x\\ \Delta\lambda \\ \Delta s\end{bmatrix}=\begin{bmatrix} -r_c \\ -r_b \\ -XSe\end{bmatrix}
因此，在求解得到相應的前進方向後，需要做的就是求解確定一個步長 α 來進行如下的更新
(x,\lambda,s)+\alpha(\Delta x, \Delta\lambda, \Delta s)
其中 α∈(0,1] 的選取要使得原始變數和對偶變數都滿足相應的約束。
看起來這種方法似乎已經可以了，但通過這種被稱為 affine scaling direction 的方向所得到的 α 往往很小。也就是說，很難在迭代中取得進展。一開始看到這個說法時，我也很難理解這句話的意思。所以在這裡我們用一個圖來說明，令 (9) 中的 c = \p\begin{bmatrix} 10 \\ \\ 1\end{bmatrix}，初始點為(0.7,2)，這裡不考慮等式約束，直接採用affine scaling direction 的方向得到的迭代點的軌跡為 


從圖中可以看出，幾乎在從第二次迭代開始，迭代點就一頭撞上了約束。後面的迭代也只能貼著約束的邊緣來走，因為要保迭代點不能違反約束，因此每次的步長 α 都只能取得很小。在這一過程中，一共進行了11次迭代才使得 duality measure μ<10−5。
因此，實際上內點法採用的是一個不那麼 aggressive 的牛頓方向，也就是控制迭代點使其徐徐想約束邊界和最優點處靠近。具體的方法是，我們在用牛頓法求解非線性系統時，在每次迭代中並不要求直接實現 xisi=0，而是令其等於一個逐漸減少的值，具體來說就是 xisi=σμ，其中 μ 是當前的 duality measure，σ∈[0,1]是用於控制下降速度的下降因子。也就是說，在每次迭代中要求解的方程組應該是
\begin{bmatrix} 0 &A^T &I \\ A & 0 & 0 \\ S & 0 &X \end{bmatrix}\begin{bmatrix} \Delta x\\ \Delta\lambda \\ \Delta s\end{bmatrix}=\begin{bmatrix} -r_c \\ -r_b \\ -XSe+\sigma\mu e\end{bmatrix} \tag{14}
其中，σ 被稱為 center parameter 。意在表示我們要通過調整 σ 使得迭代點的軌跡在 central path 附近。
Central Path
在內點法中，central path 是指滿足下面一組方程式的迭代點所組成的所組成的一條弧線
\begin{align}A^T\lambda+s&=c,\tag{15a} \\ Ax&=b,\tag{15b} \\ x_is_i&=\tau, i=1,2...n. \tag{15c} \\ (x,s)&>0 \tag{15d}\end{align}
這與 (11) 中的KKT 條件的區別就在於在第三個條件的等式右端的 0 被 τ>0 替代了。
也就是說，對偶內點法的基本思路就是依據 central path 計算出相應的方向，並控制步長 α 的選擇使得迭代點位於 central path 的某一個鄰域內。

關於步長 α 的選取，central parameter σ 的更新，以及更多的收斂性分析的內容，在這裡不作展開。

舉例
與上一張圖對應，我們同樣取 c=(10 1)，初始點為(0.7,2)，不考慮等式約束。採用對偶內點法的迭代點的軌跡為 


可以看出，引入 central path 後，迭代點能在貼近約束邊界後再次遠離約束邊界（也就是靠近 central path） ，從而保證下一次迭代能取得更大的進步。在這一過程中，一共進行了6次迭代就使得 duality measure μ<10−5。

相關文章