在Python中求解矩阵方程的有效方法答案

【问题标题】：Efficient way to solve matrix equation in Python在Python中求解矩阵方程的有效方法
【发布时间】：2018-07-02 19:46:42
【问题描述】：

现在我正在使用numpy.linalg.solve 来求解我的矩阵，但是我使用它来求解 5000*17956 矩阵的事实使得它非常耗时。它运行得很慢，我花了一个多小时才解决。求解矩阵方程的运行时间可能是 O(n^3)，但我从没想过会这么慢。有什么方法可以在 Python 中更快地解决它？

我的代码是这样的，求解方程BT * UT = BT*B a，其中m 是测试用例的数量（在我的情况下超过5000），B 是数据矩阵m*17956，并且u 是 1*m。

C = 0.005                        # hyperparameter term for regulization
I = np.identity(17956)          # 17956*17956 identity matrix
rhs = np.dot(B.T, U.T)          # (17956*m) * (m*1)     = 17956*1
lhs = np.dot(B.T, B)+C*I        # (17956*m) * (m*17956) = 17956*17956
a = np.linalg.solve(lhs, rhs)   # B.T u = B.T B a, solve for a (17956*1)

【问题讨论】：

我不是矩阵专家，但是为什么不能将每边乘以 B.T^-1 以去除 B.T 项？
求逆矩阵的计算效率很低

标签： python performance numpy matrix

【解决方案1】：

更新（2018 年 7 月 2 日）：更新后的问题询问了正则化项的影响以及矩阵中的数据类型。一般来说，这会对特定 CPU 最优化的数据类型产生很大的影响（粗略的经验法则是，当所有其他事情都是保持相等），并且大量零值的存在可以允许使用稀疏矩阵库。但在这种特殊情况下，主对角线上的变化（远低于所有考虑值的 1%）对运行时间的影响可以忽略不计。

TLDR;

一个小时是合理的（三次回归表明这在我的机器上大约需要 83 分钟 - 一台低端 chromebook）。
生成lhs 和rhs 的预处理几乎不占任何时间。
您将无法比使用numpy.linalg.solve 更快地解决那个确切的问题。
如果m 像您建议的那样小，并且如果B 是可逆的，您可以在一分钟或更短的时间内解出方程U.T=Ba。
如果这是一个更大问题的一部分，这个代价高昂的中间步骤也许可以从数学框架中简化出来。
性能瓶颈确实应该通过分析来解决，以找出导致问题的步骤。
由于这来自真实世界的数据，根据最终目标，您可能能够使用更少的特征（直接或通过 PCA、NMF 或 LLE 等缩减步骤）。
正如另一个答案中提到的，如果矩阵足够稀疏，您可以使用稀疏线性代数例程来获得很好的效果（许多自然语言处理数据源都是这样的）。
由于输出是一维向量，我将使用np.dot(U, B).T 而不是np.dot(B.T, U.T)。转置很整洁。这避免了在像B 这样的大矩阵上进行转置，但由于您将三次运算作为主要步骤，这对您的问题并不重要。
根据您是否需要原始数据以及所涉及的矩阵是否有任何其他特殊属性，您可能可以修改scipy.linalg.solve 中的参数以获得收益。李>
我在用块矩阵方程替换大型矩阵方程时取得了喜忧参半的成功，这些方程又回到了 numpy 例程上。这种方法通常比 numpy 方法节省 5-20%，并在我的系统上节省 1% 左右的 scipy 方法。我还没有完全探究造成这种差异的原因。

【讨论】：

但是现在已经过去了几个小时，它仍在计算中。
@JamesLiu 这有点奇怪。情况如何？ CPU规格？内存规格？ CPU使用率？内存使用情况？操作系统？ Python版本？ ETC...？代码很好。对于诊断，我们没有太多其他工作要做。
我使用的是python 3.7 64位，我认为它是最新版本。我使用的是 2.8 Ghz 的 7700hq，我认为这对于这个计算来说绝对绰绰有余。这与 B 和 U 的值有关吗？我的 u 仅包含 -1 和 1，而我的 U 包含 0 或 1 值。为了摆脱无法为方程提供唯一解的奇点，我在 BT*U 中添加了一个惩罚项 CI，其中 c 是一个非常小的系数，而 I 是单位矩阵。这可能与计算复杂性有关吗？我认为这无关紧要
这段代码之前或之后是否还有其他运行？在其他条件相同的情况下，我只运行了您发布的代码，无需等待几个小时。我的电脑应该比你的慢很多。它本身并不快，因为您要进行大约 6-18 个 terra 操作，但也不会花费 3 多个小时。
我之前已经测试过代码，包括读取输入只花了不到一分钟左右的时间。这真的很奇怪。我将重新运行代码以查看其余部分花费了多长时间。

【解决方案2】：

假设您的矩阵是稀疏的，scipy.sparse.linalg 模块将很有用。 Here 是整个模块的文档，here 是 spsolve 的文档。

【讨论】：