有没有一种快速的方法来在 Matlab 中反转矩阵？答案

【问题标题】：Is there a fast way to invert a matrix in Matlab?有没有一种快速的方法来在 Matlab 中反转矩阵？
【发布时间】：2011-09-24 10:57:08
【问题描述】：

我有很多大型（大约 5000 x 5000）矩阵需要在 Matlab 中求逆。我实际上需要倒数，所以我不能使用 mldivide 代替，这对于仅解决一个 b 的 Ax=b 来说要快得多。

我的矩阵来自一个问题，这意味着它们具有一些不错的属性。首先，它们的行列式是 1，所以它们绝对是可逆的。但是，它们不可对角化，或者我会尝试对它们进行对角化，反转它们，然后将它们放回原处。他们的条目都是实数（实际上是有理数）。

我正在使用 Matlab 来获取这些矩阵，为此我需要处理它们的逆矩阵，所以我更喜欢一种加快 Matlab 速度的方法。但是，如果我可以使用另一种更快的语言，请告诉我。我不懂很多其他语言（一点点C，一点点Java），所以如果它在其他语言中真的很复杂，那么我可能无法使用它。不过，请继续提出建议，以防万一。

【问题讨论】：

Matlab 的线性代数例程通常是非常优化的，因此您可能无法通过用另一种语言手动重新实现任何东西来提高速度。不过，您可能想查看 Eigen (c++)。对于 Matlab，您可能想尝试并行处理工具箱。
另外，如果您告诉我们您的矩阵具有哪些属性，或者您如何获得它们，这可能会有所帮助。也许在将其插入 MATLAB 之前进行一些线性代数操作可能比任何库都能为您提供更多帮助。
您可能想详细说明您的具体情况。你100%确定你需要逆吗？谢谢
我猜它们是稀疏的，因为大多数现实世界的问题都没有完整的 5000 x 5000 矩阵。假设它们是，你能以某种方式利用稀疏结构吗？你调查过这个吗？
矩阵有多稀疏？矩阵有什么特殊的形状，或者条目之间的关系吗？

标签： algorithm matlab matrix performance matrix-inverse

【解决方案1】：

我实际上需要逆，所以我不能用 mldivide 代替，...

这不是真的，因为您仍然可以使用mldivide 来获得相反的结果。请注意A-1 = A-1 * I。在 MATLAB 中，这相当于

invA = A\speye(size(A));

在我的机器上，5000x5000 矩阵大约需要 10.5 秒。请注意，MATLAB 确实有一个 inv 函数来计算矩阵的逆矩阵。尽管这将花费大约相同的时间，但在数值准确性方面效率较低（链接中的更多信息）。

首先，它们的行列式是 1，所以它们绝对是可逆的

而不是det(A)=1，而是condition number of your matrix 决定了逆运算的准确性或稳定性。请注意det(A)=∏i=1:n λi。所以只需设置λ1=M、λn=1/M 和λi≠1,n=1 就会给你det(A)=1。但是，如M → ∞、cond(A) = M2 → ∞ 和λn → 0，这意味着您的矩阵正在接近奇点，在计算逆时会有很大的数值误差。

我的矩阵来自一个问题，这意味着它们具有一些不错的属性。

当然，如果您的矩阵是稀疏的或具有其他有利的属性，则可以使用其他更有效的算法。但是如果没有关于您的具体问题的任何其他信息，就没有什么可以说的了。

我更喜欢加快 Matlab 速度的方法

MATLAB 使用高斯消元法计算使用 mldivide 的一般矩阵的逆矩阵（满秩、非稀疏、没有任何特殊属性），这是 Θ(n3)，其中 n 是矩阵的大小。因此，在您的情况下，n=5000 和 1.25 x 1011 浮点运算。因此，在具有大约 10 Gflops 计算能力的合理机器上，您将需要至少 12.5 秒来计算逆，并且没有办法解决这个问题，除非您利用“特殊属性”（如果它们是可利用的） )

【讨论】：

这真的有什么好处吗？
@yoda：你说你的方法需要 10.5 秒。我很想知道inv 处理相同的输入需要多长时间。
@PengOne:inv 大约需要相同的时间（约 10 秒），但不太准确。见inv，他们在哪里谈论这个。
一般指导是使用A\b而不是inv(A)*b来解决Ax=b当b是一个向量时。当b 是一个矩阵时——尤其是 yoda 建议的单位矩阵——我认为mldivide 形式没有任何好处。
@nibot, @yoda：求解程序通常执行 1) 矩阵分解（通常是某种形式的 LU 或 QR，或对称矩阵的 Cholesky）2) 反向代换（求解三角形或正交系统） 3）溶液抛光。固定矩阵求逆程序通常会跳过第 3 步（或进行完全高斯消除）。由于求解系统比矩阵求逆更频繁，因此矩阵求逆例程的优化程度往往较低，而且很可能比 1 + 2 + 3 与 n 右手边（1 必须完成一次）比一个没有人使用的通用反演程序。

