跳转至

矩阵分解法

LU分解法

我们在高斯消元法中,对增广矩阵进行初等行变换得到了系数矩阵为上三角阵的线性方程组,由此我们可以通过许多个初等矩阵左乘原方程组增广矩阵得到最终的方程组,为此我们做下面定义:

定义第 \(k\) 步消元所对应的初等矩阵 \(L_k\) 为:

\[ L_k = \begin{pmatrix} 1 & & & & & & \\ & \ddots & & & & & \\ & & 1 & & & & \\ & & -l_{(k+1)k} & 1 & & & \\ & & \vdots & & \ddots & & \\ & & -l_{nk} & & & 1 & \end{pmatrix} \]

我们直接给出其逆矩阵形式,不再赘述:

\[ L_k^{-1} = \begin{pmatrix} 1 & & & & & & \\ & \ddots & & & & & \\ & & 1 & & & & \\ & & l_{(k+1)k} & 1 & & & \\ & & \vdots & & \ddots & & \\ & & l_{nk} & & & 1 & \end{pmatrix} \]

由此我们可以得到: \(L_{n-1} L_{n-2} \cdots L_2 L_1 (A^{(0)} \mid b^{(0)}) = (A^{(n-1)} \mid b^{(n-1)})\)

根据高斯消元法过程我们知道,最后 \(A^{(n-1)}\) 是一个上三角矩阵,我们不妨设 \(A^{(n-1)} = U\)

对方程进行变换得到:\(A =A^{(0)}= L_1^{-1} L_2^{-1} \cdots L_{n-1}^{-1} U\)

\(L_1^{-1} L_2^{-1} \cdots L_{n-1}^{-1} = L\),则可以求得:

\[ L = \begin{pmatrix} 1 & & & & \\ l_{21} & 1 & & & \\ l_{31} & l_{32} & 1 & & \\ \vdots & \vdots & \vdots & \ddots & \\ l_{n1} & l_{n2} & \cdots & l_{n,n-1} & 1 \end{pmatrix} \]

则我们将系数矩阵 \(A\) 分解为一个单位下三角矩阵 \(L\) 和一个上三角矩阵 \(U\) 的乘积,即 \(A = LU\)。即 \(Ax = LUx = b\) 。关于求解过程求解 \(A\) 分解方法和 \(x\) 的过程见下:

求解步骤

求解 \(LU \to\) 求解 \(x\)

  1. 分解过程:对 \(k = 1\)\(n\),计算:
  2. \(U\) 的第 \(k\) 行:\(u_{kj} = a_{kj} - \sum_{r=1}^{k-1} l_{kr} u_{rj} \quad (j = k, \cdots, n)\)
  3. \(L\) 的第 \(k\) 列:\(l_{ik} = (a_{ik} - \sum_{r=1}^{k-1} l_{ir} u_{rk}) / u_{kk} \quad (i = k+1, \cdots, n)\) 计算顺序就是:\(U\) 一行,\(L\) 一列(利用 \(L\) 第一行已知仅首元为1其余皆为0的条件不断地用 \(L\) 第一行乘以 \(U\) 的各列,得到各列首元,即 \(U\) 的第一行,此时 \(U\) 第一列仅首元不为零,做相似操作求出 \(L\) 第一列,不断重复即可求出 \(L,U\)

  4. 求解过程:令 \(Ux=y\) ,这部分相对较为简单,只需要利用系数矩阵是上(或下)三角矩阵求解就行,U的对角线并不一定都是1,需要求解

  5. \(Ly = b\) (前代): \(y_1 = b_1, \quad y_k = b_k - \sum_{j=1}^{k-1} l_{kj} y_j \quad (k=2,\cdots,n)\)
  6. \(Ux = y\) (回代): \(x_n = y_n / u_{nn}, \quad x_k = (y_k - \sum_{j=k+1}^n u_{kj} x_j) / u_{kk} \quad (k=n-1, \cdots, 1)\)

注意保证分解是合理适用的

关于 LU 三角分解的几点说明

  1. 用直接三角分解法求线性方程组 \(Ax=b\) 的乘除法的计算量也是 \(O(n^3)\) 数量级,与高斯消元法解同阶数方程组所需计算量基本相同
  2. 从三角分解的计算公式可以看出,稀疏矩阵\(A\) 的分解与 \(b\) 无关,因此解许多个系数矩阵相同的线性方程组时会更加简便
  3. 完成 \(A=LU\) 分解后,\(\det(A) = \det(L)\times \det(U) = u_{11}u_{22}\cdots u_{nn}\)
  4. 直接三角解法一般地可采用选主元地技术进行方程组求解使算法更加稳定
  5. 另外一种分解方法(Crout分解):令 \(U\) 的对角线元素全是1,\(L\) 全未知

特殊情况(Cholesky分解法):

若系数矩阵 \(A\) 为对称正定矩阵,则存在一个对角元均为正数的下三角矩阵 \(L\in R^{n\times n}\) 使得:

\[ A = LL^T \]

其中 \(L\) 称为\(A\) 的 Cholesky 因子

复习一下什么是正定矩阵吧

一个 \(n \times n\) 的实对称矩阵 \(A\) 是正定的,如果对于任意非零向量 \(x \in \mathbb{R}^n\),都有:

\[ x^\dagger A x > 0 \]

两个判断条件:

  1. 特征值判定法:\(A\) 是正定矩阵 当且仅当 它的所有特征值都是正数: \(\(\lambda_i > 0 \quad (i = 1, 2, \cdots, n)\)\)
  2. 顺序主子式判定法:\(A\) 是正定矩阵 当且仅当 它的所有顺序主子式都大于零: $$ \begin{aligned} &\Delta_1 = a_{11} > 0 \ &\Delta_2 = \begin{vmatrix} a_{11} & a_{12} \ a_{21} & a_{22} \end{vmatrix} > 0 \ &\Delta_3 = \begin{vmatrix} a_{11} & a_{12} & a_{13} \ a_{21} & a_{22} & a_{23} \ a_{31} & a_{32} & a_{33} \end{vmatrix} > 0 \ &\vdots \ &\Delta_n = \det(A) > 0 \end{aligned} $$