R和Python中的LU分解结果不一致 [英] Inconsistent results between LU decomposition in R and Python
问题描述
我在R中具有以下矩阵A:
I have the following matrix A in 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分解如下:
I compute its LU-decomposition as follows:
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代码:
On the other hand, this is the Python code for the same job:
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).你知道我在做什么错吗?
I'm wondering why the two C and lu matrices in R and Python (respectively) are not the same. The point is that I have to get the same result as Python version (i.e. matrix lu). Do you have any idea what I am doing wrong?
推荐答案
在1.5年后意识到我的原始答案并不完全正确是令人尴尬的.尽管正确指出问题中A的秩不足是原因，但将其归因于根本原因是不正确的. 真正的问题是枢轴的非唯一选择，即使A排名较高，也有可能发生(尽管可能性较小).鉴于该帖子已被查看700多次，并且获得6分，我可能误导了许多读者.抱歉！
It is embarrassing to realize 1.5 years later that my original answer is not fully correct. While it correctly pointed out that the rank-deficiency of A in the question is the cause, it is incorrect to attribute this as the root cause. The real issue is the non-unique choice of pivot, which could happen (although less likely) even if A is full-rank. Given that this post has been viewed 700+ times and has a score of 6, I might have misled many readers. SORRY!
我发布了写了一个可跟踪的R函数，该函数模仿LAPACK的dgetrf进行LU分解，并回答了它.这个问题包含一个LU函数，用于不进行透视分解的LU分解，答案包含两个函数LUP和LUP2，用于具有行透视的版本，与LAPACK的 dgetrf 一致， Matrix::lu和R基本函数solve.的密集方法特别是LUP2函数允许人们逐步跟踪因式分解.我将在这里使用此功能进行调查.
I posted Write a trackable R function that mimics LAPACK's dgetrf for LU factorization and answered it just now. The question contains a LU function for LU factorization without pivoting, and the answer contains two functions LUP and LUP2 for a version with row pivoting that are consistent with LAPACK's dgetrf, which underlies the dense method of Matrix::lu and R base function solve. In particular, the LUP2 function allows one to track the factorization step by step. I will use this function here for investigation.
所以您要分解A的转置.
从R和Python的输出中，我们看到它们产生相同的第一个枢轴-1.1527778和第二个枢轴-1.2746988，而第三个枢轴却不同.因此，当分解完成前两列/行并进行到第三列/行时，肯定会发生一些有趣的事情.现在让我们暂停R的LU分解:
From the output of R and Python we see that they yield identical 1st pivot -1.1527778 and 2nd pivot -1.2746988, while the 3rd pivot differs. So there is definitely something interesting happening when the factorization has done the first two columns / rows and proceeds to the the 3rd column / row. Let's pause R's LU factorization at this point:
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)减小为
At this point, the Gaussian elimination has reduced t(A) to
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可以是枢轴，因为它们具有相同的绝对值.
Wow, we see something really interesting now: the 3rd and 4th rows only differ by a sign change. Then the Gaussian elimination is not unique, because either -1.2291115 or 1.2291115 can be a pivot as they have the same absolute value.
很明显，R选择了-1.2291115作为枢轴，但是Python选择了1.2291115作为枢轴. Python中将在第3行和第4行之间进行行交换.在您的问题中，您没有报告Python给您的置换索引，但是它应该1, 2, 4, 3，而不是R中的1, 2, 3, 4.
Clearly, R has chosen -1.2291115 as the pivot, but Python has chosen 1.2291115 as the pivot. A row exchange between 3rd and 4th rows will happen in Python. In your question, you didn't report what permutation index Python gives you, but it should 1, 2, 4, 3, instead of the 1, 2, 3, 4 in R.
无论哪种方式，U因数最终都在底部以一排零结尾，因此t(A)或A并非完全排名.如果要在两个软件之间看到更一致的行为，则最好尝试使用完整等级矩阵.在这种情况下，在LU分解期间，您不太可能有一个以上的枢轴选择.您可以在R中生成一个随机的全秩矩阵，
Either way, the U factor ends up with a row of zeros in the bottom, so that t(A) or A is not full-rank. If you want to see a more consistent behavior between two software, you'd better try a full-rank matrix. In that case, it is less likely that you will have more than one choices for a pivot during LU factorization. You can generate a random full-rank matrix in R by
A <- matrix(runif(16), 4, 4)
这篇关于R和Python中的LU分解结果不一致的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
