文江博客开发文档 fast.ai 数值线性代数讲义中文版 v2 文章详情

文章来源于网络收集而来，版权归原创者所有，如有侵权请及时联系！

五、LU 分解

发布于 2025-01-01 12:38:40 字数 10687 浏览 0 评论 0 收藏 0

fbpca 和我们自己的 randomized_range_finder 方法都使用 LU 分解，它将矩阵分解为下三角矩阵和上三角矩阵的乘积。

高斯消元

本节基于 Trefethen 的 20-22 讲座。

如果你不熟悉高斯消元或需要复习，请观看此可汗学院视频。

让我们手动使用高斯消元来回顾：

$A=\begin{pmatrix} 1 & -2 & -2 & -3 \\ 3 & -9 & 0 & -9 \\ -1 & 2 & 4 & 7 \\ -3 & -6 & 26 & 2 \end{pmatrix}$

答案：

$LU = \begin{bmatrix} 1 & 0 & 0 & 0\\ 3 & 1 & 0 & 0 \\ -1 & 0 & 1 & 0 \\ -3 & 4 & -2 & 1\end{bmatrix} \cdot \begin{bmatrix} 1 & -2 & -2 & -3 \\ 0 & -3 & 6 & 0 \\ 0 & 0 & 2 & 4 \\ 0 & 0 & 0 & 1 \end{bmatrix}$

以上示例来自 Trefethen 的讲座 20,21。

高斯消元通过在左侧应用线性变换，将线性方程组变换为上三角形方程组。它是三角形三角化。

$L_{m-1} \dots L_2 L_1 A = U$

L 是单位下三角形：所有对角线元素都是 1。

def LU(A):
    U = np.copy(A)
    m, n = A.shape
    L = np.eye(n)
    for k in range(n-1):
        for j in range(k+1,n):
            L[j,k] = U[j,k]/U[k,k]
            U[j,k:n] -= L[j,k] * U[k,k:n]
    return L, U

A = np.array([[2,1,1,0],[4,3,3,1],[8,7,9,5],[6,7,9,8]]).astype(np.float)

L, U = LU(A)

np.allclose(A, L @ U)

# True

LU 分解很有用！

求解 Ax = b 变为 LUx = b ：

找到 A = LU
解 Ly = b
解 Ux = y
完事

工作量

高斯消元的工作量： $2\cdot\frac{1}{3} n^3$

内存

在上面，我们创建了两个新的矩阵， L 和 U 。但是，我们可以将 L 和 U 的值存储在矩阵 A 中（覆盖原始矩阵）。由于 L 的对角线都是 1，因此不需要存储。在原地进行因式分解或计算，是数值线性代数中用于节省内存的常用技术。注意：如果你将来需要再次使用原始矩阵 A ，则不希望这样做。其中一个作业问题是重写 LU 方法来原地操作。

考虑矩阵：

$A = \begin{bmatrix} 10^{-20} & 1 \\ 1 & 1 \end{bmatrix}$

A = np.array([[1e-20, 1], [1,1]])

手动使用高斯消元法计算 L 和 U ：

# 练习：

np.set_printoptions(suppress=True)

# 练习：

L2, U2 = LU(A)

'''
[[  1.00000000e-20   1.00000000e+00]
 [  0.00000000e+00  -1.00000000e+20]]
'''

L2, U2

'''
(array([[  1.00000000e+00,   0.00000000e+00],
        [  1.00000000e+20,   1.00000000e+00]]),
 array([[  1.00000000e-20,   1.00000000e+00],
        [  0.00000000e+00,  -1.00000000e+20]]))
'''

np.allclose(L1, L2)

# True

np.allclose(U1, U2)

# True

np.allclose(A, L2 @ U2)

# False

这是使用交换主元进行 LU 分解的动机。

这也说明 LU 分解是稳定的，但不是向后稳定的。（剧透：即使部分交换主元，LU 对某些矩阵来说“爆炸性不稳定”，但在实践中稳定）

稳定性

问题 f 的算法 $\hat f$ 是稳定的，如果对于每个 x ：

$\frac{\lVert \hat{f}(x) - f(y) \rVert}{ \lVert f(y) \rVert } = \mathcal{O}(\varepsilon_{machine})$

对于一些 y ：

$\frac{\lVert y - x \rVert }{\lVert x \rVert} = \mathcal{O}(\varepsilon_{machine})$

一个稳定的算法几乎可以为几乎正确的问题提供正确的答案（Trefethen，第 104 页）。

翻译：

正确的问题： x
几乎正确的问题： y
正确答案： f
几乎正确的问题的正确答案： f(y)

向后稳定

向后稳定性比稳定性更强大，更简单。

问题 f 的算法 $\hat f$ 是向后稳定的，如果对于每个 x ，

