【學界/編碼】凸優化演算法 I: 內點法(interior point method)求解線性規劃問題

04-14

作者： @李恩志 Louisiana State University 計算物理學博士。
『運籌OR帷幄』責任編輯： @文雨之（東北大學系統工程博士生）
本篇文章是由以上作者在知乎上的優秀文章（原文鏈接:凸優化演算法 I: 內點法(interior point method)求解線性規劃問題），通過『運籌OR帷幄』責任編輯整理修改而成的。
歡迎原鏈接轉發，付費轉載請前往 @留德華叫獸 的主頁獲取信息，盜版必究。
敬請關注和擴散本專欄及同名公眾號，會邀請全球知名學者陸續發布運籌學、人工智慧中優化理論等相關乾貨、知乎Live及行業動態：

『運籌OR帷幄』大數據人工智慧時代的運籌學

0前言

內點法是一種處理帶約束優化問題的方法，其在線性規劃，二次規劃，非線性規劃等問題上都有著很好的表現。在線性規劃的問題上，相對於鼎鼎大名的單純形法，內點法是多項式演算法，而單純形法並非多項式演算法。從實際應用的效果來說，內點法也達到了足以和單純形法分庭抗衡的地步，尤其針對大規模的線性規劃問題內點法有著更大的發展潛力。

線性規劃單純形法：通過一系列迭代達到最優解，迭代點沿著可行多面體的邊界從一個頂點到另一個頂點，直到得到最優解。一般而言單純形法每次迭代的開銷相對內點法來說較小，但所需迭代次數較多。
線性規劃內點法：同樣是通過一系列迭代達到最優解，但其是從多面體內部逐漸收斂到最優解。一般而言內點法每次迭代的開銷相對單純形法來說較大，但所需迭代次數較少。

內點法並不僅僅用於線性規劃的求解，值得一提的是內點法的很多思想有著更廣泛的應用，例如障礙函數法的思想。

1導論

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

$ext{min } c^{T}x \ ext{s.t. } Ax preceq b$

這裡，目標函數為線性函數 $c^{T}x: mathbb{R}^{n} ightarrow mathbb{R}$ ，約束條件為 $A_{ij}x_{j} le b_{i}, i = 1, 2, ..., m.$ 矩陣 $A$ 為 $m imes n$ 的滿秩矩陣，其中 $m$ 為約束條件的個數， $n$ 為變數的個數。通常約束條件的個數大於變數的個數，所以有 $m > n$ .

2線性規劃問題的等價(近似)表述

這個線性規劃問題可以重新表述為計算 $ext{min } f(x)$ ，其中

$f(x) = c^{T}x + sum_{i = 1}^{m}I(A_{ij}x_{j} - b_{i})$

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

$[ I(u) = egin{cases} 0 & ext{if } x le 0 \ infty & ext{if } x > 0 end{cases} ]$

引入這個函數的意義在於可以將約束條件直接寫入到目標函數裡面，這樣我們直接求新的函數的極小值就可以了，而不必藉助於未知乘子。但是這裡有一個問題，那就是indicator函數存在不可求導的點，因此在求函數極小值的時候我們沒法通過普通的微分法來確定函數的極小值。為了規避這個問題，我們可以用一個光滑的函數來近似這個indicator函數。一個不錯的選擇是用 $I_{t}(u) = -frac{1}{t}log (-u)$ 來代替indicator函數。 $I_{t}(u)$ 只有在 $u < 0$ 的時候有定義，我們規定當 $u > 0$ 的時候 $I_{t}(u) = infty$ . 而且參數 $t > 0$ 越大，函數 $I_{t}(u)$ 就越接近於 $I(u)$ . 所以我們可以通過調節 $t$ 的值來調節這個函數的近似程度。使用這個近似的indicator函數，我們新的的目標函數可以寫作

$f(x) =t c^{T}x - sum_{i = 1}^{m} log(-A_{ij} x_j + b_i)$ .

因為線性函數是凸函數，並且 $I_{t}(u)$ 也是凸函數，所以 $f(x)$ 是凸函數，因此我們可以很容易用凸優化的經典方法得到該函數的極小值。

3計算函數的梯度和Hessian矩陣

為了求函數的極小值，根據微積分的經典結果，我們只需令函數的梯度等於零，然後計算梯度為零時對應的解 $x^{star}$ .

函數的梯度為

$frac{partial f}{partial x_{k}} = t c_{k} + sum_{i = 1}^{m}frac{-A_{ik}}{A_{ij}x_{j} - b_{i}}$

Hessian 矩陣為

$frac{partial^2 f}{partial x_{k}partial x_{l}} = sum_{i = 1}^{m}frac{A_{ik} A_{il}}{Big(A_{ij} x_{j} - b_{i} Big)^2}$

定義對角型矩陣為

$D_{ij} = delta_{ij} frac{1}{Big(A_{ik}x_{k} - b_{i}Big)^2}$

於是Hessian矩陣可以寫作

$H_{f} = A^{T} D A$

