У меня есть код, который только что был завершен. Работает как задумано. Я решил использовать dot в Numpy, поскольку, по моему ограниченному опыту, это быстрее, чем обычные формы умножения матриц, если в вашей системе установлен BLAS. Тем не менее, вы заметите, чем мне пришлось перенести много вещей. Я уверен, что в этом случае перевешивает пользу от использования dot.
. Я предоставляю математическое уравнение, найденное в статье. Обратите внимание, что это рекурсия
Вот моя реализация кода:
Сначала я предоставляю необходимые компоненты и их размеры
P = array([[0.73105858, 0.26894142],
[0.26894142, 0.73105858]]) # shape (K,K)
B = array([[6.07061629e-09, 0.00000000e+00],
[0.00000000e+00, 2.57640371e-10]]) # shape (K,K)
dP = array([[[ 0.19661193, -0.19661193],
[ 0. , 0. ]],
[[ 0. , 0. ],
[ 0.19661193, -0.19661193]],
[[ 0. , 0. ],
[ 0. , 0. ]],
[[ 0. , 0. ],
[ 0. , 0. ]]]) # shape (L,K,K)
dB = array([[[ 0.00000000e+00, 0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00]],
[[ 0.00000000e+00, 0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00]],
[[-1.16721049e-09, 0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00]],
[[ 0.00000000e+00, 0.00000000e+00],
[ 0.00000000e+00, -1.27824683e-09]]]) # shape (L,K,K)
a = array([[[ 7.60485178e-08, 5.73923956e-07]],
[[-5.54100398e-09, -8.75213012e-08]],
[[-1.25878643e-08, -1.48361081e-07]],
[[-2.73494035e-08, -1.74585971e-07]]]) # shape (L,1,K)
alpha = array([[0.11594542, 0.88405458]]) # shape (1,K)
c = 1 # scalar
Вот фактический Numpy расчет. Обратите внимание на все варианты использования транспонирования.
term1 = (a/c).dot(P).dot(B)
term2 = dP.transpose((0,2,1)).dot(alpha.T).transpose((0,2,1)).dot(B)
term3 = dB.dot( P.T.dot(alpha.T) ).transpose((0,2,1))
a = term1 + term2 + term3
Затем нужно получить:
>>> a
array([[[ 1.38388584e-10, -5.87312190e-12]],
[[ 1.05516813e-09, -4.47819530e-11]],
[[-3.76451117e-10, -2.88160549e-17]],
[[-4.06412069e-16, -8.65984406e-10]]])
Обратите внимание, что форма alpha, а также a была выбрана мной. Их можно изменить, если вы обнаружите, что они обеспечивают превосходную производительность.
Я хотел бы отметить, что я думаю, что существующий код работает быстро. На самом деле, очень быстро. Тем не менее, я продолжаю задаваться вопросом, можно ли сделать лучше. Пожалуйста, дайте ему go. Я профилировал свой код (в котором много Numpy вещания и векторизации), и это не обязательно является узким местом, так как для оценки на моей очень старой машине требуется 23 микросекунды. Тем не менее, это один шаг рекурсии. Это означает, что он оценивается N раз последовательно. Следовательно, даже малейший выигрыш будет иметь большое значение для большой последовательности.
Спасибо за ваше время.
РЕДАКТИРОВАТЬ / ОБНОВИТЬ:
Спасибо @ max9111, который предложил мне посмотреть на этот вопрос здесь , я удалось обработать некоторый код Numba, который работает быстрее, чем расчет Numpy для a. Это занимает 14 микросекунд, в отличие от исходных 23.
Вот оно:
import numba as nb
@nb.njit(fastmath=True,parallel=True,boundscheck=False)
def get_next_a(a,alpha,P,dP,B,dB,c):
N,M,_ = dP.shape
new_a = np.zeros((N,1,M),dtype=np.float64)
new_a = np.zeros((N,1,M))
entry = 0
for idx in nb.prange(N):
for i in range(M):
for j in range(M):
term1 = a[idx,0,j]*P[j,i]*B[i,i]/c
term2 = alpha[0,j]*dP[idx,j,i]*B[i,i]
term3 = alpha[0,j]*P[j,i]*dB[idx,i,i]
entry += term1 + term2 + term3
new_a[idx,0,i] = entry
entry = 0
return new_a
Видно, что get_next_a возвращает желаемый результат. Однако, когда я вызываю его в чистой функции python, которая содержит Numpy, он жалуется. Вот фрагмент моего действительного кода:
def forward_recursions(X,working_params):
# P,dP,B,dB,pi,dpi = setup(X,working_params)
# Dummy Data and Parameters instead of setup
working_params = np.random.uniform(0,2,size=100)
X = np.random.uniform(0,1,size=1000)
P = np.random.uniform(0,1,size=(10,10))
norm = P.sum(axis=1)
P = P/norm[:,None]
dP = np.random.uniform(-1,1,size=(100,10,10))
# We pretend that all 1000 of the 10 x 10 matrices
# only have diagonal entries
B = np.random.uniform(0,1,size=(1000,10,10))
dB = np.random.uniform(0,1,size=(1000,100,10,10))
pi = np.random.uniform(0,1,size=10)
norm = pi.sum()
pi = (pi/norm).reshape(1,10)
dpi = np.random.uniform(0,1,size=(100,1,10))
T = len(X)
N = len(working_params)
M = np.int(np.sqrt(N))
ones = np.ones((M,1))
alpha = pi.dot(B[0])
scale = alpha.dot(ones)
alpha = alpha/scale
ll = np.log(scale)
a = dpi.dot(B[0]) + dB[0].dot(pi.T).transpose((0,2,1))
for t in range(1,T):
old_scale = scale
alpha = alpha.dot(P).dot(B[t])
scale = alpha.dot(ones)
ll += np.log(scale)
alpha = alpha/scale
# HERE IS THE NUMBA FUNCTION
a = get_next_a(a,alpha,P,dP,B[t],dB[t],old_scale)
dll = a.dot(ones).reshape((N,1))/scale
return ll,dll,a
Я знаю, что включение моего собственного кода зависит от других функций, которые не включены, и, следовательно, означает, что forward_recursions не будет работать. Я просто надеюсь, что это, возможно, даст некоторую перспективу.
Я получаю ошибку:
TypingError: Invalid use of Function(<built-in function iadd>) with argument(s) of type(s): (Literal[int](0), array(float64, 2d, C))
Known signatures:
* (int64, int64) -> int64
* (int64, uint64) -> int64
* (uint64, int64) -> int64
* (uint64, uint64) -> uint64
* (float32, float32) -> float32
* (float64, float64) -> float64
* (complex64, complex64) -> complex64
* (complex128, complex128) -> complex128
* parameterized
In definition 0:
All templates rejected with literals.
In definition 1:
All templates rejected without literals.
In definition 2:
All templates rejected with literals.
In definition 3:
All templates rejected without literals.
In definition 4:
All templates rejected with literals.
In definition 5:
All templates rejected without literals.
In definition 6:
All templates rejected with literals.
In definition 7:
All templates rejected without literals.
In definition 8:
All templates rejected with literals.
In definition 9:
All templates rejected without literals.
In definition 10:
All templates rejected with literals.
In definition 11:
All templates rejected without literals.
In definition 12:
All templates rejected with literals.
In definition 13:
All templates rejected without literals.
In definition 14:
All templates rejected with literals.
In definition 15:
All templates rejected without literals.
This error is usually caused by passing an argument of a type that is unsupported by the named function.
[1] During: typing of intrinsic-call at <ipython-input-251-50e636317ef8> (13)
Я не понимаю, что это значит. Вы, наверное, знаете, как я могу исправить что-то подобное. Большое спасибо за ваше время.