$\hat{f}(x) = f(y)$

对于一些 y ：

$\frac{\lVert y - x \rVert }{\lVert x \rVert} = \mathcal{O}(\varepsilon_{machine})$

向后稳定的算法为几乎正确的问题提供了正确的答案（Trefethen，第 104 页）。

翻译：

正确的问题： x
几乎正确的问题： y
正确答案： f
几乎正确的问题的正确答案： f(y)

带有交换主元的 LU 分解

让我们看看矩阵：

$\hat{A} = \begin{bmatrix} 1 & 1 \\ 10^{-20} & 1 \end{bmatrix}$

A = np.array([[1,1], [1e-20, 1]])

手动使用高斯消元法计算 L 和 U ：

$\hat{L} = \begin{bmatrix} 1 & 0 \\ 10^{-20} & 1 \end{bmatrix}$

$\hat{U} = \begin{bmatrix} 1 & 1 \\ 0 & 1 - 10^{-20} \end{bmatrix}$

L, U = LU(A)

np.allclose(A, L @ U)

# True

想法：我们可以切换行的顺序，来获得更稳定的答案！这相当于乘以置换矩阵 P 。例如，

$\begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} \cdot \begin{bmatrix} 10^{-20} & 1 \\ 1 & 1 \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ 10^{-20} & 1 \end{bmatrix}$

$PA = \hat{A}$

对 PA 应用高斯消元。

在每个步骤中，选择列 k 中的最大值，并将该行移动到行 k 。

作业

def swap(a,b):
    temp = np.copy(a)
    a[:] = b
    b[:] = temp

a=np.array([1,2,3])
b=np.array([3,2,1])
swap(a,b)
a,b

# 练习：重新编写上面的 LU 分解以使用交换主元

示例

A = np.array([[2,1,1,0],[4,3,3,1],[8,7,9,5],[6,7,9,8]]).astype(np.float)

L, U, P = LU_pivot(A)

可以比较下面 Trefethen，第 159 页的答案：

A

'''
array([[ 2.,  1.,  1.,  0.],
       [ 4.,  3.,  3.,  1.],
       [ 8.,  7.,  9.,  5.],
       [ 6.,  7.,  9.,  8.]])
'''

U

'''
array([[ 8.        ,  7.        ,  9.        ,  5.        ],
       [ 0.        ,  1.75      ,  2.25      ,  4.25      ],
       [ 0.        ,  0.        , -0.28571429,  0.57142857],
       [ 0.        ,  0.        ,  0.        , -2.        ]])
'''

P

'''
array([[ 0.,  0.,  1.,  0.],
       [ 0.,  0.,  0.,  1.],
       [ 1.,  0.,  0.,  0.],
       [ 0.,  1.,  0.,  0.]])
'''

部分交换主元可以置换行。这是一种普遍的做法，这通常是 LU 分解的意思。

完全交换主元可以置换行和列。完全交换主元非常耗时，很少在实践中使用。

示例

考虑方程组：

$\begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 1 \\ -1 & 1 & 0 & 0 & 0 & 1 \\ -1 & -1 & 1 & 0 & 0 & 1 \\ -1 & -1 & -1 & 1 & 0 & 1 \\ -1 & -1 & -1 & -1 & 1 & 1 \\ -1 & -1 & -1 & -1 & -1 & 1 \end{bmatrix} \mathbf{x} = \begin{bmatrix} 1 \\ 1 \\ 1 \\ 1 \\ 2 \\ 1 \end{bmatrix}$

def make_matrix(n):
    A = np.eye(n)
    for i in range(n):
        A[i,-1] = 1
        for j in range(i):
            A[i,j] = -1
    return A 

def make_vector(n):
    b = np.ones(n)
    b[-2] = 2
    return b

make_vector(7)

# array([ 1.,  1.,  1.,  1.,  1.,  2.,  1.])

练习

练习：让我们在 5×5 方程组上使用高斯消元法。

Scipy 也有这种功能。让我们看看最后 5 个方程的解，其中 n = 10,20,30,40,50,60 。

np.set_printoptions(precision=3, suppress=True)

?scipy.linalg.solve

for n, ls in zip(range(10, 70, 10), ['--', ':', '-', '-.', '--', ':']):
    soln = scipy.linalg.lu_solve(scipy.linalg.lu_factor(make_matrix(n)), make_vector(n))
    plt.plot(soln[-5:], ls)
    print(soln[-5:])

'''
[-0.062 -0.125 -0.25   0.5    1.002]
[-0.062 -0.125 -0.25   0.5    1.   ]
[-0.062 -0.125 -0.25   0.5    1.   ]
[-0.062 -0.125 -0.25   0.5    1.   ]
[-0.062 -0.125 -0.25   0.5    1.   ]
[ 0.  0.  0.  0.  1.]
'''