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代码和Cython的差不多，只是没有cdef语句和参数的类型信息。我检查过，输入值相同的情况下，输出结果也是一样的。我的问题是，为什么scipy的fsolve在一个版本上会出问题，而在另一个版本上却没事。我猜这可能是我Cython代码的问题？

我是在Anaconda上运行Python 2.7，通过pyximport编译扩展模块。

scipy 数据类型误差分析 cython Anaconda 根查找 optimize fsolve

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

1 个回答

撰写回答