R和Python中LU分解的不一致结果
我在R中有以下矩阵A：
# [,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中的两个C矩阵和lu矩阵(分别)不是相同的。关键是，我必须获得与Python相同的结果(即矩阵lu)。你知道我做错了什么吗？
回答 1
Stack Overflow用户
发布于 2017-01-10 14:27:09
一年半后才意识到我原来的答案并不完全正确，实在令人尴尬。虽然它正确地指出了A在问题中的等级缺陷是原因，但把它作为根本原因是不正确的。真正的问题是支点的非唯一选择，即使A是全等级的，也可能发生(尽管可能性较小)。鉴于这篇文章已经被“700+时报”看过，得分为6分，我可能误导了许多读者。抱歉的!
我刚刚发了Write a trackable R function that mimics LAPACK's dgetrf for LU factorization并回复了它。该问题包含一个用于不旋转的LU分解的LU函数，而对于一个与LAPACK的dgetrf一致的版本，答案包含了LUP和LUP2两个函数，后者是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哇，我们现在看到了一些非常有趣的东西:第三行和第四排只有一个符号变化才能区别。然后高斯消除并不是唯一的，因为无论是-1.2291115还是1.2291115都可以作为支点，因为它们具有相同的绝对值。
显然，R选择了-1.2291115作为枢轴，但是Python选择了1.2291115作为枢轴。Python中将发生第3行和第4行之间的行交换。在您的问题中，您没有报告Python提供的置换索引，但是它应该是1, 2, 4, 3，而不是R中的1, 2, 3, 4。
无论哪种方式，U因子的底部都有一行零，因此t(A)或A并不是满的。如果您想在两个软件之间看到更一致的行为，最好尝试一个完整的矩阵。在这种情况下，在LU分解过程中，不太可能有一个以上的支点选择。你可以在R中生成一个随机满秩矩阵
A <- matrix(runif(16), 4, 4)https://stackoverflow.com/questions/41570787
复制相似问题
腾讯云开发者
