NumPy ma trận thủ đoạn gian trá - tổng của ma trận nghịch đảo lần

Question

Tôi đang cố gắng để làm như sau, và lặp lại cho đến khi hội tụ:

nơi mỗi X _i là n x p, và có r trong số họ trong một mảng r x n x p được gọi là samples. U là n x n, V là p x p. (Tôi nhận được MLE của một matrix normal distribution.) Các kích thước là tất cả khả năng lớn-ish; Tôi đang mong đợi những thứ ít nhất theo thứ tự của r = 200, n = 1000, p = 1000.

mã hiện tại của tôi không

V = np.einsum('aji,jk,akl->il', samples, np.linalg.inv(U)/(r*n), samples) 
U = np.einsum('aij,jk,alk->il', samples, np.linalg.inv(V)/(r*p), samples) 

này hoạt động ổn, nhưng tất nhiên bạn sẽ không bao giờ phải thực sự tìm nghịch đảo và nhân thứ bởi nó. Nó cũng sẽ tốt nếu tôi bằng cách nào đó có thể khai thác thực tế là U và V là đối xứng và xác định tích cực. Tôi rất muốn có thể tính toán hệ số Cholesky của U và V trong lần lặp lại, nhưng tôi không biết làm thế nào để làm điều đó vì tổng số.

tôi có thể tránh được nghịch đảo bằng cách làm một cái gì đó giống như

V = sum(np.dot(x.T, scipy.linalg.solve(A, x)) for x in samples) 

(hoặc một cái gì đó tương tự mà khai thác psd-Ness), nhưng sau đó có một vòng lặp Python, và điều đó làm cho các nàng tiên NumPy khóc.

tôi cũng có thể tưởng tượng được định hình lại samples theo cách như vậy mà tôi có thể nhận được một loạt các A^-1 x sử dụng solve cho mỗi x mà không cần phải làm một vòng lặp Python, nhưng điều đó làm cho một mảng lớn phụ trợ đó là một sự lãng phí bộ nhớ. Có một số đại số tuyến tính hoặc mẹo vặt tôi có thể làm để có được tốt nhất trong cả ba: không có đảo ngược rõ ràng, không có vòng lặp Python, và không có mảng aux lớn nào không? Không. Hoặc là đặt cược tốt nhất của tôi thực hiện một với một vòng lặp Python trong một ngôn ngữ nhanh hơn và gọi ra nó? (Chỉ cần chuyển trực tiếp đến Cython có thể giúp ích, nhưng vẫn liên quan đến nhiều cuộc gọi phương thức Python; nhưng có thể sẽ không quá nhiều rắc rối để thực hiện các thủ tục liên quan đến blas/lapack trực tiếp mà không gặp quá nhiều rắc rối.)

(Khi nó quay ra, tôi không thực sự cần các ma trận U và V cuối cùng - chỉ là yếu tố quyết định của chúng, hoặc thực sự chỉ là yếu tố quyết định cho sản phẩm Kronecker của chúng. Vì vậy, nếu ai đó có ý tưởng thông minh về cách làm ít công việc hơn có được yếu tố quyết định ra ngoài, mà sẽ được nhiều đánh giá cao.)

Answer 1

cho đến khi một người nào đó đi lên với một câu trả lời đầy cảm hứng hơn, nếu tôi là bạn, tôi muốn để cho các nàng tiên khóc ...

r, n, p = 200, 400, 400 

X = np.random.rand(r, n, p) 
U = np.random.rand(n, n) 

In [2]: %timeit np.sum(np.dot(x.T, np.linalg.solve(U, x)) for x in X) 
1 loops, best of 3: 9.43 s per loop 

In [3]: %timeit np.dot(X[0].T, np.linalg.solve(U, X[0])) 
10 loops, best of 3: 45.2 ms per loop 

Vì vậy, có một vòng python, và phải tổng hợp tất cả các kết quả với nhau, lấy 390 ms nhiều hơn 200 lần những gì nó cần để giải quyết từng trong số 200 hệ thống phải được giải quyết. Bạn sẽ nhận được ít hơn một cải tiến 5% nếu vòng lặp và tổng hợp được miễn phí. Có thể có một số chức năng gọi hàm python trên đầu cũng có, nhưng nó có lẽ vẫn sẽ không đáng kể so với thời gian thực giải phương trình, bất kể ngôn ngữ bạn mã hóa nó là gì.