首页/资讯中心/详情

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

从Waring到DC分解：多项式凸表示的理论与算法实践
📅 发布时间：2026/6/27 0:42:26

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，使得：

  1. g和h的Hessian矩阵是SOS矩阵（从而保证凸性）。
  2. 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的系数向量。
  • 约束条件：
    1. (等式约束) H_g(x) 的表达式等于由系数c_g构成的多项式矩阵。
    2. (SDP约束) H_g(x) 具有一个预定的平方和表示形式，这等价于存在一个半正定矩阵Q，使得与多项式的向量化基相乘后满足等式。这是一个线性矩阵不等式（LMI）。

调用SDP求解器（如MOSEK, SDPA, 或CVXPY+SCS）即可解出当前最优的凸多项式g_{k+1}。

步骤3：固定g，优化h（凸下近似子问题）。在得到g_{k+1}后，我们更新h。关系式为h(x, y) = g_{k+1}(x, y) - f(x, y)。现在，我们需要确保新得到的h是一个凸多项式。但根据等式，h已经确定。如果它恰好是凸的，那么我们就得到了一个有效的DC分解，算法可以终止。

然而，很多时候这步得到的h未必是凸的。此时，我们需要“修正”h，使其变为凸的，同时尽可能不破坏等式关系。这可以建模为另一个优化问题：寻找一个凸多项式h_new，使得h_new尽可能接近h_current = g_{k+1} - f。其建模方式与步骤2完全对称，也是一个带凸性（SOS）约束的回归问题，形成另一个SDP。

解出h_{k+1}后，等式f = g_{k+1} - h_{k+1}不再精确成立，因为我们对h进行了修正。误差为e = h_{k+1} - (g_{k+1} - f)。这个误差将被带入下一轮迭代。

步骤4：迭代与收敛。重复步骤2和步骤3。在每一步，我们交替地优化g和h，使其在满足各自凸性的前提下，相互之间尽可能满足差分关系。算法可以设置一个收敛阈值，例如当误差e的范数小于某个值，或者目标函数（g和h与理想值的差异之和）的变化很小时，停止迭代。

输出最终得到的凸多项式g和h，即为我们所需的DC分解。

3.2 算法实现中的关键细节与技巧

  1. 多项式的基与阶数选择：g和h应该选择多大的阶数？理论上，不低于f的阶数。实践中，通常选择与f同阶。如果同阶找不到，再尝试升高一阶。过高的阶数会导致SDP问题规模急剧膨胀（决策变量数与基的维度的平方相关），计算会非常昂贵。我的经验是从同阶开始尝试。

  2. 凸性的SOS松弛阶数：要求Hessian矩阵是SOS矩阵时，我们需要指定SOS多项式的阶数。这个阶数至少是Hessian矩阵元素最高阶数的一半。通常选择最小的可能阶数以控制问题规模。如果最小阶数松弛不可行，则需要提高SOS阶数，这同样会增加计算负担。

  3. 定义域的处理：上述算法默认追求全局凸性，这要求Hessian矩阵在整个R^n上半正定，非常严格。在实际优化问题中，变量往往有定义域约束（如一个紧致的盒子约束）。我们可以只要求函数在定义域Ω上是凸的，这对应要求Hessian矩阵在Ω上半正定。这可以通过Putinar或Schmüdgen的矩量定理，结合局部化的SOS约束来处理，从而得到更紧、更容易找到的分解。这是算法实用化的关键一步。

  4. 利用稀疏性：多项式的Hessian矩阵通常很稀疏。在构建SDP时，利用这种稀疏性可以极大减少矩阵变量的规模和约束数量。专业的SOS建模工具如SOSTOOLS（MATLAB）或SumOfSquares.jl（Julia）可以自动识别并利用这种结构。

  5. 初始化的艺术：好的初始化能大幅减少迭代次数。除了简单的二次函数，也可以尝试从f的凸包或凸下半近似出发来构造初始h。例如，可以先对f在一些采样点上的值进行凸回归，得到一个凸的下近似函数，以其作为初始h。

4. 从理论到实战：DC分解在优化中的应用与挑战

找到DC分解f = g - h本身不是目的，我们的目的是利用这个结构来解决非凸优化问题min f(x)。这就是著名的DC规划（D.C. Programming）算法框架。

