化三角矩阵计算行列式的算法实现

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{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} $$

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

Theory

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

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

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

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

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

根据性质 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 行，再进行消元。

还是举个例子：

第一步消元之后得到：

‍

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

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

需要注意的是，这样的交换过后，根据性质 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 行。 据说这样可以增加精度，可能是算法竞赛中的一些都市传说吧……

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

‍