简介SCASequential Convex Approximation算法示例压缩包面向凸优化与非凸优化问题研究者、算法学习者和MATLAB使用者特别适合需要处理非凸目标函数并追求全局最优解的工程场景。资源包内共2个.m文件大小仅3KB其中一个是可复用的SCA通用算法框架另一个是针对具体优化问题的示例实现两者可对照展示凸近似的构造思路、迭代更新方式以及收敛判定过程。借助两个脚本读者能快速掌握从目标函数近似、近似问题求解到解更新的完整SCA流程也能体会Taylor展开、松弛或截断等近似构造手法在实际代码中的落地方式进一步替换目标函数与约束即可适配不同领域的非凸优化任务。文件体量小、结构清晰适合在MATLAB中直接运行观察迭代效果也便于作为算法课堂演示或工程预研的基础模板。目前已有2692人学习下载适合在无线通信、信号处理、能源系统优化等领域开展算法演练与二次开发的读者。1. 为什么说 SCA 凸优化是处理非凸问题的默认解法做过多用户波束成形、相位恢复或传感器定位的工程师大概率有过这种经历目标函数写出来很好看SINR 分母里带干扰项把 log 一展开就是一堆非凸项CVX 里直接obj cp.sum(cp.log1p(...))会报错或给出毫无意义的警告。这时最常见的选择不是去套一个全局优化求解器而是打开 SCASuccessive Convex Approximation连续凸近似的迭代框架把非凸目标在当前迭代点附近凸化成子问题求解后更新再重复。SCA 凸优化的核心思想很朴素不试图一次性解决非凸问题而是在每个迭代点上用一个强凸的代理函数surrogate function去贴合原问题走一小步再重新贴合。这个“每次只走一小步”的策略使它可以处理大规模、高维的非凸问题并且能保证收敛到一阶稳定点。适合的人群包括通信系统设计、信号处理、机器学习稀疏建模、控制与机器人领域的算法工程师和研究员它不要求你有很强的凸分析基础但需要理解替代函数的构造条件和几个关键参数的含义。2. SCA 凸优化的核心框架替代函数、信任域与收敛条件2.1 从 MM 算法到 SCA全局上界条件为什么太苛刻SCA 的前身是 MMMajorization-Minimization算法。MM 要求替代函数是原函数的全局上界对任意 x替代函数值都要不小于原函数值。这个条件的优点是单调性证明非常漂亮但工程上几乎做不到。考虑一个常见的非凸项 f(x) -log(1 x²)它在 x 0 附近有凹陷两侧逐渐回平。想找一个全局上界函数其实并不难但要让它在所有 x 上都压住这个双峰形态同时还要保持凸性通常会得到一个很“松”的界迭代收敛奇慢。SCA 放弃了这个全局约束只要求替代函数在迭代点 x^k 附近与真实函数有相同函数值、相同梯度并且替代函数自身是强凸的。这个放宽的意义是革命性的你可以用最简单的一阶泰勒展开加一个二次正则项来构造替代函数而不必担心它在远处是否越界。2.2 构造凸近似的三个必要条件在迭代点 x^k为原目标 f(x) 构造凸近似\tilde{f}(x; x^k)时我一般会检查三个条件。第一个是函数值一致\tilde{f}(x^k; x^k) f(x^k)。这个条件保证代理目标与真实目标在迭代点处对齐否则你迭代过程中记录的“目标下降”都不是真实目标的变化。第二个是梯度一致\nabla \tilde{f}(x^k; x^k) \nabla f(x^k)。梯度一旦算错方向整个迭代可能看着在收敛最后落到的却是一个毫无意义的位置。第三个是强凸性\tilde{f}(\cdot; x^k)对 x 是强凸的这保证子问题有唯一解同时迭代步长被自然限制。实际中最常用的模板是一阶展开加二次正则项\tilde{f}(x; x^k) f(x^k) \nabla f(x^k)^T (x - x^k) (1 / 2\tau) \|x - x^k\|^2这里的\tau是正则化系数本质上是步长的倒数。\tau越大二次正则项越小子问题解的移动越激进\tau越小移动越保守。注意如果 f 本身是凸函数这个模板会把它变得更凸加了强凸项子问题依然好解如果 f 本身非凸一阶展开后第一项是仿射函数依然是凸的加上二次项后整个子问题严格凸。这正是 SCA 可以无视原问题凸性的原因。2.3 SCA 迭代为什么能收敛不动点视角收敛性的直观解释可以从不动点条件看。子问题的最优解 x^{k1} 满足一阶最优性条件其中代理函数的梯度在 x^k 处与真实梯度相等。当迭代接近收敛时x^{k1} 与 x^k 的距离趋近于零子问题的一阶条件就逼近原问题的 KKT 条件。只要替代函数与真实函数的误差以\|x - x^k\|^2的阶有界即不会在线性项上产生偏差那么目标函数值序列是单调不增的。实际代码里你只需要在每次迭代时打印真实目标值看到它持续下降且最终平稳基本就可以认为 SCA 在工作如果曲线出现明显上升优先怀疑替代函数梯度算错或者\tau设得太小。这个检查点比任何收敛证明都先用得上。3. 在稀疏回归上用 SCA 凸优化跑通一个真实案例3.1 问题背景用非凸正则逼近 l0 范数稀疏回归的标准做法是用 l1 正则即 LASSO。但 l1 有一个众所周知的缺陷它对大幅系数的惩罚是线性的导致估计结果相比真值系统性偏小。业界在计算资源允许时更倾向使用逼近 l0 的连续非凸正则项。这里我选一个经典形式g(x_i) 1 - exp(-\alpha x_i^2)当 |x_i| 远大于 0 时g 趋近于 1当 x_i 0 时g 0。它像一个“软开关”参数\alpha控制开关的陡峭程度。整个优化问题是min_x \|y - A x\|^2 \lambda \sum_i (1 - exp(-\alpha x_i^2))第一项是凸的第二项非凸。这个问题的规模可以做得很高维非常适合用来演示 SCA 的完整落地过程。3.2 凸子问题怎么推导一阶展开加二次正则在迭代点 x记作 x^k对非凸正则项做一阶泰勒展开g_i(x_i) \approx g_i(x_i^k) g_i(x_i^k)(x_i - x_i^k)其中g_i(x_i^k) 2\alpha x_i^k exp(-\alpha (x_i^k)^2)。把这一项代入原目标再加上二次正则项得到子问题min_x \|y - A x\|^2 \lambda \sum_i [g_i(x_i^k) g_i(x_i^k)(x_i - x_i^k)] (1 / 2\tau) \|x - x^k\|^2这是一个标准二次规划QP二次项来自最小二乘和二次正则线性项来自 g 的导数贡献。CVXPY 可以直接建模底层用 OSQP 求解。注意常数项g_i(x_i^k)对子问题最优解没有影响但如果不把它计入记录的目标值你看到的迭代曲线会与真实目标曲线错位一个常数这对调试不友好所以记录时要加回去。3.3 可运行的 Python 实现下面是完整可跑的代码使用 NumPy 和 CVXPYimport numpy as np import cvxpy as cp def sca_sparse_recovery(A, y, lam0.01, alpha1.0, tau1.0, max_iter30, tol1e-4): n A.shape[1] x np.zeros(n) hist_true [] hist_proxy [] for it in range(max_iter): # 在当前迭代点 x 处计算非凸正则项及其导数 g 1.0 - np.exp(-alpha * np.square(x)) grad_g 2.0 * alpha * x * np.exp(-alpha * np.square(x)) xv cp.Variable(n) residual cp.sum_squares(y - A xv) # 将 g 在 x 处一阶展开加二次正则项保证子问题强凸 linear_penalty lam * cp.sum(g grad_g (xv - x)) proximal_term 0.5 * tau * cp.sum_squares(xv - x) objective cp.Minimize(residual linear_penalty proximal_term) prob cp.Problem(objective) prob.solve(solvercp.OSQP, verboseFalse) x_new xv.value hist_proxy.append(prob.value) true_obj (np.linalg.norm(y - A x_new) ** 2 lam * np.sum(1.0 - np.exp(-alpha * np.square(x_new)))) hist_true.append(true_obj) if np.linalg.norm(x_new - x) tol * (1 np.linalg.norm(x)): x x_new break x x_new return x, hist_true, hist_proxy # 生成仿真数据 m, n, k 80, 200, 10 np.random.seed(0) A np.random.randn(m, n) x_true np.zeros(n) x_true[:k] np.random.randn(k) * 2.0 y A x_true 0.05 * np.random.randn(m) x_hat, hist_true, hist_proxy sca_sparse_recovery( A, y, lam0.01, alpha2.0, tau1.0 ) print(支持集恢复错误:, np.linalg.norm(x_hat - x_true) / np.linalg.norm(x_true))代码逻辑分四步先计算非凸正则项在当前点的函数值和导数然后构造子问题的三项目标分别是残差、一阶展开的惩罚项和二次正则项调用 OSQP 求解最后比较新旧 x 的范数差判断是否收敛。参数说明如下lam正则化系数控制稀疏性强度和偏差之间的平衡lam越大解越稀疏但幅值偏差越大。alpha非凸正则的陡峭度。alpha太小g 接近二次函数退化成岭回归alpha太大g 逼近硬阈值函数子问题附近的梯度变化剧烈需要更小的迭代步长。tau二次正则项系数也是有效步长的倒数。tau偏大会让每次迭代移动距离偏小迭代次数增加tau偏小则可能出现目标值振荡。在数据规模 m80、n200 时用 CVXPY 的 OSQP 求解子问题通常只需几毫秒。打印hist_true会看到目标值单调下降且大多数下降发生在前 5 次迭代。4. SCA 凸优化的参数设定与收敛性控制4.1 核心参数的含义与初始值设定SCA 凸优化在工程上可调的参数并不多但每个参数都直接影响收敛行为。下面是我自己的参数设定表可以直接作为起点参数含义推荐初始值取值偏小取值偏大tau二次正则系数步长倒数与 Lipschitz 常数同量级稀疏回归里取|A|_2^2的倒数迭代步长过大目标值振荡或上升步长过小收敛极慢容易出现“假收敛”alpha非凸正则陡峭度0.1 到 10视数据幅度而定接近凸问题失去非凸优势梯度不平滑数值不稳定lam正则强度取交叉验证后的值解过拟合支持集偏密解全为零信息丢失max_iter最大迭代次数20 到 50不满足收敛判据就退出浪费计算时间tau的初始值可以从问题的 Lipschitz 常数估计。以第 3 节的稀疏回归为例目标函数的二次项是\|y - Ax\|^2其 Hessian 为2 A^T A所以tau可以取tau 1 / (2 * np.linalg.norm(A, 2) ** 2)。通信和信号处理里很多问题可以这样粗估然后按数量级扫一遍tau从估计值的 0.1 倍到 10 倍之间取 3 个值对比收敛曲线即可。4.2 信任域约束与惩罚项的二选一SCA 凸优化里有两种方式限制迭代步长。第一种是惩罚式即目标函数里加(1 / 2\tau) \|x - x^k\|^2这是第 3 节代码里用的方式优点是不需要额外调约束半径缺点是tau的值需要根据数据分布调整。第二种是信任域式在子问题里显式加约束cp.norm(xv - x, 2) R这时二次正则项可以去掉或设得极小。我自己的经验是信任域式更容易调试因为半径R有明确的几何含义惩罚式写起来更省事适合快速原型。两者不要同时大尺度使用否则等效步长会过小。信任域式的子问题建模只需要改两行R 0.1 prob cp.Problem(objective, [cp.norm(xv - x, 2) R])其中R的意义是每次迭代允许 x 离开旧点的最大欧氏距离。R可以按迭代次数衰减比如R R0 / (1 it)这样前期大步探索、后期小步收敛。4.3 加速收敛的两条有效路线第一个是冷启动。用凸松弛解作为 SCA 的初始点比如先用 LASSO 求一个解作为x0可以显著减少 SCA 迭代次数并且避开对称的非凸局部极小值。代码改动非常小from sklearn.linear_model import Lasso x0 Lasso(alpha0.01).fit(A, y).coef_ x x0.copy()在稀疏回归的例子里这种做法通常能把迭代次数从 20 次以上降到 10 次以内。第二个是动量外插。在进入下一次迭代前对迭代点做外插x^k x^k beta * (x^k - x^{k-1})beta取 0.5 到 0.9。动量能加快收敛但对非凸问题可能引入振荡如果看到目标值上升先把beta归零。4.4 收敛判定标准我建议同时使用两个判据而不是只用一个。第一个是目标函数相对变化量|f_{k1} - f_k| / max(1, |f_k|) 1e-4。这个判据在目标函数“地势平坦”时容易过早触发因为目标变化很小但变量还在大范围移动。第二个是变量变化量\|x_{k1} - x_k\|_\infty 1e-3。这个判据防止平地假收敛。两者同时满足才退出对应的 Python 判断是delta_obj abs(hist_true[-1] - hist_true[-2]) / max(1.0, abs(hist_true[-2])) delta_x np.max(np.abs(x_new - x)) if delta_obj 1e-4 and delta_x 1e-3: break需要说明的是这两个阈值的单位依赖于问题本身的数据尺度。如果 y 的量级在 1e6目标绝对值很大max(1, |f_k|)中的 1 可能需要替换成目标量级的千分之一变量同理按1 \|x\|缩放。5. 快速验证 SCA 凸优化实现是否正确的三个检查点5.1 梯度一致性检查SCA 实现里最高频的 bug 是梯度方向写错或漏掉因子 2。可以用数值梯度做一次断言。取一个随机测试点x_test计算解析梯度和数值梯度的最大偏差应小于1e-4eps 1e-6 grad_anal 2.0 * alpha * x_test * np.exp(-alpha * np.square(x_test)) grad_num np.zeros_like(x_test) for i in range(len(x_test)): xp, xm x_test.copy(), x_test.copy() xp[i] eps xm[i] - eps gp 1.0 - np.exp(-alpha * np.square(xp)) gm 1.0 - np.exp(-alpha * np.square(xm)) grad_num[i] (gp[i] - gm[i]) / (2 * eps) print(np.max(np.abs(grad_anal - grad_num)))如果偏差超过1e-4不要往下调试 SCA先把梯度修对。梯度错误会导致目标曲线看上去正常下降但最终解完全错误。5.2 目标函数单调性与回溯对最小化问题SCA 的正确实现应该让真实目标值序列单调不增。如果出现连续两次上升最可能是tau太小。可以加一个回溯逻辑如果hist_true[-1] hist_true[-2]就把tau乘以 2 并重解当前子问题最多重试 5 次。这个策略在很多论文里叫自适应正则化实现成本只有几行却能显著提升算法的鲁棒性。注意回溯只作用于当前迭代点不要修改已经接受的历史迭代点。5.3 用退化场景对齐全局最优把非凸正则项关掉也就是设lam 0此时原问题退化为普通最小二乘最优解有闭式表达式x_opt pinv(A) y。SCA 的每次子问题也退化为带二次正则的最小二乘迭代若干次后应该收敛到与闭式解几乎一致的位置。如果这个测试不过说明 SCA 框架本身有问题和正则项无关。检查时用相对误差例如\|x_hat - x_opt\| / \|x_opt\| 1e-4。这一步通过后再把lam恢复到非零值观察非凸正则带来的稀疏性是否与预期一致。另一个退化方向是让alpha趋近零此时正则项近似二次函数问题退化为岭回归也可以用闭式解做对照。这两个退化测试能快速区分“SCA 实现 bug”和“模型选择问题”算是这个框架里性价比最高的排查手段。本文还有配套的精品资源点击获取