sympy的nsolve求解方程组不断抛出ValueError/数值奇异错误，求解决方法

你当前尝试运行如下代码，基于给定的股票价格序列，通过sympy的nsolve求解含绝对值、指数项的非线性方程组，拟合alpha、beta两个模型参数：

# %% Imports
import numpy as np
import sympy as sy
import scipy as sc

alpha = sy.Symbol('alpha', real=True) 
beta = sy.Symbol('beta', real=True) 
s_0 = sy.Symbol('s_0', real=True) 
s_1 = sy.Symbol('s_1', real=True) 
s_2 = sy.Symbol('s_2', real=True) 
t = sy.Symbol('t', real=True)


#run checks on S(t) whether S(t) > S0 for current iteration
def main():
    # A. Defining equations to constant variables to be called later
    #S(t) > S0, (S(t) + β/α)*S(t) > 0
    MODEL_1_1 = ((beta/alpha) * (1/sy.Abs((s_0 + beta/alpha)/s_0)) * sy.exp(beta*t)) / (1 - (1/sy.Abs((s_0 + beta/alpha)/s_0) * sy.exp(beta*t)))

    #S(t) > S0, (S(t) + β/α)*S(t) < 0
    MODEL_1_2 = -((beta/alpha) * (1/sy.Abs((s_0 + beta/alpha)/s_0)) * sy.exp(beta*t)) / (1 + (1/sy.Abs((s_0 + beta/alpha)/s_0) * sy.exp(beta*t)))

    #S(t) < S0, (S(t) + β/α)*S(t) > 0
    MODEL_1_3 = ((beta/alpha) * (1/sy.Abs((s_0 + beta/alpha)/s_0)) * sy.exp(-sy.Abs(beta)*t)) / (1 - (1/sy.Abs((s_0 + beta/alpha)/s_0) * sy.exp(-sy.Abs(beta)*t)))

    #S(t) < S0, (S(t) + β/α)*S(t) > 0
    MODEL_1_4 = -((beta/alpha) * (1/sy.Abs((s_0 + beta/alpha)/s_0)) * sy.exp(-sy.Abs(beta)*t)) / (1 + (1/sy.Abs((s_0 + beta/alpha)/s_0) * sy.exp(-sy.Abs(beta)*t)))

    #---------------------------------------------------------------
    # B. Initializing values
    stock_price = [25725, 25600, 24600]

    #---------------------------------------------------------------
    # C. Get 2 functions for every t = 1 and t = 2
    ## C.1. processing t1
    if (stock_price[1] > stock_price[0]):
        eq1 = MODEL_1_1.subs([(t,1),(s_0,stock_price[0])]) #(S(t) + β/α)*S(t) > 0
        eq2 = MODEL_1_2.subs([(t,1),(s_0,stock_price[0])]) #(S(t) + β/α)*S(t) < 0
    
    else:
        eq1 = MODEL_1_3.subs([(t,1),(s_0,stock_price[0])]) #(S(t) + β/α)*S(t) > 0
        eq2 = MODEL_1_4.subs([(t,1),(s_0,stock_price[0])]) #(S(t) + β/α)*S(t) < 0
    
    #print(eq1)
    #print(eq2)

    ## C.2. processing t2
    if (stock_price[2] > stock_price[0]):
        eq3 = MODEL_1_1.subs([(t,2),(s_0,stock_price[0])]) #(S(t) + β/α)*S(t) > 0
        eq4 = MODEL_1_2.subs([(t,2),(s_0,stock_price[0])]) #(S(t) + β/α)*S(t) < 0
    
    else:
        eq3 = MODEL_1_3.subs([(t,2),(s_0,stock_price[0])]) #(S(t) + β/α)*S(t) > 0
        eq4 = MODEL_1_4.subs([(t,2),(s_0,stock_price[0])]) #(S(t) + β/α)*S(t) < 0

    #print(eq3)
    #print(eq4)

    #------------------------------------------------
    # D. Run 4 different solution, get 4 sets of possible (α,β) pair
    ## D.1. eq1 vs. eq3  -- (S(1) + β/α)*S(1) > 0, (S(2) + β/α)*S(2) > 0
    print(eq1)
    print(eq3)
    sol1 = sy.solvers.nsolve([sy.Eq(eq1,stock_price[1]), sy.Eq(eq3,stock_price[2])],[alpha,beta],[-1,1])


    print(sol1)
    ## D.2. eq2 vs. eq3  -- (S(1) + β/α)*S(1) < 0, (S(2) + β/α)*S(2) > 0
    ## D.3. eq1 vs. eq4  -- (S(1) + β/α)*S(1) > 0, (S(2) + β/α)*S(2) < 0
    ## D.4. eq2 vs. eq4  -- (S(1) + β/α)*S(1) < 0, (S(2) + β/α)*S(2) < 0 

