單形法－oooooooooolichen1943的部落格

線性規劃 (四)：單形法

Posted on 09/27/2013 by ccjou

本文的閱讀等級：中級

美國數學家丹齊格 (George Dantzig) 於1947年提出第一個有效的線性規劃算法──單形法 (simplex method，或稱單純形法)，被後人譽為線性規劃之父。1946年丹齊格從加州大學柏克萊分校取得博士學位，次年他以數學顧問身分為美國空軍工作，經常他被要求解決一些與規劃相關的問題，譬如，如何配置預算、人力、飛機，及其他資源使達到最佳的成本效益。因為這些問題多少與經濟學有關，丹齊格跑去徵詢經濟學家科普斯曼 (Tjalling Koopmans，1975年諾貝爾經濟學獎得主)。丹齊格設想經濟學家應當早已發展出解線性規劃問題的技術，出乎意料之外，科普斯曼告訴他經濟學家也沒有線性規劃問題的系統化解法。於是在1947年的暑假，丹齊格決定自己著手尋找一個方法^[1]。

George Dantzig (1914-2005) From http://www.uh.edu/engines/georgedantzig.jpg

丹齊格在嘗試找出求解算法的過程中很快地發現一個重要的幾何性質：可行域 (滿足約束條件的可行解所形成的集合) 是一個 (有界或無界) 多胞形。直白地說，多胞形是一個如鑽石形狀的凸集，其邊界都是平的 (見“多胞形”)。更關鍵的是，如果一線性規劃問題存在一有限數值的最佳解，則最佳解必定是多胞形的一個端點。這是否意味我們必須算出每一個端點的目標函數值才能找到最佳解？丹齊格辯稱無此必要，他想出一個簡單的算法 (見下圖)：因為每一個端點的目標函數值通常都不相同，我們可以從任何一個端點開始，沿著多胞形的邊移動至可使目標函數下降的鄰接端點 (線性規劃標準型設定為最小化問題)。丹齊格稱這個算法為單形法。起初，丹齊格以為單形法的搜尋程序可能極度低效，不過這回他錯了。除了極少數刻意製造的特例或退化問題，在絕大多數的現實應用中，從起始端點至最佳解的移動次數與線性規劃問題的限制條件數有相同的數量級。

單形法幾何解釋

我們將單形法的演算程序整理於下：

找出可行域的一個端點作為起始端點，計算此端點的目標函數值。
檢查與目前端點鄰接的所有端點的目標函數值是否小於目前端點的目標函數值。如果是，挑選一個最具改善潛力的鄰接端點，沿著可行域的邊移動至該端點。
重複步驟2直到不存在可使目標函數下降的鄰接端點，終止端點即為最佳解。

讀者或許懷疑：單形法找到的是一個局部 (local) 最佳解 (或稱相對最佳解)，我們怎麼能夠斷言它就是全域 (global) 最佳解？線性規劃問題是一個典型的凸優化問題，意思是可行域為一凸集，目標函數為一凸函數。凸優化問題蘊含一個美好的性質：任何一個局部最佳解即為全域最佳解 (見“凸優化──凸函數的最小化”，定理一)。吃下定心丸後，現在我們可以放心開始推演單形法的理論細節。

單形法的幾何概念雖然簡單明瞭，但對應的代數運算卻不易理解。繼續往下閱讀之前，建議讀者先充實背景知識“線性規劃 (三)：最佳基可行解定理”。單形法的演算程序包含四關鍵問題：(1) 如何找出一個可行解作為起始端點？(2) 給定一端點，如何選擇使目標函數下降的鄰接端點？(3) 如何確定目前的端點是最佳解？(4) 選定一鄰接端點後，如何從目前端點移動至該端點？下面依序解釋。

考慮線性規劃標準型問題

$\begin{array}{lll} \min& z=\mathbf{c}^T\mathbf{x} \\ \hbox{subject to}& A\mathbf{x}=\mathbf{b}, ~~\mathbf{x}\ge\mathbf{0}, \end{array}$

其中 $A$ 是一 $m\times n$ 階矩陣， $m\le n$ ， $\mathbf{b}$ 是一 $m$ 維向量， $\mathbf{c}$ 是一 $n$ 維向量， $\mathbf{x}$ 是一 $n$ 維未知向量 (見“線性規劃 (一)：標準型問題”)。單形法要求 $\hbox{rank}A=m$ ，也就是說，約束條件的係數矩陣 $A$ 有線性獨立的列 (row)。令 $\mathbf{a}_j$ 是 $A=[a_{ij}]$ 的第 $j$ 行 (column)， $j=1,\ldots,n$ 。我們的目標是找出 $x_1\ge 0, x_2\ge 0,\ldots,x_n\ge 0$ 使最小化 $z=c_1x_1+\cdots+c_nx_n$ 並滿足

