使用Scipy optimize.root求解方程组非负根的问题

我尝试寻找不动点f(x)=x，将其转化为求解g(x)=f(x)-x的根，实现代码如下：

from scipy.optimize import root
import numpy as np

def get_attraction(b1,b2, b3):
    c1 = 20 * b1 + 12 * b2 + 6 * b3
    c2 = 12 * b1 + 24 * b2 + 18 * b3
    c3 = 0 * b1 + 14 * b2 + 30 * b3
    return c1,c2,c3

def get_a_tau_max(p1_tau, p2_tau, p3_tau, a):
    b1_tau_l1_a1_a1 = p1_tau * (1 - a) + a
    b2_tau_l1_a1_a1 = p2_tau * (1 - a)
    b3_tau_l1_a1_a1 = p3_tau * (1 - a)

    a1_tau_l1_a1, a2_tau_l1_a1, a3_tau_l1_a1 = get_attraction(b1_tau_l1_a1_a1, b2_tau_l1_a1_a1,
                                                                   b3_tau_l1_a1_a1)


    b1_tau_l1_a1_a2 = p1_tau * (1 - a)
    b2_tau_l1_a1_a2 = p2_tau * (1 - a) + a
    b3_tau_l1_a1_a2 = p3_tau * (1 - a)
    a1_tau_l1_a2, a2_tau_l1_a2, a3_tau_l1_a2 = get_attraction(b1_tau_l1_a1_a2, b2_tau_l1_a1_a2,
                                                                   b3_tau_l1_a1_a2)

    b1_tau_l1_a1_a3 = p1_tau * (1 - a)
    b2_tau_l1_a1_a3 = p2_tau * (1 - a)
    b3_tau_l1_a1_a3 = p3_tau * (1 - a) + a
    a1_tau_l1_a3, a2_tau_l1_a3, a3_tau_l1_a3 = get_attraction(b1_tau_l1_a1_a3, b2_tau_l1_a1_a3,
                                                                   b3_tau_l1_a1_a3)

    a_tau_max = max(a1_tau_l1_a1, a2_tau_l1_a2, a3_tau_l1_a3)

    return a_tau_max

def get_prob(x, s1, s2, s3, a, lamda,p11,p12,p13,p21,p22,p23,p31,p32,p33):
    p1, p2, p3 = x[0], x[1], x[2]

    b1 = s1 * (1-a) + p1 * a
    b2 = s2 * (1-a) + p2 * a
    b3 = s3 * (1-a) + p3 * a

    c1,c2,c3=get_attraction(b1,b2,b3)
    
    t1 = get_a_tau_max(p11, p12, p13, a)
    t2 = get_a_tau_max(p21, p22, p23, a)
    t3 = get_a_tau_max(p31, p32, p33, a)

    a1 = c1+t1
    a2 = c2+t2
    a3 = c3+t3

    nom1 = np.exp(lamda*a1)
    nom2 = np.exp(lamda*a2)
    nom3 = np.exp(lamda*a3)
    dem = nom1 + nom2 + nom3

    p1_t = nom1 / dem
    p2_t = nom2 / dem
    p3_t = nom3 / dem

    return [p1_t - p1, p2_t - p2, p3_t - p3]


s1=2.0454331980638877e-42
s2=1.4352264582396677e-21
s3=1.0
a=11/13
lamda=4.0
p11=0.999999614695216
p12=3.853047896087925e-07
p13=5.242553502503789e-19
p21=1.425164081709055e-21
p22=0.9999999992759427
p23=7.240549367895844e-10
p31=2.03109266273484e-42
p32=1.4251640827409756e-21 
p33=1.0
result = root(get_prob, x0=np.array([0.1, 0.4, 0.5]), args=(s1, s2, s3, a, lamda,p11,p12,p13,p21,p22,p23,p31,p32,p33), method="lm")
print(result)

运行结果：

