函数scipy.linalg.lu中的上三角矩阵总是行阶梯形吗？

Question

我有一个 m x n 的矩阵 A，其中 n 大于 m，我想通过这个矩阵的行阶梯形态来找出独立的行。函数 scipy.linalg.lu 返回了我的矩阵的 PLU 分解，但 U 这个部分似乎并不是行阶梯形态，也就是说，主元（pivot）并没有呈现出阶梯状的排列。根据我的了解，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.]])

Answer 1

LU分解并不能保证得到的U矩阵是行阶梯形的。可能出现这种情况的原因是矩阵是奇异的；我非常怀疑这就是原因，因为奇异矩阵的一个特征就是缺少完整的主元。

有关一般矩阵的详细信息，请查看https://en.wikipedia.org/wiki/LU_decomposition中的相关部分和提供的参考资料。

Answer 2

我以前没见过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结果：
U =

     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,] . . . . |