Introduction

行列式(Determinant) 是矩阵的重要属性。

在手动计算行列式时，我们常常使用两种方法：

• 按行/列进行拉普拉斯展开。
• 利用矩阵在任意行/列加减其他行列的任意倍后行列式不变的性质，化为三角矩阵后，计算主对角线元乘积求解。

前者的复杂度是 $O(n!)$ 级别的，在计算约 $12$ 阶的矩阵时就会需要超过 1s 的时间，而计算 $1000 \times 1000$ 的矩阵需要进行约：

402387260077093773543702433923003985719374864210714632543799910429938512398629020592044208486969404800479988610197196058631666872994808558901323829669944590997424504087073759918823627727188732519779505950995276120874975462497043601418278094646496291056393887437886487337119181045825783647849977012476632889835955735432513185323958463075557409114262417474349347553428646576611667797396668820291207379143853719588249808126867838374559731746136085379534524221586593201928090878297308431392844403281231558611036976801357304216168747609675871348312025478589320767169132448426236131412508780208000261683151027341827977704784635868170164365024153691398281264810213092761244896359928705114964975419909342221566832572080821333186116811553615836546984046708975602900950537616475847728421889679646244945160765353408198901385442487984959953319101723355556602139450399736280750137837615307127761926849034352625200015888535147331611702103968175921510907788019393178114194545257223865541461062892187960223838971476088506276862967146674697562911234082439208160153780889893964518263243671616762179168909779911903754031274622289988005195444414282012187361745992642956581746628302955570299024324153181617210465832036786906117260158783520751516284225540265170483304226143974286933061690897968482590125458327168226458066526769958652682272807075781391858178889652208164348344825993266043367660176999612831860788386150279465955131156552036093988180612138558600301435694527224206344631797460594682573103790084024432438465657245014402821885252470935190620929023136493273497565513958720559654228749774011413346962715422845862377387538230483865688976461927383814900140767310446640259899490222221765904339901886018566526485061799702356193897017860040811889729918311021171229845901641921068884387121855646124960798722908519296819372388642614839657382291123125024186649353143970137428531926649875337218940694281434118520158014123344828015051399694290153483077644569099073152433278288269864602789864321139083506217095002597389863554277196742822248757586765752344220207573630569498825087968928162753848863396909959826280956121450994871701244516461260379029309120889086942028510640182154399457156805941872748998094254742173582401063677404595741785160829230135358081840096996372524230560855903700624271243416909004153690105933983835777939410970027753472000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

次运算。

这样计算行列式的效率显然是极低的。而通过化三角矩阵，我们可以用 $O(n^3)$ 的复杂度完成行列式的求解。对于同样的矩阵，我们只需要进行 $1 \times 10^9$ 的运算。这对于中小规模的矩阵已经足够快速了。

Prerequisites

对于矩阵 $\mathbf{A}$，记 $\begin{vmatrix}\mathbf{A}\end{vmatrix}$ 为其行列式的值。

在进行算法实现之前，我们需要对用到的性质做一定了解。
本文不会对性质进行证明，可以自行参考相关教材。

行列式消元：

$$ \begin{vmatrix} a_{1,1} & a_{1,2} & \cdots & a_{1n} \ \vdots & \vdots & & \vdots \a_{i,1} + ka_{j,1} & a_{i,2} +ka_{j,2}& \cdots & a_{i,n}+ka_{j,n} \\vdots & \vdots & \vdots & \a_{n,1} & a_{n,2} & \cdots & a_{n,n}\end{vmatrix} = \begin{vmatrix} a_{1,1} & a_{1,2} & \cdots & a_{1n} \ \vdots & \vdots & & \vdots \a_{i,1} & a_{i,2}& \cdots & a_{i,n} \\vdots & \vdots & \vdots & \a_{n,1} & a_{n,2} & \cdots & a_{n,n}\end{vmatrix} \tag{1} $$

三角矩阵的行列式计算：

