R和Python中LU分解之间的结果不一致

我A在R中有以下矩阵:

           # [,1]       [,2]   [,3]       [,4]
# [1,] -1.1527778  0.4444444  0.375  0.3333333
# [2,]  0.5555556 -1.4888889  0.600  0.3333333
# [3,]  0.6250000  0.4000000 -1.825  0.8000000
# [4,]  0.6666667  0.6666667  0.200 -1.5333333

A <- structure(c(-1.15277777777778, 0.555555555555556, 0.625, 0.666666666666667, 
0.444444444444444, -1.48888888888889, 0.4, 0.666666666666667, 
0.375, 0.6, -1.825, 0.2, 0.333333333333333, 0.333333333333333, 
0.8, -1.53333333333333), .Dim = c(4L, 4L), .Dimnames = list(NULL, 
    NULL))

我按如下方式计算其LU分解:

library(Matrix)
ex <- expand(lu(t(A)))
L <- ex$L
P <- ex$P
C <- U <- ex$U
C[lower.tri(U)] <- L[lower.tri(L)]

print(C)

# 4 x 4 Matrix of class "dgeMatrix"
           # [,1]       [,2]       [,3]          [,4]
# [1,] -1.1527778  0.5555556  0.6250000  6.666667e-01
# [2,] -0.3855422 -1.2746988  0.6409639  9.236948e-01
# [3,] -0.3253012 -0.6124764 -1.2291115  9.826087e-01
# [4,] -0.2891566 -0.3875236 -1.0000000 -2.220446e-16

另一方面,这是同一作业的Python代码:

lu, piv = scipy.linalg.lu_factor(A.T, check_finite=False)

print(lu)

# [[ -1.15277778e+00   5.55555556e-01   6.25000000e-01   6.66666667e-01]
 # [ -3.85542169e-01  -1.27469880e+00   6.40963855e-01   9.23694779e-01]
 # [ -2.89156627e-01  -3.87523629e-01   1.22911153e+00  -9.82608696e-01]
 # [ -3.25301205e-01  -6.12476371e-01  -1.00000000e+00   7.69432827e-16]]

我想知道为什么R和Python(分别)中的两个C和lu矩阵不一样.关键是我必须得到与Python版本(即矩阵lu)相同的结果.你知道我做错了什么吗？

1.5年后意识到我的原始答案并不完全正确,这令人尴尬.虽然它正确地指出问题中的等级缺陷A是原因,但将其归因于根本原因是不正确的.真正的问题是枢轴的非唯一选择,即使A是满级,也可能发生(虽然不太可能).鉴于这篇文章被观看了700多次,得分为6,我可能误导了许多读者.抱歉!

我发布了一个可跟踪的R函数,它模仿LAPACK的dgetrf进行LU分解,并且刚刚回答了它.该问题包含一个LU没有旋转的LU分解函数,答案包含两个函数LUP,LUP2对于一个与LAPACK的dgetrf一致的行旋转版本,它是密集方法Matrix::lu和R基函数的基础solve.特别是,该LUP2函数允许跟踪分解是一步一步的.我将在此处使用此功能进行调查.

所以你是因素分解的转置A.

从R和Python的输出中我们看到它们产生相同的第一枢轴-1.1527778和第二枢轴-1.2746988,而第三个枢轴不同.因此,当因子分解完成前两列/行并进入第三列/行时,肯定会发生一些有趣的事情.让我们暂停R的LU分解:

oo <- LUP2(t(A), to = 2)
#           [,1]       [,2]       [,3]       [,4]
#[1,] -1.1527778  0.5555556  0.6250000  0.6666667
#[2,] -0.3855422 -1.2746988  0.6409639  0.9236948
#[3,] -0.3253012 -0.6124764 -1.2291115  0.9826087
#[4,] -0.2891566 -0.3875236  1.2291115 -0.9826087
#attr(,"to")
#[1] 2
#attr(,"pivot")
#[1] 1 2 3 4

此时,高斯消元已减少t(A)到

getU(oo)
#           [,1]       [,2]       [,3]       [,4]
#[1,] -1.1527778  0.5555556  0.6250000  0.6666667
#[2,]  0.0000000 -1.2746988  0.6409639  0.9236948
#[3,]  0.0000000  0.0000000 -1.2291115  0.9826087
#[4,]  0.0000000  0.0000000  1.2291115 -0.9826087
#attr(,"to")
#[1] 2

哇,我们现在看到一些非常有趣的东西:第3行和第4行仅因符号变化而不同.然后高斯消除不是唯一的,因为它们具有相同的绝对值,-1.2291115或者1.2291115可以是枢轴.

显然,R已选择-1.2291115作为支点,但Python已选择1.2291115作为支点.Python中将发生第3行和第4行之间的行交换.在你的问题中,你没有报告Python给你的排列索引,但它应该1, 2, 4, 3,而不是1, 2, 3, 4在R.

无论哪种方式,U因子最终会在底部有一行零,因此t(A)或者A不是全等级.如果你想在两个软件之间看到更一致的行为,你最好尝试一个全秩矩阵.在这种情况下,在LU分解期间,您不太可能有多个枢轴选择.您可以在R中生成随机满秩矩阵

A <- matrix(runif(16), 4, 4)