python - numpy matrix trickery - 逆时矩阵之和

Question

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

其中每个 X_i 是 n x p , 还有 r他们中的一个 r x n x p名为 samples 的数组. 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-ness 的东西)，但是有一个 Python 循环，这让 numpy 仙女们哭了。

我还可以想象 reshape samples这样我就可以得到一个 A^-1 x 的数组使用 solve对于每个 x无需执行 Python 循环，但这会产生一个很大的辅助数组，这会浪费内存。

我可以使用一些线性代数或 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 循环，并且必须将所有结果加在一起，所花费的时间是 390 毫秒，是解决 200 个必须解决的系统中的每一个所花费时间的 200 多倍。如果循环和求和是免费的，您将获得不到 5% 的改进。也可能有一些调用 python 函数的开销，但无论您使用哪种语言编写代码，与求解方程式的实际时间相比，它可能仍然可以忽略不计。