Диагональ матрицы numpy без вычисления всей матрицы

#python #numpy

Вопрос:

У меня есть простая алгебрическая проблема, и я хотел бы решить ее с помощью numpy (конечно, я мог бы легко решить ее с помощью numba, но дело не в этом).

Давайте рассмотрим первую случайную матрицу A с размером (m x n), с n большим значением, и вторую случайную матрицу B с размером (n x n).

 A = np.random.random((1E6, 1E2))
B = np.random.random((1E2, 1E2))

Мы хотим вычислить следующее выражение:

 np.diag(np.dot(np.dot(A,B),B.T))

Проблема в том, что вся матрица загружается в память и только затем извлекается диагональ. Можно ли выполнить эту операцию более эффективным способом?

1. Являются ли входные матрицы действительно случайными или это только примеры?

2. Это упражнение по стекированию, а не реальная проблема. Я выбираю случайные матрицы в качестве наиболее общего случая.

Ответ №1:

Вот как я бы подошел к этому с вашего начального выражения

 np.diag(np.dot(np.dot(A,B),B.T))

Вы можете начать с группировки терминов:

 np.diag(np.dot(A, np.dot(B,B.T)))

затем используйте только первую соответствующую (квадратную) часть:

 np.diag(np.dot(A[:B.shape[0], :], np.dot(B,B.T)))

а затем избегайте дополнительных умножений (которые выпадут из диагонали), выполняя умножение по элементам самостоятельно:

 np.sum( np.multiply(A[:B.shape[0], :].T, np.dot(B,B.T)), 0)

1. Уоу! Мне просто нужно было вытащить бумагу, чтобы понять, почему np.sum(A.T * B, 0) == np.diag(A @ B) . С ума сойти

2. Поскольку каждый элемент A @ B является точечным продуктом строки A и столбца B, я транспонирую A и выполняю точечные произведения строк (путем поэлементного умножения и горизонтальной суммы). Это даст вам только те элементы, которые будут падать по диагонали. Это всего лишь увеличение ~2,5 x, собственный продукт matrix-matrix хорошо оптимизирован.

3. Да. Я только что проверил… вы можете пропустить создание промежуточного массива с помощью np.einsum('ji,ij->j', ...) и сократить таким образом несколько микросекунд, но я не уверен, что это того стоит.

4. H=np.squeeze(A[:B.shape[0],None,:]@(B@B.T)[:,:,None]) это еще быстрее, но труднее понять.

5. @hpaulj :О, это действительно быстрее, удивительный трюк.

Ответ №2:

Изменено (A*B)*B.T на A*(B*B.T)
Умножьте только эту часть A ( A[:B.shape[0]] ), что приведет к диагональной части матрицы

 import numpy as np
import time

A = np.random.random((1000_000, 100))
B = np.random.random((100, 100))

start_time = time.time()
result = np.diag(np.dot(np.dot(A, B), B.T))
print('Baseline: ', time.time() - start_time)

start_time = time.time()
for i in range(100):
    result2 = np.diag(np.dot(A[:B.shape[0]], np.dot(B, B.T)))
print('Optimized: ', (time.time() - start_time) / 100)
stop = 1

assert np.allclose(result, result2)

 Baseline:  1.7957241535186768
Optimized:  0.00016015291213989258

Ответ №3:

ДА.

 N = 1E6

A = np.random.random((N, 1E2))
B = np.random.random((1E2, 1E2))

result = 0;
for i in range(N):
  result  = np.dot(np.dot(A[i,:], B[i,:])[i, :], B.T[i, :])

  # Replacing B.T[i, :] with B[:, i].T might be a little more efficient

Объяснение:

Скажем, у нас есть: K = np.dot(np.dot(A,B),B.T) .

Затем, K[0,0] = (A[0, :] * B[:,0])[0, :] * B.T[:])

Позвольте X = (A[0, :] * B[:,0]) , что является [0, 0] элементом np.dot(A,B)
Тогда X[0, :] * B.T[:, 0] является [0, 0] элементом np.dot(np.dot(A,B),B.T)
Тогда X[0, :] * B.T[:, 0] = (A[0, :] * B[:,0])[0, :] * B.T[:])

Мы также можем обобщить этот результат на: K[i,i] = (A[i, :] * B[:,i])[i, :] * B.T[:, i])