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

高斯消元

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

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

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

五、LU 分解 - 图1

答案：

五、LU 分解 - 图2

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

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

五、LU 分解 - 图3

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

  1. def LU(A):
  2.     U = np.copy(A)
  3.     m, n = A.shape
  4.     L = np.eye(n)
  5.     for k in range(n-1):
  6.         for j in range(k+1,n):
  7.             L[j,k] = U[j,k]/U[k,k]
  8.             U[j,k:n] -= L[j,k] * U[k,k:n]
  9.     return L, U
  11. A = np.array([[2,1,1,0],[4,3,3,1],[8,7,9,5],[6,7,9,8]]).astype(np.float)
  13. L, U = LU(A)
  15. np.allclose(A, L @ U)
  17. # True

LU分解很有用！

求解Ax = b变为LUx = b：

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

工作量

高斯消元的工作量：五、LU 分解 - 图4

内存

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

考虑矩阵：

五、LU 分解 - 图5

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

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

  1. # 练习：
  3. np.set_printoptions(suppress=True)
  5. # 练习：
  7. L2, U2 = LU(A)
  9. '''
  10. [[  1.00000000e-20   1.00000000e+00]
  11.  [  0.00000000e+00  -1.00000000e+20]]
  12. '''
  14. L2, U2
  16. '''
  17. (array([[  1.00000000e+00,   0.00000000e+00],
  18.         [  1.00000000e+20,   1.00000000e+00]]),
  19.  array([[  1.00000000e-20,   1.00000000e+00],
  20.         [  0.00000000e+00,  -1.00000000e+20]]))
  21. '''
  23. np.allclose(L1, L2)
  25. # True
  27. np.allclose(U1, U2)
  29. # True
  31. np.allclose(A, L2 @ U2)
  33. # False

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

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

稳定性

问题f的算法 五、LU 分解 - 图6 是稳定的，如果对于每个x：

五、LU 分解 - 图7

对于一些y：

五、LU 分解 - 图8

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

翻译：

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

向后稳定

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

问题f的算法 五、LU 分解 - 图9 是向后稳定的，如果对于每个x，

五、LU 分解 - 图10

对于一些y：

五、LU 分解 - 图11

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

翻译：

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

带有交换主元的 LU 分解

让我们看看矩阵：

五、LU 分解 - 图12

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

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

五、LU 分解 - 图13

五、LU 分解 - 图14

  1. L, U = LU(A)
  3. np.allclose(A, L @ U)
  5. # True

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

五、LU 分解 - 图15

五、LU 分解 - 图16

对PA应用高斯消元。

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

作业

  1. def swap(a,b):
  2.     temp = np.copy(a)
  3.     a[:] = b
  4.     b[:] = temp
  6. a=np.array([1,2,3])
  7. b=np.array([3,2,1])
  8. swap(a,b)
  9. a,b
  11. # 练习：重新编写上面的 LU 分解以使用交换主元

示例

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

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

  1. A
  3. '''
  4. array([[ 2.,  1.,  1.,  0.],
  5.        [ 4.,  3.,  3.,  1.],
  6.        [ 8.,  7.,  9.,  5.],
  7.        [ 6.,  7.,  9.,  8.]])
  8. '''
  10. U
  12. '''
  13. array([[ 8.        ,  7.        ,  9.        ,  5.        ],
  14.        [ 0.        ,  1.75      ,  2.25      ,  4.25      ],
  15.        [ 0.        ,  0.        , -0.28571429,  0.57142857],
  16.        [ 0.        ,  0.        ,  0.        , -2.        ]])
  17. '''
  19. P
  21. '''
  22. array([[ 0.,  0.,  1.,  0.],
  23.        [ 0.,  0.,  0.,  1.],
  24.        [ 1.,  0.,  0.,  0.],
  25.        [ 0.,  1.,  0.,  0.]])
  26. '''

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

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

示例

