5.1 LU分解

解方程原理

  LU分解，就是将矩阵化成单位下三角矩阵乘以上三角的形式，也就是分解为A=LU。其中L是单位下三角矩阵，U是上三角矩阵。LU分解有很多应用，其中最普遍的应用就是用来解方程，那具体怎么使用LU分解来解方程呢？
  如下是一个方程：
$$\left(\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{11}& a_{12}& \cdots & a_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ a_{11}& a_{12}& \cdots & a_{nn}\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right)=\left(\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right)$$
  将这个方程的三部分分别用符号$A,x,b$表示，就是：
$$Ax=b$$
  将分解结果$A=LU$代入，得到了：
$$LUx=b$$
  定义向量$y=Ux$，就有了那么$Ly=B$，这是个新方程，但是这个方程的系数矩阵L是下三角矩阵，所以很容易求出解$y$。然后再解$Ux=y$这个方程。因为$U$是一个上三角矩阵，所以这个方程的解就很容易出来。
  对于LU分解，$A=LU$，如果$L$是下三角矩阵，那么称为Doolittle分解。如果U是上三角矩阵，那么成为Crout分解。在实际应用中，一般使用Doolittle分解。Doolittle分解也是默认的LU分解，也就是说如果没有上下文，那么LU分解指的就是Doolittle分解。

分解原理

  如果学习了初等矩阵，就很容易理解LU分解，其基本原理就是初等行变换。假设矩阵$A$经历了一系列初等行变换变成了上三角矩阵U，那么相当于左乘了一系列初等矩阵，那么有：
$$L_n{L}_{n-1}\cdots L_2{L}_{1}A=U$$
  因为初等矩阵都是可逆的，所以把U同样左乘相反的初等行变换就可以变回A，于是有：
$$A=L_1^{-1}{L}_2^{-1}\cdots L_{n-1}^{-1}{L}_n^{-1}U$$
  所以L就是初等矩阵的乘积：
$$L=L_1^{-1}{L}_2^{-1}\cdots L_{n-1}^{-1}{L}_n^{-1}$$

舒尔补方法

  首先把矩阵拆分：
$$\begin{gathered} A=\left(\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ a_{n1}& a_{n2}& \cdots & a_{nn}\end{array}\right) \\ =\left(\begin{array}{cc}{a}_{11}& w^T\\ v& A^{\prime }\end{array}\right) \end{gathered}$$
  那么A可以分解为这个样子：
$$\begin{gathered} A=\left(\begin{array}{cc}{a}_{11}& w^T\\ v& A^{\prime }\end{array}\right) \\ =\left(\begin{array}{cc}1& 0\\ v/a_{11}& I_{n-1}\end{array}\right)\left(\begin{array}{cc}{a}_{11}& w^T\\ 0& A^{\prime }-v{w}^T/a_{11}\end{array}\right) \end{gathered}$$
  这个分解用初等矩阵理论其实很容易理解。其实用到的是行倍加变换。对于第一列的元素，把第一行乘以$-a_{i1}/a_{11}$加到第$i$行，就把第$i$行的第一个元素变成0了。那么这一系列初等矩阵的逆的乘积就是下面这个形式：
$$\left(\begin{array}{cc}1& 0\\ v/a_{11}& I_{n-1}\end{array}\right)$$
  右下角这个东西$A^{\prime }-v{w}^T/a_{11}$，就叫舒尔补schur’s complement。因为是第一行乘以倍数加到其他行，所以可以认为是将$w^T$乘以了$-v/a_{11}$加到了右下角的矩阵上。所以LU分解就是不断递归求舒尔补的过程。首先得到L和U的第一行和第一列和舒尔补。然后把舒尔补视为新矩阵，计算出舒尔补的L和U的第一行和第一列，以及舒尔补的舒尔补。直到最终的舒尔补变成$1\times 1$的矩阵。所以舒尔补方法得到的就是Doolittle分解。

伪代码

  我们大可不必递归，只需要按对角线循环，把剩余的矩阵变成舒尔补。定义两个矩阵，一个L，一个U。然后按圈层循环，0代表最左上角第一圈,用i代表圈层。每个循环先把L矩阵的左上角变成1，第一行其他的变成0。然后处理U矩阵。第一行其他变成$w^T$,第一列其他变成0,A矩阵剩余部分变成舒尔补。伪代码就是：

L[i][i]=1
For j to n-1 :
	L[i][j]=0; 
	L[j][i]=L[j][i]/a[i][i]
	U[i][i]=A[i][i]
For j from i to n-1 :
	U[i][j]=A[i][j]; 
	U[j][i]=0

  合并起来就是:

L[i][i]=1
U[i][i]=A[i][i]
for j from i+1 to n-1 :
	L[i][j]=0; L[j][i]=L[j][i]/a[i][i]
	U[i][j]=A[i][j]; U[j][i]=0

  A剩余部分变成舒尔补,需要使用矩阵乘法，但由于是$n\times 1$矩阵乘以$1\times n$矩阵，所以无法使用Strassen算法进行优化。只能遍历v的行$w^T$的列,在循环里累加，正常的矩阵乘法是三层循环。但是$v$每一行才1一个元素，$w^T$每一列也才1个元素，所以第三层循环只循环一次。这样的话就可以省略这个循环，改成两层循环。伪代码如下：

For j from i+1 to n-1
   For k from i+1 to n-1
      A[j][k]=A[j][k]-A[j][i]*A[i][k]/A[i][i]

  伪代码继续合并为：

For i from 0 to n-1
    L[i][i]=1
    U[i][i]=A[i][i]
    For j from i+1 to n-1 :
	    L[i][j]=0; L[j][i]=a[j][i]/a[i][i]
	    U[i][j]=A[i][j]; U[j][i]=0
	For k from i+1 to n-1
    	A[j][k]=A[j][k]-A[j][i]*A[i][k]/A[i][i]

Python实现

  有了伪代码，翻译为python代码就简单了：

def lu_decomposition(self):

    n = len(self.__lines)
    # 初始化L U为同样大小的0矩阵
    l = copy.deepcopy(self.__lines)
    u = copy.deepcopy(self.__lines)
    a = copy.deepcopy(self.__lines)

    for i in range(0, n):
        # 每一行的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] -= l[j][i] * a[i][k]
        print(a)
    return Matrix(l), Matrix(u)

  但是LU分解之后是为了解方程的。解方程就是直接代入。代码并不难：

def solve_by_lu(self, values):
    l, u = self.lu_decomposition()
    # ly = b
    # 因为l是单位下三角矩阵，所以直接代入法
    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] = y[i] / u.__lines[i][i]
        for j in range(0, i):
            y[j] -= u.__lines[j][i] * solution[i]

    return solution

  这里提供一个验证结果是否正确的向量类，这个类我简单地继承了Matrix类:

# _*_ coding:utf-8 _*_
# 向量
from com.youngthing.mathalgorithm.matrix import Matrix

class Vector(Matrix):

    def __init__(self, array):
        lines = [[num] for num in array]
        super().__init__(lines)

  测试代码：

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

  测试结果:

[21.875, -3.25, 32.0, -13.0]
[1.0]
[2.0]
[3.0]
[4.0]

  可见，结果完全正确。