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.

Answer 1

正如上面的评论中所提到的,Python和Cython版本的结果之间存在差异的原因是在Cython函数中,几个输入被声明为float,而实际的Python变量是double精度的.

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

另外,目标函数中的舍入误差阻止收敛的情况表示病态问题.您可能想要考虑是否可以重新制定模型以提高其数值稳定性.