简介面向MATLAB优化学习者和工程应用人员的二次规划QP求解源码包适合希望从代码层面理解约束优化算法实现的读者。资源围绕二次规划的标准形式展开涉及Hessian矩阵、梯度向量及不等式/等式约束的设置并展示quadprog函数与有效集法、路径跟踪法等核心求解思路可直接用于工程优化、经济模型、信号处理等典型场景。压缩包共4个文件全部为.m源文件整体大小2KB文件虽少但模块划分清晰便于逐行研读和二次开发。已有293人浏览学习。通过研读这份源码读者可掌握二次规划建模与算法迭代细节了解自定义优化选项如迭代限制、精度控制对结果的影响为后续扩展和迁移到实际问题提供扎实基础。 年初做投资组合权重优化的时候我又一次被“二次规划”卡住了。需求说简单也简单在总仓位等于1、不允许做空的条件下找一组风险最小的资产权重。说难也难因为打开MATLAB文档准备调用quadprog时发现自己手里的数学式子跟函数要求的输入形式根本对不上。后来我把二次规划的数学模型、quadprog的输入输出以及源代码转换这一步从头到尾捋了一遍才把问题彻底搞明白。这篇笔记适合刚接触MATLAB优化的人也适合已经调过quadprog但总在结果上栽跟头的工程师。核心就三件事二次规划的标准形式是什么、怎么把实际问题写成MATLAB代码、以及求解失败时该从哪里下手。1. 二次规划到底在解什么先把数学形式对上号1.1 那个 0.5 的约定最容易栽跟头二次规划的标准形式长这样min 0.5 * x * H * x f * x s.t. A * x b Aeq * x beq lb x ub其中H是二次项系数矩阵f是一次项系数向量A和b是不等式约束Aeq和beq是等式约束lb和ub是决策变量的上下界。这些都好理解真正坑人的是目标函数开头的那个0.5。quadprog约定的是最小化0.5*x*H*x但很多教材、论文里写的二次规划是min x*H*x f*x没有0.5。如果你从那些资料里直接抄一个H矩阵塞进quadprog等于把目标函数放大了两倍。对于无约束问题最优解位置可能不变但一旦加上不等式或等式约束目标函数缩放就会改变约束与目标之间的权衡结果就偏了。判断方法很简单看H的来源。如果原始问题本身是min xΣx这种形式那转成quadprog时H要写成2*Σ如果原始问题写成了min 0.5*xΣx cxH就是Σ。我建议所有代码里都在注释里写明“这里的H对应的是带0.5的标准形式还是不带0.5的形式”不然过两周回来看代码铁定犯迷糊。1.2 约束矩阵怎么对上号二次规划的应用场景很广投资组合里常见的是权重约束模型预测控制里常见的是状态量、控制量的上下限约束工程优化里常见的是资源总量限制。所有这些约束最终都要落到四种类型里不等式、等式、上下界、以及没有约束。写代码前先把下面这张对应表想清楚问题里的描述需要转换成的形式转换方式所有变量之和等于1Aeq * x beqAeq ones(1, n), beq 1某个线性组合不小于阈值A * x b两边取负号变成 -组合系数*x -阈值变量在0到1之间lb x ublb zeros(n, 1), ub ones(n, 1)没有约束传空数组对应位置写 []新手最容易出问题的就是“大于等于”约束。quadprog只接受小于等于所以遇到mu * w targetRet这种约束必须两边乘负号变成-mu * w -targetRet。这一步丢失了符号后面怎么调都白搭。维度也要反复确认。A必须是m×nb必须是m×1Aeq是p×nbeq是p×1。MATLAB对维度不匹配的报错还算友好但有些版本会把维度错误静默地转成另一种解释导致结果不对还不报错这个后面调试章节会细说。2. 源代码实战从最小二乘到带约束的组合优化2.1 最小二乘问题怎么改写成二次规划很多人不知道最普通的最小二乘拟合其实就是二次规划的特例。考虑min ||C*x - d||^2展开之后||C*x - d||^2 x * (C*C) * x - 2 * d * C * x d * d最后一项d*d是常数不影响最优解可以丢掉。对照quadprog标准形式H 2 * (C * C)f -2 * (C * d)如果原始问题里写的是min 0.5 * ||C*x - d||^2那H和f都要跟着缩放但缩放后最优解x不变。这里可以写一小段代码验证% 最小二乘转二次规划验证 C [1 2; 3 4]; d [1; 1]; H 2 * (C * C); f -2 * (C * d); % 解析解x (C*C)^{-1} * C * d x_analytic -H \ f; % quadprog 求解 options optimoptions(quadprog, Display, off); x_qp quadprog(H, f, [], [], [], [], [], [], [], options); disp([x_analytic, x_qp]);如果这两个结果不一致说明H和f构造有问题而不是求解器有问题。这个例子我每次教别人调quadprog都会先用一遍几秒钟就能定位出90%的模型转换错误。2.2 一个完整可运行的组合权重优化示例回到投资组合的例子。假设有5个资产已知期望收益向量mu和协方差矩阵Sigma想求解风险最小的权重组合约束条件有三个全仓、不允许做空、组合收益不低于目标收益率。% 模拟数据5个资产的期望收益 mu [0.10; 0.12; 0.08; 0.15; 0.11]; % 模拟协方差矩阵确保半正定 rng(0); A_rand randn(5) * 0.2; Sigma A_rand * A_rand eye(5) * 0.01; % 目标最小化 w * Sigma * w % 转成 quadprog 标准形式H 2*Sigma, f 0 n length(mu); H 2 * Sigma; f zeros(n, 1); % 全仓约束sum(w) 1 Aeq ones(1, n); beq 1; % 不允许做空w 0同时每只股票权重不超过1 lb zeros(n, 1); ub ones(n, 1); % 收益约束mu * w targetRet % 转成 A*w b -mu * w -targetRet targetRet 0.10; A -mu; b -targetRet; % 求解 options optimoptions(quadprog, ... Algorithm, interior-point-convex, ... Display, iter); [x, fval, exitflag, output] quadprog(H, f, A, b, Aeq, beq, lb, ub, [], options);注意两点。第一H为什么是2*Sigma因为原始目标函数是w*Sigma*wquadprog标准形式里自带0.5所以必须乘2。第二收益约束为什么要对A和b取负号因为quadprog只能处理小于等于不等式取负号只是移项不是改变约束本质。跑完以后x就是最优权重fval是0.5*x*H*x f*x由于f是零向量fval就等于x*Sigma*x也就是组合方差。想要组合波动率就取sqrt(fval)。这里是最容易误读的地方quadprog返回的fval不是原始目标函数w*Sigma*w本身吗严格说fval是标准形式下的目标函数值但因为H2Σ0.5*xHx xΣx恰好等于原始目标所以直接读fval没毛病。如果你把H写成Σ那fval就不对应原始风险值了。2.3 exitflag与output结果到底靠不靠谱常看exitflag的取值大概分几种情况exitflag含义常见处理方式1找到最优解满足收敛条件正常0达到最大迭代次数可能还没收敛调大MaxIterations或换算法-2找不到可行点约束之间矛盾检查约束是否过紧、是否自相矛盾-3问题无界检查lb/ub以及H是否半正定-4求解过程遇到NaN或Inf检查矩阵中是否有异常值数据是否脏output结构体里值得多看两个字段iterations和firstorderopt。前者告诉你了迭代多少步后者是一阶最优性条件的度量理论上接近0才是好的。如果exitflag0但firstorderopt已经很小说明其实离最优解很近了只是没有达到默认的OptimalityTolerance这时候把容差稍微放松一点就能拿到结果。3. 结果离谱时完整的调试排查链路3.1 第一件事用解析解验证模型而不是怀疑求解器quadprog本身是经过大量测试的绝大多数“结果不对”都是模型到代码的转换问题。所以我的调试顺序很固定先跑一个没有约束的小问题拿解析解和quadprog结果对比。比如min 0.5*xHx fx无约束最优解就是x -H\f。这一步能快速判断H和f是否构造正确。如果解析解对上了再去检查约束部分如果连解析解都对不上问题一定出在H/f转换和后端求解器没有关系。别一上来就调各种tolerance那是在错误方向上浪费时间。3.2 迭代不收敛先看算法再看迭代上限exitflag0对应的是“迭代没跑完但触顶了”。我见过很多人在这种情况下直接把MaxIterations从默认值调大但结果还是老样子。更有效的做法是先把Display设为iter看每一轮的目标函数值和一阶最优性指标有没有在下降。如果每轮都在下降只是下降得慢那调大MaxIterations有意义如果几轮之后目标函数基本不动了说明算法已经到极限这时候应该调整的是约束的数值尺度或者换一个算法。quadprog在凸二次规划上常用的算法是interior-point-convex它默认自带一些预处理。老版本或者特殊场景下也可以选active-set。active-set在中小规模问题上迭代次数通常更稳定但它要求H必须是正定的如果有零特征值或者负特征值会遇到问题。遇到迭代不收敛我可以先跑一个关闭约束的版本再把约束逐步加回来定位到底是哪个约束引起的。3.3 数值尺度不一致一个小权重一把辛酸泪实际项目里的数据不会像官方示例那么整洁。比如协方差矩阵的元素数量级可能只有10^-4而收益率约束是0.1目标函数和约束条件之间尺度差了好几个数量级。interior-point算法对尺度很敏感矩阵条件数太差时即使问题本身有解也可能算出非常离谱的权重。处理办法有两个。一是对数据做归一化比如把收益和协方差都乘同一个倍数让它们的数量级落在1附近。二是设置合理的ConstraintTolerance和OptimalityTolerance。默认的tolerance不一定适配你的问题尺度但也不要随便放大放大太多会得到一个“看着可行其实违反约束”的伪最优解。另外如果协方差矩阵是用历史数据算出来的经常会有很小的负特征值这是数值噪声导致的不是真的非凸。先eig(Sigma)看一眼如果负特征值接近机器精度可以直接用Sigma (Sigma Sigma) / 2 1e-8 * eye(n)做对称化加微扰保证半正定。3.4 无可行解约束自相矛盾的排查顺序exitflag-2是最让人头疼的因为它不是算法不好而是可行域本身就是空的。也就是说你给出的约束条件里根本没有一个点能满足所有限制。排查原则先从最简单的约束组合开始逐步加约束。比如先只保留全仓约束和上下界看有没有解再加上收益约束再看有没有解。哪一步开始报不可行问题就出在哪一步。还可以用一个最小可行性的线性规划来验证% 用 linprog 检查可行域是否为空 f_feas zeros(n, 1); [x_feas, ~, exitflag_feas] linprog(f_feas, A, b, Aeq, beq, lb, ub);如果exitflag_feas不是1说明约束本身就矛盾了。最常见的矛盾是哪几种lb和ub重叠例如要求权重在0到0.5之间又要求权重之和等于1显然不可能。还有就是“大于等于”约束取负号时没转干净比如本该是-mu*w -targetRet写成了mu*w -targetRet方向反了可行域自然可能变成空集。4. 进阶边界什么时候该换其他求解思路4.1 大规模稀疏问题quadprog的极限在哪里quadprog对于几千个变量、几千个约束的稠密问题基本能扛住但规模再往上走或者约束矩阵具有明显稀疏结构就要注意内存和计算效率了。interior-point-convex算法内部需要解一个大型线性方程组H和A如果不是稀疏存储内存会先撑不住。一个有效做法是显式把矩阵转成稀疏类型再传给quadprogH_sparse sparse(H); A_sparse sparse(A);有了稀疏格式quadprog内部会走稀疏线性代数路径能处理的规模会大很多。再往上如果变量上万约束结构又很复杂MATLAB生态里还可以接Gurobi、Mosek这类专业求解器它们在预处理、并行计算方面更强。纯MATLAB环境下也有osqp这样的ADMM求解器但需要自己处理QP标准形式的转换学习成本另算。4.2 非凸二次规划与整数变量quadprog解决不了的部分quadprog默认要求问题至少是凸的更准确地说H必须是对称半正定矩阵。判断方法是eig(H)只要出现一个明显的负特征值就不适合直接扔给quadprog。此时quadprog可能报错也可能给出某个局部解但无法保证是全局最优。对小型非凸问题可以考虑枚举或分支定界对中型问题可以试试fmincon换个非线性求解器但对于二次目标它也只是找局部解需要根据实际场景判断这个局部解接不接受。还有一类更特殊的问题变量要求取整数比如“哪些资产纳入组合用0/1变量表示”。这种带整数约束的二次规划叫MIQPMATLAB自带的intlinprog只能处理线性目标函数quadprog又不接受整数约束两边都搭不上。到了这一步就需要上商业求解器或者专门的启发式算法了。所以遇到MIQP第一反应不应该是找“matlab怎么解”而是先评估问题规模再决定是用Gurobi这类工具还是自己写启发式。最后分享一个我自己的习惯。无论问题看起来多简单我都会先用一个可以手算的小样例跑通再换真实数据。quadprog把矩阵运算和迭代细节封装得很干净也正因如此一旦结果不对问题往往出在“标准形式”的转换上而不是求解器本身。拿这篇文章里的源代码框架换成你自己的C、d、A、b参数基本半小时就能跑出第一版结果。真正的难点从来不是调用函数而是建模时把每一个约束都翻译成quadprog认识的矩阵和向量。本文还有配套的精品资源点击获取