基于Python实现的Cholesky分解与Crout分解

阅读这篇博客需要有矩阵Doolittle分解的理论基础。

---

Cholesky分解

对于对称正定矩阵，如果 $A$ 的所有顺序主子式均大于0，则该矩阵可分解为 $A=LU$ 。其中，$U$ 为上三角矩阵，并分解为 $U=DM$ , $D$ 为对角矩阵 $diag(u_{11},u_{22},\dots ,u_{nn})$，M为单位上三角矩阵。

利用矩阵的对称性和Doolittle分解的唯一性进行如下推导
$$\begin{gathered} A=LDM=A^T=(LDM{)}^T=M^{T}D{L}^T \\ M=L^T(分解的唯一性) \\ A=LD{L}^T \end{gathered}$$
因为A是正定的，所以若令 $D^{\frac{1}{2}}=diag(\sqrt{u_{11}},\sqrt{u_{22}},\dots ,\sqrt{u_{nn}})$，有
$$\begin{gathered} D=D^{\frac{1}{2}}{D^{\frac{1}{2}}}^T \\ 令G=L{D}^{\frac{1}{2}} \\ A=G{G}^T \end{gathered}$$
G为下三角矩阵，上述分解为正定矩阵的Cholesky分解

代码

#为方便示意不考虑时间空间代价
def Cholesky(matrix):
    w = matrix.shape[0]
    G = np.zeros((w,w))#实际上只用一半的空间就可以完成矩阵分解
    for i in range(w):
        G[i,i] = (matrix[i,i] - np.dot(G[i,:i],G[i,:i].T))**0.5
        for j in range(i+1,w):
            G[j,i] = (matrix[j,i] - np.dot(G[j,:i],G[i,:i].T))/G[i,i]
    return G

实验

C = np.array([[4,-2,4],[-2,5,0],[4,0,6]])
print(Cholesky(C))

结果
原矩阵

G

由于G的对称性，Cholesky分解所需的计算量和存储量仅为Doolittle分解的一半。

避免开平方的改进

在Cholesky分解的过程中，G的主元需要通过开平方运算，对于早期的计算机，开平方运算很耗费资源，所以人们在Cholesky分解的基础上进行改进，并没有引入 $G$ 而是将矩阵分解为如下的形式
$$A=LD{L}^T$$
代码

#为方便示意不考虑时间空间代价
def Cholesky_plus(matrix):
    w = matrix.shape[0]
    L = np.zeros((w,w))
    for i in range(w):
        L[i,i] = 1
    D = np.zeros((w,w))
    for i in range(w):
        D[i,i] = matrix[i,i] - np.dot(np.dot(L[i,:i],D[:i,:i]),L[i,:i].T)
        for j in range(i+1,w):
            L[j,i] = (matrix[j,i] - np.dot(np.dot(L[j,:i],D[:i,:i]),L[i,:i].T))/D[i,i]
    return L,D

实验

L,D = Cholesky_plus(C)
print(L)
print(D)
print(np.dot(np.dot(L,D),L.T))

结果
原矩阵

$L、D$ 以及 $LD{L}^T$

Crout分解

很简单的
$$\begin{gathered} A=LDM \\ 令T=LD \\ A=TM \end{gathered}$$
这个 $A=TM$ 就是Crout分解，其中T为下三角矩阵，M为单位上三角矩阵。嗯，和Doolittle分解很类似。。。。。

代码

def Crout(matrix):
    w = matrix.shape[0]
    T = np.zeros((w,w))
    M = np.zeros((w,w))
    for i in range(w):
        M[i,i] = 1
    for i in range(w):
        for j in range(i,w):
            T[j,i] = matrix[j,i] - np.dot(T[j,:i],M[:i,i])
        for j in range(i+1,w):
            M[i,j] = (matrix[i,j] - np.dot(T[i,:i],M[:i,j]))/T[i,i]
    return T,M

实验

A = np.array([[2,-1,0,0],[-1,2,-1,0],[0,-1,2,-1],[0,0,-1,2]])
T,M = Crout(A)
print(T)
print(M)

结果
原矩阵

T

M