5.2 LUP分解

大家好，欢迎来到IT知识分享网。

解方程步骤

  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]