因為 $D$ 為正定矩陣，所以Hessian矩陣至少為半正定矩陣。所以函數 $f(x)$ 是一個凸函數。而且矩陣 $D$ 為可逆矩陣，矩陣 $A$ 滿秩，所以Hessian矩陣為可逆矩陣。於是函數 $f(x)$ 為強凸函數。所以，要計算 $abla f = 0$ 的根，我們可以用高效的牛頓迭代法。

4牛頓迭代法

我們現在的目標是計算 $abla f = 0$ 的根。因為Hessian矩陣可逆，所以我們可以用牛頓迭代法求解。牛頓迭代法為

$x^{(n+1)} = x^{(n)} - H^{-1} abla f$

在這個迭代過程中，參數 $t$ 為固定的。每對應一個 $t$ ，我們都可以得到一個解 $x_{t}^{star}$ . 如果我們掃描參數 $t$ ，我們就可以得到一系列的解。其中最大的 $t$ 的對應的解應該最精確。

一個更好的演算法是選取一個比較小的初始參數 $t_{0}$ ，求出這個參數對應的解 $x_{t_0}^{star}$ . 然後增加 $t$ ，用之前得到的解 $x_{t_0}^{star}$ 來初始化當前 $t$ 所對應的牛頓迭代法的試解。這樣算出來的解應該比直接計算 $t$ 所對應的解更加精確。這樣逐步迭代從小 $t_{0}$ 到大 $t$ 可以得到一系列的解，最後得到的解 $x_{t}^{star}$ 稱作是淬火解。這個解應該可以滿足我們的精度需求。

5一個例子

為了展示該演算法的功效，我們舉一個簡單的例子。這個例子可以用簡單的幾何方法求解出來。我們要展示的是，我們可以用數值方法得到同樣的結果（在一定的精度範圍之內）。之所以要用數值方法求解這個簡單的問題，是因為數值方法對更複雜的問題同樣有效，而簡單的幾何方法對複雜問題卻已經不適用了。

現在要研究的例子為

$ext{min } x + y \ ext{s.t. }\ x + 2y le 1 \ 2x + y le 1 \ x ge 0 \ y ge 0$

也就是

$c = egin{pmatrix} 1 \ 1 end{pmatrix}$ , $A = egin{pmatrix} 1 & 2 \ 2 & 1 \ -1 & 0 \ 0 & -1 end{pmatrix}$ , $b = egin{pmatrix} 1 \ 1 \ 0 \ 0 end{pmatrix}$ .

梯度為

$abla f(x) = tc - A^{T} egin{pmatrix} frac{1}{x + 2y -1}\ frac{1}{2x + y -1} \ frac{1}{-x} \ frac{1}{-y} end{pmatrix}$

Hessian矩陣為

$H_{f}(x) = A^{T}DA$

其中，對角型矩陣 $D$ 為

$D = egin{pmatrix} frac{1}{(x + 2y - 1)^2 } & 0 & 0 & 0 \ 0 & frac{1}{(2x + y - 1)^2 } & 0 & 0 \ 0& 0 & frac{1}{x^2} & 0 \ 0 & 0 & 0 & frac{1}{y^2} end{pmatrix}$

對於一個固定的參數 $t$ ，選擇一個恰當的初始解 $x^{(0)}$ ，代入牛頓迭代公式

$x^{(n+1)} = x^{(n)} - H_{f}(x^{(n)})^{-1} abla f(x^{(n)}), n ge 0$

可以得到一個依賴於參數 $t$ 的解 $x^{star}_{t}$ .

用幾何方法很容易求得這個例子的解為 $(x^{star}, y^{star} ) = (0, 0)$ . 所以我們期待，當 $lim_{t ightarrow infty} (x_{t}^{star}, y_{t}^{star}) = (0, 0)$ .

我需要用程序求這個例子的數值解。我寫了一個Python程序實現這個演算法，程序地址為PrimerLi/linear-programming.

程序結果如下：

這裡用了四個不同的初始點來初始化程序，最終得到解都收斂到了 $(0, 0)$ . 箭頭方向為參數 $t$ 增加時 $(x_t^{star}, y_t^{star} )$ 移動的方向。這正是我們期待的結果。

在此感謝『運籌OR帷幄』審稿人對本文提出了寶貴的意見。

『運籌OR帷幄』審稿人 @覃含章，美國麻省理工學院（MIT）計算科學與工程方向博士在讀，清華大學工業工程及數學與應用數學（第二學位）本科。研究興趣主要為優化理論，機器學習演算法在運營管理中的應用。

以上『運籌OR帷幄』專欄所有文章都會同步發送至留德華叫獸的頭條主頁，以及同名微信公眾號，目前預計受眾10w +

如果你是運籌學/人工智慧碩博或在讀，請在下圖的公眾號後台留言：「加微信群」。系統會自動辨認你的關鍵字，並提示您進一步的加群要求和步驟，邀請您進全球運籌或AI學者群（群內學界、業界大佬雲集）。

運籌學|控制論|優化理論愛好者，歡迎加qq群：686387574

人工智慧愛好者，歡迎加qq群： 685839321

數據科學|分析愛好者，歡迎加qq群：130414574

最後敬請關注和擴散本專欄及同名公眾號，會陸續發布運籌學、人工智慧中優化理論相關乾貨及行業動態：『運籌OR帷幄』大數據和人工智慧時代下的運籌學 - 知乎專欄