numpy矩阵欺骗 - 逆矩阵和的总和

Question

我正在尝试执行以下操作,并重复直到收敛:

每个X _i的位置n x p,并且r它们在一个r x n x p名为的数组中samples.U是n x n,V是p x p.(我得到了矩阵正态分布的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-ness的东西),但随后有一个Python循环,这使得那些笨拙的仙女哭了.

我还可以想象samples以这样一种方式重塑:我可以为每一个A^-1 x使用数组而不必进行Python循环,但这会产生一个大的辅助数组,这会浪费内存.solvex

是否有一些线性代数或numpy技巧可以做到最好的三个:没有明确的反转,没有Python循环,没有大的辅助数组？或者我最好用更快的语言实现Python环路并调用它？(直接将它移植到Cython可能会有所帮助,但仍然会涉及很多Python方法调用;但也许直接制作相关的blas/lapack例程并不会太麻烦.)

(事实证明,我实际上并不需要矩阵U,V最后只是他们的决定因素,或者实际上只是他们的Kronecker产品的决定因素.所以,如果有人对如何减少工作并且仍然得到决定因素,非常感谢.)

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

因此,拥有一个python循环,并且必须将所有结果汇总在一起,所需的时间比解决每个必须解决的200个系统所需的时间长390毫秒.如果循环和求和是免费的,你得到的改善不到5%.可能还有一些调用python函数的开销,但无论你用什么语言编写代码,它对解决方程的实际时间可能仍然可以忽略不计.