计算矩阵行列式的最快算法是什么？

Question

我被指派研究计算矩阵行列式的最快算法，这与IT技术有关。

我已经了解到LU分解和Bareiss算法，它们都运行在O(n^3)的时间复杂度，但是通过一些挖掘，似乎存在一些时间复杂度介于n^2和n^3之间的算法。

这个来源（见第113-114页）和这个来源（见第198页）说，基于Coppersmith-Winograd算法进行矩阵乘法，存在一个运行时间为O(n^2.376)的算法。然而，我找不到任何关于这种算法的详细信息。

我的问题是：

1. 计算矩阵行列式的已知最快算法是什么？
2. 在哪里可以找到关于这个最快算法的信息？

非常感谢。

回答：

1. 已知的计算矩阵行列式的最快算法是基于Coppersmith-Winograd算法的，其运行时间为O(n^2.376)。
2. 目前没有找到有关该算法的详细信息。

- John-Luke Laue

---

矩阵有多大？你想计算多少个行列式？ - Absurd-Mind

我会假设矩阵非常大（N > 22应该足够大了吧？）。那么有多少个呢？对于给定的矩阵只有一个行列式。 输入：1个大矩阵 输出：输入矩阵的单个行列式。 - John-Luke Laue

数值稳定性也是一个问题吗？ - Henry

值得探究Ladderman乘法的问题和答案，https://dev59.com/CWgv5IYBdhLWcg3wI9Wg。 - Hooked

@Henry 抱歉，我不知道什么是数值稳定性。 - John-Luke Laue

@jlguitar287 http://en.wikipedia.org/wiki/Numerical_stability - Henry

Answer 1

我相信实践中（并且常用）最快的算法是 Strassen 算法。您可以在 维基百科 上找到解释以及示例 C 代码。

基于 Coppersmith-Winograd 乘法算法 的算法过于复杂，不适合实际使用，尽管它们具有最佳渐近复杂度。

- Sopel

---

值得注意的是，目前的上限为O(n^{2.3728639})，但正如Sopel所说，由于巨大的常数因子，追求较小的指数在实际计算中可能并不切实际。 - Hooked

@Sopel 我认为 Strassen 算法用于矩阵求逆。参见http://en.wikipedia.org/wiki/Computational_complexity_of_mathematical_operations至于基于 Coppersmith-Winograd 的算法，你是正确的 - 我问了我的教授，他说它们是最快的，但过于复杂而不实用。谢谢！ - John-Luke Laue

我相信它们是非常相似的问题（我是从Coppersmith-Winograd算法可转换为计算行列式中得出的结论），但我对此并没有很大的了解，所以我可能是错误的。 - Sopel

你能解释一下矩阵乘法算法如何用于计算行列式吗？ - HelloGoodbye

Answer 2

我知道这不是对我的问题的直接回答，但是为了完成我的研究论文，这已经足够了。我最终向我的教授提出了问题，现在我将总结他的回答：

总结：

• 最快的矩阵乘法算法（例如Coppersmith-Winograd和更近期的改进）可以使用O(n^~2.376)个算术运算，但使用复杂的数学工具，通常不实用。
• LU分解和Bareiss确实使用O(n^3)个算术运算，但是更实用。

简而言之，尽管LU分解和Bareiss不如最有效的算法快，但它们更实用，我应该把研究论文重点放在这两个方面。

感谢所有评论和帮助！

- John-Luke Laue

---

“LU分解和Bareiss算法不够快”这个说法有点误导。你想要表达的是“从渐近复杂度上来看，LU分解和Bareiss算法并不够快”。实际上，在执行时间函数中，渐近复杂度更高的算法可能会有很高的常数项，导致在实践中速度更慢的算法反而更快。例如，对于所有n < 1e5，1e6n比0.0001n^2更大。 - Gene

你可以更具体地说：“当n < 一些大常数”时，“LU分解和Bareiss方法比Coppersmith-Winograd方法更快地找到nxn矩阵的行列式”。当然，这需要一些工作来找出这个大常数。 - Stef

Answer 3

请看下面的Matlab测试脚本，它可以计算任意方阵的行列式（还包括与Matlab内置函数的比较）：

nMin = 2;  % Start with 2-by-2 matrices
nMax = 50; % Quit with 50-by-50 matrices
nTests = 10000;
detsOfL = NaN*zeros(nTests, nMax - nMin + 1);
detsOfA = NaN*zeros(nTests, nMax - nMin + 1);
disp(' ');
for n = nMin:nMax
    tStart = tic;
    for j = 1:nTests
       A = randn(n, n);
       detA1 = det(A); % Matlab's built-in function
       if n == 1
           detsOfL(j, 1) = 1;
           detsOfA(j, 1) = A;
           continue; % Trivial case => Quick return
       end
       [L, U, P] = lu(A);
       PtL = P'*L;
       realEigenvaluesOfPtL = real(eig(PtL));
       if min(prod(realEigenvaluesOfPtL)) < 0 % det(L) is always +1 or -1
           detL = -1;
       else
           detL = 1;
       end
       detU = prod(diag(U));
       detA2 = detL * detU; % Determinant of A using LU decomposition
       if detA1 ~= detA2
           error(['Determinant computation failed at n = ' num2str(n) ', j = ' num2str(j)]);
       end
       detsOfL(j, n - nMin + 1) = detL;
       detsOfA(j, n - nMin + 1) = detA2;
    end
    tElapsed = toc(tStart);
    disp(sprintf('Determinant computation was successful for n = %d!! Elapsed time was %.3f seconds', n, tElapsed));
end
disp(' ');

- vnormi75

---

只是简单地倾倒代码并不能构成一个答案。(0) 这是你写的吗？它的许可证是什么？如果不是，那么你从哪里得到了它？现在来真正回答问题：(1) 它与其他解决方案相比在性能上如何，有明确的量化数据吗？(2) 它为什么能够如此出色地表现的详细信息是什么？ - underscore_d