單純性法（Simplex Method）是“線性規劃之父”George Dantzig 最著名的成果，也是求解線性規劃最有力的演算法之一。而這一演算法在求解二次規劃（Quadratic Programming, QP）時的升級版就是有效集法（Active Set Method, ASM）。這兩種演算法的特點都是迭代點會循著約束邊界前進，直到達到問題的最優點。本文對用於求解 QP 命題的 Primal ASM 演算法作以介紹，主要內容的來源是 Nocedal 等人 2006 年的著作 Numerical Optimization 第二版。

有效集（Active Set） 要說有效集法，首先要說說什麼是有效集。有效集是指那些在最優點有效（active）的不等式約束所組成的集合。例如，考慮二次函式  min.fs.t.xy=x2+y2+2x+y≥0≥−1 min.f=x2+y2+2x+ys.t.x≥0y≥−1 函式的等高線及兩個約束的影象如下： 

通過計算或者從圖上可以看出，當沒有約束時，目標函式的最小值在 (−1,−0.5)(−1,−0.5) 處取得。當考慮上述的兩條約束時，目標函式的最小值在 (0,−0.5)(0,−0.5) 處取得。這個時候在最優點處，約束 x≥0x≥0 中的等號被啟用，這條約束就被稱為有效約束（active constraint）。如果我們記兩條約束的編號為 11 和 22，那麼在最優點處的有效集就可以記為  A∗={1} A∗={1} 而如果原命題中 y≥−1y≥−1 這條約束改為 y≥−0.5y≥−0.5 ，即要求解的優化命題變為  min.fs.t.xy=x2+y2+2x+y≥0≥−0.5 min.f=x2+y2+2x+ys.t.x≥0y≥−0.5 函式的等高線及兩個約束的影象如下： 

此時優化命題的解仍然是 (0,−0.5)(0,−0.5) ，但在最優點處兩條約束均被啟用，即此時的有效集可以記為  A∗={1,2} A∗={1,2} 有效集法 從上面可以看出，如果我們能提前知道在最優點處有效的約束，那我們就可以把那些未有效的不等式約束剔除掉並把原命題轉化成更易求解的等式約束命題。然而，在求解之前我們往往對最優點處的有效約束知之甚少。因此如何找到最優點處的有效約束也就是有效集法的主要工作。另外在這裡要提一點，其實在一些應用中，我們需要求解一系列類似的 QP 命題，這個時候我們往往對最優點處的有效約束有一個初始猜測，因此通過這種方式可以實現演算法的熱啟動（Warm Start），從而加速演算法的收斂。下面我們正式開始介紹 ASM 的理論。我們考慮的是帶有不等式約束的凸二次規劃命題：  minxq(x)subject toaTix=12xTGx+xTc≥bi,i∈I minxq(x)=12xTGx+xTcsubject toaiTx≥bi,i∈I 在 ASM 中，我們會構造一個工作集（Working Set），與有效集類似，工作集也是有效約束的集合，但只是我們認為在某次迭代中有效約束的集合，它可能與最優點處的有效集相同，也有可能不同。如果相同，我們可以通過計算對偶變數 λλ 瞭解到此時已經是最優點從而退出迭代。如果不同，我們會對工作集進行更新，從現有工作集中刪除一條約束或者增加一條新的約束到工作集中。因為工作集是隨著迭代而改變的，因此記工作集為 WkWk 。

