NumPy矩阵技巧 - 逆矩阵与矩阵相乘求和

Question

我试图做以下事情，并重复直到收敛：

其中每个 X_i 都是 n x p，在一个名为 samples 的 r x n x p 数组中有 r 个。 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性质的东西），但接下来有一个Python循环，这会让numpy小精灵哭泣。

我也可以想象将samples重塑，以便我可以使用solve为每个x获取A^-1 x数组，而无需进行Python循环，但这会产生一个浪费内存的大型辅助数组。

是否有一些线性代数或numpy技巧可以同时实现三个最佳方案：没有显式逆矩阵，没有Python循环，没有大型辅助数组？或者我的最佳选择是用一种更快的语言实现带有Python循环的方法并调用它？（直接将其移植到Cython可能有所帮助，但仍然涉及大量Python方法调用；但也许不需要太麻烦地直接制作相关的blas/lapack例程。）

（事实证明，我实际上不需要最终的矩阵U和V——只需要它们的行列式，或者实际上只需要它们的Kronecker积的行列式。因此，如果有人有一个聪明的主意，可以做更少的工作，仍然可以得到行列式，那将不胜感激。）

- Danica

---

很好的问题。我的大脑今天不太正常，但我想建议你把数学部分从开头到结尾至少发布到math.stackexchange.com，以防你错过了一个显而易见的捷径。你是对的，感觉上似乎有一种利用spd矩阵属性的方法，但我也看不出来。 - YXD

@MrE 感谢您的建议；我也在那里发布了。 - Danica

Answer 1

除非有更为聪明的答案出现，如果我是你的话，我会让仙女们哭泣……

（注：此段内容暂无IT技术相关性）

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毫秒，即该时间的约200倍。即使循环和求和是免费的，你也只会获得不到5%的改进。可能还存在一些调用Python函数的开销，但无论你在哪种编程语言中编写代码，它都可能仍然可以忽略不计，因为实际解决方程所需的时间要远大于此。

- Jaime

---

嗯...说得好。我愚蠢地在一个非常大的r和非常小的n和p的情况下计算了einsum方法与solve方法的时间，当然Python循环开销会更加重要。明天我会在我的真实数据上尝试一下，并看看比较结果如何。 - Danica

事实证明，使用scipy.linalg.cho_solve进行Python循环足够快，可以满足我的需求。我仍然很好奇是否有算法加速的方法，因此我会保留math.SE问题的开放状态。 - Danica