$\displaystyle \mathbf{a}_1x_1+\mathbf{a}_2x_2+\cdots+\mathbf{a}_nx_n=\mathbf{b}$ ，

其中

$\displaystyle \mathbf{a}_j=\begin{bmatrix} a_{1j}\\ a_{2j}\\ \vdots\\ a_{mj} \end{bmatrix},~~~\mathbf{b}=\begin{bmatrix} b_1\\ b_2\\ \vdots\\ b_m \end{bmatrix}$ ，

且 $a_{ij}, b_i, c_j$ 都是常數 (實數)。

(1) 如何找出一個可行解作為起始端點？

假設我們將未知變數 $x_1,x_2,\ldots,x_n$ 重新命名排序使得 $\mathbf{a}_1,\mathbf{a}_2,\ldots,\mathbf{a}_m$ 組成一線性獨立集，並令

$\displaystyle B=\begin{bmatrix} \mathbf{a}_1&\mathbf{a}_2&\cdots&\mathbf{a}_m \end{bmatrix}$ 。

明顯地， $B$ 是一 $m\times m$ 階可逆矩陣。在 $A\mathbf{x}=\mathbf{b}$ 左乘 $B^{-1}$ ，如下：

$\displaystyle B^{-1}A\mathbf{x}=(B^{-1}\mathbf{a}_1)x_1+(B^{-1}\mathbf{a}_2)x_2+\cdots+(B^{-1}\mathbf{a}_n)x_n=B^{-1}\mathbf{b}$ ，

或表示為

$\displaystyle \bar{\mathbf{a}}_1x_1+\bar{\mathbf{a}}_2x_2+\cdots+\bar{\mathbf{a}}_nx_n=\bar{\mathbf{b}}$ ，

其中 $\bar{\mathbf{a}}_j=B^{-1}\mathbf{a}_j$ 和 $\bar{\mathbf{b}}=B^{-1}\mathbf{b}$ 分別是 $\mathbf{a}_j$ 和 $\mathbf{b}$ 參考基底 $\{\mathbf{a}_1,\mathbf{a}_2,\ldots,\mathbf{a}_m\}$ 的座標向量。我們稱上式為約束條件的典型形式 (canonical form)。因為 $B^{-1}B=I$ ，可知 $\bar{\mathbf{a}}_i=B^{-1}\mathbf{a}_i=\mathbf{e}_i$ ， $1\le i\le m$ ，其中 $\mathbf{e}_i$ 是標準單位向量 (第 $i$ 元等於 $1$ ，其餘元等於 $0$ )。所以典型形式可簡化表示為

$\displaystyle\begin{array}{cccccccccccccl} x_1& & & & +&\bar{a}_{1,m+1}x_{m+1}& +&\bar{a}_{1,m+2}x_{m+2} & +&\cdots& +&\bar{a}_{1n}x_n&=&\bar{b}_1\\ &x_2& & &+&\bar{a}_{2,m+1}x_{m+1}& +&\bar{a}_{2,m+2}x_{m+2} & +&\cdots& +&\bar{a}_{2n}x_n&=&\bar{b}_2\\ & &\ddots& &&&&\vdots &\\ & & & x_m &+&\bar{a}_{m,m+1}x_{m+1}& +&\bar{a}_{m,m+2}x_{m+2} & +&\cdots& +&\bar{a}_{mn}x_n&=&\bar{b}_m. \end{array}$

上式中， $x_1,x_2,\ldots,x_m$ 稱為基變數 (basic variable)， $x_{m+1},x_{m+2},\ldots,x_n$ 稱為非基變數。從計算面來看，典型形式其實就是增廣矩陣 $\begin{bmatrix} A~\vert~\mathbf{b} \end{bmatrix}$ 經過基本列運算得到的簡約列梯形式 (reduced row echelon form，見“高斯─約當法”)。基變數對應線性方程的軸變數 (pivot variable)，非基變數對應自由變數 (free variable)。設非基變數為零，可得一解，稱為基解 (basic solution)：

$\displaystyle x_1=\bar{b}_1,x_2=\bar{b}_2,\ldots,x_m=\bar{b}_m,x_{m+1}=0,x_{m+2}=0,\ldots,x_n=0$ 。

如果每一 $\bar{b}_i\ge 0$ ，所有的基變數都不為負值 (非基變數等於零)，則此基解是一個可行解，故稱為基可行解 (basic feasible solution)。可行域的端點就是基可行解 (見“線性規劃 (二)：端點與基解”，基可行解─端點定理)，我們解決了第一個問題。