在第 kk 次迭代開始時，我們首先檢查當前的迭代點 xkxk 是否是當前工作集 WkWk 下的最優點。如果不是，我們就通過求解一個等式約束的 QP 命題來得到一個前進方向 pp ，在計算 pp 的時候，我只關注 WkWk 中的等式約束而忽略原命題的其他約束。為了便於計算，我們定義  p=x−xk,gk=Gxk+c p=x−xk,gk=Gxk+c 把上面的定義代入原命題有：  q(x)=q(xk+p)=12pTGp+gTkp+ρk q(x)=q(xk+p)=12pTGp+gkTp+ρk 其中 ρk=12xTkGxk+xTkcρk=12xkTGxk+xkTc 在這裡可以看作是常數因此可以從上面的優化命題中去掉。因此我們可以寫出如下所示的在第 kk 次迭代需要求解的 QP 子命題：  minqsubject to12pTGp+gTkpaTip=0,i∈Wk minq12pTGp+gkTpsubject toaiTp=0,i∈Wk 我們記上面的命題求解得到的最優解為 pkpk，首先假定求解得到的 pkpk 不為0∗∗。即我們有了一個方向使得沿著這個方向目標函式會下降，然後我們需要確定一個步長 αkαk 來決定走多遠。如果 xk+pkxk+pk 對於所有原命題的約束都是可行的，那麼此時的 αk=1αk=1，否則 αkαk 則是一個小於1的正數。定義步長 αkαk 之後我們就可以得到迭代點更新的式子：  xk+1=xk+αkpk xk+1=xk+αkpk 關於步長 αkαk 的計算，我們有如下的分析。其實，步長的計算主要就是要保證新的迭代點不要違反原命題的約束，那麼在工作集 WkWk 內的約束就不用擔心了，因為在求解 pkpk 的過程中就已經保證了工作集內的約束必須得到滿足。那麼對於非工作集中的約束 i∉Wki∉Wk，我們首先判斷 aTipkaiTpk 的符號，如果 aTipk≥0aiTpk≥0 ，那麼通過分析可知只要步長 αkαk 大於0，則該約束一定可以滿足，因此我們需要關注的主要就是 aTipk<0aiTpk<0 的那些約束。為了保證原命題中的約束 aTi(xk+αkpk)≥biaiT(xk+αkpk)≥bi 能夠滿足，我們需要使得 αkαk 滿足  αk≤bi−aTixkaTipk αk≤bi−aiTxkaiTpk 因此，我們可以把計算 αkαk 的方法顯式地寫出來：  αk=min(1,mini∉Wk,aTipk<0bi−aTixkaTipk) αk=min(1,mini∉Wk,aiTpk<0bi−aiTxkaiTpk) 要注意的是，雖然我們這裡的確把 αkαk 顯式地表示出來，但仔細看就會發現求解 αkαk 其實還是一個逐條約束去檢查的過程。比如對於我們從 αk=1αk=1 開始，逐條約束開始判斷，首先判斷如果 aTipk≥0aiTpk≥0 則直接跳到一條約束，否則可以判斷 αkαk 與 bi−aTixkaTipkbi−aiTxkaiTpk 的大小，如果前者小於後者，那麼沒問題可以去檢查下一條約束，否則直接令 αkαk 取值為後者。以我自己的經驗，當原命題的約束個數較多時，計算步長 αkαk 往往是一個比較耗時的環節。如果能有什麼好辦法可以削減這部分的計算量，那麼對於 QP 的加速求解將會是很有益的工作。

當我們在計算 αkαk 的時候，每條約束都會對應一個不違反其約束的最小的 αkαk 的值，而且在最小的 αkαk 就是這次迭代的步長值 αkαk ，而這條約束就被稱為 blocking constraint（如果 αk=1αk=1 且沒有新的約束在下一步迭代點處啟用，那麼此次迭代沒有 blocking constraint）。注意，也有可能出現 αk=0αk=0 的情況，這是因為在當前迭代點處有其他的有效約束沒有被新增到工作集中。如果 αk<1αk<1， 也就是說下降方向 pkpk 被某條不在工作集 WkWk 內的約束阻攔住了。因此我們可以通過將這條約束新增到工作集的方法來構造新的工作集 Wk+1Wk+1 。

通過上面的方法我們可以持續地向工作集內新增有效約束，直到在某次迭代中我們發現當前的迭代點已經是當期工作集下的最優點。這種情況很容易判斷，因為此時計算出來的 pk=0pk=0。接下來我們就要驗證當前的迭代點是不是原命題的最優點，驗證的方法就是判斷工作集內的約束對應的對偶變數 λλ 是否都大於等於 0，如果滿足這一點，那麼所有的 KKT 條件都成立，可以退出迭代並給出原命題的最優解（關於這部分的推導，詳細內容請見原書的 16.5 節）。λλ 的計算由下式給出  ∑i∈W^aiλi^=g=Gx^+c ∑i∈W^aiλi^=g=Gx^+c 如果通過上式計算出來的有一個或者多個 λ^iλ^i 的值小於 0。那麼就表明通過去掉工作集的某一條或幾條約束，目標函式值可以進一步下降。因此我們會從對應的 λ^iλ^i 值小於 0 的約束中選擇一條，將其從工作集 WkWk 中剔除從而構造出新的工作集 Wk+1Wk+1 。如果有多於一條的可選約束，那麼不同的剔除方法會遵循不同原則，我們在這裡不加更多說明地遵循去除對應 λ^iλ^i 值最小（絕對值最大）的那條約束。

基於上面的討論，我們可以得到一個 Primal 的 ASM 的演算法，具體的演算法可見原書 Algorithm 16.3，這裡我們給出一個演算法的流程圖。

這裡要注意，Primal ASM 演算法的迭代點都是在可行域內或者可行域的邊界上移動，這樣的好處是即使你提前終止迭代，那麼演算法得到的也是一個可行解（Dual 的演算法往往不能保證這一點）。而這樣的缺點則是在演算法啟動是需要我們給定一個可行的初始點，而求解這樣一個可行的初始點往往也不是一件簡單的事情。現在的一種方法是通過一個 Phase I 的過程來得到這樣一個可行的初始點，也就是求解一個只有約束的線性規劃命題。而初始的工作集往往可以取為空集。

