从Waring到DC分解：多项式凸表示的理论与算法实践

1. 从“和”到“差”理解多项式凸表示的核心范式转换在优化、控制乃至机器学习领域我们常常需要处理一个核心问题如何将一个复杂的非线性函数特别是多项式函数表示成更容易处理的形式一个直观的想法是把它写成一堆“好”的东西的“和”。这里的“好”通常指的是凸函数。这就是Waring分解或称“平方和”分解的朴素思想它试图将一个非负的多项式表示为若干个单项式平方的和。如果这个表示存在那么原多项式在全局上都是非负的这为证明全局最优性提供了强有力的工具。我在早期研究多项式优化时曾对Waring分解寄予厚望因为它形式优美且与半定规划SDP有着深刻的联系可以通过Lasserre松弛等框架进行求解。然而现实很快给了我一盆冷水。Waring分解的要求太强了。它要求多项式本身是全局非负的。但实践中绝大多数我们关心的目标函数无论是经济模型中的利润函数还是神经网络中的损失函数抑或是物理系统中的能量函数都并非全局凸更非全局非负。它们往往拥有多个局部极值点呈现出复杂的非凸地形。试图用“和”的范式去逼近或表示这类函数要么根本不存在这样的分解对于不定多项式要么找到的分解阶数高得离谱计算完全不可行。这就像试图只用加法来描绘一幅既有山峰又有山谷的地形图你发现你需要无数个小小的凸起平方项去拟合一个深谷效率极其低下甚至无法实现。这时DC分解Difference-of-Convex凸差分解的范式就显得尤为关键和实用。它的核心思想从“和”转向了“差”将一个可能非常非凸的函数表示为两个凸函数的差。即 f(x) g(x) - h(x)其中 g(x) 和 h(x) 都是凸函数。这个范式转换的魅力在于它一下子放宽了表示的存在性条件。理论上任何连续可微且具有利普希茨连续梯度的函数都可以表示为DC函数。而对于多项式这类光滑函数DC分解总是存在的。这为我们处理非凸多项式优化问题打开了一扇新的大门我们不再纠结于证明全局非负而是接受非凸性并通过分解将其结构显式地暴露出来进而利用凸函数的良好性质设计算法。那么具体到多项式函数如何实现这种从“和”到“差”的表示呢这就是多项式凸表示理论与算法要解决的问题。它不再追求将多项式表示为平方和这属于“和凸”表示而是追求将其表示为“凸差”。这里的“凸”指的是每一项函数是凸函数。对于多项式一个自然的思路是能否将其直接表示为两个凸多项式的差这就引出了本文要探讨的核心如何系统性地找到一个多项式的DC分解特别是保证两个分量都是凸多项式的分解。这不仅仅是理论上的存在性证明更需要构造性的算法以及对其几何意义和算法复杂度的深入理解。接下来我将深入拆解这一理论框架与关键算法。2. DC分解的理论基石为什么凸差是可行的在深入算法之前我们必须夯实理论基础理解为什么DC分解对于多项式不仅是可能的而且是可以系统化求得的。这关系到我们后续所有算法的底气和边界。首先从最一般的函数角度看DC表示之所以广泛存在源于凸分析中的一个深刻事实任何连续函数在足够小的局部都可以被一个凸函数如下方图所控制也可以被一个凹函数如上方图所控制。全局上看通过适当的“拼接”与“调整”总可以构造出全局的凸上界和凸下界其差就是原函数。对于具有利普希茨连续梯度的函数一个经典的构造是利用其梯度的利普希茨常数L。我们可以定义 g(x) (L/2) * ||x||^2那么函数 f(x) g(x) 的 Hessian 矩阵的最小特征值将非负从而成为凸函数。因此f(x) [f(x) (L/2)||x||^2] - [(L/2)||x||^2] 就构成了一个平凡的DC分解。这个分解被称为“基于利普希茨常数的分解”它虽然总是存在但往往非常保守因为L可能很大导致分解后的两个凸函数都非常“陡峭”算法效率不佳。对于多项式这一特殊而重要的函数类我们可以做得比这种平凡分解好得多。多项式的关键在于其Hessian矩阵也是一个多项式矩阵函数。一个二次可微函数是凸的当且仅当其Hessian矩阵处处半正定即所有特征值非负。对于多项式其Hessian矩阵的元素是多项式。因此判断一个多项式是否凸就等价于判断一个多项式矩阵是否全局半正定。这是一个“矩阵值多项式”的非负性问题比标量多项式的非负性问题更复杂但依然处于多项式优化理论的射程之内。那么如何构造两个凸多项式g和h呢一个核心的观察点是利用多项式的齐次展开。任何一个多项式都可以写成其各阶齐次部分的和。其中最高阶齐次部分主部的凸性对整个多项式在无穷远处的行为起主导作用。因此一个构造性思路是首先确保分解后两个多项式的主部都是凸的。这可以通过对原多项式的主部进行适当的“分裂”来实现。例如如果一个多项式的主部是H(x)我们可以寻找一个凸的多项式项P(x)使得 H(x) P(x) 和 P(x) 都是凸多项式的主部。然后再通过调整低阶项来保证全局凸性。更系统化的理论工具是半定规划SDP和平方和SOS规划。我们可以将“寻找凸多项式g”表述为一个约束g的Hessian矩阵 H_g(x) 是一个多项式矩阵且需要满足对于所有x H_g(x) ≽ 0半正定。这个“对所有x”的约束是无限维的。幸运的是对于多项式我们可以借助平方和SOS松弛来将其转化为一个有限维的、可计算的充分条件即要求 H_g(x) 是一个平方和矩阵。这意味着存在一个多项式矩阵V(x)使得 H_g(x) V(x)^T V(x)。这是一个更强的条件SOS矩阵必然是半正定的反之不一定但在许多实际情况下是可行的并且可以精确地转化为一个SDP问题来求解。因此为多项式f寻找DC分解 f g - h 的问题可以转化为一个双凸优化问题我们同时搜索两个多项式g和h使得g和h的Hessian矩阵是SOS矩阵从而保证凸性。f g - h 恒等成立。这本质上是一个在多项式系数空间中的搜索问题可以通过交替优化g和h并利用SDP求解子问题来实现。理论保证了对于给定阶数的多项式只要我们将g和h的阶数设得足够高至少是f的阶数并且允许SOS矩阵的阶数足够高这样的分解总是可以找到。这构成了多项式凸表示算法的理论基础。3. 核心算法拆解如何构造一个多项式的DC分解有了理论保证我们进入最实用的部分具体算法。我将以一个具体的二元三次多项式为例手把手拆解一个基于SOS规划的DC分解算法流程。假设我们有一个非凸多项式f(x, y) x^3 y^3 - 4*x*y x^2 y^2我们的目标是为其找到一个分解 f(x, y) g(x, y) - h(x, y)其中g和h都是凸多项式。3.1 算法框架与问题建模算法的核心思想是迭代优化。我们无法一次性同时找到最优的g和h但可以固定其中一个优化另一个。步骤1初始化。一个简单可靠的初始化是采用前述的“利普希茨”思想但进行多项式特化。我们可以计算f的Hessian矩阵并估计其最大特征值上界ρ这可以通过在单位球上采样或求解一个SDP来粗略估计。然后令初始的凸函数为h0(x, y) (ρ/2) * (x^2 y^2)。这是一个非常简单的凸二次函数。那么初始的g0 f h0。但g0不一定是凸的。所以这只是一个起点我们需要进入迭代。更常见的初始化是直接设h0(x, y) 0然后尝试寻找一个凸的g0来匹配f。这等价于直接寻找f的一个凸上近似如果f本身不是凸的这一步就会失败。因此带有一个简单凸h的初始化更为稳健。步骤2固定h优化g凸上近似子问题。假设我们有一个凸多项式h_k。我们需要找到一个凸多项式g使得g(x, y) f(x, y) h_k(x, y)。注意此时等式右边是已知的。因此问题转化为给定多项式r(x, y) f(x, y) h_k(x, y)寻找一个凸多项式g使得g在某种意义下“接近”r。这可以建模为一个回归问题但带有凸性约束。我们假设g具有预定的阶数例如与f同阶或高一阶其系数为待优化变量c_g。优化目标是最小化g与r之间的差异例如最小化二范数误差 ∫_Ω (g(x) - r(x))^2 dx 在一个定义域Ω上或者简单地要求g与r的对应项系数尽可能接近。约束条件是g的Hessian矩阵 H_g(x) 是半正定的。利用SOS松弛我们将其加强为H_g(x) 是一个平方和矩阵。这便形成了一个半定规划SDP问题决策变量g的系数向量c_g以及用于构造H_g(x)的SOS表示所需的辅助矩阵变量。目标函数Minimize ||c_g - c_r||^2其中c_r是r的系数向量。约束条件(等式约束) H_g(x) 的表达式等于由系数c_g构成的多项式矩阵。(SDP约束) H_g(x) 具有一个预定的平方和表示形式这等价于存在一个半正定矩阵Q使得与多项式的向量化基相乘后满足等式。这是一个线性矩阵不等式LMI。调用SDP求解器如MOSEK, SDPA, 或CVXPYSCS即可解出当前最优的凸多项式g_{k1}。步骤3固定g优化h凸下近似子问题。在得到g_{k1}后我们更新h。关系式为h(x, y) g_{k1}(x, y) - f(x, y)。现在我们需要确保新得到的h是一个凸多项式。但根据等式h已经确定。如果它恰好是凸的那么我们就得到了一个有效的DC分解算法可以终止。然而很多时候这步得到的h未必是凸的。此时我们需要“修正”h使其变为凸的同时尽可能不破坏等式关系。这可以建模为另一个优化问题寻找一个凸多项式h_new使得h_new尽可能接近h_current g_{k1} - f。其建模方式与步骤2完全对称也是一个带凸性SOS约束的回归问题形成另一个SDP。解出h_{k1}后等式f g_{k1} - h_{k1}不再精确成立因为我们对h进行了修正。误差为e h_{k1} - (g_{k1} - f)。这个误差将被带入下一轮迭代。步骤4迭代与收敛。重复步骤2和步骤3。在每一步我们交替地优化g和h使其在满足各自凸性的前提下相互之间尽可能满足差分关系。算法可以设置一个收敛阈值例如当误差e的范数小于某个值或者目标函数g和h与理想值的差异之和的变化很小时停止迭代。输出最终得到的凸多项式g和h即为我们所需的DC分解。3.2 算法实现中的关键细节与技巧多项式的基与阶数选择g和h应该选择多大的阶数理论上不低于f的阶数。实践中通常选择与f同阶。如果同阶找不到再尝试升高一阶。过高的阶数会导致SDP问题规模急剧膨胀决策变量数与基的维度的平方相关计算会非常昂贵。我的经验是从同阶开始尝试。凸性的SOS松弛阶数要求Hessian矩阵是SOS矩阵时我们需要指定SOS多项式的阶数。这个阶数至少是Hessian矩阵元素最高阶数的一半。通常选择最小的可能阶数以控制问题规模。如果最小阶数松弛不可行则需要提高SOS阶数这同样会增加计算负担。定义域的处理上述算法默认追求全局凸性这要求Hessian矩阵在整个R^n上半正定非常严格。在实际优化问题中变量往往有定义域约束如一个紧致的盒子约束。我们可以只要求函数在定义域Ω上是凸的这对应要求Hessian矩阵在Ω上半正定。这可以通过Putinar或Schmüdgen的矩量定理结合局部化的SOS约束来处理从而得到更紧、更容易找到的分解。这是算法实用化的关键一步。利用稀疏性多项式的Hessian矩阵通常很稀疏。在构建SDP时利用这种稀疏性可以极大减少矩阵变量的规模和约束数量。专业的SOS建模工具如SOSTOOLSMATLAB或SumOfSquares.jlJulia可以自动识别并利用这种结构。初始化的艺术好的初始化能大幅减少迭代次数。除了简单的二次函数也可以尝试从f的凸包或凸下半近似出发来构造初始h。例如可以先对f在一些采样点上的值进行凸回归得到一个凸的下近似函数以其作为初始h。4. 从理论到实战DC分解在优化中的应用与挑战找到DC分解f g - h本身不是目的我们的目的是利用这个结构来解决非凸优化问题min f(x)。这就是著名的DC规划D.C. Programming算法框架。其核心算法是DCADifference-of-Convex Algorithm。DCA的思想非常直观利用了凸函数的次梯度性质。对于问题min f(x) g(x) - h(x)其迭代步骤如下初始化选择一个初始点 x^0。次梯度步在当前点 x^k计算凸函数 h 在 x^k 处的一个次梯度 ξ^k ∈ ∂h(x^k)。凸优化步求解下一个迭代点 x^{k1}它是如下凸优化问题的解x^{k1} ∈ argmin_x { g(x) - [h(x^k) ξ^k, x - x^k] }注意方括号内的部分正是函数 h(x) 在点 x^k 处的一个仿射下近似由于h是凸函数其切线全局在下。因此原目标函数 f(x) 被替换为了一个凸的上界函数g(x) - [h(x^k) ξ^k, x - x^k]。最小化这个凸上界函数就得到了下一个迭代点。迭代重复步骤2和3直到收敛。DCA的魅力在于每一步都只需要计算一个次梯度和求解一个凸优化问题。只要凸子问题能高效求解例如g和h都是凸多项式时子问题可能是一个凸的二次或半定规划整个算法就非常可行。并且DCA具有收敛到临界点通常是局部极小值的理论保证。现在将我们前面得到的多项式DC分解与DCA结合就形成了一套处理非凸多项式优化的完整流程建模与分解对目标多项式f运行前述的DC分解算法得到凸多项式g和h。实施DCA应用DCA算法进行优化。局部精炼DCA收敛到一个临界点后可以将其作为初始点送入局部优化器如基于梯度的方法进行精炼。实战中的挑战与心得分解的质量决定算法的命运DCA的收敛速度和最终解的质量极度依赖于DC分解的质量。一个“好”的分解应该使得g和h“尽可能凸”但又不是通过引入巨大的二次项如平凡的利普希茨分解来实现。理想的分解是g和h的曲率都相对平缓这样DCA每一步的凸子问题条件数更好更容易求解且迭代路径更平滑。我们之前的分解算法旨在寻找“最近”的凸函数这在一定程度上能产生质量较好的分解。次梯度的选择与凸子问题的求解对于光滑的凸多项式h其次梯度就是梯度∇h(x)。凸子问题min_x g(x) - ∇h(x^k), x是一个无约束凸多项式最小化问题。当g是二次凸函数时子问题有解析解。当g是更高阶的凸多项式时子问题本身可能仍需用迭代法求解但由于它是凸的且通常形式比原问题简单求解依然比原非凸问题容易得多。我的经验是对于中低阶多项式如4阶以下子问题可以用内点法快速求解对于更高阶的可能需要专门的多项式优化求解器。多个局部解与初始化DCA只能保证收敛到临界点而多项式优化问题通常有多个局部极小值。因此算法对初始点x^0敏感。在实际应用中需要结合多种策略a) 从多个随机初始点运行DCAb) 利用多项式优化本身的全局方法如Lasserre的SOS层次结构提供一个下界并启发式地产生初始点c) 使用我们得到的DC分解结构设计分支定界框架以寻求全局最优解。复杂度瓶颈整个流程的复杂度瓶颈主要在第一步——DC分解。求解大规模的SDP问题尤其是高维、高阶多项式是计算昂贵的。这是将方法应用于大规模问题的首要障碍。研究更高效的、可扩展的DC分解算法例如基于一阶方法近似求解SDP或利用深度学习学习分解是当前的前沿方向。5. 超越多项式DC分解思想的延伸与工具生态虽然本文聚焦于多项式但DC分解的思想和DCA算法具有极大的通用性。任何可以表示为两个凸函数之差的函数都可以纳入这个框架。这包括许多带激活函数的神经网络、一些特殊的超越函数组合等。对于非多项式函数获得其DC分解可能没有系统性的算法但通常可以通过函数本身的特性进行“手工”分解或利用已知的凸凹性质。例如log(1exp(x))(log-sum-exp) 本身是凸函数。x^4 - 2*x^2可以分解为(x^4 a*x^2) - ((a2)*x^2)通过选择合适的a0可以使两项均为凸函数。神经网络中常用的ReLU函数max(0, x)本身就是凸函数而许多损失函数如交叉熵也可以找到其DC表示。在工具层面对于多项式凸表示和DC优化已经有一些优秀的开源库SOSTOOLS (MATLAB)最老牌、功能最全面的SOS优化工具箱可以方便地建模多项式凸性SOS矩阵约束用于实现第3节的分解算法。SumOfSquares.jl (Julia)基于JuMP建模语言与Julia的优化生态系统无缝集成是目前最活跃、性能也相当出色的工具。它支持利用多核并行和稀疏性能处理规模更大的问题。CVXPY (Python)结合cvxpy的cp.sum_of_squares特性也可以进行SOS规划但生态相对Julia版本稍弱。DCA算法实现DCA的核心步骤相对简单可以自己编码实现。也有一些专门的工具箱如虽然不限于多项式dca_algorithm等。在我自己的研究实践中我通常使用SumOfSquares.jl进行多项式DC分解的建模和求解然后将得到的g和h的系数导出再用JuMP或Convex.jl编写DCA的迭代循环。Julia语言在科学计算和优化方面的性能优势使得这套流程可以处理变量数在10个左右、阶数在4阶以下的中等问题。最后必须提的注意事项是DC分解和DCA是一个强大的框架但它不是万能的银弹。它的效力介于纯粹的局部方法如梯度下降可能陷入很差的局部点和全局方法如Lasserre层次结构计算成本极高之间。对于具有特殊结构如稀疏性、低有效维数的中等规模非凸多项式优化问题这套基于凸表示的方法往往能在可接受的时间内找到高质量的解甚至通过全局化策略证明全局最优性。理解其理论边界熟练运用现有的工具链并根据具体问题设计分解策略是发挥其最大效力的关键。这需要不断的实践和试错但一旦掌握它就成为了你处理复杂非线性问题工具箱中一件非常犀利的武器。