numpy 矩阵技巧 - 逆矩阵乘法之和

Question

我正在尝试做以下事情，并重复直到结果稳定：

这里每个 X_i 是一个 n x p 的矩阵，而总共有 r 个这样的矩阵，组成一个 r x n x p 的数组，叫做 samples。U 是一个 n x n 的矩阵，V 是一个 p x p 的矩阵。（我在计算一个 矩阵正态分布 的最大似然估计。） 这些矩阵的大小都可能比较大；我预计 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)

（或者其他利用正定特性的方式），但这样就会出现 Python 循环，这样会让 numpy 的性能下降。

我也想过重新调整 samples 的形状，这样我就可以用 solve 为每个 x 计算 A^-1 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函数也会有一些额外的时间开销，但无论你用什么编程语言，解决方程的实际时间都还是会占大头。