умножение матриц с ndarray

Question

Как выполнить умножение вектора и матрицы где v - вектор (v1, v2, v3), а A - матрица 3x3? Python жалуется на то, что фигуры не выровнены, возможно, потому что v - ndarray. Любые идеи о том, как сделать эту операцию? конечный результат должен быть скалярным в каждой точке координат (v1, v2, v3). Приведенный ниже основной код прерывается при попытке выполнить умножение.

import numpy as np
a = np.linspace(0, 10, 21)
b = np.linspace(0, 20, 41)
a, b = np.meshgrid(a,b)
v = np.array([a*b, a+b, a])
A = np.ones([3,3])
s = v.T @ A @ v     # doesn't work

Error

----> 1 s = v.T @ A @ v    
ValueError: shapes (21,41,3) and (3,41,21) not aligned: 3 (dim 2) != 41 (dim 1)

Редактировать: матричная операция должна выполняться в каждой точке v , где v обычно представляет собой большой массив (векторов). Например, возьмите куб 1 м с центром в начале координат и оцените матричную операцию в каждой точке сетки, скажем, каждые 10 см в каждой оси координат.

редактировать 2 пример для отдельной точки в (x, y, z)

A = np.zeros([3,3])
A[0][0] = 1
A[1][1] = 2
A[2][2] = 3
x,y,z = 1, 1, 0
v = np.array([x, y, z])
s = v.T @ A @ v   # should give s=3

Следующий шаг - заставить код работать для большого массива векторов v . За исключением того, что это немного сложнее, потому что координаты вектора (x, y, z) должны быть параметризованы в терминах координат (a, b). Исходный код выше пытается это сделать, но не работает и может быть не лучшим подходом. Есть другие идеи?

Answer 1

Когда вы умножаете две N-мерные матрицы на numpy, я предполагаю, что он автоматически умножает два последних измерения вместе и сохраняет первые. Умножение N-мерной матрицы более или менее похоже на умножение 2D. Ваши матрицы должны иметь одинаковую форму, за исключением двух измерений. В этих двух измерениях вы должны соблюдать правила двумерного умножения. Например, если у вас есть матрица A с формой (a,b,c, ...,d,e,f) и вы хотите умножить ее на матрицу B, форма B должна быть (a,b,c,...,d,f,g), а форма результата будет (a,b,c, ...,d,e,g).

Давайте забудем, что мы находимся в четырехмерном пространстве. Если у вас была только одна точка, v^T*A*v должен иметь форму (1,3)x(3,3)x(3,1). Мы просто хотим применить это для каждой точки в сетке (41,21). Это дает нам последние измерения каждого компонента, который нам нужно умножить. Чтобы быть последовательным, v^T*A*v должен иметь форму (41,21,1,3)x(3,3)x(41,21,3,1).

 import numpy as np
 a = np.linspace(0, 10, 21)
 b = np.linspace(0, 20, 41)
 a, b = np.meshgrid(a,b)
 a = np.expand_dims(a, axis=0)
 b = np.expand_dims(b, axis=0)
 print("Shape a = {}, Shape b = {}".format(a.shape, b.shape))
 v = np.array([a*b, a+b, a])
 print("Shape v = {}".format(v.shape))
 u1 = v.transpose((2,3,1,0))
 print("Shape u1 = {}".format(u1.shape))
 s = u1 @ A
 u2 = v.transpose((2,3,0,1))
 print("Shape u2 = {}".format(u2.shape))
 s = s @ u2
 print("{} x {} x {} = {} x {} = {}".format(u1.shape, A.shape, u2.shape, (u1 @ A).shape, u2.shape, s.shape))

возвращает:

Shape a = (1, 41, 21), Shape b = (1, 41, 21)
Shape v = (3, 1, 41, 21)
Shape u1 = (41, 21, 1, 3)
Shape u2 = (41, 21, 3, 1)
(41, 21, 1, 3) x (3, 3) x (41, 21, 3, 1) = (41, 21, 1, 3) x (41, 21, 3, 1) = (41, 21, 1, 1)

Я предлагаю вам это решение. Вы начинаете с добавления размера 1 к вашим векторам a и b. Вместо того, чтобы иметь форму (41,21), они будут иметь форму (1,41,21). Теперь, когда вы строите v, вы получаете форму (3,1,41,21). Теперь, если вы используете обычную транспозицию, вы просто переворачиваете все измерения, а это не то, что вы хотите. Вы хотите, чтобы v ^ T было умножено на A, формы (3,3). Таким образом, вы вручную определяете, как изменить размеры вектора, чтобы перейти от (3,1,41,21) к (41,21,1,3) и (41,21,3,1). В конце концов, вы можете, наконец, умножить его, и оно будет последовательным.

ПРИМЕЧАНИЕ 1 В теории вы можете умножать на другие измерения, кроме последних, при условии соблюдения правила двумерного умножения для этих измерений на те же другие измерения. Но это способ сделать это в Python.

ПРИМЕЧАНИЕ 2 Вы можете задаться вопросом, почему мы можем умножить матрицу формы (41,21,1,3) на матрицу формы (3,3). Это точно такая же механика, что и при умножении 2D матрицы на скаляр. Когда вы делаете это, вы увеличиваете размер скаляра до двух измерений (в основном это матрица со скаляром везде), и вы выполняете поэлементное умножение. Точно так же вы создаете матрицу формы (41,21,3,3) и умножаете поэлементно, или «блочно» (2D матричное умножение). Элементы дают умножение (1,3)x(3,3).

Answer 2

Кажется, упоминая вектор v из трех элементов, вы имели в виду ndarray с тремя элементами вдоль его первой оси, каждый из которых содержит данные массива n-dim. Для приведенного примера у вас есть то же самое, что и 3D-массив. Также представляется, что выходной сигнал должен быть уменьшен до скаляра для каждого вектора из трех элементов, то есть выходной сигнал будет двухмерным. Итак, чтобы решить ваш случай, нам нужно использовать тензорное умножение для первого: V.T @ A суммирование по первым осям, что дает нам трехмерный массив. Затем используйте einsum, чтобы выровнять первые две оси, и суммируйте последние последние, как показано ниже -

p1 = np.tensordot(v,A,axes=((0,0)))
out = np.einsum('jkl,ljk->jk',p1,v)

В качестве альтернативы, используя einsum, мы можем сделать все за один шаг, например, -

out = np.einsum('ijk,il,ljk->jk',v,A,v)

Возможно, мы сможем сделать это быстрее с помощью einsum's необязательный аргумент: optimize, установленный как True:

np.einsum(..., optimize=True)