NumPy와 매트릭스 속임수 - 역 번 행렬의 합

Question

나는 다음을 수행하고, 수렴 할 때까지 반복하기 위해 노력하고있어 :

각 X _내가는 n x p이며, 그 중 r가있는 경우 samples이라는 r x n x p 배열에 있습니다. U은 n x n, V은 p x p입니다. (나는 matrix normal distribution의 MLE를 얻고있다.) 크기는 모두 잠재적으로 크다. 나는 적어도 r = 200, n = 1000, p = 1000의 순서로 물건을 기대하고 있습니다.

내 현재 코드는

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) 

이 괜찮 작동하지만, 물론 당신은 실제로 역을 찾아 그것에 의해 물건을 곱하면 안되는 않습니다. U와 V가 대칭적이고 긍정적 인 사실을 어떻게 든 활용할 수 있다면 좋을 것입니다. 반복에서 U와 V의 Cholesky 팩터를 계산할 수 있기를 원하지만, 합계 때문에이를 수행하는 방법을 모르겠습니다.

나는

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

(또는 PSD-다움을 악용 그 비슷한)과 같은 작업을 수행하여 역을 피하기,하지만 거기에 파이썬 루프, 그리고 그는 NumPy와 요정 울려요 수 있습니다.

나는 또한 내가 파이썬 루프를 할 필요없이 모든 x에 대한 solve를 사용 A^-1 x의 배열을 얻을 수있는 방식으로 samples을 재편 상상할 수 있지만, 메모리의 낭비 큰 보조 배열을 만든다.

명시 적 반전 없음, 파이썬 루핑 없음, 큰 aux 배열 없음을 모두 만족하기 위해 할 수있는 선형 대수 또는 numpy 트릭이 있습니까? 아니면 더 빠른 언어로 파이썬 루프를 구현하고 호출하는 것이 최선의 방법일까요? (Cython으로 직접 이식하는 것만으로도 도움이 될 수 있지만 여전히 많은 Python 메서드 호출이 필요하지만 어쨌든 너무 많은 문제없이 관련 blas/lapack 루틴을 직접 작성하는 것은 그리 어렵지 않을 수 있습니다.)

(실제로 밝혀진 바 최종적으로 행렬 U과 V은 실제로 필요하지 않으며 단지 결정 요인이거나 크로 네커 제품의 결정 요인 일뿐입니다. 따라서 작업량을 줄이는 방법에 대한 영리한 아이디어가 있다면 결정 요인을 얻을, 그것은 많이 주시면 감사하겠습니다.) 누군가가 더 영감 대답과 함께 때까지 내가 당신이라면

Answer 1

가, 내가 요정 울고하자 것 ...

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 

그래서 파이썬 루프를 사용하고 모든 결과를 종합해야하는 경우 해결해야하는 200 개 시스템 각각을 해결하는 데 소요되는 시간이 390 밀리 초가 넘습니다. 루핑과 합계가 무료 인 경우 5 % 미만의 개선 효과를 얻습니다. 파이썬 함수 오버 헤드를 호출하는 호출이있을 수 있습니다.하지만 방정식을 푸는 실제 시간에 대해서는 무시할 수 있습니다. 코드를 작성한 언어가 무엇이든 상관 없습니다.