0%

从锥到半定约束:把凸优化模型落到可验证的工作流

优化建模里最容易造成混乱的,往往不是求解器,而是问题类别的名称。线性规划、二次规划、二阶锥规划和半定规划经常同时出现在一篇论文或一段代码中;它们看似是一串由简到繁的缩写,实际上分别对应不同的可行集结构、不同的表达能力,也对应不同的数值风险。

更可靠的工作方式不是先问“该调用哪一个算法”,而是先把问题分成三件事:决策变量是什么,允许的解构成什么集合,目标函数究竟在衡量哪一种代价。只有这三件事明确,才能判断一个二次项能否写成锥约束、一个矩阵不等式是否真的是 LMI,以及所谓“松弛”换来了什么、丢掉了什么。

本文把常见的 LP、QP、QCQP、SOCP 和 SDP 放在同一个框架中讨论,并给出从数学模型到 MATLAB/YALMIP 表达、再到结果验收的完整路径。重点不在记住每种问题的标准型,而在于建立一套可检查的建模习惯。

优化问题先由结构而不是缩写区分

一个一般的最小化问题可以写作:

其中 $x$ 是决策变量,$f_0$ 是目标,约束共同定义可行集。若 $f_0$ 与全部不等式约束都是凸函数,等式约束是仿射函数,则问题是凸优化问题。这里的关键不只是“容易求”,而是局部最优就是全局最优:求解器不必在多个局部谷底之间猜测哪一个更好。

常见类别可从这个角度排列,而不必把它们误解成严格的难度序列:

类别 典型目标与约束 主要结构 建模时最需要确认的事实
LP 仿射目标、仿射不等式与等式 多面体 系数、量纲与边界是否完整
凸 QP 凸二次目标、仿射约束 椭球等值面与多面体的组合 二次项矩阵是否半正定
凸 QCQP 二次目标与二次约束 由凸二次函数定义的集合 每个二次约束是否保持凸性
SOCP 仿射目标、仿射约束、二范数锥约束 二阶锥的仿射逆像 范数右侧是否为非负的仿射上界
SDP 仿射目标、仿射约束、半正定矩阵约束 半正定锥与仿射空间的交 矩阵是否对称,约束是否真是仿射矩阵函数

LP、SOCP 与 SDP 都是锥规划的特例。它们共享“在线性目标下与一个凸锥求交”的思路,只是所使用的锥不同。这个视角比缩写更有用:它解释了为什么许多看似非线性的约束仍可以用内点法处理,也解释了为什么有些非凸问题必须通过松弛、近似或局部算法另行处理。

凸集与锥:可行性的几何语言

集合 $C$ 若满足:对任意 $x,y\in C$ 和任意 $\theta\in[0,1]$,都有

则称 $C$ 为凸集。直观地说,集合内任意两点的连线仍留在集合内。一个集合若同时满足 $x\in K, \alpha\geq0\Rightarrow\alpha x\in K$,称为锥;若还对任意 $x,y\in K$ 有 $x+y\in K$,则称为凸锥。

这两个定义并非抽象装饰。仿射约束 $Ax=b$ 给出一个平面,线性不等式 $Gx\leq h$ 给出半空间;它们的交仍是凸集。范数球、二阶锥和半正定锥也都是凸集或凸锥。只要目标和可行集保持这种结构,求解器便可以利用对偶性、障碍函数和残差控制来给出具有明确含义的数值结果。

二阶锥(也称 Lorentz 锥)定义为

约束

等价于要求 $(Ax+b,c^\mathsf{T}x+d)$ 落在这个锥中。右侧不需要额外再写 $\geq0$:只要不等式可行,范数非负这一事实已经隐含了该条件;不过在模型推导和实现中,仍应检查右侧没有被其他约束推到不可行区域。

二次项出现时,先判断它是否凸

二次规划常写为

若 $P\succeq0$,即 $P$ 为半正定矩阵,目标函数凸;若 $P\succ0$,目标严格凸,在可行解存在时最优点通常唯一。这里的 $\succeq0$ 表示对任意向量 $z$ 均有 $z^\mathsf{T}Pz\geq0$。对称矩阵的这一性质等价于全部特征值非负。

但“式子里有二次项”并不等于“这是一个凸 QP”。考虑二次约束二次规划:

当 $P_0\succeq0$ 且每个 $P_i\succeq0$ 时,它是凸 QCQP。只要某个不等式约束的二次部分不定,或者等式约束含有非仿射二次项,可行集就可能变成非凸集。此时“用凸优化工具箱求一下”并不会凭空恢复全局最优性;工具箱可能拒绝模型,也可能转交给只保证局部性的非凸求解策略。

