5.2 LUP分解

解方程步骤

  LUP分解是一种解方程的办法，这个办法和LU分解不同的地方的地方，在于如果对角线上存在0，或者舒尔补的对角线上存在0的时候，会出现除于0的操作，这个时候是没法进行LU分解的，所以引入了置换矩阵。
  但是LUP分解的出现不仅仅是为了除以0，而是为了尽量减少误差。有一类矩阵叫做病态矩阵(深入了解的话，可以看我的博客矩阵的条件数)，这种矩阵会将LU分解的计算误差放大成百上千倍。而浮点数的运算，误差无法避免。
  什么是置换矩阵呢？就是单位矩阵进行行变换后产生的矩阵就叫做置换矩阵。置换矩阵与普通矩阵的乘法是对矩阵进行行交换或者列交换。如果置换矩阵在左，就是行变交，如果置换矩阵在右则是列交换。
  LUP分解是找出三个矩阵，满足以下等式：
$$PA=LU$$
  L是单位下三角矩阵。
  U是上三角矩阵。
  P是置换矩阵。所以加上方程组和值的向量就变成了以下形式：
$$\begin{gathered} PAx=Pb \\ LUx=Pb \end{gathered}$$
  定义向量y = Ux，有：
$$Ly=Pb$$
  计算$Pb$时没必要用矩阵乘法去做。直接进行行交换就行。
  因为L是单位下三角矩阵，所以用代入法很快能得到y这个向量。
  那么$y=Ux$，也就是$Ux=y$，因为U是上三角，代入法很快能得到x这个向量的解。
  这整个过程除了和和前面的LU分解差不多，难就难在置换。

求LUP三个矩阵的步骤

  LUP分解比LU分解稍微复杂一点。因为LUP分解不仅仅是为了处理0，而是为了解决病态矩阵问题，所以置换的条件每次把绝对值最大的值放在第一行。假设第一列绝对值最大的行为k，第k行和第一行进行交换。这个交换，是左边乘以一个置换矩阵，这个矩阵我们定义为Q。再按以下方式分块：
$$QA=\left(\begin{array}{cc}{a}_{k1}& w_T\\ v& A^{\prime }\end{array}\right)$$
最终的置换矩阵是P，而P和我们仅对第一列进行置换的矩阵Q有以下关系：
$$P=\left(\begin{array}{cc}1& 0\\ 0& P^{\prime }\end{array}\right)Q$$
  虽然有上面的公式，但是根据置换矩阵的性质，没必要在内存中存储这么大的一个完整矩阵，也没必要用稀疏数组去存储，只需要用置换的理论，存一个一维数组就够了。这个一维数组，下标表示原始的列，存储的值表示交换后的列。
  LUP分解中，分块方式变了，是先交换行再进行分块。LUP分解需要注意的地方还有舒尔补的那一部分也需要置换，也就是：
$$P^{\prime }(A^{\prime }-v{w}^T/a_{k1})=L^{\prime }{U}^{\prime }$$
  这样递归下去。如果用递归做代码倒是好写，但是和LU分解一样，我们大可不必递归，用循环就行。
  完整的证明过程如下：
$$\begin{gathered} PA=\left(\begin{array}{cc}1& 0\\ 0& P^{\prime }\end{array}\right)QA \\ =\left(\begin{array}{cc}1& 0\\ 0& P^{\prime }\end{array}\right)\left(\begin{array}{cc}1& 0\\ v/a_{k1}& I_{n-1}\end{array}\right)\left(\begin{array}{cc}{a}_{k1}& w^T\\ 0& A^{\prime }-v{w}^T/a_{k1}\end{array}\right) \\ =\left(\begin{array}{cc}1& 0\\ P^{\prime }v/a_{k1}& P^{\prime }\end{array}\right)\left(\begin{array}{cc}{a}_{k1}& w^T\\ 0& A^{\prime }-v{w}^T/a_{k1}\end{array}\right) \\ =\left(\begin{array}{cc}1& 0\\ P^{\prime }v/a_{k1}& I_{n-1}\end{array}\right)\left(\begin{array}{cc}{a}_{k1}& w^T\\ 0& P^{\prime }(A^{\prime }-v{w}^T/a_{k1})\end{array}\right) \\ =\left(\begin{array}{cc}1& 0\\ P^{\prime }v/a_{k1}& I_{n-1}\end{array}\right)\left(\begin{array}{cc}{a}_{k1}& w^T\\ 0& L^{\prime }{U}^{\prime }\end{array}\right) \\ =\left(\begin{array}{cc}1& 0\\ P^{\prime }v/a_{k1}& L^{\prime }\end{array}\right)\left(\begin{array}{cc}{a}_{k1}& w^T\\ 0& U^{\prime }\end{array}\right) \\ =LU \end{gathered}$$
   这么复杂的过程，我们看看就行了。又不是立志做数学家，数学家才深究证明过程和计算过程，我们是程序员，奉行拿来主义，把数学家的结论拿过来用就行了。
   所以LUP分解的步骤是：
   1. 找出第一列绝对值最大的一行，然后和第一行置换，记录置换；
   2. 对置换后的矩阵进行LU分解，计算舒尔补；
   3. 将舒尔补视为新的新的矩阵，重复进行这个LUP分解过程；
   4. 待矩阵变成$1\times 1$矩阵时，结束。