$$ \begin{aligned} &{\mathbf {U}}={\begin{bmatrix}u_{{1,1}}&u_{{1,2}}&u_{{1,3}}&\ldots &u_{{1,n}}\&u_{{2,2}}&u_{{2,3}}&\ldots &u_{{2,n}}\\vdots &&\ddots &\ddots &\vdots \&(0)&&\ddots &u_{{n-1,n}}\0&&\cdots &&u_{{n,n}}\end{bmatrix}} \ &\det (\mathbf{U}) = \prod \limits {i=1}^n u{i,i}
\end{aligned} \tag{2} $$

$$ \text{互换矩阵的任意两行，行列式变号。} \tag{3} $$

$$ \text{矩阵中某行或某列全为零时，行列式为零。} \tag{4} $$

如果你了解高斯消元相关的内容，那再好不过了。

Theory

通过性质 1，我们可以对矩阵进行变换，将其化为三角矩阵，从而通过性质 2 的方法求解行列式。

先从一个具体的例子入手。

$$ \mathbf{A} = \begin{bmatrix}1 & 1 & 1 \ 1 & 3 & 2 \ 1 & 5 & 7\end{bmatrix} $$

将第二、三行减去第一行，，得到：

$$ \mathbf{A} \xrightarrow{r_2 - r_1, r_3 - r_1} \mathbf{B} = \begin{bmatrix}1 & 1 & 1 \ 0 & 2 & 1 \ 0 & 4 & 6\end{bmatrix} $$

此时第一列已经满足三角矩阵的要求了，对第二列进行操作。

将第三行减去两倍的第二行，得到：

$$ \mathbf{B} \xrightarrow{r_3 - 2r_2} \mathbf{C} = \begin{bmatrix}1 & 1 & 1 \ 0 & 2 & 1 \ 0 & 0 & 4\end{bmatrix} $$

根据性质 2，$|\mathbf{C}| = 1 \times 2 \times 4 = 8$，则 $|\mathbf{A}| = |\mathbf{B}| = |\mathbf{C}| = 8$，我们就完成了对矩阵 $\mathbf{A}$ 的行列式求解。

从特殊到一般，我们可以这样描述我们的算法流程：

• 枚举 $i=1,2,\ldots,n$，选取 $a_{i,i}$，对于第 $j$ 行($j=i+1,i+2,\ldots,n)$，整行减去第 $i$ 行的 $\dfrac{a_{j,i}}{a_{i,i}}$ 倍。
• 重复此流程，直到 $i=n$.
• 计算 $\prod \limits {i=1}^n a{i,i}$，即为所求的行列式。

可以发现，第一步完成后，第 $i+1$ 行到第 $n$ 行的第 $i$ 列都为零。反复消去，就能得到一个上三角矩阵。

但这里需要注意一个 corner case：$a_{i,i} = 0$ 怎么办。
在第一步中，如果 $a_{i,i}=0$，我们就无法用第 $i$ 行消去其余行的第 $i$ 列。

一个合理的做法是：遍历第 $i+1$ 行到第 $n$ 行，找到 $a_{j,i} \neq 0$ 的行，将其交换到第 $i$ 行，再进行消元。

还是举个例子：

$$ \mathbf{A} = \begin{bmatrix}1 & 1 & 1\ 1 & 1 & 0 \ 0 & 2 & 3\end{bmatrix} $$

第一步消元之后得到：

$$ \begin{bmatrix}1 & 1 & 1\ 0 & 0 & -1 \ 0 & 2 & 3\end{bmatrix} $$

‍

此时 $a_{2,2}=0$，我们找到第三行 $a_{3,2} \neq 0$，交换到第二行，继续执行消元：

$$ \begin{bmatrix}1 & 1 & 1\ 0 & 0 & -1 \ 0 & 2 & 3\end{bmatrix} \xrightarrow{r_2\leftrightarrow r_3}\begin{bmatrix}1 & 1 & 1\ 0 & 2 & 3 \ 0 & 0 & -1 \end{bmatrix} $$

