Могу ли я ускорить этот базовый код линейной алгебры?

Мне было интересно, можно ли оптимизировать следующее с помощью Numpy или математического обмана.

def f1(g, b, dt, t1, t2):
  p = np.copy(g)
  for i in range(dt):
    p += t1*np.tanh(np.dot(p, b)) + t2*p
  return p

где g - вектор длины n, b - это матрица n x n, dt - число итераций, а t1 и t2 - скаляры.

Я быстро исчерпал идеи о том, как оптимизировать это дальше, потому что p используется в цикле во всех трех терминах уравнения: при добавлении к себе; в точечном продукте; и в скалярном умножении.

Но, возможно, есть другой способ представить эту функцию или есть другие трюки, чтобы повысить ее эффективность. Если возможно, я бы предпочел не использовать Cython и т.д., Но я был бы готов использовать его, если бы улучшения скорости были значительными. Заранее спасибо, и извиняюсь, если вопрос вне сферы видимости.

Обновление:

Представленные ответы до сих пор более сфокусированы на том, какие значения ввода/вывода могут быть во избежание ненужных операций. Теперь я обновил MWE с правильными значениями инициализации для переменных (я не ожидал, что идеи оптимизации придут с этой стороны - извинения). g будет находиться в диапазоне [-1, 1], а b будет находиться в диапазоне [-infinity, infinity]. Аппроксимация вывода не является опцией, потому что возвращенные векторы позже передаются функции оценки - аппроксимация может возвращать один и тот же вектор для довольно аналогичного ввода, поэтому это не вариант.

MWE:

import numpy as np
import timeit

iterations = 10000

setup = """
import numpy as np
n  = 100
g  = np.random.uniform(-1, 1, (n,)) # Updated.
b  = np.random.uniform(-1, 1, (n,n)) # Updated.
dt = 10
t1 = 1
t2 = 1/2

def f1(g, b, dt, t1, t2):
  p = np.copy(g)
  for i in range(dt):
    p += t1*np.tanh(np.dot(p, b)) + t2*p
  return p
"""

functions = [
  """
    p = f1(g, b, dt, t1, t2)
  """
]

if __name__ == '__main__':
  for function in functions:
    print(function)
    print('Time = {}'.format(timeit.timeit(function, setup=setup,
                                           number=iterations)))

Ответы

Ответ 1

Чтобы ускорить работу кода без cython или jit, некоторые математические обманщики могут быть более легким подходом. Мне кажется, что если мы определяем a k(g, b) = f1(g, b, n+1, t1, t2)/f1(g, b, n, t1, t2) для n в положительном N, то функция k должна иметь предел t1+t2 (пока еще нет твердого доказательства, просто ощущение кишки; - частный случай для E (g) = 0 и E (p) = 0 также.). Для t1=1 и t2=0.5, k(), как представляется, приближается к пределу довольно быстро, для N>100 оно почти константа 1.5.

Итак, я думаю, что подход с численным приближением должен быть самым простым. enter image description here

In [81]:

t2=0.5
data=[f1(g, b, i+2, t1, t2)/f1(g, b, i+1, t1, t2) for i in range(1000)]
In [82]:

plt.figure(figsize=(10,5))
plt.plot(data[0], '.-', label='1')
plt.plot(data[4], '.-', label='5')
plt.plot(data[9], '.-', label='10')
plt.plot(data[49], '.-', label='50')
plt.plot(data[99], '.-', label='100')
plt.plot(data[999], '.-', label='1000')
plt.xlim(xmax=120)
plt.legend()
plt.savefig('limit.png')

In [83]:

data[999]
Out[83]:
array([ 1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,
        1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,
        1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,
        1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,
        1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,
        1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,
        1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,
        1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,
        1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,  1.5,
        1.5])

Ответ 2

Я смущаюсь дать это как ответ, поскольку я думаю, что это может быть артефакт входных данных, которые вы нам дали. Тем не менее, обратите внимание, что tanh(x) ~ 1 для x>>1. Ваши входные данные во все времена, когда я его запускал, имеют x = np.dot(p,b) >> 1, поэтому мы можем заменить f1 на f2.

def f1(g, b, dt, t1, t2):
  p = np.copy(g)
  for i in range(dt):
      p += t1*np.tanh(np.dot(p, b)) + t2*p
  return p

def f2(g, b, dt, t1, t2):
  p = np.copy(g)
  for i in range(dt):
      p += t1 + t2*p
  return p

print np.allclose(f1(g,b,dt,t1,t2), f2(g,b,dt,t1,t2))

Что действительно показывает, что две функции численно эквивалентны. Заметим, что f2 является неоднородным линейным рекуррентным отношением, и его можно решить за один шаг, если вы решите это сделать.