精华内容
下载资源
问答
  • 矩阵三角分解

    2014-12-03 22:59:01
    基于高斯消去的非奇异矩阵三角分解,将矩阵分解成上三角矩阵U和下三角矩阵L
  • (1)编制解n阶线性方程组Ax=b的列主元Gauss消去法的通用程序;
  • 作了比较全面的统述,它包括矩阵三角分解存在唯一性的充要条件、存在性的充要条件等.并介绍了三类特殊矩阵的 三角分解,且对矩阵三角分解的计算及应用作了较详细的论述.最后,本文特别例举了比较有意思的简单多项式...
  • 实现了矩阵的LU分解,也称三角分解,输入为n阶方阵A,输出为L和U,同时会自动进行矩阵是否可以进行LU分解的判断,注释翔实,有一定基础的一定可以轻松看懂。
  • 多核与多GPU系统下的一种矩阵三角分解并行算法.pdf
  • 解线性方程组地矩阵三角分解法* 第五章解线性方程组的直接方法 计算方法 —— 矩阵三角分解法 * 本讲内容 一般线性方程组 LU 分解与 PLU 分解 对称正定线性方程组 平方根法--Cholesky 分解 对角占优三对角线性方程...

    解线性方程组地矩阵三角分解法

    * 第五章解线性方程组的直接方法 计算方法 —— 矩阵三角分解法 * 本讲内容 一般线性方程组 LU 分解与 PLU 分解 对称正定线性方程组 平方根法--Cholesky 分解 对角占优三对角线性方程组 追赶法 * LU 分解 将一个矩阵分解成结构简单的三角形矩阵的乘积 矩阵的三角分解 矩阵的 LU(Doolittle) 分解 矩阵的 LDR 分解 克洛脱 (Crout) 分解 * 计算 LU 分解 利用矩阵乘法直接计算 LU 分解 L ? U = A ?比较等式两边的第一行得: u1j = a1j 比较等式两边的第一列得: ?比较等式两边的第二行得: 比较等式两边的第二列得: ( j = 1,…, n ) ( i = 2,…, n ) ( j = 2,…, n ) ( i = 3,…, n ) U 的第一行 L 的第一列 U 的第二行 L 的第二列 * 计算 LU 分解 第 k 步:此时 U 的前 k-1 行和 L 的前 k-1 列已经求出 直到第 n 步,便可求出矩阵 L 和 U 的所有元素。 比较等式两边的第 k 行得: ( j = k, …, n ) 比较等式两边的第 k 列得: ( i = k+1, …, n ) * LU 分解算法 算法 :(LU 分解 ) for k = 1 to n end j = k, …, n i = k+1, …, n Matlab程序参见:ex51.m 乘除法运算量:(n3 - n)/3 为了节省存储空间,通常用 A 的绝对下三角部分来存放 L (对角线元素无需存储),用 A 的上三角部分来存放 U * PLU 分解 矩阵的 PLU 分解 for k = 1 to n end i = k, k+1, …, n j = 1, 2, …, n i = k+1, …, n j = k+1, …, n Matlab程序:上机练习 * Cholesky 分解 对称正定矩阵的三角分解--Cholesky 分解 定理:设 A 是对称矩阵,若 A 的所有顺序主子式都不为 0,则 A 可唯一分解为 其中 L 为单位下三角阵,D 为对角矩阵 A = LDLT 定理:(Cholesky分解)若 A 对称正定,则 A 可唯一分解为 其中 L 为下三角实矩阵,且对角元素都大于 0 A = LLT * 计算 Cholesky 分解 Cholesky 分解的计算 直接比较等式两边的元素 计算公式 * Cholesky 分解算法 for j = 1 to n end i = j +1, …, n 算法 :(Cholesky 分解 ) * 平方根法 A 对称正定 算法 :(解对称正定线性方程组的平方根法 ) 计算 A 的 Cholesky 分解 解方程:Ly = b 和 LTx = y i = 2, 3, …, n i = n-1, …, 2, 1 * 改进的 Cholesky 分解 计算公式 改进的 Cholesky 分解 * 改进的 Cholesky 分解 for j = 1 to n end i = j +1, …, n 算法 :(改进的 Cholesky 分解 ) 优点:避免开方运算 * 改进的平方根法 A 对称正定 算法 :(解对称正定线性方程组的改进的平方根法 ) 计算 改进的 Cholesky 分解 解方程:Ly = b 和 DLTx = y i = 2, 3, …, n i = n-1, …, 2, 1 * 追赶法 对角占优的三对角矩阵的 LU 分解 计算公式 i = 2, 3, …, n-1 * 追赶法 A 三对角矩阵(对角占优) 算法 :(追赶法 ) i = 2, 3, …, n i = n-1, …, 2, 1 i = 2, 3, …, n-1 运算量:5n-4 2n – 3 次 2n 次 n – 1 次

    展开全文
  • * 一、直接法概述 直接法是将原方程组化为一个或若干个三角形 方程组的方法,共有若干种. 对于线性方程组 其中 系数矩阵 未知量向量 常数项 根据Cramer...则 都是三角 形方程组 上述方法称为直接三角形分解法 §2...

    * 一、直接法概述 直接法是将原方程组化为一个或若干个三角形 方程组的方法,共有若干种. 对于线性方程组 其中 系数矩阵 未知量向量 常数项 根据Cramer(克莱姆)法则,若 determinantal 行列式的记号 若用初等变换法求解,则对其增广矩阵作行初等变换: 经过n-1次 同解 即 以上求解线性方程组的方法称为Gauss消去法 则 都是三角 形方程组 上述方法称为直接三角形分解法 §2 Matrix Factorization – Doolittle ? 道立特分解法 /* Doolittle Factorization */: —— LU 分解的紧凑格式 /* compact form */ 反复计算, 很浪费哦 …… 通过比较法直接导出L 和 U 的计算公式。 思路 §2 Matrix Factorization – Doolittle 固定 i : 对 j = i, i+1, …, n 有 lii = 1 a 固定 j ,对 i = j, j+1, …, n 有 b 上述解线性方程组的方法称为 直接三角分解法的 Doolittle法 例1. 用Doolittle法解方程组 解: 由Doolittle分解 Doolittle法在计算机上实现是比较容易的 但如果按上述流程运算仍需要较大的存储空间: 因此可按下列方法存储数据: 直接三角分解的Doolittle法可以用以下过程表示: 存储单元(位置) 紧凑格式的 Doolittle法 例2. 用紧凑格式的Doolittle法解方程组(例1) 解: 所以 Matrix Factorization – Choleski ? 平方根法 /* Choleski’s Method */: ——对称 /* symmetric */ 正定 /* positive definite */ 矩阵的分解法 定义 一个矩阵 A = ( aij )n?n 称为对称阵,如果 aij = aji 。 定义 一个矩阵 A 称为正定阵,如果 对任意非零向量 都成立。 ?回顾:对称正定阵的几个重要性质 ? A?1 亦对称正定,且 aii > 0 ? 若不然,则 存在非零解,即 存在非零解。 ? ? 对任意 , 存在 , 使得 , 即 。 ? 其中 第 i 位 ? A 的顺序主子阵 /* leading principal submatrices */ Ak 亦对 称正定 对称性显然。对任意 有 , 其中 。 ? A 的特征值 /* eigen value */ ?i > 0 设对应特征值 ? 的非零特征向量 为 ,则 。 ? A 的全部顺序主子式 det ( Ak ) > 0 因为 一、对称正定矩阵的三角分解(Cholesky分解) 记为 Diagonal:对角 因此 所以 综合以上分析, 则有 定理1. (Cholesky分解) 且该分解式唯一 这种关于对称正定矩阵的分解称为Cholesky分解 二、对称正定线性方程组的解法 线性方程组 则线性方程组(10)可化为两个三角形方程组 对称正定方程 组的平方根法 例1. 用平方根法解对称正定方程组 解: 即 三、平方根法的数值稳定性 用平方根法求解对称正定方程组时不需选取主元 由 可知 因此 平方根法是数值稳定的 事实上,对称正定方程组也可以用顺序Gauss消去法求解 而不必加入选主元步骤 §2 Matrix Factorization – Tridiagonal System ? 追赶法解三对角方程组 /* Crout Reduction for Tridiagonal Linear System */ Step 1: 对 A 作Crout 分解 直接比较等式两边的元素,可得到计算公式。 Step 2: 追——即解 : Step 3: 赶——即解 : 与G.E.类似,一旦 ?i = 0 则算法中断,故并非任何 三对角阵都可以用此方法分解。

    展开全文
  • 矩阵三角分解-LDR分解

    2012-09-25 16:01:49
    线性代数课本中给出了矩阵的LDR三角分解的定义形式以及相应的定理,但一般要求所给矩阵是奇异的方阵。该文档适当放宽条件,给予了证明。
  • 矩阵三角分解法(LU分解)

    万次阅读 2018-01-13 16:52:42
    矩阵分解法是高斯消元法的变形,它的复杂度和高斯消元法一样都是O(n^3),但是矩阵分解法在处理线性方程组系(具有相同的系数矩阵,但是右端项不同的方程组)时,运算比较方便。  下面是矩阵分解原理的原理:   ...

            矩阵分解法是高斯消元法的变形,它的复杂度和高斯消元法一样都是O(n^3),但是矩阵分解法在处理线性方程组系(具有相同的系数矩阵,但是右端项不同的方程组)时,运算比较方便。

            下面是矩阵分解原理的原理:

       

        

       

         下面是如何来求解L和U矩阵:

         

        在求L和U矩阵的时候,要注意两点:

       1>先求U矩阵中的一行,然后在求L矩阵的一列。次序不能颠倒。

       2>不论求L和U矩阵,都要用到相应的A矩阵中的数值。

       下面是求LU矩阵的python实现。

    data = [[2, 2, 3], [4, 7, 7], [-2, 4, 5]]
    
    
    #在实施LU分解的时候,所有的操作都在data矩阵上进行。因为LU分解的过程决定,它是一行一行进行分解的,所以用完的行可以被
    #LU矩阵中的值所代替。
    i = 0
    size = len(data)
    while i < size:
        if i == 0:
            j = 1
            while j < size:
                data[j][0]=data[j][0]/data[0][0]
                j += 1
        else:
            #下面是对U矩阵进行操作,操作过程是,先求一行U矩阵,后求一列L矩阵
            j = i
            while j < size:
                sum_column = 0
                flag_sum = i - 1
                while flag_sum >= 0:
                    if j == i:
                        sum_column += data[flag_sum][j]*data[j][flag_sum]
                    else:
                        sum_column += data[flag_sum][j]*data[j-1][flag_sum]
                    flag_sum -= 1
                data[i][j] = data[i][j]-sum_column
                j += 1
    
            #下面的是对L矩阵进行操作,操作过程是求一列L矩阵
            m = i+1
            while m <size:
                sum_column_L=0
                flag_sum_L = i-1
                while flag_sum_L >= 0:
                    sum_column_L += data[i][flag_sum_L]*data[m][flag_sum_L]
                    flag_sum_L -= 1
                data[m][i] = (data[m][i]-sum_column_L)/data[i][i]
                m += 1
        i += 1
    
    
    "输出LU矩阵"
    L = []
    U = []
    l = 0
    while l < size:
        r = 0
        temp=[]
        temp1=[]
        while r < size:
            if l > r:
                temp.append(data[l][r])
                temp1.append(0)
            else:
                if l == r:
                    temp.append(1)
                    temp1.append((data[l][r]))
                else:
                    temp.append(0)
                    temp1.append(data[l][r])
            r += 1
        L.append(temp)
        U.append(temp1)
        l += 1
    
    print("L矩阵为:\n")
    for x in L:
        print(x)
    print("\n")
    
    print("U矩阵为:\n")
    for x in U:
        print(x)
    
    
    


    展开全文
  • 1. 矩阵三角分解法 实现代码 2. LU分解 2.1 基本步骤 2.2 LU分解的计算公式 2.3 LU分解的结果表示 实现代码 3. 选主元的LU分解 3.1 基本原理 3.2 算法过程描述 实现代码 4. 对称正定矩阵的Cholesky分解 4.1 基本原理...

    \qquad 假设线性方程组

    { a 00 x 0 + a 01 x 1 + ⋯ + a 0 , n − 1 x n − 1 = b 0 a 10 x 0 + a 11 x 1 + ⋯ + a 1 , n − 1 x n − 1 = b 1 ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ a n − 1 , 0 x 0 + a n − 1 , 1 x 1 + ⋯ + a n − 1 , n − 1 x n − 1 = b n − 1 \qquad\qquad\left\{\begin{aligned}a_{\scriptscriptstyle00}x_{\scriptscriptstyle0}&+a_{\scriptscriptstyle01}x_{\scriptscriptstyle1}+\cdots+a_{\scriptscriptstyle0,n-1}x_{\scriptscriptstyle n-1}=b_{\scriptscriptstyle 0}\\ a_{\scriptscriptstyle10}x_{\scriptscriptstyle0}&+a_{\scriptscriptstyle11}x_{\scriptscriptstyle1}+\cdots+a_{\scriptscriptstyle1,n-1}x_{\scriptscriptstyle n-1}=b_{\scriptscriptstyle 1}\\ \cdots\cdots&\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\\ a_{\scriptscriptstyle n-1,0}x_{\scriptscriptstyle0}&+a_{\scriptscriptstyle n-1,1}x_{\scriptscriptstyle1}+\cdots+a_{\scriptscriptstyle n-1,n-1}x_{\scriptscriptstyle n-1}=b_{\scriptscriptstyle n-1}\end{aligned}\right. a00x0a10x0an1,0x0+a01x1++a0,n1xn1=b0+a11x1++a1,n1xn1=b1+an1,1x1++an1,n1xn1=bn1

    \qquad 写成矩阵形式 A X = b AX=b AX=b,也就是

    [ a 00 a 01 ⋯ a 0 , n − 1 a 10 a 11 ⋯ a 1 , n − 1 ⋮ ⋮ ⋮ a n − 1 , 0 a n − 1 , 1 ⋯ a n − 1 , n − 1 ] ⏟ A [ x 0 x 1 ⋮ x n − 1 ] ⏟ X = [ b 0 b 1 ⋮ b n − 1 ] ⏟ b \qquad\qquad\underbrace{\begin{bmatrix} a_{\scriptscriptstyle 00}&a_{\scriptscriptstyle 01}&\cdots&a_{\scriptscriptstyle 0,n-1} \\ a_{\scriptscriptstyle 10}&a_{\scriptscriptstyle 11}&\cdots&a_{\scriptscriptstyle 1,n-1} \\ \vdots&\vdots& &\vdots\\ a_{\scriptscriptstyle n-1,0}&a_{\scriptscriptstyle n-1,1}&\cdots&a_{\scriptscriptstyle n-1,n-1} \end{bmatrix}}_A\underbrace{\begin{bmatrix}x_{\scriptscriptstyle0}\\x_{\scriptscriptstyle1}\\\vdots\\x_{\scriptscriptstyle n-1}\end{bmatrix}}_X=\underbrace{\begin{bmatrix}b_{\scriptscriptstyle0}\\b_{\scriptscriptstyle1}\\\vdots\\b_{\scriptscriptstyle n-1}\end{bmatrix}}_b A a00a10an1,0a01a11an1,1a0,n1a1,n1an1,n1X x0x1xn1=b b0b1bn1
    \qquad

    1. 矩阵三角分解法

    \qquad 矩阵三角分解法,是将线性方程组的系数矩阵 A A A 分解为 A = L U A=LU A=LU,其中 L L L 为下三角矩阵, U U U 为上三角矩阵:

    [ a 00 a 01 ⋯ a 0 , n − 1 a 10 a 11 ⋯ a 1 , n − 1 ⋮ ⋮ ⋱ ⋮ a n − 1 , 0 a n − 1 , 1 ⋯ a n − 1 , n − 1 ] ⏟ A = [ 1 l 10 1 ⋮ ⋮ ⋱ l n − 1 , 0 l n − 1 , 1 ⋯ 1 ] ⏟ L [ u 00 u 01 ⋯ u 0 , n − 1 u 11 ⋯ u 1 , n − 1 ⋱ ⋮ u n − 1 , n − 1 ] ⏟ U \qquad\qquad\underbrace{\begin{bmatrix} a_{\scriptscriptstyle 00}&a_{\scriptscriptstyle 01}&\cdots&a_{\scriptscriptstyle 0,n-1} \\ a_{\scriptscriptstyle 10}&a_{\scriptscriptstyle 11}&\cdots&a_{\scriptscriptstyle 1,n-1} \\ \vdots&\vdots&\ddots&\vdots\\ a_{\scriptscriptstyle n-1,0}&a_{\scriptscriptstyle n-1,1}&\cdots&a_{\scriptscriptstyle n-1,n-1} \end{bmatrix}}_A =\underbrace{\begin{bmatrix} 1 \\ l_{\scriptscriptstyle 10}&1 \\ \vdots&\vdots&\ddots&\\ l_{\scriptscriptstyle n-1,0}&l_{\scriptscriptstyle n-1,1}&\cdots&1 \end{bmatrix}}_L \underbrace{\begin{bmatrix}u_{\scriptscriptstyle 00}&u_{\scriptscriptstyle 01}&\cdots&u_{\scriptscriptstyle 0,n-1} \\ &u_{\scriptscriptstyle 11}&\cdots&u_{\scriptscriptstyle 1,n-1} \\ & &\ddots&\vdots\\ & & &u_{\scriptscriptstyle n-1,n-1} \end{bmatrix}}_U A a00a10an1,0a01a11an1,1a0,n1a1,n1an1,n1=L 1l10ln1,01ln1,11U u00u01u11u0,n1u1,n1un1,n1

    \qquad 因此,解线性方程组 A x = b Ax=b Ax=b 的过程,可以分为两步:
    ( 1 )   L Y = b \qquad(1)\ LY=b (1) LY=b,求出向量 Y Y Y
    ( 2 )   U X = Y \qquad(2)\ UX=Y (2) UX=Y ,求出解 X X X

    \qquad 由于 L L L 为下三角矩阵, U U U 为上三角矩阵,线性方程组 L Y = b , U X = Y LY=b,UX=Y LY=b,UX=Y 的求解只要采用类似“回代”的方式求解即可(不需要消去的过程)。
    \qquad

    实现代码

    def solveLineq(L,U,b):
        # Ly = b
        rows = len(b)
        y = np.zeros(rows)
        y[0] = b[0]/L[0,0]
        for k in range(1,rows):
            y[k] = (b[k] - np.sum(L[k,:k]*y[:k]))/L[k,k]   
        # Ux = y (back substitution)     
        x = np.zeros(rows)
        k = rows-1
        x[k] = y[k]/U[k,k] 
        for k in range(rows-2,-1,-1):
            x[k] = (y[k] - np.sum(x[k+1:]*U[k,k+1:]))/U[k,k]    
        return x
    

    \qquad

    2. LU分解

    2.1 基本步骤

    第   0   步 \textbf{第 0 步}  0 

    \qquad 由于下三角矩阵 L L L 的对角线元素为 1 1 1,而 U U U 为上三角矩阵,根据矩阵乘法 A = L U A=LU A=LU,可以按以下顺序可以首先求出(如图 1 1 1):

    ( 1 ) \qquad(1) (1) L L L 的对角线元素值 :   l r r = 1 l_{rr}=1 lrr=1 ,   r = 0 , 1 , ⋯   , n − 1 ,\ r=0,1,\cdots,n-1 , r=0,1,,n1

    ( 2 ) \qquad(2) (2) U U U 的第 0 0 0 行元素值:  u 0 i = a 0 i u_{\scriptscriptstyle0\scriptstyle i}=a_{\scriptscriptstyle0\scriptstyle i} u0i=a0i ,   i = 0 , 1 , ⋯   , n − 1 ,\ i=0,1,\cdots,n-1 , i=0,1,,n1

    ( 3 ) \qquad(3) (3) L L L 的第 0 0 0 列元素值:  l r 0 = a r 0 u 00 l_{r\scriptscriptstyle0}=\dfrac{a_{r\scriptscriptstyle0}}{u_{\scriptscriptstyle00}} lr0=u00ar0 ,   r = 1 , ⋯   , n − 1 ,\ r=1,\cdots,n-1 , r=1,,n1
    在这里插入图片描述

    图1 可直接求出的部分:(1) L L L 的对角线元素;(2) U U U 的第 0 0 0 行元素;(3) L L L 的第 0 0 0 列元素

    第   1-1   步 \textbf{第 1-1 步}  1-1 

    \qquad 在求出了图 1 1 1 标识的这部分元素值之后,可以看出,此时的下三角矩阵 L L L 的第 1 1 1 行元素均已求出,也就是:

    L [ 1 , : ] = [ l 10 , 1 , 0 , ⋯ ⋯   , 0 ⏟ 2   ∼   n − 1 ] \qquad\qquad L[1,:]=[l_{10},1,\underbrace{0,\cdots\cdots,0}_{2\ \sim\ n-1}] L[1,:]=[l10,1,2  n1 0,,0]

    在这里插入图片描述

    图2 L [ 1 , : ] L[1,:] L[1,:] 可求出 U [ 1 , : ] U[1,:] U[1,:]

    \qquad 此时,可以求出上三角矩阵 U U U 的第 1 1 1 行元素 U [ 1 , : ] = [ 0 , u 11 , ⋯   , u 1 i , ⋯   , u 1 , n − 1 ] U[1,:]=[0,u_{11},\cdots,u_{1i},\cdots,u_{1,n-1}] U[1,:]=[0,u11,,u1i,,u1,n1],也就是:

    a 1 i = l 10 u 0 i + 1 ⋅ u 1 i ⟹ u 1 i = a 1 i − l 10 u 0 i \qquad\qquad a_{1i}=l_{10}u_{0i}+1\cdot u_{1i}\quad\Longrightarrow\quad u_{1i}=a_{1i}-l_{10}u_{0i} a1i=l10u0i+1u1iu1i=a1il10u0i

    第   1-2   步 \textbf{第 1-2 步}  1-2 

    \qquad 经过了前面的步骤之后,可以看出,此时的上三角矩阵 U U U 的第 1 1 1 行元素均已求出,也就是:

    U [ : , 1 ] = [ u 01 , u 11 , 0 , ⋯ ⋯   , 0 ⏟ 2   ∼   n − 1 ] \qquad\qquad U[:,1]=[u_{01},u_{11},\underbrace{0,\cdots\cdots,0}_{2\ \sim\ n-1}] U[:,1]=[u01,u11,2  n1 0,,0]

    在这里插入图片描述

    图3 U [ : , 1 ] U[:,1] U[:,1] 可求出 L [ 2 : , 1 ] L[2:,1] L[2:,1]

    \qquad 此时,可以求出下三角矩阵 L L L 的第 1 1 1 列中的未知元素 L [ 2 : , 1 ] = [ l 21 , ⋯   , u r 1 , ⋯   , u n − 1 , 1 ] L[2:,1]=[l_{21},\cdots,u_{r1},\cdots,u_{n-1,1}] L[2:,1]=[l21,,ur1,,un1,1],也就是:

    a r 1 = l r 0 u 01 + l r 1 u 11 ⟹ l r 1 = a r 1 − l r 0 u 01 u 11 \qquad\qquad a_{r1}=l_{r0}u_{01}+l_{r1}u_{11}\quad\Longrightarrow\quad l_{r1}=\dfrac{a_{r1}-l_{r0}u_{01}}{u_{11}} ar1=lr0u01+lr1u11lr1=u11ar1lr0u01

    第   2   步 ⋯ ⋯ \textbf{第 2 步}\cdots\cdots  2 

    \qquad 在这里插入图片描述

    图4 L [ 2 , : ] L[2,:] L[2,:] 可求出 U [ 2 , : ] U[2,:] U[2,:]

    类似的过程可以继续下去 ⋯ ⋯ \cdots\cdots ,直到求出 L , U L,U L,U 矩阵的所有元素值。
    \qquad

    2.2 LU分解的计算公式

    \qquad 将上述过程归纳一下,假设 L , U L,U L,U 矩阵的元素下标满足 r ≤ i   ,   i = r , r + 1 , ⋯   , n − 1 r\le i\ ,\ i=r,r+1,\cdots,n-1 ri , i=r,r+1,,n1

    1. 先直接求出:
       
      (1) L L L 的对角线元素
      (2) U U U 的第 0 0 0 行元素
      (3) L L L 的第 0 0 0 列元素
      \qquad

    2. r r r 步(第一阶段)

      由于 L [ r , : r ] = [ l r 0 , ⋯   , l r , r − 1 ] L[r,:r]=[l_{r0},\cdots,l_{r,r-1}] L[r,:r]=[lr0,,lr,r1] 已求出,也就是:
      \qquad
      L [ r , : ] = [ l r 0 , ⋯   , l r , r − 1 , 1 , 0 , ⋯ ⋯   , 0 ⏟ r + 1   ∼   n − 1 ] \qquad\qquad L[r,:]=[l_{r0},\cdots,l_{r,r-1},1,\underbrace{0,\cdots\cdots,0}_{r+1\ \sim\ n-1}] L[r,:]=[lr0,,lr,r1,1,r+1  n1 0,,0]

      根据矩阵乘法:
      \qquad
      a r i = ∑ k = 0 r − 1 l r k u k i + u r i   ⟹ \qquad a_{ri}=\sum\limits_{k=0}^{r-1}l_{rk}u_{ki}+u_{ri}\ \Longrightarrow ari=k=0r1lrkuki+uri   u r i = a r i − ∑ k = 0 r − 1 l r k u k i u_{ri}=a_{ri}-\sum\limits_{k=0}^{r-1}l_{rk}u_{ki} uri=arik=0r1lrkuki   , i = r , r + 1 , ⋯   , n − 1 \ ,\quad i=r,r+1,\cdots,n-1  ,i=r,r+1,,n1
      \qquad
      如图 5 5 5 所示,可求得:

      U [ r , r : ] = [ u r r , ⋯   , u r , n − 1 ] \qquad\qquad U[r,r:]=[u_{rr},\cdots,u_{r,n-1}] U[r,r:]=[urr,,ur,n1]
      也就是:
      U [ r , : ] = [ 0 , ⋯ ⋯   , 0 ⏟ 0   ∼   r − 1 , u r r , ⋯   , u r , n − 1 ] \qquad\qquad U[r,:]=[\underbrace{0,\cdots\cdots,0}_{0\ \sim\ r-1},u_{rr},\cdots,u_{r,n-1}] U[r,:]=[0  r1 0,,0,urr,,ur,n1]

    在这里插入图片描述

    图5 L [ r , : ] L[r,:] L[r,:] 可求出 U [ r , : ] U[r,:] U[r,:]

    \qquad
    3. 第 r r r 步(第二阶段)
     
     由于 U [ r , : r ] = [ u 0 r , ⋯   , u r r ] U[r,:r]=[u_{0r},\cdots,u_{rr}] U[r,:r]=[u0r,,urr] 已求出,也就是:
     
    U [ : , r ] = [ u 0 r , ⋯   , u r r , 0 , ⋯ ⋯   , 0 ⏟ r + 1   ∼   n − 1 ] \qquad\qquad U[:,r]=[u_{0r},\cdots,u_{rr},\underbrace{0,\cdots\cdots,0}_{r+1\ \sim\ n-1}] U[:,r]=[u0r,,urr,r+1  n1 0,,0]
     根据
    a i r = ∑ k = 0 r − 1 l i k u k r + l i r u r r   ⟹ \qquad a_{ir}=\sum\limits_{k=0}^{r-1}l_{ik}u_{kr}+l_{ir}u_{rr}\ \Longrightarrow air=k=0r1likukr+lirurr   l i r = a i r − ∑ k = 0 r − 1 l i k u k r u r r l_{ir}=\dfrac{a_{ir}-\sum\limits_{k=0}^{r-1}l_{ik}u_{kr}}{u_{rr}} lir=urrairk=0r1likukr   , i = r + 1 , ⋯   , n − 1 \ ,\quad i=r+1,\cdots,n-1  ,i=r+1,,n1
    \qquad
     如图 6 6 6 所示,可求得:
     
    L [ r + 1 : , r ] = [ l r + 1 , r , ⋯   , l n − 1 , r ] \qquad\qquad L[r+1:,r]=[l_{r+1,r},\cdots,l_{n-1,r}] L[r+1:,r]=[lr+1,r,,ln1,r]
     也就是:
    L [ : , r ] = [ 0 , ⋯ ⋯   , 0 ⏟ 0   ∼   r − 1 , 1 , l r + 1 , r , ⋯   , l n − 1 , r ] \qquad\qquad L[:,r]=[\underbrace{0,\cdots\cdots,0}_{0\ \sim\ r-1},1,l_{r+1,r},\cdots,l_{n-1,r}] L[:,r]=[0  r1 0,,0,1,lr+1,r,,ln1,r]

    在这里插入图片描述

    图6 U [ : , r ] U[:,r] U[:,r] 可求出 L [ : , r ] L[:,r] L[:,r]

    \qquad

    2.3 LU分解的结果表示

    \qquad 系数矩阵 A A A 分解为 A = L U A=LU A=LU,其中 L L L 为下三角矩阵, U U U 为上三角矩阵:

    [ a 00 a 01 ⋯ a 0 , n − 1 a 10 a 11 ⋯ a 1 , n − 1 ⋮ ⋮ ⋱ ⋮ a n − 1 , 0 a n − 1 , 1 ⋯ a n − 1 , n − 1 ] ⏟ A = [ 1 l 10 1 ⋮ ⋮ ⋱ l n − 1 , 0 l n − 1 , 1 ⋯ 1 ] ⏟ L [ u 00 u 01 ⋯ u 0 , n − 1 u 11 ⋯ u 1 , n − 1 ⋱ ⋮ u n − 1 , n − 1 ] ⏟ U \qquad\qquad\underbrace{\begin{bmatrix} a_{\scriptscriptstyle 00}&a_{\scriptscriptstyle 01}&\cdots&a_{\scriptscriptstyle 0,n-1} \\ a_{\scriptscriptstyle 10}&a_{\scriptscriptstyle 11}&\cdots&a_{\scriptscriptstyle 1,n-1} \\ \vdots&\vdots&\ddots&\vdots\\ a_{\scriptscriptstyle n-1,0}&a_{\scriptscriptstyle n-1,1}&\cdots&a_{\scriptscriptstyle n-1,n-1} \end{bmatrix}}_A =\underbrace{\begin{bmatrix} 1 \\ l_{\scriptscriptstyle 10}&1 \\ \vdots&\vdots&\ddots&\\ l_{\scriptscriptstyle n-1,0}&l_{\scriptscriptstyle n-1,1}&\cdots&1 \end{bmatrix}}_L \underbrace{\begin{bmatrix}u_{\scriptscriptstyle 00}&u_{\scriptscriptstyle 01}&\cdots&u_{\scriptscriptstyle 0,n-1} \\ &u_{\scriptscriptstyle 11}&\cdots&u_{\scriptscriptstyle 1,n-1} \\ & &\ddots&\vdots\\ & & &u_{\scriptscriptstyle n-1,n-1} \end{bmatrix}}_U A a00a10an1,0a01a11an1,1a0,n1a1,n1an1,n1=L 1l10ln1,01ln1,11U u00u01u11u0,n1u1,n1un1,n1

    \qquad 为了节省内存,可以将下三角矩阵 L L L 和上三角矩阵 U U U 合并保存(略去了 L L L 的对角线元素):

    [ 1 l 10 1 ⋮ ⋮ ⋱ l n − 1 , 0 l n − 1 , 1 ⋯ 1 ] ⏟ L [ u 00 u 01 ⋯ u 0 , n − 1 u 11 ⋯ u 1 , n − 1 ⋱ ⋮ u n − 1 , n − 1 ] ⏟ U ∼ [ u 00 u 01 ⋯ u 0 , n − 1 l n − 1 , 0 u 11 ⋯ u 1 , n − 1 ⋮ ⋮ ⋱ ⋮ l n − 1 , 0 l n − 1 , 1 ⋯ u n − 1 , n − 1 ] ⏟ L U \qquad\qquad\underbrace{\begin{bmatrix} 1 \\ l_{\scriptscriptstyle 10}&1 \\ \vdots&\vdots&\ddots&\\ l_{\scriptscriptstyle n-1,0}&l_{\scriptscriptstyle n-1,1}&\cdots&1 \end{bmatrix}}_L \underbrace{\begin{bmatrix}u_{\scriptscriptstyle 00}&u_{\scriptscriptstyle 01}&\cdots&u_{\scriptscriptstyle 0,n-1} \\ &u_{\scriptscriptstyle 11}&\cdots&u_{\scriptscriptstyle 1,n-1} \\ & &\ddots&\vdots\\ & & &u_{\scriptscriptstyle n-1,n-1} \end{bmatrix}}_U\sim\underbrace{\begin{bmatrix}\textcolor{red}{u_{\scriptscriptstyle 00}}&\textcolor{red}{u_{\scriptscriptstyle 01}}&\cdots&\textcolor{red}{u_{\scriptscriptstyle 0,n-1}} \\ \textcolor{blue}{l_{\scriptscriptstyle n-1,0}}&\textcolor{red}{u_{\scriptscriptstyle 11}}&\cdots&\textcolor{red}{u_{\scriptscriptstyle 1,n-1}} \\ \vdots&\vdots&\ddots&\vdots\\ \textcolor{blue}{l_{\scriptscriptstyle n-1,0}}&\textcolor{blue}{l_{\scriptscriptstyle n-1,1}}&\cdots&\textcolor{red}{u_{\scriptscriptstyle n-1,n-1}} \end{bmatrix}}_{LU} L 1l10ln1,01ln1,11U u00u01u11u0,n1u1,n1un1,n1LU u00ln1,0ln1,0u01u11ln1,1u0,n1u1,n1un1,n1

    L L L U U U 合并保存也是为了实现“选主元的三LU分解”

    \qquad

    实现代码

    import numpy as np
    def LUFact(A):# LU Factorization
        n = len(A)
        L = np.eye(n)
        U = np.zeros(A.shape)
        U[0,:] = A[0,:]
        L[1:,0] = A[1:,0]/U[0,0]
        for r in range(1,n):
            lt = L[r,:r].reshape(r,1)
            U[r,r:] = A[r,r:] - np.sum(lt*U[:r,r:],axis=0)
            print('\n#U%d:'%(r))
            print(U)
            if r==n-1:
                continue
            ut = U[:r,r].reshape(r,1)
            L[r+1:,r] = (A[r+1:,r] - np.sum(ut*L[r+1:,:r].T,axis=0))/U[r,r]   
            print('\n#L%d:'%(r))
            print(L)
        return L,U
    
    def Doolittle(A):# LU Factorization---Doolittle
        n = len(A)
        LU = A.copy()
        LU[1:,0] = LU[1:,0]/LU[0,0]
        for r in range(1,n):
            lt = LU[r,:r].reshape(r,1)
            LU[r,r:] = LU[r,r:] - np.sum(lt*LU[:r,r:],axis=0)
            print('\n#U%d:'%(r))
            print(LU)
            if r==n-1:
                continue
            ut = LU[:r,r].reshape(r,1)
            LU[r+1:,r] = (LU[r+1:,r] - np.sum(ut*LU[r+1:,:r].T,axis=0))/LU[r,r] 
            print('\n#L%d:'%(r))   
            print(LU)        
        U = np.triu(LU)
        L = np.tril(LU) - np.diag(np.diag(LU)) + np.eye(n)
        return LU,L,U
        
    def solveLineq(L,U,b):
    	(见第1节,略)  
        return x
        
    if __name__ == '__main__':
        A = np.array( [[  2.,   5.,   7.,   9.,   1.],
                       [  4.,  18.,  20.,  23.,   9.],
                       [  6.,  87.,  76.,  80.,  75.],
                       [  8.,  60.,  64., 112.,  95.],
                       [  2.,  29.,  32.,  89.,  97.]],dtype='float')
        b = np.array([[ 74., 237., 1103., 1243., 997.]],dtype='float').T               
        La,Ua = LUFact(A)
        LU,La,Ua = Doolittle(A)  
        x = solveLineq(La,Ua,b)
        print(x)	# x = [1. 2. 3. 4. 5.]
    

    运行结果:
    在这里插入图片描述

    3. 选主元的LU分解

    3.1 基本原理

    \qquad 高斯消去法在消去过程中一旦遇到主元为 0 0 0,消去过程就会无法继续。与之类似,LU分解的过程中有着同样的问题。图 6 6 6 中所描述的由 U [ : , r ] U[:,r] U[:,r] L [ : , r ] L[:,r] L[:,r] 的过程中,包含了除法运算

    l i r = a i r − ∑ k = 0 r − 1 l i k u k r u r r   , i = r + 1 , ⋯   , n − 1 \qquad\qquad l_{ir}=\dfrac{a_{ir}-\sum\limits_{k=0}^{r-1}l_{ik}u_{kr}}{\textcolor{red}{u_{rr}}}\ ,\quad i=r+1,\cdots,n-1 lir=urrairk=0r1likukr ,i=r+1,,n1

    \qquad 一旦出现了 u r r = 0 u_{rr}=0 urr=0 的情况,那么LU分解也无法继续下去。为了避免这种情况,同样需要采取“选主元”的策略。

    \qquad 由于在LU分解过程中,更新一次 L L L U U U 矩阵,实际包含了 2 2 2 个阶段, 2 2 2 个阶段之间又通过 u r r u_{rr} urr 这个元素值联系在一起:

    ( 1 ) \qquad(1) (1) 第一阶段如图 5 5 5 所示,先求出 U [ r , r : ] = [ u r r , ⋯   , u r , n − 1 ] U[r,r:]=[u_{rr},\cdots,u_{r,n-1}] U[r,r:]=[urr,,ur,n1],其中的 u r r u_{rr} urr

    ( 2 ) \qquad(2) (2) 第二阶段如图 6 6 6 所示, u r r u_{rr} urr 被用于更新 L [ r + 1 : , r ] = [ l r + 1 , r , ⋯   , l i r , ⋯   , l n − 1 , r ] L[r+1:,r]=[l_{r+1,r},\cdots,l_{ir},\cdots,l_{n-1,r}] L[r+1:,r]=[lr+1,r,,lir,,ln1,r]
    \qquad   为了能够计算出 L [ r + 1 : , r ] L[r+1:,r] L[r+1:,r] 中的 l i r l_{ir} lir 元素,必须要求 u r r ≠ 0 u_{rr}\neq 0 urr=0

    ( 3 ) \qquad(3) (3) 为了满足 u r r ≠ 0 u_{rr}\neq 0 urr=0 的要求,可以试着计算 u r r , u r + 1 , r , ⋯   , u n − 1 , r u_{rr},u_{r+1,r},\cdots,u_{n-1,r} urr,ur+1,r,,un1,r,选择其中的最大值作为 u r r u_{rr} urr(通过矩阵的行交换),从而可以避免出现 u r r = 0 u_{rr}=0 urr=0 的情况。

    \qquad

    3.2 算法过程描述

    \qquad 为了实现选主元策略,更适合采用“将 L , U L,U L,U 矩阵显示在一起”的方式来处理。

    \qquad 假设第 r − 1 r-1 r1 步已经完成,此时 L U LU LU 矩阵的状态为:
    \qquad
    [ u 00 u 01 ⋯ u 0 , r − 1 u 0 r u 0 , r + 1 ⋯ ⋯ ⋯ u 0 , n − 1 l 10 u 11 ⋯ u 1 , r − 1 u 1 r u 1 , r + 1 ⋯ ⋯ ⋯ u 1 , n − 1 ⋮ ⋮ ⋱ ⋮ ⋮ ⋮ ⋮ l r − 1 , 0 l r − 1 , 1 ⋯ u r − 1 , r − 1 u r − 1 , r u r − 1 , r + 1 ⋯ ⋯ ⋯ u r − 1 , n − 1 l r 0 l r 1 ⋯ l r , r − 1 a r r a r , r + 1 ⋯ ⋯ ⋯ a r , n − 1 l r + 1 , 0 l r + 1 , 1 ⋯ l r + 1 , r − 1 a r + 1 , r a r + 1 , r + 1 ⋯ ⋯ ⋯ a r + 1 , n − 1 ⋮ ⋮ ⋮ ⋮ ⋮ ⋱ ⋮ l i r , 0 l i r , 1 ⋯ l i r , r − 1 a i r , r a i r , r + 1 ⋱ a i r , n − 1 ⋮ ⋮ ⋮ ⋮ ⋮ ⋱ ⋮ l n − 1 , 0 l n − 1 , 1 ⋯ l n − 1 , r − 1 a n − 1 , r a n − 1 , r + 1 ⋯ ⋯ ⋯ a n − 1 , n − 1 ] ⏟ L U \qquad\qquad\underbrace{\begin{bmatrix}u_{\scriptscriptstyle 00}&u_{\scriptscriptstyle 01}&\cdots&u_{\scriptscriptstyle 0,r-1}&u_{\scriptscriptstyle 0r}&u_{\scriptscriptstyle 0,r+1}&\cdots&\cdots&\cdots&u_{\scriptscriptstyle 0,n-1} \\ l_{\scriptscriptstyle 10}&u_{\scriptscriptstyle 11}&\cdots&u_{\scriptscriptstyle 1,r-1}&u_{\scriptscriptstyle 1r}&u_{\scriptscriptstyle 1,r+1}&\cdots&\cdots&\cdots&u_{\scriptscriptstyle 1,n-1} \\ \vdots&\vdots&\ddots&\vdots&\vdots&\vdots&&&&\vdots \\ l_{\scriptscriptstyle r-1,0}&l_{\scriptscriptstyle r-1,1}&\cdots&u_{\scriptscriptstyle r-1,r-1}&u_{\scriptscriptstyle r-1,r}&u_{\scriptscriptstyle r-1,r+1}&\cdots&\cdots&\cdots&u_{\scriptscriptstyle r-1,n-1}\\ \textcolor{Magenta}{l_{\scriptscriptstyle r0}}&\textcolor{Magenta}{l_{\scriptscriptstyle r1}}&\cdots&\textcolor{Magenta}{l_{\scriptscriptstyle r,r-1}}&\textcolor{red}{a_{\scriptscriptstyle rr}}&\textcolor{red}{a_{\scriptscriptstyle r,r+1}}&\cdots&\cdots&\cdots&\textcolor{red}{a_{\scriptscriptstyle r,n-1}}\\ \textcolor{LimeGreen}{l_{\scriptscriptstyle r+1,0}}&\textcolor{LimeGreen}{l_{\scriptscriptstyle r+1,1}}&\cdots&\textcolor{LimeGreen}{l_{\scriptscriptstyle r+1,r-1}}&\textcolor{red}{a_{\scriptscriptstyle r+1,r}}&\textcolor{red}{a_{\scriptscriptstyle r+1,r+1}}&\cdots&\cdots&\cdots&\textcolor{red}{a_{\scriptscriptstyle r+1,n-1}}\\ \vdots&\vdots&&\vdots&\vdots&\vdots&\ddots&&&\vdots \\ \textcolor{MediumVioletRed}{l_{\scriptscriptstyle i_r,0}}&\textcolor{MediumVioletRed}{l_{\scriptscriptstyle i_r,1}}&\cdots&\textcolor{MediumVioletRed}{l_{\scriptscriptstyle i_r,r-1}}&\textcolor{red}{a_{\scriptscriptstyle i_r,r}}&\textcolor{red}{a_{\scriptscriptstyle i_r,r+1}}& &\ddots&&\textcolor{red}{a_{\scriptscriptstyle i_r,n-1}}\\ \vdots&\vdots&&\vdots&\vdots&\vdots&&&\ddots&\vdots \\ \textcolor{LimeGreen}{l_{\scriptscriptstyle n-1,0}}&\textcolor{LimeGreen}{l_{\scriptscriptstyle n-1,1}}&\cdots&\textcolor{LimeGreen}{l_{\scriptscriptstyle n-1,r-1}}&\textcolor{red}{a_{\scriptscriptstyle n-1,r}}&\textcolor{red}{a_{\scriptscriptstyle n-1,r+1}}&\cdots&\cdots&\cdots&\textcolor{red}{a_{\scriptscriptstyle n-1,n-1}} \end{bmatrix}}_{LU} LU u00l10lr1,0lr0lr+1,0lir,0ln1,0u01u11lr1,1lr1lr+1,1lir,1ln1,1u0,r1u1,r1ur1,r1lr,r1lr+1,r1lir,r1ln1,r1u0ru1rur1,rarrar+1,rair,ran1,ru0,r+1u1,r+1ur1,r+1ar,r+1ar+1,r+1air,r+1an1,r+1u0,n1u1,n1ur1,n1ar,n1ar+1,n1air,n1an1,n1

    \qquad s i = a i r − ∑ k = 0 r − 1 l i k u k r s_i=a_{ir}-\sum\limits_{k=0}^{r-1}l_{ik}u_{kr} si=airk=0r1likukr ,   i = r , r + 1 , ⋯   , n − 1 ,\ i=r,r+1,\cdots,n-1 , i=r,r+1,,n1

    ( 1 ) (1)\quad (1) 由于 a r r = ∑ k = 0 r − 1 l r k u k r + u r r a_{rr}=\sum\limits_{k=0}^{r-1}l_{rk}u_{kr}+u_{rr} arr=k=0r1lrkukr+urr,所以 u r r = a r r − ∑ k = 0 r − 1 l r k u k r u_{rr}=a_{rr}-\sum\limits_{k=0}^{r-1}l_{rk}u_{kr} urr=arrk=0r1lrkukr,正好等于 s r s_r sr 的值。

    \qquad  如果此时的 u r r = 0 u_{rr}=0 urr=0,随后的求取 L L L 矩阵元素的过程就无法继续。

    \qquad
    ( 2 ) (2)\quad (2) 假设第 r r r 步正好出现了 u r r = 0 u_{rr}=0 urr=0 的情况。

    \qquad  此时必须将 r r r 行的前 r r r 个元素 L [ r , : r ] = [ l r , 0 , ⋯   , l r , r − 1 ] L[r,:r]=[l_{r,0},\cdots,l_{r,r-1}] L[r,:r]=[lr,0,,lr,r1],与 i r i_r ir 行的前 r r r 个元素 L [ i r , : r ] = [ l i r , 0 , ⋯   , l i r , r − 1 ] L[i_r,:r]=[l_{i_r,0},\cdots,l_{i_r,r-1}] L[ir,:r]=[lir,0,,lir,r1] , r + 1 ≤ i r ≤ n − 1 ,r+1\le i_r\le n-1 ,r+1irn1 进行交换,才能改变第 r r r 步中 u r r = 0 u_{rr}=0 urr=0 的情况。

    \qquad  因此,公式 s i = a i r − ∑ k = 0 r − 1 l i k u k r s_i=a_{ir}-\sum\limits_{k=0}^{r-1}l_{ik}u_{kr} si=airk=0r1likukr ,   i = r + 1 , ⋯   , n − 1 ,\ i=r+1,\cdots,n-1 , i=r+1,,n1 就是把第 r r r 行下面的某一行 i i i 交换到第 r r r 行时、所求的 u r r u_{rr} urr 的值。

    ( 3 ) (3)\quad (3) 为了选择一个最靠谱的行换到第 r r r 行,仿照列主元消去法的思路,选择具有最大 ∣ u r r ∣ \vert u_{rr}\vert urr 的值的行进行行交换。

    \qquad  也就是选择 ∣ s i r ∣ = max ⁡ r + 1 ≤ i ≤ n − 1 ∣ s i ∣ \vert s_{i_r}\vert=\displaystyle\max_{\scriptscriptstyle r+1\le i\le n-1}\vert s_i\vert sir=r+1in1maxsi,使用第 i r i_r ir 行的前 r r r 个元素 L [ i r , : r ] = [ l i r , 0 , ⋯   , l i r , r − 1 ] L[i_r,:r]=[l_{i_r,0},\cdots,l_{i_r,r-1}] L[ir,:r]=[lir,0,,lir,r1] 算出来的 ∣ s i r ∣ \vert s_{i_r}\vert sir 值最大。这样就可以避免第 r r r 步中 u r r = 0 u_{rr}=0 urr=0 的情况。

    [ u 00 u 01 ⋯ u 0 , r − 1 u 0 r u 0 , r + 1 ⋯ u 0 , i ⋯ u 0 , n − 1 l 10 u 11 ⋯ u 1 , r − 1 u 1 r u 1 , r + 1 ⋯ u 1 , i ⋯ u 1 , n − 1 ⋮ ⋮ ⋱ ⋮ ⋮ ⋮ ⋮ ⋮ l r − 1 , 0 l r − 1 , 1 ⋯ u r − 1 , r − 1 u r − 1 , r u r − 1 , r + 1 ⋯ u r − 1 , i ⋯ u r − 1 , n − 1 l i r , 0 l i r , 1 ⋯ l i r , r − 1 a i r , r a i r , r + 1 ⋯ a i r , i ⋯ a i r , n − 1 l r + 1 , 0 l r + 1 , 1 ⋯ l r + 1 , r − 1 a r + 1 , r a r + 1 , r + 1 ⋯ a r + 1 , i ⋯ a r + 1 , n − 1 ⋮ ⋮ ⋮ ⋮ ⋮ ⋱ ⋮ l r 0 l r 1 ⋯ l r , r − 1 a r r a r , r + 1 ⋱ a r , n − 1 ⋮ ⋮ ⋮ ⋮ ⋮ ⋱ ⋮ l n − 1 , 0 l n − 1 , 1 ⋯ l n − 1 , r − 1 a n − 1 , r a n − 1 , r + 1 ⋯ a n − 1 , i ⋯ a n − 1 , n − 1 ] ⏟ L U \qquad\qquad\underbrace{\begin{bmatrix}u_{\scriptscriptstyle 00}&u_{\scriptscriptstyle 01}&\cdots&u_{\scriptscriptstyle 0,r-1}&u_{\scriptscriptstyle 0r}&u_{\scriptscriptstyle 0,r+1}&\cdots&\textcolor{blue}{u_{\scriptscriptstyle 0,i}}&\cdots&u_{\scriptscriptstyle 0,n-1} \\ l_{\scriptscriptstyle 10}&u_{\scriptscriptstyle 11}&\cdots&u_{\scriptscriptstyle 1,r-1}&u_{\scriptscriptstyle 1r}&u_{\scriptscriptstyle 1,r+1}&\cdots&\textcolor{blue}{u_{\scriptscriptstyle 1,i}}&\cdots&u_{\scriptscriptstyle 1,n-1} \\ \vdots&\vdots&\ddots&\vdots&\vdots&\vdots&&\vdots&&\vdots \\ l_{\scriptscriptstyle r-1,0}&l_{\scriptscriptstyle r-1,1}&\cdots&u_{\scriptscriptstyle r-1,r-1}&u_{\scriptscriptstyle r-1,r}&u_{\scriptscriptstyle r-1,r+1}&\cdots&\textcolor{blue}{u_{\scriptscriptstyle r-1,i}}&\cdots&u_{\scriptscriptstyle r-1,n-1}\\ \textcolor{MediumVioletRed}{l_{\scriptscriptstyle i_r,0}}&\textcolor{MediumVioletRed}{l_{\scriptscriptstyle i_r,1}}&\cdots&\textcolor{MediumVioletRed}{l_{\scriptscriptstyle i_r,r-1}}&\textcolor{red}{a_{\scriptscriptstyle i_r,r}}&\textcolor{red}{a_{\scriptscriptstyle i_r,r+1}}&\cdots&\textcolor{red}{a_{\scriptscriptstyle i_r,i}}&\cdots&\textcolor{red}{a_{\scriptscriptstyle i_r,n-1}}\\ \textcolor{LimeGreen}{l_{\scriptscriptstyle r+1,0}}&\textcolor{LimeGreen}{l_{\scriptscriptstyle r+1,1}}&\cdots&\textcolor{LimeGreen}{l_{\scriptscriptstyle r+1,r-1}}&\textcolor{red}{a_{\scriptscriptstyle r+1,r}}&\textcolor{red}{a_{\scriptscriptstyle r+1,r+1}}&\cdots&\textcolor{red}{a_{\scriptscriptstyle r+1,i}}&\cdots&\textcolor{red}{a_{\scriptscriptstyle r+1,n-1}}\\ \vdots&\vdots&&\vdots&\vdots&\vdots&\ddots&&&\vdots \\ \textcolor{Magenta}{l_{\scriptscriptstyle r0}}&\textcolor{Magenta}{l_{\scriptscriptstyle r1}}&\cdots&\textcolor{Magenta}{l_{\scriptscriptstyle r,r-1}}&\textcolor{red}{a_{\scriptscriptstyle rr}}&\textcolor{red}{a_{\scriptscriptstyle r,r+1}}& &\ddots&&\textcolor{red}{a_{\scriptscriptstyle r,n-1}}\\ \vdots&\vdots&&\vdots&\vdots&\vdots&&&\ddots&\vdots \\ \textcolor{LimeGreen}{l_{\scriptscriptstyle n-1,0}}&\textcolor{LimeGreen}{l_{\scriptscriptstyle n-1,1}}&\cdots&\textcolor{LimeGreen}{l_{\scriptscriptstyle n-1,r-1}}&\textcolor{red}{a_{\scriptscriptstyle n-1,r}}&\textcolor{red}{a_{\scriptscriptstyle n-1,r+1}}&\cdots&\textcolor{red}{a_{\scriptscriptstyle n-1,i}}&\cdots&\textcolor{red}{a_{\scriptscriptstyle n-1,n-1}} \end{bmatrix}}_{LU} LU u00l10lr1,0lir,0lr+1,0lr0ln1,0u01u11lr1,1lir,1lr+1,1lr1ln1,1u0,r1u1,r1ur1,r1lir,r1lr+1,r1lr,r1ln1,r1u0ru1rur1,rair,rar+1,rarran1,ru0,r+1u1,r+1ur1,r+1air,r+1ar+1,r+1ar,r+1an1,r+1u0,iu1,iur1,iair,iar+1,ian1,iu0,n1u1,n1ur1,n1air,n1ar+1,n1ar,n1an1,n1

    \qquad
    ( 4 ) (4)\quad (4) 矩阵 L U LU LU 经过行交换之后,会导致 A A A 矩阵中的行交换。

    \qquad   L [ i r , : ] = [ l i r , 0 , ⋯   , l i r , r − 1 , 1 , 0 , ⋯   , 0 ⏟ r + 1 ∼ n − 1 ] L[i_r,:]=[l_{i_r,0},\cdots,l_{i_r,r-1},1,\underbrace{0,\cdots,0}_{r+1\sim n-1}] L[ir,:]=[lir,0,,lir,r1,1,r+1n1 0,,0] 与任何一列 U [ : r , i ] = [ u 0 i , ⋯   , u r − 1 , i , ∗ , ⋯   , 0 ⏟ r ∼ n − 1 ] U[:r,i]=[u_{0i},\cdots,u_{r-1,i},\underbrace{*,\cdots,0}_{r\sim n-1}] U[:r,i]=[u0i,,ur1,i,rn1 ,,0] 的向量积,即为 A A A 矩阵中第 i r i_r ir A [ i r , : ] A[i_r,:] A[ir,:]

    \qquad   L [ r , : ] = [ l r , 0 , ⋯   , l r , r − 1 , 1 , 0 , ⋯   , 0 ⏟ r + 1 ∼ n − 1 ] L[r,:]=[l_{r,0},\cdots,l_{r,r-1},1,\underbrace{0,\cdots,0}_{r+1\sim n-1}] L[r,:]=[lr,0,,lr,r1,1,r+1n1 0,,0] 与任何一列 U [ : r , i ] = [ u 0 i , ⋯   , u r − 1 , i , ∗ , ⋯   , 0 ⏟ r ∼ n − 1 ] U[:r,i]=[u_{0i},\cdots,u_{r-1,i},\underbrace{*,\cdots,0}_{r\sim n-1}] U[:r,i]=[u0i,,ur1,i,rn1 ,,0] 的向量积,即为 A A A 矩阵中第 r r r A [ r , : ] A[r,:] A[r,:]

    由于到第 r − 1 r-1 r1 步为止,右下角红色位置 a i j a_{ij} aij 对应的L矩阵(除了对角线元素为 1 1 1之外)和U矩阵的元素值都为 0 0 0,因此,并不会影响此部分区域之外的 a i j a_{ij} aij 元素的值。

    \qquad  该过程相当于 P A = L U PA=LU PA=LU P = [ 1 ⋱ ⋱ 1 ⋱ 1 ⋱ ⋱ 1 ] i i r P=\begin{bmatrix}1&&&\\&\ddots&&\\&&\ddots&&1&\\&&&\ddots\\&&1&&\ddots\\&&&&&\ddots&&\\&&&&&&1\end{bmatrix}\begin{matrix}i\\\\i_r\\\end{matrix} P=1111iir
    \qquad
    \qquad
    \qquad 因此,经过选主元的操作之后,LU分解的过程仍旧采用1.2节中描述的过程。如果在整个过程中发生了 m m m 次行交换,那么最终的LU分解过程实际上是:

    P m ⋯ P 2 P 1 A = L U \qquad\qquad\qquad P_m\cdots P_2P_1A=LU PmP2P1A=LU

    \qquad

    实现代码

    import numpy as np
    def Doolittle_pivot(A):  
        n = len(A)
        LU = A.copy()
        #LU[1:,0] = A[1:,0]/LU[0,0]
        order1 = np.arange(n)
        for r in range(n):        
            ut = LU[:r,r].reshape(r,1)
            si = A[r:,r] - np.sum(ut*LU[r:,:r].T,axis=0)
            ir = np.argmax(np.abs(si))
            if ir!=0:
                LU[[r,r+ir],:] = LU[[r+ir,r],:]
                order1[[r,r+ir]] = order1[[r+ir,r]]           
            lt = LU[r,:r].reshape(r,1)
            LU[r,r:] = LU[r,r:] - np.sum(lt*LU[:r,r:],axis=0)
            print('\n#U%d:'%(r))
            print(LU)        
            if r==n-1:
                continue
            # ut = LU[:r,r].reshape(r,1)
            LU[r+1:,r] = (LU[r+1:,r] - np.sum(ut*LU[r+1:,:r].T,axis=0))/LU[r,r] 
            print('\n#L%d:'%(r))   
            print(LU)           
        U = np.triu(LU)
        L = np.tril(LU) - np.diag(np.diag(LU)) + np.eye(n)
        order0 = []
        [order0.insert(i,np.where(order1==i)[0][0]) for i in range(n)]        
        return LU,L,U,np.array(order0),order1
    
    def solveLineq(L,U,b):
    	(见第1节,略)  
        return x
        
    if __name__ == '__main__':
        A = np.array( [[  2.,   5.,   7.,   9.,   1.],
                       [  4.,  18.,  20.,  23.,   9.],
                       [  6.,  87.,  76.,  80.,  75.],
                       [  8.,  60.,  64., 112.,  95.],
                       [  2.,  29.,  32.,  89.,  97.]],dtype='float')
        b = np.array([[ 74., 237., 1103., 1243., 997.]],dtype='float').T
        LU,La,Ua,order0,order1 = Doolittle_pivot(A)
        print(La)
        print(Ua)
        print('L*U:\n',np.dot(La,Ua)[order0,:])
        x = solveLineq(La,Ua,b[order1,:])      # b行变换:Pb
        print(x)	# x = [1. 2. 3. 4. 5.]
    

    运行结果:
    在这里插入图片描述

    4. 对称正定矩阵的Cholesky分解

    4.1 基本原理

    • A A A对称矩阵,在进行 A = L U A=LU A=LU 分解后,将 U U U 矩阵写为 U = D U 0 U=DU_0 U=DU0,也就是:

    [ u 00 u 01 ⋯ u 0 , n − 1 u 11 ⋯ u 1 , n − 1 ⋱ ⋮ u n − 1 , n − 1 ] ⏟ U = [ u 00 u 11 ⋱ u n − 1 , n − 1 ] ⏟ D [ 1 u 01 u 00 ⋯ ⋯ u 0 , n − 1 u 00 1 u 12 u 11 ⋯ u 1 , n − 1 u 11 ⋱ ⋮ ⋱ ⋮ 1 ] ⏟ U 0 \qquad\qquad \underbrace{\begin{bmatrix}u_{\scriptscriptstyle 00}&u_{\scriptscriptstyle 01}&\cdots&u_{\scriptscriptstyle 0,n-1} \\ &u_{\scriptscriptstyle 11}&\cdots&u_{\scriptscriptstyle 1,n-1} \\ & &\ddots&\vdots\\ & & &u_{\scriptscriptstyle n-1,n-1} \end{bmatrix}}_U=\underbrace{\begin{bmatrix}u_{\scriptscriptstyle 00} \\ &u_{\scriptscriptstyle 11} \\ & &\ddots\\ & & &u_{\scriptscriptstyle n-1,n-1} \end{bmatrix}}_D \underbrace{\begin{bmatrix}1&\frac{u_{\scriptscriptstyle 01}}{u_{\scriptscriptstyle 00}}&\cdots&\cdots&\frac{u_{\scriptscriptstyle 0,n-1}}{u_{\scriptscriptstyle 00}} \\ &1&\frac{u_{\scriptscriptstyle 12}}{u_{\scriptscriptstyle 11}}&\cdots&\frac{u_{\scriptscriptstyle 1,n-1}}{u_{\scriptscriptstyle 11}} \\ & &\ddots&&\vdots\\ &&&\ddots&\vdots\\ & &&&1 \end{bmatrix}}_{U_0} U u00u01u11u0,n1u1,n1un1,n1=D u00u11un1,n1U0 1u00u011u11u12u00u0,n1u11u1,n11

    \qquad 其中, D D D 是对角阵, U 0 U_0 U0 是单位上三角阵(对角线元素值为 1 1 1)。

    \qquad 那么

    { A = L U = L D U 0 A = A T ⟹ A = L D U 0 = ( L D U 0 ) T = U 0 T D L T \qquad\qquad\left\{\begin{aligned}A&=LU=LDU_0 \\A&=A^T \end{aligned}\right.\Longrightarrow A=LDU_0=(LDU_0)^T=U_0^TDL^T {AA=LU=LDU0=ATA=LDU0=(LDU0)T=U0TDLT

    \qquad 由于LU分解的唯一性,可得到: U 0 = L T \quad U_0=L^T U0=LT

    \qquad 因此,对称矩阵 A A A 可以分解为: A = L D L T \quad A=LDL^T A=LDLT
    \qquad

    • 假设对称矩阵 A A A 也是正定的,那么

    [ d 0 d 1 ⋱ d n − 1 ] ⏟ D = [ d 0 1 2 d 1 1 2 ⋱ d n − 1 1 2 ] ] ⏟ D 1 2 [ d 0 1 2 d 1 1 2 ⋱ d n − 1 1 2 ] ⏟ D 1 2 \qquad\qquad\underbrace{\begin{bmatrix}d_{\scriptscriptstyle 0} \\ &d_{\scriptscriptstyle 1} \\ & &\ddots\\ & & &d_{\scriptscriptstyle n-1} \end{bmatrix}}_D= \underbrace{\begin{bmatrix}d_{\scriptscriptstyle 0}^{\frac{1}{2}} \\ &d_{\scriptscriptstyle 1}^{\frac{1}{2}} \\ & &\ddots\\ & & &d_{\scriptscriptstyle n-1}^{\frac{1}{2}} \end{bmatrix}]}_{D^{\frac{1}{2}}} \underbrace{\begin{bmatrix}d_{\scriptscriptstyle 0}^{\frac{1}{2}} \\ &d_{\scriptscriptstyle 1}^{\frac{1}{2}} \\ & &\ddots\\ & & &d_{\scriptscriptstyle n-1}^{\frac{1}{2}} \end{bmatrix}}_{D^{\frac{1}{2}}} D d0d1dn1=D21 d021d121