numpy 矩阵技巧 - 逆时间矩阵的总和

Posted 2023-03-12

【中文标题】numpy 矩阵技巧 - 逆时间矩阵的总和

【英文标题】：numpy matrix trickery - sum of inverse times matrices

【发布时间】：2013-01-24 20:34:36

【问题描述】：

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

其中每个 X_i 是 n x p，并且在名为 samples 的 r x n x p 数组中存在 r。 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 的东西），但是有一个 Python 循环，这让 numpy 仙女哭了。

我还可以想象重塑 samples 的方式，即我可以使用 solve 为每个 x 获得一个 A^-1 x 数组，而无需执行 Python 循环，但这会产生一个很大的辅助数组浪费内存。

我是否可以做一些线性代数或 numpy 技巧来充分利用这三者：没有显式逆运算、没有 Python 循环以及没有大的辅助数组？或者我最好的选择是用一种更快的语言实现一个 Python 循环并调用它？ （直接将其移植到 Cython 可能会有所帮助，但仍会涉及大量 Python 方法调用；但直接制作相关的 blas/lapack 例程可能不会太麻烦。）

（事实证明，我实际上并不需要矩阵 U 和 V 最后 - 只是它们的行列式，或者实际上只是他们的 Kronecker 乘积的行列式。所以如果有人有一个聪明的想法如何做更少的工作并仍然得到决定因素，这将不胜感激。）

【问题讨论】：

写得很好的问题。今天我的大脑运转不佳，但我只是想建议您至少将数学部分从头到尾发布到 math.stackexchange.com，以防您遗漏了明显的捷径。你是对的，它感觉就像可能是一种利用spd矩阵属性的方法，但我也看不到它。

@MrE 感谢您的建议； I've posted it there also.

【参考方案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 函数的开销，但对于求解方程的实际时间，它可能仍然可以忽略不计，无论您使用哪种语言编写代码。

【讨论】：

嗯...好点。在r 非常大且n 和p 非常小的情况下，我愚蠢地对einsum 方法与solve 方法进行了计时，当然，Python 循环开销是有意义的重要得多。明天我会在我的真实数据上尝试一下，看看比较是什么。

事实证明，使用scipy.linalg.cho_solve 执行 python 循环对我的需求来说已经足够快了。我仍然很好奇是否有算法加速，所以我把 math.SE 问题留了下来。