此时交换后恰好无须消元，算法结束。

需要注意的是，这样的交换过后，根据性质 3，行列式变号。因此在算法过程中需要在交换时额外处理一下。

进一步的 corner case：假如第 $i$ 行到第 $n$ 行的第 $j$ 列全都为零呢？

不失一般性，我们举一个例子来考虑，从第三行开始无法消元的情况：

$$ \mathbf{A}= \begin{bmatrix} a_{1,1} & a_{1,2} & a_{1,3} & \cdots & a_{1,n} \ 0 & a_{2,2} & a_{1,3} & \cdots & a_{2,n} \ 0 & 0 & 0 & \cdots & a_{3,n} \ 0 & 0 & 0 & \cdots & a_{4,n} \ \vdots & \vdots & \vdots & &\vdots \ 0 & 0 & 0 & \cdots & a_{n,} \end{bmatrix} $$

此时，我们对第一列做拉普拉斯展开，得到：

$$ |\mathbf{A}|= a_{1,1} \times \begin{vmatrix} a_{2,2} & a_{1,3} & \cdots & a_{2,n} \ 0 & 0 & \cdots & a_{3,n} \ 0 & 0 & \cdots & a_{4,n} \ \vdots & \vdots & &\vdots \ 0 & 0 & \cdots & a_{n,} \end{vmatrix} $$

再对第二列进行拉普拉斯展开，即有：

$$ |\mathbf{A}|= a_{1,1} \times a_{2,2} \times \begin{vmatrix} 0 & \cdots & a_{3,n} \ 0 & \cdots & a_{4,n} \ \vdots & &\vdots \ 0 & \cdots & a_{n,} \end{vmatrix} $$

根据性质 4，展开后的矩阵行列式为零，则 $|\mathbf{A}| = 0$.

更一般的，若从第 $i$ 行开始无法消元，则对 $\mathbf{A}$ 进行 $i-1$ 次展开后，余子式第一列必定全为零，则 $|\mathbf{A}| = 0$.

因此，如果我们在算法执行过程中遇到无法消元的困境，直接返回结果零即可。

Implementation

「Talk is cheap, show me the code.」

在实现中可以有一个小细节：我们在利用 $a_{i,i}$ 消元时，可以找到 $|a_{j,i}|$ 最大的值所在行，将其交换到第 $i$ 行。

据说这样可以增加精度，可能是算法竞赛中的一些都市传说吧……

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

#include <algorithm>
#include <math.h>
#include <stdio.h>
const int maxn = 1e3 + 10;
const double EPS = 1e-6;
double A[maxn][maxn], det = 1.0;
int n;
int main()
{
    scanf("%d", &n);
    for (int i = 1; i <= n; ++i)
        for (int j = 1; j <= n; ++j)
            scanf("%lf", A[i] + j);
    for (int r = 1; r < n; ++r)
    {
        // find the row with max a_{j,i} and swap it to row i
        int maxx = r;
        for (int i = r + 1; i <= n; ++i)
            if (fabs(A[i][r]) > fabs(A[maxx][r]))
                maxx = i;
        if (maxx != r)
        {
            det *= -1;  // swap two rows change the sign of determinant
            for (int i = 1; i <= n; ++i)
                std::swap(A[r][i], A[maxx][i]);
        }
        // exit if all a_{j,i} equals to zero
        if (fabs(A[r][r]) < EPS) break;
        // eliminate [r+1, ... , n] rows
        for (int i = r + 1; i <= n; ++i)
        {
            if (i == r) continue;
            double ratio = A[i][r] / A[r][r];
            for (int j = 1; j <= n; ++j)
                A[i][j] -= A[r][j] * ratio;
        }
    }
    for (int i = 1; i <= n; ++i) det *= A[i][i];
    printf("%.2f", det);
    return 0;
}

事实上，这个算法与高斯消元(Gauss Elimination)几乎一致。

‍