行列式(Determinant) 是矩阵的重要属性。
在手动计算行列式时,我们常常使用两种方法:
前者的复杂度是 O(n!) 级别的,在计算约 12 阶的矩阵时就会需要超过 1s 的时间,而计算 1000 \times 1000 的矩阵需要进行约:
402387260077093773543702433923003985719374864210714632543799910429938512398629020592044208486969404800479988610197196058631666872994808558901323829669944590997424504087073759918823627727188732519779505950995276120874975462497043601418278094646496291056393887437886487337119181045825783647849977012476632889835955735432513185323958463075557409114262417474349347553428646576611667797396668820291207379143853719588249808126867838374559731746136085379534524221586593201928090878297308431392844403281231558611036976801357304216168747609675871348312025478589320767169132448426236131412508780208000261683151027341827977704784635868170164365024153691398281264810213092761244896359928705114964975419909342221566832572080821333186116811553615836546984046708975602900950537616475847728421889679646244945160765353408198901385442487984959953319101723355556602139450399736280750137837615307127761926849034352625200015888535147331611702103968175921510907788019393178114194545257223865541461062892187960223838971476088506276862967146674697562911234082439208160153780889893964518263243671616762179168909779911903754031274622289988005195444414282012187361745992642956581746628302955570299024324153181617210465832036786906117260158783520751516284225540265170483304226143974286933061690897968482590125458327168226458066526769958652682272807075781391858178889652208164348344825993266043367660176999612831860788386150279465955131156552036093988180612138558600301435694527224206344631797460594682573103790084024432438465657245014402821885252470935190620929023136493273497565513958720559654228749774011413346962715422845862377387538230483865688976461927383814900140767310446640259899490222221765904339901886018566526485061799702356193897017860040811889729918311021171229845901641921068884387121855646124960798722908519296819372388642614839657382291123125024186649353143970137428531926649875337218940694281434118520158014123344828015051399694290153483077644569099073152433278288269864602789864321139083506217095002597389863554277196742822248757586765752344220207573630569498825087968928162753848863396909959826280956121450994871701244516461260379029309120889086942028510640182154399457156805941872748998094254742173582401063677404595741785160829230135358081840096996372524230560855903700624271243416909004153690105933983835777939410970027753472000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
次运算。
这样计算行列式的效率显然是极低的。而通过化三角矩阵,我们可以用 O(n^3) 的复杂度完成行列式的求解。对于同样的矩阵,我们只需要进行 1 \times 10^9 的运算。这对于中小规模的矩阵已经足够快速了。
对于矩阵 \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} $$
如果你了解高斯消元相关的内容,那再好不过了。
通过性质 1,我们可以对矩阵进行变换,将其化为三角矩阵,从而通过性质 2 的方法求解行列式。
先从一个具体的例子入手。
将第二、三行减去第一行,,得到:
此时第一列已经满足三角矩阵的要求了,对第二列进行操作。
将第三行减去两倍的第二行,得到:
根据性质 2,|\mathbf{C}| = 1 \times 2 \times 4 = 8,则 |\mathbf{A}| = |\mathbf{B}| = |\mathbf{C}| = 8,我们就完成了对矩阵 \mathbf{A} 的行列式求解。
从特殊到一般,我们可以这样描述我们的算法流程:
可以发现,第一步完成后,第 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.
因此,如果我们在算法执行过程中遇到无法消元的困境,直接返回结果零即可。
「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)几乎一致。