Как рассчитать ковариационную матрицу массивов 3d numpy? - PullRequest
Новая
       49

Как рассчитать ковариационную матрицу массивов 3d numpy?

1 голос
/ 26 февраля 2020

У меня есть матрица A, в форме (N, D, 4). Сначала я вычисляю A транспонированный, A_t. Я хочу рассчитать произведение A_t раз A. Я хочу, чтобы полученная матрица имела форму (D, D), и произведение матриц было бы таким, как если бы последний вектор из 4 компонентов был числом. (Точечное произведение двух векторов является числом.) 

import numpy as np

N = 15
D = 98
A = np.random.random((N, D, 4))
A_t = np.zeros((D, N, 4))
for i in range(N):
    A_t[:, i] = A[i]

S = np.zeros((D, D))

for i in range(D):
    row = A_t[i]
    for j in range(D):
        col = A[:, j, :]
        val = 0
        for n in range(N):
            val += np.matmul(row[n], col[n])
        S[i][j] = val

print(A.shape)
print(A_t.shape)
print(S.shape)

1 Ответ

1 голос
/ 27 февраля 2020

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

A_t = np.swapaxes(A, 0, 1)

Это эквивалентно 

A_t = np.transpose(A, [0, 1, 2])

или 

A_t = A.transpose([0, 1, 2])

. Как это ни случается, ни то, ни другое не требуется для вашего текущего приложения. Чтобы понять почему, давайте поработаем с упрощенным примером: 

np.random.seed(42)
N = 4
D = 3
K = 2
A = np.random.randint(0, 10, (N, D, K))

В вашем внешнем l oop у вас есть row = A_t[i]. Но по определению вашего транспонирования это идентично row = A[:, i, :], что делает вашу жизнь намного проще, а транспонирование излишним.

Внутренний l oop суммирует некоторые точечные произведения: 

val = 0
for n in range(N):
    val += np.matmul(row[n], col[n])

Если вы помните определение точечного произведения, вы увидите, что вы делаете эквивалент 

np.sum(np.sum(row * col, axis=1), axis=0)

Внутренняя сумма является суммой-произведением в вашем l oop, а внешняя сумма это вычисление val. Суммирование по обоим измерениям по отдельности аналогично суммированию всего буфера за раз, поэтому мы можем сразу заменить внутренний l oop просто 

for i in range(D):
    for j in range(D):
        S[i][j] = np.sum(A[:, i, :] * A[:, j, :])

Вы можете упростить это с помощью np.dot, np.tensordot, np.einsum или просто эфирное вещание. Первые два излишне сложны, потому что вы действительно умножаете сумму на два измерения одновременно. np.einsum предлагает самое прямое решение в целом, но это менее простой перевод вашего кода.

Решение 1. Трансляция

Давайте начнем с прямой версии вещания двойного l oop, прежде чем перейти к более идиоматическому c решению: 

S = (A[:, None, ...] * A[:, :, None, ...]).sum(axis=(0, -1))

или 

S = np.sum(A[:, None, ...] * A[:, :, None, ...], axis=(0, -1))

Это создает вид A в форме (N, 1, D, K) и (N, D, 1, K) соответственно. Умножение транслирует реплицированные оси D в каждом случае точно на то, что делают циклы for, поэтому итоговая сумма по осям N и K делает именно то, что делала линия S[i][j] = np.sum(A[:, i, :] * A[:, j, :]) ранее.

Решение 2: np.einsum

Это решение позволяет вам применять sum-product непосредственно к любым осям, которые вы хотите: 

S = np.einsum('ijk,ihk->jh', A, A)

Обратите внимание, что вы должны использовать другая буква для второй оси второй матрицы (j и h), чтобы указать, что вы будете не суммировать по этой оси. S симметрично, но если бы это было не так, вы могли бы транспонировать его, транспонируя в результат ->hj.

...