scipy LU分解置换矩阵

Ale*_*lec 3 python numpy linear-algebra scipy

据我所知,LU分解,这意味着对于下三角矩阵L和上三角矩阵U,矩阵A可写为A = LU.

然而,在SciPy的关于LU因式分解的功能(lu,lu_factor,lu_solve)似乎涉及第三矩阵P,使得A = PLU,P是置换矩阵(和L,U如前是).

这个排列矩阵有什么意义?如果一个"真正的"LU分解始终是可能的,为什么P不是单位矩阵?

Dom*_*ack 8

考虑高斯消除过程.如果枢轴上有零,你会怎么做?您必须切换行,这会引入P矩阵.

此外,非常小的非零枢轴值导致浮点环境中的数值不稳定.基本算法通过在pivot列中搜索具有最大绝对值的条目并使用pivot行切换相应的行来避免这种情况.

这种开关可能很昂贵,因此通常最大的绝对值输入必须比枢轴的绝对值大一些因子,例如10,才能使开关发生.这减少了交换机的数量,但保留了限制浮点错误所必需的数量.

针对该问题的任意数量的良好资源搜索"LU分解与部分旋转".

注意:由于P是置换矩阵,P ^ T = P ^( - 1).因此,Ax = b与LUx = P ^ T b具有相同的解决方案(一些实现返回你所谓的P,而其他实现返回你所谓的P ^ T并称之为P - 确保你知道它是哪一个它是'PA = LU'和'A = PLU'之间的差异 - 在每种情况下P都不相同).