Tính "v^TA v" cho một ma trận của vectơ v

Question

Tôi có một k*n ma trận X, và một k*k ma trận A. Đối với mỗi cột của X, tôi muốn để tính toán vô hướng

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

(hoặc, toán học, Xi' * A * Xi).

Hiện nay, tôi có một vòng lặp for:

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

nhưng vì n là lớn, tôi muốn thực hiện điều này nhanh hơn nếu có thể (ví dụ: sử dụng một số chức năng NumPy thay vì một vòng lặp).

Answer 1

Điều này dường như làm nó độc đáo: (X.T.dot(A)*X.T).sum(axis=1)

Chỉnh sửa: Đây là nhanh hơn một chút. np.einsum('...i,...i->...', X.T.dot(A), X.T). Cả hai hoạt động tốt hơn nếu X và A là Fortran liền kề nhau.

Answer 2

Bạn không thể thực hiện nhanh hơn trừ khi bạn song song toàn bộ nội dung: Một chuỗi cho mỗi cột. Bạn vẫn sẽ sử dụng vòng lặp - bạn không thể thoát khỏi điều đó.

Giảm bản đồ là một cách hay để xem xét vấn đề này: bội số bước bản đồ, giảm tổng số bước.

Answer 3

Bạn có thể sử dụng numpy.einsum:

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

này sẽ nhận được kết quả tương tự. Chúng ta hãy xem nếu nó là nhanh hơn nhiều:

Hình như dot vẫn là lựa chọn nhanh hơn, đặc biệt là vì nó sử dụng ren BLAS, như trái ngược với einsum mà chạy trên một lõi.

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