4.6 高斯消元法与Doolittle算法

理论基础

  行列式有这么一个性质：
$$|AB|=|A||B|$$
  这个性质可以快速的计算矩阵的行列式。具体怎么算呢?前面我说过初等行变换相当于左乘一个初等矩阵，假设对矩阵进行一系列初等行变换变成了上三角矩阵$U$，就相当于左乘了$L_n,L_{n-1},\cdots ,L_2,L_1$。那么矩阵的行列式就可以这么计算：
$$\begin{gathered} L_n{L}_{n-1}\cdots L_2{L}_{1}A=U \\ |L_n||L_{n-1}|\cdots |L_2||L_1||A|=|U| \end{gathered}$$
  接下来，再研究初等矩阵的行列式，经过前人的研究发现，三种初等行变换有以下性质：

1. 行交换的初等矩阵行列式为-1;
2. 行倍加的初等矩阵行列式为1；
3. 行倍数的初等矩阵行列式为倍数k。

  再看看上三角矩阵，因为左下角全部是0，容易观察出，上三角矩阵的行列式等于对角线元素的乘积。所以算法就简单了：

1. 进行初等行变换，如果交换了两行，则记录一个倍数-1，如果倍乘了一行，则记录倍数k；
2. 计算获得的上三角矩阵的对角线元素的乘积；
3. 将乘积除于第一步记录的倍数的乘积，结果就是矩阵的行列式。

  举个例子，以下矩阵：
$$\left[\begin{array}{ccc}1& 1& 1\\ 2& 1& 1\\ 2& 2& 1\end{array}\right]$$
  变成上三角形式：
$$\left[\begin{array}{ccc}1& 1& 1\\ 0& 1& 1\\ 0& 0& -1\end{array}\right]$$
  把对角线乘起来，行列式就是$1\times 1\times -1=-1$。

Doolittle法

  这种方法，如果不使用行倍数和行交换也叫Doolittle法。在这里，我解释清楚，因为这是我带头的一个错误。错误根源在于一个外文网站，它的错误在于混淆了Crout分解与Doolittle分解，先看原文：

The Kraut algorithm finds decomposition of matrix  $A$  as  $A=LU$  where  $L$  is lower triangular and  $U$  is upper triangular matrix. Without loss of generality, we can assume that all the diagonal elements of  $L$  are equal to 1. Once we know these matrices, it is easy to calculate the determinant of  $A$ : it is equal to the product of all the elements on the main diagonal of the matrix  $U$ .

  什么意思呢？LU分解为L是单位下三角矩阵，行列式就等于U矩阵对角线元素的乘积，这不是Doolittle分解吗？另外，怎么名字是Kraut呢？而且主流文章，用的是Crout，而不是拼写错误的Kraut。我之前没认真考证这篇文章，还写了一篇博文以讹传讹，现在我已经将那篇博文删除了。
  此外，LU分解求行列式，本质上就是只使用行倍加的高斯消元，为什么呢？前面我写过这样的公式：
$$L_n{L}_{n-1}\cdots L_2{L}_{1}A=U$$
  再进行改造，两边乘以初等矩阵的逆矩阵，消掉左边：
$$\begin{gathered} (L_1^{-1}{L}_2^{-1}\cdots L_{n-1}^{-1}{L}_n^{-1})L_n{L}_{n-1}\cdots L_2{L}_{1}A=L_1^{-1}{L}_2^{-1}\cdots L_n^{-1}U \\ \Rightarrow A=L_1^{-1}{L}_2^{-1}\cdots L_n^{-1}U \\ L=L_1^{-1}{L}_2^{-1}\cdots L_n^{-1} \\ A=LU \end{gathered}$$
  把所有初等矩阵的逆矩阵乘起来，记作L，这就是LU分解。我们知道以下事实：

1. 初等矩阵一定有逆矩阵，它的逆矩阵还是初等矩阵；
2. 行倍加初等矩阵的逆矩阵还是行倍加初等矩阵
3. 行倍加的初等矩阵对角线元素全部为1；
4. 多个行倍加的初等矩阵相乘对角线元素依旧全部为1。

  基于以上四点，所以只用行倍加的高斯消元，只要记录行变换过程，就成了Doolittle分解。因为计算行列式只需要用到U矩阵，所以只要严格只用行倍加，连记录行变换过程都不需要了。

代码实现

  所以我们要写一个新的转换上三角矩阵的方法。这个新的方法，不会把要变换的行乘以k倍。

public Matrix<T> toBaseUpperTriangular() {
    // 大循环，先用第一行合并以后所有行
    final T[][] ts = copyArray();
    for (int i = 0; i < ts.length - 1; i++) {
        // 循环遍历其后的所有行,其后所有行进行改变
        for (int j = i + 1; j < ts.length; j++) {
            rowOperator(ts[i], ts[j], divide(ts[j][i],ts[i][i]), oneValue(), i, ts[i].length);
        }
    }
    return createMatrix(ts);
}

  然后就是行列式了

public T determinant() {
    final Matrix<T> matrix = this.toBaseUpperTriangular();
    T result = oneValue();
    for (int i = 0; i < matrix.array.length; i++) {
        result = multiply(result, matrix.array[i][i]);
    }
    return result;
}

  python代码也很简单：

def to_base_upper_triangle(self):
    # 将矩阵变成上三角矩阵
    lines = copy.deepcopy(self.__lines)
    for i in range(len(lines) - 1):
        for j in range(i + 1, len(lines)):
            Matrix.row_operate(lines[i], lines[j], lines[j][i]/lines[i][i], 1, i, len(lines[i]))
    return Matrix(lines)

  行列式代码

def determinant(self):
    matrix = self.to_base_upper_triangle()
    result = 1
    for i in range(0, len(matrix.__lines)):
        result *= matrix.__lines[i][i]
    return result

  python源码地址：https://e.coding.net/buildt/python-study/math-algorithm.git

结语

  至此，行列式计算的六种主流算法就介绍完了，接下来，就该讨论行列式与线性方程组的关系了，也就是著名的克拉默法则。