(2) 給定一端點，如何選擇使目標函數下降的鄰接端點？

為簡化問題，假設所有的基解都是非退化 (nondegenerate) 基解，意思是基解的全部基變數值都不等於零 (這個漏洞有補救辦法，在此從略)。假定我們已經得到了一個非退化基可行解 $(x_1,\ldots,x_m,x_{m+1},\ldots,x_n)$ ，其中 $x_j=\bar{b}_j>0$ ， $1\le j\le m$ ，且 $x_j=0$ ， $m+1\le j\le n$ 。令這個基可行解的目標函數值為

$\displaystyle z_0=c_1x_1+c_2x_2+\cdots+c_mx_m=\mathbf{c}_B^T\bar{\mathbf{b}}$ ，

其中 $\displaystyle \mathbf{c}_B^T=\begin{bmatrix} c_1&c_2&\cdots&c_m \end{bmatrix}^T$ 。單形法的思路是選擇一個非基變數 $x_k$ ， $m+1\le k\le n$ ，讓它進入基變數集，使新產生的基可行解的目標函數值小於原本解的目標函數值。為了達到這個目的，我們運用一些代數技巧設法將 $z_0$ 引進目標函數 $z$ 。寫出典型形式的向量表達式

$\displaystyle \mathbf{e}_1x_1+\mathbf{e}_2x_2+\cdots+\mathbf{e}_mx_m=\bar{\mathbf{b}}-\bar{\mathbf{a}}_{m+1}x_{m+1}-\bar{\mathbf{a}}_{m+2}x_{m+2}-\cdots-\bar{\mathbf{a}}_nx_n$ ，

左乘 $\mathbf{c}_B^T$ ，可得

$\displaystyle c_1x_1+c_2x_2+\cdots+c_mx_m=\mathbf{c}_B^T\bar{\mathbf{b}}-\sum_{j=m+1}^nz_jx_j$ ，

其中 $z_j=\mathbf{c}_B^T\bar{\mathbf{a}}_j$ 。上式等號兩邊同時加上 $\sum_{j=m+1}^nc_jx_j$ ，

$\displaystyle z=z_0+\sum_{j=m+1}^n(c_j-z_j)x_j$ 。

注意，上式可以用來計算任一可行解的目標函數值。根據這個目標函數表達式，下面的定理提供了一個選擇有效的非基變數 $x_k$ 的方法。

基可行解改善定理：給定一非退化基可行解 (端點)，對應的目標函數值為 $z_0$ 。如果存在 $m+1\le k\le n$ 使得 $c_{k}-z_k<0$ ，則存在一基解其目標函數值為 $z<z_0$ 。若非基變數 $x_k$ 納入基變數集後可得一基可行解，此基可行解即為鄰接端點。

證明於下。考慮基可行解 $(x_1,x_2,\ldots,x_m,0,\ldots,0)$ ，目標函數值等於 $z_0$ 。在不失一般性的原則下，假定 $c_{m+1}-z_{m+1} < 0$ 。將新可行解 $(x'_1,x'_2,\ldots,x'_m,x'_{m+1},0,\ldots,0)$ ， $x'_{m+1}>0$ ，代入前面導出的目標函數表達式，立得

$\displaystyle z-z_0=(c_{m+1}-z_{m+1})x'_{m+1}<0$ ，

即不論新可行解為何，必定有 $z<z_0$ 。明顯地， $x'_{m+1}$ 的數值應當越大越好，也就是沿著可行域的邊盡可能地移動使目標函數越小越好，直到觸碰一個新基可行解 (端點)，證畢。

如果存在不只一個 $k$ 使得 $c_k-z_k<0$ ，我們可以選擇

$\displaystyle k=\arg\min_{m+1\le j\le n}c_j-z_j$ 。

不過這僅表示沿著對應非基變數 $x_k$ 變動的該邊移動具有最陡的目標函數下降率，最後獲得的目標函數減少量仍由 $(c_k-z_k)x_k$ 決定，下面 (4) 將探討這個問題。

(3) 如何確定目前的端點是最佳解？

根據基可行解改善定理，如果每一 $c_j-z_j\ge 0$ ， $m+1\le j\le n$ ，則目前的基可行解即為最佳解。因為任何其他的可行解受限於 $x_j\ge 0$ ， $1\le j\le n$ ，由上面的目標函數表達式可知目前的目標函數值滿足 $z-z_0\ge 0$ 。

(4) 選定一鄰接端點後，如何從目前端點移動至該端點？