Python实现

# LUP分解
def lup_decomposition(self):

    n = len(self.__lines)
    # 初始化L U为同样大小的0矩阵
    l = copy.deepcopy(self.__lines)
    u = copy.deepcopy(self.__lines)
    a = copy.deepcopy(self.__lines)
    p = [i for i in range(0, n)]
    # i代表处理的行和列，因为每次都是从左上到右下不断缩小这个矩阵
    for i in range(0, n):
        # 先求出自己的置换
        # 注意不能用是否等于0判断，而是考虑浮点数精度进行判断
        # 找最大绝对值
        swap_line = i

        # 进行置换
        # 找出不为0的一行，进行置换
        for k in range(i + 1, n):
            if abs(a[k][i]) > abs(a[swap_line][i]):
                swap_line = k
        if abs(a[swap_line][i]) <= MAX_ERROR:
            raise Exception("矩阵无解")
        # 交换两行
        a[i], a[swap_line] = a[swap_line], a[i]
        l[i], l[swap_line] = l[swap_line], l[i]
        p[i], p[swap_line] = p[swap_line], p[i]

        # 剩下的就是LU分解了
        # 每一行的l进行赋值
        # 这个循环其实是对角线方向
        l[i][i] = 1
        u[i][i] = a[i][i]
        for j in range(i + 1, n):
            l[i][j] = 0
            l[j][i] = a[j][i] / a[i][i]
            u[i][j] = a[i][j]
            u[j][i] = 0
            # 处理A比较麻烦，需要计算舒尔补
            for k in range(i + 1, n):
                # 舒尔补是 a[j][k] -= a[j][i] * a[i][k] / a[i][i]
                # 因为前面计算l[j][i] = a[j][i] / a[i][i]
                # 所以可以直接替换
                a[j][k] = round(a[j][k] - l[j][i] * a[i][k], 2)
    return Matrix(l), Matrix(u), Matrix([[1 if j == i else 0 for j in range(0, n)] for i in p])

  解方程可以把LU分解中的代入法抽出一个新方法:

def solve(l, u, values):
    n = len(values)
    y = [0] * n
    for i in range(0, n):
        #
        # 遍历其后的所有行,让values的值减掉
        y[i] = values[i]
        for j in range(i + 1, n):
            values[j] -= l.lines[j][i] * values[i]
    # 再解Ux=y
    solution = [0] * n
    for i in range(n - 1, -1, -1):
        solution[i] = round(y[i] / u.lines[i][i], 2)
        for j in range(0, i):
            y[j] -= u.lines[j][i] * solution[i]
    return solution

  然年LUP分解解方程代码就简单了：

def solve_by_lup(self, values):
    l, u, p = self.lup_decomposition()
    # lux=pb
    pb = p * Matrix([[num] for num in values])
    pb_values = [i[0] for i in pb.__lines]
    # 因为l是单位下三角矩阵，所以直接代入法
    return solve(l, u, pb_values)

  测试代码：

if __name__ == '__main__':
    a = Matrix([[2, 0, 2, 0.6], [3, 3, 4, -2], [5, 5, 4, 2], [-1, -2, 3.4, -1]])
    l, u, p = a.lup_decomposition()
    print(l)
    print(u)
    print(p)
    print(l*u)
    print(p*a)

    solutions = a.solve_by_lup([1, 2, 3, 4])
    print(solutions)
    print(a*Vector(solutions))

  测试结果完全正确：

[1, 0, 0, 0]
[0.4, 1, 0, 0]
[-0.2, 0.5, 1, 0]
[0.6, -0.0, 0.4, 1]
[5, 5, 4, 2]
[0, -2.0, 0.4, -0.2]
[0, 0, 4.0, -0.5]
[0, 0, 0, -3.0]
[0, 0, 1, 0]
[1, 0, 0, 0]
[0, 0, 0, 1]
[0, 1, 0, 0]
[5, 5.0, 4.0, 2.0]
[2.0, 0.0, 2.0, 0.6]
[-1.0, -2.0, 3.4, -1.0]
[3.0, 3.0, 4.0, -2.0]
[5, 5, 4.0, 2.0]
[2, 0, 2.0, 0.6]
[-1, -2, 3.4, -1.0]
[3, 3, 4.0, -2.0]
[-0.91, 0.29, 1.24, 0.56]
[1.0]
[1.98]
[2.98]
[3.99]