numpy matrix population very slow

Question

Я пишу программу для выполнения численного расчета с помощью матрицы Гессиана. Матрица Гессиана 500 х 500, и мне нужно многократно заполнять ее сотни раз. Я каждый раз заполняю его двумя циклами. Моя проблема заключается в том, что она медленно замедляется. Вот мой код:

#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, которые дают вам свои показатели.

Я пытаюсь дать вам инструменты и хотя процессы для решения этого с помощью нескольких векторизованных операций, а не итерации по всем элементам 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.