其核心算法是DCA（Difference-of-Convex Algorithm）。DCA的思想非常直观，利用了凸函数的次梯度性质。对于问题min f(x) = g(x) - h(x)，其迭代步骤如下：

  1. 初始化：选择一个初始点 x^0。
  2. 次梯度步：在当前点 x^k，计算凸函数 h 在 x^k 处的一个次梯度 ξ^k ∈ ∂h(x^k)。
  3. 凸优化步：求解下一个迭代点 x^{k+1}，它是如下凸优化问题的解：x^{k+1} ∈ 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>]。最小化这个凸上界函数，就得到了下一个迭代点。
  4. 迭代：重复步骤2和3，直到收敛。

DCA的魅力在于，每一步都只需要计算一个次梯度和求解一个凸优化问题。只要凸子问题能高效求解（例如，g和h都是凸多项式时，子问题可能是一个凸的二次或半定规划），整个算法就非常可行。并且，DCA具有收敛到临界点（通常是局部极小值）的理论保证。

现在，将我们前面得到的多项式DC分解与DCA结合，就形成了一套处理非凸多项式优化的完整流程：

  1. 建模与分解：对目标多项式f，运行前述的DC分解算法，得到凸多项式g和h。
  2. 实施DCA：应用DCA算法进行优化。
  3. 局部精炼：DCA收敛到一个临界点后，可以将其作为初始点，送入局部优化器（如基于梯度的方法）进行精炼。

实战中的挑战与心得：

  1. 分解的质量决定算法的命运：DCA的收敛速度和最终解的质量极度依赖于DC分解的质量。一个“好”的分解应该使得g和h“尽可能凸”，但又不是通过引入巨大的二次项（如平凡的利普希茨分解）来实现。理想的分解是g和h的曲率都相对平缓，这样DCA每一步的凸子问题条件数更好，更容易求解，且迭代路径更平滑。我们之前的分解算法旨在寻找“最近”的凸函数，这在一定程度上能产生质量较好的分解。

  2. 次梯度的选择与凸子问题的求解：对于光滑的凸多项式h，其次梯度就是梯度∇h(x)。凸子问题min_x g(x) - <∇h(x^k), x>是一个无约束凸多项式最小化问题。当g是二次凸函数时，子问题有解析解。当g是更高阶的凸多项式时，子问题本身可能仍需用迭代法求解，但由于它是凸的，且通常形式比原问题简单，求解依然比原非凸问题容易得多。我的经验是，对于中低阶多项式（如4阶以下），子问题可以用内点法快速求解；对于更高阶的，可能需要专门的多项式优化求解器。

  3. 多个局部解与初始化：DCA只能保证收敛到临界点，而多项式优化问题通常有多个局部极小值。因此，算法对初始点x^0敏感。在实际应用中，需要结合多种策略：a) 从多个随机初始点运行DCA；b) 利用多项式优化本身的全局方法（如Lasserre的SOS层次结构）提供一个下界，并启发式地产生初始点；c) 使用我们得到的DC分解结构，设计分支定界框架，以寻求全局最优解。

  4. 复杂度瓶颈：整个流程的复杂度瓶颈主要在第一步——DC分解。求解大规模的SDP问题（尤其是高维、高阶多项式）是计算昂贵的。这是将方法应用于大规模问题的首要障碍。研究更高效的、可扩展的DC分解算法（例如，基于一阶方法近似求解SDP，或利用深度学习学习分解）是当前的前沿方向。

5. 超越多项式：DC分解思想的延伸与工具生态

虽然本文聚焦于多项式，但DC分解的思想和DCA算法具有极大的通用性。任何可以表示为两个凸函数之差的函数，都可以纳入这个框架。这包括许多带激活函数的神经网络、一些特殊的超越函数组合等。

对于非多项式函数，获得其DC分解可能没有系统性的算法，但通常可以通过函数本身的特性进行“手工”分解或利用已知的凸凹性质。例如：

  • log(1+exp(x))(log-sum-exp) 本身是凸函数。
  • x^4 - 2*x^2可以分解为(x^4 + a*x^2) - ((a+2)*x^2)，通过选择合适的a>0，可以使两项均为凸函数。
  • 神经网络中常用的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层次结构，计算成本极高）之间。对于具有特殊结构（如稀疏性、低有效维数）的中等规模非凸多项式优化问题，这套基于凸表示的方法往往能在可接受的时间内，找到高质量的解，甚至通过全局化策略证明全局最优性。理解其理论边界，熟练运用现有的工具链，并根据具体问题设计分解策略，是发挥其最大效力的关键。这需要不断的实践和试错，但一旦掌握，它就成为了你处理复杂非线性问题工具箱中一件非常犀利的武器。