建模前至少应做两项检查:

  1. 将二次项矩阵显式对称化为 $\frac{1}{2}(P+P^\mathsf{T})$ 后检查特征值或 Cholesky 分解条件;数值上接近零的负特征值应结合缩放与容差判断,不能机械截断。
  2. 区分“凸函数不等式”与“凸函数等式”。$f(x)\leq0$ 在 $f$ 凸时保留凸性;$f(x)=0$ 通常不保留,除非 $f$ 为仿射函数。

这一点对功率控制、波束形成、几何估计和低秩恢复都很重要。许多推导从物理或统计模型出发自然得到 QCQP,但是否可凸化,取决于二次型的符号、变量替换和约束方向,而不是任务名称。

SOCP:把长度、能量和最坏扰动写成锥约束

SOCP 的标准形式可写为:

范数约束的价值在于,它把许多非线性的物理量保留为凸结构,而不是强行线性化。下面三类情形尤其常见。

最小二乘与残差上界

最小二乘问题

与

有相同的最优 $x$。后一种表述把平方目标替换成线性目标加 SOCP 约束,便于同其他锥约束合并。若改用 $\lVert Ax-b\rVert_\infty\leq t$,问题则落入 LP;不同范数对应不同的误差风险模型。

鲁棒线性约束

设约束系数存在欧氏球不确定性:

要求它对所有允许的 $\delta$ 都成立,等价于

这个等价关系来自 Cauchy—Schwarz 不等式,最坏扰动会沿着 $x$ 的方向出现。原本的无限约束族因此压缩成一个二阶锥约束。若不确定集采用 $L_\infty$ 球或多面体,得到的对偶范数和最终模型也会不同;不确定集的形状本身就是风险假设,不能只把 $\rho$ 当作可随意调节的“鲁棒参数”。

二次约束的锥表达

若 $Q\succeq0$ 且 $Q=R^\mathsf{T}R$,则

可写为 $\lVert Rx\rVert_2\leq t$。这说明一部分凸二次约束可以转为 SOCP。但这个变换的前提是二次型可分解为平方范数,且右侧的变量与符号满足锥约束的要求。对不定二次型,形式上凑出“平方”并不能改变其非凸本质。

SDP:用矩阵半正定性表达整体耦合

半定规划的基本形式是

其中 $F_i$ 都是实对称矩阵,$F(x)\succeq0$ 表示其所有特征值非负。它不是逐元素不小于零,而是一种整体谱性质:允许矩阵包含负的非对角元素,也允许部分元素为负,只要求任意方向的二次型 $z^\mathsf{T}F(x)z$ 不为负。

这一点恰好使 SDP 适合描述变量之间的全局耦合。控制中的 Lyapunov 稳定性、相关矩阵、最大特征值界、波束形成的协方差矩阵以及许多组合问题的松弛,都可以写成 LMI 或半定约束。

例如,对称矩阵 $M$ 的最大特征值上界可写为

因此最小化 $t$ 并施加 $tI-M(x)\succeq0$,就是在凸框架下最小化最大特征值。相对于逐个计算特征值后再做外层搜索,LMI 直接把谱界纳入模型,便于与其他线性约束共同处理。

SOCP 可嵌入 SDP。一个常用等价式是

这并不意味着应优先把所有 SOCP 都改写成 SDP。SOCP 的结构更特殊,专用求解器往往更省内存、更高效。表达能力更强的模型类别并不必然是更好的数值选择;应选取能够准确表达问题的最小凸类别。

从非凸 QCQP 到 SDP 松弛:被放松的究竟是什么

很多问题具有 $x^\mathsf{T}Qx$ 这样的二次形式。令

则有 $x^\mathsf{T}Qx=\operatorname{tr}(QX)$。原来的二次项变成了关于 $X$ 的线性函数,但同时引入了三个关键条件:

半正定条件是凸的,秩为一条件通常不是凸的。所谓半定松弛,正是保留线性关系和 $X\succeq0$,删除或放宽 $\operatorname{rank}(X)=1$。得到的 SDP 通常更容易求解,但其最优值是原问题的界,而不自动就是原问题的可实现解。

这一过程至少要区分三种结果:

松弛后的现象 可以得出的结论 仍然需要做的事
$X^\star$ 的数值秩接近 1 松弛可能紧,主特征向量可用于恢复候选 $x$ 将候选代回原约束并比较目标值
$X^\star$ 明显高秩 得到的是界或相关结构,不是直接的原变量解 设计随机化、舍入或问题特定的恢复步骤
松弛不可行 原问题也不可行,前提是松弛确实包含原问题 复核建模与数值缩放,定位冲突约束

