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

Question

我在 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 和 Python（分别）中的两个 C 和 lu 矩阵不相同。关键是我必须得到与 Python 版本相同的结果（即矩阵 lu）。你知道我做错了什么吗？

Answer 1

1.5 年后才发现我原来的答案并不完全正确，这让我很尴尬。虽然它正确地指出问题中 A 的排名不足是原因，但将其归因于根本原因是不正确的。 真正的问题是主元 的非唯一选择，即使 A 是满秩也可能发生（尽管不太可能）。鉴于此 post 已被浏览 700+ 次且评分为 6，我可能误导了很多读者。抱歉！

我post编辑了，刚才回答了。该问题包含一个 LU 函数用于 LU 分解而不旋转，答案包含两个函数 LUP 和 LUP2 用于与 LAPACK 的 dgetrf 一致的行旋转版本，它是 Matrix::lu 和 R 基函数 solve. 的密集方法的基础，特别是 LUP2 函数允许人们逐步跟踪因式分解。我将在这里使用此功能进行调查。

所以你正在分解 A 的转置。

从 R 和 Python 的输出我们看到它们产生相同的第一个枢轴 -1.1527778 和第二个枢轴 -1.2746988，而第三个枢轴不同。因此，当分解完成前两列/行并进行到第 3 列/行时，肯定会发生一些有趣的事情。让我们在此时暂停 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都可以作为一个pivot，因为它们的绝对值是一样的。

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

中的 1, 2, 3, 4

无论如何，U 因子在底部以一行零结尾，因此 t(A) 或 A 不是满秩。如果您想在两个软件之间看到更一致的行为，您最好尝试使用满秩矩阵。在这种情况下，您不太可能在 LU 分解期间有多个主元选择。您可以通过

在 R 中生成随机满秩矩阵

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

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

Inconsistent results between LU decomposition in R and Python

python

r

matrix

matrix-decomposition