numpy: трансляция в несколько внутренних продуктов и инверсий

У меня есть массивы e (форма q на l), f (форма n на l) и w (форма n на l), и я хочу создать массив M, где M[s,i,j] = np.sum(w[s, :] * e[i, :] * e[j, :]), и массив F, где F[s,j] = np.sum(w[s, :] * f[s, :] * e[j, :]).

И то, и другое достаточно легко сделать, например, перебирая элементы M, но я хочу быть более эффективным (мои реальные данные содержат что-то вроде 1M записей длиной 5k). Для F я могу использовать F = np.inner(w * f, e) (который, как я убедился, дает тот же ответ, что и цикл). M сложнее, поэтому первый шаг — прокрутить нулевое измерение с пониманием списка, говоря, что M = np.stack([np.inner(r[:] * e, e) for r in w]) (я убедился, что это также работает так же, как и цикл). np.inner() не принимает никаких аргументов осей, поэтому мне не ясно, как заставить массивы просто транслировать все строки w.

Наконец, мне нужно использовать элементы M и F для создания матрицы A, где A[s,i] = np.sum(np.linalg.inv(M[s, :, :])[i, :] * F[i, :]). Это также выглядит как внутренний продукт, но выполнение множества отдельных инверсий отнимает много времени, так что есть ли способ вычислить инверсии срезов, кроме цикла?

Некоторые тестовые значения в моих массивах следующие:

e = array([[-0.9840087 , -0.17812043],
           [ 0.17812043, -0.9840087 ]])

w = array([[  1.12545297e+01,   1.48690140e+02],
           [  7.30718244e+00,   4.07840612e+02],
           [  2.62753065e+02,   2.27085711e+02],
           [  1.53045364e+01,   5.63025281e+02],
           [  8.00555079e+00,   2.16207407e+02],
           [  1.54070190e+01,   1.87213209e+06],
           [  2.71802081e+01,   1.06392902e+02],
           [  3.46300255e+01,   1.29404438e+03],
           [  7.77638140e+00,   4.18759293e+03],
           [  1.12874849e+01,   5.75023379e+02]])

f = array([[ 0.48907404,  0.06111084],
           [-0.21899297, -0.02207311],
           [ 0.58688524,  0.05156326],
           [ 0.57407751,  0.10004592],
           [ 0.94172351,  0.03895357],
           [-0.7489003 , -0.08911183],
           [-0.7043736 , -0.19014227],
           [ 0.58950925,  0.16587887],
           [-0.35557142, -0.14530267],
           [ 0.24548714,  0.03221844]])

python numpy broadcasting matrix-inverse inner-product
person DathosPachy    schedule 01.07.2016    source источник
comment
В M[s,i,j] = np.sum(w[s] * e[i] * e[j]) какая сумма закончилась? Я вижу i,j,s с обеих сторон. Это (w[s,k]*e[i,k]*e[j,k]).sum(axis=k)?   -  person hpaulj    schedule 01.07.2016
comment
np.einsum — это быстрый и простой инструмент для вычисления такой суммы произведений. Уравнения s,i,j почти пишутся сами собой.   -  person hpaulj    schedule 01.07.2016
comment
строки w, f и e являются векторами длины l, поэтому поэлементное произведение трех суммируется по его единственной оси   -  person DathosPachy    schedule 01.07.2016
comment
Иногда понятнее написать e[j,:], чем просто e[j], даже если интерпретатор трактует оба одинаково.   -  person hpaulj    schedule 01.07.2016

Ответы (1)


arrow_upward
3
arrow_downward 

M[s,i,j] = np.sum(w[s, :] * e[i, :] * e[j, :])

переводится как

M = np.einsum('sk,ik,jk->sij',w,e,e)

а также

F[s,j] = np.sum(w[s, :] * f[s, :] * e[j, :])
F = np.einsum('sk,sk,jk->sj', w, f, e)

Я не тестировал их с вашими образцами, но перевод достаточно прост.

С действительно большими массивами вам, возможно, придется разбить выражения на части. С 4 переменными итерации общее пространство итераций может быть очень большим. Но сначала посмотрите, работают ли эти выражения с массивами скромного размера.

Что касается

A[s,i] = np.sum(np.linalg.inv(M[s, :, :])[i, :] * F[i, :])

Похоже, np.linalg.inv(M) работает, выполняя s i x i инверсию

Если так, то

 IM = np.linalg.inv(M)
 A = np.einsum('skm,ik,im->si', IM, F)

Я предполагаю больше здесь.

Опять же, размер может стать слишком большим, но сначала попробуйте его уменьшить.

Обычно решения линейных уравнений рекомендуются вместо прямых обратных, например

  A = F/M
  A = np.linalg.solve(M, F)

так как вы, вероятно, хотите A таких, что M@A=F (@ матричный продукт). Но я немного заржавел в этих делах. Также проверьте tensorsolve и tensorinv.

person hpaulj    schedule 01.07.2016
comment
A = np.einsum('skm,ik,im->si', IM, F) имеет слишком много индексов, но A = np.linalg.solve(M, F) создает массив с правильными размерами, что является половиной дела. - person DathosPachy; 01.07.2016
comment
sij,ij-›si` лучше? - person hpaulj; 01.07.2016
comment
Нет, эта корректировка индекса вызывает ошибку с повторной широковещательной передачей операнда [original->remapped]: (50,2,2)->(50,2,2) (50,2)->(50,2). Однако с реализацией np.linalg.solve() этот метод дает тот же ответ, что и f.dot(e.T), когда w везде один. Это ожидаемое поведение, поэтому я уверен, что оно работает. - person DathosPachy; 01.07.2016