Есть ли «улучшенный» метод «тупой / тупой»?

Question

Задача

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

Y = A**T * Q * A

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

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

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

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

1. Y = np.dot(np.dot(A.T, np.diag(Q)), A) и
2. 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 вместе с ним?

Answer 1

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

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

>>> 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, передайте массив и установите соответствующий флаг order в 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

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

Но почему массивы копируются перед передачей в dot ?

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

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

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

Таким образом, чтобы вычислить точку (А.Т., А) без создания дополнительной копии:

>>> 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)

Answer 2

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

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

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

Answer 3

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

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

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