“SDR 解出来了”只表示放松后的凸问题求解结束。发布结果时还应报告恢复后的原问题可行性、目标差距和对随机初始化或阈值的敏感性。若恢复步骤含随机化,应至少报告多次运行的分位数,而不是只保留一次最好的样本。

YALMIP 的职责是建模,不是替代求解器

YALMIP 的作用是把接近数学记号的 MATLAB 表达转换为某个求解器接受的标准形式。它负责识别线性、二次、二阶锥、半定或整数结构,完成规范化与接口调用;真正的数值算法由后端求解器执行。把二者混为一谈,会导致两个常见误判:一是将建模错误归咎于求解器,二是以为安装了 YALMIP 就自动具备了所有问题类别的求解能力。

一个清晰的建模流程可固定为五步:

  1. 用 sdpvar、binvar、intvar 等声明变量,并明确矩阵是否必须为对称或满矩阵。
  2. 用一个约束数组累积物理、统计或业务约束;不要把约束分散在多个未检查的临时变量中。
  3. 单独定义标量目标函数,并核对最小化或最大化的方向。
  4. 通过 sdpsettings 选择适配的后端与输出级别,再调用 optimize。
  5. 读取诊断状态、约束残差和原始变量值,最后在工具箱之外复算关键量。

对于声明,P = sdpvar(n,n,'symmetric') 表示实对称决策矩阵;P = sdpvar(n,n,'full') 则允许所有元素独立变化。半正定约束只对对称矩阵有意义,因此若本意是 LMI,就不应误用 full。反过来,若一个一般矩阵被无意声明为对称,模型会悄悄少掉自由度,得到的“最优解”可能对应的是另一个问题。

一个 SOCP 例子:带幅度限制的拟合

下面的例子在拟合 $Ax\approx b$ 时限制决策向量的能量,并把残差的二范数显式写为锥约束:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
% 已知数据:A (m x n), b (m x 1), rho > 0
n = size(A, 2);
x = sdpvar(n, 1);
t = sdpvar(1, 1);

Constraints = [norm(A*x - b, 2) <= t, ...
norm(x, 2) <= rho];
Objective = t;

options = sdpsettings('solver', 'mosek', 'verbose', 1);
diagnostics = optimize(Constraints, Objective, options);

if diagnostics.problem ~= 0
error(yalmiperror(diagnostics.problem));
end

x_hat = value(x);
residual_norm = norm(A*x_hat - b, 2);

这里 $t$ 不是可有可无的辅助变量。它把目标中的非线性残差转为线性目标与一个锥成员关系,使求解器可以识别为 SOCP。求解完成后,residual_norm 应与 value(t) 在数值容差内一致;若明显不一致,优先检查求解状态、变量值是否为空、尺度是否极端,而不是先怀疑 norm 的数学含义。

若环境中没有指定的后端,应该明确更换为本机已安装且支持 SOCP 的求解器,而不是删去 solver 字段后假设自动选择一定符合预期。可用后端、许可证、问题规模和锥类型都会影响实际选择。

一个 SDP 例子:用 Lyapunov 不等式验证稳定性

连续时间线性系统 $\dot z=Az$ 的渐近稳定性可由如下条件保证:存在对称正定矩阵 $P$,使

数值优化不能直接处理严格不等式,通常用明确的安全裕量替代:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
% 已知系统矩阵 A,epsilon 应依据数据尺度设定
n = size(A, 1);
P = sdpvar(n, n, 'symmetric');
epsilon = 1e-6;

