在Python中拟合全局单调递增多项式至数据集

核心前提 #

要让n阶单变量多项式在全体实数范围内单调递增，其导数（n-1阶多项式）必须全局非负。而只有当n为奇数时，n-1阶多项式是偶次多项式，才有可能满足全局非负（奇次多项式在x→±∞时符号相反，无法全局非负）。若n为偶数，不存在这样的多项式，直接无解。

实现思路 #

利用平方和（SOS）多项式的性质：一个实系数多项式全局非负的充要条件是它可以表示为若干实系数多项式的平方和。通过CVXPY的SOS约束，直接从多项式结构上保证导数全局非负，进而确保原多项式全局单调递增，再结合最小二乘目标拟合数据。

Python代码实现 #

import cvxpy as cp
import numpy as np
import sympy as sp

# 1. 准备数据（替换为你的数据集）
np.random.seed(42)
x_data = np.linspace(-5, 5, 50)
y_data = 2 + 3*x_data + 0.5*x_data**3 + np.random.normal(0, 2, size=x_data.shape)

# 2. 定义多项式阶数（必须为奇数）
n = 3

# 3. 定义多项式系数变量
a = cp.Variable(n + 1)

# 4. 用Sympy构建多项式及其导数
x_sym = sp.symbols('x')
p_sym = sum(a[i] * x_sym**i for i in range(n + 1))
p_prime_sym = sp.diff(p_sym, x_sym)

# 5. 转换为CVXPY表达式并添加SOS约束（保证导数全局非负）
p_prime_cp = cp.from_sympy(p_prime_sym)
constraints = [cp.SOS(p_prime_cp)]

# 6. 定义最小二乘目标函数
x_data_cp = cp.Parameter(x_data.shape, value=x_data)
y_data_cp = cp.Parameter(y_data.shape, value=y_data)
p_x = sum(a[i] * x_data_cp**i for i in range(n + 1))
objective = cp.Minimize(cp.sum_squares(p_x - y_data_cp))

# 7. 求解问题（使用支持SOS约束的求解器，如SCS）
prob = cp.Problem(objective, constraints)
prob.solve(solver=cp.SCS)

# 8. 输出结果并验证
print("拟合得到的多项式系数（a0, a1, ..., an）：")
print(a.value.round(4))

# 验证导数全局非负
test_x = np.linspace(-100, 100, 200)
p_prime_func = sp.lambdify(x_sym, p_prime_sym.subs({a[i]: a.value[i] for i in range(n+1)}), 'numpy')
deriv_values = p_prime_func(test_x)
print(f"导数在[-100, 100]范围内的最小值：{np.min(deriv_values):.4f}")

关键说明 #

• 求解器选择：必须使用支持SOS约束的求解器，如SCS（开源）、MOSEK（商业，精度更高）。
• 阶数限制：务必确保n为奇数，否则约束无法满足，求解会失败。
• 复杂度：多项式阶数n越高，求解速度越慢，可根据拟合需求选择合适的阶数。
• 与局部约束的区别：之前仅在数据区间内取点加导数非负约束，只能保证局部单调；而SOS约束从多项式结构上确保全局非负，彻底解决全局单调的要求。

内容的提问来源于stack exchange，提问作者user999605