化問題的工程實(shí)踐與避坑指南)
簡介這份資源聚焦SCASequential Convex Approximation凸優(yōu)化算法的實(shí)現(xiàn)面向具備一定數(shù)學(xué)優(yōu)化基礎(chǔ)、希望將非凸問題轉(zhuǎn)化為連續(xù)凸近似求解的學(xué)習(xí)者與工程人員可用于無線通信、信號處理、能源系統(tǒng)優(yōu)化等場景的算法驗(yàn)證。壓縮包共2個文件均為MATLAB腳本.m整體約3KB體量輕便便于直接閱讀與二次修改。其中一份腳本給出特定功率問題的示例建模另一份則承載通用SCA算法框架二者配合可幫助讀者理解凸函數(shù)與凸集、凸優(yōu)化問題形式、近似構(gòu)造、迭代更新及收斂性分析等關(guān)鍵環(huán)節(jié)。目前已有2699人學(xué)習(xí)下載說明其在同類算法資料中具備一定參考價值。讀者可借此掌握Taylor展開、松弛等近似技術(shù)并對照源碼梳理從近似、求解到更新的完整流程進(jìn)而定制適配自身問題的SCA實(shí)現(xiàn)。1. SCA 與凸優(yōu)化從一個「非凸卡死」的現(xiàn)場說起如果你調(diào)過通信里的波束成形、搞過 IRS 相控陣、或者碰過機(jī)器學(xué)習(xí)里的稀疏正則大概率遇到過這種場面目標(biāo)函數(shù)寫出來挺漂亮一求導(dǎo)發(fā)現(xiàn)非凸CVX 直接報Disciplined convex programming error梯度下降跑一夜 loss 在幾個局部極小值之間反復(fù)橫跳。這時候老工程師通常會甩給你三個字母——SCASuccessive Convex Approximation逐次凸近似。它干的事很樸素既然原問題非凸那我不直接解你我在當(dāng)前點(diǎn)附近構(gòu)造一個凸的「替身」解替身拿解當(dāng)新起點(diǎn)再構(gòu)造再解直到收斂。凸優(yōu)化在這里不是目的是工具SCA 是把非凸問題「翻譯」成一系列凸問題的翻譯官。這套方法適合誰適合手里已經(jīng)有凸優(yōu)化求解器CVX、CVXPY、MOSEK、SCS 都行、但被非凸約束或非凸目標(biāo)卡住的從業(yè)者。下面我按自己踩過的順序把 SCA 從原理到能跑通的代碼、參數(shù)、翻車點(diǎn)講一遍。2. SCA 憑什么能收斂凸近似、代理函數(shù)與三個前提2.1 非凸問題為什么不能直接丟給求解器先把話說死凸優(yōu)化求解器的「凸」不是建議是硬約束。CVX 這類工具用的是 disciplined convex programming它靠一套規(guī)則判斷你寫的表達(dá)式是否滿足凸性組合一旦出現(xiàn)變量相乘、變量做分母、log 里套非凹函數(shù)直接拒絕。這不是求解器弱而是凸問題的全局最優(yōu)性保證依賴于可行域凸、目標(biāo)凸這兩個條件破壞任何一個KKT 點(diǎn)就不再等價于全局最優(yōu)。現(xiàn)實(shí)里的問題偏偏愛非凸。舉幾個我常碰到的波束成形里功率約束下的和速率最大化速率是 log(1SINR)SINR 里分子分母都含優(yōu)化變量整體非凸稀疏感知里的 L0 范數(shù)組合性質(zhì)天然非凸IRS 場景里反射相移和信道耦合雙線性結(jié)構(gòu)。這些問題的共同點(diǎn)是——局部看近似凸全局看不是。SCA 的思路就是承認(rèn)「全局我搞不定」轉(zhuǎn)而做「局部可靠」。它在第 k 次迭代點(diǎn) x_k 處找一個凸函數(shù) g(x; x_k)滿足兩個條件一是 g(x_k; x_k) f(x_k)在當(dāng)前點(diǎn)值和原函數(shù)相等二是 g(x; x_k) ≥ f(x)或 ≤取決于最大化還是最小化即代理函數(shù)是原函數(shù)的全局上界或下界。滿足這兩條解代理函數(shù)得到的解一定不會讓原目標(biāo)變差這就是單調(diào)性來源。2.2 代理函數(shù)怎么造一階泰勒、二次上界與 MM 框架造代理函數(shù)是 SCA 的核心手藝常見三條路。第一條是一階泰勒展開。對凸函數(shù)泰勒展開是全局下界對凹函數(shù)泰勒展開是全局上界。所以最大化一個凹函數(shù)時直接在當(dāng)前點(diǎn)做一階泰勒得到的就是合法的凹下界代理。這是最省事的一類很多速率最大化問題里 log(1SINR) 對某些變量是凹的直接泰勒就能用。第二條是二次上界典型場景是處理 log 和分式。比如 log(1x) 這種凹函數(shù)除了泰勒還可以用二次函數(shù)在展開點(diǎn)處構(gòu)造上界收斂性質(zhì)更好但計算稍重。分式結(jié)構(gòu)常用二次變換quadratic transform把分子分母解耦成可交替優(yōu)化的形式。第三條是 MM 框架Majorization-Minimization。它不要求代理函數(shù)是泰勒只要求是原函數(shù)的 majorizer上界且在當(dāng)前點(diǎn)緊。MM 和 SCA 經(jīng)常被混著叫區(qū)別在于 MM 更強(qiáng)調(diào)「majorize」這個構(gòu)造動作SCA 更強(qiáng)調(diào)「逐次凸化」這個流程。實(shí)操里我基本把它們當(dāng)一回事用。提示代理函數(shù)必須滿足「在當(dāng)前點(diǎn)緊」這個條件否則單調(diào)性不成立。我見過有人隨手寫個上界但當(dāng)前點(diǎn)不相等結(jié)果迭代震蕩查了半天以為是步長問題。2.3 收斂性依賴的三個前提缺一個就翻車SCA 的收斂保證不是白給的它依賴三個前提我按重要性排。第一代理函數(shù)在當(dāng)前點(diǎn)必須緊且是原函數(shù)的全局上界最大化時或下界最小化時。這條破了單調(diào)性沒了收斂無從談起。第二每次子問題要解到足夠精度。理論上要求解到全局最優(yōu)實(shí)操里如果子問題只解了個近似外層迭代可能停在偽收斂點(diǎn)。CVXPY 里我會把solver的eps調(diào)緊一點(diǎn)別用默認(rèn)的粗精度。第三迭代序列要有界或者目標(biāo)函數(shù)要有下界。如果問題本身無界SCA 會一路發(fā)散。這個在功率控制里要特別注意功率上界約束必須寫死。這三條里第一條是設(shè)計問題第二條是工程問題第三條是建模問題。我踩過的坑基本都落在這三類的某一類里。3. 用 Python CVXPY 跑通一個 SCA 最小例子3.1 選一個能體現(xiàn)非凸性的最小問題為了讓你能直接抄我選一個足夠小但確實(shí)非凸的問題最大化 sum(log(1 x_i * a_i))約束是 sum(x_i) ≤ Px_i ≥ 0。這里 a_i 是給定正系數(shù)x_i 是優(yōu)化變量。這個問題本身其實(shí)是凹的log 是凹復(fù)合仿射保持凹所以它不算真非凸。為了制造非凸我把目標(biāo)改成 sum(log(1 x_i * a_i)) - c * sum(x_i^2)加一個凹的負(fù)二次項(xiàng)整體就非凹了。這個結(jié)構(gòu)在能效優(yōu)化里很常見速率減去功耗懲罰。原問題maximize sum(log(1 a_i * x_i)) - c * sum(x_i^2) subject to sum(x_i) P x_i 0非凸來源是 -c * sum(x_i^2)它是凹函數(shù)但前面是負(fù)號整體目標(biāo)變成凹減凹不保證凹。SCA 的處理是把 -c * sum(x_i^2) 在當(dāng)前點(diǎn)做凹函數(shù)的泰勒展開得到全局上界然后最大化這個上界。3.2 代理函數(shù)的代碼實(shí)現(xiàn)import cvxpy as cp import numpy as np np.random.seed(0) N 8 a np.random.rand(N) 0.5 c 0.1 P 5.0 x cp.Variable(N, nonnegTrue) # 原目標(biāo)非凸不能直接丟給 CVXPY # obj cp.sum(cp.log(1 cp.multiply(a, x))) - c * cp.sum_squares(x) # SCA 外層迭代 x_k np.ones(N) * (P / N) # 初始點(diǎn)均勻分配 max_iter 50 tol 1e-4 for it in range(max_iter): # 構(gòu)造代理-c*sum(x^2) 在 x_k 處的凹泰勒上界 # f(x) -c*x^2, f(x) f(x_k) f(x_k)*(x - x_k) # -c*x_k^2 - 2*c*x_k*(x - x_k) # 常數(shù)項(xiàng)不影響 argmax可省略 linear_term -2 * c * x_k # 梯度系數(shù) proxy cp.sum(cp.log(1 cp.multiply(a, x))) linear_term x constraints [cp.sum(x) P] prob cp.Problem(cp.Maximize(proxy), constraints) prob.solve(solvercp.ECOS, verboseFalse) x_new x.value # 用原目標(biāo)評估真實(shí)進(jìn)展 true_obj np.sum(np.log(1 a * x_new)) - c * np.sum(x_new ** 2) if it % 5 0: print(fiter {it:3d} true_obj {true_obj:.6f}) if np.linalg.norm(x_new - x_k) tol: print(fconverged at iter {it}) break x_k x_new print(final x , np.round(x_new, 4))這段代碼的邏輯分三層。第一層原目標(biāo)里的 log 項(xiàng)是凹的保留不動二次項(xiàng) -csum(x^2) 是凹的但我們要最大化它凹函數(shù)最大化本身沒問題問題在于它和 log 項(xiàng)加在一起后整體不保證凹——實(shí)際上 log 是凹-x^2 也是凹兩個凹相加還是凹等等這里我得糾正自己凹加凹確實(shí)是凹所以這個例子其實(shí)還是凹的。為了真正制造非凸得讓二次項(xiàng)帶正號即 csum(x^2)這樣凹加凸整體非凹。我重新調(diào)整目標(biāo)改成 sum(log(1 a_i * x_i)) - c * sum(x_i^2) 里把 -c 改成 c即 sum(log(1a_i x_i)) c * sum(x_i^2)最大化它。此時 log 凹x^2 凸凹加凸非凹。SCA 處理凸項(xiàng) x^2 時因?yàn)橐畲蠡购瘮?shù)的最大化不能直接做需要對凸函數(shù)做全局下界凸函數(shù)的泰勒展開是全局下界即 x^2 ≥ x_k^2 2 x_k (x - x_k)。用這個下界替換代理函數(shù)變成凹的可以最大化。# 修正后的代理構(gòu)造 for it in range(max_iter): # 對凸項(xiàng) c*sum(x^2) 做全局下界凸函數(shù)泰勒展開 # x^2 x_k^2 2*x_k*(x - x_k) # 常數(shù)項(xiàng)省略線性系數(shù)為 2*c*x_k linear_term 2 * c * x_k proxy cp.sum(cp.log(1 cp.multiply(a, x))) linear_term x constraints [cp.sum(x) P] prob cp.Problem(cp.Maximize(proxy), constraints) prob.solve(solvercp.ECOS) x_new x.value true_obj np.sum(np.log(1 a * x_new)) c * np.sum(x_new ** 2) if np.linalg.norm(x_new - x_k) tol: break x_k x_new3.3 參數(shù)怎么設(shè)初始點(diǎn)、步長與停止條件初始點(diǎn) x_k 的選擇直接影響收斂速度和能不能收斂到好點(diǎn)。我一般用均勻分配或者可行域內(nèi)的隨機(jī)點(diǎn)跑幾次取最好的。均勻分配在功率分配問題里通常不差因?yàn)閷ΨQ性。停止條件我用兩個變量變化量 norm(x_new - x_k) tol或者連續(xù)兩次原目標(biāo)變化小于某個閾值。tol 取 1e-4 到 1e-6 之間看問題尺度。如果變量量級是 1e3tol 要相應(yīng)放大。求解器選擇上ECOS 適合小規(guī)模二階錐問題SCS 適合大規(guī)模但精度粗MOSEK 最穩(wěn)但要 license。CVXPY 里可以指定solvercp.ECOS如果報 solver 不支持換cp.SCS試試。注意子問題求解精度別用默認(rèn)值。ECOS 默認(rèn)abstol1e-8還行SCS 默認(rèn)精度很粗外層迭代會被子問題的誤差帶偏表現(xiàn)為原目標(biāo)曲線鋸齒狀。3.4 怎么驗(yàn)證你真的在收斂光看變量變化不夠要看原目標(biāo)序列。我習(xí)慣把每次迭代的真實(shí)目標(biāo)打出來畫一條曲線。健康的 SCA 曲線是單調(diào)上升最大化問題然后趨于平緩。如果出現(xiàn)下降說明代理函數(shù)構(gòu)造錯了或者子問題沒解到最優(yōu)。如果曲線一直上升不收斂檢查約束是不是沒寫全問題可能無界。另一個驗(yàn)證手段是跑多個初始點(diǎn)看收斂到的目標(biāo)值是否接近。如果差異很大說明問題有多個局部最優(yōu)SCA 只能保證局部這時候要么換更好的初始點(diǎn)要么考慮全局化策略。4. 避坑與排查SCA 落地時最容易翻車的五件事4.1 代理函數(shù)方向搞反迭代直接發(fā)散現(xiàn)象原目標(biāo)曲線一路下降或者震蕩幅度越來越大。原因最大化問題里代理函數(shù)必須是原函數(shù)的全局上界我見過有人把凸函數(shù)的泰勒展開當(dāng)上界用實(shí)際上凸函數(shù)泰勒是下界方向反了。解決每次構(gòu)造完代理在當(dāng)前點(diǎn)驗(yàn)證 g(x_k) f(x_k)再隨機(jī)取一個點(diǎn)驗(yàn)證 g(x) f(x)最大化或 g(x) f(x)最小化。這個檢查寫成一個 assert能省掉大量調(diào)試時間。4.2 子問題求解精度不夠偽收斂現(xiàn)象變量變化量很小看起來收斂了但換一個求解器或者調(diào)緊精度后目標(biāo)還能再漲一截。原因子問題沒解到全局最優(yōu)外層迭代停在了子問題誤差造成的偽駐點(diǎn)。解決把子問題求解器的精度調(diào)緊ECOS 用abstol1e-9, reltol1e-9SCS 用eps1e-6并增加max_iters。如果子問題規(guī)模大考慮換 MOSEK。4.3 約束非凸但被忽略解出來不可行現(xiàn)象求解器返回成功但把解代回原問題約束被違反。原因SCA 只凸化了目標(biāo)約束里的非凸部分沒處理或者處理時用了錯誤的近似方向。解決約束的非凸部分同樣需要凸化比如 x*y b 這種雙線性約束在當(dāng)前點(diǎn)對其中一個變量做泰勒固定另一個。檢查方法是把最終解代入所有原始約束逐條驗(yàn)證。4.4 初始點(diǎn)不可行第一步就報 infeasible現(xiàn)象CVXPY 報Problem status: infeasible。原因初始點(diǎn)不滿足約束而代理子問題的可行域雖然包含原可行域但如果初始點(diǎn)離可行域太遠(yuǎn)子問題可能無解。解決初始點(diǎn)必須選在可行域內(nèi)。功率分配里就是均勻分配相移優(yōu)化里就是全零相位總之先找一個滿足所有約束的點(diǎn)。4.5 收斂判據(jù)只看變量忽略目標(biāo)尺度現(xiàn)象變量變化量小于 tol 但目標(biāo)還在緩慢改善或者變量變化量很大但目標(biāo)幾乎不變。原因變量尺度和目標(biāo)尺度不匹配單一判據(jù)不可靠。解決同時監(jiān)控變量變化和目標(biāo)變化兩個都小于閾值才停。目標(biāo)變化閾值取1e-6 * abs(true_obj)這種相對量比絕對量穩(wěn)。5. 進(jìn)階技巧把 SCA 嵌進(jìn)交替優(yōu)化與加速收斂5.1 塊坐標(biāo)下降 SCA多變量耦合時的標(biāo)準(zhǔn)打法很多問題有多個變量塊比如波束成形里的發(fā)射波束和 IRS 相移兩者耦合導(dǎo)致聯(lián)合非凸。標(biāo)準(zhǔn)做法是塊坐標(biāo)下降BCD固定相移優(yōu)化波束固定波束優(yōu)化相移每塊內(nèi)部用 SCA。這樣每塊都是凸子問題交替求解。收斂性上BCD 加 SCA 能保證目標(biāo)單調(diào)但收斂速度取決于塊之間的耦合強(qiáng)度。耦合強(qiáng)的時候交替次數(shù)會很多我一般設(shè)最大交替次數(shù) 20 到 50配合目標(biāo)變化閾值提前停。5.2 用外推加速Nesterov 和 Anderson 加速的取舍SCA 的收斂速度通常是次線性的迭代次數(shù)多。加速手段有兩類一是 Nesterov 外推在代理函數(shù)里加動量項(xiàng)二是 Anderson 加速用歷史迭代點(diǎn)做線性組合。Nesterov 實(shí)現(xiàn)簡單但步長參數(shù)不好調(diào)調(diào)不好會破壞單調(diào)性。Anderson 加速更穩(wěn)但需要存歷史點(diǎn)內(nèi)存開銷大。我的經(jīng)驗(yàn)是問題規(guī)模小變量少于 100用 Anderson規(guī)模大用 Nesterov 或者干脆不加速因?yàn)榧铀賻淼氖找婵赡鼙幻坎降念~外計算抵消。5.3 收斂性驗(yàn)證清單跑完 SCA我習(xí)慣做三件事確認(rèn)結(jié)果可信。第一把最終解代入原問題檢查所有約束滿足程度違反量應(yīng)該在求解器精度范圍內(nèi)。第二從不同初始點(diǎn)跑三到五次看目標(biāo)值分布如果方差很小說明局部最優(yōu)解質(zhì)量穩(wěn)定。第三把原目標(biāo)曲線畫出來確認(rèn)單調(diào)性如果有下降段回去查代理函數(shù)。下面這張表是我常用的參數(shù)配置按問題規(guī)模分問題規(guī)模求解器子問題精度停止 tol最大迭代變量 50ECOS1e-91e-610050 ~ 500SCS1e-61e-5200 500MOSEK1e-81e-4300這張表不是金科玉律但能讓你少走彎路。我早期用 SCS 默認(rèn)精度跑小問題結(jié)果外層迭代 200 次還沒收斂換成 ECOS 后 30 次就停了血淚經(jīng)驗(yàn)。5.4 一個我常犯的錯誤最后說個我自己的教訓(xùn)。有次做 IRS 相移優(yōu)化SCA 跑出來目標(biāo)值比預(yù)期低很多查了兩天以為是代理函數(shù)錯了最后發(fā)現(xiàn)是初始相位設(shè)成了全零而全零相位在這個場景里恰好是個很差的局部點(diǎn)SCA 從那兒出發(fā)就再也沒爬出來。后來改成隨機(jī)相位跑五次取最好目標(biāo)直接漲了 15%。SCA 是局部方法初始點(diǎn)的重要性怎么強(qiáng)調(diào)都不過分別在初始點(diǎn)上省事。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取