LU分解使用R

k4emjkb1  于 2023-09-27  发布在  其他
关注(0)|答案(1)|浏览(118)

我正在尝试使用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)
xzv2uavs

xzv2uavs1#

在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是单位的(不需要行交换)。

相关问题