Cython与Python在scipy.optimize.fsolve中产生不同的结果

Question

我在代码中调用了一个函数，我将其进行了Cython优化。Cython版本和原始Python代码给出的答案相同（误差小于1e-7，这可能与Cython和Python类型有关，但不是本问题的重点）。

我尝试使用scipy.optimize.fsolve()找到该函数的根。Python版本可以正常工作，但Cython版本发散了。

由于代码比较复杂，并且有一个大的外部文件来准备一些参数，因此我无法发布全部内容。我会发布Cython代码。完整代码在这里。

def euler_outside(float b_prime, int index_b,
                  np.ndarray[np.double_t, ndim=1] b_grid, int index_y,
                  np.ndarray[np.double_t, ndim=1] y_grid,
                  np.ndarray[np.double_t, ndim=1] y_vec,
                  np.ndarray[np.double_t, ndim=2] pol_mat_b, float q,
                  np.ndarray[np.double_t, ndim=2] pol_mat_q,
                  np.ndarray[np.double_t, ndim=2] P, float beta,
                  int n_ygrid, int check=0):
    '''
    b_prime - the variable of interest. want to find b_prime that solves this
    function
    '''
    cdef double b, y, c, uc, e_ucp, eul_val
    cdef int i
    cdef np.ndarray[np.float64_t, ndim=1] uct, c_prime = np.zeros((n_ygrid,))

    b = b_grid[index_b]
    y = y_grid[index_y]

    # Get value of consumption today
    c = b + y - b_prime/q

    # Get possible values of consumption tomorrow
    if check:
        c_prime = b_prime + y_vec - b_grid[0]/q
    else:
        for i in range(n_ygrid):
            c_prime[i] = (b_prime + y_vec[i] -
                         (np.interp(b_prime, b_grid, pol_mat_b[:,i]) /
                          np.interp(b_prime, b_grid, pol_mat_q[:,i])))

    if c<0:
        return 1e10

    uc = utility_prime(c)
    uct = utility_prime(c_prime)

    e_ucp = np.inner( uct, P[index_y,:] )
    eul_val = uc - beta*q * e_ucp

    return eul_val

这个 Python 代码与带有 cdef 语句和参数类型信息的代码相同。我已经检查过，对于相同的输入值，输出结果是相同的。我的问题是为什么 scipy 的 fsolve 对其中一个函数失效，而对另一个函数不失效。我猜测这是因为我的 Cython 出了问题？

我使用 Anaconda 中的 Python 2.7 运行该代码。通过 pyximport 编译扩展模块。

- benten

---

你正在声明一些输入变量为 float。它们真的是浮点数吗，还是实际上是双精度浮点数？将双精度浮点数向下转换为浮点数可能会导致 Python 和 Cython 结果之间的差异。如果你的系统状况不佳，在梯度中额外的舍入误差可能会导致无法收敛。 - ali_m

感谢 @ali_m。将浮点数更改为双精度浮点数，使得该函数的输出与原始Python代码完全相同，并且 fsolve() 收敛于相同的答案。请将其发布为答案，我会接受它。 - benten

Answer 1

如上方评论所述，Python和Cython版本结果不一致的原因在于，在Cython函数中，多个输入被声明为float，而实际的Python变量是双精度double。

这导致了Cython函数舍入误差的增加，似乎是fsolve无法收敛的原因 - 当将这些输入声明为double时，Python和Cython版本产生完全相同的结果，并且fsolve对两者都正确地收敛。

---

作为一个旁注，目标函数中的舍入误差阻止收敛的情况表明问题病态条件。您可能需要考虑是否有可能重新构建模型以提高其数值稳定性。

- ali_m