函数 scipy.linalg.lu 中的上三角矩阵是否始终为行梯形？答案

【问题标题】：Is the upper triangular matrix in function scipy.linalg.lu always in row echelon form?函数 scipy.linalg.lu 中的上三角矩阵是否始终为行梯形？
【发布时间】：2014-09-30 22:24:24
【问题描述】：

我有一个 m x n 矩阵 A，n > m，我试图通过它的行梯形来识别独立的行。函数 scipy.linalg.lu 返回我的矩阵的 PLU 分解，但 U 因子似乎不是梯形形式，即枢轴不是阶梯模式。据我所知，U 因子应该总是呈阶梯状。

考虑以下示例：

from numpy import array
from scipy.linalg import lu

A = array([[1, 1, 1, 1, 0, 1, 1, 1, 1, 0],
           [1, 1, 0, 0, 1, 0, 1, 0, 1, 1],
           [1, 1, 0, 0, 0, 1, 1, 0, 0, 0],
           [0, 1, 0, 1, 1, 0, 1, 0, 0, 1],
           [1, 1, 0, 0, 1, 1, 1, 1, 1, 1]])

P, L, U = lu(A)

U 因子不是行梯形。对于每一行 k，枢轴应始终位于第 k-1 行中枢轴的右侧。看到第五行的枢轴不在第四行的枢轴右侧：

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

【问题讨论】：

请提供示例代码。
顺便说一句，我已经解决了通过 QR 分解来识别独立行的问题，虽然我的第一次尝试是按照线性代数教科书通过高斯消去法来识别它们。
scipy.linalg.lu (docs.scipy.org/doc/scipy/reference/generated/…) 并没有说 U 将采用行梯形 (en.wikipedia.org/wiki/Row_echelon_form)。它说它将是上三角形（mathworld.wolfram.com/UpperTriangularMatrix.html）。

标签： python numpy matrix scipy linear-algebra

【解决方案1】：

我以前从未见过 LUP 分解，但我怀疑它并没有完全按照您的想法进行。我在 Matlab 中进行了与您在 Python 中所做的相同的分解，并且得到了完全相同的结果。所以我不认为 Python LU 函数有问题。

更新：我也在 R（下面的代码）中实现了这一点，并且 U 矩阵再次给出了与 Matlab 和 Python 相同的结果。与其他人不同，它给出以下警告：

Warning message:
In .local(x, ...) :
  Exact singularity detected during LU decomposition: U[i,i]=0, i=4.

Matlab：

代码：

A = ([[1, 1, 1, 1, 0, 1, 1, 1, 1, 0],
           [1, 1, 0, 0, 1, 0, 1, 0, 1, 1],
           [1, 1, 0, 0, 0, 1, 1, 0, 0, 0],
           [0, 1, 0, 1, 1, 0, 1, 0, 0, 1],
           [1, 1, 0, 0, 1, 1, 1, 1, 1, 1]]);
[L,U,P] = lu(A);
U

Matlab 结果：
你=

     1     1     1     1     0     1     1     1     1     0
     0     1     0     1     1     0     1     0     0     1
     0     0    -1    -1     0     0     0    -1    -1     0
     0     0     0     0     1    -1     0     0     1     1
     0     0     0     0     1     0     0     1     1     1

Python 结果（来自您的代码）：

U = 

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

R：

代码：

library(Matrix)
A <- Matrix( c( 1, 1, 1, 1, 0, 1, 1, 1, 1, 0 , 1, 1, 0, 0, 1, 0, 1, 0, 1, 1 , 1, 1, 0, 0, 0, 1, 1, 0, 0, 0 , 0, 1, 0, 1, 1, 0, 1, 0, 0, 1 , 1, 1, 0, 0, 1, 1, 1, 1, 1, 1 ),nrow=5,ncol=10,byrow=TRUE )
expand(mylu <-lu(A))

结果：

class "dtrMatrix" (unitriangular)
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    .    .    .    .
[2,]    0    1    .    .    .
[3,]    1    0    1    .    .
[4,]    1    0    1    1    .
[5,]    1    0    1    0    1

$U
5 x 10 Matrix of class "dgeMatrix"
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    1    1    1    1    0    1    1    1    1     0
[2,]    0    1    0    1    1    0    1    0    0     1
[3,]    0    0   -1   -1    0    0    0   -1   -1     0
[4,]    0    0    0    0    1   -1    0    0    1     1
[5,]    0    0    0    0    1    0    0    1    1     1

$P
5 x 5 sparse Matrix of class "pMatrix"

[1,] | . . . .
[2,] . . . | .
[3,] . . | . .
[4,] . | . . .
[5,] . . . . |

【讨论】：

这很奇怪。根据 U 因子 (STRANG, 1988, p.72) 的以下属性，第五列中枢轴下方的元素应为零，而不是 1：“(i) 非零行首先出现 - 否则将有行交换 - 枢轴是这些行中的第一个非零条目。（ii）每个枢轴下方是一列零，通过消除获得。（iii）每个枢轴位于上一行中枢轴的右侧。这会产生楼梯图案。” STRANG, G. 线性代数及其应用。第三版，汤姆森，1988 年。
事实上，我应用 LU 分解的矩阵比上面示例中给出的矩阵要大得多，而阶梯模式中的 U 因子远小于上面的。我目前的假设是 Python 和 Matlab 中的 LU 函数计算的 U 因子中的行没有正确排序。
我刚刚也在 R 中进行了分解（请参阅更新的答案），它也给出了相同的结果。
R 给出的奇点警告只对方阵有意义。矩形矩阵在消除后很可能在对角线位置有零。与方阵不同，对角线中的零并不意味着矩阵秩不足或没有逆矩阵。实际上，示例中的矩阵 A 是满秩的。
有趣的是，如果我取消项 U[5,5]，通过减去第 4 行和第 5 行，我得到了预期的行梯形，但是乘以因子 PLU 不会恢复 A。这似乎违反了 Strang (1998, p.72) 中的 (iii)。我认为这更像是一个线性代数问题，而不是计算问题。我将在数学交流中发布这个问题。

【解决方案2】：

不能保证 LU 分解会产生行梯形的 U 矩阵。可能在这种情况下失败的原因是因为矩阵是奇异的；我强烈怀疑这就是原因，因为奇异矩阵的一个特征是没有完整的枢轴集。

请参阅https://en.wikipedia.org/wiki/LU_decomposition 上的通用矩阵部分以及其中提供的参考以了解详细信息。

【讨论】：