2.1 高斯消元

数学原理

  雷米兹算法的基础是N元一次方程的求解。此外，离散傅里叶变换也用到了N元一次方程的求解。当然，实际应用中求解多元一次方程组的例子非常多，我就不一一列举。求解的方法主要有两种，一种是高斯消元，一种是递归法。
  递归法不推荐，已经被我删除。而高斯消元不仅仅是用来解方程，还可以用来求行列式、逆矩阵。
  高斯消元是利用矩阵表示N元一次方程组。然后将矩阵变成上三角矩阵或下三角矩阵。
  所谓上三角矩阵就是左下角全部为0的矩阵，下三角反之。
  如下，就是转换过程：
$$\left[\begin{array}{ccccc}1& 1& 1& 1& 10\\ 2& 1& 1& 1& 11\\ -1& 2& 1& 1& 10\\ -3& 2& 1& 2& 14\end{array}\right]$$
  将第一列变成0：
$$\left[\begin{array}{ccccc}1& 1& 1& 1& 10\\ 0& 1& 1& 1& 9\\ 0& -3& -2& -2& -20\\ 0& -5& -4& -5& -44\end{array}\right]$$
  将第二列变成0：
$$\left[\begin{array}{ccccc}1& 1& 1& 1& 10\\ 0& 1& 1& 1& 9\\ 0& 0& -1& -1& -7\\ 0& 0& -1& 0& -1\end{array}\right]$$
  将第三列变成0：
$$\left[\begin{array}{ccccc}1& 1& 1& 1& 10\\ 0& 1& 1& 1& 9\\ 0& 0& -1& -1& -7\\ 0& 0& 0& 1& 6\end{array}\right]$$
  这个转换过程，就是先用第一行和后续的所有行进行计算，让后续所有行以0开始。
  然后将右下角视为一个子矩阵，继续这个运算。这个算法可以用递归实现，也可以用循环实现。但是最好实用循环实现，也就是高斯消元方法。
  然后使用代入法就可以解开，但是这时候千万不要去再转换为下三角矩阵，并且进一步转为单位矩阵。因为那样会有很多不必要的运算。

转换为三角矩阵

  这个核心的行变换方法其实就是行倍加。假设第一行的系数是$a$,第二行的系数是$b$,那么就是把第一行乘以$-\frac{b}{a}$加到第二行上。
  行例如第一行是$x+2y+3z=6$，第二行是$3x-3y+9z=9$，那么第一行就乘以$-\frac{3}{1}=-3$，然后加到第二行。
  代码如下：

public Matrix<T> toUpperTriangular() {
    // 大循环，先用第一行合并以后所有行
    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], ts[j][i], ts[i][i], i, ts[i].length);
        }
    }
    return createMatrix(ts);
}

/**
 * 行变换
 *
 * @param line1
 * @param line2       要变的行
 * @param coefficent1 第一行乘以的系数
 * @param coefficent2 第二行乘以的系数
 * @return
 */
public void rowOperator(T[] line1, T[] line2, T coefficent1, T coefficent2, int startColumn, int endColumn) {
    if (line1 == null || line2 == null) {
        throw new NullPointerException("合并的方程不能为空");
    }
    if (line1.length != line2.length) {
        throw new NullPointerException("两个方程参数不一致");
    }
    for (int i = startColumn; i < endColumn; i++) {
        final T ai = line1[i];
        final T bi = line2[i];
= subtract(multiply(ai, coefficent1), multiply(bi, coefficent2));
    }
}

  而python代码则如此：

def to_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], i, len(lines[i]))
    return Matrix(lines)

行变换代码

@staticmethod
def row_operate(line1, line2, coefficient1, coefficient2, start_column, end_column):
    for i in range(start_column, end_column):
        line2[i] = line1[i] * coefficient1 - line2[i] * coefficient2

代入法

  而获取了上三角矩阵之后。
  有两种方法可以解开方程，一种办法是代入法，另一种方法是从下到上使用行倍加变换。代入法性能高一点，性能高的原因是避免了不必要的乘以0的运算。
  就是套进去解开了。解开的方法代码如下：
  其具体的思路是用最后一行代进去得到一个数字，比如如下的上三角矩阵:
$$\left[\begin{array}{ccccc}1& 1& 1& 1& 10\\ 0& 1& 1& 1& 9\\ 0& 0& -1& -1& -7\\ 0& 0& 0& 1& 6\end{array}\right]$$
  从最后一行开始，得到$x_4=6$，然后$x_4=6$代入上面的每一行，矩阵变成了这样
$$\left[\begin{array}{ccccc}1& 1& 1& 0& 4\\ 0& 1& 1& 0& 3\\ 0& 0& 1& 0& 1\\ 0& 0& 0& 1& 6\end{array}\right]$$
  然后把$x_3$求出来，代入上面的每一行，最终，整个矩阵变成
$$\left[\begin{array}{ccccc}1& 0& 0& 0& 1\\ 0& 1& 0& 0& 2\\ 0& 0& 1& 0& 1\\ 0& 0& 0& 1& 6\end{array}\right]$$
  所以这个方程的解就是：
$$\begin{gathered} x_1=1 \\ {x}_2=2 \\ {x}_3=1 \\ {x}_4=6 \end{gathered}$$

最终代码

  代入法计算性能更高，如果再转为下三角会有很多非必要的预算，再转为单位矩阵，又是一堆没必要的运算。高斯消元总体代码如下：

public T[] gaussianReduction() {
    final Matrix<T> tMatrix = this.toUpperTriangular();
    // 进行代入法变换。
    T[] solution = newArray(this.array.length);
    final T[][] ts = tMatrix.copyArray();
    final int n = ts.length;
    System.out.println();
    System.out.println(tMatrix.toBeautifulString());
    // 从最后一行出来-
    for (int i = n - 1; i >= 0; i--) {
        // 首先计算出自己
        solution[i] = divide(ts[i][n], ts[i][i]);
        // 遍历之上的所有行
        for (int j = 0; j < i; j++) {
            ts[j][n] = subtract(ts[j][n], multiply(ts[j][i], solution[i]));
        }
    }
    return solution;
}

  python代码如下：

def gaussian_reduction(self):
    matrix = self.to_upper_triangle()
    solution = [0 for _ in matrix.__lines]
    n = len(matrix.__lines)
    for i in range(n - 1, -1, -1):
        solution[i] = matrix.__lines[i][n] / matrix.__lines[i][i]
        for j in range(0, i):
            matrix.__lines[j][n] = matrix.__lines[j][n] - matrix.__lines[j][i] * solution[i]
    return solution