考虑方程组：

五、LU 分解 - 图17

  1. def make_matrix(n):
  2.     A = np.eye(n)
  3.     for i in range(n):
  4.         A[i,-1] = 1
  5.         for j in range(i):
  6.             A[i,j] = -1
  7.     return A 
  9. def make_vector(n):
  10.     b = np.ones(n)
  11.     b[-2] = 2
  12.     return b
  14. make_vector(7)
  16. # array([ 1.,  1.,  1.,  1.,  1.,  2.,  1.])

练习

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

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

  1. np.set_printoptions(precision=3, suppress=True)
  3. ?scipy.linalg.solve
  5. for n, ls in zip(range(10, 70, 10), ['--', ':', '-', '-.', '--', ':']):
  6.     soln = scipy.linalg.lu_solve(scipy.linalg.lu_factor(make_matrix(n)), make_vector(n))
  7.     plt.plot(soln[-5:], ls)
  8.     print(soln[-5:])
  10. '''
  11. [-0.062 -0.125 -0.25   0.5    1.002]
  12. [-0.062 -0.125 -0.25   0.5    1.   ]
  13. [-0.062 -0.125 -0.25   0.5    1.   ]
  14. [-0.062 -0.125 -0.25   0.5    1.   ]
  15. [-0.062 -0.125 -0.25   0.5    1.   ]
  16. [ 0.  0.  0.  0.  1.]
  17. '''

五、LU 分解 - 图18

当n = 60时会发生什么？

定理：让矩阵A的因式分解PA = LU通过高斯消元和部分交换主元来计算。 所得矩阵（由计算机使用浮点算术） 五、LU 分解 - 图19，五、LU 分解 - 图20 和 五、LU 分解 - 图21 满足：

五、LU 分解 - 图22

其中ρ是增长因子。

五、LU 分解 - 图23

对于我们上面的矩阵，五、LU 分解 - 图24。

理论上不稳定，实际上稳定

大多数算法（例如 QR）的稳定性很简单。 具有部分交换主元的高斯消元不是这种情况。 只有当L和/或U相对于A的大小较大时，才会出现高斯消元（有或没有交换主元）的不稳定性。

Trefethen：“尽管有（22.4）这样的例子，部分交换主元的高斯消元在实践中是完全稳定的……在计算的五十年中，在自然环境下不会出现产生爆炸性不稳定性的矩阵问题。”【虽然人为的例子很容易构造】

虽然有些矩阵会导致不稳定，但由于统计原因，占所有矩阵的比例非常小，因此“从不”出现。 “如果你随机挑选十亿个矩阵，你几乎肯定找不到高斯消元不稳定的矩阵。”

扩展阅读

  • 高斯消元/ LU 分解 - Trefethn 讲座 20
  • 交换主元 - Trefethn 讲座 21
  • 高斯消除的稳定性 - Trefethn 讲座 22

随机投影发生了什么？

我们在下面的矩阵中采用线性组合（带有随机权重）：

  1. plt.figure(figsize=(12, 12))
  2. plt.imshow(M, cmap='gray')
  4. # <matplotlib.image.AxesImage at 0x7f601f315fd0>

五、LU 分解 - 图25

这就像一个随机加权平均值。 如果你选取其中的几个，你最终会得到彼此不高度相关的列（大致正交）。

Johnson-Lindenstrauss 引理：（来自维基百科）高维空间中的一小组点可以嵌入到更低维度的空间中，使得点之间的距离几乎保持不变。

我们期望，能够以保留相关结构的方式，减少数据的维度。 Johnson-Lindenstrauss 引理是这种类型的经典结果。

高斯消元的历史

有趣的事实：高斯并没有发明高斯消元，但可能在 Cholesky 之前发现了 Cholesky 因子分解

— Rachel Thomas (@math_rachel) ) 2017 年 6 月 6 日

