Doolittle三角矩阵分解

有关Doolittle三角矩阵分解的思路和代码实现。

什么是Doolittle分解

首先对一个行列式不为零的矩阵来说，一定可以通过初等行变换 将 该矩阵转化成一个初等矩阵和一个上三角矩阵相乘的形式。
$$A=L\cdot U$$
此时如果L 是一个单位下三角矩阵 ，表明矩阵A 可以进行Doolittle 分解 ，当然，Doolittle分解的三角形 并不唯一；
设矩阵D是 可逆对角矩阵，则有
$$D\cdot D^{-1}=E$$,故
$$A=L\cdot D\cdot D^{-1}\cdot U$$
同时，又可变化为

$$A=(L\cdot D)\cdot (D^{-1}\cdot U)$$
就相当于，又产生了一组 L 和 U。
Doolittle 分解唯一的充要条件是 N-1 阶顺序主子式大于零，如果在n-1阶子式中有不为零的项，可通过初等行变换，消除这种现象。

使用程序进行矩阵分解。

首先将一个矩阵分解为两个矩阵的乘积可以通过待定系数法 进行实现，矩阵如下所示。

$$\left\{\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right\}=\left\{\begin{array}{ccc}1& 0& 0\\ l_{21}& 1& 0\\ l_{31}& l_{32}& 1\end{array}\right\}\cdot \left\{\begin{array}{ccc}{r}_{11}& r_{12}& r_{13}\\ 0& r_{22}& r_{23}\\ 0& 0& r_{33}\end{array}\right\}$$
容易求得 矩阵L的第一列和矩阵U的第一行比较容易计算，即 满足如下公式

{ r 1 , j = a 1 , j j ⊂ { 1 , 2 , 3... n } l i , 1 = a i , 1 a 1 , 1 j \begin{cases} r_{1 ,j} =a_{1,j} & j \subset \{1,2,3...n\} \\ l_{i,1}= \frac{a_{i,1}}{a_{1,1}} & j \end{cases} {r1,j​=a1,j​li,1​=a1,1​ai,1​​​j⊂{1,2,3...n}j​
同时此公式可以作为矩阵的初始化 使用。
下面开始计算非特殊位置
{ r i , j = a i , j − ∑ m = 1 i − 1 l i , m ⋅ r m , j i ⊂ { 2 , 3 , . . n } , j ⊂ { i , i + 1... n } l j , i = a j , i − ∑ m = 1 i − 1 l j , m ⋅ r m , i i ⊂ { 2 , 3 , . . n − 1 } , j ⊂ { i + 1... n } \begin{cases} r_{i,j}=a_{i,j}-\sum_{m=1}^{i-1}l_{i,m}\cdot r_{m,j} &i \subset\{2,3,..n\},j \subset\{i,i+1...n\}\\ \\ l_{j,i}=a_{j,i}-\sum_{m=1}^{i-1}l_{j,m}\cdot r_{m,i} &i \subset\{2,3,..n-1\},j \subset\{i+1...n\}\\ \end{cases} ⎩⎪⎨⎪⎧​ri,j​=ai,j​−∑m=1i−1​li,m​⋅rm,j​lj,i​=aj,i​−∑m=1i−1​lj,m​⋅rm,i​​i⊂{2,3,..n},j⊂{i,i+1...n}i⊂{2,3,..n−1},j⊂{i+1...n}​
经过分析可知 ，求元素l或r的任意元素，只和其所在行列的上层元素有关，如
将l和u矩阵合并在一起，就会有如下结果
$$\left\{\begin{array}{ccc}{r}_{11}& r_{12}& r_{13}\\ l_{21}& r_{22}& r_{23}\\ l_{31}& l_{32}& r_{33}\end{array}\right\}$$
如果求$r_{22}$ 就需要知道$r_{12}$和$l_{21}$的值，即
所以说，我们可以用迭代的方式求出两个矩阵的所有元素。

同时，i和j，k的元素有重叠的部分，可以放在一个for循环之下执行。
{ i ⊂ { 2 , 3 , . . n } j ⊂ { i , i + 1... n } \begin{cases} i \subset\{2,3,..n\}\\j \subset\{i,i+1...n\} \end{cases} {i⊂{2,3,..n}j⊂{i,i+1...n}​
{ i ⊂ { 2 , 3 , . . n − 1 } , j ⊂ { i + 1... n } \begin{cases} i \subset\{2,3,..n-1\},\\j \subset\{i+1...n\} \end{cases} {i⊂{2,3,..n−1},j⊂{i+1...n}​

取重合部分 ，$i\subset [2,3,..n-1],j\subset [i+1,i+2.....n]$ 进行大循环，同时我们可以使用if 判断 来去除一些元素

#初始化 第一行 和第一列
lu[0]=A[0]
for i in range(1,n):# n 是矩阵的阶数
    lu[i][0]=A[i][0]/A[0][0]
# 初始化结束，求其他数值
for  i in range(1,n):
    for j in range(i,n):
        lu[i][j]=A[i][j]-sum([lu[i][k]*lu[k][j] for k in range(0,i)])
        if i< n-1 and j!=i:# 通过判断排除 这里是由ij 代替了k
            lu[j][i]=(A[j][i]-sum([lu[j][k]*lu[k][i] for k in range(0,i)]))/lu[i][i]
     

完整代码

# Doolittle.py
A=[[2,2,-5],[1,-3,1],[1,5,2]]
A=[[2,-1,0,1],[1,3,0,1],[0,1,-1,2],[1,0,0,1]]
n=len(A)
# //创建一个和矩阵A同阶的矩阵
lu=[[0 for i in range(n)] for j in range(n)]
# 初始化
lu[0]=A[0]
for i in range(1,n):
    lu[i][0]=A[i][0]/A[0][0]
# 初始化结束
for  i in range(1,n):
    for j in range(i,n):
        lu[i][j]=A[i][j]-sum([lu[i][k]*lu[k][j] for k in range(0,i)])
        if i< 4-1 and j!=i:
            lu[j][i]=(A[j][i]-sum([lu[j][k]*lu[k][i] for k in range(0,i)]))/lu[i][i]


#将大矩阵可视化分解
L=[[0 for i in range(n)] for j in range(n)]
U=[[0 for i in range(n)] for j in range(n)]
for  i in range(0,n):
    for j in range(0,n):
        if i==j:
            L[i][j]=1
            U[i][j]=lu[i][j]
        if i>j:
            L[i][j]=lu[i][j]
        else:
            U[i][j]=lu[i][j]
print("原矩阵")
print("\n".join(["".join(str(i)) for i in A]))
print("矩阵L")
print("\n".join(["".join(str(i)) for i in L]))
print("矩阵U")
print("\n".join(["".join(str(i)) for i in U]))

测试用例

原矩阵
[2, 2, -5]
[1, -3, 1]
[1, 5, 2]
矩阵L
[1, 0, 0]
[0.5, 1, 0]
[0.5, -1.0, 1]
矩阵U
[2, 2, -5]
[0, -4.0, 3.5]
[0, 0, 8.0]

原矩阵
[2, -1, 0, 1]
[1, 3, 0, 1]
[0, 1, -1, 2]
[1, 0, 0, 1]
矩阵L
[1, 0, 0, 0]
[0.5, 1, 0, 0]
[0.0, 0.2857142857142857, 1, 0]
[0.5, 0.14285714285714285, -0.0, 1]
矩阵U
[2, -1, 0, 1]
[0, 3.5, 0.0, 0.5]
[0, 0, -1.0, 1.8571428571428572]
[0, 0, 0, 0.4285714285714286]