不完全LU分解

不完全 LU 分解（incomplete LU factorization，ILU）在高斯消去中舍弃部分填充或小元素，构造近似因子 $A\approx LU$。这些因子通常用于预条件：每轮通过两个三角系统近似求解 $Az=r$，而不是把不完全分解误当成原系统的精确分解。舍弃规则决定存储、求解成本与近似质量。

一个被舍弃的填充位置

考虑对称正定矩阵 $$A=\left(\begin{array}{ccc}4& -1& -1\\ -1& 4& 0\\ -1& 0& 4\end{array}\right).$$ 消去第一列时，第二、三行的首乘子都为 $-1/4$，两个剩余对角元变为 $15/4$，而原来为零的 $(2,3)$、$(3,2)$ 位置会出现 $-1/4$。这类新非零元称为填充。

ILU(0) 只保留原矩阵非零位置和对角位置的因子元素，丢掉这些新位置，得到 $$L=\left(\begin{array}{ccc}1& 0& 0\\ -1/4& 1& 0\\ -1/4& 0& 1\end{array}\right),U=\left(\begin{array}{ccc}4& -1& -1\\ 0& 15/4& 0\\ 0& 0& 15/4\end{array}\right).$$ 但乘积为 $$M=LU=\left(\begin{array}{ccc}4& -1& -1\\ -1& 4& 1/4\\ -1& 1/4& 4\end{array}\right)\ne A.$$ “零填充”约束的是因子的存储模式，并不要求乘积 $LU$ 在所有被舍弃的位置也为零。这个区别正是理解 ILU(0) 的关键。

矩阵位置示意。两个非对角空位原本为零，消去第一行后会产生耦合；ILU(0)保留原稀疏模式，换取近似而非精确分解。

保留位置上的分解条件

设预定模式为 $𝒮$，包含所有对角位置。取 $L_{ii}=1$，并在不属于允许模式的位置设相应因子元素为零。对保留的上三角位置，按 $$U_{ij}=A_{ij}-\sum _{k<i}{L}_{ik}{U}_{kj},i\le j,$$ 对保留的下三角位置，按 $$L_{ij}=\frac{A_{ij}-\sum _{k<j}{L}_{ik}{U}_{kj}}{U_{jj}},i>j,$$ 依消去顺序计算，前提是所需主元非零。这些等式保证保留位置上 $(LU{)}_{ij}=A_{ij}$；舍弃位置的误差不受这条等式控制。ILU(k) 根据填充层级保留更多位置，阈值 ILU 则根据数值大小和存储限制丢弃元素，名字相近但规则不同。

通过三角求解应用预条件子

对残差 $r$，先解 $Ly=r$，再解 $Uz=y$，便得 $z=M^{-1}r$。无需形成逆矩阵。若 $A=M-E$，则左预条件矩阵为 $$M^{-1}A=I-M^{-1}E.$$ 这说明近似质量取决于误差经过 $M^{-1}$ 后的作用，而不只是 $A-LU$ 的条目看起来小。对预条件 Richardson 迭代，误差递推为 $e_{new}=M^{-1}Ee$；若在某个一致矩阵范数下 $\Vert M^{-1}E\Vert <1$，重复取界即可证明几何收敛。这个充分条件不是所有 ILU 必然满足的保证。

主例取精确解为全一向量，右端为 $b=(2,3,3{)}^T$。一次近似求解先得 $y=(2,7/2,7/2{)}^T$，再得 $$z=(29/30,14/15,14/15{)}^T.$$ 直接代回得到真实残差 $b-Az=(0,7/30,7/30{)}^T$，因此不能报告已经精确解出系统。有解变式将矩阵改成三对角 $\left(\begin{array}{ccc}4& -1& 0\\ -1& 4& -1\\ 0& -1& 4\end{array}\right)$；无主元消去不产生模式外填充，ILU(0)便与精确 LU 相同。

breakdown、排序与算法选择

即使原矩阵非奇异，舍弃后的某个主元仍可能为零，算法因而失败。重排序、对角扰动和允许更多填充可改善表现，但也改变了近似与成本。对称正定矩阵使用 CG 时，还应保持对称正定的预条件作用；一般带置换的 ILU 不自动满足这个要求，常改用不完全 Cholesky 或选 GMRES、BiCGSTAB 等相容方法。

下面直接实现不带选主元的 ILU(0)，用稠密列表展示模式检查，再提供前后代入。它返回两个新因子，不修改输入；稠密三重循环的最坏时间为 $O(n^3)$、存储 $O(n^2)$。真正的稀疏实现只遍历允许的非零及相应交集，成本取决于模式，不能由本代码声称线性时间。

Python 3

import math


def ilu0(A):
    W = [list(map(float, row)) for row in A]; n = len(W)
    if not n or any(len(row) != n for row in W) or not all(math.isfinite(x) for row in W for x in row):
        raise ValueError("finite nonempty square matrix required")
    pattern = [[i == j or W[i][j] != 0 for j in range(n)] for i in range(n)]
    for i in range(n):
        for j in range(i):
            if pattern[i][j]:
                if W[j][j] == 0:
                    raise ArithmeticError("zero incomplete pivot")
                W[i][j] /= W[j][j]
                for k in range(j+1, n):
                    if pattern[i][k]:
                        W[i][k] -= W[i][j]*W[j][k]
        if W[i][i] == 0 or not all(math.isfinite(x) for x in W[i]):
            raise ArithmeticError("zero pivot or nonfinite factor")
    L = [[1.0 if i == j else W[i][j] if i > j else 0.0 for j in range(n)] for i in range(n)]
    U = [[W[i][j] if i <= j else 0.0 for j in range(n)] for i in range(n)]
    return L, U


def apply_ilu(L, U, rhs):
    n = len(rhs); y = [0.0]*n; x = [0.0]*n
    for i in range(n):
        y[i] = float(rhs[i])-sum(L[i][j]*y[j] for j in range(i))
    for i in range(n-1, -1, -1):
        x[i] = (y[i]-sum(U[i][j]*x[j] for j in range(i+1, n)))/U[i][i]
    if not all(math.isfinite(v) for v in x):
        raise ArithmeticError("nonfinite triangular solution")
    return x


A = [[4, -1, -1], [-1, 4, 0], [-1, 0, 4]]
L, U = ilu0(A)
print(apply_ilu(L, U, [2, 3, 3]))  # [29/30, 14/15, 14/15] up to rounding

辅助求解函数要求输入确由同一 ilu0 返回的因子，以及匹配维数的有限右端。若将它封装成面向任意调用者的库接口，还应补上独立的维数、主元和输入类型检查。