Как ускорить вычисление векторного перекрестного произведения
Привет, я относительно новичок здесь и пытаюсь сделать некоторые расчеты с NumPy. Я испытываю длительное время, затрачиваемое на один конкретный расчет, и не могу найти более быстрый способ добиться того же.
По сути, это часть алгоритма пересечения треугольника с лучами, и мне нужно вычислить все произведения векторных кроссов из двух матриц разных размеров.
Код, который я использовал, был:
allhvals1 = numpy.cross( dirvectors[:,None,:], trivectors2[None,:,:] )
где dirvectors это массив n* vectors (xyz) а также trivectors2 это массив m*vectors(xyz), allhvals1 массив перекрестных произведений размера n*M*vector (xyz), Это работает, но очень медленно. По сути, это матрица n*m каждого вектора из каждого массива. Надеюсь, что вы понимаете. Размеры каждого из них варьируются от 1 до 4000 в зависимости от параметров (я в основном делю на куски в зависимости от размера).
Любой совет приветствуется. К сожалению, моя матричная математика несколько ненормальная.
2 ответа
Если вы посмотрите на исходный код np.cross, это в основном перемещает xyz Размер до начала кортежа формы для всех массивов, а затем вычисление каждого из компонентов изложено следующим образом:
x = a[1]*b[2] - a[2]*b[1]
y = a[2]*b[0] - a[0]*b[2]
z = a[0]*b[1] - a[1]*b[0]
В вашем случае каждый из этих продуктов требует выделения огромных массивов, поэтому общее поведение не очень эффективно.
Давайте настроим некоторые тестовые данные:
u = np.random.rand(1000, 3)
v = np.random.rand(2000, 3)
In [13]: %timeit s1 = np.cross(u[:, None, :], v[None, :, :])
1 loops, best of 3: 591 ms per loop
Мы можем попытаться вычислить его, используя символы Леви-Чивита и np.einsum следующее:
eijk = np.zeros((3, 3, 3))
eijk[0, 1, 2] = eijk[1, 2, 0] = eijk[2, 0, 1] = 1
eijk[0, 2, 1] = eijk[2, 1, 0] = eijk[1, 0, 2] = -1
In [14]: %timeit s2 = np.einsum('ijk,uj,vk->uvi', eijk, u, v)
1 loops, best of 3: 706 ms per loop
In [15]: np.allclose(s1, s2)
Out[15]: True
Так что, пока он работает, он имеет худшую производительность. Дело в том, что np.einsum возникают проблемы, когда имеется более двух операндов, но оптимизированы пути для двух или менее. Таким образом, мы можем попытаться переписать его в два этапа, чтобы посмотреть, поможет ли это:
In [16]: %timeit s3 = np.einsum('iuk,vk->uvi', np.einsum('ijk,uj->iuk', eijk, u), v)
10 loops, best of 3: 63.4 ms per loop
In [17]: np.allclose(s1, s3)
Out[17]: True
Бинго! Близко к порядку улучшения...
Некоторые показатели производительности для NumPy 1.11.0 с a=numpy.random.rand(n,3), b=numpy.random.rand(n,3):
Вложенные einsum примерно в два раза быстрее cross для самого большого n испытания.
При написании динамических симуляций для подводных аппаратов я нашел этот метод для быстрого перекрестного произведения:
Что хорошо работает, написано в Matlab, но код очень прост. Просто прочитайте комментарии вверху.