Numpy矩阵填充非常慢

Question

我正在写一个程序来执行数值计算与黑森矩阵。黑森矩阵是500 × 500我需要填充它上百次。我每次都用两个for循环填充它。我的问题是速度太慢了。下面是我的代码:

#create these outside function
hess = np.empty([500,500])
b = np.empty([500])
def hess_h(x):
    #create these first so they aren't calculated every iteration
    for k in range(500):
        b[k] = (1-np.dot(a[k],x))**2
    for i in range(500):
        for j in range(500):
            if i == j:
                #these are values along diagonal
                hess[i,j] = float(2*(1-x[i])**2 + 4*x[i]**2)/(1-x[i]**2)**2 
                            - float(a[i,j]*sum(a[i]))/b[i]
            #the matrix is symmetric so only calculate upper triangle
            elif j > i :
                hess[i,j] = -float(a[i,j]*sum(a[i]))/b[i]
            elif i > j:
                hess[i,j] = hess[j,i]
    return hess

我计算hess_h(np.zeros(500))需要10.2289998531秒来运行。这太长了，我需要另想办法

Answer 1

在你的计算中寻找模式，特别是那些你可以在整个i和j范围内计算的东西。

例如，我看到i==j 的对角线

hess[i,j] = float(2*(1-x[i])**2 + 4*x[i]**2)/(1-x[i]**2)**2 
            - float(a[i,j]*sum(a[i]))/b[i]

你能把它改成一次性的表达式吗，比如:

2*(1-x)**2 + 4*x**2)/(1-x**2)**2 - np.diagonal(a)*sum(a)/b

其他部分使用上下三角形元素。有像np.triu这样的函数可以给你它们的索引。

我试着给你一些工具和过程来解决这个问题，用一些numpy的矢量化操作，而不是迭代i和j的所有元素。

看起来像

 -a[i,j]*sum(a[i])/b[i]

用于每个元素。我假设a是一个(500,500)数组。你能用

吗?

 -a*a.sum(axis=?)/b

b可以'矢量化'

 b[k] = (1-np.dot(a[k],x))**2

，例如:

 (1 - np.dot(a, x))**2

或

 (1 - np.einsum('kj,ji',a,x))**2

在较小的a上测试细节