使用 R 进行 LU 分解

Question

我正在尝试使用 R 运行 LU 分解。这是可重现的代码。我不明白为什么我的排列矩阵与解决方案不同。L 和 U 矩阵是正确的。但对于置换矩阵，第 1 行和第 2 行以及第 3 行和第 4 行互换。因此，我没有得到线性方程组的正确解。将不胜感激您的帮助。

install.packages("Matrix")
library(Matrix)
(A <- matrix(c(4, 3, -2, 5, 2, -4, 6, 1, -1, 2, -5, 6, 3, 5, -2, -3), nrow = 4))
(B <- matrix(c(16.9, -14, 25, 9.4), nrow = 4))

luA <- lu(A)
elu <- expand(luA)
(L <- elu$L)
(U <- elu$U)
(P <- elu$P)

(Y <- solve(L) %*% P %*% B)
(X <- solve(U) %*% Y)

Answer 1

通过 R 实现，我们似乎有A = PLU（而不是PA = LU）。因此，进行以下工作：

all.equal(Matrix(A), with(elu, P %*% L %*% U))
# TRUE

(Y <- solve(L, solve(P) %*% B)) # solve LY = inv(P).B instead of LY = PB
(X <- solve(U, Y))
 
X
#4 x 1 Matrix of class "dgeMatrix"
#     [,1]
#[1,]  4.5
#[2,]  1.6
#[3,] -3.8
#[4,] -2.7 

all.equal(X, Matrix(solve(A, B)))
# TRUE

为了进行双重检查，对于给定的和，我们可以使用高斯消去法计算LU分解：AB

# Gaussian elimination steps
E21 <- matrix(c(1,0,0,0,-3/4,1,0,0,0,0,1,0,0,0,0,1), nrow=4, byrow=T)
E31 <- matrix(c(1,0,0,0,0,1,0,0,1/2,0,1,0,0,0,0,1), nrow=4, byrow=T)
E41 <- matrix(c(1,0,0,0,0,1,0,0,0,0,1,0,-5/4,0,0,1), nrow=4, byrow=T)
E32 <- matrix(c(1,0,0,0,0,1,0,0,0,14/11,1,0,0,0,0,1), nrow=4, byrow=T)
E42 <- matrix(c(1,0,0,0,0,1,0,0,0,0,1,0,0,-3/11,0,1), nrow=4, byrow=T)
E43 <- matrix(c(1,0,0,0,0,1,0,0,0,0,1,0,0,0,13/4,1), nrow=4, byrow=T)

U <- E43 %*% E42 %*% E32 %*% E41 %*% E31 %*% E21 %*% A
L <- solve(E43 %*% E42 %*% E32 %*% E41 %*% E31 %*% E21)

all.equal(L %*% U, A)
# [1] TRUE

all.equal(solve(U, solve(L, B)), solve(A, B))
# [1] TRUE

这里的置换矩阵 P 是恒等的（不需要行交换）。