#check for contradictions against initial condition
#return values that satisfies all condition
#run t = 3 based on the iteration
#repeat
main()

报错信息 #

运行时持续抛出两类错误：
第一类收敛失败报错：

Could not find root within given tolerance. (655317900.776688019361 > 2.16840434497100886801e-19) Try another starting point or tweak arguments.

调整初始迭代点后抛出第二类矩阵奇异报错：

matrix is numerically singular

完整报错回溯如下：

Traceback (most recent call last):
  File "D:\SKRIPSI ANDREE\Code\model1.py", line 76, in <module>
    main()
  File "D:\SKRIPSI ANDREE\Code\model1.py", line 64, in main
    sol1 = sy.solvers.nsolve([sy.Eq(eq1,stock_price[1]), sy.Eq(eq3,stock_price[2])],[alpha,beta],[-1,1])
  File "C:\Users\Andre\AppData\Local\Programs\Python\Python39\lib\site-packages\sympy\utilities\decorator.py", line 88, in func_wrapper
    return func(*args, **kwargs)
  File "C:\Users\Andre\AppData\Local\Programs\Python\Python39\lib\site-packages\sympy\solvers\solvers.py", line 2954, in nsolve
    x = findroot(f, x0, J=J, **kwargs)
  File "C:\Users\Andre\AppData\Local\Programs\Python\Python39\lib\site-packages\mpmath\calculus\optimization.py", line 969, in findroot
    for x, error in iterations:
  File "C:\Users\Andre\AppData\Local\Programs\Python\Python39\lib\site-packages\mpmath\calculus\optimization.py", line 660, in __iter__
    s = self.ctx.lu_solve(Jx, fxn)
  File "C:\Users\Andre\AppData\Local\Programs\Python\Python39\lib\site-packages\mpmath\matrices\linalg.py", line 226, in lu_solve
    A, p = ctx.LU_decomp(A)
  File "C:\Users\Andre\AppData\Local\Programs\Python\Python39\lib\site-packages\mpmath\matrices\linalg.py", line 136, in LU_decomp
    raise ZeroDivisionError('matrix is numerically singular')
ZeroDivisionError: matrix is numerically singular

错误原因 #

• 收敛失败问题：你的模型包含绝对值、指数项、分段规则，属于非光滑函数，nsolve默认采用的牛顿类迭代法高度依赖初始值，且对函数可导性要求高。你给定的初始值[-1,1]和股票价格的万级量级差距过大，很容易导致迭代步长跑飞，无法达到默认收敛精度；同时仅用3个样本点的情况下，方程组大概率不存在精确解析根，自然无法满足nsolve的残差要求。
• 矩阵奇异问题：牛顿迭代每一步需要对雅可比矩阵求逆，当迭代过程中雅可比矩阵行列式趋近于0时，就会出现数值奇异无法求逆的问题。你的模型包含大量除法、指数运算，参数轻微偏移就可能导致雅可比矩阵的行/列线性相关，触发该报错。

替代求解方案 #

• 改用数值最小二乘拟合：放弃符号求解逻辑，直接将方程组残差平方和作为损失函数，调用scipy.optimize.least_squares做优化，不需要方程组存在精确根，对非光滑函数的兼容性更强，还可以给alpha、beta设置合理的取值边界，避免出现除零等数值问题。
• 先做变量降维：引入中间变量k = beta/alpha，先消解模型里的除法和冗余参数，降低方程组的非线性程度，再通过网格搜索先锁定参数的大致取值区间，再带入nsolve迭代，可大幅提升收敛概率。
• 更换数值求解器：调用scipy.optimize.root，选择method='lm'（Levenberg-Marquardt算法）或broyden1等拟牛顿法，这类求解器对雅可比奇异的容忍度远高于nsolve默认的牛顿法，适配性更强。

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