scipy LU分解置换矩阵

Question

根据我的理解，LU分解是指一个矩阵A可以写成A = LU，其中L是一个下三角矩阵，U是一个上三角矩阵。

不过，在scipy中与LU分解相关的函数（lu、lu_factor、lu_solve）似乎涉及到一个第三个矩阵P，使得A = PLU，而P是一个置换矩阵（L和U还是之前的样子）。

那么这个置换矩阵有什么用呢？如果“真正的”LU分解总是可以实现，那为什么P会是除了单位矩阵以外的其他东西呢？

Answer 1

补充一下@DomJack的观点：改变排列顺序（也就是重新排序）也会影响L和U这两个部分中非零元素的数量。所以，重新排序可能会让计算变得更高效，尤其是在内存使用方面。

Answer 2

并不是所有的矩阵都有LU分解。不过，每个方阵至少可以通过调整行的顺序，找到一种LU分解的方式。

Answer 3

考虑一下高斯消元的过程。如果在主元位置上出现了零，你该怎么办呢？这时候你需要交换行，这就引入了一个P矩阵。

而且，如果主元的值非常小但不为零，这在浮点数计算中会导致数值不稳定。为了避免这种情况，基本的算法会在主元列中寻找绝对值最大的数，并将对应的行与主元行交换。

不过，这种交换可能会比较耗时，所以通常情况下，只有当绝对值最大的数比主元的绝对值大一个倍数，比如10，才会进行交换。这样可以减少交换的次数，但又能保留那些为了限制浮点数误差而必须进行的交换。

如果你想了解更多，可以搜索“带部分主元的LU分解”，会有很多不错的资源。

注意：因为P是一个置换矩阵，所以P的转置等于P的逆。也就是说，方程Ax = b和LUx = P^T b有相同的解（有些实现会返回你称之为P的矩阵，而其他的则返回你称之为P^T的矩阵并称之为P——一定要搞清楚是哪一个。这就是'PA = LU'和'A = PLU'之间的区别——在这两种情况下，P的含义是不同的）。