這一節課講解了線性規划中的原始對偶方法(primal-dual method),並以最短路問題為例說明該方法的應用。
原始對偶方法
原始對偶方法利用的就是上一節課中講到的互補松弛定理。我們首先找到對偶問題的一個可行解 $y$,並嘗試找到一個原問題的可行解 $x$,使得 $x$ 和 $y$ 滿足互補松弛定理。如果我們找到了這樣的 $x$,那么 $x$ 和 $y$ 就分別是原問題和對偶問題的最優解;否則我們就需要調整 $y$,讓它變得更好,繼續嘗試,直到找到最優解為止。
尋找對偶可行解
考慮以下線性規划問題作為原問題 $$\begin{matrix} \min\limits_x & c^Tx \\ \text{s.t.} & Ax = b \ge 0 \\ & x \ge 0 \end{matrix}$$ 它的對偶問題為 $$\begin{matrix} \max\limits_x & b^Ty \\ \text{s.t.} & A^Ty \le c \end{matrix}$$ 我們首先需要獲得一個對偶問題的可行解。如果 $c \ge 0$,顯然 $y = 0$ 就是對偶問題的一個可行解,否則我們可以使用如下方法尋找可行解。
給原問題增加一個變量與一條約束 $x_1 + x_2 + \dots + x_n + x_{n+1} = b_{m+1}$,其中 $b_{m+1}$ 需要足夠大,至少要大等於 $x_1 + x_2 + \dots + x_n$ 可能的最大值。這樣的約束就不會改變原問題。
加入新約束后,我們寫出新的對偶問題 $$\begin{matrix} \max\limits_y & b^Ty + b_{m+1}y_{m+1} \\ \text{s.t.} & A^Ty + y_{m+1}\begin{bmatrix} 1 & 1 & \dots & 1 \end{bmatrix}^T \le c \\ & y_{m+1} \le 0 \end{matrix}$$ 我們很容易看出這個新對偶問題的可行解為 $y_i = 0 \quad \forall i \in \{1, 2, \dots, m\}$ 以及 $y_{m+1} = \min\limits_i \quad c_i$。
限制的原問題 (RP) 與限制的對偶問題 (DRP)
假設我們找到了對偶問題的可行解 $y$,現在我們需要尋找一個原問題的可行解 $x$ 滿足互補松弛定理。
設 $A_j$ 表示矩陣 $A$ 的第 $j$ 列,定義 $J = \{ j | A_j^Ty = c_j \}$(稱 $J$ 為允許指標集,簡單來說 $J$ 就是以 $y$ 作為對偶可行解時,對偶問題中較緊的那些限制的編號)。根據原問題的定義和互補松弛定理,我們有 $$Ax = b \\ x_j = 0 \quad \forall j \not\in J \\ x_j \ge 0 \quad \forall j \in J$$ 如果我們能找到一個 $x$ 滿足上面三個條件,$x$ 和 $y$ 就能滿足互補松弛定理。
為了獲得一個可行的 $x$,我們使用一個類似於兩階段法的思想(雖然兩階段法和原始對偶方法可能沒啥聯系- -),構造一個優化問題,稱為限制的原問題(restricted primal,RP)$$\begin{matrix} \max\limits_x & \sum\limits_{i=1}^m \bar{x}_i \\ \text{s.t.} & \sum\limits_{j\in J} a_{i,j}x_j + \bar{x}_i = b_i \\ & x_j \ge 0, \bar{x}_i \ge 0 \end{matrix}$$ 很顯然,若目標函數值取 0,那么我們就找到了滿足互補松弛定理的 $x$。
我們再來寫出 RP 的對偶問題,稱為限制的對偶問題(DRP)$$\begin{matrix} \max\limits_x & b^Ty \\ \text{s.t.} & A^T_jy \le 0 & \forall j \in J \\ & y_i \le 1 & \forall i \in \{1, 2, \dots, m\} \end{matrix}$$ 設 $\bar{y}$ 為 DRP 的最優解,根據強對偶定理,若 $b^T\bar{y} = 0$ 則 RP 的目標函數值也可以取到 0,那么當前對偶可行的 $y$ 就是對偶問題的最優解;否則我們有 $b^T\bar{y} \ge 0$,因為至少 $\bar{y} = 0$ 是 DRP 的可行解。
一開始要求原問題的 $b \ge 0$ 可能是為了保證 RP 問題有可行解(如果 $b \ge 0$ 那么 RP 問題至少有可行解 $x = 0$ 和 $\bar{x} = b$)。不過這可能也不是很有必要,如果發現 RP 問題無可行解,或者 DRP 問題無有限最優解,那就說明了原問題無可行解。
改進當前的對偶可行解
如果 DRP 的最優解讓 DRP 的目標函數值超過 0,說明當前的 $y$ 還不是最優的。我們來想辦法用 $\bar{y}$ 改進 $y$。我們讓新的對偶可行解 $\hat{y} = y + \theta\bar{y}$(其中 $\theta \ge 0$),顯然有 $b^T\hat{y} = b^Ty + \theta b^T\bar{y} > b^Ty$,說明對偶問題的目標函數值的確改進了。但別忘了,$\hat{y}$ 仍然需要是對偶可行的,也就是說,$A^T\hat{y} = A^Ty + \theta A^T\hat{y} \le c$ 仍要滿足。
我們來考慮 $\theta$ 的取值范圍。對於 $j \in J$,因為 $A^T_j\bar{y} \le 0$,無論 $\theta$ 取多大,都不會超過 $c$ 的限制,所以這些項不會限制 $\theta$ 的取值;對於 $j \not\in J$,我們選擇 $$\theta = \min_{j \not\in J, A^T_j\bar{y} > 0} \quad \frac{c_j - A^T_jy}{A^T_j\bar{y}}$$ 這樣就會讓 $j \not\in J$ 中的一條限制變緊,其它限制則不會超過,仍然是對偶可行的。當然啦,如果 $\theta$ 可以無限增大,說明對偶問題沒有有限最優解,那么原問題無可行解。
將 $y$ 調整為 $\hat{y}$ 之后,就進入下一輪迭代繼續調整(當然我們不需要再從原問題開始,慢慢寫出 RP 和 DRP 了,直接確定新的 $J$ 和新的 DRP 即可),直到 DRP 的最優解讓目標函數值為 0,此時的 $y$ 就是對偶問題的最優解。
算法的時間復雜度
事實上,原始對偶問題就是通過枚舉 $J$ 來獲得對偶問題的最優解。由於 $J$ 至多有 $2^n$ 個,所以這個方法是一定會終止的。當然,這個方法的最差復雜度也是指數級的。如果限制條件在進入 $J$ 后不會出來,那么算法就能在線性步數內結束。
這個方法看起來比單純形法麻煩多了。單純形法只要解一個優化問題就能得到最優解,這個方法一下子變出來四個優化問題。但是這個方法對特定問題是非常有效的,對於一些特定問題來說,我們可以一下子看出(或者在很短的時間內方便地算出)DRP 的最優解。下面就來舉例說明原始對偶方法的應用。
應用:最短路問題
用線性規划表示最短路問題
我們來考慮有向圖上的最短路問題,起點為 $s$,終點為 $t$。對有向圖定義點 - 弧關聯矩陣,這個矩陣中每列對應一條邊,每行對應一個頂點。若一條邊是一個頂點的出邊,那么矩陣對應元素為 1;若一條邊是一個頂點的入邊,那么矩陣對應元素為 -1。設 $w_i$ 表示第 $i$ 條邊的長度,$x_i$ 表示第 $i$ 條邊是否在最短路上,那么最短路問題就是解如下線性規划問題 $$\begin{matrix} \min\limits_x & w^Tx \\ \text{s.t.} & Ax = v \\ & x \ge 0 \end{matrix}$$ 其中 $$v_i = \begin{cases} 1 & i = s \\ -1 & i = t \\ 0 & \text{otherwise} \end{cases}$$ 雖然這個問題是一個線性規划問題,但是最優解中每個變量的值不是 0 就是 1(最短路問題中不會有半條邊在最短路上吧- -)。最優解中每個變量的取值范圍在 0~1 之間還好理解(如果有一個變量大於 1,就必須有變量取負值,這樣不滿足 $x \ge 0$ 的條件),為什么每個變量的取值還是整數呢?
我們使用單純形法求解時,得到的解是 $x_B = A_B^{-1}b$。如果 $A_B^{-1}$ 和 $b$ 的元素都是整數,得到的解自然是整數。我們引入全幺模矩陣(total unimodular matrix)的定義:一個整數矩陣 $A$ 稱為全幺模矩陣,如果 $A$ 的任何一個子方陣的行列式取值為 0,1 或 -1。如果我們能證明有向圖的點 - 弧關聯矩陣是全幺模矩陣,那么 $A_B$ 的行列式取值就為 0,1 或 -1。由於 $A_B^{-1} = A_B^* / |A_B|$,而由於 $A_B$ 的元素均為整數,伴隨矩陣 $A_B^*$ 的元素也為整數,那么 $A_B^{-1}$ 的元素自然也都是整數了。
以下對方陣的邊長 $n$ 使用數學歸納法,證明任意一個有向圖的點 - 弧關聯矩陣的子方陣行列式取值為 0,1 或 -1。
起始步驟:對於任意有向圖 $n = 1$ 的子方陣,根據點 - 弧關聯矩陣的定義,子方陣的行列式取值為 0,1 或 -1。
推遞步驟:假設對於任意有向圖 $n \le k$ 的子方陣,都有行列式取值為 0,1 或 -1。
對於任意有向圖 $n = k + 1$ 的子方陣,根據點 - 弧關聯矩陣的定義,每一列至多有兩個非 0 元素,且這些元素中至多一個為 1,至多一個為 -1。
如果該子方陣存在一列沒有非 0 元素,那么該子方陣的行列式取值為 0;
如果該子方陣存在一列只有一個非 0 元素,由於該元素為 1 或 -1,那么該子方陣行列式的絕對值等於該元素余子式的絕對值。將原有向圖去掉該元素對應的點和邊后,這個余子陣可以看作是新有向圖的子方陣. 根據歸納法假設,該元素余子式的絕對值為 0 或 1,那么該子方陣的行列式取值為 0,1 或 -1;
如果該子方陣的每一列都有兩個非 0 元素,那么根據點 - 弧關聯矩陣的定義,對該子方陣的每一行求和,將會得到一個元素都是 0 的行向量,所以該子方陣的行列式取值為 0。
綜上所述,對於 $n = k + 1$ 的子方陣,其行列式的取值仍為 0,1 或 -1。
根據上述數學歸納法,有向圖的點 - 弧關聯矩陣是全幺模矩陣。
最短路問題的對偶問題和 DRP
根據點 - 弧關聯矩陣的定義很容易發現,把 $A$ 的每一行加起來,最后會獲得都是 0 的一行,也就是說 $A$ 中的行向量是線性相關的。那么不妨從 $A$ 中去掉代表終點 $t$ 的那一行,得到新矩陣 $\bar{A}$。這樣,最短路問題就變為 $$\begin{matrix} \min\limits_x & w^Tx \\ \text{s.t.} & \bar{A}x = v \\ & x \ge 0 \end{matrix}$$ 其中 $$v_i = \begin{cases} 1 & i = s \\ 0 & \text{otherwise} \end{cases}$$ 我們寫出這個問題的對偶問題 $$\begin{matrix} \max\limits_{\pi} & \pi_s \\ \text{s.t.} & \pi_i - \pi_j \le w_e & \forall e = (i, j) \in E \\ & \pi_t = 0 \end{matrix}$$ 其中 $(i, j)$ 代表一條從 $i$ 連到 $j$ 的邊,$E$ 是有向圖的邊集。這個對偶問題的可行解很容易得到,只要令$\pi = 0$ 即可。另外,雖然 $\bar{A}$ 矩陣比 $A$ 矩陣少了代表 $t$ 的那一行,但是只要我們定義 $\pi_t = 0$,就又能把對偶問題寫成上面的形式了。這個形式就是差分約束系統,以前做算法競賽的時候一直不能理解為什么差分約束求最大值反而要跑一個最短路,現在終於明白了,原來差分約束是最短路問題的對偶問題...
我們再來寫出最短路問題的 DRP $$\begin{matrix} \max\limits_{\pi} & \pi_s \\ \text{s.t.} & \pi_i - \pi_j \le 0 & \forall (i, j) \in J \\ & \pi_i \le 1 \\ & \pi_t = 0 \end{matrix}$$ 這個 DRP 問題的解是非常容易看出的。如果從點 $s$ 到某個點 $p$ 的最短路已經算出來了,而邊 $j$ 就在這條最短路上,根據最短路的三角不等式,這條邊一定在 $j$ 中。我們把所有已知最短路的點構成的連通塊稱為 $C$,為了讓 $\pi_s$ 的取值最大並滿足 DRP 的限制條件,$C$ 中所有點的變量值都必須和 $\pi_s$ 相同。顯然,若 $t \not\in C$,那么 $\pi_s = 1$;若 $t \in C$,那么由於 $\pi_t = 0$ 的條件,$\pi_s = 0$,說明我們找到了從 $s$ 到 $t$ 的最短路。
我們用一張有向圖來模擬一下最短路的原始對偶算法。
$$\begin{array}{c|cccccc} \text{iter}. & \pi_s & \pi_{v2} & \pi_{v3} & \pi_{v4} & \pi_{v5} & \pi_t \\ \hline 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 & 0 & 0 & 0 \\ 2 & 2 & 0 & 1 & 0 & 0 & 0 \\ 3 & 4 & 2 & 3 & 0 & 2 & 0 \\ 4 & 6 & 5 & 4 & 2 & 4 & 0 \end{array}$$ 其實這就是我們熟悉的 Dijkstra 算法的過程。