python - cython 与 python 在 scipy.optimize.fsolve 中的不同结果

Question

我用 cython 化了一个在我的代码中多次调用的函数。 cython 版本和原始 python 代码给了我相同的答案(在 1e-7 内，我理解这与 cython 与 python 类型有关......不是这里的问题，但可能很重要) .

我尝试使用 scipy.optimize.fsolve() 找到函数的根。 python 版本工作正常，但 cython 版本不同。

代码非常复杂，并且有一个很大的外部文件来准备一些参数，所以我无法发布所有内容。我发布了cython代码。完整代码是 here .

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 编译扩展模块。

Answer 1

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

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

---

顺便说一句，目标函数中的舍入误差阻止收敛的情况表明 ill-conditioned problems .您可能想考虑是否可以重新制定模型以提高其数值稳定性。