用于最小二乘拟合的 Cython/numpy 与纯 numpy答案

【问题标题】：Cython/numpy vs pure numpy for least square fitting用于最小二乘拟合的 Cython/numpy 与纯 numpy
【发布时间】：2016-07-15 03:24:27
【问题描述】：

学校的一位 T.A 向我展示了这段代码作为最小二乘拟合算法的示例。

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))

【问题讨论】：

最好的办法是保存一些您多次使用的中间结果（例如M.T*V.I）。此外，当矩阵求逆立即与另一个矩阵相乘时，通常有一些方法可以避免矩阵求逆，这对速度和数值精度都有好处。

标签： python numpy cython least-squares

【解决方案1】：

除了[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)

做同样的事情。

所以不，这看起来不像是一个好的 cython 候选者——除非你准备好以可编译的细节写出所有这些矩阵产品。

我也会跳过 asmatrix 的内容并使用常规的 dot、@ 和 einsum，但这更多的是风格问题而不是速度问题。

【讨论】：