換一個問法，選定一個非基變數 $x_k$ ， $m+1\le k\le n$ ，如何得到與目前的基可行解鄰接的新基可行解 (如果存在的話)？約束條件式可以僅用基變數表示如下：

$\displaystyle x_1\mathbf{a}_1+x_2\mathbf{a}_2+\cdots+x_m\mathbf{a}_m=\mathbf{b}$ 。

如果要把非基變數 $x_k$ 納入基變數集，有一個基變數必須轉為非基變數 (因為總共有 $m$ 個基變數)。對應 $x_k$ 的係數向量 $\mathbf{a}_k$ 可唯一表示為 $\mathbf{a}_1,\mathbf{a}_2,\ldots,\mathbf{a}_m$ 的線性組合

$\displaystyle \mathbf{a}_k=\alpha_{1k}\mathbf{a}_1+\alpha_{2k}\mathbf{a}_2+\cdots+\alpha_{mk}\mathbf{a}_m$ 。

將上式通乘 $\epsilon\ge 0$ ，並與前述約束條件式相減，

$\displaystyle (x_1-\epsilon\alpha_{1k})\mathbf{a}_1+(x_2-\epsilon\alpha_{2k})\mathbf{a}_2+\cdots+(x_m- \epsilon\alpha_{mk})\mathbf{a}_m+\epsilon\mathbf{a}_k=\mathbf{b}$ 。

對於任一 $\epsilon\ge 0$ ，上式給出由 $m+1$ 個向量組合的約束式，除了 $x_k$ ，其他非基變數仍固定為零。若 $\epsilon=0$ ，此即原本的基可行解。當 $\epsilon$ 由零開始增大， $\mathbf{a}_k$ 的係數隨之增加，微小的 $\epsilon$ 仍產生可行解，但已不再是基解，因為當 $\epsilon$ 增大時， $\mathbf{a}_1,\ldots,\mathbf{a}_m$ 的係數也將以線性方式隨之增大或減小。如果存在至少一 $\alpha_{ik}>0$ ，設

$\displaystyle \epsilon=\min_{1\le i\le m}\left\{\frac{x_i}{\alpha_{ik}}~\vline~\alpha_{ik}>0\right\}$ 。

根據非退化假設， $x_i>0$ ， $i=1,2,\ldots,m$ ，可知 $\epsilon>0$ 。若 $p$ 是使上式最小化的指標，則 $x_p-\epsilon\alpha_{pk}=0$ 。這麼一來，原本的基變數 $x_p$ 便被非基變數 $x_k$ 取代。如果每一 $\alpha_{ik}\le 0$ ，也就是說，只要 $\epsilon$ 增大， $x_i-\epsilon\alpha_{ik}$ 也跟著增大，這意味新基可行解 (端點) 不存在。若 $x_k$ 加入基變數集後無法得到一基可行解，表示可行域是一無界集，我們大可沿著可行域的邊盡情地移動可行解，也就是說，目標函數可以任意減小 (趨於 $-\infty$ )。

從1947年以來，不論問題尺度如何增大，電腦技術如何演變，單形法是最常被採用的線性規劃算法，直到近年出現了內點法 (interior point method)，單形法獨霸的局面才受到挑戰。丹齊格在線性規劃的開創性貢獻赫赫有名，不過多數人可能從未聽聞青年丹齊格的一件真實事件。1939年丹齊格是加州大學柏克萊分校的博士生。有一天丹齊格上課遲到，鐘響不久內曼 (Jerzy Neyman) 教授在黑板上寫了兩個著名的未解統計學問題，丹齊格進教室後以為這兩個問題是指定作業便將它們抄下。據丹齊格說，那些問題「看來比平常的難一點」，幾天後他遞交了兩個問題的完整解答，還誤認自己遲交作業。過了六個星期，內曼教授興奮地探訪丹齊格，告訴他其中一題的解答即將送交一份數學期刊發表。1986年丹齊格於接受訪問時透露：「一年後，我正在為論文題目發愁，內曼只是聳聳肩跟我說，把這兩個問題裝訂好，他願意接受它們作為我的論文。」多年以後，另一位研究者瓦爾德 (Abraham Wald) 解出第二個問題並準備發表論文，當他獲知丹齊格之前的解答，便大方地將丹齊格列為合著者。丹齊格的事蹟很快地散播開來，並作為啟發教材展示正面思考的力量。隨著時間過去，故事主角丹齊格的名字已被人遺忘，實際發生的細節也被更改，但基本故事仍以都會傳奇的形式流傳^[2]。

註解：
[1] John L. Casti，Five Golden Rules，1996年出版，頁189。
[2] 維基百科：George Dantzig