Cython/numpy vs純numpy最小二乘擬合

T.A在學校向我展示了這個代碼作爲最小二乘擬合算法的一個例子。Cython/numpy vs純numpy最小二乘擬合

import numpy as np 
#return the coefficients (a0,..aN) of the fit y=a0+a1*x+..an*x^n 
#with associated sigma dy 
#x,y,dy are all np.arrays with dtype= np.float64 
def fit_poly(x,y,dy,n): 
    V = np.asmatrix(np.diag(dy**2)) 
    M = [] 

    for k in range(n+1): 
     M.append(x**k) 
    M = np.asmatrix(M).T 
    theta = (M.T*V.I*M).I*M.T*V.I*np.asmatrix(y).T 
    cov_t = (M.T*V.I*M).I 

    return np.asarray(theta.T)[0], np.asarray(cov_t)

我試圖用cython優化他的代碼。我得到這個代碼

cimport numpy as np 
import numpy as np 
cimport cython 

@cython.boundscheck(False) 
@cython.wraparound(False) 
cpdef poly_c(np.ndarray[np.float64_t, ndim=1] x , 
    np.ndarray[np.float64_t, ndim=1] y np.ndarray[np.float64_t,ndim=1]dy , np.int n): 

    cdef np.ndarray[np.float64_t, ndim=2] V, M 

    V=np.asmatrix(np.diag(dy**2),dtype=np.float64) 
    M=np.asmatrix([x**k for k in range(n+1)],dtype=np.float64).T 

    return ((M.T*V.I*M).I*M.T*V.I*(np.asmatrix(y).T))[0],(M.T*V.I*M).I

但運行時似乎是這兩個程序是相同的，我沒有使用「斷言」，其中相同的，以確保輸出。我錯過了什麼/做錯了什麼？

謝謝你的時間，希望你能幫助我。

PS：這是代碼即時分析與（不知道我是否可以調用此配置，但W/E）

import numpy as np 
from polyC import poly_c 
from time import time 
from pancho_fit import fit_poly 
#pancho's the T.A,sup pancho 
x=np.arange(1,1000) 
x=np.asarray(x,dtype=np.float64) 
y=3*x+np.random.random(999) 
y=np.asarray(y,dtype=np.float64) 
dy=np.array([y.std() for i in range(1,1000)],dtype=np.float64) 
t0=time() 
a,b=poly_c(x,y,dy,4) 
#a,b=fit_poly(x,y,dy,4) 
print("time={}s".format(time()-t0))

來源

2016-07-15 cnmoya

最好的辦法是保存一些中間結果，您可以使用多次（例如'M.T * V.I'）。當它們立即乘以另一個對速度和數字精度有利的矩陣時，通常還有避免矩陣求逆的方法。 – DavidW

除了[x**k for k in range(n+1)]我沒有看到任何迭代cython改善。大部分行動都在矩陣產品中。這些已經通過編譯代碼完成（np.dot爲ndarrays）。

和n只有4個，迭代次數不多。

但爲什麼迭代呢？

In [24]: x=np.arange(1,1000.) 
In [25]: M1=x[:,None]**np.arange(5) 
# np.matrix(M1)

做同樣的事情。

所以不，這看起來不像一個好的候選人 - 除非你準備把所有這些矩陣產品編寫成可編譯的細節。

我也跳過asmatrix東西，並使用常規的dot,@和einsum，但這更多的是風格而不是速度。

來源

2016-07-15 03:52:11 hpaulj

Cython/numpy vs純numpy最小二乘擬合

回答

相關問題