Constraints = [P >= epsilon*eye(n), ...
A'*P + P*A <= -epsilon*eye(n)];
Objective = trace(P);

options = sdpsettings('solver', 'mosek', 'verbose', 1);
diagnostics = optimize(Constraints, Objective, options);

if diagnostics.problem ~= 0
error(yalmiperror(diagnostics.problem));
end

P_hat = value(P);
min_eig_P = min(eig(P_hat));
max_eig_L = max(eig(A'*P_hat + P_hat*A));

代码中的矩阵不等式表示半正定序,而不是逐元素比较。P >= epsilon*eye(n) 要求 $P-\epsilon I\succeq0$;第二个约束要求 Lyapunov 矩阵具有负定裕量。trace(P) 只是为了给可行解选择一个温和的尺度,并非稳定性判据本身。实际验收仍需检查 min_eig_P 和 max_eig_L,并把 $ε$ 与系统矩阵的量级放在一起解释。

SeDuMi 标准型与高层建模之间的边界

直接调用 SeDuMi 一类锥求解器时,使用者需要自行将模型装配成锥标准型,例如构造系数矩阵、右端项、目标向量和锥维度信息。线性自由变量、非负正交锥、二阶锥块和半定锥块在内部具有不同的向量化与排列规则。矩阵变量的对称向量化、内积缩放、符号方向和块边界只要有一处不一致,求解器仍可能返回一个数值结果,但它未必对应原始数学模型。

这正是高层建模工具的价值:它把符号表达、标准型转换与求解器接口分开。不过高层接口不会替使用者判断模型是否表达了原意。以下几个边界应始终保留:

工作 高层建模工具可以协助完成 使用者仍须负责
变量与约束翻译 识别多数凸结构,生成标准型 确认变量维度、对称性、单位与约束方向
求解器调用 传递模型、参数和状态码 选择支持目标锥与问题规模的后端
数值诊断 返回问题代码、残差和日志 判断容差、病态性与解是否符合原问题
松弛与恢复 表达半定松弛 验证秩、构造候选解并审查原约束

在旧代码中还会遇到 solvesdp 等历史接口。维护时不宜只做字符串替换:应先确认工具箱版本、后端版本以及模型是否依赖旧默认行为,再迁移到明确的 optimize(Constraints,Objective,options) 调用。这能把约束、目标和求解参数分开,也更容易做回归检查。

数值结果不是证明:至少检查四类残差

凸性保证的是精确数学问题的全局性质;浮点求解面对的却是有限精度、停止准则和缩放后的模型。因此,diagnostics.problem == 0 是必要条件,但不是验收终点。

原始可行性

将求得的 $x^\star$ 代回每类约束,分别计算等式残差、线性不等式最大违反量、锥剩余量和最小特征值。例如 SOCP 约束应检查

LMI 则应检查最小特征值而非每个元素的符号。若变量经过缩放,残差也应在原始物理单位下重新报告。

对偶与最优性

对于凸问题,求解器的原始残差、对偶残差和对偶间隙共同说明停止点质量。不要只看目标值的小数位数。两个候选方案的目标差若远小于残差、数据测量误差或建模误差,就不应据此宣称某一个方案有实际优势。

结构一致性

半定松弛需要看矩阵秩谱;稀疏建模需要看非零项和阈值的关系;整数模型则要检查变量距离整数的偏差。求解器输出的连续近似值是否能够安全转换为业务需要的离散决定,是另一层问题。

对数据与尺度的敏感性

若 $A$ 的部分列量级为 $10^{-8}$、另一些为 $10^8$,即使数学模型没有错误,数值条件也会很差。优先对变量、约束和目标做有物理含义的归一化,再设置容差。盲目把 epsilon 设得更小或把迭代上限调得更大,通常无法修复尺度失衡,反而会把不可靠的解伪装成收敛。

一套更稳妥的建模检查表

在交给求解器之前,建议依次完成以下检查:

  1. 定义决策含义。 每个变量的单位、上下界、是否连续或离散都应在模型旁明确写出。
  2. 逐条确认凸性。 不只看目标函数;二次约束、等式、乘积项、分式和秩约束都需要单独判断。
  3. 选择最小充分结构。 能用 LP 不用 SOCP,能用 SOCP 不无故升级为 SDP;这既改善效率,也减少数值负担。
  4. 保留原问题的验收式。 即使经过变量替换或松弛,也要保留可以直接代回的物理、统计或业务约束。
  5. 区分界与解。 松弛得到的是界、候选或精确解,必须明确标注;不要把 SDP 的可行矩阵直接当成原变量的可行决定。
  6. 记录求解环境。 工具箱版本、求解器、主要选项、缩放方式和容差是结果的一部分,尤其在复现实验时。

结语

从 LP 到 SOCP、再到 SDP,变化的不是缩写的长度,而是可行集能够表达的几何结构。范数约束适合描述长度、能量和最坏扰动;半正定约束适合描述谱界、相关性与全局耦合;半定松弛则以删除秩条件为代价,把难解的二次结构变成可计算的凸界。

YALMIP 能够降低从数学公式到锥标准型的实现门槛,却不能替代结构判断与结果验收。一个可靠的优化工作流,应先确认模型是不是凸的、松弛是否仍保留所需含义,再选择与结构匹配的求解器,并在返回数值解后重新检查可行性、谱性质、间隙和恢复误差。这样,优化工具才是把问题讲清楚的手段,而不是掩盖建模不确定性的黑箱。

------ 本文结束感谢您的阅读 ------