基于Python实现的矩阵Doolittle分解及应用

若方阵A的各阶主子式不为0，则矩阵$A$可分解为$A=LU$，$L$为下三角矩阵，$U$为上三角矩阵，这种方法称为矩阵的Doolittle分解。

主元直接分解

原理

设 $A=LU$，则有
$$\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \dots & a_{2n}\\ ⋮& ⋮& & ⋮\\ a_{n1}& a_{n2}& \dots & a_{nn}\end{array}\right]=\left[\begin{array}{cccc}1& 0& \cdots & 0\\ l_{21}& 1& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ l_{n1}& l_{n2}& \dots & 1\end{array}\right]\left[\begin{array}{cccc}{u}_{11}& u_{12}& \cdots & u_{1n}\\ 0& u_{22}& \dots & u_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & u_{nn}\end{array}\right]$$
根据矩阵乘法可得
$$a_{ij}=\sum _{m=1}^n{l}_{im}{u}_{mj},i,j=1,2,\dots ,n$$
可以发现 $l_{ii}=1，$且L为下三角矩阵，U为上三角矩阵，利用这些信息可以得到
$$\begin{gathered} a_{ij}=\sum _{m=1}^{i-1}{l}_{im}{u}_{mj}+u_{ij} \\ {a}_{ij}=\sum _{m=1}^{j-1}{l}_{im}{u}_{mj}+l_{ij}{u}_{jj} \end{gathered}$$

就可以计算L和U了，具体公式如下

计算 $u_{ij}$
$$u_{ij}=a_{ij}-\sum _{m=1}^{i-1}{l}_{im}{u}_{mj}$$
计算 $l_{ij}$
$$l_{ij}=\frac{a_{ij}-{\sum }_{m=1}^{j-1}{l}_{im}{u}_{mj}}{u_{jj}}$$

代码

def Dolittle(Matrix):
#事实上，w、L、U可以用一个矩阵存储，为了方便理解，写成了三个矩阵
    w = Matrix.shape[0]
    L = np.zeros((w,w))
    U = np.zeros((w,w))

    #U的第一行
    U[0,:] = Matrix[0,:]
    #L的主元和第一列
    for i in range(w):
        L[i,i] = 1
        L[i,0] = Matrix[i,0] / U[0,0]
    for i in range(1,w):
        for j in range(i,w):
            #U的第i行
            U[i,j] = Matrix[i,j] - np.dot(L[i,:i],U[:i,j])
            #L的第i列
            if(j+1<w):
                L[j+1,i] = (Matrix[j+1,i] - np.dot(L[j+1,:i],U[:i,i]))/U[i,i]
    return L,U

实验

Matrix = np.array([
    [2,1,1],[4,4,3],[6,7,7]
])
L,U=Dolittle(Matrix)

结果

原矩阵：

L(分解得到的下三角矩阵)：

U(分解得到的上三角矩阵)：

主元直接分解解方程

直接上公式
$$AX=LUX=b$$令
$$\{\begin{array}{l}Ly=b\\ Ux=y\end{array}$$
解方程组就完事了

代码

def SolveEquations(L,U,b):
    w = L.shape[0]
    x = np.zeros(w)
    y = np.zeros(w)
    for i in range(w):
        y[i] = (b[i] - np.dot(L[i,:],y.T))
    for i in range(w-1,-1,-1):
        x[i] = (y[i] - np.dot(U[i,:],x.T))/U[i,i]
    return x

列主元直接三角分解

上述三角消元法源于顺序Gauss消元法，没有选择主元，可能产声数值不稳定的现象，而列主元直接三角分解解决了这个问题。列主元三角分解的思想与列主元Guass消元法类似，需要选主元。

事实上，计算机在进行数值计算的过程中由于数字位数限制，计算的结果都是近似值，存在一定的误差。而在进行直接三角分解的过程中，$U$ 的主元 $U_{kk}$ 会在计算 $l_{ik}$ 时作为分母，如果 $U_{kk}$ 太小，会使 $l_{ik}$ 很大，从而导致整个计算过程中的累计误差很大。

所以，计算矩阵 $U$ 的第 $k$ 行时，我们考察 $U_{kk},U_{k+1k},\dots ,U_{nk}$等元素，将该值最大的一行与第 $k$ 行交换，然后计算其他元素的值。

在计算过程中矩阵会发生变换，此时不满足 $A=LU$，如果仍要写成等式的样子需要引入行变换矩阵 $P$ 。即
$$PA=LU$$
代码

def Dolittle_Col(Matrix):
#事实上，w、L、U可以用一个矩阵存储，为了方便理解，写成了三个矩阵
    w = Matrix.shape[0]
    L = np.zeros((w,w))
    U = np.zeros((w,w))

    #U的第一行
    U[0,:] = Matrix[0,:]
    #L的主元和第一列
    for i in range(w):
        L[i,i] = 1
        L[i,0] = Matrix[i,0] / U[0,0]
    for i in range(1,w):
        #交换主元最大的行到第i行
        maxrow = i
        maxValue = float('-inf')
        for k in range(i,w):
            value = Matrix[k,i] - np.dot(L[k,:i],U[:i,i])
            if(maxValue<value):
                maxrow = k
                maxValue = value
        Matrix[[i,maxrow],:] = Matrix[[maxrow,i],:] #交换
        L[[i,maxrow],:i] = L[[maxrow,i],:i] #交换
        for j in range(i,w):
            #U的第i行
            U[i,j] = Matrix[i,j] - np.dot(L[i,:i],U[:i,j])
            #L的第i列
            if(j+1<w):
                L[j+1,i] = (Matrix[j+1,i] - np.dot(L[j+1,:i],U[:i,i]))/U[i,i]
    return L,U

实验

Matrix = np.array([
    [2,1,1],[4,4,3],[6,7,7]
])
L,U=Dolittle_Col(Matrix)

结果
原矩阵：

L(分解得到的下三角矩阵)：

U(分解得到的上三角矩阵)：

上面结果并不满足 $A=LU$，但是满足 $PA=LU$， L与U相乘的结果是A第二行与第三行互换后的样子。

主元直接分解求逆矩阵

设
$$X=A^{-1}$$
则有
$$AX=E$$
令
$$\begin{gathered} X=(x_1,x_2,\dots ,x_n) \\ E=(e_1,e_2,\dots ,e_n) \end{gathered}$$
有
$$A{x}_i=e_i$$
对系数矩阵A使用直接三角分解，然后利用三角分解易解同系数矩阵的性质，解得每一个 $x_i$ ，组成逆矩阵。

代码

def GetInverseMatrix(matrix):
    w = matrix.shape[0]
    L,U = Dolittle(matrix)
    inverseMatrix = np.zeros((w,w))
    for i in range(w):
        e = np.zeros(w)
        e[i]=1
        inverseMatrix[:,i] = SolveEquations(L,U,e)
    return inverseMatrix

实验
随便在百度找了一道求逆矩阵的题

A = np.array([[2,3,1],[0,1,3],[1,2,5]])
inverseMatrix = GetInverseMatrix(A)
print(inverseMatrix)

结果

答案

很好用有木有。