根据维基百科，Stigler 的 Eponymy 定律：“没有任何科学发现以它的原始发现者命名。例子包括哈勃定律，它是由 Georges Lemaître 在 Edwin Hubble 两年之前得到的，毕达哥拉斯定理在毕达哥拉斯之前为巴比伦数学家所知，哈雷彗星是自公元前 240 年以来天文学家观察到的彗星。Stigler 本人将社会学家 Robert K. Merton 命名为 Stigler 定律的发现者，表明它遵循自己的法令，尽管这一现象之前曾被其他人注意到。”

迷人的高斯消元的历史。一些亮点：

  • 公元前 20 0年左右，高斯消元的第一个书面记录在中文书籍“九章算术”中。
  • 古代中国人使用彩色竹棒放在“计数板”的列中。
  • 日本数学家 Seki Kowa（1643-1708）在 1683 年之前推进了中国的淘汰消元，并发明了行列式。大约在同一时间，莱布尼兹独立地发现了相似的发现，但是 Kowa 和莱布尼兹都没有因为他们的发现而受到赞扬。
  • 高斯称消元方法是“众所周知的”并且从未声称已经发明了它，尽管他可能已经发明了 Cholesky 分解。

这里有更多历史

加速高斯消元

并行 LU 分解：LU 分解可以完全并行化

随机化 LU 分解（2016 年文章）：随机 LU 完全为在标准 GPU 上运行而实现，无需任何 GPU-CPU 数据传输。

scipy.linalg vs lu_solve

  1. n = 60
  2. A = make_matrix(n)
  3. b = make_vector(n)

这个问题有很大的增长因子= 259。 我们使用scipy.linalg.lu_solve获得了错误的答案，但使用scipy.linalg.solve得到了正确的答案。什么是scipy.linalg.solve呢？

  1. print(scipy.linalg.lu_solve(scipy.linalg.lu_factor(A), b)[-5:])
  2. print(scipy.linalg.solve(A, b)[-5:])
  4. '''
  5. [ 0.  0.  0.  0.  1.]
  6. [-0.062 -0.125 -0.25   0.5    1.   ]
  7. '''
  9. %%timeit
  10. soln = scipy.linalg.lu_solve(scipy.linalg.lu_factor(A), b)
  11. soln[-5:]
  13. # 91.2 µs ± 192 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
  15. %%timeit
  16. soln = scipy.linalg.solve(A, b)
  17. soln[-5:]
  19. # 153 µs ± 5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

查看scipy的源代码，我们看到它正在调用LAPACK例程gesvx。 这是sgesvx的 Fortran 源代码（s指的是单个，也有用于浮点的dgesvx和复数的cgesvx）。 在注释中，我们看到它正在计算 reciprocal 主元增长因子，因此它考虑了这个增长因子，并做了一些比普通的部分主元 LU 分解更复杂的事情。

分块矩阵

经典的矩阵乘法

问题：计算两个n×n的矩阵A×B = C的矩阵乘法的计算复杂度（大O）是多少？

你可以在 Codecademy 学习（或复习）大O。

它的样子是：

  1. for i=1 to n
  2.     {read row i of A into fast memory}
  3.     for j=1 to n
  4.         {read col j of B into fast memory}
  5.         for k=1 to n
  6.             C[i,j] += A[i,k] x B[k,j]
  7.         {write C[i,j] back to slow memory}

问题：进行了多少次读写操作？

分块矩阵相乘

将A，B，C分成大小为N/n × N/n的N×N个块。

五、LU 分解 - 图26

来源

它的样子是：

  1. for i=1 to N
  2.     for j=1 to N
  3.         for k=1 to N
  4.             {read block (i,k) of A}
  5.             {read block (k,j) of B}
  6.             block (i,j) of C += block of A times block of B
  7.         {write block (i,j) of C back to slow memory}

问题 1：这个的大O是什么？

问题 2：进行了多少次读写操作？