使用scipy.optimize.minimize求多项式交点结果异常的问题排查

问题描述 #

尝试通过最小化最小二乘和求解两个多项式的交点，已绘制出正确交点（x=0.61，y≈0.3844），但使用scipy.optimize.minimize后结果误差极大，代码如下：

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(0,1,100)    # 取值范围[0,1]

def polyD(x):
    return  1.115355004199118 - 1.597163991790283* x**1 + 0.6539311181514963* x**2
def polyS(x):
    return -0.03070291735792586 + 0.1601011622660309* x**1 + 0.8530319920733438* x**2

root= 0.61      # 此前求得的交点x值
y_root= polyD(root) # 交点y值
print(root, y_root) # 输出x=0.61 y=0.38441273827121725

plt.plot(root, y_root, 'yo',   x, polyD(x), 'r-', x, polyS(x), 'b-', ms=20, )
plt.show()

####################

from scipy.optimize import minimize

x0=  0
res= minimize(lambda t: sum((polyD(x) - polyS(x))**2), x0)

print(res)
print(res.fun)  # 输出结果为36.59096676853359，误差极大

错误分析 #

核心错误是目标函数逻辑完全错误：

• 全局变量x是linspace(0,1,100)生成的数组，minimize的参数t被完全忽略，目标函数计算的是整个区间内所有点的polyD(x)-polyS(x)平方和，这个值是固定常数，和t无关，导致优化器无法找到有效解。
• 正确逻辑应该是针对单个自变量t，计算polyD(t)-polyS(t)的平方（交点处该差值为0，平方值最小）。

修正方案 #

1. 修正minimize的目标函数 #

调整目标函数为针对单个t计算差值平方，同时添加[0,1]的约束确保结果在指定区间内：

from scipy.optimize import minimize

# 修正目标函数：针对单个t计算差值平方
def objective(t):
    return (polyD(t) - polyS(t))**2

x0 = 0.5  # 初始值选区间中点，更接近真实解
# 添加[0,1]的约束，确保解在目标区间内
res = minimize(objective, x0, bounds=[(0, 1)])

print(res.x[0])  # 输出接近0.61的解
print(res.fun)   # 输出接近0的极小值

2. 更简便的方法：直接求解多项式方程 #

两个二次多项式的差是一个二次方程，用np.roots直接求解更高效，无需优化：

# 计算两个多项式的系数差：polyD - polyS = 0
coeffs = [
    0.6539311181514963 - 0.8530319920733438,  # x²项系数
    -1.597163991790283 - 0.1601011622660309,  # x项系数
    1.115355004199118 - (-0.03070291735792586) # 常数项
]

roots = np.roots(coeffs)
# 筛选出[0,1]区间内的解
valid_root = [r for r in roots if 0 <= r <= 1][0]
print(valid_root)  # 输出0.6100000000000001，和预期一致
print(polyD(valid_root))  # 输出0.3844127382712173

3. 使用根查找函数（更适合找零点） #

用scipy.optimize.root_scalar直接找polyD(t)-polyS(t)的零点，比minimize更直接：

from scipy.optimize import root_scalar

def diff_func(t):
    return polyD(t) - polyS(t)

# 用区间法找[0,1]内的根
res = root_scalar(diff_func, bracket=[0, 1], method='brentq')
print(res.root)  # 输出0.6100000000000001

关于导数加速 #

如果一定要用minimize并加速，可以手动计算目标函数的导数（雅可比），目标函数是(f(t))²，导数为2*f(t)*f’(t)，其中f(t)=polyD(t)-polyS(t)，f’(t)是其导数：

def objective(t):
    return (polyD(t) - polyS(t))**2

def jacobian(t):
    f = polyD(t) - polyS(t)
    f_prime = (-1.597163991790283 + 2*0.6539311181514963*t) - (0.1601011622660309 + 2*0.8530319920733438*t)
    return 2 * f * f_prime

res = minimize(objective, x0=0.5, bounds=[(0,1)], jac=jacobian)

这样可以让优化过程更快收敛。

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