舉例 我們用下面這個二維的例子來說明 ASM 演算法的過程。  minxq(x)=(x1−1)2+(x2−2.5)2ssubject tox1−2x2+2−x1−2x2+6−x1+2x2+2x1x2≥0,≥0,≥0,≥0,≥0. minxq(x)=(x1−1)2+(x2−2.5)2ssubject tox1−2x2+2≥0,−x1−2x2+6≥0,−x1+2x2+2≥0,x1≥0,x2≥0. 示意圖如下： 

假設此時迭代初始點為 x0=(2,0)Tx0=(2,0)T，初始的工作集為 W0={3,5}W0={3,5} 。注意此時初始工作可以只有約束 3 或 5，也可以是空集，這裡我們只是舉一個 W0={3,5}W0={3,5} 的例子。當這兩條約束點都啟用時，迭代點 x0=(2,0)x0=(2,0) 顯然是當前工作集下的最優點，即 p0=0p0=0。因此我們有 x1=(2,0)Tx1=(2,0)T。當 p0=0p0=0 時，我們需要檢查對偶變數 λλ 的符號。計算可以得到 (λ^3,λ^5)=(−2,−1)(λ^3,λ^5)=(−2,−1)。因為不是所有的對偶變數都是正值，說明當前迭代點不是最優點，因此我們需要從工作集中刪除一條約束，這裡我們選擇刪除對應的對偶變數的絕對值最大的值所對應的約束，即第 3 條約束。

即此時工作集內只有第 5 條約束， W1={5}W1={5} 。在當前的工作集下計算得到的 p1=(−1,0)Tp1=(−1,0)T 。因為在前進的路上沒有 blocking constraint 的阻礙，因此計算得到的步長 α1=1α1=1。同時我們可以得到新的迭代點 x2=(1,0)Tx2=(1,0)T。

因為步長 α1=1α1=1，因此我們在工作集不變的情況下，即 W2={5}W2={5}，再求解一次等式約束命題，此時我們會得到前進方向 p2=0p2=0，因此更新後的迭代點為 x3=(1,0)Tx3=(1,0)T。因此我們去檢查各個對偶變數的正負號，通過計算可得 λ^5=−5λ^5=−5，說明當前點不是最優點，同時我們將約束 5 從工作集中刪除。

此時工作集為空集，即 W3=∅W3=∅。此時求解下降方向就相當於在求解一個無約束命題，從圖中可以看出，下降方向為 p3=(0,2.5)p3=(0,2.5)。然而此時，在下降的方向上有 blocking constraint 的存在，因此計算得到的步長為 α3=0.6α3=0.6，同時我們得到新的迭代點為 x4=(1,1.5)Tx4=(1,1.5)T，同時我們將 blocking constraint 新增到工作集中得到更新後的工作集為 W4={1}W4={1}。

此時進入下一次迭代，求解得到的下降方向為 p4=(0.4,0.2)p4=(0.4,0.2)。並且在前進的方向上沒有 blocking constraint，因此求得的步長為 1 。得到的新的迭代點為 x5=(1.4,1.7)Tx5=(1.4,1.7)T。因為步長為 1，因此我們在不改變工作集的情況下再計算一次下降方向，得到 p5=0p5=0，此時我們檢查對偶變數的值，發現 λ^1=0.8λ^1=0.8 是正值，說明當前的迭代點已經是最優點。因此我們退出迭代，並得到問題的最優解 x∗=(1.4,1.7)Tx∗=(1.4,1.7)T。

通過上面這張圖，我們可以較為清楚地瞭解 ASM 演算法的計算流程。

後記 最後要提一點，關於演算法的收斂性，原書中的定理 16.6 證明了，只要每次迭代的前進方向 pkpk 是從上文所提到的等式約束命題求解出來的，那麼就可以證明，目標函式的值沿著該方向一定會減少，而這就保證了ASM 演算法的迭代可以在有限次數內終止，從而證明的演算法的收斂性。那麼其實演算法的收斂性不是靠在新增約束和刪除約束的對約束的選擇來保證的，那就說明，在具體的應用中我們可以根據情況選擇新增和刪除約束的策略，從而實現演算法的快速收斂。

∗∗ 所謂判斷一個向量是否為0在程式碼實現時是一個很微妙的問題，根據我自己的經驗。當採用 double 精度實現時，如果所有元素的絕對值小於 10−610−6 就可以認為 pkpk 是 0，而對於 float 精度的實現，這一閾值約為 10−310−3。 原文：https://blog.csdn.net/dymodi/article/details/50358278