message: The relative error between two consecutive iterates is at most 0.000000
 success: True
  status: 2
     fun: [ 1.773e-33  8.479e-28  0.000e+00]
       x: [-1.773e-33  7.462e-26  1.000e+00]
   cov_x: [[ 1.000e+00  0.000e+00  0.000e+00]
           [ 0.000e+00  1.000e+00 -3.309e-24]
           [ 0.000e+00 -3.309e-24  1.000e+00]]
    nfev: 13
    fjac: [[-1.000e+00 -3.309e-24 -1.000e+00]
           [-3.309e-24  1.000e+00 -3.309e-24]
           [ 0.000e+00  0.000e+00  1.000e+00]]
    ipvt: [3 2 1]
     qtf: [ 4.441e-16  3.888e-16  8.130e-22]

补充说明：

• s1、s2、s3是通过logistic函数计算的概率
• p11、p12等矩阵元素因精度问题未严格求和为1
• 因代码中使用exp函数，期望得到非负根，但依赖初始值x0有时会得到极小负根

疑问：

1. 这是否是root方法的精度问题？
2. 如何获取非负根？（注：求根已是我最小化问题的一部分，无法改用optimize.minimize）

1. 极小负根的成因 #

确实和数值精度有关：

• 目标函数中p1_t、p2_t是softmax输出，理论上完全非负，但当数值极小时，浮点数计算误差会导致结果出现极小量级的负数（比如1e-33）。
• root方法的终止条件是迭代相对误差足够小，这种极小负数已经满足g(x)≈0的要求，因此会被作为有效根返回。

2. 获取非负根的方法 #

方法一：结果后处理 #

负根的绝对值极小，可直接截断非负后归一化：

• 将结果中小于0的分量设为0，再调整三个分量使其和为1。
• 比如你的结果x: [-1.773e-33, 7.462e-26, 1.000e+00]，处理后得到[0, 0, 1]，代入get_prob计算会发现fun值几乎为0，完全满足不动点条件。

方法二：修改目标函数加入非负惩罚 #

在get_prob中对负分量添加惩罚项，引导优化器避开负数解：

def get_prob(x, s1, s2, s3, a, lamda,p11,p12,p13,p21,p22,p23,p31,p32,p33):
    p1, p2, p3 = x[0], x[1], x[2]
    
    # 添加惩罚项：对负分量施加惩罚
    penalty = 1e6  # 惩罚系数可根据情况调整
    # 先修正输入的负分量
    p1_clamped = max(p1, 0)
    p2_clamped = max(p2, 0)
    p3_clamped = max(p3, 0)
    
    b1 = s1 * (1-a) + p1_clamped * a
    b2 = s2 * (1-a) + p2_clamped * a
    b3 = s3 * (1-a) + p3_clamped * a

    c1,c2,c3=get_attraction(b1,b2,b3)
    
    t1 = get_a_tau_max(p11, p12, p13, a)
    t2 = get_a_tau_max(p21, p22, p23, a)
    t3 = get_a_tau_max(p31, p32, p33, a)

    a1 = c1+t1
    a2 = c2+t2
    a3 = c3+t3

    nom1 = np.exp(lamda*a1)
    nom2 = np.exp(lamda*a2)
    nom3 = np.exp(lamda*a3)
    dem = nom1 + nom2 + nom3

    p1_t = nom1 / dem
    p2_t = nom2 / dem
    p3_t = nom3 / dem

    # 对原始输入的负分量添加惩罚
    res1 = p1_t - p1_clamped
    res2 = p2_t - p2_clamped
    res3 = p3_t - p3_clamped
    if x[0] < 0:
        res1 += penalty * abs(x[0])
    if x[1] < 0:
        res2 += penalty * abs(x[1])
    if x[2] < 0:
        res3 += penalty * abs(x[2])
        
    return [res1, res2, res3]

方法三：优化初始值 #

根据参数特征选择更贴合预期的初始值：

• 你的s3=1，计算出的a3远大于a1和a2，softmax输出后p3_t会非常接近1，直接将初始值设为[0, 0, 1]，root会更快收敛到非负的稳定解，避免出现负根。

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