简介:面向需要掌握凸优化理论并在MATLAB中动手实践的工程师、科研人员与学生,这份代码包以凸优化算法实例为核心,同时提供Python和Julia实现,便于跨语言对照解决实际建模与求解中的常见问题;适用于课程作业、毕业设计及实际工程参考,能够显著降低凸优化算法的入门门槛。凸优化在机器学习、信号处理、控制系统设计等领域应用广泛,资源内容紧贴这些场景,从基础建模到算法实现逐步展开,配合示例与测试数据,帮助读者快速搭起理论到实践的桥梁。压缩包共199个文件,包含103个.m脚本、46个.py脚本和41个.jl脚本,以及少量.png示意图、.mat数据和.txt说明,整体仅595KB,轻量而聚焦;不同语言文件适合对比各自语法与求解器调用方式。目前已有441人学习/下载,是经过一定验证的入门与进阶资料。内容覆盖基础凸优化问题建模、梯度下降法、内点法、ADMM等常见算法的实现,并配套最小二乘、支持向量机等典型应用实例;从基础例子到测试验证脚本均有涉及,通过运行代码,读者可以直观理解不同算法的收敛行为与适用场景,同时提升MATLAB、Python、Julia的多语言编程能力;目录按功能模块划分,便于按需查阅和快速上手。
1. 拿到 Code.zip 之后,先别急着 run
解压一个名为 Code.zip 的凸优化代码包,通常看到的是这样几种东西:一叠.m脚本、几个被反复调用的函数文件、一份写得像论文附录的 README,以及藏在某个子目录里的cvx_setup.m或运行日志。这个标题把 convex optimization 和 MATLAB 放在一起,本质上是在划一条边界:MATLAB 里能解优化问题的工具很多,但凸优化值得单独拿出来讲,是因为这一类问题的局部最优解就是全局最优解,求解器可以放心地走内点法、梯度下降这类局部迭代路线,而不必担心陷在鞍点或局部极小值里出不来。所以这篇文章直接从建模、选求解器、调参数三层往下讲,适合手里已经有一堆数据、想把问题转成标准凸形式、并且希望在 MATLAB 里跑出稳定结果的人。我默认你已经装了 Optimization Toolbox,也知道linprog和quadprog的大概用途,下面只讲怎么写对、怎么调稳,以及拿到陌生代码包时怎么最快看懂它在干什么。
2. 在 MATLAB 里建立凸优化模型的三个前置动作
2.1 先判断问题是不是凸的:用hessian做数值探针
很多人在 MATLAB 里把fmincon当作万能求解器,丢进去一个目标函数就开始迭代。这样在凸问题上通常也能收敛,但局面会变得不可控:你不知道求解器走了多少步、不知道结果是否稳定、也不知道如果再改一个约束会不会直接发散。所以在写任何调用代码之前,我习惯先做一个十几行的数值探针,检查目标函数和约束函数在可行域内的凸性。这个探针不追求数学严格,只用来快速排除明显非凸的情况。
% check_convexity.m % 用有限差分计算目标函数在若干随机点处的 Hessian 特征值 rng(42); x0 = randn(10, 1); h = 1e-6; H = zeros(10, 10); for i = 1:10 for j = 1:10 e_i = zeros(10,1); e_i(i) = 1; e_j = zeros(10,1); e_j(j) = 1; H(i,j) = (f(x0 + h*e_i + h*e_j) - f(x0 + h*e_i - h*e_j) ... - f(x0 - h*e_i + h*e_j) + f(x0 - h*e_i - h*e_j)) / (4*h^2); end end eigs_H = eig((H + H')/2); fprintf('最小特征值: %.6e, 最大特征值: %.6e\n', min(eigs_H), max(eigs_H));这段代码用中心差分离散求 Hessian,在多个随机点上取特征值。如果所有特征值都大于负的容差(比如-1e-8),说明目标函数在局部是凸的。注意这里的关键点是中心差分需要调用 4 次目标函数,维度高了之后成本会翻倍,但作为建模前的快速探针完全够用。特征值的结果如果出现明显的负值,就要回到公式层面检查是哪里引入了非凸项。
常见误用是用fminunc的输出去反推问题是否凸,这是不成立的:梯度为零只说明找到了一个平稳点,平稳点可能是鞍点,也可能是极大值点。数值探针配合eig是直接看曲率,比看收敛曲线可靠得多。
2.2 建模工具的选型:linprog、quadprog、fmincon、CVX、YALMIP 怎么选
确定了问题是凸的,接下来就是选建模语言。这个问题在 MATLAB 生态里通常有三条路线。第一条是直接用 Optimization Toolbox 的底层 API,linprog对应线性规划,quadprog对应二次规划,fmincon对应一般非线性规划。这条路的好处是不需要装任何第三方包,受 MathWorks 官方维护,调用方式稳定,适合生产环境。第二条是 CVX,用的是类似 cvx_begin 的声明式建模语法,把约束和目标函数以几乎数学原样的形式写出来,内部自动做对偶转换和锥规划映射,不需要手推标准形。第三条是 YALMIP,本身不自带求解器,它提供一个统一建模接口,后端可以挂 Gurobi、Mosek、SeDuMi 等求解器,灵活性最高。
从我的使用习惯来说,建模阶段会用 CVX 快速验证模型正确性,等模型定稿之后再用底层 API 重写一遍做性能优化。原因是 CVX 的可读性高,约束写错了一眼能看出来,但它在某些大规模稀疏问题上会引入额外的预处理开销;quadprog在稠密中小规模问题上通常更快。下面这组对比表格可以快速决策。
| 场景 | 推荐工具 | 原因 |
|---|---|---|
| 目标函数和约束都是线性 | linprog | 自带单纯形法和内点法,零配置 |
| 二次目标 + 线性约束 | quadprog | 支持稀疏矩阵,内存占用可控 |
| 带非光滑项(如 L1 范数) | CVX 或 YALMIP | 自动做变量分裂和锥规划转换 |
| 需要频繁更换求解器做对比 | YALMIP | 后端求解器可插拔,方便做交叉验证 |
| 问题规模超过 1e5 个变量 | quadprog或 Mosek | CVX 的预处理可能成为性能瓶颈 |
选工具的关键不只在语法,而在于你打算怎么修改模型。如果约束条件在一个项目中会被反复调整,CVX 的声明式语法能让你在五分钟内完成从“加一个约束”到“重新求解”的整个循环;如果你已经知道模型结构不会再变,quadprog最直接的收益就是少一层封装,启动时间和每次迭代的开销都会更低。
2.3 从数学形式到 MATLAB 表达:一个 L1 范数最小化的转换示例
真实项目中很少会遇到恰好就是A*x <= b这种现成形式的问题。更多的情况是需要把数学上的绝对值、范数、最大最小操作改写成 MATLAB 能识别的约束。以最经典的 L1 范数最小化为例:min ||x||_1看起来不能用linprog直接解,但通过引入辅助变量t,可以等价转换成线性规划:min sum(t),约束是-t <= x <= t。这个转换的好处是把不可微的目标函数变成了光滑的线性目标加线性约束,求解效率会提高一个数量级。
% l1_minimization.m % 用 linprog 求解 min ||x||_1, subject to A*x = b n = 10; m = 6; A = randn(m, n); x_true = zeros(n, 1); x_true([2 5 8]) = [1.5 -2 3]; % 稀疏真值 b = A * x_true; % 定义 linprog 的系数矩阵,x 扩增为 [x; t] f = [zeros(n,1); ones(n,1)]; % 目标只对 t 求和 Aeq = [A, zeros(m, n)]; % A*x = b beq = b; Aineq = [eye(n), -eye(n); -eye(n), -eye(n)]; % t >= x 和 t >= -x bineq = zeros(2*n, 1); lb = [-inf(n,1); zeros(n,1)]; % t 必须非负 [x_opt, fval, exitflag] = linprog(f, Aineq, bineq, Aeq, beq, lb); x_opt = x_opt(1:n);这里的核心逻辑是把|x_i|拆成两个不等式:x_i - t_i <= 0和-x_i - t_i <= 0。因为t_i被目标函数最小化,在最优解处必然有t_i = |x_i|。这种转换在有linprog时特别实用,能规避对不可微点的特殊处理。参数上要注意,Aineq矩阵的分块顺序不能写反,否则约束会被错误地松绑;lb里只对t部分设置非负下界,x本身允许正负。用稀疏矩阵构造这几个分块矩阵时,可以用sparse函数节省内存,尤其当n到1e4以上时,直观的全矩阵写法会让内存先于迭代崩溃。
3. 凸优化核心求解器在 MATLAB 里的标准用法
3.1 用linprog跑通线性规划的最小完整例子
线性规划是凸优化里最基础也最常用的一类。调用linprog的完整语法是[x, fval, exitflag, output, lambda] = linprog(f, A, b, Aeq, beq, lb, ub, options)。这里f是目标函数的系数向量,A和b定义不等式约束,Aeq和beq定义等式约束,lb和ub是决策变量的边界。一个常被忽略的细节是:如果问题里没有不等式约束,必须传A = []和b = [],而不是直接跳过这两个参数。
% lp_basic.m % min f'*x, s.t. A*x <= b, Aeq*x = beq f = [-3; -2]; % 最大化 3x1+2x2 的等价写法 A = [1, 1; 2, 1; -1, 0; 0, -1]; b = [10; 16; 0; 0]; options = optimoptions('linprog', 'Algorithm', 'dual-simplex', 'Display', 'iter'); [x, fval, exitflag, output] = linprog(f, A, b, [], [], [], [], options);dual-simplex算法在有热启动需求时会比默认的内点法更有优势,因为单纯形法可以复用一个可行基解做增量迭代;内点法每次从零开始,适合大规模稀疏问题。Display设为'iter'后,命令行会输出每次迭代的对偶可行性和互补间隙,这是判断算法是否健康最直接的窗口。如果输出的exitflag不是 1,不要只看错误码,要配合output.dualfeasibility和output.iterations一起判断是数值问题还是模型不可行。
3.2 用quadprog实现带约束的最小二乘估计
带线性约束的最小二乘是凸优化里最常被“线性化误用”的问题。很多人遇到min ||Ax - b||^2就直接用A\b,一旦加上Cx = d这样的硬约束,就不知道该怎么办了。正确方案是用quadprog:目标函数的二次项系数矩阵是A'*A,一次项是-2*A'*b,约束照写。关键是要确认A'*A是半正定的,这在大多数实际数据下成立,但如果A是列亏秩的,quadprog可能会报错说 Hessian 不正定。
% qp_constrained_LS.m % min ||Ax-b||^2 subject to Cx = d A = randn(20, 8); b = randn(20, 1); C = [1, -1, 0, 0, 0, 0, 0, 0; 0, 0, 1, -1, 0, 0, 0, 0]; d = [0; 0]; H = A' * A; f = -2 * A' * b; options = optimoptions('quadprog', 'Algorithm', 'interior-point-convex', ... 'OptimalityTolerance', 1e-10); [x, fval, exitflag, output] = quadprog(H, f, [], [], C, d, [], [], [], options); fprintf('约束残差: %.3e\n', norm(C*x - d));interior-point-convex算法对半正定 Hessian 的预处理做了专门优化,比默认的算法在大型问题上更稳定。OptimalityTolerance直接控制一阶最优性条件的满足程度,如果你发现约束残差总是在1e-6附近徘徊,通常不是迭代没收敛,而是这个容差给大了。需要注意的是quadprog只接受凸二次规划,如果 Hessian 有负特征值,无论容差怎么调都不会得到全局最优解,需要回到建模阶段修正。
3.3 用 CVX 写 LASSO:建模可读性优先
虽然quadprog可以解决 L2 正则化的最小二乘,但在 L1 正则化(即 LASSO)问题上,quadprog需要引入额外的扩展维度,代码可读性会显著下降。CVX 的优势在这一类问题上体现得最明显,因为它的声明式语法允许你直接把数学式写出来。
% lasso_cvx.m rng(7); X = randn(100, 20); w_true = zeros(20,1); w_true([1 5 10]) = [2 -3 1]; y = X*w_true + 0.1*randn(100,1); lambda = 0.1; cvx_begin quiet variable w(20) minimize( 0.5 * sum_square(X*w - y) + lambda * norm(w, 1) ) cvx_end % 对比: 不考虑稀疏性的最小二乘 w_ls = X \ y; fprintf('CVX 求解器状态: %s\n', cvx_status); fprintf('w 的非零个数: %d (真实为 3)\n', nnz(abs(w) > 1e-4));sum_square等价于sum((X*w-y).^2),norm(w, 1)就是 L1 范数。CVX 内部会把这个问题自动转换为锥规划,后端默认调用 SDPT3 或 SeDuMi,不需要关心底层细节。关于lambda的选择,网格搜索一个对数等距序列通常能找到统计意义上更合理的正则化强度;值得注意的是cvx_status要确认是Solved,如果显示Inaccurate/Solved,通常说明数值条件已经接近求解器能力边界,需要回头查看数据的尺度是否严重不平衡。
3.4fmincon在凸优化中的定位与局限
许多工程人员习惯用fmincon处理所有非线性优化问题。在凸优化的语境下,fmincon是可以胜任的,特别是在约束中包含非线性函数的情况下。但它的默认算法interior-point处理非光滑函数时会止步于不可微点,所以当目标函数包含 L1 范数、绝对值、分段线性项时,更稳妥的路线是先做等价变换(如前面章节提到的线性化),而不是直接在fmincon里传一个包含abs函数的句柄。
% fmincon_demo.m % min (x1-2)^2 + (x2-1)^2, subject to x1^2 + x2^2 <= 1 fun = @(x) (x(1)-2)^2 + (x(2)-1)^2; nonlcon = @(x) deal(x(1)^2 + x(2)^2 - 1, []); % 不等式约束 c(x) <= 0 opts = optimoptions('fmincon', 'SpecifyConstraintGradient', false, ... 'Display', 'final', 'Algorithm', 'sqp'); x0 = [0; 0]; [x, fval, exitflag, output] = fmincon(fun, x0, [], [], [], [], [], [], nonlcon, opts);sqp算法在处理中等规模非线性约束时通常比内点法收敛更快,且每一步都保持可行解。SpecifyConstraintGradient设为false是让求解器做数值差分,但注意这会把约束梯度的误差带进 KKT 条件的判定里,所以如果条件数不好,建议手写gradient函数并置为true。这类问题并不是严格意义上的凸问题(约束集是凸集,目标是凸函数),所以fmincon的收敛结论是可靠的,但只有当目标函数是凸函数时,局部解才等于全局解。
4. 凸优化求解器的参数选择与排错方法
4.1 解读 output 结构体:迭代次数与一阶梯度的门道
linprog、quadprog、fmincon的第四个输出参量是output结构体,里面包含的字段在不同求解器中略有不同,但iterations、algorithm、firstorderopt、constrviolation这四样是通用的。firstorderopt是求解器判断是否达到最优的标准,它本质上是对 KKT 条件的残差度量。当这个值接近你设置的OptimalityTolerance时,求解器会判定收敛,而不论迭代步数还剩下多少。所以如果你发现求解器很少迭代几步就停了,先看firstorderopt是否符合预期,而不是急着加迭代次数上限。
% inspect_output.m % 以 quadprog 为例,检查收敛指标 [~, ~, ~, output] = quadprog(H, f, [], [], C, d); disp(output); fprintf('firstorderopt = %.6e\n', output.firstorderopt); fprintf('constraint violation = %.6e\n', output.constrviolation); if output.firstorderopt > 1e-6 warning('KKT 残差偏大,结果可能不是精确最优解'); endconstrviolation是约束违反量的指标。如果它大于1e-8,说明最终解在数值上突破了某个约束边界,可能是由于模型矩阵条件数过大,或是ConstraintTolerance设置得过松。这里的定位方法是:把约束矩阵做一次奇异值分解,查看最小奇异值是否接近机器精度。如果是,那就不是求解器的问题,而是约束本身线性相关。
4.2 内点法与有效集法的适用范围:Algorithm 参数怎么选
现在的问题是:quadprog支持三种算法,interior-point-convex、trust-region-reflective,以及active-set。不同算法面对同一问题的稳定性差异巨大。interior-point-convex是大规模问题的首选,它利用稀疏结构,迭代次数与问题维度近似对数关系,适合变量数量在万级以上的场景。但它的缺点是解可能不是严格意义上的基本可行解,边界上的值可能是“擦着边”的。trust-region-reflective要求目标函数必须是凸的,且只支持边界约束和线性等式约束,不适合带不等式的问题。active-set在小规模稠密问题上表现优秀,因为它显式维护工作集,每一步求解的子问题规模小,且天然返回极点解。
% compare_algorithms.m opts1 = optimoptions('quadprog', 'Algorithm', 'interior-point-convex'); opts2 = optimoptions('quadprog', 'Algorithm', 'active-set'); [~, ~, ~, out1] = quadprog(H, f, [], [], C, d, [], [], opts1); [~, ~, ~, out2] = quadprog(H, f, [], [], C, d, [], [], opts2); fprintf('interior-point: %d 次迭代\n', out1.iterations); fprintf('active-set: %d 次迭代\n', out2.iterations);从实际工程经验看,当变量规模小于 500 且约束数量不多时,active-set的数值精度往往更好。当变量规模超过 5000 或约束矩阵是稀疏的时候,interior-point-convex是默认选择。如果有人问“为什么两种算法的结果差一点”,通常不是算法实现的问题,而是两种算法收敛到了不同精度的解。
4.3 五类常见报错与排查顺序
第一类是最常见的错误提示"Solver stopped prematurely",在 CVX 中常伴随cvx_status为Failed。排查方向是先看问题是否是数值病态的:把目标函数的系数都归一化到相近尺度,再重新求解。很多这类问题都与系数矩阵量级差异过大有关。
% diagnose_scaling.m % 检查约束矩阵的尺度差异 col_norms = sqrt(sum(Aineq.^2, 1)); fprintf('列范数范围: %.3e ~ %.3e\n', min(col_norms), max(col_norms)); if max(col_norms) / min(col_norms) > 1e4 warning('矩阵尺度不平衡,建议先做行/列归一化'); end第二类是"Hessian is not positive semidefinite"。这个错误的意思是模型本身被写成了非凸的形式,而不是求解器有问题。解决办法是检查H的特征值是否出现负值,若出现则说明问题的二次型本质上是非凸的。第三类是"Number of iterations exceeded options.MaxIterations"。这往往不是步数不够,而是目标函数梯度与约束梯度之间存在尺度冲突,先检查梯度的数值量级。第四类是exitflag = 0表示迭代次数达到上限,但output.firstorderopt仍然很大,这时应当做变量标准化。第五类是内存错误,常见于linprog在稠密矩阵配合内点法时出现,解决办法是检查矩阵是否已用sparse存储。
4.4 用optimoptions控制内部迭代参数
optimoptions是控制求解器内部参数的最直接手段。常用的几个参数分别是OptimalityTolerance(一阶最优性容差)、ConstraintTolerance(约束可行性容差)、StepTolerance(步长最小阈值)、MaxIterations(迭代次数上限,默认通常为 2000)、FunctionTolerance(目标函数值变化阈值)。需要特别注意的是,OptimalityTolerance和ConstraintTolerance不是越小越好。太紧的容差会让内点法在某一步之后进入缓慢的“磨”阶段,反复迭代却几乎不前进。经验上,对大多数工程问题,将这两个容差设置在1e-8到1e-10之间已足够,过小的容差只会让求解器浪费计算资源在测试稀疏因式分解精度上。
% set_tolerances.m options = optimoptions('quadprog', ... 'OptimalityTolerance', 1e-10, ... 'ConstraintTolerance', 1e-8, ... 'MaxIterations', 3000, ... 'Display', 'iter');如果在这组参数下仍无法收敛,通常意味着模型结构需要调整,而不是继续放宽容差。例如增加一个正则化项lambda * norm(x)^2可以改善病态条件下的收敛性,这在工程上是常见的手段,即使原始模型并不需要正则化,加入一个极小的lambda(比如1e-8)也能显著提高数值稳定性。
5. 快速拆解陌生凸优化代码包的五个位置与数值验证技巧
5.1 先扫描这五处,能省下一小时
拿到一个陌生代码包,不要从头到尾读每个.m文件。我一般按顺序检查五个位置:主入口脚本、目标函数的定义、约束函数的定义、求解器选择(linprog/quadprog/cvx/fmincon)、最后一个位置的配置参数。入口脚本通常决定了整个流程的顺序;目标函数定义里如果看到norm、abs、sum_square,说明这是用 CVX 或 YALMIP 建模的;约束函数里如果出现非线性的function句柄,说明它走的是fmincon路线。除了注意版权信息,更重要的是寻找一个通向模型标准的“桥”:看看解压目录下是否有cvx_setup.m,如果有,整个包大概率是围着 CVX 写的。
% scan_package.m % 在当前目录及一级子目录中搜索关键调用 files = dir('**/*.m'); keywords = {'cvx_begin', 'linprog', 'quadprog', 'fmincon', 'yalmip', 'optimproblem'}; for k = 1:length(keywords) count = 0; for i = 1:length(files) text = fileread(fullfile(files(i).folder, files(i).name)); if contains(text, keywords{k}) count = count + 1; end end if count > 0 fprintf('%s: %d 个文件引用\n', keywords{k}, count); end end这种方式可以快速知道一个代码包的技术路线。看到optimproblem,说明用的是 Optimization Toolbox 的面向对象接口,这比较容易转回quadprog或linprog。看到cvx_begin则说明写法是声明式,要改的话必须安装 CVX。
5.2 用数值梯度校验解析梯度
凸优化代码包里经常出现手写解析梯度的函数,这是性能的关键,也是 bug 的高发区。写错梯度的现象很常见,最稳妥的验证方式是数值梯度交叉检验:找一个随机点,对比解析梯度与中心差分离散求导的结果,二者应相对一致到1e-5的量级。
% check_gradient.m x_test = randn(10, 1); analytic_grad = grad_fun(x_test); % 代码包提供的解析梯度 h = 1e-6; numeric_grad = zeros(size(x_test)); for i = 1:length(x_test) e = zeros(size(x_test)); e(i) = 1; numeric_grad(i) = (f_fun(x_test + h*e) - f_fun(x_test - h*e)) / (2*h); end rel_err = norm(analytic_grad - numeric_grad) / max(1, norm(numeric_grad)); fprintf('梯度相对误差: %.3e\n', rel_err); if rel_err > 1e-4 error('解析梯度与数值梯度不一致,检查梯度推导'); end如果解析梯度的表达式推导是正确的,相对误差应该在1e-6到1e-8的范围内。误报通常来自有限差分步长h的选择,1e-6是一个常见的安全区间,太小会让舍入误差占主导,太大又会引入截断误差。若相对误差在1e-4量级,不建议继续使用这个函数,因为fmincon在梯度残差较大时,KKT 条件的判定会失灵,导致提前终止或收敛到错误方向。
5.3 一个能改写的终极验证技巧:无约束极值核对
拿到一个凸优化代码包后,最快的正确性验证方法不是直接跑很大规模的数据,而是把问题缩小到二维或三维,在图上画出目标函数的等值线,同时画出约束边界,然后把求解器的解标在图上。如果数值解明显位于等值线的最低处并与约束边界相切,那么代码包的基本逻辑是通的。这个操作在 MATLAB 里用fcontour加plot就能完成,而且对于自查代码模型错误很有用,尤其是约束方向写反、目标函数系数正负号写反这类隐蔽错误。
然后针对高维数据,另加一条验证路径:构造一个已知全局最优解的问题,再丢进求解器,比较输出解和真实解之间的误差。这条思路简单却非常有效——很多代码包在真实数据上跑了很久才发现问题,而一个精心构造的合成数据实验能在几秒内定位算法或建模的实现缺陷。用randn生成稀疏权重做 LASSO 的对照实验时,注意要固定随机种子,对比w非零位置是否与真实非零位置一致,这才是真正有意义的检验。
本文还有配套的精品资源,点击获取