Вычислить «v^TA v» для матрицы векторов v

Question

У меня есть k*n матрица X, и k*k матрицы A. Для каждого столбца X, я хотел бы, чтобы вычислить скалярное

X[:, i].T.dot(A).dot(X[:, i]) 

(или, математически, Xi' * A * Xi).

В настоящее время у меня есть for цикла:

out = np.empty((n,)) 
for i in xrange(n): 
    out[i] = X[:, i].T.dot(A).dot(X[:, i]) 

но так n велик, я хотел бы сделать это быстрее, если это возможно (то есть, используя некоторые функции Numpy вместо цикла).

Answer 1

Это, кажется, сделать это красиво: (X.T.dot(A)*X.T).sum(axis=1)

Edit: Это немного быстрее. np.einsum('...i,...i->...', X.T.dot(A), X.T). Оба работают лучше, если X и A являются фортранами непрерывными.

Answer 2

Вы не можете сделать это быстрее, если вы не распараллелите все это: один поток за столбец. Вы по-прежнему будете использовать циклы - вы не можете уйти от этого.

Сокращение карты - это хороший способ взглянуть на эту проблему: умножить кратные шаги, уменьшить суммы шага.

Answer 3

Вы можете использовать numpy.einsum:

np.einsum('ji,jk,ki->i',x,a,x) 

Это позволит получить тот же результат. Давайте посмотрим, если это намного быстрее:

Похоже dot еще быстрее вариант, особенно потому, что он использует резьбовую BLAS, в отличие от einsum, который работает на одном ядре.

import perfplot 
import numpy as np 


def setup(n): 
    k = n 
    x = np.random.random((k, n)) 
    A = np.random.random((k, k)) 
    return x, A 


def loop(data): 
    x, A = data 
    n = x.shape[1] 
    out = np.empty(n) 
    for i in range(n): 
     out[i] = x[:, i].T.dot(A).dot(x[:, i]) 
    return out 


def einsum(data): 
    x, A = data 
    return np.einsum('ji,jk,ki->i', x, A, x) 


def dot(data): 
    x, A = data 
    return (x.T.dot(A)*x.T).sum(axis=1) 


perfplot.show(
    setup=setup, 
    kernels=[loop, einsum, dot], 
    n_range=[2**k for k in range(10)], 
    logx=True, 
    logy=True, 
    xlabel='n, k' 
    )