Существует ли "улучшенный" метод numpy/scipy dot?

Проблема

Я хотел бы вычислить следующее с помощью numpy или scipy:

Y = A**T * Q * A

где A - матрица m x n, A**T - это транспонирование A, а Q - диагональная матрица m x m.

Так как Q - диагональная матрица, я сохраняю только ее диагональные элементы в качестве вектора.

Способы решения для `Y`

В настоящее время я могу представить два способа расчета Y:

Y = np.dot(np.dot(A.T, np.diag(Q)), A) и
Y = np.dot(A.T * Q, A).

Ясно, что вариант 2 лучше, чем вариант 1, так как никакая реальная матрица не должна создаваться с помощью diag(Q) (если это то, что действительно делает numpy...)
Тем не менее, оба метода страдают от недостатка необходимости выделять больше памяти, чем это действительно необходимо, поскольку A.T * Q и np.dot(A.T, np.diag(Q)) необходимо сохранить вместе с A, чтобы вычислить Y.

Вопрос

Есть ли метод numpy/scipy, который бы устранил ненужное выделение дополнительной памяти, где вы могли бы пропускать только две матрицы A и B (в моем случае B есть A.T) и весовой вектор Q вместе с ним?

Ответ 1

(w/r/t последнее предложение OP: я не знаю такого метода numpy/scipy, но w/r/t Вопрос в заголовке OP (то есть, улучшая производительность NumPy dot), что ниже с другой стороны, мой ответ направлен на повышение производительности большинства шагов, включающих вашу функцию для Y).

Во-первых, это должно дать вам заметный импульс по методу NumPy dot ванили:

>>> from scipy.linalg import blas as FB
>>> vx = FB.dgemm(alpha=1., a=v1, b=v2, trans_b=True)

Обратите внимание, что два массива, v1, v2 находятся в порядке C_FORTRAN

Вы можете получить доступ к байтовому порядку массива NumPy через атрибут массива flags, например:

>>> c = NP.ones((4, 3))
>>> c.flags
      C_CONTIGUOUS : True          # refers to C-contiguous order
      F_CONTIGUOUS : False         # fortran-contiguous
      OWNDATA : True
      MASKNA : False
      OWNMASKNA : False
      WRITEABLE : True
      ALIGNED : True
      UPDATEIFCOPY : False

чтобы изменить порядок одного из массивов, чтобы оба они были выровнены, просто вызовите конструктор массива NumPy, пройдите в массив и установите соответствующий флаг порядка на True

>>> c = NP.array(c, order="F")

>>> c.flags
      C_CONTIGUOUS : False
      F_CONTIGUOUS : True
      OWNDATA : True
      MASKNA : False
      OWNMASKNA : False
      WRITEABLE : True
      ALIGNED : True
      UPDATEIFCOPY : False

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

Но почему массивы копируются перед тем, как их передать точка?

Точечный продукт полагается на операции BLAS. Эти операции требуют массивов, хранящихся в C-смежном порядке - это ограничение, которое вызывает копирование массивов.

С другой стороны, транспонирование не влияет на копию, но, к сожалению, возвращает результат в порядке Fortran:

Следовательно, чтобы устранить узкое место производительности, вам необходимо устранить шаг предикатного массива; для этого просто требуется, чтобы оба массива располагались в C-смежном порядке *.

Итак, чтобы вычислить точку (A.T., A) без дополнительной копии:

>>> import scipy.linalg.blas as FB
>>> vx = FB.dgemm(alpha=1.0, a=A.T, b=A.T, trans_b=True)

В сумме выражение выше (вместе с инструкцией импорта предикатов) может заменить точку, чтобы обеспечить ту же функциональность, но более высокую производительность

вы можете связать это выражение с такой функцией:

>>> super_dot = lambda v, w: FB.dgemm(alpha=1., a=v.T, b=w.T, trans_b=True)

Ответ 2

Я просто хотел поставить это на SO, но этот запрос на перенос должен быть полезным и удалить необходимость отдельной функции для numpy.dot https://github.com/numpy/numpy/pull/2730 Это должно быть доступно в numpy 1.7

Тем временем я использовал приведенный выше пример для записи функции, которая может заменить точку numpy, независимо от порядка ваших массивов, и сделать правильный вызов fblas.dgemm. http://pastebin.com/M8TfbURi

Надеюсь, что это поможет,

Ответ 3

numpy.einsum - это то, что вы ищете:

numpy.einsum('ij, i, ik -> jk', A, Q, A)

Это не требует дополнительной памяти (хотя обычно einsum работает медленнее, чем операции BLAS)