劃(SDP)核心原理、應用與求解器實戰(zhàn))
整理2021年11月8日的筆記時翻到一頁紙標題只寫了半行字半定規(guī)劃SDP。那天我在調一個圖劃分問題的松弛折騰到很晚隨手記下來的SDP要點如今回頭看竟然成了之后好幾個項目的解題地基。半定規(guī)劃這個名字聽起來很數(shù)學但本質其實不難它就是把線性規(guī)劃里的“非負向量”約束換成“半正定矩陣”約束。它之所以值得花時間去啃是因為它能做很多線性規(guī)劃做不了的事——給組合優(yōu)化提供強松弛、把非凸的秩問題轉換成凸問題、把控制理論里的穩(wěn)定性條件寫成可求解的矩陣不等式。這篇文章適合正在跟整數(shù)規(guī)劃、非線性約束、低秩恢復打交道或者被某個看似簡單卻死活解不動的小模型卡住的人讀。我會把核心原理、應用范式、求解器選型以及實際踩坑的過程都捋一遍盡量說人話。1. 線性松弛為什么擋不住這種問題從最大割說起1.1 一個讓線性松弛當場翻車的例子最大割問題是SDP最經(jīng)典的出場場景。問題本身很直白給你一張無向圖每條邊帶權重把頂點分成兩組讓“兩端不在同一組”的邊權重總和盡量大。這是個NP-hard問題聚類、電路劃分、網(wǎng)絡設計里都藏著它的變體。第一次接觸這類問題時大多數(shù)人會嘗試整數(shù)規(guī)劃建模。對每條邊(i,j)引入變量表示兩端是否同組最直接的做法是用x_i±1編碼組別然后最大化∑ w_ij · (1 - x_i x_j) / 2這里x_i x_j1表示同組該邊不貢獻x_i x_j-1表示異組貢獻w_ij。理想情況是讓x_i只取±1但那是整數(shù)規(guī)劃不好解。所以經(jīng)典的松弛思路是允許x_i落在[-1,1]區(qū)間把問題變成線性規(guī)劃。這個松弛在理論上是成立的——可行域變大了最優(yōu)值只會更好。但問題是它好得過頭了。我在一個模擬項目X里測過一張單位權隨機圖整數(shù)最優(yōu)割大約能拿到總邊權的84%線性松弛給出來的上界卻在92%以上。這個上界虛高到?jīng)]法用來做分支定界的剪枝幾乎等于沒給信息。原因在于[-1,1]這個盒子太松松弛解可以讓每條邊都“半貢獻”把這個割的邊數(shù)同時掙滿這在真正的離散解里根本做不到。換成SDP之后情況完全不一樣。做法是把標量x_i升級成向量v_i約束張成Gram矩陣X其中X_ii1X_ijv_i^T v_j。允許X_ij在[-1,1]之間連續(xù)變化但要求整個Gram矩陣半正定。同樣是那張隨機圖SDP給出來的上界只比整數(shù)最優(yōu)高兩三個百分點。這就是SDP在組合優(yōu)化里最核心的價值松弛質量緊邊界有用。1.2 為什么SDP松弛更緊一個幾何直覺很多人第一次看到“把標量升級成向量”會覺得莫名其妙這不是把簡單問題搞復雜了嗎但如果從幾何角度理解就順了。線性松弛相當于把{-1,1}^n的頂點集合放進一個n維立方體。立方體是凸的好優(yōu)化但它的體積遠大于離散點真正的凸包尤其在組合結構復雜時立方體的很多角落在組合意義上根本不可達。SDP則不同它把每個頂點映射為單位球上的一個向量可行域變成所有Gram矩陣的集合——你可以讓某些邊的內積給到部分的貢獻但這些內積必須來自一組真實存在的歐氏空間向量不能每個方向都撒謊。舉個例子三條邊的三角形如果想讓三條邊同時貢獻接近滿分線性松弛可以說“每條邊都切一半”三個x_i_i的對都取接近0目標值虛高。但SDP里三個向量兩兩夾角不能同時接近180°Gram矩陣半正定這個約束會把這種不可能直接判死。所謂“半正定約束沒有時間先后它管的是全局一致性”這就是SDP比LP緊的幾何本質。2. 半定錐、LMI與對偶SDP的數(shù)學骨架2.1 標準形式從頭到尾只有三個要素SDP的標準形式寫出來非常簡短min tr(C X) s.t. tr(A_i X) b_i, i 1,...,m X ? 0且X為對稱矩陣這里tr(AX)是矩陣內積等于把兩個矩陣對應元素相乘再求和。X ? 0表示X是半正定矩陣所有特征值非負。三個要素分別是線性目標、線性等式約束、半定矩陣錐約束。跟線性規(guī)劃相比唯一的變化就是約束從“x ≥ 0”變成了“X ? 0”。但這一變表達力就完全不同了。標量版本里x_i≥0是一些獨立的半射線矩陣版本里X ? 0是整個半定錐它對矩陣整體結構施加約束能把“這些向量之間滿足某種內積關系”這種全局性質直接表達出來。很多看起來完全非線性的問題放進這個框架里反而變成凸問題。2.2 LMI與Schur補把二次約束翻譯進半定世界線性矩陣不等式LMI是SDP最常見的表述形式。它一般寫成F(x) A_0 ∑ x_i A_i ? 0其中A_i是常數(shù)對稱矩陣x是決策向量。只要目標函數(shù)是線性的這整個問題就是一個SDP。之所以強調LMI是因為大量工程約束都不是天然寫成SDP的標準形式而是能通過LMI等價轉化。最常用的翻譯工具是Schur補。比如我要表達“t ≥ x^T x”這個二次約束它本身是凸的但沒法寫進線性規(guī)劃。通過Schur補可以寫成[ I, x ] [ x^T, t ] ? 0這個分塊矩陣半正定等價于t ≥ x^T x??雌饋碇皇菗Q了個寫法但在求解上是從“非線性約束”變成了“矩陣錐約束”直接落入SDP的射程。我做傳感器定位項目時這個技巧幾乎天天用——把距離約束的平方項塞進矩陣的角落問題就從“不知道從哪下手”變成了“標準的凸優(yōu)化”。2.3 對偶與Slater條件求解器真正在解什么SDP的拉格朗日對偶寫出來跟原始問題長得非常對稱max b^T y s.t. C - ∑ y_i A_i ? 0弱對偶始終成立任意原始可行點的目標值都不小于任意對偶可行點。但真正讓求解器能高效工作的是強對偶也就是原始最優(yōu)值等于對偶最優(yōu)值而這需要在原始問題里存在一個嚴格可行點某個X滿足等式約束并且X嚴格正定所有特征值都大于0。這個條件叫Slater條件。大多數(shù)SDP求解器同時維護原始路徑和對偶路徑靠兩者之間的對偶間隙來判斷收斂。如果問題里沒有嚴格可行點對偶路徑會一直不穩(wěn)求解器要么慢得要命要么直接報不可行。我在實際項目里吃過這方面的虧——模型的理論形式看著沒問題等式約束也都對得上但變量被約束在低秩流形上整個可行集里找不到一個嚴格正定的點求解器就是不給面子。解決辦法通常是在目標里加一個很小的正則項比如εtr(X)或者把某些等式約束放寬到帶ε的區(qū)間把問題推向內點存在的情形。2.4 多項式時間可解為什么SDP能上工程SDP是凸優(yōu)化最經(jīng)典的算法是內點法。它的思路是在目標函數(shù)上疊加一個對數(shù)障礙項比如min tr(CX) - μ log det(X)然后沿著μ遞減的路徑做牛頓迭代。每步牛頓法要分解一個大型矩陣代價大約O(n^3)但迭代次數(shù)只需要約O(√n log(1/ε))——這里的n是矩陣維度ε是目標精度。這意味著求解時間隨問題規(guī)模是多項式增長而不是指數(shù)爆炸這是SDP能用于實際工程的根本保證。不要被“矩陣”這個外殼嚇住從計算復雜度框架看SDP和LP非常接近只是每一步要多拆一個更大的矩陣。實際表現(xiàn)上幾千維以內的SDP在主流求解器上都能在幾秒到幾分鐘內收斂到1e-6左右的精度完全夠用。3. 四個能直接抄的SDP應用范式松弛、核范數(shù)、李雅普諾夫與定位3.1 最大割的經(jīng)典套路SDP求解隨機舍入最大割的完整流程值得認真走一遍它是理解SDP后處理的標準模板。第一步建模目標。設X是Gram矩陣約束X_ii1最大化∑ w_ij (1 - X_ij)/2。求解SDP得到最優(yōu)的X*。第二步特征分解。把X*寫成VΛV^T取特征向量矩陣V每一行就是原問題頂點對應的高維向量v_i。第三步隨機舍入。生成一個標準高斯隨機向量r對所有頂點計算s_isign(r^T v_i)s_i為1的頂點歸入A組-1歸入B組。第四步重復多次。采樣幾百上千個隨機方向取切割權重最大的一次作為最終解。這個流程里最關鍵的理解點是SDP解出來的不是一組離散標簽而是一組連續(xù)向量。標簽是在后處理階段“舍入”出來的舍入質量直接決定最終結果。理論上有保證期望切割權重不低于SDP最優(yōu)值的0.878倍這是經(jīng)典結果。實際跑下來通常更好尤其在冗余采樣幾百組隨機向量后經(jīng)常能拿到0.9以上的近似比。3.2 核范數(shù)最小化把非凸秩問題變成凸問題推薦系統(tǒng)里的矩陣補全圖像修復里的低秩重構本質都在做同一件事有一個大矩陣只知道部分元素想恢復整個矩陣而且相信真實矩陣是低秩的。直接最小化rank(X)是非凸的幾乎沒法解。但rank函數(shù)有一個凸包絡核范數(shù)也就是所有奇異值之和。把目標換成min ||X||_*問題變成凸的。更有意思的是核范數(shù)可以寫成SDPmin (1/2)(tr(U) tr(V)) s.t. [ U, X; X^T, V ] ? 0這個等價關系在工程上意味著什么意味著低秩恢復可以被一套標準的凸優(yōu)化工具鏈處理不需要給每個問題手寫專用算法。這在以前是不可想象的。我在一個圖像去噪的demo里試過噪聲污染了大約60%的像素點核范數(shù)SDP恢復出來的矩陣在視覺上幾乎看不出缺失痕跡。3.3 李雅普諾夫穩(wěn)定性控制里的SDPSDP不只是離散優(yōu)化和機器學習的工具控制理論里它幾乎是標配。一個線性系統(tǒng)dx/dt Ax怎么判斷它是否穩(wěn)定可以找一個李雅普諾夫函數(shù)V(x)x^T P x如果P正定并且V沿系統(tǒng)軌跡單調下降系統(tǒng)就穩(wěn)定。這個條件寫成矩陣不等式就是P ? 0, A^T P P A ? 0兩個都是矩陣不等式變量是P這正是標準SDP。求解之后不僅能判斷穩(wěn)定性還能拿P去構造控制器或者在滿足LMI約束的前提下進一步優(yōu)化某個性能指標比如最小化某個衰減率這些追加約束都還是SDP。這個應用提醒我SDP的覆蓋面遠比單純“組合優(yōu)化松弛”要寬它本身就是處理矩陣層面約束的通用語言。3.4 傳感器定位提升半定松弛丟棄秩約束傳感器定位問題是一個特別好的范式平面上布了N個節(jié)點只知道部分兩兩距離測量值想估計每個節(jié)點的真實坐標。約束長這樣(p_i - p_j)^T (p_i - p_j) d_ij^2這是非凸的二次等式直接解非常困難。SDP的做法是把這個非線性問題“提升”到矩陣空間。定義P為坐標矩陣d×n構造分塊矩陣Y [I, P; P^T, Q]其中Q是對稱矩陣。注意如果Q恰好等于P^T P那么Y自動半正定而且秩等于d。所以原始問題等價于Y ? 0、Y的某些固定塊已知、外加rank(Y)d。丟掉秩約束保留半正定和固定塊約束就是一個SDP。這個“提升到矩陣 半定松弛 丟秩”的配方可以遷移到很多幾何問題上包括相機位姿估計、分子構象恢復等。我自己的體會是看到“二次等式約束”先別慌先想能不能把它塞進一個分塊矩陣的Schur補里往往答案就出來了。4. 求解器實戰(zhàn)場內點法、一階法與建模接口怎么選4.1 兩類主流求解器精度還是規(guī)模半定規(guī)劃的求解器大體可以分兩類選錯了會非常難受。類型每步代價迭代次數(shù)適用規(guī)模典型表現(xiàn)內點法型O(n^3)約O(√n log(1/ε))幾千維以內穩(wěn)健gap能到1e-8一階法型O(nnz)依賴稀疏度幾百到幾萬受精度影響大超大規(guī)模稀疏矩陣快但高精度收斂慢內點法型的看家本領是精度和穩(wěn)健性適合驗證模型、產(chǎn)出基準解、處理中小規(guī)模問題。一階法型更適合超大規(guī)模、稀疏結構明顯、只需要近似解做快速探索的場景。我現(xiàn)在的習慣是首次接觸一個問題永遠先用內點法在小規(guī)模數(shù)據(jù)上跑出目標值的基準再切到一階法去逼近。否則你會分不清誤差是來自松弛模型本身還是來自算法沒收斂。4.2 建模接口三行代碼聲明一個半定變量SDP建模接口已經(jīng)做得很成熟核心抽象就三個元素半定矩陣變量、等式約束、線性目標。以某通用凸優(yōu)化建模接口為例最大割的SDP松弛寫起來大約是這樣的概念代碼# 概念代碼最大割SDP松弛 X Semidefinite(n) # 聲明n×n半定矩陣變量 constraints [X[i][i] 1 for i in range(n)] objective sum(w[i][j] * (1 - X[i][j]) / 2 for i j) prob Minimize(-objective) # 等價于最大化 prob.solve()這段代碼幾乎和數(shù)學公式一一對應不用手動推導對偶形式不用管內點法的具體步進策略。聲明變量、寫約束、設目標、調求解器四步走完。這也是我推薦每個人都應該學一點SDP建模的原因數(shù)學底子要求不高收益卻非常實在。4.3 求解參數(shù)與終止判據(jù)調什么、什么時候調拿回一個SDP模型求解參數(shù)通常不需要大動但有幾個點值得了解。第一是初始點給X傳一個合理的正定矩陣比如單位陣的倍數(shù)往往能顯著減少收斂時間。第二是終止間隙gap建模階段如果只是探索解的行為gap設在1e-3到1e-4足夠如果要拿去做后續(xù)精確驗證壓到1e-7更穩(wěn)。第三是迭代上限大規(guī)模一階法給一個較大的上限然后盯著gap下降曲線看而不是干等它結束。如果返回狀態(tài)是原始不可行或無界不要急著改約束先跑一個小規(guī)模代理問題把目標系數(shù)縮小十倍看問題形態(tài)有沒有變化。無界問題通常表現(xiàn)為目標函數(shù)在某個可行方向上能無限下降不可行問題則要檢查是否有冗余等式約束或者你手動構造一個明顯可行的X代入約束驗證一下。這兩件事做下來大部分建模初期的問題都能定位。5. 我踩過的三個SDP數(shù)值坑從報錯到修復的完整鏈路5.1 坑一約束系數(shù)尺度懸殊gap卡死不動現(xiàn)象很典型求解器迭代了上千次gap停在1e-3附近死活不下來目標值卻在微小的范圍內抖動。我的排查鏈路是先打印約束矩陣A_i的元素絕對值和分布。結果發(fā)現(xiàn)有些約束的系數(shù)集中在1e7量級另一些只有1e-7。這樣量級懸殊的系數(shù)放在同一個SDP里等價于把一個橢球拉伸得極扁牛頓步在數(shù)值上非常不穩(wěn)定。SDP的內點法對尺度是敏感的盡管數(shù)學上同比例縮放不影響最優(yōu)解但在浮點數(shù)世界里等價不等于可計算。修復方式很簡單把所有約束矩陣和目標矩陣統(tǒng)一除以最大系數(shù)的絕對值也就是整體縮放到0.1到10的區(qū)間內。做完之后同一求解器gap很快降到1e-7。驗證方式也簡單兩次求解的目標值一致但收斂路徑和精度完全不同。這個坑是隱藏的因為目標值看著沒變如果把解拿去做特征分解特征向量會被數(shù)值噪聲污染得一塌糊涂后續(xù)舍入和恢復步驟全崩。5.2 坑二模型看似可行其實沒有Slater點有一次在模擬項目X里做傳感器定位求解器報原始不可行。我很確定約束方程組有解因為所有測量值都是真實的不存在矛盾。仔細排查后發(fā)現(xiàn)問題出在“幾何退化”。所有傳感器的真實坐標落在一個低維子空間里比如全部共線這導致任何滿足等式約束的Y矩陣在對應方向上的特征值自動變成0整個可行集里找不到一個嚴格正定的內點。也就是說雖然可行集非空但Slater條件不滿足對偶問題失效求解器只能判定原始不可行。修復的方法是加正則項給目標函數(shù)加上δ·tr(X)δ取1e-6左右相當于把可行集朝內點方向“推”了一下。這會對最優(yōu)值引入δ量級的擾動但對大多數(shù)工程應用來說可以接受。這個坑的教訓是SDP數(shù)學上寫著“半正定”求解器數(shù)值上要的是“正定”這兩個狀態(tài)之間的差距往往就是一次求解成功和失敗的分界線。5.3 坑三隨機舍入前忘了清洗特征向量最大割流程走到最后一步SDP解出來了gap已經(jīng)到1e-6隨機超平面舍入的結果卻慘不忍睹近似比只有0.6左右遠低于0.878的理論保障。一開始我以為是隨機種子問題多跑幾組還是沒有改善。后來查了解矩陣本身發(fā)現(xiàn)問題出在求解器默認容差下X*帶著微小的負特征值。特征分解之后負特征值對應在特征向量里的小尾巴會讓向量長度漂移有些向量的模長變成了0.985或者1.03而舍入算法假設它們嚴格落在單位球上符號判斷的準確率自然掉得厲害。修復方式是在舍入之前做清洗把特征值做clip負值截斷為0重新構造一個盡量貼近原始解的半正定矩陣再對所有向量歸一化到單位長度。清洗之后重新跑幾千組隨機超平面近似比穩(wěn)定回到了0.9附近。這一步的代碼其實只有幾行但效果立竿見影現(xiàn)在它已經(jīng)成為我所有SDP后處理流程的固定環(huán)節(jié)。5.4 建模完成后必做的三問檢查法踩過這么多坑之后我給自己定了一個SDP模型提交前的三問檢查第一問約束里有沒有冗余等式冗余約束會讓對偶尺度退化有時還會在Slater點附近制造數(shù)值陷阱能刪就刪。第二問是否存在嚴格正定的可行點如果X_iI或者某個合理構造的正定矩陣代入都不滿足約束就要考慮加正則項。第三問目標函數(shù)有沒有無用的常數(shù)項常數(shù)項不影響最優(yōu)解但會影響對偶gap的讀數(shù)干擾你看收斂曲線。這三問花不了五分鐘但能避免大多數(shù)“求解器莫名其妙不工作”的局面。最后再分享一個我最近一直在用的習慣每個SDP模型求解完成后別急著看目標值先打印最優(yōu)X的特征值譜。如果譜在某個k之后迅速衰減說明這個松弛的解本質上是低秩的后處理階段完全可以把尾部噪聲權重丟掉轉成普通的結構化問題如果譜從頭到尾都拖得很長說明模型大概率太松得回去補約束。這個檢查十秒鐘就能做完卻能省下大量在舍入階段試錯的時間。