基于LU分解的行列式计算算法
题目描述
我们考虑如何利用矩阵的LU分解高效地计算其行列式(determinant)。给定一个 \(n \times n\) 的可逆方阵 \(A\),我们需要计算 \(\det(A)\)。直接使用行列式的定义(莱布尼茨公式)或拉普拉斯展开计算复杂度为 \(O(n!)\) 或 \(O(n^3)\),对于大型矩阵不可行。LU分解将矩阵 \(A\) 分解为一个下三角矩阵 \(L\) 和一个上三角矩阵 \(U\) 的乘积(有时包含行交换矩阵 \(P\)),即 \(PA = LU\)。利用三角矩阵行列式的性质,可以快速计算出 \(\det(A)\)。本题目要求阐述基于LU分解的行列式计算算法的原理、步骤、复杂度及数值稳定性。
解题过程
第一步:理解行列式的基本性质
在介绍算法前,需要回顾几个关键性质:
- 三角矩阵的行列式:上三角矩阵 \(U\) 或下三角矩阵 \(L\) 的行列式等于其主对角线元素的乘积。
\[ \det(U) = u_{11} u_{22} \cdots u_{nn}, \quad \det(L) = l_{11} l_{22} \cdots l_{nn}. \]
- 乘积的行列式:对于任意两个 \(n \times n\) 矩阵 \(B, C\),有
\[ \det(BC) = \det(B) \det(C). \]
- 置换矩阵的行列式:如果 \(P\) 是一个置换矩阵(代表行交换),则 \(\det(P) = +1\) 或 \(-1\),具体取决于交换次数的奇偶性(奇数次为 -1,偶数次为 +1)。
第二步:回顾LU分解(带部分选主元)
对于一般的方阵 \(A\),为了数值稳定性,通常采用部分选主元的LU分解:
\[PA = LU, \]
其中:
- \(P\) 是一个置换矩阵,记录行交换操作。
- \(L\) 是单位下三角矩阵(主对角线元素均为1)。
- \(U\) 是上三角矩阵。
如果不采用选主元(例如矩阵是严格对角占优或对称正定),则分解为 \(A = LU\),此时 \(P = I\)。
第三步:推导行列式计算公式
由 \(PA = LU\),两边取行列式:
\[\det(P) \det(A) = \det(L) \det(U). \]
因为 \(L\) 是单位下三角矩阵,其行列式 \(\det(L) = 1\)。所以:
\[\det(A) = \frac{\det(U)}{\det(P)}. \]
又因为 \(\det(P) = \pm 1\),并且 \(\det(U) = \prod_{i=1}^{n} u_{ii}\),所以最终公式为:
\[\det(A) = \left( \prod_{i=1}^{n} u_{ii} \right) \cdot \det(P), \]
其中 \(\det(P) = (-1)^k\),\(k\) 是在LU分解过程中执行的行交换次数。
特别情况:如果未使用行交换(\(P = I\)),则 \(\det(A) = \prod_{i=1}^{n} u_{ii}\)。
第四步:算法步骤
- 输入:\(n \times n\) 矩阵 \(A\)。
- 执行带部分选主元的LU分解:
- 初始化置换向量 \(p = [1, 2, \dots, n]^T\) 或直接记录交换次数 \(k = 0\)。
- 对于每一列 \(j = 1, \dots, n-1\):
a. 在列 \(j\) 的第 \(j\) 行到第 \(n\) 行中,找到绝对值最大的元素所在行 \(i_{\text{max}}\)。
b. 如果 \(i_{\text{max}} \neq j\),则交换 \(A\) 的第 \(j\) 行和第 \(i_{\text{max}}\) 行,并更新置换记录(\(k\) 增加1)。
c. 对于 \(i = j+1, \dots, n\):- 计算乘子 \(\ell_{ij} = a_{ij} / a_{jj}\) 并存储在 \(A(i,j)\)(作为 \(L\) 的元素)。
- 更新 \(A(i, j+1:n) = A(i, j+1:n) - \ell_{ij} \cdot A(j, j+1:n)\)。
- 最终矩阵 \(A\) 的上三角部分(包括主对角线)即为 \(U\),下三角部分(主对角线以下)存储了 \(L\) 的乘子(主对角线为1,不显式存储)。
- 计算行列式:
- 提取 \(U\) 的主对角线元素 \(u_{11}, u_{22}, \dots, u_{nn}\)。
- 计算乘积 \(d = u_{11} u_{22} \cdots u_{nn}\)。
- 根据行交换次数 \(k\),令 \(\det(A) = (-1)^k \cdot d\)。
- 输出:\(\det(A)\)。
第五步:算法复杂度与数值稳定性
- 时间复杂度:LU分解本身需要约 \(\frac{2}{3} n^3\) 次浮点运算(加法和乘法)。计算行列式只需要 \(O(n)\) 次乘法(计算对角线乘积),因此总复杂度为 \(O(n^3)\),远低于直接计算的 \(O(n!)\)。
- 数值稳定性:部分选主元(选择列中最大元素作为主元)可以避免除零和小主元导致的数值不稳定,是实践中标准做法。即使 \(A\) 是奇异或接近奇异的,算法仍能执行(但行列式可能为零或非常小,需要注意浮点误差)。
第六步:示例演示
考虑矩阵
\[A = \begin{bmatrix} 2 & 1 & 1 \\ 4 & 1 & 0 \\ -2 & 2 & 1 \end{bmatrix}. \]
步骤1:执行LU分解(带选主元)
- 第一列:主元为4(第2行),交换第1行和第2行,\(k=1\)。
交换后:
\[ A' = \begin{bmatrix} 4 & 1 & 0 \\ 2 & 1 & 1 \\ -2 & 2 & 1 \end{bmatrix}. \]
计算乘子:\(\ell_{21} = 2/4 = 0.5, \ell_{31} = -2/4 = -0.5\)。
更新第2、3行:
\[ \begin{aligned} \text{第2行} &:= [2,1,1] - 0.5 \cdot [4,1,0] = [0,0.5,1], \\ \text{第3行} &:= [-2,2,1] - (-0.5) \cdot [4,1,0] = [0,2.5,1]. \end{aligned} \]
此时矩阵为:
\[ \begin{bmatrix} 4 & 1 & 0 \\ 0.5 & 0.5 & 1 \\ -0.5 & 2.5 & 1 \end{bmatrix}. \]
- 第二列:主元为2.5(第3行),交换第2行和第3行,\(k=2\)。
交换后:
\[ \begin{bmatrix} 4 & 1 & 0 \\ -0.5 & 2.5 & 1 \\ 0.5 & 0.5 & 1 \end{bmatrix}. \]
计算乘子:\(\ell_{32} = 0.5/2.5 = 0.2\)。
更新第3行:
\[ \text{第3行} := [0.5,0.5,1] - 0.2 \cdot [ -0.5,2.5,1] = [0.5,0.5,1] - [-0.1,0.5,0.2] = [0.6,0,0.8]. \]
最终得到:
\[ U = \begin{bmatrix} 4 & 1 & 0 \\ 0 & 2.5 & 1 \\ 0 & 0 & 0.8 \end{bmatrix}, \quad L = \begin{bmatrix} 1 & 0 & 0 \\ -0.5 & 1 & 0 \\ 0.5 & 0.2 & 1 \end{bmatrix}. \]
步骤2:计算行列式
- \(U\) 的对角线乘积:\(4 \times 2.5 \times 0.8 = 8\)。
- 行交换次数 \(k = 2\),所以 \((-1)^k = 1\)。
- 因此 \(\det(A) = 8\)。
可以验证直接计算 \(\det(A) = 2(1-0) -1(4-0) + 1(8-(-2)) = 2 -4 +10 = 8\),结果一致。
总结
基于LU分解的行列式计算算法将复杂的行列式计算转化为高效的矩阵分解问题,利用三角矩阵行列式的简便性,在 \(O(n^3)\) 时间内得到结果。通过部分选主元保证数值稳定性,适用于大多数可逆矩阵。该算法是数值线性代数中计算行列式的标准方法,广泛应用于科学计算和工程问题中。