线性系统

5.31 线性系统

高斯消元法

主元：pivot 归一化，不断向下消除

主元为0时，与当前列最大值的行交换（避免误差）

缺点：化成阶梯矩阵后还需要回代求解（繁琐）

高斯-约旦消元法（Gauss-Jordan Elimination）

没b用，就是再反向进行一遍高斯消元法，化成行最简矩阵

实现G-J

我的实现：

   def gauss_jordan_elimination(self):
        """返回Gauss-Jordan elimination的结果"""
        T = self
        # TODO 增广矩阵判断，rank实现以后需要重写
        assert T.col_num() == T.row_num() + 1, \
            "Error in Gauss-Jordan elimination! Not an augmented-matrix."
        for i in range(T.row_num()):
            if T[i, i] == 0:
                for j in range(i + 1, T.row_num()):
                    if T[j, i] != 0:
                        tmp = T[j]
                        T[j] = T[i]
                        T[i] = tmp
                        break
            T[i] /= T[i, i]
        for i in range(T.row_num - 2,-1,-1):
            if T[i, i+1] != 0:
                T[i] -= T[i+1] * T[i, i+1]
        return T

出现问题：

Traceback (most recent call last):
  File "E:/LinearAlgebra/main_matrix.py", line 55, in <module>
    print("{}.gauss_jordan_elimination() is {}".format_map(T64, T64.gauss_jordan_elimination()))
  File "E:\LinearAlgebra\playLA\Matrix.py", line 139, in gauss_jordan_elimination
    T[i] /= T[i, i]
  File "E:\LinearAlgebra\playLA\Matrix.py", line 25, in __getitem__
    r, c = pos
TypeError: cannot unpack non-iterable int object

为了不打破Vector不可变的性质：

添加underlying_list: 返回底层的list，用来操作

vec：返回vec的引用

vec[:] 或 vec.copy()：返回vec的副本

行最简矩阵 & 解的结构

阶梯矩阵 --Gauss-Jordan Elimination–> 行最简矩阵（RREF：Reduced Row Echelon Form）

$\{\begin{array}{l}无解\\ 唯一解\\ 无数解\end{array}$

n个未知数，n个约束（方程），才可能有唯一解

行最简矩阵：系数矩阵的非零行不可能大于未知数个数 (可以从秩或空间的角度考虑)

系数矩阵非方阵的Gauss-Jordan Elimination

主元列：pivot colum

自由列：freedom colum

• 我的实现：

    def _get_max_pivot(self, r_num, c_num, n):
        """获第index列最大元素所在的行号"""
        max, ret = self.Ab[r_num][c_num], r_num

        for i in range(r_num + 1, n):
            if self.Ab[i][c_num] > max:
                max, ret = self.Ab[i][c_num], i

        return ret

    def _forward(self):
        """前向化简：Gauss Elimination"""
        n = self._row_num
        # now_index = 0   # 当前列号
        for i in range(n):
            now_index = i   # TODO 这里可以优化：维持一个全局的pivot的列坐标
            while now_index < self._col_num and self.Ab[i][now_index] == 0:
                max_row = self._get_max_pivot(r_num=i, c_num=now_index, n=n)
                if self.Ab[max_row][now_index] == 0:
                    now_index += 1
                    continue
                self.Ab[i], self.Ab[max_row] = self.Ab[max_row], self.Ab[i]
                break
            if self.Ab[i][now_index] != 0:
                self.Ab[i] = self.Ab[i] / self.Ab[i][now_index] # TODO self.Ab[i][i] == 0
                for j in range(i+1, n):
                    self.Ab[j] = self.Ab[j] - self.Ab[i] * self.Ab[j][now_index]
                self._backward(i, now_index)

    def _backward(self, r_num, c_num):
        """将第index列，后上化简"""
        for i in range(r_num - 1, -1, -1):
            # self.Ab[r_num][c_num] 为 pivot
            self.Ab[i] = self.Ab[i] - self.Ab[r_num] * self.Ab[i][c_num]

    def gauss_jordan_elimination(self):
        """返回Gauss-Jordan Elimination的结果"""
        self._forward()

• 测试：

    # 6-8 实现更一般的Gauss-Jordan Elimination
    A7 = Matrix([
        [1, -1, 2, 0, 3],
        [-1, 1, 0, 2, -5],
        [1, -1, 4, 2, 4],
        [-2, 2, -5, -1, -3]
    ])
    b7 = Vector([1, 5, 13, -1])
    ls7 = LinearSystem(A7, b7)
    ls7.gauss_jordan_elimination()
    ls7.fancy_print()
    print()


    A8 = Matrix([
        [2, 2],
        [2, 1],
        [1, 2]
    ])
    b8 = Vector([3, 1.5, 7])
    ls8 = LinearSystem(A8, b8)
    ls8.gauss_jordan_elimination()
    ls8.fancy_print()
    print()
    
        # TODO Error: 化简过程中某一行全为0
    A9 = Matrix([
        [1, 2, 3],
        [2, 4, 6],
        [3, 5, 7],
        [5, 8, 10]
    ])
    b9 = Vector([1, 2, 3, 4])
    ls9 = LinearSystem(A9, b9)
    ls9.gauss_jordan_elimination()
    ls9.fancy_print()
    print()