【解决方案2】：

无论您使用什么语言，对任意 5000 x 5000 矩阵求逆在计算上都不容易。我建议研究近似值。如果你的矩阵是低秩的，你可能想尝试一个低秩近似 M = USV'

这里有一些来自 math-overflow 的更多想法：

https://mathoverflow.net/search?q=matrix+inversion+approximation

【讨论】：

嘘。仅仅因为矩阵具有 det(A)==1 并不意味着低秩近似不足以用于工程目的。参见例如A = 诊断（[1e5,1,1e-5]）。 A 的行列式仍然是 1，但如果您将 A 用于某些目的，它也可能是 diag([1e5 0 0])。
det A = 1这一事实并不意味着A在数值上是可逆的。对于大型矩阵，条件数很大的可能性大到足以使大多数矩阵不可逆（除非它们来自例如来自可逆运算符的离散化）。 SVD 方法可能是唯一明智的方法，除非 OP 为我们提供了问题的背景。
@AlexandreC.：什么是 SVD 方法？ SVD 代表什么？
@hellogoodbye ：这是帖子中的 M=USV' 方法（奇异值分解）

【解决方案3】：

首先假设特征值都是1。让A 成为矩阵的 Jordan 规范形式。然后您可以仅使用矩阵乘法和加法计算A^{-1}

A^{-1} = I + (I-A) + (I-A)^2 + ... + (I-A)^k

在哪里k < dim(A)。为什么这行得通？因为生成函数很棒。回忆一下展开

(1-x)^{-1} = 1/(1-x) = 1 + x + x^2 + ...

这意味着我们可以使用无限和来反转(1-x)。你想反转一个矩阵A，所以你想取

A = I - X

求解X 得到X = I-A。因此，通过替换，我们有

A^{-1} = (I - (I-A))^{-1} = 1 + (I-A) + (I-A)^2 + ...

在这里，我刚刚使用单位矩阵I 代替了数字1。现在我们要处理收敛问题，但这实际上不是问题。通过假设A 是Jordan 形式并且所有特征值等于1，我们知道A 是上三角形，所有1s 在对角线上。因此I-A 是上三角形，所有0s 在对角线上。因此I-A的所有特征值都是0，所以它的特征多项式是x^dim(A)，对于一些k < dim(A)，它的最小多项式是x^{k+1}。由于矩阵满足其最小（和特征）多项式，这意味着(I-A)^{k+1} = 0。因此，上述系列是有限，最大的非零项是(I-A)^k。所以它收敛了。

现在，对于一般情况，把你的矩阵变成 Jordan 形式，这样你就有一个块三角矩阵，例如：

A 0 0
0 B 0
0 0 C

每个块在对角线上都有一个值。如果A 的值是a，则使用上述技巧反转1/a * A，然后将a 相乘。由于完整矩阵是块三角形，因此逆矩阵是

A^{-1} 0      0
0      B^{-1} 0
0      0      C^{-1}

拥有三个块并没有什么特别之处，所以无论你有多少块都可以。

请注意，只要您有一个 Jordan 形式的矩阵，这个技巧就会起作用。在这种情况下，逆矩阵的计算在 Matlab 中会非常快，因为它只涉及矩阵乘法，您甚至可以使用技巧来加快计算速度，因为您只需要单个矩阵的幂。但是，如果将矩阵转换为 Jordan 形式真的很昂贵，这可能对您没有帮助。

【讨论】：

如果行列式是 1 并且条目是真实的，为什么特征值应该是 ±1？试试A=[1,3;2,7]。特征多项式为1 - 8*x + x^2。它的根源是 4±sqrt(15) 而不是 ±1。
@Yoda：不知怎的，我在想整数……好点。我将编辑我的答案。
@Alexandre C：它的哪一部分？可能我对一般情况很草率，但是当所有特征值都是1 时，这个解决方案确实有效。
@PengOne：您同时编辑了答案。您说“由于 det(A) = 1，特征值为 +-1”。除此之外，计算 Jordan 形式比计算逆矩阵要慢得多。我不知道 Jordan 形式的稳定性，但我怀疑它在很大程度上取决于 A 的条件数，因为矩阵很大，这可能会非常高。但是，如果所有特征值都是 1，则 A = I + N 具有 N 幂零，并且级数技巧有效（尽管在最坏的情况下使用 5000 项，这给出了与普通反演相同的 O(N^3) 复杂度）。
@Alexandre C：您的评论是在我编辑之后发表的，并指出整个答案（“this”）是错误的。是的，计算 Jordan 形式的成本可能很高，但基于“它们不可对角化，或者我会尝试对它们进行对角化，反转它们，然后将它们放回去”这一行。在这个问题中，我认为这种方法值得一提。